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

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.

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) 