FFT implementation in ImageJ

[Kenneth R Sloan]

Read a bit more!  If RI is a real-valued image, then FT(RI) is symmetric - which means you
only have to store half of the (complex) result.  So…it’s a wash (but requires knowledge about
the “scrambling” of your representation.)

You can think of this from an information theory point of view: the FT does not create any new information, so it should be obvious that you don’t need more bits to represent the result.

Perhaps trickier is the other way around - given a symmetric spectrum,
the inverseFT gives a real image.  Depending on the datatype you are using, this
might require special handling (coercing tiny imaginary values to 0).

This is one of the reasons why you might prefer to use an “industrial strength” FFT implementation if you want to USE the FT (instead of wanting to implement from first principles so that you understand it).

Of course, if memory is free (or your images are small), you might not find this space-saving complication as a worthwhile tradeoff.  When the FFT was invented, memory was NOT free!
Today, it might be more correct to implement FT and inverseFT as strictly complex<->complex, and
convert images real<->complex.

On the third hand, for large enough images, you should also consider memory access patterns.  Memory might APPEAR to be “free” - but USING all that memory may not be.

Are we having fun, yet?

--
Kenneth Sloan
[hidden email]<mailto:[hidden email]>
"La lutte elle-même vers les sommets suffit à remplir un coeur d'homme; il faut imaginer Sisyphe heureux."


On Jun 19, 2014, at 11:06 , Student1 <[hidden email]<mailto:[hidden email]>> wrote:

Hi Burger,

I kind of wanted to implement the FFT by hand instead of using a library because I really want to understand how the FFT works and make sure I got the implementation correct. I know that combining the complex and real array is doubling the size of the array and I don’t know exactly how to fix that. I don’t thinks theres a problems in the rest of the code I implemented.

Thanks.

Mariam Dost


On Jun 19, 2014, at 9:46 AM, Burger Wilhelm <[hidden email]<mailto:[hidden email]>> wrote:

Hello Miriam,

have you considered the FFT implementation in Apache Commons Math?:
http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/transform/FastFourierTransformer.html

Works "out of the box" and results are correct.

--Wilhelm Burger

________________________________________
From: ImageJ Interest Group [[hidden email]] On Behalf Of Student1 [[hidden email]]
Sent: Wednesday, June 18, 2014 23:58
To: [hidden email]
Subject: Re: FFT implementation in ImageJ

Hi,

First of all thanks for everybody’s input. I implemented an FFT algorithm by basically looking it up online and I am tracing through the calculation by hand to get a better understanding of the FFT. The FFT I implemented does not produce the image that should be produced by ImageJ. I was wondering if anyone can see where I am going wrong with my calculation. My work is here:

http://stackoverflow.com/questions/24291896/fft-image-is-not-the-same-as-the-fft-image-produced-in-imagej

Thanks.

Mariam Dost

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html


--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Daniel Sage]

FFT is the important block of many image processing algorithms. So, it is crucial to get an efficient FFT. The problem here is that the efficiency is related to the available resources of the client machine: memory, number of cores, GPU, OS, C-language.

In the Java world, we can find many FFT libraries:
- Several translation of the Numerical Recipes FFT
- FHT of ImageJ
- Colt (Parallel Colt) and JTransform
- Apache Commons Math
- Java Wrapper for the C library FFTW (Fast Fourier Transform in the West)
- Java wrapper for FFTPACK
- JCUFFT (Cuda / GPU)


The choice of the "good" FFT library is difficult because several criteria have to taken into account:

- execution time (without any specific hardware, Matlab has the fastest one!)
- multithreading
- efficient for every size of signals, not only the the power of 2
- exactness for every size (padding?)
- energy conservation
- data type: float or double and data dimension: 1D, 2D, and 3D
- data structure: a complex object for every pixel, real and imaginary arrays, or interleave format. Surprisingly, the last structure is use in many libraries.
- compactness: using the Hermitian symmetry to save memory
- providing binaries (DLL) for several OS when the library is written in C/C++
- licensing
- checking the presence of GPU to run on it
- and probably the most important for the user: the availability of complex-number methods to process the coefficients of the Fourier transform avoiding multiple copies from one data structure to a other data structure.

Piotr Wendykier has some benchmarking: http://dl.acm.org/citation.cfm?id=1824809

Best,

Daniel

> Read a bit more!  If RI is a real-valued image, then FT(RI) is symmetric - which means you
> only have to store half of the (complex) result.  So…it’s a wash (but requires knowledge about
> the “scrambling” of your representation.)
>
> You can think of this from an information theory point of view: the FT does not create any new information, so it should be obvious that you don’t need more bits to represent the result.
>
> Perhaps trickier is the other way around - given a symmetric spectrum,
> the inverseFT gives a real image.  Depending on the datatype you are using, this
> might require special handling (coercing tiny imaginary values to 0).
>
> This is one of the reasons why you might prefer to use an “industrial strength” FFT implementation if you want to USE the FT (instead of wanting to implement from first principles so that you understand it).
>
> Of course, if memory is free (or your images are small), you might not find this space-saving complication as a worthwhile tradeoff.  When the FFT was invented, memory was NOT free!
> Today, it might be more correct to implement FT and inverseFT as strictly complex<->complex, and
> convert images real<->complex.
>
> On the third hand, for large enough images, you should also consider memory access patterns.  Memory might APPEAR to be “free” - but USING all that memory may not be.
>
> Are we having fun, yet?
>
> --
> Kenneth Sloan
> [hidden email]<mailto:[hidden email]>
> "La lutte elle-même vers les sommets suffit à remplir un coeur d'homme; il faut imaginer Sisyphe heureux."
>
>
> On Jun 19, 2014, at 11:06 , Student1 <[hidden email]<mailto:[hidden email]>> wrote:
>
> Hi Burger,
>
> I kind of wanted to implement the FFT by hand instead of using a library because I really want to understand how the FFT works and make sure I got the implementation correct. I know that combining the complex and real array is doubling the size of the array and I don’t know exactly how to fix that. I don’t thinks theres a problems in the rest of the code I implemented.
>
> Thanks.
>
> Mariam Dost
>
>
> On Jun 19, 2014, at 9:46 AM, Burger Wilhelm <[hidden email]<mailto:[hidden email]>> wrote:
>
> Hello Miriam,
>
> have you considered the FFT implementation in Apache Commons Math?:
> http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/transform/FastFourierTransformer.html
>
> Works "out of the box" and results are correct.
>
> --Wilhelm Burger
>
> ________________________________________
> From: ImageJ Interest Group [[hidden email]] On Behalf Of Student1 [[hidden email]]
> Sent: Wednesday, June 18, 2014 23:58
> To: [hidden email]
> Subject: Re: FFT implementation in ImageJ
>
> Hi,
>
> First of all thanks for everybody’s input. I implemented an FFT algorithm by basically looking it up online and I am tracing through the calculation by hand to get a better understanding of the FFT. The FFT I implemented does not produce the image that should be produced by ImageJ. I was wondering if anyone can see where I am going wrong with my calculation. My work is here:
>
> http://stackoverflow.com/questions/24291896/fft-image-is-not-the-same-as-the-fft-image-produced-in-imagej
>
> Thanks.
>
> Mariam Dost
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Kenneth R Sloan]

At some point, it helps to use precise terminology.

FT  - Fourier Transform - a mathematical function mapping complex signals to complex spectra
FFT - Fast Fourier Transform - a particular algorithm which implements the FT

Most implementations of the FT depend, more or less, on the FFT.  Different choices for representation,
etc.  These might still be called “implementations of the FFT”.

Others may differ enough that they deserve their own name - when they do, it’s proper to say
that they implement the FT and not the FFT!

Still others are implementation of something “similar, but not quite the same as the FT”.  They should probably NOT ever be called either FT or FFT.

I am reminded of another system I use regularly which has a primitive called “cube” - which takes
three dimensions (height, width, depth).  It always hurts my brain to type it.

I am also reminded of a question I once actually heard asked at an imaging conference: “has anyone
considered implementing the Fast Fourier Transform using optics?”

--
Kenneth Sloan
[hidden email]<mailto:[hidden email]>
"La lutte elle-même vers les sommets suffit à remplir un coeur d'homme; il faut imaginer Sisyphe heureux."


On Jun 20, 2014, at 01:46 , Daniel Sage <[hidden email]<mailto:[hidden email]>> wrote:

FFT is the important block of many image processing algorithms. So, it is crucial to get an efficient FFT. The problem here is that the efficiency is related to the available resources of the client machine: memory, number of cores, GPU, OS, C-language.

In the Java world, we can find many FFT libraries:
- Several translation of the Numerical Recipes FFT
- FHT of ImageJ
- Colt (Parallel Colt) and JTransform
- Apache Commons Math
- Java Wrapper for the C library FFTW (Fast Fourier Transform in the West)
- Java wrapper for FFTPACK
- JCUFFT (Cuda / GPU)


The choice of the "good" FFT library is difficult because several criteria have to taken into account:

- execution time (without any specific hardware, Matlab has the fastest one!)
- multithreading
- efficient for every size of signals, not only the the power of 2
- exactness for every size (padding?)
- energy conservation
- data type: float or double and data dimension: 1D, 2D, and 3D
- data structure: a complex object for every pixel, real and imaginary arrays, or interleave format. Surprisingly, the last structure is use in many libraries.
- compactness: using the Hermitian symmetry to save memory
- providing binaries (DLL) for several OS when the library is written in C/C++
- licensing
- checking the presence of GPU to run on it
- and probably the most important for the user: the availability of complex-number methods to process the coefficients of the Fourier transform avoiding multiple copies from one data structure to a other data structure.

Piotr Wendykier has some benchmarking: http://dl.acm.org/citation.cfm?id=1824809

Best,

Daniel

Read a bit more!  If RI is a real-valued image, then FT(RI) is symmetric - which means you
only have to store half of the (complex) result.  So…it’s a wash (but requires knowledge about
the “scrambling” of your representation.)

You can think of this from an information theory point of view: the FT does not create any new information, so it should be obvious that you don’t need more bits to represent the result.

Perhaps trickier is the other way around - given a symmetric spectrum,
the inverseFT gives a real image.  Depending on the datatype you are using, this
might require special handling (coercing tiny imaginary values to 0).

This is one of the reasons why you might prefer to use an “industrial strength” FFT implementation if you want to USE the FT (instead of wanting to implement from first principles so that you understand it).

Of course, if memory is free (or your images are small), you might not find this space-saving complication as a worthwhile tradeoff.  When the FFT was invented, memory was NOT free!
Today, it might be more correct to implement FT and inverseFT as strictly complex<->complex, and
convert images real<->complex.

On the third hand, for large enough images, you should also consider memory access patterns.  Memory might APPEAR to be “free” - but USING all that memory may not be.

Are we having fun, yet?

--
Kenneth Sloan
[hidden email]<mailto:[hidden email]><mailto:[hidden email]>
"La lutte elle-même vers les sommets suffit à remplir un coeur d'homme; il faut imaginer Sisyphe heureux."


On Jun 19, 2014, at 11:06 , Student1 <[hidden email]<mailto:[hidden email]><mailto:[hidden email]>> wrote:

Hi Burger,

I kind of wanted to implement the FFT by hand instead of using a library because I really want to understand how the FFT works and make sure I got the implementation correct. I know that combining the complex and real array is doubling the size of the array and I don’t know exactly how to fix that. I don’t thinks theres a problems in the rest of the code I implemented.

Thanks.

Mariam Dost


On Jun 19, 2014, at 9:46 AM, Burger Wilhelm <[hidden email]<mailto:[hidden email]><mailto:[hidden email]>> wrote:

Hello Miriam,

have you considered the FFT implementation in Apache Commons Math?:
http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/transform/FastFourierTransformer.html

Works "out of the box" and results are correct.

--Wilhelm Burger

________________________________________
From: ImageJ Interest Group [[hidden email]] On Behalf Of Student1 [[hidden email]]
Sent: Wednesday, June 18, 2014 23:58
To: [hidden email]
Subject: Re: FFT implementation in ImageJ

Hi,

First of all thanks for everybody’s input. I implemented an FFT algorithm by basically looking it up online and I am tracing through the calculation by hand to get a better understanding of the FFT. The FFT I implemented does not produce the image that should be produced by ImageJ. I was wondering if anyone can see where I am going wrong with my calculation. My work is here:

http://stackoverflow.com/questions/24291896/fft-image-is-not-the-same-as-the-fft-image-produced-in-imagej

Thanks.

Mariam Dost

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html


--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html


--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[lechristophe]

An alternative meaning to the FFT acronym may be emerging though

https://twitter.com/mc_hankins/status/479721275842371584


2014-06-20 11:24 GMT+02:00 Kenneth R Sloan <[hidden email]>:

> At some point, it helps to use precise terminology.
>
> FT  - Fourier Transform - a mathematical function mapping complex signals
> to complex spectra
> FFT - Fast Fourier Transform - a particular algorithm which implements the
> FT
>
> Most implementations of the FT depend, more or less, on the FFT.
>  Different choices for representation,
> etc.  These might still be called “implementations of the FFT”.
>
> Others may differ enough that they deserve their own name - when they do,
> it’s proper to say
> that they implement the FT and not the FFT!
>
> Still others are implementation of something “similar, but not quite the
> same as the FT”.  They should probably NOT ever be called either FT or FFT.
>
> I am reminded of another system I use regularly which has a primitive
> called “cube” - which takes
> three dimensions (height, width, depth).  It always hurts my brain to type
> it.
>
> I am also reminded of a question I once actually heard asked at an imaging
> conference: “has anyone
> considered implementing the Fast Fourier Transform using optics?”
>
> --
> Kenneth Sloan
> [hidden email]<mailto:[hidden email]>
> "La lutte elle-même vers les sommets suffit à remplir un coeur d'homme; il
> faut imaginer Sisyphe heureux."
>
>
> On Jun 20, 2014, at 01:46 , Daniel Sage <[hidden email]<mailto:
> [hidden email]>> wrote:
>
> FFT is the important block of many image processing algorithms. So, it is
> crucial to get an efficient FFT. The problem here is that the efficiency is
> related to the available resources of the client machine: memory, number of
> cores, GPU, OS, C-language.
>
> In the Java world, we can find many FFT libraries:
> - Several translation of the Numerical Recipes FFT
> - FHT of ImageJ
> - Colt (Parallel Colt) and JTransform
> - Apache Commons Math
> - Java Wrapper for the C library FFTW (Fast Fourier Transform in the West)
> - Java wrapper for FFTPACK
> - JCUFFT (Cuda / GPU)
>
>
> The choice of the "good" FFT library is difficult because several criteria
> have to taken into account:
>
> - execution time (without any specific hardware, Matlab has the fastest
> one!)
> - multithreading
> - efficient for every size of signals, not only the the power of 2
> - exactness for every size (padding?)
> - energy conservation
> - data type: float or double and data dimension: 1D, 2D, and 3D
> - data structure: a complex object for every pixel, real and imaginary
> arrays, or interleave format. Surprisingly, the last structure is use in
> many libraries.
> - compactness: using the Hermitian symmetry to save memory
> - providing binaries (DLL) for several OS when the library is written in
> C/C++
> - licensing
> - checking the presence of GPU to run on it
> - and probably the most important for the user: the availability of
> complex-number methods to process the coefficients of the Fourier transform
> avoiding multiple copies from one data structure to a other data structure.
>
> Piotr Wendykier has some benchmarking:
> http://dl.acm.org/citation.cfm?id=1824809
>
> Best,
>
> Daniel
>
> Read a bit more!  If RI is a real-valued image, then FT(RI) is symmetric -
> which means you
> only have to store half of the (complex) result.  So…it’s a wash (but
> requires knowledge about
> the “scrambling” of your representation.)
>
> You can think of this from an information theory point of view: the FT
> does not create any new information, so it should be obvious that you don’t
> need more bits to represent the result.
>
> Perhaps trickier is the other way around - given a symmetric spectrum,
> the inverseFT gives a real image.  Depending on the datatype you are
> using, this
> might require special handling (coercing tiny imaginary values to 0).
>
> This is one of the reasons why you might prefer to use an “industrial
> strength” FFT implementation if you want to USE the FT (instead of wanting
> to implement from first principles so that you understand it).
>
> Of course, if memory is free (or your images are small), you might not
> find this space-saving complication as a worthwhile tradeoff.  When the FFT
> was invented, memory was NOT free!
> Today, it might be more correct to implement FT and inverseFT as strictly
> complex<->complex, and
> convert images real<->complex.
>
> On the third hand, for large enough images, you should also consider
> memory access patterns.  Memory might APPEAR to be “free” - but USING all
> that memory may not be.
>
> Are we having fun, yet?
>
> --
> Kenneth Sloan
> [hidden email]<mailto:[hidden email]><mailto:[hidden email]>
> "La lutte elle-même vers les sommets suffit à remplir un coeur d'homme; il
> faut imaginer Sisyphe heureux."
>
>
> On Jun 19, 2014, at 11:06 , Student1 <[hidden email]<mailto:
> [hidden email]><mailto:[hidden email]>> wrote:
>
> Hi Burger,
>
> I kind of wanted to implement the FFT by hand instead of using a library
> because I really want to understand how the FFT works and make sure I got
> the implementation correct. I know that combining the complex and real
> array is doubling the size of the array and I don’t know exactly how to fix
> that. I don’t thinks theres a problems in the rest of the code I
> implemented.
>
> Thanks.
>
> Mariam Dost
>
>
> On Jun 19, 2014, at 9:46 AM, Burger Wilhelm <
> [hidden email]<mailto:[hidden email]
> ><mailto:[hidden email]>> wrote:
>
> Hello Miriam,
>
> have you considered the FFT implementation in Apache Commons Math?:
>
> http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/transform/FastFourierTransformer.html
>
> Works "out of the box" and results are correct.
>
> --Wilhelm Burger
>
> ________________________________________
> From: ImageJ Interest Group [[hidden email]] On Behalf Of Student1 [
> [hidden email]]
> Sent: Wednesday, June 18, 2014 23:58
> To: [hidden email]
> Subject: Re: FFT implementation in ImageJ
>
> Hi,
>
> First of all thanks for everybody’s input. I implemented an FFT algorithm
> by basically looking it up online and I am tracing through the calculation
> by hand to get a better understanding of the FFT. The FFT I implemented
> does not produce the image that should be produced by ImageJ. I was
> wondering if anyone can see where I am going wrong with my calculation. My
> work is here:
>
>
> http://stackoverflow.com/questions/24291896/fft-image-is-not-the-same-as-the-fft-image-produced-in-imagej
>
> Thanks.
>
> Mariam Dost
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Gabriel Landini]

On Friday 20 Jun 2014 09:24:37 you wrote:
> I am also reminded of a question I once actually heard asked at an imaging
> conference: “has anyone considered implementing the Fast Fourier Transform
> using optics?”

Doesn't the cover of Pink Floyd's Dark side of the Moon show one of such
devices? :-)

Cheers
Gabriel

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Herbie-4]

Dear Kenneth,

I widely agree with your comments, but...

 > Most implementations of the FT depend, more or less, on the FFT.

One could state even more precisely "Most digital implementations of the
DFT...", where DFT is the discrete FT that assumes what is called a
finite signal, something that implies compromises: Signals that have
limited support (most often even equal support) in _both_ domains.
According to the FT this is impossible, i.e. such signal representations
imply errors per se.

I fully agree that, e.g. DFT algorithms that are not limited to a power
of two support should not be called FFT. Their results, while being
certainly useful, often cannot be compared with those obtained from
("true") FFT implementations. That said, and with regard to ImageJ-n, it
is desirable to have access to a ("true") FFT implementation as well.

BTW, I can not realize any significant discrepancies between results
obtained with FFTJ [using float]--a "true" FFT implementation--and with
the ImageJ-1 FHT-approach.

 > “has anyone considered implementing the Fast
 > Fourier Transform using optics?”

I know what you mean (FT by diffraction of coherent light), but there is
also a field of digital optical processing...

Best regards

Herbie

:::::::::::::::::::::::::::::::::::::::::
On 20.06.14 11:24, Kenneth R Sloan wrote:

> At some point, it helps to use precise terminology.
>
> FT  - Fourier Transform - a mathematical function mapping complex
> signals to complex spectra FFT - Fast Fourier Transform - a
> particular algorithm which implements the FT
>
> Most implementations of the FT depend, more or less, on the FFT.
> Different choices for representation, etc.  These might still be
> called “implementations of the FFT”.
>
> Others may differ enough that they deserve their own name - when they
> do, it’s proper to say that they implement the FT and not the FFT!
>
> Still others are implementation of something “similar, but not quite
> the same as the FT”.  They should probably NOT ever be called either
> FT or FFT.
>
> I am reminded of another system I use regularly which has a primitive
> called “cube” - which takes three dimensions (height, width, depth).
> It always hurts my brain to type it.
>
> I am also reminded of a question I once actually heard asked at an
> imaging conference: “has anyone considered implementing the Fast
> Fourier Transform using optics?”
>
> -- Kenneth Sloan [hidden email]<mailto:[hidden email]> "La lutte
> elle-même vers les sommets suffit à remplir un coeur d'homme; il faut
> imaginer Sisyphe heureux."
>
>
> On Jun 20, 2014, at 01:46 , Daniel Sage
> <[hidden email]<mailto:[hidden email]>> wrote:
>
> FFT is the important block of many image processing algorithms. So,
> it is crucial to get an efficient FFT. The problem here is that the
> efficiency is related to the available resources of the client
> machine: memory, number of cores, GPU, OS, C-language.
>
> In the Java world, we can find many FFT libraries: - Several
> translation of the Numerical Recipes FFT - FHT of ImageJ - Colt
> (Parallel Colt) and JTransform - Apache Commons Math - Java Wrapper
> for the C library FFTW (Fast Fourier Transform in the West) - Java
> wrapper for FFTPACK - JCUFFT (Cuda / GPU)
>
>
> The choice of the "good" FFT library is difficult because several
> criteria have to taken into account:
>
> - execution time (without any specific hardware, Matlab has the
> fastest one!) - multithreading - efficient for every size of signals,
> not only the the power of 2 - exactness for every size (padding?) -
> energy conservation - data type: float or double and data dimension:
> 1D, 2D, and 3D - data structure: a complex object for every pixel,
> real and imaginary arrays, or interleave format. Surprisingly, the
> last structure is use in many libraries. - compactness: using the
> Hermitian symmetry to save memory - providing binaries (DLL) for
> several OS when the library is written in C/C++ - licensing -
> checking the presence of GPU to run on it - and probably the most
> important for the user: the availability of complex-number methods to
> process the coefficients of the Fourier transform avoiding multiple
> copies from one data structure to a other data structure.
>
> Piotr Wendykier has some benchmarking:
> http://dl.acm.org/citation.cfm?id=1824809
>
> Best,
>
> Daniel
>
> Read a bit more!  If RI is a real-valued image, then FT(RI) is
> symmetric - which means you only have to store half of the (complex)
> result.  So…it’s a wash (but requires knowledge about the
> “scrambling” of your representation.)
>
> You can think of this from an information theory point of view: the
> FT does not create any new information, so it should be obvious that
> you don’t need more bits to represent the result.
>
> Perhaps trickier is the other way around - given a symmetric
> spectrum, the inverseFT gives a real image.  Depending on the
> datatype you are using, this might require special handling (coercing
> tiny imaginary values to 0).
>
> This is one of the reasons why you might prefer to use an “industrial
> strength” FFT implementation if you want to USE the FT (instead of
> wanting to implement from first principles so that you understand
> it).
>
> Of course, if memory is free (or your images are small), you might
> not find this space-saving complication as a worthwhile tradeoff.
> When the FFT was invented, memory was NOT free! Today, it might be
> more correct to implement FT and inverseFT as strictly
> complex<->complex, and convert images real<->complex.
>
> On the third hand, for large enough images, you should also consider
> memory access patterns.  Memory might APPEAR to be “free” - but USING
> all that memory may not be.
>
> Are we having fun, yet?
>
> -- Kenneth Sloan
> [hidden email]<mailto:[hidden email]><mailto:[hidden email]> "La lutte
> elle-même vers les sommets suffit à remplir un coeur d'homme; il faut
> imaginer Sisyphe heureux."
>
>
> On Jun 19, 2014, at 11:06 , Student1
> <[hidden email]<mailto:[hidden email]><mailto:[hidden email]>>
> wrote:
>
> Hi Burger,
>
> I kind of wanted to implement the FFT by hand instead of using a
> library because I really want to understand how the FFT works and
> make sure I got the implementation correct. I know that combining the
> complex and real array is doubling the size of the array and I don’t
> know exactly how to fix that. I don’t thinks theres a problems in the
> rest of the code I implemented.
>
> Thanks.
>
> Mariam Dost
>
>
> On Jun 19, 2014, at 9:46 AM, Burger Wilhelm
> <[hidden email]<mailto:[hidden email]><mailto:[hidden email]>>
> wrote:
>
> Hello Miriam,
>
> have you considered the FFT implementation in Apache Commons Math?:
> http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/transform/FastFourierTransformer.html
>
>  Works "out of the box" and results are correct.
>
> --Wilhelm Burger
>
> ________________________________________ From: ImageJ Interest Group
> [[hidden email]] On Behalf Of Student1
> [[hidden email]] Sent: Wednesday, June 18, 2014 23:58 To:
> [hidden email] Subject: Re: FFT implementation in ImageJ
>
> Hi,
>
> First of all thanks for everybody’s input. I implemented an FFT
> algorithm by basically looking it up online and I am tracing through
> the calculation by hand to get a better understanding of the FFT. The
> FFT I implemented does not produce the image that should be produced
> by ImageJ. I was wondering if anyone can see where I am going wrong
> with my calculation. My work is here:
>
> http://stackoverflow.com/questions/24291896/fft-image-is-not-the-same-as-the-fft-image-produced-in-imagej
>
>  Thanks.
>
> Mariam Dost
>
> -- ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>
>
> -- ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>
> -- ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>
>
> -- ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Dimiter Prodanov (imec)]

Dear Mariam,

It is a good idea to always check against Matlab when you do processing in the Fourier domain.
The implementations there are tested and they are numerically correct.
If you want to understand FFT you have to start with the papers of  Cooley  Tukey and Lanczos.

As fat as ImageJ is concerned there is no standard for handling complex float or double numbers which makes using FFT-based algorithms more complicated.
Daniel pointed out all the issues that have to be taken on account when deciding which algorithm to use.
I would say that the non-pow2 issue is less important as confocal microscopic images are power of 2.
What is would add to the list are also windowing functions and boundary conditions.
The usual zero-padding for example in most of the cases is not a good idea for images.  
Dimensionality is also not that important as FFT is separable.

 
Best regards,


Dimiter Prodanov

 

 

 

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Kenneth R Sloan]

Non power of 2 should not really be an issue (I was using an FFT using arbitrary factoring in 1978), but when it is, remember that “zero-padding” is a kind of mis-translation from the original signal processing terminology.  What you really want (if you are going to pad) is “DC-padding”.
Other methods are also useful - the idea is to make the finite signal cyclic without introducing spurious energy at the edges.

--
Kenneth Sloan
[hidden email]<mailto:[hidden email]>
"La lutte elle-même vers les sommets suffit à remplir un coeur d'homme; il faut imaginer Sisyphe heureux."


On Jun 23, 2014, at 07:24 , Dimiter Prodanov (imec) <[hidden email]<mailto:[hidden email]>> wrote:

Dear Mariam,

It is a good idea to always check against Matlab when you do processing in the Fourier domain.
The implementations there are tested and they are numerically correct.
If you want to understand FFT you have to start with the papers of  Cooley  Tukey and Lanczos.

As fat as ImageJ is concerned there is no standard for handling complex float or double numbers which makes using FFT-based algorithms more complicated.
Daniel pointed out all the issues that have to be taken on account when deciding which algorithm to use.
I would say that the non-pow2 issue is less important as confocal microscopic images are power of 2.
What is would add to the list are also windowing functions and boundary conditions.
The usual zero-padding for example in most of the cases is not a good idea for images.
Dimensionality is also not that important as FFT is separable.


Best regards,


Dimiter Prodanov







--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html


--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Piotr Wendykier-2]

I have just released JTransforms 3.0:
http://sites.google.com/site/piotrwendykier/software/jtransforms
The main new feature is a support for 1D arrays of size up to 2^63 via
JLargeArrays library (http://github.com/IcmVis/JLargeArrays/).
Benchmarks are also updated.

Piotr Wendykier

On Fri, Jun 20, 2014 at 8:46 AM, Daniel Sage <[hidden email]> wrote:

> FFT is the important block of many image processing algorithms. So, it is crucial to get an efficient FFT. The problem here is that the efficiency is related to the available resources of the client machine: memory, number of cores, GPU, OS, C-language.
>
> In the Java world, we can find many FFT libraries:
> - Several translation of the Numerical Recipes FFT
> - FHT of ImageJ
> - Colt (Parallel Colt) and JTransform
> - Apache Commons Math
> - Java Wrapper for the C library FFTW (Fast Fourier Transform in the West)
> - Java wrapper for FFTPACK
> - JCUFFT (Cuda / GPU)
>
>
> The choice of the "good" FFT library is difficult because several criteria have to taken into account:
>
> - execution time (without any specific hardware, Matlab has the fastest one!)
> - multithreading
> - efficient for every size of signals, not only the the power of 2
> - exactness for every size (padding?)
> - energy conservation
> - data type: float or double and data dimension: 1D, 2D, and 3D
> - data structure: a complex object for every pixel, real and imaginary arrays, or interleave format. Surprisingly, the last structure is use in many libraries.
> - compactness: using the Hermitian symmetry to save memory
> - providing binaries (DLL) for several OS when the library is written in C/C++
> - licensing
> - checking the presence of GPU to run on it
> - and probably the most important for the user: the availability of complex-number methods to process the coefficients of the Fourier transform avoiding multiple copies from one data structure to a other data structure.
>
> Piotr Wendykier has some benchmarking: http://dl.acm.org/citation.cfm?id=1824809
>
> Best,
>
> Daniel
>
>> Read a bit more!  If RI is a real-valued image, then FT(RI) is symmetric - which means you
>> only have to store half of the (complex) result.  So…it’s a wash (but requires knowledge about
>> the “scrambling” of your representation.)
>>
>> You can think of this from an information theory point of view: the FT does not create any new information, so it should be obvious that you don’t need more bits to represent the result.
>>
>> Perhaps trickier is the other way around - given a symmetric spectrum,
>> the inverseFT gives a real image.  Depending on the datatype you are using, this
>> might require special handling (coercing tiny imaginary values to 0).
>>
>> This is one of the reasons why you might prefer to use an “industrial strength” FFT implementation if you want to USE the FT (instead of wanting to implement from first principles so that you understand it).
>>
>> Of course, if memory is free (or your images are small), you might not find this space-saving complication as a worthwhile tradeoff.  When the FFT was invented, memory was NOT free!
>> Today, it might be more correct to implement FT and inverseFT as strictly complex<->complex, and
>> convert images real<->complex.
>>
>> On the third hand, for large enough images, you should also consider memory access patterns.  Memory might APPEAR to be “free” - but USING all that memory may not be.
>>
>> Are we having fun, yet?
>>
>> --
>> Kenneth Sloan
>> [hidden email]<mailto:[hidden email]>
>> "La lutte elle-même vers les sommets suffit à remplir un coeur d'homme; il faut imaginer Sisyphe heureux."
>>
>>
>> On Jun 19, 2014, at 11:06 , Student1 <[hidden email]<mailto:[hidden email]>> wrote:
>>
>> Hi Burger,
>>
>> I kind of wanted to implement the FFT by hand instead of using a library because I really want to understand how the FFT works and make sure I got the implementation correct. I know that combining the complex and real array is doubling the size of the array and I don’t know exactly how to fix that. I don’t thinks theres a problems in the rest of the code I implemented.
>>
>> Thanks.
>>
>> Mariam Dost
>>
>>
>> On Jun 19, 2014, at 9:46 AM, Burger Wilhelm <[hidden email]<mailto:[hidden email]>> wrote:
>>
>> Hello Miriam,
>>
>> have you considered the FFT implementation in Apache Commons Math?:
>> http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/transform/FastFourierTransformer.html
>>
>> Works "out of the box" and results are correct.
>>
>> --Wilhelm Burger
>>
>> ________________________________________
>> From: ImageJ Interest Group [[hidden email]] On Behalf Of Student1 [[hidden email]]
>> Sent: Wednesday, June 18, 2014 23:58
>> To: [hidden email]
>> Subject: Re: FFT implementation in ImageJ
>>
>> Hi,
>>
>> First of all thanks for everybody’s input. I implemented an FFT algorithm by basically looking it up online and I am tracing through the calculation by hand to get a better understanding of the FFT. The FFT I implemented does not produce the image that should be produced by ImageJ. I was wondering if anyone can see where I am going wrong with my calculation. My work is here:
>>
>> http://stackoverflow.com/questions/24291896/fft-image-is-not-the-same-as-the-fft-image-produced-in-imagej
>>
>> Thanks.
>>
>> Mariam Dost
>>
>> --
>> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>>
>>
>> --
>> ImageJ mailing list: http://imagej.nih.gov/ij/list.html
>
> --
> ImageJ mailing list: http://imagej.nih.gov/ij/list.html

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html

[Michael Schmid]

Hi everyone,

in the recent daily builds, Wayne has added two convenient ways to access the built-in 1D FHT of ImageJ.

   FHT.transform1D(float[] x)
   FHT.inverseTransform1D(float[] fht)

The FHT (Fast Hartley transform) is closely related to the FFT, and it also needs array sizes that are a power of 2.
http://en.wikipedia.org/wiki/Discrete_Hartley_transform

Actually, it was possible to use the built-in FHT already in previous versions of ImageJ, but that was not so obvious. In the daily build, the 1D FHT also works with array sizes >2^16 (for very large arrays, note that the accuracy is limited to maybe 5 digits or so because it's all 32-bit float, not double-precision arithmetics).

From the FHT, one can easily determine the Fourier (FFT) coefficients for either sine and cosine, or the coefficients of the complex FFT.  Of course, you can also get the power spectrum, see the API documentation of FHT.transform1D.

Michael

--
ImageJ mailing list: http://imagej.nih.gov/ij/list.html