DGtal  1.5.beta
RealFFT.h
1 
17 #pragma once
18 
29 #if defined(RealFFT_RECURSES)
30 #error Recursive header files inclusion detected in RealFFT.h
31 #else // defined(RealFFT_RECURSES)
33 #define RealFFT_RECURSES
34 
35 #if !defined RealFFT_h
37 #define RealFFT_h
38 
39 #ifndef WITH_FFTW3
40  #error You need to have activated FFTW3 (WITH_FFTW3) to include this file.
41 #endif
42 
44 // Inclusions
45 #include <cstddef> // std::size_t
46 
47 #include <complex> // To be included before fftw: see http://www.fftw.org/doc/Complex-numbers.html#Complex-numbers
48 #include <type_traits>
49 #include <fftw3.h>
50 
51 #include <boost/math/constants/constants.hpp>
52 
53 #include "DGtal/kernel/domains/HyperRectDomain.h"
54 #include "DGtal/images/ArrayImageAdapter.h"
55 #include "DGtal/kernel/PointVector.h"
56 
57 namespace DGtal
58 {
59 
60 // Implementation details.
61 namespace detail
62 {
63 
65 template <typename TFFTW>
66 struct FFTWComplexCast
67  {
68  static inline
69  typename TFFTW::complex* apply( typename TFFTW::complex* ptr ) { return ptr; }
70 
71  static inline
72  typename TFFTW::complex* apply( std::complex<typename TFFTW::real>* ptr ) { return reinterpret_cast<typename TFFTW::complex*>(ptr); }
73  };
74 
75 
82 #define FFTW_WRAPPER_GEN(suffix) \
83  using size_t = std::size_t; \
84  using complex = fftw ## suffix ## _complex; \
85  using plan = fftw ## suffix ## _plan; \
86  using self = FFTWWrapper<real>; \
87  \
88  static inline void* malloc( size_t n ) noexcept { return fftw ## suffix ## _malloc(n); } \
89  static inline void free( void* p ) noexcept { fftw ## suffix ## _free(p); } \
90  static inline void execute( const plan p ) noexcept { fftw ## suffix ## _execute(p); } \
91  static inline void destroy_plan( plan p ) noexcept { fftw ## suffix ## _destroy_plan(p); } \
92  \
93  template < typename C > \
94  static inline \
95  plan plan_dft_r2c( int rank, const int* n, real* in, C* out, unsigned flags ) noexcept \
96  { \
97  return fftw ## suffix ## _plan_dft_r2c(rank, n, in, FFTWComplexCast<self>::apply(out), flags); \
98  } \
99  \
100  template < typename C > \
101  static inline \
102  plan plan_dft_c2r( int rank, const int* n, C* in, real* out, unsigned flags ) noexcept \
103  { \
104  return fftw ## suffix ## _plan_dft_c2r(rank, n, FFTWComplexCast<self>::apply(in), out, flags); \
105  } \
106  \
107  template < typename C > \
108  static inline \
109  void execute_dft_r2c( const plan p, real* in, C* out ) noexcept \
110  { \
111  fftw ## suffix ## _execute_dft_r2c(p, in, FFTWComplexCast<self>::apply(out)); \
112  } \
113  \
114  template < typename C > \
115  static inline \
116  void execute_dft_c2r( const plan p, C* in, real* out ) noexcept \
117  { \
118  fftw ## suffix ## _execute_dft_c2r(p, FFTWComplexCast<self>::apply(in), out); \
119  } \
120  \
121  \
125  template < typename C > \
126  static inline \
127  plan plan_dft( int rank, const int* n, real* in, C* out, int way, unsigned flags ) noexcept \
128  { \
129  if ( way == FFTW_FORWARD ) \
130  return plan_dft_r2c( rank, n, in, out, flags ); \
131  else \
132  return plan_dft_c2r( rank, n, out, in, flags ); \
133  } \
134  \
135  \
139  template < typename C > \
140  static inline \
141  void execute_dft( const plan p, real* in, C* out, int way ) noexcept \
142  { \
143  if ( way == FFTW_FORWARD ) \
144  execute_dft_r2c( p, in, out ); \
145  else \
146  execute_dft_c2r( p, out, in ); \
147  } \
148 
150 template <typename Real = double>
151 struct FFTWWrapper
152  {
153  // Static errors when using wrong precision or when adapted FFTW3 library is missing.
154  static_assert( ! std::is_same<Real, float>::value, "[DGtal::RealFFT] The FFTW3 library for float precision (libfftw3f) has not been found." );
155  static_assert( ! std::is_same<Real, double>::value, "[DGtal::RealFFT] The FFTW3 library for double precision (libfftw3) has not been found." );
156  static_assert( ! std::is_same<Real, long double>::value, "[DGtal::RealFFT] The FFTW3 library for long double precision (libfftw3l) has not been found." );
157  static_assert( ! std::is_same<Real, Real>::value, "[DGtal::RealFFT] Value type not supported." );
158  };
159 
160 #ifdef WITH_FFTW3_DOUBLE
164 template <>
165 struct FFTWWrapper<double>
166  {
167  using real = double;
168  FFTW_WRAPPER_GEN()
169  };
170 #endif
171 
172 #ifdef WITH_FFTW3_FLOAT
176 template <>
177 struct FFTWWrapper<float>
178  {
179  using real = float;
180  FFTW_WRAPPER_GEN(f)
181  };
182 #endif
183 
184 #ifdef WITH_FFTW3_LONG
188 template <>
189 struct FFTWWrapper<long double>
190  {
191  using real = long double;
192  FFTW_WRAPPER_GEN(l)
193  };
194 #endif
195 
196 } // detail namespace
197 
199 
205 template <
206  class TDomain,
207  typename T = double
208 >
209 class RealFFT;
211 
216 template <typename TSpace, typename T>
217 class RealFFT< HyperRectDomain<TSpace>, T >
218  {
219  private:
220  using FFTW = detail::FFTWWrapper<T>;
221 
222  // ----------------------------- Aliases ----------------------------------
223  public:
224  using Space = TSpace;
225  using Real = T;
226  using Domain = HyperRectDomain<Space>;
227  using Point = typename Space::Point;
228  using RealPoint = PointVector<Space::dimension, Real>;
229  using Dimension = typename Space::Dimension;
230  using Complex = std::complex<Real>;
231  using Self = RealFFT< Domain, T >;
232 
233  using SpatialImage = ArrayImageAdapter< Real*, Domain >;
234  using ConstSpatialImage = ArrayImageAdapter< const Real*, Domain >;
235  using FreqImage = ArrayImageAdapter< Complex*, Domain >;
236  using ConstFreqImage = ArrayImageAdapter< const Complex*, Domain >;
237 
238  static const constexpr Dimension dimension = Domain::dimension;
239 
240  // ----------------------- Standard services ------------------------------
241  public:
242 
252  RealFFT( Domain const& aDomain );
253 
263  RealFFT( Domain const& aDomain, RealPoint const& aLowerBound, RealPoint const& anExtent );
264 
266  RealFFT( Self const & /* other */ ) = delete;
267 
269  RealFFT( Self && /* other */ ) = delete;
270 
272  Self & operator= ( Self const & /* other */ ) = delete;
273 
275  Self & operator= ( Self && /* other */ ) = delete;
276 
278  ~RealFFT();
279 
280  // ----------------------- Interface --------------------------------------
281  public:
282 
286 
288 
294  std::size_t getPadding() const noexcept;
296 
298 
305  Real* getSpatialStorage() noexcept;
306 
307  const Real* getSpatialStorage() const noexcept;
309 
311 
323  SpatialImage getSpatialImage() noexcept;
324 
325  ConstSpatialImage getSpatialImage() const noexcept;
327 
330  Domain const& getSpatialDomain() const noexcept;
331 
333  Point const& getSpatialExtent() const noexcept;
335 
337 
341 
343 
348  Complex* getFreqStorage() noexcept;
349 
350  const Complex* getFreqStorage() const noexcept;
351 
353 
355 
370  FreqImage getFreqImage() noexcept;
371 
372  ConstFreqImage getFreqImage() const noexcept;
374 
376 
381  Domain const& getFreqDomain() const noexcept;
382 
388  Point const& getFreqExtent() const noexcept;
390 
392 
396 
414  void createPlan( unsigned flags, int way );
415 
442  void doFFT( unsigned flags = FFTW_ESTIMATE, int way = FFTW_FORWARD, bool normalized = false );
443 
459  void forwardFFT( unsigned flags = FFTW_ESTIMATE );
460 
484  void backwardFFT( unsigned flags = FFTW_ESTIMATE, bool normalized = true );
485 
487 
491 
494  RealPoint getScaledSpatialExtent() const noexcept;
495 
499  void setScaledSpatialExtent( RealPoint const& anExtent ) noexcept;
501 
504  RealPoint getScaledSpatialLowerBound() const noexcept;
505 
509  void setScaledSpatialLowerBound( RealPoint const& aPoint ) noexcept;
511 
513 
517 
519 
527  RealPoint calcScaledSpatialCoords( Point const& aPoint ) const noexcept;
528 
542  RealPoint calcScaledFreqCoords( Point const& aPoint ) const noexcept;
543 
554  Complex calcScaledFreqValue( RealPoint const& aScaledPoint, Complex const& aValue ) const noexcept;
556 
558 
566  Point calcNativeSpatialCoords( RealPoint const& aScaledPoint ) const noexcept;
567 
583  Point calcNativeFreqCoords( RealPoint const& aScaledPoint, bool & applyConj ) const noexcept;
584 
595  Complex calcNativeFreqValue( RealPoint const& aScaledPoint, Complex const& aScaledValue ) const noexcept;
597 
599 
603 
605 
611  inline
612  Real getScaledSpatialValue( RealPoint const& aScaledPoint ) const noexcept
613  {
614  return getSpatialImage()( getNativeSpatialCoords( aScaledPoint ) );
615  }
616 
623  inline
624  void setScaledSpatialValue( RealPoint const& aScaledPoint, Real aValue ) noexcept
625  {
626  getSpatialImage().setValue( getNativeSpatialCoords( aScaledPoint ), aValue );
627  }
629 
630 
632 
639  inline
640  Complex getScaledFreqValue( RealPoint const& aScaledPoint ) const noexcept
641  {
642  bool apply_conj;
643  const auto aPoint = calcNativeFreqCoords( aScaledPoint, apply_conj );
644  const Complex aScaledValue = calcScaledFreqValue( aScaledPoint, getFreqImage()( aPoint ) );
645  return apply_conj ? std::conj( aScaledValue ) : aScaledValue;
646  }
647 
655  inline
656  void setScaledFreqValue( RealPoint const& aScaledPoint, Complex aScaledValue ) noexcept
657  {
658  bool apply_conj;
659  const auto aPoint = calcNativeFreqCoords( aScaledPoint, apply_conj );
660  const Complex aValue = calcNativeFreqValue( aScaledPoint, aScaledValue );
661 
662  if ( apply_conj )
663  getFreqImage().setValue( aPoint, std::conj( aValue ) );
664  else
665  getFreqImage().setValue( aPoint, aValue );
666  }
668 
670 
674 
678  bool isValid() const noexcept;
679 
684  void selfDisplay ( std::ostream & out ) const;
685 
687 
688  // ------------------------- Public Datas --------------------------------
689  public:
690  const Real pi = boost::math::constants::pi<Real>();
691 
692  // ------------------------- Private Datas --------------------------------
693  private:
694  const Domain mySpatialDomain;
695  const Point mySpatialExtent;
696  const Point myFreqExtent;
697  const Domain myFreqDomain;
698  const Domain myFullSpatialDomain;
699  void* myStorage;
700  RealPoint myScaledSpatialExtent;
701  RealPoint myScaledSpatialLowerBound;
702  Real myScaledFreqMag;
703 
704  };
705 
712  template <
713  class TDomain,
714  typename T
715  >
716  std::ostream&
717  operator<< ( std::ostream & out, const RealFFT<TDomain,T> & object );
718 
719 } // namespace DGtal
720 
722 // Includes inline functions.
723 #include "DGtal/math/RealFFT.ih"
724 
725 // //
727 
728 #endif // !defined RealFFT_h
729 
730 #undef RealFFT_RECURSES
731 #endif // else defined(RealFFT_RECURSES)
732 
DGtal::HyperRectDomain
Aim: Parallelepidec region of a digital space, model of a 'CDomain'.
Definition: HyperRectDomain.h:100
DGtal::RealFFT< HyperRectDomain< TSpace >, T >::ConstSpatialImage
ArrayImageAdapter< const Real *, Domain > ConstSpatialImage
Constant spatial image type.
Definition: RealFFT.h:228
DGtal::RealFFT< HyperRectDomain< TSpace >, T >::Self
RealFFT< Domain, T > Self
Self type.
Definition: RealFFT.h:225
DGtal::RealFFT< HyperRectDomain< TSpace >, T >::Point
typename Space::Point Point
Point type.
Definition: RealFFT.h:221
DGtal::RealFFT< HyperRectDomain< TSpace >, T >::Dimension
typename Space::Dimension Dimension
Space dimension type.
Definition: RealFFT.h:223
DGtal::RealFFT< HyperRectDomain< TSpace >, T >::SpatialImage
ArrayImageAdapter< Real *, Domain > SpatialImage
Mutable spatial image type.
Definition: RealFFT.h:227
DGtal::RealFFT< HyperRectDomain< TSpace >, T >::Complex
std::complex< Real > Complex
Complex value type.
Definition: RealFFT.h:224
DGtal::RealFFT< HyperRectDomain< TSpace >, T >::RealFFT
RealFFT(Self &&)=delete
Move constructor. Deleted.
DGtal::RealFFT< HyperRectDomain< TSpace >, T >::ConstFreqImage
ArrayImageAdapter< const Complex *, Domain > ConstFreqImage
Constant frequency image type.
Definition: RealFFT.h:230
DGtal::RealFFT< HyperRectDomain< TSpace >, T >::FreqImage
ArrayImageAdapter< Complex *, Domain > FreqImage
Mutable frequency image type.
Definition: RealFFT.h:229
DGtal::SpaceND::Dimension
DGtal::Dimension Dimension
Copy of the type used for the dimension.
Definition: SpaceND.h:129
DGtal
DGtal is the top-level namespace which contains all DGtal functions and types.
boost
Definition: Boost.dox:28
DGtal::detail::FFTWComplexCast
Facility to cast to the complex type used by fftw.
Definition: RealFFT.h:67
DGtal::detail::FFTWComplexCast::apply
static TFFTW::complex * apply(typename TFFTW::complex *ptr)
Definition: RealFFT.h:69
DGtal::detail::FFTWComplexCast::apply
static TFFTW::complex * apply(std::complex< typename TFFTW::real > *ptr)
Definition: RealFFT.h:72
DGtal::detail::FFTWWrapper
Wrapper to fftw functions depending on value type.
Definition: RealFFT.h:146
Point
MyPointD Point
Definition: testClone2.cpp:383
aPoint
const Point aPoint(3, 4)
RealPoint
PointVector< 3, double > RealPoint
Definition: testTriangulatedSurface.cpp:51