Download Digital Image Processing - Assignment 2 Solutions | ECE 533 and more Assignments Digital Signal Processing in PDF only on Docsity!
Digital Image Processing ECE 533
Solutions Assignment 2
Department of Electrical and Computing Engineering, University of New Mexico.
Professor Majeed Hayat, [email protected]
February 27, 2008
Problem 1
- From the notes we know that Fs(u, v) = (^) ∆^1 x∆^1 y
∑^ ∞
−∞
−∞
F (u − m∆u, v − n∆v) , (1)
where ∆x = ∆y = 1 mm and ∆u = ∆v = 1 cycle/mm. Therefore, aliasing will occur in the u spatial frequency. In Fig. 1 we can see the resulting Fourier transform and the aliasing effect.
- Again, from the notes we know Fs(u, v) = (^) ∆^1 x∆^1 y
∑^ ∞
−∞
−∞
P (m∆u, n∆v)F (u − m∆u, v − n∆v) , (2)
where ∆x = ∆y = 1 mm and ∆u = ∆v = 1 cycle/mm. However, P (m∆u, n∆v) = 1 if m^2 + n^2 ≤ 1 and P (m∆u, n∆v) = 0 otherwise. Therefore, Fs(u, v) =
(m,n)∈{(0,0),(1,0),(0,1),(− 1 , 0 ,),(0,−1)}
F (u − m∆u, v − n∆v).
Hence, aliasing occurs once again.
Figure 1: From top to bottom: The spatial Fourier transform of the image F (u, v), the spatial spatial Fourier transform of the train of impulses used to sample F (u, v), and the sampled version of F (u, v). The images on the left correspond to the 3D visualization of the spatial Fourier transforms while the right images are 2D representations in gray-scale.
- No, since horizontal duplicates of F (u, v) overlap as in part 1.
- Let fˆ be a 45o^ rotation of f , and let Fˆ (u, v) = F { fˆ (u, v)}. Observation: a vertical projection of fˆ is equivalent to a 45o^ line projection of f. Why? We can think of the problem in the following way: i) Rotate f by 45o^ as shown in Fig. 3. ii) Now the observation about a vertical projection of fˆ being equivalent to a 45o^ line projection of f becomes clear. iii) Look for the horizontal slice in Fˆ. iv) By using the result
Figure 3: A 45-degree line projection of f and its resulting 1D Fourier transform.
- The vectors b 1 = ˆi + 12 ˆj and b 2 = ˆj, with |ˆi| = |ˆj| = 1 and < ˆi,ˆj >= 0, are the vectors generating the skew symmetry in the Fourier domain. Note that sin θ = √^25 , where θ is the angle between b 1 and b 2. Then, the vectors a 1 and a 2 in the domain of the image are given by a 1 = ˆi and a 2 = − 12 ˆi + ˆj.
- In the previous part we have (^) ∆x^1 ∆y = 11 = 1 sample/mm^2. In part 1 we have again 1 sample/mm^2. Although the sampling rate is the same, for the skew sampling we don’t have aliasing.
Problem 2
- The result of the convolution is given in the following array
(x ∗ h)(m, n) =
where the indexes m and n take the following values: m = − 2 , − 1 ,... , 5 (increasing going from left to right) and n = 0, 1 ,... , 6 (increasing going from bottom to top).
- If we take M ×N DFTs then we must have M ≥ 6+3−1 = 8 and N ≥ 3+5−1 = 7.
- The following code does what it is required in parts c) and d) clear all; clc; x5=[0 0 0 0 0 0; %<- n= 0 0 0 0 0 0; 1 1 1 1 0 0; 0 1 1 1 1 0; 0 0 1 1 1 1]; %<- n= x=x5(3:end,:); h5=[0 0 1 0 0 0; %<- n= 0 0 1 0 0 0; 0 0 1 1 0 0; 0 0 1 1 0 0; 0 0 1 1 1 0;]; %<- n= h=h5(:,3:5); Conv=conv2(x,h) % Padding xp=x; xp(7,8)=0; % zero padding (MATLAB takes care) Xp=fft2(xp); hp=h; hp(7,8)=0; Hp=fft2(hp); % Correct result Conv=ifft2(Xp.Hp) % WRONG RESULT Conv=ifft2(fft2(x5).fft2(h5))
Problem 3
Note that we think of the region A(x, y) as the region A(0, 0) shifted to the point (x, y). Let us denote the area of the region A by A¯. The PSF can be obtained using the definition
h(x, y) = (^) A^1 ¯
A(x,y)^ δ(α, β)dαdβ^ =
A^ ¯−^1 , if (0, 0) ∈ A(x, y) 0 , otherwise = (^) A^1 ¯ IA(x,y)(0, 0) = A^1 ¯ IA(0,0)(x, y), (5)
(˜x(0), 0) to (˜x(T ), 0) in T seconds. Define the parametric curve (˜x(t), 0), t ∈ [0, T ], as the trajectory followed by the camera, then
g(x, y) = O(f )(x, y) =
∫ T
0 f^ (x^ −^ x˜(t), y)dt. When ˜x(t) = vt, the PSF becomes
h(x, y) =
∫ T
0 δ(x^ −^ vt, y)dt^ =^ v
− 1 ∫^ vT 0 δ(x^ −^ u, y)du^ =^ v
− 1 (u(x) − u(x − vT ))δ(y).
Namely, the PSF is a segment of a horizontal delta sheet, scaled by 1/v, extending from x = 0 to x = vT.
Problem 5
The following code is an example of what you can do. Also, see the analytical solutions below the code.
clear all; clc; M=256; N=256; [X Y]=meshgrid(1:M,1:N); % Note that X and Y are increasing % i.e., the Y axis is inverted compared to the traditional axis for i=9: switch i case 1 % image 1 hor. line at the bottom Img=zeros(M,N); Img(M,:)=1; case 2 % image 2: ver. line at the bottom Img=zeros(M,N); Img(:,N/2)=1; case 3 % image 3: hor. lines at the bottom Img=zeros(M,N); Img((M-1):M,:)=1; case 4 % image 4: diagonal from left-top to right-bottom Img=double(X==Y); case 5 % image 5: f(x,y)=2+cos(2pi/256(2m)) Img=2+cos(2pi/256(2X)); case 6 % image 6: f(x,y)=2+cos(2pi/256(2m+3n)) Img=2+cos(2pi/256(2X+3*Y)); case 7 % image 7: f(x,y)=2+cos((5m+7n)/25)
Img=2+cos((5X+7Y)/25); case 8 % image 8: 15x15 square in the middle Img=zeros(M,N); Idxy=floor((256-15)/2):(floor((256-15)/2)+15); Img(Idxy,Idxy)=1; case 9 % image 8: circle of radius 10 in the middle Img=zeros(M,N); R=10; Img=double(((X-128).^2+(Y-128).^2)<R^2); end % Compute the FFT in dB IMG=fftshift(fft2(Img)/(MN)); MAGIMG=10log10(abs(IMG+eps).^2);
figure(1) subplot(2,2,1); %Image imagesc(Img); xlabel(’m’); ylabel(’n’); title(’f(m,n)’,’FontName’,’Times’,’FontSize’,18); subplot(2,2,2); %Image in 3D surf(X,Y,Img,’FaceColor’,’interp’,’EdgeColor’,’none’); xlabel(’m’); ylabel(’n’); zlabel(’f(m,n)’,’FontName’,’Times’,’FontSize’,18); subplot(2,2,3); %Image of the FFT imagesc(MAGIMG); xlabel(’u’); ylabel(’v’); title(’F(u,v)’,’FontName’,’Times’,’FontSize’,18); subplot(2,2,4); %Image in 3D surf(X,Y,MAGIMG,’FaceColor’,’interp’,’EdgeColor’,’none’); xlabel(’u’); ylabel(’v’); zlabel(’F(u,v)’,’FontName’,’Times’,’FontSize’,18); pause end
a. f (m, n) = δ(n); F (k, ) = N −^1 δ(k). b. f (m, n) = δ(m); F (k,) = N −^1 δ(). c. f (m, n) = δ(n) + δ(n − 1); F (k,) = N −^1 (1 + e−j^2 π/N^ )δ(k). d. f (m, n) = ∑N r=0^ −^1 δ(m − r, n − r); F (k,) = N −^1 δ(k + ), a diagonal line perpen- dicular to the original diagonal line. e. F (k,) = 2δ(k, ) + 0. 5 δ(k + 2) + 0. 5 δ(k − 2). f. F (k,) = 2δ(k, ) + 0. 5 δ(k + 2, + 3) + 0. 5 δ(k − 2 , ` − 3).
subplot(3,3,2); image(ImgDown); % colormap(gray); xlabel(’Downsampled’) title(strcat(’Downsample factor \Deltax=\Deltay=’,num2str(DS)));
% Define the downsampled grid [Xd Yd]=meshgrid(1:DS:N,1:DS:M); % Loop along the interpolation methods for j=1:size(Str,1) % Interpolate the image ImgInterp=interp2(Xd,Yd,ImgDown,X,Y,Str(j,:)); % Compute the MSE. Depending on the interpolation used some values % can be returned as NaN, so it is safer (but slower) check pixel % by pixel to compute the error. error2=0; for ii=1:M for jj=1:N if isfinite(ImgInterp(ii,jj)) error2=error2+(Img(ii,jj)-ImgInterp(ii,jj)).^2; else % if the interpolated value turns out to be inf or NaN % then we assume the mean error of the gray scale: 128 error2=error2+(128).^2; end end end Error2(k,j)=error2/(M*N);
subplot(3,3,2+j); image(ImgInterp); % colormap(gray); xlabel(strcat(upper(Str(j,:)),’, MSE=’,num2str(Error2(k,j),5))); end
% Define the interpolation methods and the downsampling factors Str=strvcat(’nearest’,’bilinear’,’bicubic’); ImgInterp=imresize(ImgDown,[M N],’nearest’);
% Compute the MSE. Depending on the interpolation used some values % can be returned as NaN, so it is safer (but slower) check pixel % by pixel to compute the error. error2=0; for ii=1:M for jj=1:N if isfinite(ImgInterp(ii,jj)) error2=error2+(Img(ii,jj)-ImgInterp(ii,jj)).^2; else % if the interpolated value turns out to be inf or NaN % then we assume the mean error of the gray scale: 128 error2=error2+(128).^2; end end end Error2(k,5)=error2/(M*N); subplot(3,3,7); image(ImgInterp); colormap(gray); xlabel(strcat(’NEAREST2’,’, MSE=’,num2str(Error2(k,5),5)));
ImgInterp=imresize(ImgDown,[M N],’bilinear’); % Compute the MSE. Depending on the interpolation used some values % can be returned as NaN, so it is safer (but slower) check pixel % by pixel to compute the error. error2=0; for ii=1:M for jj=1:N if isfinite(ImgInterp(ii,jj)) error2=error2+(Img(ii,jj)-ImgInterp(ii,jj)).^2; else % if the interpolated value turns out to be inf or NaN % then we assume the mean error of the gray scale: 128 error2=error2+(128).^2; end end end Error2(k,6)=error2/(M*N);
% Read image and calculate its FFT-2D % The images available are: IRImage.gif, IRImage2.gif, and IRImage3.gif Img2Use=2; % Select the image to use SaveGif=0; % Save the images as a .gif file switch Img2Use case 1 Img=double(imread(’IRImage.gif’))/255; case 2 Img=double(imread(’IRImage2.gif’))/255; case 3 Img=double(imread(’IRImage3.gif’))/255; end N=size(Img,2); M=size(Img,1);
NxFFT=240; NyFFT=320; % Number of points of the DFT IMG=fftshift(fft2(Img,NxFFT,NyFFT)/(MN)); MAGIMG=10log10(abs(IMG+eps).^2);
% Parameters Ax=19/255; % Amplitude of the Stripe Dx=12; % Width in pixels Ay=17/255; % Amplitude of the Stripe Dy=16; % Width in pixels
% Show image, histogram and FFT-2D PlotImages(Img,MAGIMG,kf); kf=kf+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % V-STRIPES %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Create vertical stripes and compute its FFT-2D Stripes=0.1+[Ayones(1,Dy) zeros(1,Dy)]; Stripes=repmat(Stripes,1,N/(2Dy)); Stripes=repmat(Stripes,M,1); STRIPES=fftshift(fft2(Stripes,NxFFT,NyFFT)/(MN)); MAGSTRIPES=10log10(abs(STRIPES+eps).^2);
% Show vertical stripes, its histogram and its FFT-2D PlotImages(Stripes,MAGSTRIPES,kf); kf=kf+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % V-STRIPES + IMAGE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Add vertical stripes of amplitude A to the image, and % compute the FFT-2D Noisy=Img+Stripes; NOISY=fftshift(fft2(Noisy,NxFFT,NyFFT)/(MN)); MAGNOISY=10log10(abs(NOISY+eps).^2); if SaveGif== imwrite(Noisy,’IRImageVStripes’,’gif’); end
% Show image + vertical stripes, its histogram and its FFT-2D PlotImages(Noisy,MAGNOISY,kf); kf=kf+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % FREQUENCY ENHANCEMENT: LOCATE THE POSITION OF THE FIRST DELTA AND THEN % USE THE KNOWLEDGE ABOUT THE FREQUENCY: 20 LINES (OBTAINED SIMPLY % COUNTING) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ENHANCED=NOISY; Min=min(min(abs(NOISY))); ENHANCED(121,11:20:320)=Min; % <- it can be zero also Enhanced=ifft2(ifftshift(ENHANCED)); MAGENHANCED=10*log10(abs(ENHANCED+eps).^2);
% Show resulting image PlotImages(Enhanced,MAGENHANCED,kf); kf=kf+1;
A MATLAB code used to generate Figs. 1 and 2
clear all; clc; kf=1; % Parameters Du=0.01; Dv=0.01;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Sampling with impulses %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
surf(Xd,Yd,Rs,’FaceColor’,’interp’,’EdgeColor’,’none’); xlabel(’u’); ylabel(’v’); zlabel(’F_{s}(u,v)’,’FontName’,’Times’,’FontSize’,18); subplot(3,2,6) imshow(Rs,[]); xlabel(’u’); ylabel(’v’); %Title(’F_{s}(u,v)’,’FontName’,’Times’,’FontSize’,18); print -deps Solutions/Delta.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Sampling with circ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Create a circ function R=3; IdxX=(x-R):(x+R); IdxY=(y-R):(y+R); D0(IdxX,IdxY)=double(not(CreateShapes(’circle’,R))); % D=[D0(x:end,y:end) D0(x:end,:) D0(x:end,:) D0(x:end,:) D0(x:end,1:y); ... % D0(:,y:end) D0 D0 D0 D0(:,1:y);... % D0(:,y:end) D0 D0 D0 D0(:,1:y);... % D0(:,y:end) D0 D0 D0 D0(:,1:y);... % D0(1:x,y:end) D0(1:x,:) D0(1:x,:) D0(1:x,:) D0(1:x,1:y)]; D=[Z0(x:end,y:end) D0(x:end,:) Z0(x:end,1:y); ... D0(:,y:end) D0 D0(:,1:y);... Z0(1:x,y:end) D0(1:x,:) Z0(1:x,1:y)]; Rs=conv2(Z,D); [Xd Yd]=meshgrid(linspace(-5,5,size(Rs,2)),linspace(-5,5,size(Rs,1)));
% Plot the results figure(kf); kf=kf+1; colormap(hot) subplot(3,2,1) surf(X,Y,Z,’FaceColor’,’interp’,’EdgeColor’,’none’); xlabel(’u’); ylabel(’v’); zlabel(’F(u,v)’,’FontName’,’Times’,’FontSize’,18); subplot(3,2,2) imshow(Z,[]) xlabel(’u’); ylabel(’v’); %Title(’F(u,v)’,’FontName’,’Times’,’FontSize’,18); MyFigStyle(1) subplot(3,2,3) surf(linspace(-5,5,size(D,1)),linspace(-5,5,size(D,1)),... D,’FaceColor’,’interp’,’EdgeColor’,’none’);
xlabel(’u’); ylabel(’v’); zlabel(’\delta_{c}(u,v)’,’FontName’,’Times’,’FontSize’,18); subplot(3,2,4) imshow(D,[]); xlabel(’u’); ylabel(’v’); %Title(’\delta_{c}(u,v)’,’FontName’,’Times’,’FontSize’,18); subplot(3,2,5) surf(Xd,Yd,Rs,’FaceColor’,’interp’,’EdgeColor’,’none’); xlabel(’u’); ylabel(’v’); zlabel(’F_{s}(u,v)’,’FontName’,’Times’,’FontSize’,18); subplot(3,2,6) imshow(Rs,[]); xlabel(’u’); ylabel(’v’); %Title(’F_{s}(u,v)’,’FontName’,’Times’,’FontSize’,18); print -deps Solutions/Circle.eps
