Download Digital Image Processing - Homework 5 Solutions | ECE 533 and more Assignments Digital Signal Processing in PDF only on Docsity!
Digital Image Processing ECE 533
Solutions to Assignment 5
Department of Electrical and Computing Engineering, University of New Mexico.
Professor Majeed Hayat, [email protected]
April 4, 2008
1 Frequency Domain Design
- To generate and display f use for example: % This uses the RGB values of the image FullLena=sum(double(imread(’Lena.tif’)),3)/(3*255); Lena=FullLena(192:319,192:319); imshow(Lena);
- If h(m, n) is real then H(u, v) is conjugate symmetric, i.e., H(u, v) = H∗(−u, −v). Hence, we can exploit the symmetry of the Fourier transform about the pixel (N/2 = 64 , N/2 = 64) to complete the specification of the filter. Then, it is not hard to realize that H(u, v) = 1, if (u, v) ∈ A, and H(u, v) = 0, if (u, v) ∈/ A, where A = ({ 0 ,... , 31 } × { 0 ,... , 31 }) ∪ ({ 97 ,... , 127 } × { 97 ,... , 127 }) ∪ ({ 0 ,... , 31 } × { 97 ,... , 127 })∪({ 97 ,... , 127 }×{ 0 ,... , 31 }). The symmetry about the middle point imposes the rectangle { 97 ,... , 127 } × { 97 ,... , 127 } from the rectangle { 0 ,... , 31 } × { 0 ,... , 31 }. Note that the size of the new rectangle changes because of the effect of an even number of points. Similarly, the rectangle { 0 ,... , 31 } × { 97 ,... , 127 } imposes the rectangle { 97 ,... , 127 } × { 0 ,... , 31 }. See Figs. 1(a)–(c).
f(m,n)
m
n
|F(u,v)|dB
u
v
|H(u,v)|dB
u
v
|H(u,v)|dB centered
u
v
(a) (b)
h(m,n)
m
n
h(m,n) centered
m
n
g(m,n) via conv
m
n
m
n
g(m,n) via DFT
m
n
|G(u,v)|dB
u
v
(c) (d)
Fig. 1: Processing of Lena, f (m, n), using the filter H(u, v). (a) Lena and the magnitude square, in dB, of its spatial Fourier transform. (b) Magnitude square, in dB, of H(u, v). Spectrum with symmetry in the middle point and centered. (c) Spatial representation of h(m, n) = F −^1 {H(u, v)}. (d) Result of convolving f (m, n) and h(m, n) in spatial (upper-left) and frequency (upper-right) domains. The magnitude square, in dB, of G(u, v).
%######################################### % Filter specification % Pixels [0,31]x[97,127] H(1:32,1:32)=1; % Pixels [0,31]x[0,31] H(1:32,98:128)=1; % Completing the filter: symmetry point (64.5,64.5) H(98:128,98:128)=1; % symmetric part of [0,31]x[97,127]
g=ifft2(fft2(Lenap).*fft2(hp));
- By comparing the images and the spectra in Figs. 1(a) and (d), we see that the filter H is a low-pass filter and the resulting image is blurred.
- In Figs. 2(a)–(d), we see the result of processing the image Lena.tif using the filter H 2 (u, v) = 1 − H(u, v). By comparison, we see that the filter is a high-pass and the resulting image looks sharper than the original image.
- In MATLAB:
%Create C matrix n0=9; C=zeros(N^2,n0^2); Hv=reshape(H,N^2,1); % Very inefficient way to do it for u=0:(N-1) for v=0:(N-1) for m=0:(n0-1) for n=0:(n0-1) C(u+Nv+1,m+n0n+1)=exp(-j2pi(um+vn)/N); end end end end hv=reshape((C’C)(C’)*Hv,n0,n0);
- In Fig. 3(a) we show the result of the mask filter, with size 9 × 9, applied to the image of Lena. Also, in Fig. 3(b) it is shown the restriction of h(m, n) to the first 9 × 9 pixels and the designed mask filter ˜h(m, n). Finally, the last comparison done is in the frequency domain, in the same figure, in which the DFT of h(m, n) is shown as well as the N × N point 2D DFT of ˜h(m, n).
- To see that ˜h(m, n) is the restriction of h(m, n) to a 9 × 9 subset of pixels we will consider two approches. In the first approach, we know that the mask filter ˜h(m, n) is obtained from h˜ = (CH^ C)−^1 CH^ H, where h˜ and H are stacked versions of ˜h(m, n) and h(m, n), respec-
tively, C is an N 2 ×n^20 DFT matrix and H stands for conjugate-transpose (Hermitian operation). In the notes we saw that the columns of C are linearly independent, and moreover, we calculated the inner product between any two columns, say i and j, of C: i) CHi Cj = 0, if i 6 = j; and ii) CHi Cj = N 2 , if i = j. Therefore, CH^ C = N 2 I, and (CH^ C)−^1 = N −^2 I, with I the identity matrix of dimension n^20 × n^20. Then, it is easy to realize that h˜ = N −^2 CH^ H, which corresponds to the first n^20 terms of a 2D inverse DFT of H with N 2 points. Recalling that we have to undo the stacking the claim is proved. In the second approach (as done in class), consider the optimization criteria used to construct the mask filter: ε^2 = || H˜ − H||^22 =
N∑ − 1
i=
N∑ − 1
j=
| H˜(i, j) − H(i, j)|^2
ε^2 ∼
N∑ − 1
i=
N∑ − 1
j=
|˜h(i, j) − h(i, j)|^2 (from Parserval identity)
=
n∑ 0 − 1 i=
n∑ 0 − 1 j=
|˜h(i, j) − h(i, j)|^2 +
N∑ − 1
i=n 0
N∑ − 1
j=n 0
|h(i, j)|^2 (by mask size).
It is therefore seen right-hand side can be minimized by setting its first term to zero, or equivalently, by setting ˜h = h for all (m, n) ∈ { 0 ,... , n 0 − 1 } × { 0 ,... , n 0 − 1 }.
- Zero padding, resulting in a 255 × 255 image, is performed so that circular convo- lution becomes equivalent to (linear) convolution in the range m, n = 0, 1 ,... , 254. This is so because zero padding prevents aliasing that results from circular convo- lution.
2 Frequency Domain Enhancement and Interactive Restora-
tion
- In this problem we used interactive image enhancement in order to eliminate unde- sired frequency components, where undesired means any frequency content that the
3 Wiener Filtering
- It is easy to see that the images follow the binomial distribution with p = 12 , so
P {f (m, n) = r 1 } =
(N 2 − 1
K− 1
(N 2
K
) =^ NK 2 ,
where (nk^ )^ is the binomial coefficient. E [f (m, n)] = r 0 (1 − P {f (m, n) = r 1 }) + r 1 P {f (m, n) = r 1 } = r 0
1 − NK 2
- r (^1) N^ K 2 , μ Also, the second moment is given by E [f (m, n)^2 ]^ = r 02 (1 − P {f (m, n) = r 1 }) + r 12 P {f (m, n) = r 1 } = r 02
1 − NK 2
- Using the definition of conditional probabilities we can calculate the joint pmf. P {f (m, n) = r 1 , f (m′, n′) = r 0 } = P^ {f^ (m, n) =^ r^1 |f^ (m
′, n′) = r 0 } P {f (m′, n′) = r 0 }. But,
P {f (m, n) = r 1 |f (m′, n′) = r 0 } =
(N 2 − 2
K− 1
(N 2 − 1
K
) =^ N 2 K − 1
P {f (m, n) = r 1 |f (m′, n′) = r 1 } =
(N 2 − 2
K− 2
(N 2 − 1
K− 1
) =^ NK 2 − −^1 .
So we obtain P {f (m, n) = r 1 , f (m′, n′) = r 0 } = (^) N 2 K (^) − 1 N^
2 − K
N 2 =^
K
N 2 − 1
1 − NK 2
P {f (m, n) = r 1 , f (m′, n′) = r 1 } = (^) NK 2 − (^) −^1 1 N^ K 2 P {f (m, n) = r 0 , f (m′, n′) = r 0 } = (1 − P {f (m, n) = r 1 |f (m′, n′) = r 0 }) P {f (m′, n′) = r 0 } =
1 − N 2 K − 1
1 − NK 2
In addition, by symmetry we have P {f (m, n) = r 0 , f (m′, n′) = r 1 } = P {f (m, n) = r 1 , f (m′, n′) = r 0 } = (^) N 2 K (^) − 1
1 − NK 2
Therefore we obtain E [f (m, n)f (m′, n′)] =
∑^1
i=
∑^1
j=
rirj P {f (m, n) = ri, f (m′, n′) = rj }
= r^20
1 − N 2 K − 1
1 − NK 2
- 2r 0 r (^1) N 2 K (^) − 1
1 − NK 2
- r^21 N^ K 2 − (^) −^1 1 N^ K 2 , ρ.
- The elements, Rf f (m, n), of autocorrelation matrix, Rf f , are given by
Rf f (0, 0) = E [f (m, n)f (m, n)] = E [f (m, n)^2 ]^ = μ 2 Rf f (m, n) = E [f (i, i)f (i + m, j + n)] = ρ, m, n = 1, 2 ,... , N − 1.
- In the case of the first and the second moment of η(m, n) we have
E [η(m, n)] =
∫ (^) a −a
2 a xdx^ = 0 E [η(m, n)^2 ]^ =
∫ (^) a −a
2 a x
(^2) dx = a^2 3 =^ ση^ , while E [η(m, n)η(m′, n′)] = E [η(m, n)] E [η(m′, n′)] = 0, because of the indepen- dence between the entries of the noise in the noise matrix.
- From the previous calculation we see that Sηη (0, 0) = E [η(m, n)^2 ]^ = ση and Sηη (m, n) = 0, m, n = 1, 2 ,... , N − 1.
- From notes, we know that the Wiener filter is given by
Hr (u, v) = H
∗(u, v)Sf f (u, v) |H(u, v)|^2 Sf f (u, v) + Rηη (u, v) , and the restored image is Fˆ (u, v) = G(u, v)Hr (u, v) = G(u, v) (^) |H(u, vH)|∗ 2 (Su, v)Sf f^ (u, v) f f (u, v) +^ Rηη (u, v)^ , where G(u, v) = F {g(m, n)}.
- The inverse filtering was implemented using a threshold value to avoid dividing by zero. The results obtained for three different threshold values are shown in Fig. 4. From notes we know that the inverse filtering is extremely sensitive to noise; therefore, we will obtain a poor performance in cases of low signal-to-noise ratio (SNR) as it is shown in Fig. 4.
k=900; L=15; r0=5; r1=10; % Blurring model h(1:L,1:L)=L^(-2); % Wiener Filter H=fft2(h,M,N);
for i=1: switch RestorationType case 1 % Inverse Filtering % Threshold the min. values Thr=0.05; % Try: 0.01, 0.05 and 0. Idx=abs(H)<=Thr; H(Idx)=Thr; % Calculate the FFT of the degraded image G=fft2(g(:,:,i),M,N); % Convolve in frequency FHAT=G./H; fhat=ifft2(FHAT); case 2 % Wiener Filtering switch Quantities case 0 % Use theoretical values % Calculating rho and mu for the autocorrelation functions rho=r0^2(1-k/(MN-1))(1-k/(MN))+2r0r1(k/(MN-1))(1-k/(MN))+... r1^2k(k-1)/(MN(MN-1)); mu2=r0^2(1-k/(MN))+r1^2k/(M*N); % Autocorrelation of the image class and PSD rff(1:M,1:N)=rho; rff(1,1)=mu2; Sff=fft2(rff); % Autocorrelation of the noise and PSD
rnn=zeros(M,N); rnn(1,1)=A(i)^2/3; Snn=fft2(rnn); Hr=conj(H).Sff./((H.conj(H)).Sff+Snn); case 1 % Try to approximate the theoretical values % Experimental autocorrelation of the images: generate % a sample image and compute its autocorrelation Pr1=k/(MN); SampImg=rand(M,N); SampImg=double(r0(SampImg>Pr1))+double(r1(SampImg<=Pr1)); rff=xcorr2(SampImg)/(MN); rff=rff(M:(2M-1),N:(2N-1)); Sff=fft2(rff); % Experimental autocorrelation of the noise, similarly eta=2A(i)(rand(M,N)-0.5); rnn=xcorr2(eta)/(MN); rnn=rnn(M:(2M-1),N:(2N-1)); Snn=fft2(rnn); Hr=conj(H).Sff./((H.conj(H)).Sff+Snn); case 2 SNR=10^(12/10); % Try others (rff./rnn Theoretical value) Hr=conj(H)./((H.conj(H))+1./SNR); end FHAT=Hr.*fft2(g(:,:,i)); fhat=real(ifft2(FHAT)); end
figure(i); imshow(g(:,:,i),[]); print(i,’-deps’,strcat(’../Solutions/Prob03InvF’,num2str(i),’.eps’));
Table 1: The average and the standard deviation, in seconds, of the computing-time of the convolution between a Wiener filter and a degraded image. Calculations were performed in spatial-domain (SD) as well as in the frequency-domain (FD). Average of convolution’s computing time. Mask: 5× 5 Mask: 11× 11 Mask: 15× 15 Mask: 128× 128 SD FD SD FD SD FD SD FD fdeg001 0.0115 0.0055 0.0120 0.0085 0.0205 0.0080 1.3825 0. fdeg01 0.0070 0.0065 0.0130 0.0060 0.0195 0.0045 1.3825 0. fdeg2 0.0045 0.0080 0.0115 0.0065 0.0190 0.0075 1.3935 0. Standard deviation of convolution’s computing time. fdeg001 0.0237 0.0051 0.0052 0.0037 0.0022 0.0041 0.0055 0. fdeg01 0.0047 0.0049 0.0047 0.0050 0.0022 0.0051 0.0044 0. fdeg2 0.0051 0.0041 0.0037 0.0049 0.0031 0.0055 0.0075 0.
Additional Problems on Image Restoration
4 Problem
The type of noise corrupting Lincoln image is periodic. This is a typical example of elec- trical or electromechanical interference arising during the image acquisition step. In such scenarios, the periodic noise can be reduced significantly via frequency-domain filtering. Hence, we employ to different types of notch filters: i) a Butterworth notch filter of order n, and ii) a Gaussian notch filter. The transfer function of such filters are
HB (u, v) =
1 1+
( (^) D (^20) D 1 (u,v)D 2 (u,v)
)n , if ||Di|| 2 ≤ D 0 , i = 1, 2
1 , otherwise
HG(u, v) =
1 − exp
( (^) D 1 (u,v)D 2 (u,v) D^20
, if ||Di|| 2 ≤ D 0 , i = 1, 2 1 , otherwise
where D = (u 0 , v 0 ) is the location of the center of a band of frequencies to reject, D 0 is the radius of such band, and Di(u, v) = √(u − M/2 + (−1)iu 0 )^2 + (v − N/2 + (−1)iv 0 )^2 , i = 1 , 2, formulation that considers the symmetry of the DFT. The results of applying a Gaussian filter and two Butterworth filters (of order n = 2
and n = 10) are shown in Fig. 7. The frequency bands rejected by the filters are shown in Table 2. Clearly, the poriodic noise is reduced in all cases. Note that for a Butterworth filter of degree n = 10 the decay is much more abrupt than in the case of the same type of filter with degree n = 2 or for the Gaussian filter. Note also that ringing effects does not appear, to the naked eye at least, because the filters are not ideal. Notes: i) Recall that if we want to perform a linear convolution operation, then we have to zero pad the image and the filter properly before multiplying in the frequency domain. Note that, in this type of restoration we do not zero pad the image. You should be able to answer yourself why this does not affect the restored image. ii) The symmetry in the spatial frequency domain is important because the images in the spatial domain are real. iii) Some of you still use the metrics Q-index and MSE in this problem, so it is important that you understand that those are not useful in this problem because we don’t have the original image. Therefore, only visual inspection and roughness coefficient are valid metrics. A MATLAB function implementing the filters is the following:
function H=NotchFilters(D0,SizeImg,FilterType,d) % % NotchFilters: function to create a Gaussian or Butterworth notch filter % % Inputs: % D0: matrix of size Mx3 where M is the center of the frequencies to % eliminate. The first two columns give the position and the last one % the radius of the notch % SizeImg: [M N] the dimension of the filter % FilterType: ’Gauss ’ Gaussian notch filter; ’Butter’ Butterworth % filter % d: order of the Butterworth filter % % Output: % H: the notch filter
M=SizeImg(1); N=SizeImg(2);
Table 2: Location of the center of the rejection frequency bands, D, as well as the radius, D 0 , of bands to reject using the notch filters in Problem 1. Use symmetry to eliminate the conjugate symmetric frequencies. D(u 0 , v 0 ) u 0 v 0 D 0 -147 -155 6 -25 -38 7 -160 -91 7 -118 -130 6 -162 -130 2 -174 130 10
5 Problem
It is not hard to show that the restoration filter Hr (u, v), in the frequency domain, that minimizes (^127) ∑ m=
∑^127
n=
|(q ⊗ fˆ )(m, n)|^2
subject to the constraint
∑^127 m=
∑^127
n=
|(h ⊗ fˆ )(m, n) − g(m, n)|^2 = (128)^2 ση ,
is given by Hr (u, v) = H
∗(u, v) |H(u, v)|^2 + γ|Q(u, v)|^2 , where Hr = F {h(m, n)}, Q = F {q(m, n)}, and γ is a parameter adaptively calculated using an iterative procedure, which tries to satisfy the constraint. Fig. 8 shows the results obtained for the restoration of the figures fdeg001, fdeg01, and fdeg2. The significance of the minimization is that it reduces the roughness of the image, while the constraint guarantees that h ⊗ fˆ − g is equivalent to the noise η, in a second-moment sense. Below you can find a MATLAB script to calculate the restoring filter using constrained optimization. However, you are encouraged to take a look to the following MATLAB’s Image Processing Toolbox functions: deconvwnr and deconvreg. Also take a look at: deconvblind and deconvlucy.
%######################################### % Problem 02 HW #5 ECE 533 Spring 07 %######################################### clear all; clc; load wiener.mat %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Data g=cat(3,fdeg001,fdeg01); g=cat(3,g,fdeg2); [M N]=size(g(:,:,1)); A=[0.001 0.01 2]; k=900; L=15; r0=5; r1=10; % Blurring model h(1:L,1:L)=L^(-2); % Wiener Filter H=fft2(h,M,N); % Laplacian q=[0 0 0; 0 -1 0 -1 4 - 0 -1 0]; Q=fft2(q,M,N); % Parameter gamma=[0.5 0.7 2]; Dgamma=[10^-6 10^-6 10^-6]; MaxIters=5000; G=zeros(3,MaxIters); for i=1: % Power of the noise Sigma=A(i)^2/3; % Tolerance (in %) a=20/100Sigma; % Stop condition NoiseP=MNSigma; % Counter for number of iterations ctr=0; while (1) % Calculate the restoration filter Hr=conj(H)./((H.conj(H))+gamma(i)(Q.conj(Q)));
7 Problem
The general equation of the Wiener filter is
Hr (u, v) = H
∗(u, v)Sf f (u, v) |H(u, v)|^2 Sf f (u, v) + Rηη (u, v).
If the noise is negligible, then Rηη (u, v) ≈ 0. Hence, the previous equation transforms in
Hr (u, v) ≈ (^) H(u, v^1 ) ∴ Hr (u, v) = exp{(u^2 + v^2 )/(2σ^2 )} ,
i.e., an inverse filter.
8 Problem
Assume that E[N (n, m)] = 0 and that the noise and the signal are uncorrelated, then
E[|G(u, v)|^2 ] = E[|H(u, v)F (u, v) + N (u, v)|^2 ] = E[|H(u, v)|^2 |F (u, v)|^2 + 2Re{H(u, v)F (u, v)N (u, v)} + |N (u, v)|^2 ] = E[|H(u, v)|^2 |F (u, v)|^2 ] + 2Re{E[H(u, v)F (u, v)N (u, v)]}
- E[|N (u, v)|^2 ] (by linearity of expectation) = |H(u, v)|^2 E[|F (u, v)|^2 ] + 2Re{H(u, v)E[F (u, v)N (u, v)]}
- E[|N (u, v)|^2 ] (H(u, v) is not random) = |H(u, v)|^2 E[|F (u, v)|^2 ] + 2Re{H(u, v)E[F (u, v)]E[N (u, v)]}
- E[|N (u, v)|^2 ] (noise and signal are uncorrelated) = |H(u, v)|^2 E[|F (u, v)|^2 ] + E[|N (u, v)|^2 ] (noise is zero-mean).
Therefore, using the definitions Sf f (u, v) , E[|F (u, v)|^2 ], Sηη (u, v) , E[|N (u, v)|^2 ], and Sgg (u, v) , E[|G(u, v)|^2 ] we obtain
Sgg (u, v) = |H(u, v)|^2 Sf f (u, v) + Sηη (u, v).
9 Problem
- Since S (^) fˆ fˆ (u, v) = |R(u, v)|^2 |Sgg (u, v)|^2 , from the previous problem we obtain
S (^) fˆ fˆ (u, v) = |R(u, v)|^2 |Sgg (u, v)|^2 = |R(u, v)|^2 (|H(u, v)|^2 Sf f (u, v) + Sηη (u, v))^.
If we assume that S (^) fˆ fˆ (u, v) = Sf f (u, v) and then solve for R(u, v) we obtain
R(u, v) =
|H(u, v)|^2 + S Sηηf f^ ((u, vu, v))
- When R(u, v) is as in equation (1) then
Fˆ (u, v) =
|H(u, v)|^2 + S Sf fηη^ ((u, vu, v))
G(u, v).
