RealFFT.h

 
#pragma once
 
#if defined(RealFFT_RECURSES)
#error Recursive header files inclusion detected in RealFFT.h
#else // defined(RealFFT_RECURSES)
#define RealFFT_RECURSES
 
#if !defined RealFFT_h
#define RealFFT_h
 
#ifndef DGTAL_WITH_FFTW3
  #error You need to have activated FFTW3 (DGTAL_WITH_FFTW3) to include this file.
#endif
 
// Inclusions
#include <cstddef>    // std::size_t
 
#include <complex>    // To be included before fftw: see http://www.fftw.org/doc/Complex-numbers.html#Complex-numbers
#include <type_traits>
#include <fftw3.h>
 
#include <boost/math/constants/constants.hpp>
 
#include "DGtal/kernel/domains/HyperRectDomain.h"
#include "DGtal/images/ArrayImageAdapter.h"
#include "DGtal/kernel/PointVector.h"
 
namespace DGtal
{
 
// Implementation details.
namespace detail
{
 
template <typename TFFTW>
struct FFTWComplexCast
  {
    static inline
    typename TFFTW::complex* apply( typename TFFTW::complex* ptr ) { return ptr; }
 
    static inline
    typename TFFTW::complex* apply( std::complex<typename TFFTW::real>* ptr ) { return reinterpret_cast<typename TFFTW::complex*>(ptr); }
  };
 
 
#define FFTW_WRAPPER_GEN(suffix)                                                                                         \
    using size_t  = std::size_t;                                                                                         \
    using complex = fftw ## suffix ## _complex;                                                                          \
    using plan    = fftw ## suffix ## _plan;                                                                             \
    using self    = FFTWWrapper<real>;                                                                                   \
                                                                                                                         \
    static inline void*   malloc( size_t n )      noexcept { return fftw ## suffix ## _malloc(n); }                      \
    static inline void    free( void* p )         noexcept { fftw ## suffix ## _free(p); }                               \
    static inline void    execute( const plan p ) noexcept { fftw ## suffix ## _execute(p); }                            \
    static inline void    destroy_plan( plan p )  noexcept { fftw ## suffix ## _destroy_plan(p); }                       \
                                                                                                                         \
    template < typename C >                                                                                              \
    static inline                                                                                                        \
    plan plan_dft_r2c( int rank, const int* n, real* in, C* out, unsigned flags ) noexcept                               \
      {                                                                                                                  \
        return fftw ## suffix ## _plan_dft_r2c(rank, n, in, FFTWComplexCast<self>::apply(out), flags);                   \
      }                                                                                                                  \
                                                                                                                         \
    template < typename C >                                                                                              \
    static inline                                                                                                        \
    plan plan_dft_c2r( int rank, const int* n, C* in, real* out, unsigned flags ) noexcept                               \
      {                                                                                                                  \
        return fftw ## suffix ## _plan_dft_c2r(rank, n, FFTWComplexCast<self>::apply(in), out, flags);                   \
      }                                                                                                                  \
                                                                                                                         \
    template < typename C >                                                                                              \
    static inline                                                                                                        \
    void execute_dft_r2c( const plan p, real* in, C* out ) noexcept                                                      \
      {                                                                                                                  \
        fftw ## suffix ## _execute_dft_r2c(p, in, FFTWComplexCast<self>::apply(out));                                    \
      }                                                                                                                  \
                                                                                                                         \
    template < typename C >                                                                                              \
    static inline                                                                                                        \
    void execute_dft_c2r( const plan p, C* in, real* out ) noexcept                                                      \
      {                                                                                                                  \
        fftw ## suffix ## _execute_dft_c2r(p, FFTWComplexCast<self>::apply(in), out);                                    \
      }                                                                                                                  \
                                                                                                                         \
                                                                                                                  \
    template < typename C >                                                                                              \
    static inline                                                                                                        \
    plan plan_dft( int rank, const int* n, real* in, C* out, int way, unsigned flags ) noexcept                          \
      {                                                                                                                  \
        if ( way == FFTW_FORWARD )                                                                                       \
          return plan_dft_r2c( rank, n, in, out, flags );                                                                \
        else                                                                                                             \
          return plan_dft_c2r( rank, n, out, in, flags );                                                                \
      }                                                                                                                  \
                                                                                                                         \
                                                                                                                  \
    template < typename C >                                                                                              \
    static inline                                                                                                        \
    void execute_dft( const plan p, real* in, C* out, int way ) noexcept                                                 \
      {                                                                                                                  \
        if ( way == FFTW_FORWARD )                                                                                       \
          execute_dft_r2c( p, in, out );                                                                                 \
        else                                                                                                             \
          execute_dft_c2r( p, out, in );                                                                                 \
      }                                                                                                                  \
 
template <typename Real = double>
struct FFTWWrapper
  {
    // Static errors when using wrong precision or when adapted FFTW3 library is missing.
    static_assert( ! std::is_same<Real, float>::value, "[DGtal::RealFFT] The FFTW3 library for float precision (libfftw3f) has not been found." );
    static_assert( ! std::is_same<Real, double>::value, "[DGtal::RealFFT] The FFTW3 library for double precision (libfftw3) has not been found." );
    static_assert( ! std::is_same<Real, long double>::value, "[DGtal::RealFFT] The FFTW3 library for long double precision (libfftw3l) has not been found." );
    static_assert( ! std::is_same<Real, Real>::value, "[DGtal::RealFFT] Value type not supported." );
  };
 
#ifdef DGTAL_WITH_FFTW3_DOUBLE
template <>
struct FFTWWrapper<double>
  {
    using real = double;
    FFTW_WRAPPER_GEN()
  };
#endif
 
#ifdef DGTAL_WITH_FFTW3_FLOAT
template <>
struct FFTWWrapper<float>
  {
    using real = float;
    FFTW_WRAPPER_GEN(f)
  };
#endif
 
#ifdef DGTAL_WITH_FFTW3_LONG
template <>
struct FFTWWrapper<long double>
  {
    using real = long double;
    FFTW_WRAPPER_GEN(l)
  };
#endif
 
} // detail namespace
 
 
template <
  class TDomain,
  typename T = double
>
class RealFFT;
 
template <typename TSpace, typename T>
class RealFFT< HyperRectDomain<TSpace>, T >
  {
  private:
    using FFTW = detail::FFTWWrapper<T>;
 
    // ----------------------------- Aliases ----------------------------------
  public:
    using Space   = TSpace;                       
    using Real    = T;                            
    using Domain  = HyperRectDomain<Space>;       
    using Point     = typename Space::Point;      
    using RealPoint = PointVector<Space::dimension, Real>;  
    using Dimension = typename Space::Dimension;  
    using Complex = std::complex<Real>;           
    using Self    = RealFFT< Domain, T >;         
 
    using SpatialImage      = ArrayImageAdapter<       Real*, Domain >; 
    using ConstSpatialImage = ArrayImageAdapter< const Real*, Domain >; 
    using FreqImage      = ArrayImageAdapter<       Complex*, Domain >; 
    using ConstFreqImage = ArrayImageAdapter< const Complex*, Domain >; 
 
    static const constexpr Dimension dimension = Domain::dimension; 
 
    // ----------------------- Standard services ------------------------------
  public:
 
    RealFFT( Domain const& aDomain );
 
    RealFFT( Domain const& aDomain, RealPoint const& aLowerBound, RealPoint const& anExtent );
 
    RealFFT( Self const & /* other */ ) = delete;
 
    RealFFT( Self && /* other */ ) = delete;
 
    Self & operator= ( Self const & /* other */ ) = delete;
 
    Self & operator= ( Self && /* other */ ) = delete;
 
    ~RealFFT();
 
    // ----------------------- Interface --------------------------------------
  public:
 
 
 
    std::size_t getPadding() const noexcept;
 
 
          Real* getSpatialStorage()       noexcept;
 
    const Real* getSpatialStorage() const noexcept;
 
 
    SpatialImage      getSpatialImage()       noexcept;
 
    ConstSpatialImage getSpatialImage() const noexcept;
 
    Domain const& getSpatialDomain() const noexcept;
 
    Point  const& getSpatialExtent() const noexcept;
 
 
 
 
          Complex* getFreqStorage()       noexcept;
 
    const Complex* getFreqStorage() const noexcept;
 
 
 
    FreqImage       getFreqImage()       noexcept;
 
    ConstFreqImage  getFreqImage() const noexcept;
 
 
    Domain const& getFreqDomain()    const noexcept;
 
    Point  const& getFreqExtent()    const noexcept;
 
 
 
    void createPlan( unsigned flags, int way );
 
    void doFFT( unsigned flags = FFTW_ESTIMATE, int way = FFTW_FORWARD, bool normalized = false );
 
    void forwardFFT( unsigned flags = FFTW_ESTIMATE );
 
    void backwardFFT( unsigned flags = FFTW_ESTIMATE, bool normalized = true );
 
 
 
    RealPoint getScaledSpatialExtent() const noexcept;
 
    void setScaledSpatialExtent( RealPoint const& anExtent ) noexcept;
 
    RealPoint getScaledSpatialLowerBound() const noexcept;
 
    void setScaledSpatialLowerBound( RealPoint const& aPoint ) noexcept;
 
 
 
 
    RealPoint calcScaledSpatialCoords( Point const& aPoint ) const noexcept;
 
    RealPoint calcScaledFreqCoords( Point const& aPoint ) const noexcept;
 
    Complex calcScaledFreqValue( RealPoint const& aScaledPoint, Complex const& aValue ) const noexcept;
 
 
    Point calcNativeSpatialCoords( RealPoint const& aScaledPoint ) const noexcept;
 
    Point calcNativeFreqCoords( RealPoint const& aScaledPoint, bool & applyConj ) const noexcept;
 
    Complex calcNativeFreqValue( RealPoint const& aScaledPoint, Complex const& aScaledValue ) const noexcept;
 
 
 
 
    inline
    Real getScaledSpatialValue( RealPoint const& aScaledPoint ) const noexcept
      {
        return getSpatialImage()( getNativeSpatialCoords( aScaledPoint ) );
      }
 
    inline
    void setScaledSpatialValue( RealPoint const& aScaledPoint, Real aValue ) noexcept
      {
        getSpatialImage().setValue( getNativeSpatialCoords( aScaledPoint ), aValue );
      }
 
 
 
    inline
    Complex getScaledFreqValue( RealPoint const& aScaledPoint ) const noexcept
      {
        bool apply_conj;
        const auto aPoint = calcNativeFreqCoords( aScaledPoint, apply_conj );
        const Complex aScaledValue = calcScaledFreqValue( aScaledPoint, getFreqImage()( aPoint ) );
        return apply_conj ? std::conj( aScaledValue ) : aScaledValue;
      }
 
    inline
    void setScaledFreqValue( RealPoint const& aScaledPoint, Complex aScaledValue ) noexcept
      {
        bool apply_conj;
        const auto aPoint = calcNativeFreqCoords( aScaledPoint, apply_conj );
        const Complex aValue = calcNativeFreqValue( aScaledPoint, aScaledValue );
 
        if ( apply_conj )
          getFreqImage().setValue( aPoint, std::conj( aValue ) );
        else
          getFreqImage().setValue( aPoint, aValue );
      }
 
 
 
    bool isValid() const noexcept;
 
    void selfDisplay ( std::ostream & out ) const;
 
 
    // ------------------------- Public Datas --------------------------------
  public:
    const Real pi = boost::math::constants::pi<Real>(); 
 
    // ------------------------- Private Datas --------------------------------
  private:
    const Domain  mySpatialDomain;  
    const Point   mySpatialExtent;  
    const Point   myFreqExtent;     
    const Domain  myFreqDomain;     
    const Domain  myFullSpatialDomain;  
          void*   myStorage;        
        RealPoint myScaledSpatialExtent;      
        RealPoint myScaledSpatialLowerBound;  
        Real      myScaledFreqMag;  
 
  };
 
  template <
    class TDomain,
    typename T
  >
  std::ostream&
  operator<< ( std::ostream & out, const RealFFT<TDomain,T> & object );
 
} // namespace DGtal
 
// Includes inline functions.
#include "DGtal/math/RealFFT.ih"
 
//                                                                           //
 
#endif // !defined RealFFT_h
 
#undef RealFFT_RECURSES
#endif // else defined(RealFFT_RECURSES)