5
$\begingroup$

I have an image that spurious waves texture appears:

Imag test

I want to use 2D FFT to remove this texture, so I applied the following code from Calculation of 2D FFT for an image :

data = ImageData[imag];(*get data*)
{nRow, nCol, nChannel} = Dimensions[data];
d = data[[All, All, 2]];
d = d*(-1)^Table[i + j, {i, nRow}, {j, nCol}];
fw = Fourier[d, FourierParameters -> {1, 1}];

(*adjust for better viewing as needed*)
fudgeFactor = 50;
abs = fudgeFactor*Log[1 + Abs@fw];
Labeled[a = Image[abs/Max[abs], ImageSize -> 300], Style["Magnitude spectrum", 18]]

And the result is:

Magnitude Spectrum

How do I transform the image back using inversefourier? So I could eliminate the spurious texture.

I read lot of questions but no one is very clear about how code inverse transform back to "original" image.

$\endgroup$
7
$\begingroup$

You have to reverse all the steps you took,

ifw = Chop@InverseFourier[fw, FourierParameters -> {1, 1}];
ifw = ifw*(-1)^Table[i + j, {i, nRow}, {j, nCol}];
idata = Map[{#, #, #} &, ifw, {2}];
GraphicsRow[{Image[idata],
  Image[data]}, ImageSize -> 800]

enter image description here

I noticed that the three channels of your data aren't exactly equivalent, though they are for most of the pixels. Just evaluate Length@Select[Flatten[data, 1], UnsameQ @@ # &] and you'll see that there are 78,825 pixels with some differences. Of course, their magnitude spectrum looks visually very similar,

enter image description here

but they aren't numerically identical. So really, you should split the three channels into different matrices, do a FFT on each, perform the subtraction you are going for, then do the inverse FFT, then combine them. Fairly straightforward to generalize the code to do that.

Edit This wasn't in your question, but I thought I'd take a crack at removing the unwanted texture (never done anything like this before, so this is bound to be far from the most efficient way to do it).

It's clear from the magnitude spectrum that the problem comes from the zero frequency terms. Your image has 560 rows and 720 columns. I don't know if the zero comes at n/2 or n/2+1, so I just check the data,

Last@Ordering[Abs@fw[[100]]]
Last@Ordering[Abs@fw[[All, 100]]]
(* 361 *)
(* 281 *)

So we can just replace the 361st column with the average of the two columns nearest to it, and do the same for the 281st row. Perhaps it would be better to just remove these columns and rows, but this is my first attempt.

fw2 = fw;
fw2[[All, 361]] = Mean /@ fw2[[All, {360, 362}]];
fw2[[281]] = Mean@fw2[[{280, 282}]];

Look a the spectrum without these terms,

fudgeFactor = 50;
abs = fudgeFactor*Log[1 + Abs@fw2];
Labeled[a = Image[abs/Max[abs], ImageSize -> 300], 
 Style["Magnitude spectrum", 18]]

enter image description here

Now take the inverse FFT and compare to the original image,

ifw = Chop@InverseFourier[fw2, FourierParameters -> {1, 1}];
ifw = ifw*(-1)^Table[i + j, {i, nRow}, {j, nCol}];
idata = Map[{#, #, #} &, ifw, {2}];
GraphicsRow[{Image[idata], Image[data]}, ImageSize -> 1200]

enter image description here

And here is the code for keeping all three channels,

{nRow, nCol, nChannel} = Dimensions[data];
d = data[[All, All, #]] & /@ Range[nChannel];
idata = Module[{d, fw, fw2, ifw},
     d = #*(-1)^Table[i + j, {i, nRow}, {j, nCol}];
     fw = Fourier[#, FourierParameters -> {1, 1}];
     fw2 = fw;
     fw2[[All, 361]] = Mean /@ fw2[[All, {360, 362}]];
     fw2[[281]] = Mean@fw2[[{280, 282}]];
     ifw = Chop@InverseFourier[fw2, FourierParameters -> {1, 1}];
     ifw = ifw*(-1)^Table[i + j, {i, nRow}, {j, nCol}];
     ifw
     ] & /@ d;

idata = Transpose[idata, {3, 1, 2}];

Image[idata]

enter image description here

Edit A question was asked about how to tell which column or row to remove, and hopefully this should help in that regard:

img = Import["http://i.stack.imgur.com/XvkI1.jpg"];
data = ImageData[img];(*get data*){nRow, nCol, nChannel} = 
 Dimensions[data];
d = data[[All, All, 2]];
d = d*(-1)^Table[i + j, {i, nRow}, {j, nCol}];
fw = Fourier[d, FourierParameters -> {1, 1}];
Manipulate[
 Module[{fw2, fudgeFactor, abs, ifw, idata, message},
  fw2 = fw;
  If[
   remcol,
   fw2[[All, column]] = Mean /@ fw2[[All, {column - 1, column + 1}]]
   ];
  If[
   remrow,
   fw2[[row]] = Mean@fw2[[{row - 1, row + 1}]]
   ];
  fudgeFactor = 50;
  abs = fudgeFactor*Log[1 + Abs@fw2];
  ifw = InverseFourier[fw2, FourierParameters -> {1, 1}];
  message = "";
  If[Max@Abs@Im@ifw/Max@Abs@Re@ifw > 10^-10, 
   message = 
    "Warning, Inverse FFT returned complex results,\n proceeding with \
real part only"];
  ifw = Re@ifw;
  ifw = ifw*(-1)^Table[i + j, {i, nRow}, {j, nCol}];
  idata = Map[{#, #, #} &, ifw, {2}];
  Row[{
    Labeled[Image[abs/Max[abs], ImageSize -> 500],
     Style["Magnitude spectrum", 18]],
    Labeled[Image[idata, ImageSize -> 500],
     Style["Reconstructed image", 18]],
    message}
   ]
  ],
 {{remcol, False, "Remove Column?"}, {True, False}},
 {{remrow, False, "Remove Row?"}, {True, False}},
 {{row, Floor[nRow/2]}, 2, nRow - 1, 1, Appearance -> "Open"},
 {{column, Floor[nCol/2]}, 2, nCol - 1, 1, Appearance -> "Open"},
 ContinuousAction -> False]

enter image description here

$\endgroup$
  • $\begingroup$ It worked very well, thanks Jason B! $\endgroup$ – locometro Jan 21 '16 at 12:27
  • $\begingroup$ How do you now the exactly, columns and lines? $\endgroup$ – locometro Jan 22 '16 at 15:25
  • 1
    $\begingroup$ @locometro, see the edit $\endgroup$ – Jason B. Jan 22 '16 at 16:25
  • 1
    $\begingroup$ @locometro, so when I looked at the spectrum image it was clear there was a large stripe down the middle of both axes (and therefore corresponds to zero frequency). I assumed this was what you were trying to remove, a general noise term. To find out exactly which column corresponded to this stripe I used trial and error, doing listplots of the carious rows or columns. $\endgroup$ – Jason B. Jan 22 '16 at 20:22
  • 1
    $\begingroup$ Then when I took those away I noticed that the image looked like a general noise term had been taken out. This is literally my first time with image processing as well. Try the last bit of code I put up, it can help to decide which terms to remove. You can email me at Jason.biggs@mpsd.mpg.de if you like $\endgroup$ – Jason B. Jan 22 '16 at 20:26

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.