2

I am getting familiarized with Octave and the function fft2(). In this toy example, I am aiming at producing the 2D DFT of the following 256 x 256 png image:

enter image description here

To be able to understand the output easily, I try to convert this image into a 256 x 256 image, eliminating color information:

Im = imread('circ.png');
pkg load image
Im = rgb2gray(Im);
figure, imshow(Im)

After this bookkeeping preliminaries I run:

A = fft2(double(Im));

enter image description here

OK. Now I take the same image, and analyze it with ImageJ, checking the output at point (157, 96), like this:

enter image description here

So the magnitude is going to be sqrt(7.448^2 + 10.458^2) = 12.83

And the phase, arctan(-10.458 / 7.448) = 54.54 degrees.

The question is, How can I get these values out of the fft2() output?

In case it makes a difference, this is how I plotted the Octave output 2D DFT:

subplot(132);
h = imshow(abs(fftshift(A)),[24 100000]);
h2 = get(h,'Parent');
set(h2,'YDir','Normal');
axis equal tight;
title("2D FFT Magnitude");

subplot(133);
h = imshow(angle(fftshift(A)),[-pi pi]);
h2 = get(h,'Parent');
set(h2,'YDir','Normal');
axis equal tight;
title("2D FFT Phase");

and this is the process in ImageJ:

enter image description here

Antoni Parellada
  • 4,253
  • 6
  • 49
  • 114
  • Using indexing? `A(157,96)`? – Suever Apr 03 '17 at 19:29
  • Sorry, I meant `A(96,157)` – Suever Apr 03 '17 at 19:33
  • @Suever Still far... Probably I don't know how to read the output correctly... I don't know... `>> A(96,157) ans = -472.89 - 2050.46i` – Antoni Parellada Apr 03 '17 at 19:34
  • Do you have ImageJ computing the FFT for you? – Suever Apr 03 '17 at 19:35
  • @Suever The output I get is in the image provided. After that I have to do the calculations manually, as included right after the image. – Antoni Parellada Apr 03 '17 at 19:40
  • But do you save the output from MATLAB and then load it in ImageJ or do you just load your original image in ImageJ and then compute the FFT in ImageJ? – Suever Apr 03 '17 at 19:42
  • @Suever The latter. I load the image on ImageJ independently. – Antoni Parellada Apr 03 '17 at 19:43
  • Well first of all it looks like they are displaying the log transform of the magnitude image – Suever Apr 03 '17 at 19:43
  • As I said, ImageJ is showing you a transformed magnitude/phase image and not the raw data that you have in `A` in MATLAB – Suever Apr 03 '17 at 19:52
  • 1
    `Ashift = fftshift(A); rad2deg(angle(Ashift(96, 157)));` gives almost the same angle output (-53.0599°). As @Suever suggested, the magnitude is probably transformed. – m7913d Apr 03 '17 at 19:54
  • @Toni As I said, it looks like a `log` transform. Also, if you want the magnitude you don't want to transform it... – Suever Apr 03 '17 at 20:01
  • I guess the transformation is a scaled log, i.e. `a * log(...)`. You can try to compare the logarithm of the magnitude at different points and see if there is a constant scale factor. – m7913d Apr 03 '17 at 20:07
  • @Suever The `log(sqrt(472.89^2 + 2050.46^2)) = 7.65`, as opposed to `12.83`. – Antoni Parellada Apr 03 '17 at 20:07
  • @Toni The values you read in MATLAB are correct. You don't know what transformation ImageJ applies. Also, take into account that MATLAB uses 1-based indexing, and I presume ImageJ uses 0-based indexing. Make sure you are looking at the same pixel! – Cris Luengo Apr 04 '17 at 01:15

1 Answers1

2

Here are a few observations which should clarify the scaling used:

  • ImageJ's X and Y positions are 0-based whereas Matlab's indexing is 1-based.
  • Incrementing ImageJ's X position corresponds to incrementing the column index in Matlab, and incrementing ImageJ's Y position corresponds to incrementing the row index in Matlab, so that an (X,Y) coordinate pair in ImageJ will be found at index (Y+1,X+1) in Matlab.
  • ImageJ displays the image with the 0 frequency component in the middle of the image, so measurements are effectively made on the equivalent of fftshift(Im)
  • ImageJ scales the 0-255 grayscale values to floating point values in the 0.0-1.0 range (i.e. scaling all values by dividing by 255)

So with that in mind, we have:

>> Ashifted = fftshift(A);
>> Ashifted(97,158)/255
ans =   7.4484 - 10.4582i
>> Ashifted(93,165)/255
ans =  12.1928 -  4.9850i

which correspond exactly to your illustrated measurements of the real and imaginary parts at positions (X,Y) = (157,96) and (X,Y) = (164,92) respectively.

Note that by the linearity property of the FFT, you could also divide the input and get the same results:

A = fft2(double(Im)/255.0);
>> Ashifted = fftshift(A);
>> Ashifted(97,158)
ans =   7.4484 - 10.4582i
>> Ashifted(93,165)
ans =  12.1928 -  4.9850i
SleuthEye
  • 14,379
  • 2
  • 32
  • 61
  • Fabulous! I have a follow-up question regarding the extraction of the frequency and direction (theta) at any point in Matlab. Unless it's a very quick answer, I'll make another question out of it. – Antoni Parellada Apr 04 '17 at 01:49