0

I am getting familiarized with Matlab 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

A is a 256 x 256 matrix from which amplitude and phase can be calculated.

The question is how to extract the direction (theta) and frequency (e.g. pixels/cycle)?

EXAMPLE COMPARING MATLAB AND IMAGEJ OUTPUTS AFTER SLEUTHEYE'S ANSWER:

enter image description here

with ImageJ:

Frequency = 10.24 pixels/cycle (25 cycles)
Theta (direction) = 16.26 degrees
Real part = -1.255
Imaginary part = 10.142
Phase = arctan(10.142 / -1.255) = -82.95 degrees
Magnitude = sqrt(10.142^2 + 1.255^2) = 10.2194

and with Matlab:

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

A = fft2(double(Im) / 255);
Ashifted = fftshift(A);
Ashifted(121,153)

i = 121;
j = 153;

center = size(A) / 2 + 1;
dx = (j - center(2)) / size(A,2);
dy = (center(1) - i - 1) / size(A,1);

direction = (atan2(dy, dx))
dir_degrees = direction * (360 / (2*pi)) 
frequency = 1 /sqrt(dx*dx + dy*dy)

The output:

ans =  -1.2553 + 10.1425i
direction =  0.28379
dir_degrees =  16.260
frequency =  10.240
Antoni Parellada
  • 4,253
  • 6
  • 49
  • 114
  • @Suever I think it's a misunderstanding because I'm using the same example: it is a different question, and I got a perfect answer to it by the same person who wrote a perfect answer to my prior question, which you are linking as identical. It wouldn't make any sense for me to post again the same question. – Antoni Parellada Apr 04 '17 at 03:32
  • How is it different? It looks the same to me. You're trying to get the magnitude and the angle at a particular point – Suever Apr 04 '17 at 12:41
  • 1
    @Suever It is my understanding that there are 4 values describing a 2D DFT: the magnitude and phase (object of the previous question); and the direction (theta angle) and frequency (sometimes described as a vectorial quantity). These two last values are the object of the present question. – Antoni Parellada Apr 04 '17 at 12:59
  • Right but you literally computed the direction and frequency in your previous question yourself. – Suever Apr 04 '17 at 13:02
  • @Suever Please check the answers by SleuthEye, and you'll see they are different. – Antoni Parellada Apr 04 '17 at 13:04
  • I'm not saying they aren't different, I'm saying you didnt' need to ask another quesiton because you clearly already knew from your previous question. – Suever Apr 04 '17 at 13:06
  • @Suaver I did not, and that why I asked. I am a serious learner, and don't want to waste anybody's time. Please feel free to do whatever you feel you have to do. I don't care. – Antoni Parellada Apr 04 '17 at 13:07
  • The title says Matlab, the first part of the question says Octave, and later example says Matlab. It might be interchangeable for this example, but I'd recommend you be consistent regarding what tool you're using. There have been a number of cases of solutions being unique to one or the other. – Nick J Apr 06 '17 at 14:20

1 Answers1

1

I assume this is a follow up on this question which describe your use of ImageJ which provides direct readings of the direction and frequency parameters.

Let say you have a particular pixel A(i,j), then the direction and frequency in pixels/cycle (as similarly obtained by ImageJ) can be obtained using the following:

center = size(A)/2 + 1;
dx = (j-center(2))/size(A,2);
dy = (center(1)-i-1)/size(A,1);

direction = atan2(dy, dx); % in radians
frequency = 1/sqrt(dx*dx + dy*dy);
Community
  • 1
  • 1
SleuthEye
  • 14,379
  • 2
  • 32
  • 61
  • Fabulous! Thank you very much for all your help. To be clear, the direction is in radians, and the `i` and `j` need to be entered manually before running your code, e.g. `i = 2` and `j = 3` to retrieve the direction an frequency at point `[2,3]` of `A`, correct? – Antoni Parellada Apr 04 '17 at 03:29
  • I tried running the code on an example, and the output is different than ImageJ. I bet I'm doing something wrong, but I wonder if you wouldn't mind taking a look (I edited the OP with the output and code with both software platforms). – Antoni Parellada Apr 04 '17 at 18:08
  • Actually, it's my bad. It looks like I was off by 1 on the `dy` computation. – SleuthEye Apr 04 '17 at 21:32
  • @Toni Is that still a problem? – SleuthEye Apr 04 '17 at 21:39
  • Not at all... Just fending off texts from volunteers... It works great. Ty! Where is the `- i - 1` on your edited formula coming from? – Antoni Parellada Apr 04 '17 at 21:44
  • Let us [continue this discussion in chat](http://chat.stackoverflow.com/rooms/139909/discussion-between-sleutheye-and-toni). – SleuthEye Apr 04 '17 at 21:47