How to filter in frequency domain (multiplication after fft2 of an image)

Hi there.

I am using fft2 and fftshift on a grayscale image. The results of this fft2 contain complex values.

After that I am trying to calculate a human-like image, i.e. the fft of an image as the eye would perceive it.

I want to accomplish this by using the well known contrast sensitivity function.

This function looks like A(f) = a * (b + c*f)*exp(-d*f)^e https://www.cg.tuwien.ac.at/research/theses/matkovic/node20.html

Inserting the results from fft2 for f does not make sense to me.

How do I apply such a filter to the results of the fft?

Answer 1

[Image Analyst]

You need to use meshgrid() to create a list of all x and y coordinates of your fft2. Make sure your origin is in the center so the max distance will be sqrt(2)*(image width/2).

X = linspace(-columns/2, columns/2, columns);
Y = linspace(-rows/2, rows/2, rows);
[x, y] = meshgrid(X, Y);

Then create a distance image by using Pythagorean theorem.

distanceImage = sqrt(x.^2 + y.^2);

These distances are the "f" in the equation. Then use that to create an A image (what a bad name!):

A = a * (b + c*distanceImage) .* exp(-d*distanceImage ) .^ e;

So now you have an A image. That is a 2-D image that is the attenuation at every frequency. Make sure you used fftshift to shift the center of your spectrum to the middle of the image. Then simply multiply and inverse fft:

filteredFFT = originalFFT .* A;
filteredSpatialDomainImage = ifft2(filteredFFT);
imshow(filteredSpatialDomainImage, []);

[Stefan Mayer]

Thanks for your answer. I do not understand (yet) why I have to use meshgrid, but it does not matter for now; I will get to the bottom of it, I'm sure.

With your help I was able to apply this filter function to my fft2, but the result is odd. I applied it to the fft2 of a gaussian impulse, and transformed it back to spatial domain.

MATLAB complains that the inverse transform contains complex values and that it only displays the real part.

I was expecting some overshoot and/or undershoot at the boundaries of contrast, because that is the reason why I am using this filter. I want to visualize the eye's perception at brightness boundaries.

The output image looks almost the same as the input image, although it is slightly brighter in the centre.

I added my code, and left out the imshow.

clear; close all; clc;

wDeg = 1;    %size of image (in degrees)
nPix = 200;    %resolution of image (pixels);
[x,y] = meshgrid(linspace(-wDeg/2, wDeg/2, nPix+1));
x = x(1:end-1, 1:end-1);
y = y(1:end-1, 1:end-1);
sigma = 0.25;  %width of Gaussian (1/e half-width)
Gaussian = exp(-(x.^2+y.^2)/sigma^2);
imwrite(Gaussian, 'generated.png');

image = imread('generated.png');
[nHeight, nWidth] = size(image);

FFT = fft2(image);
FFT = fftshift(FFT);
FFT_org = FFT;

FFT = FFT/(nHeight*nWidth*255);
FFT = abs(FFT);
FFT = log(FFT + 1);

%%create csf
csfA = 2.6;
csfB = 0.0192;
csfC = 0.114;
csfD = -0.114;
csfE = 1.1;

columns = nWidth;
rows  = nHeight;

X = linspace(-columns/2, columns/2, columns);
Y = linspace(-rows/2, rows/2, rows);
[x, y] = meshgrid(X, Y);

distanceImage = sqrt(x.^2 + y.^2);
f = distanceImage;

CSF = csfA .* (csfB + csfC.*f).*(exp(csfD.*f).^csfE);

%%Apply filter
A_f = FFT_org.*CSF;
A_if = ifft2(A_f);

Can you tell me if I did something wring there, or if the error lies withing the CSF?

[Image Analyst]

I'll see if I have time tonight to try it. Of course the FFT is complex. As you know, (or see http://www.cv.nrao.edu/course/astr534/FourierTransforms.html if you don't, the FFT of a real, non-symmetric image is hermitian.