Find wavelength in an image

Hello everybody,

I have a picture where I can see the wavelength with my eye quite clearly, but I want matlab to find it automatically.

I tried that by doing a fft of the center column but I couldn't see a high peak in the result besides the one at x = 0.

I know that 700 pixels equals 16mm.

The code:

grayImage = rgb2gray(myImage); %myImage is the Image I've attatched 
Y = fft(grayImage(:,450));
P2 = abs(Y/675);
P1 = P2(1:675/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = (0:(675/2))/675;
plot(f,P1) 

I know this probably is a very simple problem, but I'm new to matlab.

Regards, Hendrik Sorhage

[Rik]

Did you plot the data you are putting in to fft? You need to make sure that the visual information translates well into your data.

[Hendrik]

Yes, I've attaches an image of the plot.

[Rik]

Does this pattern match your expectation? If it doesn't match what you visually see, fft will not give you a valid result. It is also possible that the wavelength is not clear enough in the line you selected to notice it in the Fourier transform.

[Hendrik]

First of all, thanks for the fast answeres. Yes, the pattern fits the image and you can see the maxima and minima have the same wavelength as in the image. I've also selected 10 other columns in the image, the results are a bit different, but none of them shows a really high peak.

[Rik]

You can remove the 0 peak by subtracting the mean:

myImage=imread('FirstShot_klein.jpg');
grayImage = rgb2gray(myImage); %myImage is the Image I've attatched 
data=grayImage(:,450);
data=double(data);data=data-mean(data);
Y = fft(data);
P2 = abs(Y/675);
P1 = P2(1:floor(675/2+1));
P1(2:end-1) = 2*P1(2:end-1);
f = (0:(675/2))/675;
plot(f,P1)

This still doesn't result in very pronounced peaks, but you could try this for every column of your image (be sure to only select the rows with relevant data and calculate a weighted average, using the number of relevant rows as weight factor).

[Hendrik]

Thanks to both of you, with your help I've managed to get get a significant peak. I've calculated the average over 20 columns in the middle. But I still have one problem, as I said 16mm equals 700 pixels, but how do I scale the x axis with that information correctly?

Answer 1

[Bjorn Gustavsson]

Ok, there are a couple of points to make.

1, Your selected column is not entirely filled with you object but also contains a bit of bright background. Since the fft is a circular fft that means that if the array is not periodic the intensity-jump from the last element to the first will give the fft a lot of high-frequency components.

2, your function contains a lot of noise at high frequencies, therefore you might benefit from smoothing your image. Something like this:

I = imfilter(I,conv2(ones(3)/9,ones(5)/25,'full'));

3, when you plot the column, you'll see that there is a lot of "almost-periodic" variations. This will give you a wide spectrum. I'd suggest that you find the row-indices for local maxima and minima. Then you can plot the difference between these indices - if your

wave has a wavelength-modulation along the column you might pick up that way.

HTH

[Hendrik]

Thanks to both of you, with your help I've managed to get get a significant peak. I've calculated the average over 20 columns in the middle. But I still have one problem, as I said 16mm equals 700 pixels, but how do I scale the x axis with that information correctly?

[Rik]

Your x-axis doesn't have either unit, because it is a Fourier transform. It will have a unit like 1/pixel.

[Hendrik]

Yes I know, the x-axis is the frequency. But is it exactly 1\pixel? Because then it would be really easy to calculate the wavelength.

[Bjorn Gustavsson]

It should be something like this:

n_pix = numel(img_column);
dpix = 1;
k = [0:n_pix/2-1,-n_pix/2:-1] / (dpix*n_pix); % if you have an even number of pixels
k = [0:n_pix/2,-n_pix/2:-1] / (dpix*n_pix);   % if you have an odd number of pixels

That is max wavenumber should be "2*pi/(2 pixels)" which is the Nyquist-limit.