DSPRelated.com
Forums

Confusion about 2D DFTs

Started by Michel Rouzic August 3, 2008
Here's my problem. I have an image to be processed (well, deconvolved
pretty much), and a star-shaped impulse response I want to deconvolve
it with. The impulse response is symmetric in both the vertical and
horizontal axis and pretty much centred, but not very precisely. Of
course both the image and the impulse response are appropriately zero-
padded as to avoid any circular convolution.

My idea is to transform the IR into a frequency-domain (magnitude
only) representation using something like a 2D FFT, hard limit the
frequency-domain magnitude bins under a threshold value to make sure
that the weaker frequency components don't get amplified more than
they should, and then perform the deconvolution.

Well at this point I guess I could simply do a f(x) = 1/x on all the
magnitude bins, zero the phase (as to make the IR centered around 0,
0), IFFT and convolve in the space-domain, or I could FFT the image to
be deconvolved and multiply (or rather divide) with the frequency-
domain magnitude of the hard-limited impulse response. The problem is
that I don't understand how to get the magnitude of the impulse
response FFT (or similar).

There's something I don't understand about 2D DFTs. Just like the DFT
of a 1D real signal is represented by an array of complex bins
representing from DC to the Nyquist frequency (omitting the redundant
negative frequencies) I would expect a 2D DFT of a (real) image to be
a complex array ranging from each axis' DC to each axis' Nyquist
frequency (although I realise this is surely not possible since it
would take twice less space to represent a real image under that form)

Instead, it seems that 2D DFTs are achieved by performing 1D DFTs on
each row, and then performing 1D DFTs on each column. If I understand
correctly, the result is a complex array representing from -Nyquist to
+Nyquist on each axis, and is symmetrical on one axis due to the fact
that the first series of 1D DFTs is performed on real data, as the
rest is on complex data and therefore has non-redundant negative
frequency components. My problem with that is that I don't understand
why we should "privilege" an axis by processing it first, and
therefore obtain a different result that if we have processed the
other axis first.

Also, I don't understand what to do with that. If I'm not mistaken, to
obtain the magnitude of my IR in the frequency domain, I just need to
do a sqrt(real� + imaginary�) on each of my complex bins, correct? So
then I'm left with.. a magnitude that is non-symmetrical around DC on
only one axis? What do I do with that, I mean, why should it have
negative frequency components that are not the same as the positive
ones, and most importantly, why on only one (arbitrarily chose) axis??

I should also probably mention that I'm implementing that using FFTW 3
(no Matlab to toy around with), and that its documentation is pretty
poor on the very topic of my problem. It seemingly doesn't tell what
the output matches to depending on the transform. For example I
noticed that the output is quite different between a 2D Discrete
Hartley Transform and a 2D DFT (and also a 2D DCT, but I know why that
one is different, I at least understand that, but the documentation
won't say anything about it), why? Disclaimer : I don't really know
what the DHT does, but I find it striking that its output (as a real
image, so I take it it's in their "half complex" format that goes DC
Re1 Re2 Re3  Nyquist Im3 Im2 Im1) is made of only two parts (which I
guess would come as one if it were to be put as a complex array as it
may (or may not? the documentation doesn't seem to discuss this) be in
"half complex" format).

Thanks in advance
Michel Rouzic wrote:
> [snip] > I should also probably mention that I'm implementing that using FFTW 3 > (no Matlab to toy around with),[snip]
Can't help with the 2D aspect but have you looked at Scilab www.scilab.org ?
On Aug 3, 2:52&#4294967295;pm, Richard Owlett <rowl...@atlascomm.net> wrote:
> Michel Rouzic wrote: > > [snip] > > I should also probably mention that I'm implementing that using FFTW 3 > > (no Matlab to toy around with),[snip] > > Can't help with the 2D aspect but have you looked at Scilabwww.scilab.org?
A bit, but I'm not looking for anything like Matlab really ;-) maybe I could benefit from getting it and learning it but so far I like doing even prototyping in C. Which is convenient if you can't be bothered to rewrite your fully working Matlab code into C ;-)
On 3 Aug, 15:17, Michel Rouzic <Michel0...@yahoo.fr> wrote:
> Here's my problem. I have an image to be processed (well, deconvolved > pretty much), and a star-shaped impulse response I want to deconvolve > it with. The impulse response is symmetric in both the vertical and > horizontal axis and pretty much centred, but not very precisely. Of > course both the image and the impulse response are appropriately zero- > padded as to avoid any circular convolution. > > My idea is to transform the IR into a frequency-domain (magnitude > only) representation using something like a 2D FFT, hard limit the > frequency-domain magnitude bins under a threshold value to make sure > that the weaker frequency components don't get amplified more than > they should, and then perform the deconvolution.
Deconvolution is tricky business at the best of times. You might want to consult the 2008 edition of the Gozales and Woods book on image processing to see a number of approaches and some comparisions of results.
> Well at this point I guess I could simply do a f(x) = 1/x on all the > magnitude bins, zero the phase (as to make the IR centered around 0, > 0), IFFT and convolve in the space-domain, or I could FFT the image to > be deconvolved and multiply (or rather divide) with the frequency- > domain magnitude of the hard-limited impulse response. The problem is > that I don't understand how to get the magnitude of the impulse > response FFT (or similar).
OK. The 2D FFT is the problem here. Let's focus on that and forget about deconvolution for now.
> There's something I don't understand about 2D DFTs. Just like the DFT > of a 1D real signal is represented by an array of complex bins > representing from DC to the Nyquist frequency (omitting the redundant > negative frequencies) I would expect a 2D DFT of a (real) image to be > a complex array ranging from each axis' DC to each axis' Nyquist > frequency (although I realise this is surely not possible since it > would take twice less space to represent a real image under that form) > > Instead, it seems that 2D DFTs are achieved by performing 1D DFTs on > each row, and then performing 1D DFTs on each column. If I understand > correctly, the result is a complex array representing from -Nyquist to > +Nyquist on each axis, and is symmetrical on one axis due to the fact > that the first series of 1D DFTs is performed on real data, as the > rest is on complex data and therefore has non-redundant negative > frequency components. My problem with that is that I don't understand > why we should "privilege" an axis by processing it first, and > therefore obtain a different result that if we have processed the > other axis first.
If we start with the computations, then you are right in that the concept is to first perform 1D FFTs in one direction and then 1D FFTs in the other direction (although it may be possible to implement the 2d FFT in other ways; I don't know). The end result is, however, the same regardless of the order the FFTs are computed. As for the structure of the spectrum, it will become a little bit different than you expect. Assume, for the sake of argument, that you compute the FFTs in the x direction first: f(x,y) --[FFT_x]--> f(u,y) Here f(x,y) is the real-valued image, FFT_x is the FFT along the x dimension, and u is the 'spatial x frequency'. Since the image f(x,y) is real-valued, each line f(u) in the image is complex-valued and conjugate symmetric. No surprises. Now, once we compute the FFT in the other dimension, f(u,y) --[FFT_y]--> f(u,v) the input to the y-direction FFT, the columns of f(u,y), are complex-valued, not real-valued. This means that the columns f(v) in the end spectrum are no longer symmetric. Remember, Nyquist's theorem f < Fs/2 only applies to real-valued signals because of the conjugate symmetry of Euler's equations, cos(x) = 1/2 (exp(jx)+exp(-jx)) sin(x) = 1/j2(exp(jx)-exp(-jx)) (are the signs correct here? i don't remember...) Since complex-valued signals don't have the conjugate symmetric spectrum, the complex valued form of Nyquist's theorem becomes f < Fs. The net effect is that the 2D spectrum is cylindrical conjugate symmetric around origo (view with fixed-width font): ^ v | +----------------+----------------+ Fy/2 | | | | | | | | x (u,v) | 0 -+----------------+----------------+--> | x | | u | (-u,-v) | | | | | +----------------+----------------+ -Fy/2 -Fx/2 0 Fx/2 The conjugate mirror of a sinusoidal at (u,v) will occur at (-u.-v) in a 2D spectrum centered at (0,0).
> Also, I don't understand what to do with that. If I'm not mistaken, to > obtain the magnitude of my IR in the frequency domain, I just need to > do a sqrt(real&#4294967295; + imaginary&#4294967295;) on each of my complex bins, correct? So > then I'm left with.. a magnitude that is non-symmetrical around DC on > only one axis? What do I do with that, I mean, why should it have > negative frequency components that are not the same as the positive > ones, and most importantly, why on only one (arbitrarily chose) axis??
These things are explained in some detail in the Gonzales and Woods book. The trick is to use circular symmetry around origo. Rune

Rune Allnor wrote:
> If we start with the computations, then you are right in that > the concept is to first perform 1D FFTs in one direction and > then 1D FFTs in the other direction (although it may be possible > to implement the 2d FFT in other ways; I don't know). > > The end result is, however, the same regardless of the order the > FFTs are computed.
Ha, I don't understand that part. How can it be the same thing considered that, if we start from a real-only 2D image, and that therefore the result should be symmetrical for the first processed axis?
> As for the structure of the spectrum, it will become a little bit > different than you expect. Assume, for the sake of argument, that > you compute the FFTs in the x direction first: > > f(x,y) --[FFT_x]--> f(u,y) > > Here f(x,y) is the real-valued image, FFT_x is the FFT along > the x dimension, and u is the 'spatial x frequency'. > > Since the image f(x,y) is real-valued, each line f(u) in the > image is complex-valued and conjugate symmetric. No surprises. > > Now, once we compute the FFT in the other dimension, > > f(u,y) --[FFT_y]--> f(u,v) > > the input to the y-direction FFT, the columns of f(u,y), > are complex-valued, not real-valued. This means that the > columns f(v) in the end spectrum are no longer symmetric. > Remember, Nyquist's theorem > > f < Fs/2 > > only applies to real-valued signals because of the > conjugate symmetry of Euler's equations, > > cos(x) = 1/2 (exp(jx)+exp(-jx)) > sin(x) = 1/j2(exp(jx)-exp(-jx)) > > (are the signs correct here? i don't remember...) > > Since complex-valued signals don't have the conjugate > symmetric spectrum, the complex valued form of Nyquist's > theorem becomes > > f < Fs. > > The net effect is that the 2D spectrum is cylindrical > conjugate symmetric around origo (view with fixed-width > font): > ^ v > | > +----------------+----------------+ Fy/2 > | | | > | | | > | | x (u,v) | > 0 -+----------------+----------------+--> > | x | | u > | (-u,-v) | | > | | | > +----------------+----------------+ -Fy/2 > -Fx/2 0 Fx/2 > > The conjugate mirror of a sinusoidal at (u,v) will > occur at (-u.-v) in a 2D spectrum centered at (0,0). > > > Also, I don't understand what to do with that. If I'm not mistaken, to > > obtain the magnitude of my IR in the frequency domain, I just need to > > do a sqrt(real&#65533; + imaginary&#65533;) on each of my complex bins, correct? So > > then I'm left with.. a magnitude that is non-symmetrical around DC on > > only one axis? What do I do with that, I mean, why should it have > > negative frequency components that are not the same as the positive > > ones, and most importantly, why on only one (arbitrarily chose) axis?? > > These things are explained in some detail in the Gonzales and Woods > book. > The trick is to use circular symmetry around origo.
Ha. thanks for the explanation, that clears up a few things about the theoretical background. But obviously I'm not going to order a $120 book to solve this particular problem (although I'll surely consider the book again if it appears to cover the image processing problems I have yet to cover) so would you mind elaborating a bit on that please? Here's what I figure should be the way to do it, based on your explanations of 2D FFTs : I should do my hard limiting on all magnitude bins alike regardless of which frequencies they represent, then divide the image's FFT's complex bins with the aforementioned magnitude. Considered both arrays would be the same size to begin with (since they would appropriately padded) their spatial frequencies will map so I don't even really have to wonder what I'm looking at and simply process bins all the same. Does it sound correct?
On 3 Aug, 19:29, Michel Rouzic <Michel0...@yahoo.fr> wrote:
> Rune Allnor wrote: > > If we start with the computations, then you are right in that > > the concept is to first perform 1D FFTs in one direction and > > then 1D FFTs in the other direction (although it may be possible > > to implement the 2d FFT in other ways; I don't know). > > > The end result is, however, the same regardless of the order the > > FFTs are computed. > > Ha, I don't understand that part. How can it be the same thing > considered that, if we start from a real-only 2D image, and that > therefore the result should be symmetrical for the first processed > axis?
The 2D DFT rotational symmetric around origo, not around any line in (u,v) domain.
> > As for the structure of the spectrum, it will become a little bit > > different than you expect. Assume, for the sake of argument, that > > you compute the FFTs in the x direction first: > > > f(x,y) --[FFT_x]--> f(u,y) > > > Here f(x,y) is the real-valued image, FFT_x is the FFT along > > the x dimension, and u is the 'spatial x frequency'. > > > Since the image f(x,y) is real-valued, each line f(u) in the > > image is complex-valued and conjugate symmetric. No surprises. > > > Now, once we compute the FFT in the other dimension, > > > f(u,y) --[FFT_y]--> f(u,v) > > > the input to the y-direction FFT, the columns of f(u,y), > > are complex-valued, not real-valued. This means that the > > columns f(v) in the end spectrum are no longer symmetric. > > Remember, Nyquist's theorem > > > f < Fs/2 > > > only applies to real-valued signals because of the > > conjugate symmetry of Euler's equations, > > > cos(x) = 1/2 (exp(jx)+exp(-jx)) > > sin(x) = 1/j2(exp(jx)-exp(-jx)) > > > (are the signs correct here? i don't remember...) > > > Since complex-valued signals don't have the conjugate > > symmetric spectrum, the complex valued form of Nyquist's > > theorem becomes > > > f < Fs. > > > The net effect is that the 2D spectrum is cylindrical > > conjugate symmetric around origo (view with fixed-width > > font): > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;^ v > > &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;| > > &#4294967295; &#4294967295; +----------------+----------------+ Fy/2 > > &#4294967295; &#4294967295; | &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;| &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;| > > &#4294967295; &#4294967295; | &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;| &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;| > > &#4294967295; &#4294967295; | &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;| &#4294967295; &#4294967295;x (u,v) &#4294967295; &#4294967295; | > > &#4294967295;0 -+----------------+----------------+--> > > &#4294967295; &#4294967295; | &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;x &#4294967295; &#4294967295; | &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;| &#4294967295;u > > &#4294967295; &#4294967295; | &#4294967295; &#4294967295; &#4294967295;(-u,-v) &#4294967295; | &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;| > > &#4294967295; &#4294967295; | &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;| &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;| > > &#4294967295; &#4294967295; +----------------+----------------+ -Fy/2 > > &#4294967295; -Fx/2 &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295;0 &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; Fx/2 > > > The conjugate mirror of a sinusoidal at (u,v) will > > occur at (-u.-v) in a 2D spectrum centered at (0,0). > > > > Also, I don't understand what to do with that. If I'm not mistaken, to > > > obtain the magnitude of my IR in the frequency domain, I just need to > > > do a sqrt(real + imaginary ) on each of my complex bins, correct? So > > > then I'm left with.. a magnitude that is non-symmetrical around DC on > > > only one axis? What do I do with that, I mean, why should it have > > > negative frequency components that are not the same as the positive > > > ones, and most importantly, why on only one (arbitrarily chose) axis?? > > > These things are explained in some detail in the Gonzales and Woods > > book. > > The trick is to use circular symmetry around origo. > > Ha. thanks for the explanation, that clears up a few things about the > theoretical background. But obviously I'm not going to order a $120 > book to solve this particular problem (although I'll surely consider > the book again if it appears to cover the image processing problems I > have yet to cover) so would you mind elaborating a bit on that please?
I can't. It's one of those things that has become second nature to me; it's so obvious in my mind (I know, it's not to others, but I spent 6 years working with 2D spectra every day), so I just don't know what to explain.
> Here's what I figure should be the way to do it, based on your > explanations of 2D FFTs : I should do my hard limiting on all > magnitude bins alike regardless of which frequencies they represent, > then divide the image's FFT's complex bins with the aforementioned > magnitude. Considered both arrays would be the same size to begin with > (since they would appropriately padded) their spatial frequencies will > map so I don't even really have to wonder what I'm looking at and > simply process bins all the same. Does it sound correct?
The algorithm seems more or less correct. It's just that the $120 book shows a compariion between that approach and a few others where this gives the by far worst result. If possoble, pop into your local university bookstore to browse a copy. Depending on what you want to do with the processing, it might just save you a lot of time and grief. Rune
On 2008-08-03 14:29:27 -0300, Michel Rouzic <Michel0528@yahoo.fr> said:

> > > Rune Allnor wrote: >> If we start with the computations, then you are right in that >> the concept is to first perform 1D FFTs in one direction and >> then 1D FFTs in the other direction (although it may be possible >> to implement the 2d FFT in other ways; I don't know). >> >> The end result is, however, the same regardless of the order the >> FFTs are computed. > > Ha, I don't understand that part. How can it be the same thing > considered that, if we start from a real-only 2D image, and that > therefore the result should be symmetrical for the first processed > axis?
No! Only along y=0. A special case only of the inversion symmetry about the origin. It is also symmetric along y for x=0. Again another special case. To lower your confusion you should drop the "optimization" of using the real values only. It has confused you no end. Work it through without the "optimization" and only after it is "working" try to follow the effects of real values. I use "working" as deconvolution is a noise enhancing process. You will almost surely end up needing some form of regualrization.
>> As for the structure of the spectrum, it will become a little bit >> different than you expect. Assume, for the sake of argument, that >> you compute the FFTs in the x direction first: >> >> f(x,y) --[FFT_x]--> f(u,y) >> >> Here f(x,y) is the real-valued image, FFT_x is the FFT along >> the x dimension, and u is the 'spatial x frequency'. >> >> Since the image f(x,y) is real-valued, each line f(u) in the >> image is complex-valued and conjugate symmetric. No surprises. >> >> Now, once we compute the FFT in the other dimension, >> >> f(u,y) --[FFT_y]--> f(u,v) >> >> the input to the y-direction FFT, the columns of f(u,y), >> are complex-valued, not real-valued. This means that the >> columns f(v) in the end spectrum are no longer symmetric. >> Remember, Nyquist's theorem >> >> f < Fs/2 >> >> only applies to real-valued signals because of the >> conjugate symmetry of Euler's equations, >> >> cos(x) = 1/2 (exp(jx)+exp(-jx)) >> sin(x) = 1/j2(exp(jx)-exp(-jx)) >> >> (are the signs correct here? i don't remember...) >> >> Since complex-valued signals don't have the conjugate >> symmetric spectrum, the complex valued form of Nyquist's >> theorem becomes >> >> f < Fs. >> >> The net effect is that the 2D spectrum is cylindrical >> conjugate symmetric around origo (view with fixed-width >> font): >> ^ v >> | >> +----------------+----------------+ Fy/2 >> | | | >> | | | >> | | x (u,v) | >> 0 -+----------------+----------------+--> >> | x | | u >> | (-u,-v) | | >> | | | >> +----------------+----------------+ -Fy/2 >> -Fx/2 0 Fx/2 >> >> The conjugate mirror of a sinusoidal at (u,v) will >> occur at (-u.-v) in a 2D spectrum centered at (0,0). >> >>> Also, I don't understand what to do with that. If I'm not mistaken, to >>> obtain the magnitude of my IR in the frequency domain, I just need to >>> do a sqrt(real&#65533; + imaginary&#65533;) on each of my complex bin > s, correct? So >>> then I'm left with.. a magnitude that is non-symmetrical around DC on >>> only one axis? What do I do with that, I mean, why should it have >>> negative frequency components that are not the same as the positive >>> ones, and most importantly, why on only one (arbitrarily chose) axis?? >> >> These things are explained in some detail in the Gonzales and Woods >> book. >> The trick is to use circular symmetry around origo. > > Ha. thanks for the explanation, that clears up a few things about the > theoretical background. But obviously I'm not going to order a $120 > book to solve this particular problem (although I'll surely consider > the book again if it appears to cover the image processing problems I > have yet to cover) so would you mind elaborating a bit on that please? > > Here's what I figure should be the way to do it, based on your > explanations of 2D FFTs : I should do my hard limiting on all > magnitude bins alike regardless of which frequencies they represent, > then divide the image's FFT's complex bins with the aforementioned > magnitude. Considered both arrays would be the same size to begin with > (since they would appropriately padded) their spatial frequencies will > map so I don't even really have to wonder what I'm looking at and > simply process bins all the same. Does it sound correct?
You might look up ImageJ, which is public domain image processing
software and has facilities to directly manipulate 2D FFT data. It
also, by the way, has several implementations of deconvolution (as
plug-ins).

Chris
==================
Chris Bore
BORES Signal Processing
www.bores.com

On Aug 3, 2:17&#4294967295;pm, Michel Rouzic <Michel0...@yahoo.fr> wrote:
> Here's my problem. I have an image to be processed (well, deconvolved > pretty much), and a star-shaped impulse response I want to deconvolve > it with. The impulse response is symmetric in both the vertical and > horizontal axis and pretty much centred, but not very precisely. Of > course both the image and the impulse response are appropriately zero- > padded as to avoid any circular convolution. > > My idea is to transform the IR into a frequency-domain (magnitude > only) representation using something like a 2D FFT, hard limit the > frequency-domain magnitude bins under a threshold value to make sure > that the weaker frequency components don't get amplified more than > they should, and then perform the deconvolution. > > Well at this point I guess I could simply do a f(x) = 1/x on all the > magnitude bins, zero the phase (as to make the IR centered around 0, > 0), IFFT and convolve in the space-domain, or I could FFT the image to > be deconvolved and multiply (or rather divide) with the frequency- > domain magnitude of the hard-limited impulse response. The problem is > that I don't understand how to get the magnitude of the impulse > response FFT (or similar). > > There's something I don't understand about 2D DFTs. Just like the DFT > of a 1D real signal is represented by an array of complex bins > representing from DC to the Nyquist frequency (omitting the redundant > negative frequencies) I would expect a 2D DFT of a (real) image to be > a complex array ranging from each axis' DC to each axis' Nyquist > frequency (although I realise this is surely not possible since it > would take twice less space to represent a real image under that form) > > Instead, it seems that 2D DFTs are achieved by performing 1D DFTs on > each row, and then performing 1D DFTs on each column. If I understand > correctly, the result is a complex array representing from -Nyquist to > +Nyquist on each axis, and is symmetrical on one axis due to the fact > that the first series of 1D DFTs is performed on real data, as the > rest is on complex data and therefore has non-redundant negative > frequency components. My problem with that is that I don't understand > why we should "privilege" an axis by processing it first, and > therefore obtain a different result that if we have processed the > other axis first. > > Also, I don't understand what to do with that. If I'm not mistaken, to > obtain the magnitude of my IR in the frequency domain, I just need to > do a sqrt(real&#4294967295; + imaginary&#4294967295;) on each of my complex bins, correct? So > then I'm left with.. a magnitude that is non-symmetrical around DC on > only one axis? What do I do with that, I mean, why should it have > negative frequency components that are not the same as the positive > ones, and most importantly, why on only one (arbitrarily chose) axis?? > > I should also probably mention that I'm implementing that using FFTW 3 > (no Matlab to toy around with), and that its documentation is pretty > poor on the very topic of my problem. It seemingly doesn't tell what > the output matches to depending on the transform. For example I > noticed that the output is quite different between a 2D Discrete > Hartley Transform and a 2D DFT (and also a 2D DCT, but I know why that > one is different, I at least understand that, but the documentation > won't say anything about it), why? Disclaimer : I don't really know > what the DHT does, but I find it striking that its output (as a real > image, so I take it it's in their "half complex" format that goes DC > Re1 Re2 Re3 &#4294967295;Nyquist Im3 Im2 Im1) is made of only two parts (which I > guess would come as one if it were to be put as a complex array as it > may (or may not? the documentation doesn't seem to discuss this) be in > "half complex" format). > > Thanks in advance
Michel Rouzic wrote:
(snip)

> Instead, it seems that 2D DFTs are achieved by performing 1D DFTs on > each row, and then performing 1D DFTs on each column. If I understand > correctly, the result is a complex array representing from -Nyquist to > +Nyquist on each axis, and is symmetrical on one axis due to the fact > that the first series of 1D DFTs is performed on real data, as the > rest is on complex data and therefore has non-redundant negative > frequency components. My problem with that is that I don't understand > why we should "privilege" an axis by processing it first, and > therefore obtain a different result that if we have processed the > other axis first.
I am pretty sure that it doesn't depend on the order of the transforms. The Fourier transform is a linear operator. One way to look at it is through the normal modes of a square drum head. (Well, that would by 2D DST, but the symmetry is the same.) In the 2D transform array the upper left corner is the DC component, or, for a drum head, the mode with no nodal lines. As you go across the matrix in one direction you increase the number of vertical node lines. In the other direction horizontal node lines. Do the transform by hand for a 2x2 or 3x3 matrix. You will then see the symmetry. -- glen
On Aug 3, 9:17 am, Michel Rouzic <Michel0...@yahoo.fr> wrote:
> Here's my problem. I have animageto be processed (well, deconvolved > pretty much), and a star-shaped impulse response I want to deconvolve > it with. The impulse response is symmetric in both the vertical and > horizontal axis and pretty much centred, but not very precisely. Of > course both theimageand the impulse response are appropriately zero- > padded as to avoid any circular convolution. >
On the PSF topic, you (and others) may want to look at Tria at http://www.quarktet.com/Tria.html. It has built-in centering and cleaning functions for PSF. If your PSF is not centered before deconvolution, your processed image will shift by the difference between the center of the image and the center of the PSF. This may or may not be important. Tria also uses a non-iterative pseudo- inverse filter (PIF) (sometimes called the Wiener filter) that,a s you can see from the image gallery does not produce edge effects when applied correctly. The ImageJ plug-ins, as far as I can tell, are iterative processes that some people feel are more effective than the PIF. However, there are a few papers that demonstrate that this is not the case. Tria is free to try, and costs more than an image processing book, but much less than MatLab. Best Regards, Jim C