FFT phase ?

[imaggie]

Hi,

I am having some doubts as to whether -fft is producing the correct phase data. (Or perhaps documentation is incomplete as to what it is _supposed_ to produce).

I'm looking at the DFT of simple straight line segment that can be used to simulate a motion blur:



I reproduce the above image using IM then do the DFT in phase/magnitude and get the two layers as png files.

Code: Select all

length=15;
name=motion$length;
convert -size  ${ww}x${ww}  xc:black -fill white \
-draw "line  $(($half-$length/2)),$half $(($half+$length/2)),$half" -alpha off \
 ${name}.png

mean=`convert ${name}.png -format "%[mean]" info:`

mask=$name.png

ww=`identify -format "%w" $mask`
half=$(($ww/2))

# same as above to get four images Real,Imaginary
convert \( $mask  -roll -${half}-${half} -fft -level 0x$mean \) \
 \( cameraman.png -fft \) \
  cameraman_${name}-fftMP.png

convert \( $mask  -roll -${half}-${half} -fft  \) \
  ${name}-fftMP.png

The magnitude looks credible as basically sinc-like of expected proportions. However I'm rather unsure about the phase:


Latest on-line doc tells me the following:

http://www.imagemagick.org/script/comma ... 0na4pidop4
Phase values nominally range from 0 to 2*π, but for non-HDRI compilations of ImageMagick, the phase image is scaled to span the full dynamic range.

I'm using a HDRI build , so this statement (or it's implied inverse) is not true. The above plot is not in the range 0-2*pi !

Code: Select all

 convert -version
Version: ImageMagick 6.7.4-0 2012-01-03 Q16 http://www.imagemagick.org
Copyright: Copyright (C) 1999-2011 ImageMagick Studio LLC
Features:   HDRI  

I have several problems here. First it seems clear that both non_HDRI and HDRI are being scaled , or perhaps the layer is not being scaled but is scaled in writing, I don't know.

Second , this image has only three values: black and white and 50%grey. This would seem an improbably result unless this is just near zero noise being scaled up , even then it's dubious.

Finally , unless I am confusing thing I have other sources that indicate the imaginary part is a cosine and that makes the phase a smoothly rotating angle not this staccato all or nothing.

http://mathworld.wolfram.com/SincFunction.html says : {sinc(x),cos(x)}
also I have a locally built DFT routine released by arachnoid, that work very well on 1-D (audio) signals.
http://www.arachnoid.com/signal_processing/fft.html

Using the latter to transform a rectangular step I get pure cosine in the imaginary.

Using IM , I get a black square.

Code: Select all

convert \( $mask  -roll -${half}-${half} +fft  \)   ${name}+fft.png

All this suggests that some assumptions about the result are being made as "shortcuts" in the code. I'm wondering if they are correct or some undocumented modifications are being imposed of the DFT.

This may be "small but acceptable" errors and maybe not . When chasing out ringing in FFT processes you need to be nit-pickingly precise. Undocumented features (if that is the case) and incorrect statements in the doc do not help help diagnose problems.

Any clarifications as to what extent all this is a bug or a feature would be helpful.

regards.

[fmw42]

No roll is necessary for the image, unless you are using it as the filter function. Then it has to be rolled so as centered at 0,0. In my scripts such as fftconvol fftdeconvol, I do the roll of the spatial filter (only), not the image, internal to the script. So it is not necessary to do the roll externally as most people do not understand that.


I get the following for the phase of your image (without doing the roll). But you have to remember that the phase is arctan(im/re). Wolfram's extension of sinc to complex plane is extending to real imaginary and that is not the definition of phase. You are assuming that sinc is the magnitude but it is not. It is the real part of the transform. The magnitude would be sqrt(re^2+im^2). So when the phase is scaled from -pi ...pi to 0 ... 2pi, they both are positive with no negative values. Wofram's extension of the sinc to complex domain has both positive and negative values. Now I have to believe that his extension to the complex domain correct, but does not apply to the magnitude phase representation of the motion15.png straight line image when doing a real<->complex (mag phase) transform. It would appear to apply to a real<->complex (real imaginary) transform. However as you know Wolfram also suggests that the imaginary component adds little to the issue of motion blur deblur using the real imaginary format for the filter, where the sinc is the real part. In this case then the magnitude would be sqrt(re^2)=abs(re). However it would be interesting to try adding the cosine imaginary part to the sinc real part for the blur/deblur filter in the frequency domain.

With regard to the IM FFT, it is based on the well known standard FFTW code (used as a delegate library) and use extensively by many people. So I have doubts it is incorrect. The only additions that IM added were to convert the code from the half transform provided by FFTW (as half is all that is really needed) to unfold it to the full transform and to scale the phase from -pi ... pi to 0 ... pi and then to 0 ... 1 as all IM non-HDRI images are processed in that range and then multiplied by quantum range for saving. The decision was made a long time ago to make the HDRI phase the same scaling rather than leave it as -pi ... pi. It was that way for a while and then changed to be consistent.

See
http://www.fftw.org/
http://paulbourke.net/miscellaneous/dft/
http://paulbourke.net/miscellaneous/imagefilter/


convert motion15.png -fft motion15_fft_mp.png

Magnitude after applying log (log of absolute value of sinc):

convert motion15_fft_mp-0.png -evaluate log 1000000 motion15_fft_mp-0_l1000000.png




Phase:

[imaggie]

No roll is necessary for the image, unless you are using it as the filter function. Then it has to be rolled so as centered at 0,0. In my scripts such as fftconvol fftdeconvol, I do the roll of the spatial filter (only), not the image, internal to the script. So it is not necessary to do the roll externally as most people do not understand that.

A phase shift won't affect the real or magnitude of FT but it does shift the phase around a bit. The chequer board effect seen in the phase is a bit odd. I also note that your result is black on white only, you've no grey. I'll have to think about why all that should be.

I get the following for the phase of your image (without doing the roll). But you have to remember that the phase is arctan(im/re). Wolfram's extension of sinc to complex plane is extending to real imaginary and that is not the definition of phase. You are assuming that sinc is the magnitude but it is not. It is the real part of the transform. The magnitude would be sqrt(re^2+im^2).

No I'm not thinking the mag is sync, sinc has negative lobes. Clearly mag is not the same as real.

So when the phase is scaled from -pi ...pi to 0 ... 2pi, they both are positive with no negative values. Wofram's extension of the sinc to complex domain has both positive and negative values. Now I have to believe that his extension to the complex domain correct, but does not apply to the magnitude phase representation of the motion15.png straight line image when doing a real<->complex (mag phase) transform. It would appear to apply to a real<->complex (real imaginary) transform.

Well what is this real<->complex (mag phase) transform as opposed to real<->complex (real imaginary) transform ?
Whatever way you go about it the result has to be the same thing expressed in different coordinates. I'm not suggesting Wolfram {sinc(x), cos(x)} can be read as mag/phase but it MUST be the same once converted. And I'm not convinced that phase relationship between sinc and cos can be this odd binary result we are seeing.

However as you know Wolfram also suggests that the imaginary component adds little to the issue of motion blur deblur using the real imaginary format for the filter, where the sinc is the real part. In this case then the magnitude would be sqrt(re^2)=abs(re). However it would be interesting to try adding the cosine imaginary part to the sinc real part for the blur/deblur filter in the frequency domain.

Since you seem to think I'm not misreading this I'll give it a try. The " imaginary component adds little to the issue" may be the little that is causing defects.

With regard to the IM FFT, it is based on the well known standard FFTW code (used as a delegate library) and use extensively by many people. So I have doubts it is incorrect. The only additions that IM added were to convert the code from the half transform provided by FFTW (as half is all that is really needed) to unfold it to the full transform and to scale the phase from -pi ... pi to 0 ... pi and then to 0 ... 1 as all IM non-HDRI images are processed in that range and then multiplied by quantum range for saving. The decision was made a long time ago to make the HDRI phase the same scaling rather than leave it as -pi ... pi. It was that way for a while and then changed to be consistent.

So you are able to confirm what I thought , that the doc is broken on this issue and never got updated when this situation was unified.

That also has implications for deblur. I have a draft email for you about that but I was holding off to be more certain of what was going in here.

See
http://www.fftw.org/
http://paulbourke.net/miscellaneous/dft/
http://paulbourke.net/miscellaneous/imagefilter/

[fmw42]

What part of the doc is broken -- where specifically?

[fmw42]

I don't want to say that IM FFT cannot be broken. I have had problems before. But I will have to check my standard test image to see. Also you might like to look at your image results using other software. ImageJ has some good FFT tools and also see some of their plugins for more complex deblurring techniques than I have developed.

Here are my standard tests:

Input:


Mag/Phase transform:

convert square31.png -fft square31_hdri_fft_mp.png
convert square31_hdri_fft_mp-0.png -evaluate log 10000 square31_hdri_fft_mag_l1000000.png

Log mag:


Log mag profile (absolute value of sinc):
profile -r 128 -f 1 -l off square31_hdri_fft_mag_l1000000.png square31_hdri_fft_mag_l1000000_profile.gif


Phase:



The above are the same results I get from using FFTJ or ImageJ


Real/Imag transform:

convert square31.png +fft square31_hdri_fft_ri.pfm

Note values contain negatives. So the best way to check is to convert to mag/phase


convert square31_hdri_fft_ri.pfm -fx "hypot(u,v)" square31_hdri_fft_ri2mag.png
convert square31_hdri_fft_ri2mag.png -evaluate log 1000000 square31_hdri_fft_ri2mag_l1000000.png
convert square31_hdri_fft_ri.pfm -fx "0.5*atan2(v,u)+0.5" -auto-level square31_hdri_fft_ri2phase.png


Log mag from real/imaginary:



Phase





So apart from some left/right edge "noise" in the phase, the mag and phase from real/imaginary look just like the mag and phase from the direct transform.

So I would have to say that basically the FFT transfrom either mag/phase or real/imag is working correctly.

I am not sure where the "noise" in the phase from real/imag is coming from unless some kind of wrap around or precision issues?

All the above were tested using IM 6.7.4.4 Q16 HDRI on Mac OSX Snow Leopard.

Please feel free to do testing against ImageJ or other programs.

[imaggie]

What part of the doc is broken -- where specifically?

specifically the bit I noted above about non-HDRI being different. In fact non-HDRI seems to be as stated but the implication that HDRI is different has apparently not been the case " for a long time".

When I went to the doc the try to understand what was going on , being told to expect something that was not happening only added to the confusion.

It's a detail that obviously got missed when the code changed.

I still don't quite get why there's no imaginary part is not being produced. Broken or by design ?

[fmw42]

What part of the doc is broken -- where specifically?

specifically the bit I noted above about non-HDRI being different. In fact non-HDRI seems to be as stated but the implication that HDRI is different has apparently not been the case " for a long time".

I still don't quite get why there's no imaginary part is not being produced. Broken or by design ?

1) What exact page are you looking at for HDRI or where you think something needs changing/adding?

2) The imaginary part is produced by +FFT along with the real component. -FFT produces magnitude and phase.

[imaggie]

quote and link in first post. According to what I see and what you confirmed, it seems HDRI and non_HDRI are now scaled the same and it's not 0-2*pi.

[fmw42]

quote and link in first post. According to what I see and what you confirmed, it seems HDRI and non_HDRI are now scaled the same and it's not 0-2*pi.

That is correct. They are both scaled to 0 to quantumrange so that the phase if easily visible without scaling. Internally they are used properly with phase between 0 and 2*pi (or -pi to +pi). Magnitude, real and imaginary components are never scaled.

I have updated the docs on the options page, but it may take a few hours or perhaps a day until the changes become visible.

[fmw42]

Well what is this real<->complex (mag phase) transform as opposed to real<->complex (real imaginary) transform ?
Whatever way you go about it the result has to be the same thing expressed in different coordinates. I'm not suggesting Wolfram {sinc(x), cos(x)} can be read as mag/phase but it MUST be the same once converted. And I'm not convinced that phase relationship between sinc and cos can be this odd binary result we are seeing.

The fft converts a real image to a complex image in the form of real and imaginary components (the definition of complex). IM can output that by using +fft in HDRI mode. IM can output mag and phase where the real and imaginary are processed to form mag=sqrt(re^2+im^2) and phase =arctan2(im,re).

I am still puzzled by Wolframs complex sinc,cos representation.

I took my square31.png image and converted it to real and imaginary using +fft. The magnitude and phase are shown above, where the mag is the absolute value of the sinc. As I said above, converting the real and imaginary to mag phase reproduces the results from using -fft.

However, I tried various ways to scale and log process the real and imaginary components, but surprisingly the real component looks nothing like a 2D sinc. It has almost no shape whatever that I can tell but a bit brighter in the middle. The imaginary component looks nothing like a cos either and again has no obvious shape. Perhaps I am making some mistake in trying to get a viewable image. However, the real and imaginary components must be correct as they produce the correct mag/phase pair and the mag is the absolute value of the sinc. Also I get the same results for mag/phase pair using other software.

The other issue I am still puzzled about is how these real and imaginary components relate to the fact that in the 1D (analog) fourier transform, the rect and sinc are transform pairs. No one ever states anything about the imaginary part. See table in http://en.wikipedia.org/wiki/Fourier_transforms.

All I can say is that somehow(???) the real and imaginary components of the FFT of the square in my image is somehow equivalent to having a real sinc and zero imaginary component when using them as deblurring filters for motion blur. This equivalence works in practice as shown by my ability to deblur lens defocus (the jinc analogy to the sinc) in real images (not computer simulated blur). This is also not something I claim to have invented. Using the real (not-complex) sinc as the filter in the frequency domain has also been used by others to correct for motion blur. However, at this time I do not understand the mathematical connection between the real sinc and the real/imaginary components of the FFT of a rectangle (square) as to how they both seem to be different but produce the same filtering results.

[imaggie]

Code: Select all

The other issue I am still puzzled about is how these real and imaginary components relate to the fact that in the 1D (analog) fourier transform, the rect and sinc are transform pairs. No one ever states anything about the imaginary part. See table in http://en.wikipedia.org/wiki/Fourier_transforms.

Oh please, don't quote WP as a reference for anything. It's like watching a Hollywood film "based on a true story" and imagining you are learning history.

For the 1D case I don't think you are correct. As I said elsewhere , using the arachnoid DFT software I get precisely that in transforming a rect : as pure cosine imaginary component. " No one ever states " is hearsay , not science. I think this is a popular misconception. I was taught that the FT of a rect was a sinc but was never explicitly told that there is an imaginary part that is zero. Until I started using FT and looking more closely, I had assumed that was the case.

In fact I doubt that any real only function will generate a real only FT.

However, the real and imaginary components must be correct as they produce the correct mag/phase pair

I am not aware of any fast DFT that actually operates in polar co-ordinates (not to say it can't be done). My guess is that either IM converts when you use -fft or FFTW does it internally. By using +fft or -fft you are simply stating in what form you wish to have the output. That the two are consistent does prove or confirm anything is correct.

However, at this time I do not understand the mathematical connection between the real sinc and the real/imaginary components of the FFT of a rectangle (square) as to how they both seem to be different but produce the same filtering results.

Could you post a code snip that demonstrates that?

[imaggie]

quote and link in first post. According to what I see and what you confirmed, it seems HDRI and non_HDRI are now scaled the same and it's not 0-2*pi.

That is correct. They are both scaled to 0 to quantumrange so that the phase if easily visible without scaling. Internally they are used properly with phase between 0 and 2*pi (or -pi to +pi). Magnitude, real and imaginary components are never scaled.

OH What?! More undocumented "features" . No wonder I can't follow what IM is doing. Whenever I try to pull out an intermediary step to see what is going on it gets "converted" so I can see it.

This is like trying to debug a program but when you print out a value the interpreter decides that in some cases it would be better to scale to the range 0-65335 before printing it.

So how can I access this "internal", not mucked about with, transform?

Thx

[imaggie]

now If I take a x=0 or y=0 section through that plot , it's about what I would expect: an even cyclic variation of phase. The light / dark discontinuity is in fact the -pi , +pi wrap around.

I was expecting something more like that from the FT of a 1x15 rect.

[fmw42]

Oh please, don't quote WP as a reference for anything. It's like watching a Hollywood film "based on a true story" and imagining you are learning history.

So are you saying that it is wrong. A lot of experts have contributed to Wikipedia. So there must be something to it. I cannot explain to you but using the real sinc or jinc to recover the image from motion blur and defocus (without an imaginary component) works.

I am not aware of any fast DFT that actually operates in polar co-ordinates (not to say it can't be done). My guess is that either IM converts when you use -fft or FFTW does it internally. By using +fft or -fft you are simply stating in what form you wish to have the output. That the two are consistent does prove or confirm anything is correct.

FFT or DFT etc, do not use polar coordinates and I never said so. The fact that the jinc is circular symmetric comes from the fact that the spatial transform or PSF is a circle for defocus. This has nothing to do with the workings of the FFT.

Could you post a code snip that demonstrates that?

No, my conclusions are based on getting the same kind of results from simulations and real defocus images. If I use the fftdeconvol using a circle filter or I use the cameradeblur using a jinc, they both correct the defocus. The fftdeconvol has to deal however with the fact that the digital circle is not perfect, especially for small circles. Nevertheless, the both deblur the image. Both methods use real-imaginary FFT components of the image. fftdeconvol use real/imaginary components for the fft of the circle filter. cameradeblur uses on the jinc as the real part.

So how can I access this "internal", not mucked about with, transform?

As far as I know, you cannot in mag phase, only from real/imaginary in HDRI mode. You could use those to regenerate mag phase without any scaling if you want.

OH What?! More undocumented "features" . No wonder I can't follow what IM is doing. Whenever I try to pull out an intermediary step to see what is going on it gets "converted" so I can see it.

I have edited the docs to bring them up to speed with the current IM FFT. If you find any other issues or have questions, then please notify us and we will try to fix them as quickly as possible.

If you are unhappy with IM's approach to FFT, you are free to use other software. I have already suggested two -- ImageJ and FFTJ. I am sure there are lots of other tools out there.

[fmw42]

Here are the center row profiles of the mag/phase images from the square31.png image.

mag:





phase:





But I really don't know how to interpret the phase profile. It seems to have very high frequency oscillations within slower osciallations, the later of which seem to match the cycles in the mag profile.