Designing a 2d low pas filter in fourier domain for an image

Question

Firstly this is not a homework question. I am working on a project and my professor gave me this task so I can also learn about filters and fourier domain.

I have a picture that is 2800 x 4500 pixels. I need to make a filter that will only give an arbitrary number of "frequencies".

My thought process is this:
1. Transfer the picture in fourier domain using

F = fftshift(fft2(img)); 

2: Multiply it by a filter that has ones around the centre and zeroes elsewhere.

3: Transfer the image back into x-y domain

The problem is that I don't know how to make the filter since the picture is non-square. I tried making a circle but that didn't work. How do I make an ellipse shape?

In 1-d it is easy because the filter is a vector with ones around the centre.

I tried looking for answer to this but the pictures in question are always square pictures so I cannot apply the same thing. Also nowhere does it say how it works but rather it is just done.

Answer 1

It's been said so many times before that multiplying the Fourier transform coefficients of an image with a mask of ones and zeros in the frequency domain in order to achieve a short-cut to the ideal brickwall filtering in the time-domain is not a recommended method; and almost never preferred due to unavoidable, and unforeseen, errors it would introduce.

That being said, there are algorithms in computer geometry that describe how to approximate a cirle (or ellipse) on a discrete grid like that of a computer monitor.

The following is a circle algorithm that I remember having taken from wikipedia a decade ago. I have no idea on how it works, but just put it here for your sake.

% MIDPOINT CIRCLE ALGORITHM
Nx=256;          % image pixel size
Ny=256;
I=zeros(Ny,Nx);  % image onto which the circle will be drawn

x0 = 128;     % circle centers x-y
y0 = 128;
r = 50;       % circle radius

% ALGORITHM :
% -----------
f = 1-r;
ddfx = 1;
ddfy = -2*r;

x = 0;
y = r;

I(Ny-y0-r+1,x0) = 1;
I(Ny-y0+r+1,x0) = 1;
I(Ny-y0+1,x0+r) = 1;
I(Ny-y0+1,x0-r) = 1;

while x<y

    if f >= 0        
       y = y-1;
       ddfy = ddfy + 2;
       f = f + ddfy;
    end%if

    x = x + 1;

    ddfx = ddfx + 2;
    f = f + ddfx;

    I(Ny-y0-y+1,x0 + x) = 1;
    I(Ny-y0-y+1,x0 - x) = 1;
    I(Ny-y0+y+1,x0 + x) = 1;
    I(Ny-y0+y+1,x0 - x) = 1;
    I(Ny-y0-x+1,x0 + y) = 1;
    I(Ny-y0-x+1,x0 - y) = 1;
    I(Ny-y0+x+1,x0 + y) = 1;
    I(Ny-y0+x+1,x0 - y) = 1;

end%for

figure,imshow(I)

The output will be your mask :

I think you can modify it to work for ellipse and some fill in effects...

Answer 2

Keeping in mind @Fat32's warnings, based on this answer of mine to another question, you can construct the mask by testing whether:

$$\left(\frac{u-\operatorname{floor}(W/2)}{W}\right)^2 + \left(\frac{v-\operatorname{floor}(H/2)}{H}\right)^2 <\left( \frac{f_c}{f_s}\right)^2\tag{1},$$

where $u$ and $v$ are zero-based indexes to the fftshifted frequency domain data and $W$ and $H$ are the image width and height, $f_c$ is the cutoff frequency and $f_s$ is the sampling frequency in the horizontal and the vertical directions. $\operatorname{floor}(W/2)$ and $\operatorname{floor}(H/2)$ are the shift amounts by fftshift.