Real-data in-place Fast Fourier Transform using FFTW3. #1185

rolanddenis · 2016-06-01T20:22:26Z

PR Description

This adds the RealFFT class that features real(spatial) to complex(frequency) forward and backward Discrete Fourier Transform for any dimension, using FFTW3 library (http://www.fftw.org/).

Real to complex transformations allows us to make the transformation in-place in a memory area allocated by fftw3. This memory area is accessible through a real ArrayImageAdapter for the spatial image, and a complex ArrayImageAdapter for the frequency image.

The class also features:

float, double and long double value precisions if the corresponding fftw3 library files are found (fftw3f, fftw3 and fftw3l respectively).
multi-threaded FFT using (for example) pthread (appropriate fftw3 initialization is needed, see http://www.fftw.org/fftw3_doc/Usage-of-Multi_002dthreaded-FFTW.html#Usage-of-Multi_002dthreaded-FFTW).
possibility to choose the planner algorithm (see http://www.fftw.org/fftw3_doc/Planner-Flags.html#Planner-Flags).
some tools to convert digital coordinates into properly scaled spatial and frequency coordinates.
normalized or unnormalized backward transformation (see http://www.fftw.org/fftw3_doc/The-1d-Real_002ddata-DFT.html#The-1d-Real_002ddata-DFT).

Future plans: classes to make convolutions and to apply diffusion operator on an image.

Checklist

Unit-test of your feature with Catch.
Doxygen documentation of the code completed (classes, methods, types, members...)
~~Documentation module page added or updated.~~
New entry in the ChangeLog.md added.
~~Small benchmark of the various planner algorithms (as an example of the importance of this choice).~~
Travis update to test this class.
No warning raised in Debug cmake mode (otherwise, Travis C.I. will fail).
All continuous integration tests pass (Travis & appveyor)

And the FindFFTW3 script, etc. The FFTW3 available precisions are detected and associated macro are defined: WITH_FFTW3_FLOAT, WITH_FFTW3_DOUBLE and WITH_FFTW3_LONG. Threaded FFTW3 is also checked.

JacquesOlivierLachaud · 2016-06-01T20:30:14Z

Looks so cool !

rolanddenis · 2016-06-01T20:39:26Z

😎

But should be a static member (test).

dcoeurjo · 2016-06-06T06:39:03Z

WITH_FFTW3 must also be added to the DGtalConfig.cmake.in when defined.

dcoeurjo · 2016-06-06T06:43:25Z

I'm having an issue with this simple example: https://gist.github.com/dcoeurjo/4209b140b0e5d7ba290b48dc356d507e

(just adding uniform samples and computing the spatial integral contained in the (0,0,0) of the spectrum).

Here is my output:

     In file included from /Users/davidcoeurjolly/Projects/BlueLDS/src/DGtal-example/fourier3D.cpp:8:
In file included from /Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.h:471:
/Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.ih:62:18: error: implicit instantiation of undefined template
      'DGtal::detail::FFTWWrapper<double>'
    , myStorage( FFTW::malloc( sizeof(Complex) * myFreqDomain.size() ) )
                 ^
/Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.ih:73:7: note: in instantiation of member function
      'DGtal::RealFFT<DGtal::HyperRectDomain<DGtal::SpaceND<3, int> >, double>::RealFFT' requested here
    : RealFFT( aDomain
      ^
/Users/davidcoeurjolly/Projects/BlueLDS/src/DGtal-example/fourier3D.cpp:20:19: note: in instantiation of member function
      'DGtal::RealFFT<DGtal::HyperRectDomain<DGtal::SpaceND<3, int> >, double>::RealFFT' requested here
  RealFFT<Domain> fft( domain );
                  ^
/Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.h:150:8: note: template is declared here
struct FFTWWrapper;
       ^
In file included from /Users/davidcoeurjolly/Projects/BlueLDS/src/DGtal-example/fourier3D.cpp:8:
In file included from /Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.h:471:
/Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.ih:62:22: error: incomplete definition of type
      'DGtal::detail::FFTWWrapper<double>'
    , myStorage( FFTW::malloc( sizeof(Complex) * myFreqDomain.size() ) )
                 ~~~~^~
/Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.ih:86:3: error: implicit instantiation of undefined template
      'DGtal::detail::FFTWWrapper<double>'
  FFTW::free( myStorage );
  ^
/Users/davidcoeurjolly/Projects/BlueLDS/src/DGtal-example/fourier3D.cpp:20:19: note: in instantiation of member function
      'DGtal::RealFFT<DGtal::HyperRectDomain<DGtal::SpaceND<3, int> >, double>::~RealFFT' requested here
  RealFFT<Domain> fft( domain );
                  ^
/Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.h:150:8: note: template is declared here
struct FFTWWrapper;
       ^
In file included from /Users/davidcoeurjolly/Projects/BlueLDS/src/DGtal-example/fourier3D.cpp:8:
In file included from /Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.h:471:
/Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.ih:86:7: error: incomplete definition of type
      'DGtal::detail::FFTWWrapper<double>'
  FFTW::free( myStorage );
  ~~~~^~
/Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.ih:336:12: error: implicit instantiation of undefined template
      'DGtal::detail::FFTWWrapper<double>'
  typename FFTW::plan p;
           ^
/Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.ih:391:3: note: in instantiation of member function
      'DGtal::RealFFT<DGtal::HyperRectDomain<DGtal::SpaceND<3, int> >, double>::doFFT' requested here
  doFFT( flags, FFTW_FORWARD, false );
  ^
/Users/davidcoeurjolly/Projects/BlueLDS/src/DGtal-example/fourier3D.cpp:32:7: note: in instantiation of member function
      'DGtal::RealFFT<DGtal::HyperRectDomain<DGtal::SpaceND<3, int> >, double>::forwardFFT' requested here
  fft.forwardFFT();
      ^
/Users/davidcoeurjolly/Sources/DGtal/src/DGtal/math/RealFFT.h:150:8: note: template is declared here
struct FFTWWrapper;
       ^
5 errors generated.
make[2]: *** [DGtal-example/CMakeFiles/fourier3D.dir/fourier3D.cpp.o] Error 1
make[1]: *** [DGtal-example/CMakeFiles/fourier3D.dir/all] Error 2
make[1]: *** Waiting for unfinished jobs....

any idea ?

dcoeurjo · 2016-06-06T13:01:13Z

BTW, would it be possible you add a documentation (new module Dox) ?
For instance, I'm not familiar with FFTW and I'm not sure to clearly see what the normalisation and scaling things mean.

rolanddenis · 2016-06-06T15:27:36Z

About your error, I think that you do not have (or my cmake script does not find) the fftw3 library for double precision. The libraries name are as follow:

libfftw3f for float precision,
libfftw3 for double precision,
libfftw3l for long double precision.

Check your cmake variable FFTW3_DOUBLE_LIBRARIES if it is empty.

It is related to the FFTWWrapper class that is specialized for double precision, line 157 ( https://github.com/DGtal-team/DGtal/pull/1185/files#diff-ed9c28e96a4137a0b02ec3199b6da6aeR157 ) only if WITH_FFTW3_DOUBLE is defined (depends on FFTW3_DOUBLE_LIBRARIES content).

However, I will see how to display a more explicit message when someone use a precision that is not installed.

Concerning your example, I suggest your to initialize the spatial image since RealFFT doesn't do it for you (to avoid unnecessary performance penalty). An efficient way to do this is:
std::fill( spatial.begin(), spatial.end(), 0. );

Finally, it is already planned to make a module documentation dedicated to this class (a simple class description will be to short to explain how to use fftw efficiently).

rolanddenis · 2016-06-06T15:53:06Z

I have tried to update DGtalConfig.cmake.in but I am not sure of what is needed in this file ... could you check it @dcoeurjo ?

Since most users will probably use RealFFT for only a few transformations and thus ESTIMATE will be faster in that case.

rolanddenis · 2016-06-07T15:10:05Z

I have added clear messages when using RealFFT with an unsupported value type or when the associated FFTW3 library is missing.

dcoeurjo · 2016-06-08T06:52:58Z

cmake/DGtalConfig.cmake.in

@@ -157,6 +157,25 @@ IF(@CGAL_FOUND_DGTAL@)
  SET(WITH_CGAL 1)
 ENDIF(@CGAL_FOUND_DGTAL@)

+IF(@FFTW3_FOUND_DGTAL@)
+  FIND_PACKAGE(FFTW3 REQUIRED)


When you have "found" the libraries/includes, you shouldn't retry a FIND_PACKAGE in the DGtalConfig (see MAGICK section or Cairo).
Path are known and are not supposed to change between the DGtal build and the external project you're compiling.

In my system, there is no FindFFTW3.cmake in the cmake module repository so you may need to ship if with DGtal (which is not what you want).

For some deps, we force a FIND_PACKAGE because the package set specific variables..

For short: remove the find, remove the include_directories and that's it :)

dcoeurjo · 2016-06-08T07:03:03Z

Ok... Here is my FFT3W libs:

   /usr/local/lib/libfftw3f.dylib;/usr/local/lib/libfftw3.dylib;/usr/local/lib/libfftw3l.dylib;/usr/local/lib/libfftw3f_threads.dylib;/usr/local/lib/libfftw3_threads.dylib;/usr/local/lib/libfftw3l_threads.dylib.

(default hombrew build).

When I instantiate FFTW on long double, I still get the assert error "[DGtal::RealFFT] The FFTW3 library for long double precision (libfftw3l) has not been found."

and the same for "double" or "float"

rolanddenis · 2016-06-08T07:45:39Z

I think it is also related to the DGtalConfig.cmake.in: the FFTWWrapper specializations are available only if the corresponding WITH_FFTW3_(FLOAT|DOUBLE|LONG) is defined. However, they are defined in FindFFTW3.cmake depending on which fftw3 libraries he found.

So, I see the following alternatives:

we ship the FindFFTW3.cmake file in order to detect which value precisions are available at compilation.
we add @FFTW3_FLOAT_FOUND_DGTAL@ (and the two others) variables in order to define the WITH_FFTW3_FLOAT macros in `DGtalConfig.cmake'. However, DGtal needs to be re-compiled after installing a missing fftw3 library.
we remove the conditional specialization of FFTWWrapper depending on WITH_FFTW3_(FLOAT|DOUBLE|LONG) and we rely on linking errors if a library is missing.
we require all the fftw3 libraries.

What do you think ?

dcoeurjo · 2016-06-08T07:59:30Z

The easiest way would be to propagate the WITH_FFTW_xxxx defines to the DGtalConfig.cmake (solution 2).

dcoeurjo · 2016-06-24T12:04:13Z

Hi @rolanddenis I'm planing a release soon. Can we merge this PR (and postpone the moduleRealFFT.dox) ?

rolanddenis · 2016-06-24T18:10:16Z

Ok to postpone the module page but I also want to slightly change the interface in order to ease some usages. Should I also postpone it ?

dcoeurjo · 2016-06-25T11:30:06Z

Maybe let's postpone the new API. This feature is not critical at this point and we can schedule a quick 0.9.3 release.

dcoeurjo · 2016-06-27T06:50:19Z

@rolanddenis can you please add a changelog entry ?

rolanddenis · 2016-06-27T07:00:05Z

Done

dcoeurjo · 2016-06-27T07:00:47Z

thanks merging.

rolanddenis added 7 commits June 1, 2016 10:38

Copying file from DigitalSnow repository.

171915f

Updating RealFFT to meet DGtal source showing.

9ff398c

Adding WITH_FFTW3 option in cmake file.

41c70c0

And the FindFFTW3 script, etc. The FFTW3 available precisions are detected and associated macro are defined: WITH_FFTW3_FLOAT, WITH_FFTW3_DOUBLE and WITH_FFTW3_LONG. Threaded FFTW3 is also checked.

Normalization option for backward FFT.

e0426ba

Enable distributed doxygen comments for grouped members.

f551b6c

Adding scaling for spatial and frequency domain.

9fb81fe

Draft test-file for RealFFT.

e23ae96

dcoeurjo added the Image label Jun 2, 2016

dcoeurjo added this to the 0.9.2 milestone Jun 2, 2016

rolanddenis added 8 commits June 3, 2016 14:02

RealFFT test file current status.

0a177b8

Fixing RealFFT class and test file.

e1795bb

Merge branch 'master' into pr_RealFFT

ed5e836

Fixing RealFFT and his test file. Adding translation test.

c72d9ae

Making RealPoint alias dependant of the value type in RealFFT.

3f715c7

Adding FFTW3 support in Travis.

8061d38

Making pi a RealFFT public member.

323dadb

But should be a static member (test).

Fixing operator << overload for RealFFT.

5d99dfb

dcoeurjo self-assigned this Jun 6, 2016

rolanddenis added 2 commits June 6, 2016 17:37

Fixing clang warnings with Catch REQUIRE macro in testRealFFT.

18a7d13

Try to update DGtalConfig.cmake.in.

2d483cf

rolanddenis added 3 commits June 7, 2016 12:23

Changing default planner to ESTIMATE.

d3c0c1a

Since most users will probably use RealFFT for only a few transformations and thus ESTIMATE will be faster in that case.

Documenting RealFFT.

567ccce

Documentatin and adding createPlan methods in RealFFT.

c9b7bef

Clear message when FFTW3 library is missing for a given value precision.

8c2fb98

dcoeurjo reviewed Jun 8, 2016
View reviewed changes

rolanddenis added 2 commits June 8, 2016 19:04

Modifying cmake files.

2a03519

Documenting.

57379db

dcoeurjo added the waiting author feedbacks label Jun 13, 2016

rolanddenis added 4 commits June 13, 2016 13:57

Adding scale->native calculators in RealFFT.

36b064b

Fixing scaling calculators.

ea7aacc

Merge branch 'master' into pr_RealFFT

8eeb65d

Adding scaled value accessors.

03f61c3

rolanddenis removed the waiting author feedbacks label Jun 13, 2016

rolanddenis added 3 commits June 13, 2016 21:13

Fixing test file.

bef0650

Fixing wrong freq coordinates when using Hermitian symetry.

4ea853e

Fixing the previous fix.

57301a2

dcoeurjo added the waiting author feedbacks label Jun 24, 2016

Updating ChangeLog.md

8ab38e9

dcoeurjo merged commit 2be6936 into DGtal-team:master Jun 27, 2016

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Real-data in-place Fast Fourier Transform using FFTW3. #1185

Real-data in-place Fast Fourier Transform using FFTW3. #1185

rolanddenis commented Jun 1, 2016 •

edited

Loading

JacquesOlivierLachaud commented Jun 1, 2016

rolanddenis commented Jun 1, 2016

dcoeurjo commented Jun 6, 2016

dcoeurjo commented Jun 6, 2016

dcoeurjo commented Jun 6, 2016

rolanddenis commented Jun 6, 2016 •

edited

Loading

rolanddenis commented Jun 6, 2016 •

edited

Loading

rolanddenis commented Jun 7, 2016

dcoeurjo Jun 8, 2016

dcoeurjo Jun 8, 2016

dcoeurjo commented Jun 8, 2016

rolanddenis commented Jun 8, 2016

dcoeurjo commented Jun 8, 2016

dcoeurjo commented Jun 24, 2016

rolanddenis commented Jun 24, 2016

dcoeurjo commented Jun 25, 2016

dcoeurjo commented Jun 27, 2016

rolanddenis commented Jun 27, 2016

dcoeurjo commented Jun 27, 2016

Real-data in-place Fast Fourier Transform using FFTW3. #1185

Real-data in-place Fast Fourier Transform using FFTW3. #1185

Conversation

rolanddenis commented Jun 1, 2016 • edited Loading

PR Description

Checklist

JacquesOlivierLachaud commented Jun 1, 2016

rolanddenis commented Jun 1, 2016

dcoeurjo commented Jun 6, 2016

dcoeurjo commented Jun 6, 2016

dcoeurjo commented Jun 6, 2016

rolanddenis commented Jun 6, 2016 • edited Loading

rolanddenis commented Jun 6, 2016 • edited Loading

rolanddenis commented Jun 7, 2016

dcoeurjo Jun 8, 2016

Choose a reason for hiding this comment

dcoeurjo Jun 8, 2016

Choose a reason for hiding this comment

dcoeurjo commented Jun 8, 2016

rolanddenis commented Jun 8, 2016

dcoeurjo commented Jun 8, 2016

dcoeurjo commented Jun 24, 2016

rolanddenis commented Jun 24, 2016

dcoeurjo commented Jun 25, 2016

dcoeurjo commented Jun 27, 2016

rolanddenis commented Jun 27, 2016

dcoeurjo commented Jun 27, 2016

rolanddenis commented Jun 1, 2016 •

edited

Loading

rolanddenis commented Jun 6, 2016 •

edited

Loading

rolanddenis commented Jun 6, 2016 •

edited

Loading