problem with fourier spectrum

I need to show the Fourier spectrum of the image and determine one dimensional spectrum in the form: S(theta)=sum(Sr(theta)), where the limits of the sum are r=1 and r=Ro where Ro is the radius of a circle centered at the origin. Can someone help me? Thanks

[Walter Roberson]

What is Sr in this? Or does the r in Sr refer to the r=1 to r=Ro ? If so I still do not understand.

This sounds like you might be trying to do a radon transform of a circular subset of ... of the 2D fft of the image??

[Guus Hiddink]

For each frequency r, Sr(theta) is a 1-D function.

image

:

[Walter Roberson]

What is the equation for S subscript r applied at theta?

Still sounds like a radon transform to me.

[Guus Hiddink]

S(r,theta): S is the spectrum function and r and theta are the variables in polar coordinate system. For each frequency r, Sr(theta) is a 1-D function. Analyzing Sr(theta) for a fixed value of r yields the behavior along the circle centered on the origin.

[Walter Roberson]

Please post a link to Wikipedia for the defintion of the "spectrum function" that you are using. Or any other reference work.

[Guus Hiddink]

Digital Image Processing [Gonzales, Woods]

page:671

[Walter Roberson]

That does not define the Spectrum function, only refers to it. The definition must be somewhere else, perhaps in the same PDF. However that version of the PDF is not searchable, and I cannot be bothered to read 650+ pages of textbook to find the definition.

[Guus Hiddink]

I'm sorry, this is a radon function probably. I didn't know for function with that name. Is it now clearer?

[Walter Roberson]

I don't know, they might also be referring to the power spectrum function, also known as the spectral density, which is the square of the magnitude of the fourier transform; see http://www.ejbiotechnology.info/index.php/ejbiotechnology/article/view/v8n3-1/476

[Image Analyst]

Guus, what did you think of my demo below?

[Guus Hiddink]

In the book they are called magnitude as Fourier spectrum, and power spectrum is the square of the magnitude.

[Image Analyst]

OK, right. So???

[Guus Hiddink]

The previos comment was for Walter.

Answer 1

[Image Analyst]

Why not just use the trivial brute force approach?

1. Compute the spectrum
2. Scan every pixel computing r and theta
3. Sum spectrum and increment counts
4. divide sum by counts
5. plot

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;

%===============================================================================
% Read in a demo image.
folder = pwd;
baseFileName = 'slika 214.jpg';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
  % 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.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.  
% numberOfColorBands should be = 1.
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
  % It's not really gray scale like we expected - it's color.
  % Convert it to gray scale by taking only the green channel.
  grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize, 'Interpreter', 'None');

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') 
drawnow;

% Compute the 2D fft.
frequencyImage = fftshift(fft2(grayImage));

% Ask the user to select an option.
selectedOptionNumber = menu('Select an option','Magnitude', 'Real Part', 'Imaginary Part', 'Quit');
% Do stuff that depends on what option they selected.
switch selectedOptionNumber
  case 1
    % Take log magnitude so we can see it better in the display.
    amplitudeImage = abs(frequencyImage);
  case 2
    % Take log magnitude so we can see it better in the display.
    amplitudeImage = real(frequencyImage);
  case 3
    % Take log magnitude so we can see it better in the display.
    amplitudeImage = imag(frequencyImage);
  case 4
    return;
end
% Get rid of DC spike
amplitudeImage = amplitudeImage - mean2(amplitudeImage);

% Take log magnitude so we can see it better in the display.
ImageToDisplay = log(abs(amplitudeImage));

minValue = min(min(amplitudeImage))
maxValue = max(max(amplitudeImage))
subplot(2, 2, 2);
imshow(ImageToDisplay, []);
caption = sprintf('Magnitude');
title(caption, 'FontSize', fontSize);
axis on;
% zoom(10)

xCenter = columns / 2;
yCenter = rows / 2;
sumTheta = zeros(1, 361);
countsTheta = zeros(1, 361);
numRBins = ceil(sqrt(rows^2+columns^2));
sumR = zeros(1, numRBins);
countsR = zeros(1, numRBins);
for row = 1 : rows
  for col = 1 : columns
    r = round(sqrt((col - xCenter)^2 + (row - yCenter)^2)) + 1;
    theta = round(atan2d(row-yCenter, col-xCenter)) + 181;
    sumTheta(theta) = sumTheta(theta) + amplitudeImage(row, col);
    countsTheta(theta) = countsTheta(theta) + 1;
    sumR(r) = sumR(r) + amplitudeImage(row, col);
    countsR(r) = countsR(r) + 1;
  end
end
Stheta = sumTheta ./ countsTheta;
Sr = sumR ./ countsR;
subplot(2, 2, 3);
plot(Stheta, 'b-', 'LineWidth', 2);
xlabel('Angle in degrees', 'FontSize', fontSize);
ylabel('S \theta', 'FontSize', fontSize);
grid on;
xlim([0,360]);
title('S \theta', 'FontSize', fontSize);
axis on;

subplot(2, 2, 4);
plot(Sr, 'b-', 'LineWidth', 2);
xlabel('Radius', 'FontSize', fontSize);
ylabel('Sr', 'FontSize', fontSize);
grid on;
xlim([0,360]);
title('Sr', 'FontSize', fontSize);
axis on;

[Guus Hiddink]

I'm sorry for my late reply. Thank you for demo, but something isn't right. I know what graphics I should get for some images, and I try your code with these images and results are different, not the same.

[Guus Hiddink]

images and results above:

[Image Analyst]

What do you want to do with those images? I don't think it makes too much sense to deal with polar representation of spectra when you have a perpendicular pattern like the first one. And I don't know what your plans are for the second one, like do you want to measure isotropy or something? Does each image want a unique thing measured, or are you wanting to measure or process both images in the same way?

[Guus Hiddink]

I want to process images in the same way. I know what results I should get for images primer.jpg and primer1.jpg, but I need to process slika 214.jpg and I don't know results for that image. But the code need to be the same, as with other images.

[Image Analyst]

Guus, I still claim my code is correct. Let's try to visualize this. Imagine you took that spectrum on a sheet of paper and made a cut from the center to the right along the x axis. Now take that end at 360 degrees and move it along the circumference half way around the circle, from 360 to 270 to finally at 180 degrees. And then you took the center and pulled it to where the right side used to be. So now you've morphed that full plane into a half plane where the vertical direction is radius and the horizontal direction is angle. So now you can get Sr by taking the average vertical profile (averaging across columns), and you can get Stheta by taking the average horizontal profile (gotten by averaging along rows). This is essentially what I did, except that I didn't morph it. But the plots I compute are right, though your source may be displaying them differently.

But what I mean is, that's not the end of the analysis. You don't just compute those and say "OK, that's interesting but I'm done now." You actually use those plots to do something further. And what I asked (but didn't get an answer to) is what is the "further" actions that you want to do.

[Guus Hiddink]

Nothing. I need to solve that problem and learn how to do that in Matlab. This is just a task for educational purposes. Thank you very much, you are open my eyes.