1

Here is a good starting to fiter image in the frequency domain

lena = ColorSeparate[ImageResize[ExampleData[{"TestImage", "Lena"}],{121,121}]][[1]];

(Create a padding image to avoid the wraparound problem)

lenapad = ImagePad[lena, {{0, 122}, {122, 0}}]

padimage

data = ImageData[lenapad];

centereddata = Table[(-1)^(x + y)*data[[x, y]], {x, 1, 243}, {y, 1, 243}];

Image[centereddata]

enter image description here

(compute the DFT of image)

dft = Fourier[centereddata];

(Gaussian filter to blur image: spatial domain)

filter = ImageData@
   ImageAdjust[ImagePad[Image@GaussianMatrix[3], (243 - 7 + 1)/2]];

(DFT of the filter frequency domain)

huv = Im[Fourier[
   Table[filter[[x, y]]*(-1)^(x + y), {x, 1, 243}, {y, 1, 243}]]];
huvcenter=Table[Complex[0, huv[[x, y]]]*(-1)^(x + y), {x, 1, 243}, {y, 1, 243}];

(Apply the filter to the image)

hfp = huvcenter dft;

(Recuperate the filtering version of lena)

refinversehfp = Re[InverseFourier[hfp]];
refinverse = Table[refinversehfp[[x, y]]*(-1)^(x + y), {x, 1, 243}, {y, 1, 243}];

Image[refinverse] // ImageAdjust

enter image description here

The result is not good compared to the spatial filter

ImageConvolve[lena, GaussianMatrix[3]]

enter image description here

BetterEnglish
  • 2,026
  • 13
  • 19
  • 1
    Possible duplicate? http://mathematica.stackexchange.com/questions/110914/how-to-use-2d-fourier-analysis-to-clean-the-noise-in-an-image?lq=1 – dr.blochwave May 26 '16 at 06:39

1 Answers1

3

I believe that you are not converting your filter properly from the spatial domain to the Fourier domain. The process has three steps:

  1. Pad the spatial filter to the size of the padded image
  2. Multiply this new matrix by $(-1)^{x+y}$
  3. Compute the DFT

In code it would be written like this:

spatialToFourier[spatial_, image_] := Module[{padded, centered},
  padded = PadRight[spatial, 2 ImageDimensions[image]];
  centered = MapIndexed[centerPixel, padded, {2}];
  Fourier[centered]
  ]

where spatial in your case is GaussianMatrix[3] and image is lena. centerPixel is given below.

Multiplying by $(-1)^{x+y}$ is actually not necessary but is often done to center the fourier transform. It makes the visualization of the Fourier transform easier to interpret.

Example

lena = ImageResize[ExampleData[{"TestImage", "Lena"}], {121, 121}];
lena = ColorConvert[lena, "Grayscale"]

Mathematica graphics

padImage[image_] := Module[{width, height},
  {width, height} = ImageDimensions[image];
  ImagePad[image, {{0, width}, {height, 0}}]
  ]
padImage[lena]

Mathematica graphics

centerPixel[pixel_, {y_, x_}] := (-1)^(x + y) pixel
centerImage[image_] := Image[MapIndexed[centerPixel, ImageData[image], {2}]]
centeredLena = centerImage[padImage[lena]]

Mathematica graphics

transform[filter_, image_] := Module[{centered, fourier, inverse},
  centered = centerImage[padImage[lena]];
  fourier = Fourier[ImageData[centered]];
  inverse = InverseFourier[filter fourier];
  centerImage[Image[Re[inverse]]]
  ]
spatialFilter = GaussianMatrix[3];
transformed = transform[spatialToFourier[spatialFilter, lena], lena]

Mathematica graphics

cropped = ImageCrop[transformed, ImageDimensions[lena] + Dimensions[spatialFilter], {Right, Bottom}];
cropped = ImageCrop[cropped, ImageDimensions[lena], Center];
GraphicsRow[{
  cropped // ImageAdjust,
  ImageConvolve[lena, spatialFilter, Padding -> Black]
  }, ImageSize -> 2 121]

Mathematica graphics

Community
  • 1
C. E.
  • 70,533
  • 6
  • 140
  • 264
  • Why are you using Re[] instead of Abs[]? – J. M.'s missing motivation May 26 '16 at 10:08
  • @J.M. The IDFT is supposed to be real, the imaginary part is due to numerical inaccuracies. – C. E. May 26 '16 at 10:18
  • In MMA 11.3.0.0 the "transformed" image has Max@ImageData ~ 0.00347, so it is only visible if I use //ImageAdjust. Is MMA 11.3.0.0 doing something differently? Should Max of the transformed image be roughly the same as Max of the Convolved image, or do we lose the scaling by this method? (I ask because I'm trying to see the impact on absolute signal levels after filtering, so really the Total signal should be conserved). – DrBubbles Feb 21 '19 at 19:01
  • @DrBubbles You are right, thank you for pointing this out. I edited the answer accordingly. Unfortunately, I do not remember how this worked in older versions. It seems likely that I made a mistake when originally writing this post. I'm not aware of anything that would have made the scaling work then but not now. – C. E. Feb 21 '19 at 21:25
  • 1
    @C.E Thanks, I can normalize to the total signal in the original data I guess, although there's likely some trick to doing it correctly, if I find it I'll follow up. – DrBubbles Feb 21 '19 at 21:38