Search code examples
matlabfftoctavedft

Reading out the direction and frequency of a point in the Matlab 2D DFT fft2() function

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
Solution

  • 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);