LAB 9, 04/02/14 ________________ 1. An image can be imported into MATLAB using the IMREAD function. This produces an array of unsigned integers. Conversion to a numerical matrix with range [0,1] is via the function DOUBLE (for double precision float): I = imread('baseball.png'); % source: scroll to bottom whos I help uint8 I = double(I)/255; imshow(I); size(I) In this lab, we study the representation of a NxN image (where N is a power of 2) using Haar wavelets introduced in Lab 7. 2. In Lab 7, we generated a NxN Haar wavelet basis matrix V using the function NHAAR222. The projection of a signal onto different scales (each scale consisting of a set of columns of NHAAR222) is a form of band-pass filtering: the coarsest scale corresponds to low-pass filtering, the finest scale to high-pass filtering. The same projection can be performed iteratively using a pair of low-pass and high-pass filtering operations. Thus c = V'*s is equivalent to c = Pk*...*P2*P1*s (k = log2(N)) where the matrices P1, P2, ..., Pk are all variants of a matrix P, generated by the function HAARIT: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function P = haarit(N) E = sqrt(1/2)*eye(N/2); P = reshape([E; E; E; -E], N, N); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Assuming haarit.m is in the working directory, type N = 8; P = haarit(N) The top half of P performs low-pass filtering; the bottom half high-pass filtering. In the case N=8, the iteration goes through three steps, with P1 = haarit(8) ; % same as P above P2 = eye(8) ; P2(1:4,1:4) = haarit(4) ; P3 = eye(8) ; P3(1:2,1:2) = haarit(2) ; Thus Pi differs from the NxN identity in the top left 2^(k-i+1) x 2^(k-i+1) block only, which equals the matrix produced by haarit(2^(k-i+1)). 3. To illustrate the iterative Haar wavelet transform, we consider a (N=64)-point sinusoid. On the bar graphs: - blue bars indicate final results, i.e., wavelet coefficients for scales that have already been processed; - red bars indicate coefficients that will undergo further iteration. Expanded to full length N, the red bars are the projection of s onto scales yet to be processed. N = 64; t = 1:N ; s = sin(5*t/N).' ; plot(t,s), grid; pause bar(t,s,'r') pause c = s ; m = N; for i = 0:log2(N)-1 ; P = haarit(m) c(1:m) = P*c(1:m) ; clf hold on bar(1:m/2,c(1:m/2),'r'), axis([0 N+1 -5 5]); bar(m/2+1:N,c(m/2+1:N),'b'); hold off pause m=m/2; end hold on bar(1:N,c,'b'), axis([0 N+1 -5 5]); 4. The Haar wavelet transform can be also applied to images, i.e., any image can be expressed as a linear combination of orthonormal basis images which are generated using the Haar wavelet. In a nutshell: If v_i and v_j are the ith and jth columns of an orthonormal NxN matrix V, then % A_ij = (v_i)*(v_j).' is also a NxN matrix. There are N^2 such matrices A_ij. It is not difficult to show that if each of these matrices is reshaped into a single (N^2)x1 column vector, the resulting vectors are orthogonal, each with unit norm. Thus any NxN matrix X can be expressed as a linear combination of the matrices A_ij: % X = sum (c_ij)*(A_ij) = sum (c_ij)*(v_i)*(v_j).' % (sum over i,j) If C is the NxN matrix of coefficients c_ij, the above sum can be written compactly in matrix form as % X = V*C*V' (synthesis equation) C is easily obtained from X using the orthonormality of V, namely V'*V = V*V' = I: % C = V'*X*V (analysis equation) 5. In the case of the Haar wavelet transform, we have the factorization V' = Pk*...*P2*P1 The analysis equation C = V'*X*V then becomes C = Pk*...*(P1*(P1*X*P1')*P2')*...*Pk' i.e, each factor P is applied both as a column and as a row transform (the latter amounting to right-multiplication by P'). We illustrate this on the given image, by performing the first four row-column transform iterations: clf N = 512; C = I; m = N ; for i=1:4 % max. # iterations = k P = haarit(m); C(1:m,1:m) = C(1:m,1:m)*P'; imshow(C,[]); pause C(1:m,1:m) = P*C(1:m,1:m); imshow(C,[]); pause m=m/2; end 6. The 32x32 low-frequency version of the image is stored in the top left corner of the coefficient matrix C. The remaining entries of C are higher-frequency wavelet coefficients. For better contrast, we repeat using two IMSHOW displays per stage, for the low and high-frequency portions of C. Recall that the option [] in IMSHOW maps min -> 0 and max -> 1. clf N = 512; C = I; mask = zeros(N,N) ; m = N ; for i=1:4 P = haarit(m); C(1:m,1:m) = C(1:m,1:m)*P'; mask(1:N, m/2+1:N) = 1; imshow(C.*(1-mask),[]); pause imshow(abs(C).*mask,[]); pause C(1:m,1:m) = P*C(1:m,1:m); mask(m/2+1:N, 1:N) = 1; imshow(C.*(1-mask),[]); pause imshow(abs(C).*mask,[]); pause m=m/2; end clf imshow(C(1:m, 1:m), []) 7. Reconstruction of the image from its coefficients involves the synthesis equation X = V*C*V', or in iterative form: % X = Pk'*...*(P2'*(P1'*X*P1)*P2)*...*Pk We illustrate this by (effectively) reversing the iteration in 5. above. for i=1:4 m = 2*m; P = haarit(m); C(1:m,1:m) = P'*C(1:m,1:m) ; imshow(C,[]); pause C(1:m,1:m) = C(1:m,1:m)*P ; imshow(C,[]); pause end 8. We now apply the wavelet transform to a noisy version of the given image, generated as follows: I_clean = I ; SD = 20/255 ; % noise standard deviation noise = SD*randn(size(I)) ; % white Gaussian noise I_noisy = I_clean + noise ; I_noisy = max(0, min(I_noisy,1)) ; % limit to [0,1] imshow(I_noisy, []) A clean image has negligible high-frequency wavelet coefficients except where high-contrast edges occur; there, the resulting coefficients can be quite high. White noise, on the other hand, will produce coefficients that also have a white-noise like appearance. These different profiles suggest using a thresholding rule at high frequencies, whereby coefficients that are absolutely less than a certain factor of the noise standard deviation are deemed spurious (i.e., due to noise) and nulled out. We must therefore choose the following parameters: L = lowest scale to which thresholding will apply delta = threshold value (as a multiple of SD) 9. We use the same steps as in 5. to iterate down to scale L: L = 6 ; % thus low-frequency portion spans scales 0-5, % i.e., is 32x32 (as in 5. above) clf N = 512; C = I_noisy; m = N ; while m > 2^(L-1) P = haarit(m); C(1:m,1:m) = C(1:m,1:m)*P'; imshow(C,[]); pause C(1:m,1:m) = P*C(1:m,1:m); imshow(C,[]); pause m=m/2; end Apply thresholding (with a slight twist - can you see it?): delta = 1.0*SD ; C_temp = sign(C).*(abs(C)-delta).*(abs(C) >= delta) ; C_temp(1:m,1:m) = C(1:m,1:m) ; C = C_temp ; Reconstruct: while m < N m = 2*m; P = haarit(m); C(1:m,1:m) = P'*C(1:m,1:m) ; imshow(C,[]); pause C(1:m,1:m) = C(1:m,1:m)*P ; imshow(C,[]); pause end Compare to the noisy original: figure(2); imshow(I_noisy,[]); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Cropped monochrome version of photo by Johannes Gilger, from http://www.flickr.com/photos/heipei/2430910754/sizes/l/in/photostream/