频域滤波

Question

我必须以 domain.Here 的频率对图像应用 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 的 window 因为 prewitt 过滤器是 3x3 还是我目前使用过滤器的方式正确？（即通过将过滤器值放在 2Nx2N 过滤器矩阵的中心并将所有其他索引值设置为 0）。如果两者都不是，那么filter应该如何形成才能与图像的dft相乘。

Answer 1

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

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

1) 乘以 (-1)^(x+y) 只是平移 DFT，不需要。（请参阅 Foundations of Signal Processing Table 3.7 "Circular shift in frequency" 一维情况。在该表示法中，您让 k_0 成为 N/2，因此 W_N 项在左栏中只是在 -1 和 1 之间切换。）

2) 因为 Prewitt 过滤器只有 3x3 非零支持，所以你的输出只需要大小为 N+2 x N+2。这里要记住的公式是 length(signal) + length(filter) - 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,[]);

Answer 2

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

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

2) 计算图像矩阵的DFT

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

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

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

6) 对其进行反DFT计算并提取结果的实部

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

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

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

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

频域滤波

Filtering in Frequency Domain

matlab

image-processing

dft