Fourrier transform on image

How can we compute numerically the features orientation and spacing from a FFT2 applied on image with oriented grid of objects?

Answer 1

[Meg Noah]

Here's an example - the grid of features is space every 20 meters, the frequency found at 0.05/m corresponds to the frequency of dots.

% Online references for FFT's
% https://www.gaussianwaves.com/2015/11/interpreting-fft-results-complex-dft-frequency-bins-and-fftshift/
% https://blogs.uoregon.edu/seis/wiki/unpacking-the-matlab-fft/

clc
close all
clear all

% generate spatial frame data and coordinates
dx_m = 1; % [m]
dy_m = 1; % [m]
nx = 101;
ny = nx;

% a grid of values
imgData = zeros(1,nx);
imgData(10:20:nx-9) = 1;
imgData = imgData.*imgData';

if (mod(nx,2) == 0)
    X1D = dx_m.*[-nx/2:1:nx/2-1];
else
    X1D = dx_m.*[-(nx-1)/2:1:(nx-1)/2];
end

if (mod(ny,2) == 0)
    Y1D = dy_m.*[-ny/2:1:ny/2-1];
else
    Y1D = dy_m.*[-(ny-1)/2:1:(ny-1)/2];
end


% visualize spatial data
figure('Color','white');
subplot(2,1,1)
imagesc(X1D,Y1D,imgData);
title({'Image Data'},'fontsize',12);
axis equal
axis tight
colorbar
set(gca,'fontweight','bold');
xlabel('X [m]'); ylabel('Y [m]');

% now to show the power spectrum
% with frequency space grid
dfy = 1/(ny*dy_m);
dfx = 1/(nx*dx_m);
if (mod(nx,2) == 0)
    FX = dfy.*[-nx/2:1:nx/2-1];
else
    FX = dfy.*[-(nx-1)/2:1:(nx-1)/2];
end

if (mod(ny,2) == 0)
    FY = dfy.*[-ny/2:1:ny/2-1];
else
    FY = dfy.*[-(ny-1)/2:1:(ny-1)/2];
end

subplot(2,1,2)
% imagesc(FX,FY,20*log10(abs(fftshift(fft2(imgData)))));
imagesc(FX,FY,(abs(fftshift(fft2(imgData)))));
axis equal; axis tight; colorbar
set(gca,'fontweight','bold');
xlabel('Frequency [1/m]'); ylabel('Frequency [1/m]');
title('FFT2D Output','fontsize',12);

Adding the ability to rotate the image

clc
close all
clear all

% generate spatial frame data and coordinates
dx_m = 1; % [m]
dy_m = 1; % [m]
nx = 101;
ny = nx;
rotation = 30;

% a grid of values
imgData = zeros(1,nx);
imgData(10:20:nx-9) = 1;
imgData = imgData.*imgData';
imgData = imrotate(imgData,rotation,'crop');

if (mod(nx,2) == 0)
    X1D = dx_m.*[-nx/2:1:nx/2-1];
else
    X1D = dx_m.*[-(nx-1)/2:1:(nx-1)/2];
end

if (mod(ny,2) == 0)
    Y1D = dy_m.*[-ny/2:1:ny/2-1];
else
    Y1D = dy_m.*[-(ny-1)/2:1:(ny-1)/2];
end


% visualize spatial data
figure('Color','white');
subplot(2,1,1)
imagesc(X1D,Y1D,imgData);
title({'Image Data'},'fontsize',12);
axis equal
axis tight
colorbar
set(gca,'fontweight','bold');
xlabel('X [m]'); ylabel('Y [m]');

% now to show the power spectrum
% with frequency space grid
dfy = 1/(ny*dy_m);
dfx = 1/(nx*dx_m);
if (mod(nx,2) == 0)
    FX = dfy.*[-nx/2:1:nx/2-1];
else
    FX = dfy.*[-(nx-1)/2:1:(nx-1)/2];
end

if (mod(ny,2) == 0)
    FY = dfy.*[-ny/2:1:ny/2-1];
else
    FY = dfy.*[-(ny-1)/2:1:(ny-1)/2];
end

subplot(2,1,2)
% imagesc(FX,FY,20*log10(abs(fftshift(fft2(imgData)))));
imagesc(FX,FY,(abs(fftshift(fft2(imgData)))));
axis equal; axis tight; colorbar
set(gca,'fontweight','bold');
xlabel('Frequency [1/m]'); ylabel('Frequency [1/m]');
title('FFT2D Output','fontsize',12);

[Chantal Hajjar]

Thank you for your answer. What I really need is to compute numerically (without looking at the graph) the grid orientation angle. Plus, my grid is not 100% rectangular. I must be able to extract both orientation angles. Is it possible to do it ?

[Meg Noah]

It should be possible. Do you have an example image?

[Chantal Hajjar]

Please find attached the image that I wish to analyze. In fact, I would like to know numerically both orientations of the grid using Fourrier Transform or using another method.

Thank you in advance

[Meg Noah]

What is the spatial dimension of the x-axis and y-axis - in other words, how many meters is this image from pixel center to pixel center along a row and along a column?