4
$\begingroup$

I am trying to calculate the amplitude of FFT vs wavelength for some experiments. I had put a question here previously on how to compute the FFT of an image. Thanks to everyone who provided the direction. I have prepared the program with some modifications and have computed the FFT. First here's the link to the previous question: Previous question on calculation of FFT

Here's the raw data: raw data

I have computed the following code to calculate the FFT:

data = MeanFilter[img,0]
fft1 = Fourier[ImageData[data][[All, All, 2]]];
Floor[Dimensions[fft1]/2]];
fft1 = RotateLeft[fft1, Floor[Dimensions[fft1]/2]];

I have analyzed the image and found it to be a RGBA image.I have used this code to evaluate the FFT for other two channels 1 and 3 and then added up the absolute part of FFT. Once we get the FFT pattern, which looks like this FFT image. My target is to compute the amplitude vs wavelength which will effectively give me the period of these patterns. What I wish to do is to average the FFT over circles. As from my previous post you can see that I tried to do that defining a circle of a certain radius with the no. of pixels being the radius. However, I want to know is it a good method? The result is now what I expected.

If you want to exactly know what I am aiming for here's an example. (ref:) And the portion which matches well with what I want to do.desired result

$\endgroup$

1 Answer 1

4
$\begingroup$

I'm not sure how you get your FFT image. If I apply a discrete Fourier transform to your image:

img = Import["https://i.stack.imgur.com/lquy0.jpg"];
pixels = ImageData[img][[All, All, 2]];    
pixels = pixels - Mean[Flatten@pixels];

(Setting the data to 0 mean removes the bright spot in the center of the FFT - the bin for frequency 0.)

wnd1d = N@Array[HammingWindow, #, {-.5, .5}] & /@ Dimensions[pixels];    
wnd = Outer[Times, wnd1d[[1]], wnd1d[[2]]];
Image[wnd]

(I'm using a Hamming window here - if you don't, the image borders will produce slight artifacts in the result, because the Fourier transform is periodic.)

fft = Fourier[wnd*pixels];
fft = RotateLeft[fft, Floor[Dimensions[fft]/2]];    
Image[Rescale@Abs[fft]]

enter image description here

This looks more or less as you expected, right? The bright area in the center are low frequencies - maybe an artifact of the image acquisition process? (E.g. if you scanned the image from a book, it would be slightly darker to one side of the page, that would produce low-frequency peaks like that. Or maybe it's the green text in the image? It does make the image brighter to the lower right, on average.)

Now if you use the polar transform code from my last question:

radius = Sqrt[Total[Dimensions[fft]^2]]/2.;
fftPolar = 
 ImageTransformation[
  Image[Abs[fft]], {Cos[#[[1]]], Sin[#[[1]]]} #[[2]] &, 
  Round[{\[Pi]*radius, radius}], PlotRange -> {{0, 2 \[Pi]}, {0, 1}}, 
  DataRange -> {{-1, 1}, {-1, 1}}]

ListLogPlot[Mean /@ ImageData[fftPolar, DataReversed -> True], 
 DataRange -> {0, radius}, Joined -> True, PlotRange -> All]

You get this:

enter image description here

So what does that plot mean? To be honest, I can never remember myself. So I'll generate a simple test image with a known frequency pattern:

testFrequency = 25;
pixels = Array[
   Sin[Norm[{#1, #2} - 0.5 Dimensions[pixels]]*2 π/testFrequency] &, 
   Dimensions[pixels]];
Image[pixels]

enter image description here

which gives this Fourier transform:

enter image description here

and the polar developed view looks like this:

enter image description here

where the peak of the mean plot is at:

meanData = Mean /@ ImageData[fftPolar, DataReversed -> True];
maxIndex = Position[meanData, Max[meanData]][[1, 1]]

51

Which can be understood as "there are 51 full sine waves in 2*radius", i.e. the estimated frequency would be 2*radius / maxIndex == 25.1, which is quite close to the original frequency 25 I used to generate the test pattern.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.