%% Exercises from Linear algebra, signal processing, and wavelets. A unified approach.\\ Matlab version
% 
% Author: \O yvind Ryan
% 
% Date: Mar 7, 2017
% 
% % Externaldocuments: applinalg
% % Mapping from exercise labels to numbers: label2numbers = {'example:extractcolours': '9.2', 'sec:rgb-grey': '9.3', 'example:secondorder_derivative': '9.11', 'ex:contrastimages3': '9.8', 'ex:chesspatternmoluecules': '9.12', 'ex:transform2jpeg': '9.28', 'example:edge_detection': '9.10', 'ex:genimgs': '9.14', 'smoothingex': '9.9', 'example:negative': '9.4', 'example:contrast': '9.5', 'ex:impl_twodim_DFT_DCT': '9.27', 'ex:contrastex': '9.7'}
% 
% 
% Table of contents:
% 
%  Sound and Fourier series 
%  Digital sound and Discrete Fourier analysis 
%  Operations on digital sound: digital filters 
%  Symmetric filters and the DCT 
%  Motivation for wavelets and some simple examples 
%  The filter representation of wavelets 
%  Constructing interesting wavelets 
%  The polyphase representation and wavelets 
%  Digital images 
%      Example 9.2: Extracting the different colors 
%      Example 9.3: Converting from color to grey-level 
%      Example 9.4: Computing the negative image 
%      Example 9.5: Increasing the contrast 
%      Exercise 9.7: Adjust contrast in images 
%      Exercise 9.8: Adjust contrast with another function 
%      Example 9.9: Smoothing an image 
%      Example 9.10: Edge detection 
%        Partial derivative in $x$-direction 
%        Partial derivative in $y$-direction 
%        The gradient 
%      Example 9.11: Second-order derivatives 
%      Example 9.12: Chess pattern image 
%      Exercise 9.14: Generate images 
%      Exercise 9.27: Implement two-dimensional FFT and DCT 
%      Exercise 9.28: Zeroing out DCT coefficients 
%  Using tensor products to apply wavelets to images 
%  Appendix: Basic Linear Algebra 
%  Appendix: Signal processing and linear algebra: a translation guide 
% 
% 
% 
%% Sound and Fourier series
%% Digital sound and Discrete Fourier analysis
%% Operations on digital sound: digital filters
%% Symmetric filters and the DCT
%% Motivation for wavelets and some simple examples
%% The filter representation of wavelets
%% Constructing interesting wavelets
%% The polyphase representation and wavelets
%% Digital images
% 
% 
% % --- begin exercise ---
% 
%% Example 9.2: Extracting the different colors
% 
% % keywords = ipynbimageschap; mimageschap
% 
% If we have a color image 
% $$P=(r_{i,j},g_{i,j},b_{i,j})_{i,j=1}^{m,n},$$ 
% it is often useful to manipulate the three color components separately as the three images
% $$P_r=(r_{i,j})_{i,j=1}^{m,n},\quad P_r=(g_{i,j})_{i,j=1}^{m,n},\quad P_r=(b_{i,j})_{i,j=1}^{m,n}.$$
% As an example, let us first see how we can produce three separate images, showing the R,G, and B color components, respectively. 
% Let us take the image |lena.png| used in Figure 9.1: 
% 
img = double(imread('lena.png'));
% The returned object has three dimensions. The first
% two dimensions represent the spatial directions (the row-index and column-index). The third dimension represents the color component.
% One can therefore view images representing the different color components with the help of the following code: 
% 
X1 = zeros(size(img));
X1(:,:,1) = img(:,:,1);
X2 = zeros(size(img));
X2(:,:,2) = img(:,:,2);
X3=zeros(size(img));
X3(:,:,3) = img(:,:,3);
% 
imshow(uint8(X1))
imshow(uint8(X2))
imshow(uint8(X3))
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 9.3: Converting from color to grey-level
% 
% % keywords = ipynbimageschap; mimageschap
% 
% If we have a color image we can convert it to a grey-level image. This means that at each point in the image we have to replace the three color values 
% $(r,g,b)$ by a single value $p$ that will represent the grey level. If we want the grey-level image to be a reasonable representation of the color image, 
% the value $p$ should somehow reflect the intensity of the image at the point. There are several ways to do this.
% 
% It is not unreasonable to use the largest of the three color components as a measure of the intensity, i.e, to set $p=\max(r,g,b)$. 
% An alternative is to use the sum of the three values as a measure of the total intensity at the point. This corresponds to setting $p=r+g+b$. 
% Here we have to be a bit careful with a subtle point. We have required each of the $r$, $g$ and $b$ values to lie in the range $[0,1]$, 
% but their sum may of course become as large as 3. We also require our grey-level values to lie in the range $[0,1]$ so after having computed all the sums we must 
% normalise as explained above.
% A third possibility is to think of the intensity of $(r,g,b)$ as the length of the color vector, in analogy with points in space, and set $p=\sqrt{r^2+g^2+b^2}$. 
% Again, we may end up with values in the range $[0,\sqrt{3}]$ so we have to normalise like we did in the second case. 
% 
% Let us sum this up as follows: A color image $P=(r_{i,j},g_{i,j},b_{i,j})_{i,j=1}^{m,n}$ can be converted to a grey level image $Q=(q_{i,j})_{i,j=1}^{m,n}$ by one of the following three operations:
% 
% * Set $q_{i,j}=\max(r_{i,j},g_{i,j},b_{i,j})$ for all $i$ and $j$.
% 
% * Compute $\hat q_{i,j}=r_{i,j}+g_{i,j}+b_{i,j}$ for all $i$ and $j$.
% 
% * Transform all the values to the interval $[0,1]$ by setting $$q_{i,j}=\frac{\hat q_{i,j}}{\max_{k,l} \hat q_{k,l}}.$$
% 
% * Compute $\hat q_{i,j}=\sqrt{r_{i,j}^2+g_{i,j}^2+b_{i,j}^2}$ for all $i$ and $j$.
% 
% * Transform all the values to the interval $[0,1]$ by setting $$q_{i,j}=\frac{\hat q_{i,j}}{\max_{k,l} \hat q_{k,l}}.$$
% 
% In practice one of the last two methods are preferred, perhaps with a preference for the last method, but the actual choice depends on the application.
% These can be implemented as follows. 
% 
X1 = max(img, [], 3);
X2 = img(:, :, 1) + img(:, :, 2) + img(:, :, 3);
X2 = 255*mapto01(X2);
X3 = sqrt(img(:,:,1).^2 + img(:,:,2).^2 + img(:,:,3).^2);
X3 = 255*mapto01(X3);
% 
imshow(uint8(X1))
imshow(uint8(X2))
imshow(uint8(X3))
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 9.4: Computing the negative image
% 
% % keywords = ipynbimageschap; mimageschap
% 
% In film-based photography a negative image was obtained when the film was developed, and then a positive image was created from the negative. 
% We can easily simulate this and compute a negative digital image.
% 
% Suppose we have a grey-level image $P=(p_{i,j})_{i,j=1}^{m,n}$ with intensity values in the interval $[0,1]$. Here intensity value 0 corresponds to black and 1 
% corresponds to white. To obtain the negative image we just have to replace an intensity $p$ by its 'mirror value' $1-p$.
% This is also easily translated to code as above. 
% 
X1 = 255 - X1;
X2 = 255 - X2;
X3 = 255 - X3;
imshow(uint8(X1))
imshow(uint8(X2))
imshow(uint8(X3))
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 9.5: Increasing the contrast
% 
% % keywords = ipynbimageschap; mimageschap
% 
% A common problem with images is that the contrast often is not good enough. 
% This typically means that a large proportion of the grey values are concentrated in a rather small subinterval of $[0,1]$. 
% The obvious solution to this problem is to somehow spread out the values. 
% This can be accomplished by applying a monotone function $f$ which maps $[0,1]$ onto $[0,1]$. 
% If we choose $f$ so that its derivative is large in the area where many intensity values are concentrated, we obtain the desired effect. 
% We will consider two such families of functions:
% 
% $$f_n(x)        &= \frac{\arctan(n(x-1/2))}{2\arctan(n/2)}+ \frac{1}{2}  \\  g_\epsilon(x) &= \frac{\ln(x+\epsilon)-\ln \epsilon}{\ln(1+\epsilon)-\ln \epsilon}.$$
% The first type of functions have quite large derivatives near $x=0.5$ and will therefore increase the contrast in images with a concentration of intensities with value around 0.5. 
% The second type of functions have a large derivative near $x=0$ and will therefore increase the contrast in images with a large proportion of small intensity values, i.e., very dark images. 
% The following code plots the functions $f_4$, $f_{10}$, and $f_{100}$. 
% 
x=linspace(0, 1, 100);
n=4;
f4 = atan(n*(x-1/2))/(2*atan(n/2)) + 1/2;
n=10;
f10=atan(n*(x-1/2))/(2*atan(n/2)) + 1/2;
n=100;
f100=atan(n*(x-1/2))/(2*atan(n/2)) + 1/2;
figure()
plot(x,f4,'r',x,f10,'g',x,f100,'b');
axis([0 1 0 1]);
legend('n=4','n=10','n=100');
% The following code plots the functions $g_{0.1}$, $g_{0.01}$, and $g_{0.001}$. 
% 
eps=0.1;
g01=(log(x+eps)-log(eps))./(log(1+eps)-log(eps));
eps=0.01;
g001=(log(x+eps)-log(eps))./(log(1+eps)-log(eps));
eps=0.001;
g0001=(log(x+eps)-log(eps))./(log(1+eps)-log(eps));
figure();
plot(x,g01,'r',x,g001,'g',x,g0001,'b');
axis([0 1 0 1]);
legend(['$\epsilon$=0.1','$\epsilon$=0.01','$\epsilon$=0.001'])
% The following code applies $f_{10}$ to the image in the right part of Figure 9.6. 
% 
X1 = contrastadjust0(img, 10);
imshow(uint8(X1))
% The following code applies $g_{0.01}$ to the image in the right part of Figure 9.6. 
% 
X2 = contrastadjust(img, 0.01);
imshow(uint8(X2))
% 
% Since the image was quite well balanced, $f_{10}$ made the dark areas too dark and the bright areas too bright.
% $g_{0.01}$ on the other hand has made the image as a whole too bright.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Exercise 9.7: Adjust contrast in images
% 
% % keywords = ipynbimageschap; mimageschap
% 
% 
% *a)* Write a function |contrastadjust0| which instead uses the function from Equation (9.1) 
% to increase the contrast. $n$ should be a parameter to the function.
% 
% *b)* Generate the left and right images in Figure 9.9 
% on your own by writing code which uses the two functions |contrastadjust0| and |contrastadjust|.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The following code can be used.
% 
X1 = contrastadjust0(img, 10);
imshow(uint8(X1))
% 
X2 = contrastadjust(img, 0.01);
imshow(uint8(X2))
% 
% % --- end solution of exercise ---
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Exercise 9.8: Adjust contrast with another function
% 
% % keywords = ipynbimageschap; mimageschap
% 
% In this exercise we will look at another function for increasing the contrast of a picture.
% 
% 
% *a)* Show that the function $f: \R \rightarrow \R$ given by
% $$f_n(x) = x^n,$$
% for all $n$ maps the interval $[0,1] \rightarrow [0,1]$, and that $f^\prime (1) \rightarrow \infty$ as $n \rightarrow \infty$.
% 
% *b)* The color image \verb secret.jpg, shown in Figure 9.10, 
% contains some information that is nearly invisible to the naked eye on most computer monitors. Use
% the function $f(x)$, to reveal the secret message.
% 
% 
% <<images/secret.jpg>>
% 
% % --- begin hint in exercise ---
% 
% *Hint.* You will first need to convert the image to a greyscale image. You can then use the function |contrastadjust| as a starting point for your own program.
% 
% % --- end hint in exercise ---
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The following code can be used
% 
n = 100;
X = double(imread('images/secret.jpg', 'jpg'));
X = (X(:,:,1) + X(:,:,2) + X(:,:,3))/3;
X = mapto01(X);
X = X.^n;
X = 255*X;
imshow(uint8(X))
% 
% % --- end solution of exercise ---
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 9.9: Smoothing an image
% 
% % keywords = ipynbimageschap; mimageschap
% 
% When we considered filtering of digital sound, low-pass filters dampened high frequencies. 
% We will here similarly see that an image can be smoothed by applying a low-pass filters to the rows and the columns.
% Let us consider such computational molecules. In particular, let us as before take filter coefficients taken from Pascal's triangle. 
% If we use the filter $S=\frac{1}{4}\{1,\underline{2},1\}$ (row 2 from Pascal's triangle), 
% Theorem 9.9 gives the computational molecule
% 
% $$A= \frac{1}{16} \begin{pmatrix} 1 & 2 & 1\\  2 & \underline{4} & 2\\  1 & 2 & 1 \end{pmatrix}.$$
% If the pixels in the image are $p_{i,j}$, this means that we compute the new pixels by
% 
% $$\begin{multline*} \hat p_{i,j}=\frac{1}{16}\bigl(4 p_{i,j}+2(p_{i,j-1}+p_{i-1,j}+p_{i+1,j}+p_{i,j+1})\\  +p_{i-1,j-1}+p_{i+1,j-1}+p_{i-1,j+1}+p_{i+1,j+1}\bigr). \end{multline*}$$
% If we instead use the filter $S=\frac{1}{64}\{1,6,15,\underline{20},15,6,1\}$ (row 6 from Pascal's triangle), we get the computational molecule
% 
% 
% $$\frac{1}{4096} \begin{pmatrix} 1 & 6 & 15 & 20 & 15 & 6 & 1 \\   6 & 36 & 90 & 120 & 90 & 36 & 6 \\   15 & 90 & 225 & 300 & 225 & 90 & 15 \\   20 & 120 & 300 & \underline{400} & 300 & 120 & 20 \\   15 & 90 & 225 & 300 & 225 & 90 & 15 \\   6 & 36 & 90 & 120 & 90 & 36 & 6 \\   1 & 6 & 15 & 20 & 15 & 6 & 1 \end{pmatrix}.$$
% We anticipate that both molecules give a smoothing effect, but that the second molecules provides more smoothing. 
% The result of applying the two molecules in (9.6) and (9.7) to our greyscale-image 
% 
X1=CreateExcerpt();
imshow(uint8(X1))
% is shown below
% 
shortmolecule = @(x) filterS([1 2 1]/4, x, 1);
X2 = tensor_impl(X1, shortmolecule, shortmolecule);
imshow(uint8(X2))
% 
S2 = conv([1 3 3 1]/8, [1 3 3 1]/8);
longmolecule = @(x) filterS(S2, x, 1);
X3 = tensor_impl(X1, longmolecule, longmolecule);
imshow(uint8(X3))
% To make the smoothing effect visible, we have zoomed in on the face in the image
% . The smoothing effect is best visible in the second image.
% Smoothing effects are perhaps more visible if we use the following simple image.
% 
N=16;
img=255*ones(N);
img(1:(N/2),1:(N/2)) = 0;
img((N/2+1):N,(N/2+1):N) = 0;
X1 = img;
imshow(uint8(X1))
% The image can be smoothed out horizontally as follows,
% 
donothing = @(x) x;
X2 = img;
X2 = tensor_impl(X2, donothing, shortmolecule);
X2 = mapto01(X2);
X2 = X2*255;
imshow(uint8(X2))
% vertically as follows,
% 
X3 = img;
X3 = tensor_impl(X3, shortmolecule, donothing);
X3 = mapto01(X3);
X3 = X3*255;
imshow(uint8(X3))
% and both horizontally and vertically as follows.
% 
X4 = img;
X4 = tensor_impl(X4, shortmolecule, donothing);
X4 = tensor_impl(X4, donothing, shortmolecule);
X4 = mapto01(X4);
X4 = X4*255;
imshow(uint8(X4))
% Again we have used the filter $S=\frac{1}{4}\{1,\underline{2},1\}$. 
% Here we also have shown what happens if we only smooth the image in one of the directions.
% In the right image we have smoothed in both directions. 
% We clearly see the union of the two one-dimensional smoothing operations then.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 9.10: Edge detection
% 
% % keywords = ipynbimageschap; mimageschap
% 
% Another operation on images which can be expressed in terms of computational molecules is _edge detection_. 
% An edge in an image is characterised by a large change in intensity values over a small distance in the image. 
% For a continuous function this corresponds to a large derivative. 
% An image is only defined at isolated points, so we cannot compute derivatives, 
% but since a grey-level image is a scalar function of two variables, 
% we have a perfect situation for applying numerical differentiation techniques.
% 
% 
%% Partial derivative in $x$-direction
% Let us first consider computation of the partial derivative $\partial P/\partial x$ at all points in the image.
% Note first that it is the second coordinate in an image which refers to the $x$-direction used when plotting functions. This means that the familiar symmetric Newton quotient
% approximation for the partial derivative [1] takes the form
% 
% $$\frac{\partial P}{\partial x}(i,j)\approx \frac{p_{i,j+1}-p_{i,j-1}}{2},$$
% where we have used the convention $h=1$ which means that the derivative is measured in terms of 'intensity per pixel'. 
% This corresponds to applying the bass-reducing filter $S=\frac{1}{2}\{1,\underline{0},-1\}$ to all the rows 
% (alternatively, applying the tensor product $I\otimes S$ to the image). 
% We can thus express this in terms of the computational molecule.
% 
% $$\begin{equation*} \frac{1}{2} \begin{pmatrix} 0 & 0 & 0\\  1 & \underline{0} & -1\\  0 & 0 & 0 \end{pmatrix}. \end{equation*}$$
% 
% We have included the two rows of 0s just to make it clear how the computational molecule is to be interpreted when we place it over the pixels. 
% Let us first apply this molecule to the usual excerpt of the Lena image. 
% 
diffxmolecule = @(x) filterS([1 0 -1]/2, x, 0);
excerpt=CreateExcerpt();
X1=tensor_impl(excerpt, donothing, diffxmolecule);
imshow(uint8(X1))
% This images shows many artefacts since the pixel values lie outside the legal range: 
% many of the intensities are in fact negative. 
% More specifically, the intensities turn out to vary in the interval $[-0.424,0.418]$. 
% Let us therefore normalise and map all intensities to $[0,1]$.
% 
X2 = tensor_impl(excerpt, donothing, diffxmolecule);
X2 = mapto01(X2);
X2 = X2*255;
imshow(uint8(X2))
% The predominant color of this image is an average grey, i.e. an intensity of about 0.5. 
% To get more detail in the image we therefore try to increase the contrast by applying the function $f_{50}$ in equation (9.1) to each intensity value. 
% 
X3 = tensor_impl(excerpt, donothing, diffxmolecule);
X3 = mapto01(X3);
X3 = X3*255;
X3 = contrastadjust0(X3,50);
imshow(uint8(X3))
% This does indeed show more detail.
% It is important to understand the colors in these images. We have computed the derivative in the $x$-direction, and we recall that the computed values varied 
% in the interval $[-0.424,0.418]$. The negative value corresponds to the largest average decrease in intensity from a pixel $p_{i-1,j}$ to a pixel $p_{i+1,j}$. 
% The positive value on the other hand corresponds to the largest average increase in intensity. A value of 0 in the left image in Figure 9.13 corresponds to no change in 
% intensity between the two pixels.
% 
% When the values are mapped to the interval $[0,1]$ in the second image, the small values are mapped to something close to 0 (almost black), 
% the maximal values are mapped to something close to 1 (almost white), and the values near 0 are mapped to something close to $0.5$ (grey). 
% In the third image these values have just been emphasised even more.
% 
% The third image tells us that in large parts of the image there is very little variation in the intensity. However, there are some small areas where 
% the intensity changes quite abruptly, and if you look carefully you will notice that in these areas there is typically both black and white pixels close together, 
% like down the vertical front corner of the bus. This will happen when there is a stripe of bright or dark pixels that cut through an area of otherwise quite 
% uniform intensity.
% 
% 
%% Partial derivative in $y$-direction
% The partial derivative $\partial P/\partial y$ can be computed analogously to $\partial P/\partial x$,
% i.e. we apply the filter $-S=\frac{1}{2}\{-1,\underline{0},1\}$ to all columns of the image (alternatively, apply the tensor product $-S\otimes I$ to the image), where $S$ is the filter which we
% used for edge detection in the $x$-direction. Note that the positive direction of this axis in an image is opposite to the direction of the $y$-axis we use when plotting functions.
% We can express this in terms of the computational molecule
% 
% $$\begin{equation*} \frac{1}{2} \begin{pmatrix} 0 & 1 & 0\\  0 & \underline{0} & 0\\  0 & -1 & 0 \end{pmatrix}. \end{equation*}$$
% Let us compare the partial derivatives in both directions.
% First we reproduce the partial derivative in $x$-direction.
% 
excerpt=CreateExcerpt();
X1 = tensor_impl(excerpt, donothing, diffxmolecule);
X1 = mapto01(X1);
X1 = X1*255;
X1 = contrastadjust0(X1, 50);
imshow(uint8(X1))
% Then we compute the partial derivative in $y$-direction.
% 
diffymolecule = @(x) filterS([-1 0 1]/2, x, 0);
X2 = tensor_impl(excerpt, diffymolecule, donothing);
X2 = mapto01(X2);
X2 = X2*255;
X2 = contrastadjust0(X2, 50);
imshow(uint8(X2))
% The intensities have been normalised and the contrast enhanced by the function $f_{50}$ from Equation (9.1).
% 
%% The gradient
% The gradient of a scalar function is often used as a measure of the size of the first derivative. The gradient is defined by the vector
% $$
% \nabla P=\Biggl(\frac{\partial P}{\partial x}, \frac{\partial P}{\partial y}\Biggr),
% $$
% so its length is given by
% $$
% \vert \nabla P\vert = \sqrt{\Biggl(\frac{\partial P}{\partial x}\Biggr)^2 + \Biggl(\frac{\partial P}{\partial y}\Biggr)^2}.
% $$
% When the two first derivatives have been computed it is a simple matter to compute the gradient vector and its length. 
% Note that, as for the first order derivatives, it is possible for the length of the gradient to be outside the legal range of values. 
% Let us first compute and show the values of the gradient.
% 
excerpt=CreateExcerpt();
X1 = tensor_impl(excerpt, donothing, diffxmolecule);
X2 = tensor_impl(excerpt, diffymolecule, donothing);
X1 = X1.^2 + X2.^2;
imshow(uint8(X1))
% Then we map the intensities to the legal range.
% 
X2 = X1;
X2 = mapto01(X2);
X2 = X2*255;
imshow(uint8(X2))
% Then we use the function $g_{0.01}$ to adjust the contrast.
% 
X3 = X2;
X3 = contrastadjust0(X3, 50);
imshow(uint8(X3))
% The image of the gradient looks quite different from the images of the two partial derivatives. The reason is that the numbers that represent the length of 
% the gradient are (square roots of) sums of squares of numbers.  This means that the parts of the image that have virtually constant intensity 
% (partial derivatives close to 0) are colored black. In the images of the partial derivatives these values ended up in the middle of the range of intensity values, 
% with a final color of grey, since there were both positive and negative values.
% To enhance the contrast for this image we should thus do something different from what was done in the other images, since we now have a large number of intensities near 0. 
% The solution was to apply a function like the ones shown in the right plot in Figure 9.8. Here we have used the function $g_{0.01}$.
% 
% Figure 9.14 shows the two first-order partial derivatives and the gradient. If we compare the two partial derivatives we see that the 
% $x$-derivative seems to emphasise vertical edges while the $y$-derivative seems to emphasise horizontal edges. This is precisely what we must expect. 
% The $x$-derivative is large when the difference between neighbouring pixels in the $x$-direction is large, which is the case across a vertical edge. 
% The $y$-derivative enhances horizontal edges for a similar reason.
% 
% The gradient contains information about both derivatives and therefore emphasises edges in all directions. It also gives a simpler image since the sign of 
% the derivatives has been removed.
% 
% 
% 
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 9.11: Second-order derivatives
% 
% % keywords = ipynbimageschap; mimageschap
% 
% To compute the three second order derivatives we can combine the two computational molecules which we already have described. For the mixed second order derivative we get
% $(I\otimes S)((-S)\otimes I)=-S\otimes S$. For the last two second order derivative $\frac{\partial^2P}{\partial x^2}$, $\frac{\partial^2P}{\partial y^2}$, we can also use the
% three point approximation to the second derivative [1]
% 
% $$\frac{\partial P}{\partial x^2}(i,j)\approx p_{i,j+1}-2p_{i,j}+p_{i,j-1}$$
% (again we have set $h=1$). This gives a smaller molecule than if we combine the two molecules for order one differentiation 
% (i.e. $(I\otimes S)(I\otimes S)=(I\otimes S^2)$ and $((-S)\otimes I)((-S)\otimes I)=(S^2\otimes I)$), since 
% $S^2=\frac{1}{2}\{1,\underline{0},-1\}\frac{1}{2}\{1,\underline{0},-1\}=\frac{1}{4}\{1,0,\underline{-2},0,1\}$. 
% The second order derivatives of an image $P$ can thus be computed by applying the computational molecules
% 
% $$\frac{\partial^2P}{\partial x^2}:&\qquad \begin{pmatrix} 0 & 0 & 0\\  1 & \underline{-2} & 1\\  0 & 0 & 0 \end{pmatrix},     \\  \frac{\partial^2P}{\partial y \partial x}:&\qquad \frac{1}{4} \begin{pmatrix} -1 & 0 & 1\\  0 & \underline{0} & 0\\  1 & 0 & -1 \end{pmatrix},    \\  \frac{\partial^2P}{\partial y^2}:&\qquad \begin{pmatrix} 0 & 1 & 0\\  0 & \underline{-2} & 0\\  0 & 1 & 0 \end{pmatrix}.   \end{align}$$
% With these molecules it is quite easy to compute the second-order derivatives.
% First the $xx$-derivative.
% 
excerpt=CreateExcerpt();
X1 = tensor_impl(excerpt, donothing, diffxmolecule);
X1 = tensor_impl(X1, donothing, diffxmolecule);
X1 = mapto01(X1);
X1 = X1*255;
X1 = contrastadjust0(X1, 100);
imshow(uint8(X1))
% Then the $xy$-derivative.
% 
X2 = tensor_impl(excerpt, donothing, diffxmolecule);
X2 = tensor_impl(X2, diffymolecule, donothing);
X2 = mapto01(X2);
X2 = X2*255;
X2 = contrastadjust0(X2, 100);
imshow(uint8(X2))
% Then the $yy$-derivative.
% 
X3 = tensor_impl(excerpt, diffymolecule, donothing);
X3 = tensor_impl(X3, diffymolecule, donothing);
X3 = mapto01(X3);
X3 = X3*255;
X3 = contrastadjust0(X3, 100);
imshow(uint8(X3))
% The computed derivatives were first normalised and then the contrast enhanced with the function $f_{100}$ in each image, see equation (9.1).
% 
% As for the first derivatives, the $xx$-derivative seems to emphasise vertical edges and the $yy$-derivative horizontal edges. However, we also see that the second derivatives are
% more sensitive to noise in the image (the areas of grey are less uniform). The mixed derivative behaves a bit differently from the other two, and not surprisingly it seems to pick
% up both horizontal and vertical edges.
% 
% This procedure can be generalized to higher order derivatives also. To apply $\frac{\partial^{k+l}P}{\partial x^k\partial y^l}$ to an image we can compute $S_l\otimes S_k$ where
% $S_r$ corresponds to any point method for computing the $r$'th order derivative. We can also compute $(S^l)\otimes(S^k)$, where we iterate the
% filter $S=\frac{1}{2}\{1,\underline{0},-1\}$ for the first derivative, but this gives longer filters.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 9.12: Chess pattern image
% 
% % keywords = ipynbimageschap; mimageschap
% 
% Let us apply the molecules for differentiation to a chess pattern test image. 
% The chess pattern image is produced as follows.
% 
N = 128;
imgstr=zeros(N);
for x=0:(N-1)
for y=0:(N-1)
imgstr(x+1,y+1) = 255*( (mod(x-1,64)>=32) == (mod(y-1,64)>=32) );
end
end
X11 = imgstr;
imshow(uint8(X11))
% Then we apply $S\otimes I$.
% 
X12 = imgstr;
X12 = tensor_impl(X12, diffymolecule, donothing);
X12 = mapto01(X12);
X12 = X12*255;
imshow(uint8(X12))
% Then we apply $I\otimes S$.
% 
X13 = imgstr;
X13 = tensor_impl(X13, donothing, diffxmolecule);
X13 = mapto01(X13);
X13 = X13*255;
imshow(uint8(X13))
% Then we apply $S\otimes S$.
% 
X21 = imgstr;
X21 = tensor_impl(X21, donothing, diffxmolecule);
X21 = tensor_impl(X21, diffymolecule, donothing);
X21 = mapto01(X21);
X21 = X21*255;
imshow(uint8(X21))
% Then we apply $I\otimes S^2$.
% 
X22 = imgstr;
X22 = tensor_impl(X22, donothing, diffxmolecule);
X22 = tensor_impl(X22, donothing, diffxmolecule);
X22 = mapto01(X22);
X22 = X22*255;
imshow(uint8(X22))
% The we apply $S^2\otimes I$.
% 
X23 = imgstr;
X23 = tensor_impl(X23, diffymolecule, donothing);
X23 = tensor_impl(X23, diffymolecule, donothing);
X23 = mapto01(X23);
X23 = X23*255;
imshow(uint8(X23))
% 
% These images make it is clear that $S\otimes I$ detects all horizontal edges, 
% $I\otimes S$ detects all vertical edges, and that $S\otimes S$ detects all points where abrupt changes appear in both directions. 
% We also see that the second order partial derivative detects exactly the same edges which the first order
% partial derivative found. Note that the edges detected with $I\otimes S^2$ are wider than the ones detected with $I\otimes S$. The reason is that the filter $S^2$
% has more filter coefficients than $S$. Also, edges are detected with different colors. This reflects whether the difference between the neighbouring pixels
% is positive or negative. The values after we have applied the tensor product may thus not lie in the legal range of pixel values (since they may be negative). The figures
% have taken this into account by mapping the values back to a legal range of values, as we did in Chapter 9. Finally, we also see additional edges
% at the first and last rows/edges in the images. The reason is that the filter $S$ is defined by assuming that the pixels repeat periodically (i.e. it is a
% circulant Toeplitz matrix). Due to this, we have additional edges at the first/last rows/edges. 
% This effect can also be seen in Chapter 9, although
% there we did not assume that the pixels repeat periodically.
% 
% Defining a two-dimensional filter by filtering columns and then rows is not the only way we can define a two-dimensional filter. Another possible way is to let the $MN\times
% MN$-matrix itself be a filter. Unfortunately, this is a bad way to define filtering of an image, since there are some undesirable effects near the boundaries between rows:
% in the vector we form, the last element of one row is followed by the first element of the next row. These boundary effects are unfortunate when a filter is applied.
% 
% 
% 
% 
% 
% 
% 
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Exercise 9.14: Generate images
% 
% % keywords = ipynbimageschap; mimageschap
% 
% Write code which calls the function |tensor_impl| with appropriate filters and which generate the following images:
% 
% 
% *a)* The right image in Figure 9.11.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The following code can be used
% 
S2 = conv([1 3 3 1]/8, [1 3 3 1]/8);
longmolecule = @(x) filterS(S2, x, 1);
X3 = tensor_impl(X1, longmolecule, longmolecule);
imshow(uint8(X3))
% 
% % --- end solution of exercise ---
% 
% *b)* The right image in Figure 9.13.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The following code can be used
% 
X3 = tensor_impl(excerpt, donothing, diffxmolecule);
X3 = mapto01(X3);
X3 = X3*255;
X3 = contrastadjust0(X3,50);
imshow(uint8(X3))
% 
% % --- end solution of exercise ---
% 
% *c)* The images in Figure 9.15.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The following code can be used
% 
excerpt=CreateExcerpt();
X1 = tensor_impl(excerpt, donothing, diffxmolecule);
X2 = tensor_impl(excerpt, diffymolecule, donothing);
X1 = X1.^2 + X2.^2;
imshow(uint8(X1))
% 
X2 = X1;
X2 = mapto01(X2);
X2 = X2*255;
imshow(uint8(X2))
% 
X3 = X2;
X3 = contrastadjust0(X3, 50);
imshow(uint8(X3))
% 
% % --- end solution of exercise ---
% 
% *d)* The images in Figure 9.14.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The following code can be used
% 
excerpt=CreateExcerpt();
X1 = tensor_impl(excerpt, donothing, diffxmolecule);
X1 = mapto01(X1);
X1 = X1*255;
X1 = contrastadjust0(X1, 50);
imshow(uint8(X1))
% 
diffymolecule = @(x) filterS([-1 0 1]/2, x, 0);
X2 = tensor_impl(excerpt, diffymolecule, donothing);
X2 = mapto01(X2);
X2 = X2*255;
X2 = contrastadjust0(X2, 50);
imshow(uint8(X2))
% 
% % --- end solution of exercise ---
% 
% *e)* The images in Figure 9.16.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The following code can be used
% 
excerpt=CreateExcerpt();
X1 = tensor_impl(excerpt, donothing, diffxmolecule);
X1 = tensor_impl(X1, donothing, diffxmolecule);
X1 = mapto01(X1);
X1 = X1*255;
X1 = contrastadjust0(X1, 100);
imshow(uint8(X1))
% 
X2 = tensor_impl(excerpt, donothing, diffxmolecule);
X2 = tensor_impl(X2, diffymolecule, donothing);
X2 = mapto01(X2);
X2 = X2*255;
X2 = contrastadjust0(X2, 100);
imshow(uint8(X2))
% 
X3 = tensor_impl(excerpt, diffymolecule, donothing);
X3 = tensor_impl(X3, diffymolecule, donothing);
X3 = mapto01(X3);
X3 = X3*255;
X3 = contrastadjust0(X3, 100);
imshow(uint8(X3))
% 
% % --- end solution of exercise ---
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Exercise 9.27: Implement two-dimensional FFT and DCT
% 
% % keywords = ipynbimageschap; mimageschap
% 
% Implement functions |DFTImplFull|, |IDFTImplFull|, |DCTImplFull|, and |IDCTImplFull| which applies the DFT, IDFT, DCT, and IDCT, to the entire vector, and use these to implement. 
% FFT2, IFFT2, DCT2, and IDCT2 on an image, with the help of the function |tensor_impl|.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The following code can be used.
% 
DFTImplFull = @(x) fft(x(:,:));
IDFTImplFull = @(x) ifft(x(:,:));
% 
DCTImplFull = @(x) dct(x(:,:));
IDCTImplFull = @(x) idct(x(:,:));
% 
X = tensor_impl(img, DFTImplFull, DFTImplFull);
X = tensor_impl(img, IDFTImplFull, IDFTImplFull);
X = tensor_impl(img, DCTImplFull, DCTImplFull);
X = tensor_impl(img, IDCTImplFull, IDCTImplFull);
% 
% % --- end solution of exercise ---
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Exercise 9.28: Zeroing out DCT coefficients
% 
% % keywords = ipynbimageschap; mimageschap
% 
% The function |forw_comp_rev_DCT2| in the module |forw_comp_rev| applies the DCT to a part of our sample image, 
% and sets DCT coefficients below a certain threshold to be zero.
% This is very similar to what the JPEG standard does.
% 
% 
% Run |forw_comp_rev_DCT2| for different threshold parameters, 
% and with the functions |DCTImpl8|, |IDCTImpl8|, |DCTImplFull|, and |IDCTImplFull| as parameters. 
% Check that this reproduces the DCT test images of this section, 
% and that the correct numbers of values which have been neglected (i.e. which are below the threshold) are printed on screen.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* Once 
% 
X1 = forw_comp_rev_DCT2(@DCTImpl8, @IDCTImpl8, 30);
imshow(uint8(X1))
% 
X2 = forw_comp_rev_DCT2(@DCTImpl8, @IDCTImpl8, 50);
imshow(uint8(X2))
% 
X3 = forw_comp_rev_DCT2(@DCTImpl8, @IDCTImpl8, 100);
imshow(uint8(X3))
% 
X1 = forw_comp_rev_DCT2(DCTImplFull, IDCTImplFull, 30);
imshow(uint8(X1))
% 
X2 = forw_comp_rev_DCT2(DCTImplFull, IDCTImplFull, 50);
imshow(uint8(X2))
% 
X3 = forw_comp_rev_DCT2(DCTImplFull, IDCTImplFull, 100);
imshow(uint8(X3))
% 
% % --- end solution of exercise ---
% 
% % --- end exercise ---
% 
% 
%% Using tensor products to apply wavelets to images
%% Appendix: Basic Linear Algebra
%% Appendix: Signal processing and linear algebra: a translation guide

