fftshift of an image

Hi all,

I'm trying out to find the wavelength of the wave seen in this image: http://www.flickr.com/photos/66468892@N07/6052890522/

I did the fftshift of the red channel [ im(:,:,1) ] and got this image: http://www.flickr.com/photos/66468892@N07/6052339583/

I did the fftshift of the grayscale [ rgb2gray ] and got this image: http://www.flickr.com/photos/66468892@N07/6052365079/

I think the first pair of lateral peaks on either side of the DC peak, in second and third images, corresponds to the grid lines (a grid pattern can be seen in the first image). How could I correlate the co-ordinates of the next pair of lower peaks to the image? Does those peaks correspond to some characteristic of the wave seen in the first image?

Thanks.

Answer 1

[Rick Rosson]

In order to use the DFT to determine the wavelength of the pattern in the image, you will need one key piece of information: some way to relate the size of the pixels in the image to an actual length in a physical unit of measure (for example: inches, centimeters, or millimeters). More specifically, you will need to know (or at least assume) a value for the number of pixels per unit of physical distance.

I will assume for now that you would like to know the wavelength in units of centimeters. Then you would need a value for the number of pixels per centimeter in the x-direction, and likewise in the y-direction (these two values could be the same or different, depending on the scaling of the image).

Let's call these parameters the pixel sampling rates for the X-direction and the Y-direction, respectively, and use the MATLAB variables Fs_x and Fs_y to represent them. The units of measure for these variables are pixels per centimeter. I will assume the following values:

   Fs_x = 50;       % pixels per centimeter
   Fs_y = 40;

We can then compute the pixel sizes in each direction as the reciprocal of the pixel sampling rates:

   dx = 1/Fs_x;     % centimeters per pixel
   dy = 1/Fs_y;

We can then compute a linear distance scale for the X- and Y-axes of the cropped image:

   [ M, N, ~ ] = size(cimg);      % pixels

   x = dx*(0:N-1)';               % centimeters
   y = dy*(0:M-1)';

Next, we compute the frequency increments for both X and Y:

   dFx = Fs_x/N;              % cycles per centimeter
   dFy = Fs_y/M;

Finally, we can create a frequency domain array for both X and Y:

   Fx = (-Fs_x/2:dFx:Fs_x/2-dFx)';     % cycles per centimeter
   Fy = (-Fs_y/2:dFy:Fs_y/2-dFy)';

Almost there (whew!). The last piece is to use fft2 and fftshift to compute the 2D DFT of the image, and then to display the resulting spectral images using Fx and Fy to represent the X- and Y- axes, respectively. From there, you should be able to find the peak values in the spectrum in both directions. The wavelength in each direction is then the reciprocal of the frequencies associated with the peak values.

One other thing: it probably makes sense to remove the DC component of the cropped image prior to computing the DFT. That way, the peaks in the spectrum will show up much more clearly (since they will not have to compete for attention with the DC value).

HTH.

Best,

Rick

>> See also this follow-up question

[Paul T John]

Thanks a lot Rick. That helped very much.

However, I would like to get two doubts cleared:

1) I did the 2D FFT and cleared the DC component. [Image is seen here: http://www.flickr.com/photos/66468892@N07/6062522739/

Now, how can I relate the peaks to the colour waves and also to the grids? Which peak corresponds to which?

2) Also, once I find out the Fx and Fy of the peak and get the wavelength in each direction, how can I find the resultant wavelength of the propagating wave?

Answer 2

[Paul T John]

I did 2D-FFT prior to fftshift. Also, the FFT-image is actually a 2D spectrum, and it appears as 1D because I used 'view(0,0)' command. I am adding the code below:

close all;
clear all;
clc

img = imread('f75p0.png');   %importing pixel information
figure(1); imshow(img);

cimg = img(1:485,1:216,1:3); %cropping the pixels
figure(2); imshow(cimg);

imgred = cimg(:,:,1);        %just the red channel
imggrey = rgb2gray(cimg);    %greyscale of the cropped image
figure(3); imshow(imgred);

dred = double(imgred);
dgrey = double (imggrey);
figure(4); mesh(dred);

imgfftred = fftshift(fft2(dred));
figure(5); surf(abs(imgfftred));
view(0,0);
shading interp;

imgfftgrey = fftshift(fft2(dgrey));
figure(6); surf(abs(imgfftgrey));
view(0,0);
shading interp;

I hope now you would be able to reproduce the code. Sorry for the confusion.

[Mehri Mehrnia]

why fftshift is necessary?

[Walter Roberson]

@Mehri Mehrnia

fft() returns back the 2-sided fft. By convention, the order of frequency bins returns back from fft() of real input signals is like

0, 1, 2, 3, 4, ... n-1, n, conj(-n), conj(-(n-1)), ... conj(-(4)), conj(-(3)), conj(-(2)), conj(-1)

but people often prefer to visualize it in the order

conj(-(n-1)), ... conj(-(4)), conj(-(3)), conj(-(2)), conj(-1), 0, 1, 2, 3, 4, ... n-1, n

The change is not mathematically significant: it is just easier for people to look at the plots.

Answer 3

[Joshua Canfield]

How do I plot the 2D DFT against Fx or Fy to obtain the usefull spacial frequency information

Could I sum columns and rows and then plot against Fx or Fy respetively?

Answer 4

[Rick Rosson]

The result of the DFT of an image should also be an image, but the graphs are showing a 1D spectrum. How did you compute the DFT of the image? Prior to calling fftshift, did you use the fft or fft2 function (or some other method)? Can you show a more complete set of MATLAB code that is reproducible by others?

[Paul T John]

Hi Rick, I'm attaching the code (as another answer). Please see below.

Answer 5

[Rick Rosson]

Hi Joseph,

Walter is correct. A single line or line segment, as in the linked image, contains an infinite number of frequencies (or wavelengths).

To expand a bit on this concept:

The DFT of an image is a transform from the space domain (in mm or cm) to the spatial frequency domain (in radians per cm or radians per mm), also known as the wave number. If we use x and y to denote the space domain, then the spatial frequency would be kx and ky respectively. A line segment or rectangle in space does transform into a sinc function in the spatial frequency domain.

It is a bit misleading or confusing to characterize it as sinc( x ) versus x. Rather, it should be sinc( kx ) versus kx. This distinction is crucial in this discussion, because there is a one-to-one relationship between kx and the x-component of the wavelength (and likewise for y). That is,

   Lx = 2*pi/kx;  
   Ly = 2*pi/ky;

So it is actually a function of wavelength, and it includes an infinite number of wavelengths, as Walter points out.

HTH.

Rick

[Gova ReDDy]

Thankyou walter and Rick...

Answer 6

[dustin Gallegos]

Hi there - Can you help me do something similar: I am trying to determine the wavelengths of the different colors in a given image. Ultimately, I would like to plot lambda(nm) vs Intensity from the picture.

Thanks !

[Walter Roberson]

I suggest you start a new Question with that, and that you post a sample image; See http://www.mathworks.com/matlabcentral/answers/7924-where-can-i-upload-images-and-files-for-use-on-matlab-answers

Also maybe you would get the information you need from http://www.mathworks.com/matlabcentral/answers/17011-color-wave-length-and-hue