Filtering in Frequency Domain

Question

I have to apply prewit filter to an image in the frequency domain.Here is the procedure I am following.

1) Convert the NxN matrix of image to 2*Nx2*N matrix by padding zeros

2) Center the image transform by multiplying image with (-1)^(x+y)

3) Compute DFT of image matrix

4) Create the filter of dimensions 2Nx2N and the center at coordinates (N,N)

5) Multiply image matrix with filter matrix

6) Calculate inverse DFT of it and extract the real part of result.

7) Decentralize the result by multiplying with (-1)^(x+y)

8) Finally extract the upper left NxN part of the resultant matrix

My code is below:

% 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));

The method above is from an image processing book. What I am confused about is should I be using a window of 3x3 as prewitt filter is of 3x3 or is my current way of using the filter correct? (i.e. by placing the filter values at centre of 2Nx2N filter matrix and setting all other index values to 0) . If not either of them, then how should the filter be formed to be multiplied with the dft of image.

score 1 · Answer 1 · answered Apr 03 '15 at 17:48

Your current way of padding the filter to be the same size as the image is basically correct. We often speak loosely about filtering a length M signal with a length 3 filter, but the implicit assumption is that we are padding both to length M, or maybe length M+3-1.

Some details of your approach complicate things:

1) The multiplication by (-1)^(x+y) just translates the DFT and isn't needed. (See Foundations of Signal Processing Table 3.7 "Circular shift in frequency" for the 1D case. In that notation, you are letting k_0 be N/2, so the W_N term in the left column is just toggling between -1 and 1.)

2) Because the Prewitt filter only has a 3x3 non-zero support, your output only needs to be of size N+2 by N+2. The formula to remember here is length(signal) + length(filter) - 1.

Here's how I would approach this:

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,[]);

score 1 · Accepted Answer · answered Apr 04 '15 at 20:50

Here is how I have solved my problem. I first removed step 2 and 7 from the algorithm. Then centered the transform by swapping the first half of the indices with the second half, in both horizontal and vertical direction. I did this to center the transform of the image. Then I undid this after calculating the inverse DFT of the resultant matrix. I am not sure why my above method does not work but it does so now.

1) Convert the NxN matrix of image to 2*Nx2*N matrix by padding zeros

2) Compute DFT of image matrix

3) Centre the transform of the image by swapping the first and second half of rows and columns.

4) Create the filter of dimensions 2Nx2N and the center at coordinates (N,N)

5) Multiply image matrix with filter matrix

6) Calculate inverse DFT of it and extract the real part of result.

7) Decentralize the result by reapplying step 3 on the resultant matrix

8) Finally extract the upper left NxN part of the resultant matrix

The above is the modified version of steps that I have followed when applying my filtering. Here is my code (edited/new version)

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

2 Answers2