问题 1 通过计算一维傅里叶变换实现图像二维快速傅里叶变换

实现一个函数F=dft2D(f),其中f是一个灰度源图像，F是其对应的二维快速傅里叶变换(FFT)图像.具体实现要求按照课上的介绍通过两轮一维傅里叶变换实现。也就是首先计算源图像每一行的一维傅里叶变换，然后对于得到的结果计算其每一列的一维傅里叶变换。

如果实现采用MATLAB,可以直接调用函数fft计算一维傅里叶变换。如果采用其他语言，请选择并直接调用相应的一维傅里叶变换函数。

（1）思路

首先计算源图像每一行的一维傅里叶变换，结果保存在矩阵temp中，然后对于得到的结果temp计算其每一列的一维傅里叶变换，通过两个for循环实现。

（2）代码

function [F] = dft2D(f)

%F=fft(fft(f')') 

width=size(f,1);

height=size(f,2);

for i=1:height

    temp(i,1:width)=fft(f(i,1:width));

end

for i=1:width

    temp(1:height,i)=fft(temp(1:height,i));

end

F=temp;

end

调用的主函数

（3）实现效果

（a）自定义dft2D函数实现效果（b）调用matlab自带fft2实现效果

图1 问题1实现效果

（c）自定义dft2D函数得到的矩阵（d）调用matlab自带fft2得到的矩阵

图1 问题1实现效果（续）

问题 2 图像二维快速傅里叶逆变换（10分）

实现一个函数 f=idft2D(F),其中F是一个灰度图像的傅里叶变换，f是其对应的二维快速傅里叶逆变换(IFFT)图像，也就是灰度源图像.具体实现要求按照课上的介绍通过类似正向变换的方式实现。

（1）思路（摘自《数字图像处理第3版 P189页》）

（2）代码

function [f] = idft2D(F)

%F是一个灰度图像的傅里叶变换，f是其对应的二维快速傅里叶逆变换(IFFT)图像，也就是灰度源图像

M=size(F,1);

N=size(F,2);

%取F的复共轭

F= conj(F); %调用matlab自带的求复共轭函数conj

% 将F的共轭作为离散二维傅里叶正向变换的输入，调用dft2D(F)，再次取复共轭，并除以M*N

temp=conj(dft2D(F))/(M*N);

f=temp;

end

问题 3 测试图像二维快速傅里叶变换与逆变换

对于给定的输入图像 rose512.tif, 首先将其灰度范围通过归一化调整到[0,1]. 将此归一化的图像记为f.首先调用问题1下实现的函数dft2D计算其傅里叶变换，记为F。然后调用问题2下的函数idft2D计算F的傅里叶逆变换，记为g.计算并显示误差图像d=f-g.

（1）代码

img=imread("rose512.tif");

imshow(img)

M=size(img,1);

N=size(img,2);

%归一化灰度范围

f=im2double(img)

f=f/M*N;

F=dft2D(f);

g=idft2D(F);

d=f-g

result=uint8(d)

imshow(g)

（2）效果

（a）原图（b）经过傅里叶变换与反变换后的图像

（c）二者做差得到的d矩阵与d取整数的result矩阵

图2 问题3实现效果

问题 4 计算图像的中心化二维快速傅里叶变换与谱图像

我们的目标是复现下图中的结果。首先合成矩形物体图像，建议图像尺寸为 512×512，矩形位于图像中心，建议尺寸为60像素长，10像素宽，灰度假设已归一化设为1.对于输入图像f计算其中心化二维傅里叶变换F。然后计算对应的谱图像S=log(1+abs(F)).显示该谱图像。

图 3. 矩形物体图像的傅里叶变换 (a) 源图像. (b)谱图像，注意四个角的明亮区. (c) 中心化的谱图. (d) 对数变换后的谱图。

（1）思路与代码说明

首先生成一张512*512的图片，直接把长和宽中间部分的值赋为1。再生成用于x和y方向的核，即傅里叶变换公式中的部分。

接下来是两个方向的循环，生成临时矩阵Gm和Gn，再生成傅里叶谱中心化所需要的（-1）^（x+y），对应的矩阵是C，在生成谱图像时调用了matlab自带的abs函数，它可计算数组中每个元素的幅度。

（2）代码

问题4

% 合成一幅二值图像

img=zeros(512,512);

img(227:287,252:262)=1;

imshow(img)

img = im2double(img);

subplot(1,4,1);

imshow(img,[]);

title('原始图像')


M = size(img,1);  % 长

N = size(img,2);  % 宽

Gm = zeros(M)+exp(-1i*2*pi/M);

Gn = zeros(N)+exp(-1i*2*pi/N);  % G是计算时要用的矩阵


%二维离散傅里叶变换所需要的Gm和Gn

for row = 0:M-1

    for col = 0:M-1

        Gm(row+1,col+1) = Gm(row+1,col+1)^(row * col);

    end

end

% 计算Gn

for row = 0:N-1

    for col = 0:N-1

        Gn(row+1,col+1) = Gn(row+1,col+1)^(row * col);

    end

end


% 傅里叶谱中心化

C = zeros(M,N)+(-1);%生成全为-1的矩阵

for row = 1:M

    for col = 1:N

        C(row,col) = C(row,col)^(row+col);%（-1）^(x+y)

    end

end

% 谱图像

F = abs(Gm*img*Gn)% abs函数计算数组中每个元素的幅度

subplot(1,4,2);

imshow(F,[]);

title('谱图像，注意四个角的明亮区');

F = real(Gm*(img.*C)*Gn);  % F = Gm* f *Gn，只取实部

subplot(1,4,3);

imshow(abs(F),[]);

title('中心化的谱图');

subplot(1,4,4);

imshow(log(abs(F) + 1),[]);

title('对数变换后的谱图');

（3）实现效果

图4 问题4实现效果

选做题测试更多图像的二维快速傅里叶变换

计算其他5幅图像的二维快速傅里叶变换：house.tif,house02.tif, lena_gray_512.tif, lunar_surface.tif, characters_test_pattern.tif。注意，有些图像的尺寸不是 2 的整数次幂，需要进行相应的像素填补处理。如果图像有多个通道可以选择其中的一个通道进行计算。

思路

编程思路主要分为3部分，首先读取图像，判断是否有多个通道，默认取第一个。其次取图像尺寸中较长的边，与2的整数次幂进行比较，保存在edge变量中。设置一个布尔型变量flag，用于判断图像的尺寸是否为2的整数次幂，当为true时将进行最接近图像尺寸的2的整数次幂进行填充。最后调用第4题写好的函数进行傅里叶变换（带有中心化+对数变换）。

代码

主函数代码

%问题5

img=imread("characters_test_pattern.tif");

image_size=size(img)

dimension=numel(image_size);

if dimension==3

    img=img(:,:,1);% 取一个通道

end

M=size(img,1);

N=size(img,2);

subplot(2,2,1);

imshow(img,[]);

title('原图');

%取图像尺寸中较长的边，与2的整数次幂进行比较

if M>N

    edge=M;

else

    edge=N;

end

for i=4:20 %这里在将2的4次方到20次方与edge做一个比较，确定最后要填充的数

    if 2^i>edge

        k=i;

        flag=true;

        break;

    end

    if i==20 %循环完了都没找到，说明不用填充

        flag=false;

    end

end

if flag %当图像不是2的整数次幂时

    temp=zeros(2^k,2^k);

    temp(1:M,1:N)=img;

    temp(M+1:2^k,N+1:2^k)=0;

    subplot(2,2,2);

    imshow(img,[]);

    title('填充后');

    imshow(temp);

    img=temp;

end

temp=uint8(temp);

imshow(temp)

[img,revertclass]=tofloat(img);

[F, LogF]=dft_center(img);

subplot(2,2,3);

imshow(F,[]);

title('二维傅里叶中心化');

subplot(2,2,4);

imshow(LogF,[]);

title('对数变换后');

封装好的dft_center函数代码

function [outputArg1,outputArg2] = dft_center(img)

%DFT_CENTER 此处显示有关此函数的摘要

%   此处显示详细说明

M = size(img,1);  % 长

N = size(img,2);  % 宽

Gm = zeros(M)+exp(-1i*2*pi/M);

Gn = zeros(N)+exp(-1i*2*pi/N);  % G是计算时要用的矩阵


%二维离散傅里叶变换

for row = 0:M-1

    for col = 0:M-1

        Gm(row+1,col+1) = Gm(row+1,col+1)^(row * col);

    end

end

% 计算Gn

for row = 0:N-1

    for col = 0:N-1

        Gn(row+1,col+1) = Gn(row+1,col+1)^(row * col);

    end

end

% 傅里叶谱中心化

C = zeros(M,N)+(-1);%生成全为-1的矩阵

for row = 1:M

    for col = 1:N

        C(row,col) = C(row,col)^(row+col);%（-1）^(x+y)

    end

end

F = real(Gm*(img.*C)*Gn);  % F = Gm* f *Gn，只要实部

outputArg1 =  abs(F);

outputArg2 =  log(abs(F) + 1);


end

运行效果

（1）house.tif

测试效果