频域滤波

Question

我必须将 prewit 滤波器应用于频域中的图像。这是我正在遵循的过程。

1) 通过补零将图像的NxN矩阵转换为2*Nx2*N矩阵

2) 通过将图像乘以 (-1)^(x+y) 来使图像变换居中

3) 计算图像矩阵的 DFT

4) 创建尺寸为 2Nx2N 且中心位于坐标 (N,N) 的过滤器

5) 图像矩阵与滤波器矩阵相乘

6) 计算它的逆DFT并提取结果的实部。

7) 通过乘以 (-1)^(x+y) 来分散结果

8) 最后提取结果矩阵的左上NxN部分

我的代码如下：

% mask=[-1,0,1;-1,0,1;-1,0,1];

%read image
signal=imread('cman.pgm');
signal=double(signal);

% image has NxN dimensions
l=size(signal,1);

pad_signal=zeros(2*l,2*l);

pad_signal(1:l,1:l)=signal;

m=size(mask,1);

mask_f=zeros(2*l,2*l);

for i=-1:1
   mask_f(l+i,l-1)=-1;

   mask_f(l+i,l+1)=1;
end

x=1:2*l;

[x y]=meshgrid(x,x);

% Multiply each pixel f(x,y) with (-1)*(x+y)
pad_signal=pad_signal.*((-1).^(x+y));
mask_f=myDFT(mask_f);

%find the DFT of image
signal_dft=myDFT(pad_signal);

%multiply the filter with image
res=mask_f*signal_dft;

% find the inverse DFT of real values of result
res=real(myIDFT(res));

res=res.*((-1).^(x+y));

%extract the upper left NxN portion of the result
res=res(1:l,1:l);

imshow(uint8(res)); 

上面的方法来自一本图像处理的书。我感到困惑的是我应该使用 3x3 的窗口，因为 prewitt 过滤器是 3x3 还是我当前使用过滤器的方式正确？ （即通过将过滤器值放置在 2Nx2N 过滤器矩阵的中心并将所有其他索引值设置为 0）。 如果两者都不是，那么如何形成与图像的dft相乘的滤波器。

Answer 1

这是我解决问题的方法。我首先从算法中删除了第 2 步和第 7 步。然后通过在水平和垂直方向上交换索引的前半部分和后半部分来使变换居中。我这样做是为了使图像的变换居中。然后我在计算结果矩阵的逆 DFT 后取消了这个。我不确定为什么我的上述方法不起作用，但现在可以了。

1) 通过补零将图像的NxN矩阵转换为2*Nx2*N矩阵

2) 计算图像矩阵的 DFT

3) 通过交换行和列的前半部分和后半部分来使图像的变换居中。

4) 创建尺寸为 2Nx2N 且中心位于坐标 (N,N) 的过滤器

5) 图像矩阵与滤波器矩阵相乘

6) 计算它的逆DFT并提取结果的实部。

7) 通过在结果矩阵上重新应用步骤 3 来分散结果

8) 最后提取结果矩阵的左上NxN部分

以上是我在应用过滤时所遵循的步骤的修改版本。 这是我的代码（编辑/新版本）

function res=myFreqConv(signal,mask)
    signal=double(signal);

    l=size(signal,1);

    % padding the image matrix with zeros and making it's size equal to
    % 2Nx2N
    pad_signal=zeros(2*l,2*l);
    pad_signal(1:l,1:l)=signal;
    m=size(mask,1);
    mask_f=zeros(2*l,2*l);

    % Creating a mask of 2Nx2N dims where the prewitt filter values are  
    at 
    % the center of the mask i.e. the indices are like this 
    % [(N-1,N-1), (N-1,N), (N-1,N+1);(N,N-1), (N,N), (N,N+1); (N+1,N-1), 
    (N+1,N), (N+1,N+1)] 
    for i=-1:1
       mask_f(l+i,l-1)=-1;
       mask_f(l+i,l+1)=1;
    end

    % calculate DFT of mask
    mask_f=myDFT(mask_f);

    signal_dft=myDFT(pad_signal);

    % shifting the image transform to center
    indices=cell(1,2);
    indices{1}=[2*l/2+1:2*l 1:2*l/2];
    indices{2}=[2*l/2+1:2*l 1:2*l/2];
    signal_dft=signal_dft(indices{:});

    %multiply mask with image
    res=mask_f.*signal_dft;
    res=real(myIDFT(res));

    % shifting the image transform back to original
    res=res(indices{:});
    res=res(1:l,1:l);
end

Answer 2

您当前将过滤器填充为与图像相同大小的方法基本上是正确的。我们经常说用长度为 3 的滤波器来过滤长度为 M 的信号，但隐含的假设是我们将两者都填充到长度 M，或者可能是长度 M+3-1。

您的方法的一些细节使事情变得复杂：

1) 乘以 (-1)^(x+y) 只是转换 DFT，不需要。 （参见Foundations of Signal Processing 表 3.7“频率的循环偏移”对于一维情况。在这种表示法中，您让 k_0 为 N/2，因此左列中的 W_N 项只是在 -1 和 1 之间切换。）

2) 由于 Prewitt 过滤器仅支持 3x3 非零，因此您的输出只需为 N+2 x N+2 大小。这里要记住的公式是长度（信号）+长度（滤波器）- 1。

以下是我的处理方法：

clear
x = im2double(imread('cameraman.tif'));
[M, N] = size(x);

h = [-1 0 1;
    -1 0 1;
    -1 0 1];

P = M + size(h,1) - 1;
Q = N + size(h,2) - 1;

xPadded = x;
xPadded(P, Q) = 0;

hPadded = h;
hPadded(P,Q) = 0;

hShifted = circshift(hPadded, [-1 -1]);

H = fft2(hShifted);
X = fft2(xPadded);

Y = H .* X;
y = ifft2(Y);

yCropped = y(1:M, 1:N);
imshow(yCropped,[]);