4
$\begingroup$

I have the following 8bit grey scale image:

enter image description here

The 2d FFT of this image, showing the color coded Abs[fft] values, is:

enter image description here

The code to obtain the FFT image is:

img = Import["https://i.stack.imgur.com/CA0nv.png"];

dimimg = ImageDimensions[img];
rdimimg = Reverse[dimimg];

fft = Fourier[ImageData[img]];
fftRotated = RotateLeft[fft, Floor[Dimensions[fft]/2]];

fftAbsData = Abs[fftRotated];

minc = 140;
myColorTable = 
  Flatten@{Table[{Blend[{Blue, Green, Yellow, Orange}, x]}, {x, 
      1/minc, 1, 1/minc}], 
    Table[{Blend[{Orange, Red, Darker@Red}, x]}, {x, 1/(256 - minc), 
      1, 1/(256 - minc)}]};

g = Colorize[
   ImageResize[Image[fftAbsData], {rdimimg[[1]], rdimimg[[2]]}], 
   ColorFunction -> (Blend[myColorTable, #] &)];

xfrequencies = (Range[rdimimg[[1]]] - Round[rdimimg[[1]]/2])/
   rdimimg[[1]];
yfrequencies = (Range[rdimimg[[2]]] - Round[rdimimg[[2]]/2])/
   rdimimg[[2]];

minmaxxf = MinMax[xfrequencies];
minmaxyf = MinMax[yfrequencies];

dy = minmaxyf[[2]] - minmaxyf[[1]];
dx = minmaxxf[[2]] - minmaxxf[[1]];

scaleFactor = 600;

imagefft = ImageResize[g, scaleFactor*{dx, dy}]

Questions:

Now I would like to cut out the red ring and make a backward FFT to see which objects of the original image belong to the high amplitude fft data, seen in red.

How can I cout out a circular region in the fft image?

Even without cutting out a part I am not able to reproduce the original image from fft:

inverse=Image[InverseFourier[fft]]

gives me:

enter image description here

Why does InverseFourier not reproduce the original image?

$\endgroup$
  • 1
    $\begingroup$ You may ctry to create a ellipsoidic annulus-shaped mask by something along the lines of DiskMatrix[{R1, R2}, {n1, n2}] - DiskMatrix[{r1, r2}, {n1, n2}]. Here, n1 and n2 are the image dimensions, R1, R2 are the outer radii and r1, r2` are the inner radii of the annulus. And of course, by InverseFourier but I expect you knew it already ;) $\endgroup$ – Henrik Schumacher Apr 16 at 11:06
  • $\begingroup$ @Henrik Schumacher: Thank you for the hint about how to cut out the annulus. $\endgroup$ – mrz Apr 16 at 11:44
  • $\begingroup$ You're welcome. =) $\endgroup$ – Henrik Schumacher Apr 16 at 11:51
  • $\begingroup$ @Henrik Schumacher: I have a problem with `InverseFourier´ (see above). Do you have an idea where I make an error? $\endgroup$ – mrz Apr 17 at 18:34
  • 3
    $\begingroup$ Try Chop, you can see that it works by checking that img - Image[Chop@InverseFourier[fft]] is a completely black image. (Alternately, you can call Max on the result and check that the largest value is on the order of machine precision, 10^(-16)) $\endgroup$ – C. E. Apr 17 at 20:58
3
$\begingroup$

Using your image, you can first check that InverseFourier does indeed reproduce the original image:

inverse = Image[Chop@InverseFourier[fft]];
ImageDistance[inverse, img]

returns 3.90437*10^-6

Now I would like to cut out the red ring and make a backward FFT to see which objects of the original image belong to the high amplitude fft data, seen in red.

Well, don't expect "objects" - the Fourier transform is a global operation, so you'll basically see the result of a linear (bandpass) filter by masking in the FFT domain.

But you can try easily enough. Simply create a binary mask:

mask = Array[Boole[.15 < Norm[{##}] < .25] &, 
   Dimensions[fftRotated], {{-1., 1.}, {-1., 1.}}];

HighlightImage[Image[Rescale@Log[Abs@fftRotated]], Image[mask]]

enter image description here

Multiply it with the FFT and apply inverse transform:

fftMasked = RotateRight[mask*fftRotated, Floor[Dimensions[fft]/2]];

And display the real part of the result:

Image[Rescale[Chop[Re@InverseFourier[fftMasked]]]]

enter image description here

$\endgroup$
  • $\begingroup$ thank you for the solution. I found out the the frequencies of the fft image are ALWAYS from -0.5 to 0.5, on both axes. Can you explain how I can interpret this value? $\endgroup$ – mrz Apr 23 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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