How to detect wavelength patterns using fft2?

Hello,

I would like to detect these wavelengths (arrow) :

Find my image here : https://amubox.univ-amu.fr/s/2zb7Q3eRdDCyGRD

My idea was to first calculate the fft2 of my image and then use a high-pass filter to cut off the high wavelength (background). Then find the maximum magnitude of my spectrum and get my wavelength. But it didn't work.

I know that the wavelength should be approximately equal to 2-3mm on my sensor and 10-13mm in reality (due to an optical setup that focuses my image on my camera's sensor).

clear all;
clc;
close all;

% 1 load image
image_grayscale = im2gray(im2double(imread("path"))); % Charger votre image ici

Error using imread>get_full_filename (line 579)
Unable to find file "path".

Error in imread (line 372)
fullname = get_full_filename(filename);

% 2. Remove salt & Paper
med_grayImage = medfilt2(image_grayscale, [20 20]);

% 3. clean and smooth image using gaussian filter
sigma = 1.5; % STD gaussian filter
filteredImage = imgaussfilt(med_grayImage, sigma);

% FFT2
fft_image = fftshift(fft2(filteredImage));

figure(1)
subplot(1,3,1);
imshow(image_grayscale);
subplot(1,3,2)
imshow(filteredImage);
subplot(1,3,3)
imshow(log(1 + abs(fft_image)), []);

% size l'image
[m, n] = size(image_grayscale);

% scaling pixel/mm
echelle_pixel_par_mm = 5504/24; 

%%
H=6.875;
lambda_c = 1.2*H* 70/333;
frequence_critique = 2*pi/lambda_c;

% Set the cut-off frequency of the high-pass filter
frequence_coupure_x = 24/5504 *frequence_critique;
frequence_coupure_y = 24/5504 *frequence_critique ; 

% Créer le masque de filtre passe-haut
masque_filtre = ones(m, n);
centre_x = round(n/2);
centre_y = round(m/2);
for i = 1:m
    for j = 1:n
        distance = sqrt((i - centre_y)^2 + (j - centre_x)^2);
        if distance <= frequence_coupure_x * n && distance <= frequence_coupure_y * m
            masque_filtre(i, j) = 0;
        end
    end
end

% Apply the high-pass filter by multiplying the spectrum by the mask
fft_image_filtre = fft_image .* masque_filtre;

% recontructed image
new_image = real(ifft2(ifftshift(fft_image_filtre)));

figure(2)
subplot(2,2,1)
imshow(masque_filtre)
subplot(2,2,2)
imshow(fft_image_filtre)
subplot(2,2,3)
imshow(new_image)
imcontrast

% Calculer la magnitude du spectre filtré
magnitude_spectre_filtre = abs(fft_image_filtre);

% Prendre le carré de la magnitude pour obtenir la PSD
psd = magnitude_spectre_filtre.^2;

% Afficher la PSD en utilisant une échelle logarithmique
figure;
imagesc(log(1 + psd));
colormap('jet');
colorbar;
title('Densité Spectrale de Puissance (PSD)');
xlabel('Fréquence spatiale');
ylabel('Fréquence spatiale');

% Trouver les coordonnées du maximum de la magnitude du spectre filtré
[max_mag, index] = max(magnitude_spectre_filtre(:));
[indice_ligne, indice_colonne] = ind2sub(size(magnitude_spectre_filtre), index);

% Convertir les coordonnées en fréquences spatiales
frequence_x_pixel = (indice_colonne - 1 - n/2) / n;
frequence_y_pixel = (indice_ligne - 1 - m/2) / m;

% Afficher la fréquence spatiale dominante en millimètres
disp(['La fréquence spatiale pour laquelle vous avez le maximum d''énergie est : (', num2str(frequence_x_pixel), ' mm^-1, ', num2str(frequence_y_pixel), ' mm^-1)']);

Do you know how to correctly get my wavelenght ?

Do you have any advices ?

Thank you

Answer 1

[Image Analyst]

Your hypothesis is wrong. The "wavelengths" indicated by your arrows are not high frequencies - they are low frequencies.

Anyway I would not use a Fourier transform to determine those distances when it's so straightforward to do it in the spatial domain. One simple way is to use use improfile to get an intensity profile going across the light and dark bands. Then use thresholding to determine the locations where it transitions from light to dark and dark to light.

There are other methods depending on what you want to do, like measure the band thickness all along it, or find the average width of the light or dark bands. I'm attaching an old demo of mine.

[Alexandre]

Hello,

Thank you for your answer.

Unfortunately, I can't use improfile because I've obtained thousands of images like this one and my "waves" turn in all directions. In this case, they're almost horizontal, but maybe 30 seconds later, they move from top to left.

For example :

That's why I try with fft2 because I don't know any better way to do it.

Best regards

[Image Analyst]

fft won't do it. The frequency is far too low. It would be in the first few elements of the spectrum and you won't be able to get enough precision. It would essentially fit it to a set of sine waves. The frequency of the most dominant sine wave would not really be a good estimate of the thickness of those bands.

I think you're best off doing it in the spatial domain. Try finding the number of bright bands and dark bands inside the circle. Then find the max feret distance with bwferet. Then find the areas with regionprops. Divide the area by the max feret diameter to get the average width of each bright band. You might also use bwareafilt and imfill to clean up the bright blobs.