Spatial Filtering Faces Using Cycles per Face

I'm interested in writing a code that can essentially recreate the images from these papers, but one that I can manipulate to my own cycles/face.

(https://medweb4.unige.ch/labnic/papers/PV_NN2003.pdf) (<http://files.face-categorization-lab.webnode.com/200000600-6772a6b5e6/Goffaux_2005_Perception.pdf>)

Currently I have a Gaussian Bandpass filter that I am able to use to recreate the Low Spatial Frequency faces, however it takes some manipulation that reduces my confidence in its effectiveness (For example if I want a 6 cycle/face image, I don't input "6" I have to play with numbers like 8 and 9 to get images that look similar to the ones in the papers above). Another issue with the code I currently have is that on the high spatial frequency images that it produces, is has a black background. I would like all of the images to be on a 50% (or close to it) grey background. Here is the code:

function filtered_image = gaussianbpf(im_name,d0,d1)
% Butterworth Bandpass Filter
% This simple  function was written for my Digital Image Processing course
% at Eastern Mediterranean University taught by
% Assoc. Prof. Dr. Hasan Demirel
% for the 2010-2011 Spring Semester
% for the complete report:
% http://www.scribd.com/doc/51981950/HW4-Frequency-Domain-Bandpass-Filtering
%
% Written By:
% Leonardo O. Iheme ([e-mail])
% 24th of March 2011
%
%   I = The input grey scale image
%   d0 = Lower cut off frequency
%   d1 = Higher cut off frequency
%
% The function makes use of the simple principle that a bandpass filter
% can be obtained by multiplying a lowpass filter with a highpass filter
% where the lowpass filter has a higher cut off frquency than the high pass filter.
%
% Usage GAUSSIANBPF(I,DO,D1)
% Example
% ima = imread('grass.jpg');
% ima = rgb2gray(ima);
% filtered_image = gaussianbpf(ima,30,120);
% Gaussian Bandpass Filter

image = imread(im_name);
%image = rgb2gray(image);

f = double(image);
[nx ny] = size(f);
f = uint8(f);
fftI = fft2(f,2*nx-1,2*ny-1);
fftI = fftshift(fftI);

% Initialize filter.
filter1 = zeros(2*nx-1,2*ny-1);
filter2 = zeros(2*nx-1,2*ny-1);
filter3 = zeros(2*nx-1,2*ny-1);

for i = 1:2*nx-1
    for j =1:2*ny-1
        dist = ((i-(nx+1))^2 + (j-(ny+1))^2)^.4;
        % Use Gaussian filter.
        filter1(i,j) = exp(-dist^2/(2*d1^2));
        filter2(i,j) = exp(-dist^2/(2*d0^2));
        filter3(i,j) = 1.0 - filter2(i,j);
        filter3(i,j) = filter1(i,j).*filter3(i,j);
    end
end
% Update image with passed frequencies
filtered_image = filter3.*fftI;

filtered_image = real(ifft2(filtered_image));

Fcf = fftshift(filtered_image);

S1 = log(1+abs(Fcf));
figure, imshow(S1,[])

%ORIGINAL CODE STARTS HERE
% filtered_image = ifftshift(filtered_image);
% filtered_image = ifft2(filtered_image,2*nx-1,2*ny-1);
% filtered_image = real(filtered_image(1:nx,1:ny));
% filtered_image = uint8(filtered_image);
%ENDS HERE, BELOW SHOULD BE KEPT FOR FILEOUTPUT

% lower_range = int2str(d0);
% upper_range = int2str(d1);
% 
% filename = strcat(im_name,'_',lower_range,'_',upper_range,'.png');
% 
% imwrite(filtered_image,filename)

Minimally I would like a code that could produce images consisting of the 1-6 cycles/face, 7-29 c/f and 30-infinity c/f, ideally I would like to manipulate the code to produce a wide range of images.

One of the major points of confusion for me at this stage is the difference between cycles per image and frequency as defined by the MatLab fft function. In my mind these are equivalent, can someone explain?

[Alejandro Estudillo]

Hi,

Actually facing the same problem. Was wondering whether you found a solution.