How to average the 2D spectrum of an image from fft2 to get 1D spectrum?

Hi,

Say I have a 2D image with spatial resolution r=0.17 mm, which upon reading generates a matrix (e.g.) C=randn(140,450). I would like to calculate the the 1D radially averaged spectrum of this matrix along with kx and ky (wavenumbers in the horizontal and lateral directions).

Here is my code:

r=0.17; % spatial resolution of image

C=randn(140,450); % random matrix for illustration

[n,m]=size(C); % get dimensions

Cmean=mean(mean(C)); % get the mean value of all pixels (note that the image is in 
                     % rgb)

C=C-Cmean; % get the deviations from the mean (this is what I am interested in)

Cvar=sum(sum(C.^2))/m/n ; % Get the variance of the image values (the fft2 should 
                          % preserve the variance when integrated) 

F=(1/m/n)*fft2(C); % normalized fast Fourier transform

Fa=F.*conj(F);  % get the amplitude (Fa will have dimensions nxm)

FFT_VAR=sum(sum(Fa)); % Check if the variance computed from Cvar above and from 
                      % FFT_VAR are the same (it works)

% From here on I'm not sure how to proceed but this is what I did

spech_f=mean(Fa); % spectrum in the horizontal direction

specv_f=mean(Fa'); % spectrum in the lateral direction

FFT_spech=2*spech_f(1:m/2); % one-sided horizontal spectrum (even function)

FFT_specv=2*specv_f(1:n/2); % one-sided lateral spectrum (even function)

kx=[1:m/2]/m;  % construct wavenumbers (not sure where the resolution r should come 
               % in here)

ky=[1:n/2]/n;

My questions are then (i) how to construct the wavenumbers kx and ky (should have units of inverse resolution, i.e mm^(-1)) with the resolution r, (ii) how to average the spectra FFT_spech and FFT_specv to get one spectrum over each wavenumber? Ideally k should be k=sqrt(kx^2 + ky^2) and for each k value we need the radially averaged spectrum. Note that I understand that the dimensions of FFT_spech and FFT_specv are not the same and different wavenumbers are resolved in each, but I can cut the image to be of three equal sizez.

Help is really appreciated!

Thanks, Khaled

Answer 1

[Image Analyst]

Try this:

% 2D FFT Demo to get average radial profile of the spectrum of an image.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures.
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format longg;
format compact;
fontSize = 20;

% Change the current folder to the folder of this m-file.
if(~isdeployed)
  cd(fileparts(which(mfilename)));
end
  
% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
  % User does not have the toolbox installed.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end

% Read in a standard MATLAB gray scale demo image.
folder = fileparts(which('cameraman.tif')); % Determine where demo folder is (works with all versions).
baseFileName = 'cameraman.tif';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  % File doesn't exist -- didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
% Read in image.
grayImage = imread('cameraman.tif');
[rows, columns, numberOfColorChannels] = size(grayImage)
if numberOfColorChannels > 1
  grayImage = rgb2gray(grayImage);
end

% Display original grayscale image.
subplot(2, 3, 1);
imshow(grayImage)
axis on;
title('Original Gray Scale Image', 'FontSize', fontSize)
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);

% Perform 2D FFTs
fftOriginal = fft2(double(grayImage));
% Move center from (1,1) to (129, 129) (the middle of the matrix).
shiftedFFT = fftshift(fftOriginal);
subplot(2, 3, 2);
scaledFFTr = 255 * mat2gray(real(shiftedFFT));
imshow(log(scaledFFTr), []);
title('Log of Real Part of Spectrum', 'FontSize', fontSize)
subplot(2, 3, 3);
scaledFFTi = mat2gray(imag(shiftedFFT));
imshow(log(scaledFFTi), []);
axis on;
title('Log of Imaginary Part of Spectrum', 'FontSize', fontSize)

% Display magnitude and phase of 2D FFTs
subplot(2, 3, 4);
shiftedFFTMagnitude = abs(shiftedFFT);
imshow(log(abs(shiftedFFTMagnitude)),[]);
axis on;
colormap gray
title('Log Magnitude of Spectrum', 'FontSize', fontSize)

% Get the average radial profile
midRow = rows/2+1
midCol = columns/2+1
maxRadius = ceil(sqrt(129^2 + 129^2))
radialProfile = zeros(maxRadius, 1);
count = zeros(maxRadius, 1);
for col = 1 : columns
  for row = 1 : rows
    radius = sqrt((row - midRow) ^ 2 + (col - midCol) ^ 2);
    thisIndex = ceil(radius) + 1;
    radialProfile(thisIndex) = radialProfile(thisIndex) + shiftedFFTMagnitude(row, col);
    count(thisIndex) = count(thisIndex) + 1;
  end
end
% Get average
radialProfile = radialProfile ./ count;
subplot(2, 3, 5:6);
plot(radialProfile, 'b-', 'LineWidth', 2);
grid on;
title('Average Radial Profile of Spectrum', 'FontSize', fontSize)

[Khaled Ghannam]

Perfect. Thanks!

[Priyank Pathak]

@Image Analyst, Could you please explain what are the differences between in midrows or midcolumns and centre of matrix?

It seems same.

midRow = rows/2+1

midCol = columns/2+1

maxRadius = ceil(sqrt(129^2 + 129^2)) % Since you have matrix of 255 rows and 255 %columns

thanks

[Image Analyst]

@Priyank Pathak because of the way the ifftshift rearranges the matrix, that is where the zero frequency location is. It's at 129, not 128 or 128.5 (half way between 128 and 129)

[Yin Zhang]

in the last figure, how transfer the x unit Hz to cycle per degree

[Image Analyst]

@Yin Zhang I guess you'd have to figure it out knowing the cone angle of the image when you snapped the figure. Divide that by the number of pixels along that direction.

Answer 2

[Phuong Phan Hoai]

It's really interested. But what should I do if I want the horizontal axis show the spatial frequency (mm.^-1). Please help

Answer 3

[Phuong Phan Hoai]

hi, I really need help so please help me