PROJECT 04-01 [Multiple Uses] Two-Dimensional Fast Fourier Transform

The purpose of this project is to develop a 2-D FFT program “package” that will be used in several other projects that follow. Your implementation must have the capabilities to:

1. Multiply the input image by $(-1)^{x+y}$ to center the transform for filtering.
2. Multiply the resulting (complex) array by a real function (in the sense that the the real coefficients multiply both the real and imaginary parts of the transforms). Recall that multiplication of two images is done on pairs of corresponding elements.
3. Compute the inverse Fourier transform.
4. Multiply the result by $(-1)^{x+y}$ and take the real part.
5. Compute the spectrum.

Basically, this project implements Fig. 4.5. If you are using MATLAB, then your Fourier transform program will not be limited to images whose size are integer powers of 2. If you are implementing the program yourself, then the FFT routine you are using may be limited to integer powers of 2. In this case, you may need to zoom or shrink an image to the proper size by using the program you developed in Project 02-04.

An approximation: To simplify this and the following projects (with the exception of Project 04-05), you may ignore image padding (Section 4.6.3). Although your results will not be strictly correct, significant simplifications will be gained not only in image sizes, but also in the need for cropping the final result. The principles will not be affected by this approximation.

原理

1. 因为$\Im[f(x,y)(-1)^{x+y}]=F(u-M/2,v-N/2)$，因此频率域的中心平移可以由原图乘上$(-1)^{x+y}$实现。
2. 直接在得到的复数数组上乘上实函数即可。
3. 可以调用 matlab 的 ifft2
4. 同（1），获取实部可以用real函数
5. 复数的模可以用abs函数

代码

f=imread('Fig0418(a).tif');%输入图片

[M,N]=size(f);
P=2*M;Q=2*N;

fp=zeros(P,Q);
fp(1:M,1:N)=f(1:M,1:N);%图像填充

[Y,X]=meshgrid(1:Q,1:P);
ones=(-1).^(X+Y);

F=fft2(ones.*fp);%用（-1）^(x+y)乘以输入图像,来实现中心化变换
%(-1)^(x+y)*f(x,y)的傅里叶变换等于fftshift(fft2(f,P,Q));
% 因此上面可用一句F=fftshift(fft2(f,P,Q)); 代替

G=H.*F;%频率域相乘，这里H是实函数

g=real(ifft2(G));%反变换并取结果的实部
g=ones.*g;
g=g(1:M,1:N);%在g（x,y）提取M*N左上象限

I=abs(F);% 取频谱，即模长，直接调用abs即可


PROJECT 04-02 Fourier Spectrum and Average Value

2. Display the spectrum.
3. Use your result in (a) to compute the average value of the image.

代码

close all;clear all;clc;
I = double(I);
[m,n] = size(I);
s = sum(sum(I));
avg = s / (m * n);
subplot(1,2,1),imshow(uint8(I)),title(avg);
sp = abs(fftshift(fft2(I,2*m,2*n)));
subplot(1,2,2),imshow(uint8(log(1 + sp)),[]);


PROJECT 04-03 Lowpass Filtering

1. Implement the Gaussian lowpass filter in Eq. (4.3-7). You must be able to specify the size, M x N, of the resulting 2D function. In addition, you must be able to specify where the 2D location of the center of the Gaussian function.
2. Download Fig. 4.11(a) [this image is the same as Fig. 4.18(a)] and lowpass filter it to obtain Fig. 4.18(c).

代码

close all;clear all;clc;
f = mat2gray(f);
[M,N] = size(f);
P = 2 * M; Q = 2*N;

%生成实的对称的滤波函数H，中心在（M,N）
alf=100;
for i=1:P
for j=1:Q
H(i,j) =exp(-((i-P/2)^2+(j-Q/2)^2)/(2*alf^2));
end
end

fp=zeros(P,Q);
fp(1:M,1:N)=f(1:M,1:N);%图像填充

[Y,X]=meshgrid(1:Q,1:P);
ones=(-1).^(X+Y);

F=fft2(ones.*fp);%用（-1）^(x+y)乘以输入图像,来实现中心化变换
%(-1)^(x+y)*f(x,y)的傅里叶变换等于fftshift(fft2(f,P,Q));
% 因此上面可用一句F=fftshift(fft2(f,P,Q)); 代替

G=H.*F;%频率域相乘，这里H是实函数

g=real(ifft2(G));%反变换并取结果的实部
g=ones.*g;
g=g(1:M,1:N);%在g（x,y）提取M*N左上象限

subplot(1,3,1),imshow(f),title('原图');
subplot(1,3,2),imshow(H),title('滤波函数');
subplot(1,3,3),imshow(g),title('高斯低通滤波后的图像');


PROJECT 04-04 Highpass Filtering Using a Lowpass Image

1. Subtract your image in Project 04-03(b) from the original to obtain a sharpened image, as in Eq. (4.4-14). You will note that the resulting image does not resemble the Gaussian highpass results in Fig. 4.26. Explain why this is so.
2. Adjust the variance of your Gaussian lowpass filter until the result obtained by image subtraction looks similar to Fig. 4.26(c). Explain your result.

原理

$f_{hp}(x,y)=f(x,y)-f_{lp}(x,y)$

代码

close all;clear all;clc;
f = mat2gray(f);
[M,N] = size(f);
P = 2 * M; Q = 2*N;

%生成实的对称的滤波函数，中心在（M,N）
alf=100;%第二问把这里改成200
for i=1:P
for j=1:Q
H(i,j) =exp(-((i-P/2)^2+(j-Q/2)^2)/(2*alf^2));
end
end

fp=zeros(P,Q);
fp(1:M,1:N)=f(1:M,1:N);%图像填充

[Y,X]=meshgrid(1:Q,1:P);
ones=(-1).^(X+Y);

F=fft2(ones.*fp);%用（-1）^(x+y)乘以输入图像,来实现中心化变换
%(-1)^(x+y)*f(x,y)的傅里叶变换等于fftshift(fft2(f,P,Q));
% 因此上面可用一句F=fftshift(fft2(f,P,Q)); 代替

G=H.*F;%频率域相乘，这里H是实函数

g=real(ifft2(G));%反变换并取结果的实部
g=ones.*g;
g=g(1:M,1:N);%在g（x,y）提取M*N左上象限

subplot(2,2,1),imshow(f),title('原图');
subplot(2,2,2),imshow(H),title('滤波函数');
subplot(2,2,3),imshow(g),title('高斯低通滤波后的图像');



PROJECT 04-05 Correlation in the Frequency Domain

Download Figs. 4.41(a) and (b) and duplicate Example 4.11 to obtain Fig. 4.41(e). Give the (x,y) coordinates of the location of the maximum value in the 2D correlation function. There is no need to plot the profile in Fig. 4.41(f).

代码

close all;clc;clear all;
[M1, N1] = size(f1);
[M2, N2] = size(f2);
P = 298;
Q = 298;
[Y,X]=meshgrid(1:Q,1:P);
ones=(-1).^(X+Y);
fp1 = zeros(P, Q);
fp2 = zeros(P, Q);
fp1(1:M1, 1:N1) = f1(1:M1, 1:N1);
fp2(1:M2, 1:N2) = f2(1:M2, 1:N2);
Fp1 = fft2(ones.*fp1);
Fp2 = fft2(ones.*fp2);
Fp = Fp2 .* conj(Fp1);% 求共轭
fp = ifft2(Fp);
fp = real( ones.* fp);
fp = mat2gray(fp);
max_value = max(max(fp));
[row,col] = find(fp == max_value);
subplot(1,3,1),imshow(f1),title(['max value is : ', num2str(max_value)]);
subplot(1,3,2),imshow(f2),title(['(', num2str(row), ',', num2str(col),')']);
subplot(1,3,3),imshow(fp);


以 figure4.4 为例检验谱平面随着图像旋转而旋转的性质

代码

close all;clear all;clc;
I = double(I);
[m,n] = size(I);
s = sum(sum(I));
avg = s / (m * n);
subplot(2,2,1),imshow(uint8(I)),title('原图');
sp = abs(fftshift(fft2(I,2*m,2*n)));
subplot(2,2,2),imshow(uint8(log(1 + sp)),[]),title('原图频谱');
I=imrotate(I,45);
sp = abs(fftshift(fft2(I,2*m,2*n)));
subplot(2,2,3),imshow(uint8(log(1 + sp)),[]),title('原图旋转45°的频谱');
I=imrotate(I,45);
sp = abs(fftshift(fft2(I,2*m,2*n)));
subplot(2,2,4),imshow(uint8(log(1 + sp)),[]),title('原图旋转90°的频谱');