%% 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 = {'chessex': '10.3', 'piecewiseconsttwodim': '10.1', 'ex:exptestimage': '10.10', 'example:spline53exp': '10.6', 'example:pwl2_plot_functions': '10.2', 'showalldwt': '10.5'}
% 
% 
% 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 
%  Using tensor products to apply wavelets to images 
%      Example 10.1: Piecewise constant functions 
%      Example 10.2: Piecewise linear functions 
%      Example 10.3: Applying the Haar wavelet to a very simple example image 
%      Example 10.5: Detail and low-resolution approximations for different wavelets 
%      Example 10.6: The Spline 5/3 wavelet and removing bands in the detail spaces 
%      Exercise 10.10: Experiments on a test image 
%  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
%% Using tensor products to apply wavelets to images
% 
% 
% % --- begin exercise ---
% 
%% Example 10.1: Piecewise constant functions
% 
% % keywords = ipynbtensorwavelet; mtensorwavelet
% 
% If $V_m$ is the vector space of piecewise constant functions on any interval of the form $[k2^{-m},(k+1)2^{-m})$ (as in the piecewise constant
% wavelet), $V_m\otimes V_m$ is the vector space of functions in two variables which are constant on any square of the form 
% $[k_12^{-m},(k_1+1)2^{-m})\times[k_22^{-m},(k_2+1)2^{-m})$. Clearly 
% $\phi_{m,k_1}\otimes\phi_{m,k_2}$ is constant on such a square and $0$ elsewhere, and these functions are a basis for $V_m\otimes V_m$. 
% 
% Let us compute the orthogonal projection of $\phi_{1,k_1}\otimes\phi_{1,k_2}$ onto $V_0\otimes V_0$. Since the Haar wavelet is orthonormal, the basis functions 
% in (10.4) are orthonormal, and we can thus use the orthogonal decomposition formula to find this projection.  
% Clearly $\phi_{1,k_1}\otimes\phi_{1,k_2}$ has different support from all except
% one of $\phi_{0,n_1}\otimes\phi_{0,n_2}$. Since
% 
% $$\begin{equation*} \langle \phi_{1,k_1}\otimes\phi_{1,k_2},\phi_{0,n_1}\otimes\phi_{0,n_2}\rangle = \langle \phi_{1,k_1},\phi_{0,n_1}\rangle \langle \phi_{1,k_2},\phi_{0,n_2}\rangle = \frac{\sqrt{2}}{2}\frac{\sqrt{2}}{2} = \frac{1}{2} \end{equation*}$$
% when the supports intersect, we obtain
% 
% $$\begin{equation*}  \text{proj}_{V_0\otimes V_0} (\phi_{1,k_1}\otimes\phi_{1,k_2})  =  \begin{cases}  \frac{1}{2}(\phi_{0,k_1/2}\otimes\phi_{0,k_2/2})         & \text{when $k_1,k_2$ are even} \\   \frac{1}{2}(\phi_{0,k_1/2}\otimes\phi_{0,(k_2-1)/2})     & \text{when $k_1$ is even, $k_2$ is odd} \\   \frac{1}{2}(\phi_{0,(k_1-1)/2}\otimes\phi_{0,k_2/2})     & \text{when $k_1$ is odd, $k_2$ is even} \\   \frac{1}{2}(\phi_{0,(k_1-1)/2}\otimes\phi_{0,(k_2-1)/2}) & \text{when $k_1,k_2$ are odd}  \end{cases} \end{equation*}$$
% So, in this case there were $4$ different formulas, since there were $4$ different combinations of even/odd. Let us also compute the projection onto the orthogonal
% complement of $V_0\otimes V_0$ in $V_1\otimes V_1$, and let us express this in terms of the $\phi_{0,n},\psi_{0,n}$, like we did in the one-variable case. Also
% here there are $4$ different formulas. When $k_1,k_2$ are both even we obtain
% 
% $$\begin{align*}    \lefteqn{\phi_{1,k_1}\otimes\phi_{1,k_2} - \text{proj}_{V_0\otimes V_0}(\phi_{1,k_1}\otimes\phi_{1,k_2})} \\          &= \phi_{1,k_1}\otimes\phi_{1,k_2} - \frac{1}{2}(\phi_{0,k_1/2}\otimes\phi_{0,k_2/2}) \\          &= \left(\frac{1}{\sqrt{2}}(\phi_{0,k_1/2}+\psi_{0,k_1/2})\right)\otimes\left(\frac{1}{\sqrt{2}}(\phi_{0,k_2/2}+\psi_{0,k_2/2})\right)             - \frac{1}{2}(\phi_{0,k_1/2}\otimes\phi_{0,k_2/2}) \\          &= \frac{1}{2}(\phi_{0,k_1/2}\otimes\phi_{0,k_2/2}) + \frac{1}{2}(\phi_{0,k_1/2}\otimes\psi_{0,k_2/2}) \\          &  + \frac{1}{2}(\psi_{0,k_1/2}\otimes\phi_{0,k_2/2})             + \frac{1}{2}(\psi_{0,k_1/2}\otimes\psi_{0,k_2/2}) - \frac{1}{2}(\phi_{0,k_1/2}\otimes\phi_{0,k_2/2}) \\          &= \frac{1}{2}(\phi_{0,k_1/2}\otimes\psi_{0,k_2/2}) + \frac{1}{2}(\psi_{0,k_1/2}\otimes\phi_{0,k_2/2}) + \frac{1}{2}(\psi_{0,k_1/2}\otimes\psi_{0,k_2/2}). \end{align*}$$
% Here we have used the relation $\phi_{1,k_i}=\frac{1}{\sqrt{2}}(\phi_{0,k_i/2}+\psi_{0,k_i/2})$, which we have from our first analysis of the Haar wavelet. 
% Checking the other possibilities we find similar formulas for the projection onto the orthogonal complement of $V_0\otimes V_0$ in $V_1\otimes V_1$ when either
% $k_1$ or $k_2$ is odd. In all cases, the formulas use the basis functions for $W_0^{(0,1)}$, $W_0^{(1,0)}$, $W_0^{(1,1)}$. 
% Let us plot these. First we define $\phi$ and $\psi$.
% 
t1 = linspace(-0.4, 1.4, 300);
t2 = linspace(-0.4, 1.4, 300);
[x, y] = meshgrid(t1,t2);
phidef = @(t) (t >= 0).*(t <= 1);
psidef = @(t) (t >= 0).*(t <= 0.5) - (t >= 0.5).*(t <= 1);
% Then $\phi\otimes\phi$.
% 
fig = figure();
mesh(x, y, phidef(x).*phidef(y));
% Then $\phi\otimes\psi$.
% 
figure();
mesh(x, y, phidef(x).*psidef(y));
% Then $\psi\otimes\phi$.
% 
figure();
mesh(x, y, psidef(x).*phidef(y));
% Then $\psi\otimes\psi$.
% 
figure();
mesh(x, y, psidef(x).*psidef(y));
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 10.2: Piecewise linear functions
% 
% % keywords = ipynbtensorwavelet; mtensorwavelet
% 
% If we instead use any of the wavelets for piecewise linear functions, the wavelet basis functions are not orthogonal anymore, just as in the one-dimensional case.
% Let us also plot these. First we define the new $\phi$ and $\psi$.
% 
t1 = linspace(-1, 2, 300);
t2 = linspace(-1, 2, 300);
[x, y] = meshgrid(t1, t2);
phidef = @(t) (1-abs(t)).*(abs(t)<=1);
psidef = @(t) (1-2*abs(t-0.5)).*(t>=0).*(t<=1) - 0.25*phidef(t) - 0.25*phidef(t-1);
% Then $\phi\otimes\phi$.
% 
figure();
mesh(x, y, phidef(x).*phidef(y));
% Then $\phi\otimes\psi$.
% 
figure();
mesh(x, y, phidef(x).*psidef(y));
% Then $\psi\otimes\phi$.
% 
figure();
mesh(x, y, psidef(x).*phidef(y));
% Then $\psi\otimes\psi$.
% 
figure();
mesh(x, y, psidef(x).*psidef(y));
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 10.3: Applying the Haar wavelet to a very simple example image
% 
% % keywords = ipynbtensorwavelet; mtensorwavelet
% 
% Let us apply the Haar wavelet to the sample chess pattern example image from Figure 9.17. The lowpass filter of the Haar wavelet was
% essentially a smoothing filter with two elements. Also, as we have seen, the highpass filter essentially computes an approximation to the partial derivative.
% Clearly, abrupt changes in the vertical and horizontal directions appear here only at the edges in the chess pattern, and abrupt changes in both directions appear
% only at the grid points in the chess pattern. Due to Observation 10.9, after a DWT2 we expect to see the following:
% 
% * In the upper left corner, we should see a low-resolution version of the image.
% 
% * In the upper right corner, only the vertical edges in the chess pattern should be visible.
% 
% * In the lower left corner, only the horizontal edges in the chess pattern should be visible.
% 
% * In the lower right corner, only the grid points in the chess pattern should be visible.
% 
% These effects are seen clearly if we apply one level of the DWT2 to the chess pattern example image. 
% This can be achieved with the following code.
% 
img=zeros(128);
for x=0:127
for y=0:127
img(x+1,y+1)=255*( (mod(x+1,64)>=32) == (mod(y+1,64)>=32) );
end
end
img = DWT2Impl(img,1,'Haar');
img = mapto01(img);
img = img*255;
imshow(uint8(img))
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 10.5: Detail and low-resolution approximations for different wavelets
% 
% % keywords = ipynbtensorwavelet; mtensorwavelet
% 
% Let us take a closer look at the images generated when we use different wavelets. 
% Above we viewed the low-resolution approximation as a smaller image. 
% Let us compare with the image resulting from setting the wavelet detail coefficients to zero, and viewing the result as an image of the same size. 
% In particular, let us neglect the wavelet coefficients as pictured in Figure 10.4 and Figure 10.5.
% We should expect that the lower order resolution approximations from $V_0$ are worse when $m$ increase.
% 
% We will use the function |forw_comp_rev_DWT2| to generate the first four low-resolution approximations or detail for the usual excerpt of the sample image.
% The low-level resolution approximations look as follows when we use the Haar wavelet.
% 
X1 = forw_comp_rev_DWT2(1, 'Haar');
imshow(uint8(X1))
% 
X2 = forw_comp_rev_DWT2(2, 'Haar');
imshow(uint8(X2))
% 
X3 = forw_comp_rev_DWT2(3, 'Haar');
imshow(uint8(X3))
% 
X4 = forw_comp_rev_DWT2(4, 'Haar');
imshow(uint8(X4))
% The details look as follows when we use the Haar wavelet.
% 
X1 = forw_comp_rev_DWT2(1, 'Haar', 0);
imshow(uint8(X1))
% 
X2 = forw_comp_rev_DWT2(2, 'Haar', 0);
imshow(uint8(X2))
% 
X3 = forw_comp_rev_DWT2(3, 'Haar', 0);
imshow(uint8(X3))
% 
X4 = forw_comp_rev_DWT2(4, 'Haar', 0);
imshow(uint8(X4))
% The low-level resolution approximations look as follows when we use the CDF 9/7 wavelet.
% 
X1 = forw_comp_rev_DWT2(1, 'cdf97');
imshow(uint8(X1))
% 
X2 = forw_comp_rev_DWT2(2, 'cdf97');
imshow(uint8(X2))
% 
X3 = forw_comp_rev_DWT2(3, 'cdf97');
imshow(uint8(X3))
% 
X4 = forw_comp_rev_DWT2(4, 'cdf97');
imshow(uint8(X4))
% The details look as follows when we use the CDF 9/7 wavelet.
% 
X1 = forw_comp_rev_DWT2(1, 'cdf97', 0);
imshow(uint8(X1))
% 
X2 = forw_comp_rev_DWT2(2, 'cdf97', 0);
imshow(uint8(X2))
% 
X3 = forw_comp_rev_DWT2(3, 'cdf97', 0);
imshow(uint8(X3))
% 
X4 = forw_comp_rev_DWT2(4, 'cdf97', 0);
imshow(uint8(X4))
% The black color indicates values which are close to $0$. In other words, most of the coefficients are close to $0$.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 10.6: The Spline 5/3 wavelet and removing bands in the detail spaces
% 
% % keywords = ipynbtensorwavelet; mtensorwavelet
% 
% Since the detail components in images are split into three bands, another thing we can try is to neglect only parts of the detail components (i.e.e some of $W_m^{(1,1)}$, $W_m^{(1,0)}$, 
% $W_m^{(0,1)}$), contrary to the one-dimensional case. Let us use the Spline 5/3 wavelet. 
% 
% We will use the function |mmsubbands| to generate the low-resolution approximations obtained 
% by zeroing out various numbers of the detail bands at level $m$.
% The following code shows the resulting images when the bands on the first level indicated in Figure 10.4 are removed.
% 
[X0,X1,X2] = mmsubbands(1);
% Let us also try this out at level 2.
% 
[X0,X1,X2] = mmsubbands(2);
% The image is seen still to resemble the original one, even after two levels of wavelets coefficients have been neglected. 
% This in itself is good for compression purposes, since we may achieve compression simply by dropping the given coefficients. 
% However, if we continue to neglect more levels of coefficients, the result will look poorer. 
% 
% With the following code we can neglect also the third and fourth levels of detail.
% 
X3 = forw_comp_rev_DWT2(3, 'cdf53');
imshow(uint8(X3))
% Although we still can see details in the image, the quality in the image is definitely poorer.
% Although the quality is poorer when we neglect levels of wavelet coefficients, 
% all information is kept if we additionally include the detail/bands. 
% 
% The following code displays the detail we obtain for DWT's with $1$, $2$, $3$, and $4$ levels.
% 
X1 = forw_comp_rev_DWT2(1, 'Haar', 0);
imshow(uint8(X1))
% Clearly, more detail can be seen in the image when more of the detail is included. 
% 
% 
% 
% 
% 
% 
% 
% As mentioned, the procedure developed in this section for applying a wavelet transform to an image with the help of the tensor product construction, is adopted in
% the JPEG2000 standard. This lossy (can also be used as lossless) image format was developed by the Joint Photographic Experts Group and published in 2000. After significant
% processing of the wavelet coefficients, the final coding with JPEG2000 uses an advanced version of arithmetic coding. At the cost of increased encoding and
% decoding times, JPEG2000 leads to as much as 20 \% improvement in compression ratios for medium compression rates, possibly more for high or low compression rates.
% The artefacts are less visible than in JPEG and appear at higher compression rates. Although a number of components in JPEG2000 are patented, the patent holders
% have agreed that the core software should be available free of charge, and JPEG2000 is part of most Linux distributions. 
% However, there appear to be some further, rather obscure, patents that have not been licensed, and this may be the reason why JPEG2000 is not used more. 
% The extension of JPEG2000 files is \verb/.jp2/.
% 
% 
% 
% 
% 
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Exercise 10.10: Experiments on a test image
% 
% % keywords = ipynbtensorwavelet; mtensorwavelet
% 
% Consider the following image, which is very similar to that in Example 10.3.
% 
img=zeros(128);
for x=0:127
for y=0:127
img(x+1,y+1)=255*( (mod(x,64)>=32) == (mod(y,64)>=32) );
end
end
imshow(uint8(img))
% Then apply the DWT2 with the Haar wavelet. 
% 
img2 = img;
img2 = DWT2Impl(img2,1,'Haar');
img2 = mapto01(img2);
img2 = img2*255;
imshow(uint8(img2))
% You see here, however, that there seems to be no detail components, which is very different from 
% what you saw in Example 10.3, even though the images are very similar. 
% Attempt to explain what causes this to happen.
% 
% % --- begin hint in exercise ---
% 
% *Hint.* Compare with Exercise 5.21.
% 
% % --- end hint in exercise ---
% 
% 
% % --- begin solution of exercise ---
% *Solution.* In Example 10.3, the borders in the chess pattern was chosen so that they occur at odd numbers. 
% This means that the image can not be represented
% exactly in $V_{m-1}\otimes V_{m-1}$, so that there is detail present in the image at all the borders in the chess pattern. 
% In the new image, the borders in the chess pattern was chosen so that they occur at even numbers. This means that the image can be represented
% exactly in $V_{m-1}\otimes V_{m-1}$, so that there is no detail components present in the image.
% 
% % --- end solution of exercise ---
% 
% % --- end exercise ---
% 
% 
%% Appendix: Basic Linear Algebra
%% Appendix: Signal processing and linear algebra: a translation guide

