DGtal 1.4.0
Loading...
Searching...
No Matches
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
57namespace DGtal
58{
59
60// Implementation details.
61namespace detail
62{
63
65template <typename TFFTW>
66struct 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
124 \
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
138 \
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
150template <typename Real = double>
151struct 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
164template <>
165struct FFTWWrapper<double>
166 {
167 using real = double;
168 FFTW_WRAPPER_GEN()
169 };
170#endif
171
172#ifdef WITH_FFTW3_FLOAT
176template <>
177struct FFTWWrapper<float>
178 {
179 using real = float;
180 FFTW_WRAPPER_GEN(f)
181 };
182#endif
183
184#ifdef WITH_FFTW3_LONG
188template <>
189struct FFTWWrapper<long double>
190 {
191 using real = long double;
192 FFTW_WRAPPER_GEN(l)
193 };
194#endif
195
196} // detail namespace
197
199
205template <
206 class TDomain,
207 typename T = double
208>
209class RealFFT;
211
216template <typename TSpace, typename T>
217class 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 >::FreqImage
ArrayImageAdapter< Complex *, Domain > FreqImage
Mutable frequency image type.
Definition RealFFT.h:229
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 >::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 >::SpatialImage
ArrayImageAdapter< Real *, Domain > SpatialImage
Mutable spatial image type.
Definition RealFFT.h:227
DGtal::RealFFT< HyperRectDomain< TSpace >, T >::ConstFreqImage
ArrayImageAdapter< const Complex *, Domain > ConstFreqImage
Constant frequency image type.
Definition RealFFT.h:230
DGtal::SpaceND::Point
PointVector< dim, Integer > Point
Points in DGtal::SpaceND.
Definition SpaceND.h:110
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.
Definition ClosedIntegerHalfPlane.h:49
std
STL namespace.
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(std::complex< typename TFFTW::real > *ptr)
Definition RealFFT.h:72
DGtal::detail::FFTWComplexCast::apply
static TFFTW::complex * apply(typename TFFTW::complex *ptr)
Definition RealFFT.h:69
DGtal::detail::FFTWWrapper
Wrapper to fftw functions depending on value type.
Definition RealFFT.h:146
aPoint
const Point aPoint(3, 4)
RealPoint
PointVector< 3, double > RealPoint
Definition testTriangulatedSurface.cpp:51