[Numpy-discussion] FFT's & IFFT's on images

[Robert Kern]

On Wed, Jul 2, 2008 at 16:33, Mike Sarahan <msarahan at gmail.com> wrote:
> Hi,
>
> I'm trying to do phase reconstruction on images which involves switching
> back and forth between Fourier space and real space.  I'm trying to test
> numpy (& scipy, for that matter) just to see if I can go back and forth.
> After an FFT/iFFT, the resulting image is garbage.  I'm using
> numpy.fft.fftn, but I've also tried fft2, rfftn, rfft2, and the
> corresponding inverse FFT's.
>
> >From looking at the matrices, it appears to be creating complex
> components that aren't in the matrix prior to any FFT's.

Are the all of the order 1e-14 or so? These are expected.
Double-precision floating point operations each have about 1e-16
relative error. For an FFT on an image with values in [0..255], you
will get absolute errors about 1e-16 to 1e-13 reconstructing the
original array.

> Real fft's
> seem to add some small component to each value (<1).

Exactly how large are these? If they are the same 1e-14 or so, then
everything is working as expected.

For example, here is what I get when doing this with fft2 and rfft2 on
Images/lena.gif from the PIL's source tree.

In [23]: from scipy.misc import imread

In [24]: !ls
IPython system call: ls
courB08.bdf  courB08.pbm  courB08.pil  lena.gif  lena.jpg  lena.ppm

In [25]: lena = imread('lena.gif')

In [26]: lena.shape, lena.dtype
Out[26]: ((128, 128), dtype('uint8'))

In [27]: from numpy import fft

In [28]: flena = fft.fft2(lena)

In [29]: lena2 = fft.ifft2(flena)

In [30]: lena
Out[30]:
array([[100, 100, 100, ..., 197, 100, 104],
       [124,  62,  58, ...,  65,  21,  24],
       [207, 124, 104, ...,  24,  59,  59],
       ...,
       [216,  23,  84, ...,  84, 216, 193],
       [169, 216, 169, ..., 216, 118,  26],
       [ 59,  59,  59, ..., 193,  87, 128]], dtype=uint8)

In [31]: lena2
Out[31]:
array([[ 100. -1.02418074e-14j,  100. -4.73237530e-15j,
         100. -1.80966353e-14j, ...,  197. -1.54082828e-14j,
         100. +1.42941214e-14j,  104. -3.92740898e-14j],
       [ 124. +2.68025506e-14j,   62. -1.36449897e-14j,
          58. +6.41674693e-15j, ...,   65. -3.95411363e-14j,
          21. +4.25895005e-14j,   24. -1.08962931e-14j],
       [ 207. -4.84334794e-14j,  124. +4.98905975e-14j,
         104. -2.94209102e-14j, ...,   24. -5.80408469e-14j,
          59. +1.02973186e-14j,   59. +2.53408902e-14j],
       ...,
       [ 216. +5.40427922e-14j,   23. -1.46100124e-13j,
          84. -3.71689877e-15j, ...,   84. +8.35485795e-14j,
         216. +1.27404625e-13j,  193. +1.54715424e-14j],
       [ 169. -3.69981823e-14j,  216. +1.79161744e-14j,
         169. -3.98847622e-14j, ...,  216. +7.45197822e-14j,
         118. -1.14630527e-14j,   26. -7.45791820e-14j],
       [  59. -3.23731472e-16j,   59. -5.38406684e-14j,
          59. +7.23899628e-15j, ...,  193. -1.76964932e-14j,
          87. -1.37792130e-14j,  128. +4.76010179e-14j]])

In [32]: (abs(lena2 - lena) < 1e-12).all()
Out[32]: True

In [34]: rflena = fft.rfft2(lena)

In [35]: rlena2 = fft.irfft2(rflena)

In [36]: (abs(rlena2 - lena) < 1e-12).all()
Out[36]: True

In [37]: rlena2
Out[37]:
array([[ 100.,  100.,  100., ...,  197.,  100.,  104.],
       [ 124.,   62.,   58., ...,   65.,   21.,   24.],
       [ 207.,  124.,  104., ...,   24.,   59.,   59.],
       ...,
       [ 216.,   23.,   84., ...,   84.,  216.,  193.],
       [ 169.,  216.,  169., ...,  216.,  118.,   26.],
       [  59.,   59.,   59., ...,  193.,   87.,  128.]])


Note that when printing arrays, floating point numbers are rounded to
11 decimal places or so for convenience, so the teeny-tiny errors
aren't seen in the real parts of the reconstructed images.

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless
enigma that is made terrible by our own mad attempt to interpret it as
though it had an underlying truth."
 -- Umberto Eco