FFT2D.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.6.5
// ----------------------------------------------------------------------------
// pcl/FFT2D.h - Released 2024-01-13T15:47:58Z
// ----------------------------------------------------------------------------
// This file is part of the PixInsight Class Library (PCL).
// PCL is a multiplatform C++ framework for development of PixInsight modules.
//
// Copyright (c) 2003-2024 Pleiades Astrophoto S.L. All Rights Reserved.
//
// Redistribution and use in both source and binary forms, with or without
// modification, is permitted provided that the following conditions are met:
//
// 1. All redistributions of source code must retain the above copyright
//    notice, this list of conditions and the following disclaimer.
//
// 2. All redistributions in binary form must reproduce the above copyright
//    notice, this list of conditions and the following disclaimer in the
//    documentation and/or other materials provided with the distribution.
//
// 3. Neither the names "PixInsight" and "Pleiades Astrophoto", nor the names
//    of their contributors, may be used to endorse or promote products derived
//    from this software without specific prior written permission. For written
//    permission, please contact info@pixinsight.com.
//
// 4. All products derived from this software, in any form whatsoever, must
//    reproduce the following acknowledgment in the end-user documentation
//    and/or other materials provided with the product:
//
//    "This product is based on software from the PixInsight project, developed
//    by Pleiades Astrophoto and its contributors (https://pixinsight.com/)."
//
//    Alternatively, if that is where third-party acknowledgments normally
//    appear, this acknowledgment must be reproduced in the product itself.
//
// THIS SOFTWARE IS PROVIDED BY PLEIADES ASTROPHOTO AND ITS CONTRIBUTORS
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
// TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PLEIADES ASTROPHOTO OR ITS
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, BUSINESS
// INTERRUPTION; PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; AND LOSS OF USE,
// DATA OR PROFITS) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
// ----------------------------------------------------------------------------
 
#ifndef __PCL_FFT2D_h
#define __PCL_FFT2D_h
 
 
/*
 * The FFT routines used in this PCL version are based on the KISS FFT library
 * by Mark Borgerding: http://sourceforge.net/projects/kissfft/
 *
 * KISS FFT LICENSE INFORMATION
 * ============================
 *
 * Copyright (c) 2003-2009 Mark Borgerding
 *
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * * Redistributions of source code must retain the above copyright notice,
 *   this list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice,
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the author nor the names of any contributors may be used to
 *   endorse or promote products derived from this software without specific
 *   prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 * POSSIBILITY OF SUCH DAMAGE.
 */
 
#include <pcl/Defs.h>
 
#include <pcl/FFT1D.h>
#include <pcl/Matrix.h>
#include <pcl/ParallelProcess.h>
#include <pcl/StatusMonitor.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
class FFT2DBase
{
protected:
 
   static void Transform( int, int, fcomplex*, const fcomplex*, int, StatusMonitor*, bool, int );
   static void Transform( int, int, dcomplex*, const dcomplex*, int, StatusMonitor*, bool, int );
   static void Transform( int, int, fcomplex*, const float*, StatusMonitor*, bool, int );
   static void Transform( int, int, dcomplex*, const double*, StatusMonitor*, bool, int );
   static void Transform( int, int, float*, const fcomplex*, StatusMonitor*, bool, int );
   static void Transform( int, int, double*, const dcomplex*, StatusMonitor*, bool, int );
};
 
// ----------------------------------------------------------------------------
 
template <typename T>
class PCL_CLASS AbstractFFT2D : public FFT2DBase, public ParallelProcess
{
public:
 
   using scalar = T;
 
   using complex = Complex<T>;
 
   using vector = GenericVector<scalar>;
 
   using complex_vector = GenericVector<complex>;
 
   using matrix = GenericMatrix<scalar>;
 
   using complex_matrix = GenericMatrix<complex>;
 
   using transform = complex_matrix;
 
   AbstractFFT2D( int rows, int cols )
      : m_rows( rows )
      , m_cols( cols )
   {
   }
 
   AbstractFFT2D( int rows, int cols, StatusMonitor& monitor )
      : m_rows( rows )
      , m_cols( cols )
      , m_monitor( &monitor )
   {
   }
 
   ~AbstractFFT2D() override
   {
      Release();
   }
 
   int Rows() const
   {
      return m_rows;
   }
 
   int Cols() const
   {
      return m_cols;
   }
 
   int NumberOfElements() const
   {
      return m_rows*m_cols;
   }
 
   transform DFT() const
   {
      return m_dft;
   }
 
   transform& DFT()
   {
      return m_dft;
   }
 
   virtual void Release()
   {
      m_dft = transform();
   }
 
private:
 
   // Transform dimensions
   int m_rows;
   int m_cols;
 
protected:
 
   transform      m_dft;   // DFT of complex or real data
   StatusMonitor* m_monitor = nullptr;
};
 
// ----------------------------------------------------------------------------
 
#define m_dft           this->m_dft
#define m_rows          this->Rows()
#define m_cols          this->Cols()
#define m_monitor       this->m_monitor
#define m_parallel      this->m_parallel
#define m_maxProcessors this->m_maxProcessors
 
// ----------------------------------------------------------------------------
 
template <typename T>
class PCL_CLASS GenericFFT2D : public AbstractFFT2D<T>
{
public:
 
   using base = AbstractFFT2D<T>;
 
   using scalar = typename base::scalar;
 
   using complex = typename base::complex;
 
   using vector = typename base::vector;
 
   using complex_vector = typename base::complex_vector;
 
   using matrix = typename base::matrix;
 
   using complex_matrix = typename base::complex_matrix;
 
   using transform = typename base::transform;
 
   GenericFFT2D( int rows, int cols )
      : base( rows, cols )
   {
   }
 
   GenericFFT2D( int rows, int cols, StatusMonitor& status )
      : base( rows, cols, status )
   {
   }
 
   ~GenericFFT2D() override
   {
   }
 
   GenericFFT2D& operator <<( const complex* x )
   {
      if ( m_dft.IsEmpty() )
         m_dft = transform( m_rows, m_cols );
      this->Transform( m_rows, m_cols, *m_dft, x, PCL_FFT_FORWARD, m_monitor, m_parallel, m_maxProcessors );
      return *this;
   }
 
   GenericFFT2D& operator >>( complex* y ) const
   {
      if ( m_dft.IsEmpty() )
         throw Error( "Invalid out-of-place inverse FFT: No FFT has been performed." );
      this->Transform( m_rows, m_cols, y, *m_dft, PCL_FFT_BACKWARD, m_monitor, m_parallel, m_maxProcessors );
      return *this;
   }
 
   GenericFFT2D& operator <<( const complex_matrix& x )
   {
      PCL_PRECONDITION( x.Rows() == m_rows )
      PCL_PRECONDITION( x.Cols() == m_cols )
      if ( x.Rows() != m_rows || x.Cols() != m_cols )
         throw Error( "Invalid FFT input matrix dimensions." );
      return operator <<( *x );
   }
 
   GenericFFT2D& operator >>( complex_matrix& y ) const
   {
      PCL_PRECONDITION( y.Rows() == m_rows )
      PCL_PRECONDITION( y.Cols() == m_cols )
      if ( y.Rows() != m_rows || y.Cols() != m_cols )
         throw Error( "Invalid FFT output matrix dimensions." );
      return operator >>( *y );
   }
 
   GenericFFT2D& operator ()( complex* y, const complex* x, int dir = PCL_FFT_FORWARD ) const
   {
      this->Transform( m_rows, m_cols, y, x, dir, m_monitor, m_parallel, m_maxProcessors );
      return const_cast<GenericFFT2D&>( *this );
   }
 
   static int OptimizedLength( int n )
   {
      return GenericFFT<T>::OptimizedLength( n );
   }
};
 
// ----------------------------------------------------------------------------
 
template <typename T>
class PCL_CLASS GenericRealFFT2D : public AbstractFFT2D<T>
{
public:
 
   using base = AbstractFFT2D<T>;
 
   using scalar = typename base::scalar;
 
   using complex = typename base::complex;
 
   using vector = typename base::vector;
 
   using complex_vector = typename base::complex_vector;
 
   using matrix = typename base::matrix;
 
   using complex_matrix = typename base::complex_matrix;
 
   using transform = typename base::transform;
 
   GenericRealFFT2D( int rows, int cols )
      : base( rows, cols )
   {
   }
 
   GenericRealFFT2D( int rows, int cols, StatusMonitor& status )
      : base( rows, cols, status )
   {
   }
 
   ~GenericRealFFT2D() override
   {
   }
 
   GenericRealFFT2D& operator <<( const scalar* x )
   {
      if ( m_dft.IsEmpty() )
         m_dft = transform( m_rows, m_cols/2 + 1 );
      this->Transform( m_rows, m_cols, *m_dft, x, m_monitor, m_parallel, m_maxProcessors );
      return *this;
   }
 
   GenericRealFFT2D& operator >>( scalar* y ) const
   {
      if ( m_dft.IsEmpty() )
         throw Error( "Invalid out-of-place inverse FFT: No FFT has been performed." );
      this->Transform( m_rows, m_cols, y, *m_dft, m_monitor, m_parallel, m_maxProcessors );
      return *this;
   }
 
   GenericRealFFT2D& operator <<( const matrix& x )
   {
      PCL_PRECONDITION( x.Rows() == m_rows )
      PCL_PRECONDITION( x.Cols() == m_cols )
      if ( x.Rows() != m_rows || x.Cols() != m_cols )
         throw Error( "Invalid FFT input matrix dimensions." );
      return operator <<( *x );
   }
 
   GenericRealFFT2D& operator >>( matrix& y ) const
   {
      PCL_PRECONDITION( y.Rows() == m_rows )
      PCL_PRECONDITION( y.Cols() == m_cols )
      if ( y.Rows() != m_rows || y.Cols() != m_cols )
         throw Error( "Invalid FFT output matrix dimensions." );
      return operator >>( *y );
   }
 
   GenericRealFFT2D& operator()( complex* y, const scalar* x ) const
   {
      this->Transform( m_rows, m_cols, y, x, m_monitor, m_parallel, m_maxProcessors );
      return const_cast<GenericRealFFT2D&>( *this );
   }
 
   GenericRealFFT2D& operator()( scalar* y, const complex* x ) const
   {
      this->Transform( m_rows, m_cols, y, x, m_monitor, m_parallel, m_maxProcessors );
      return const_cast<GenericRealFFT2D&>( *this );
   }
 
   static int OptimizedLength( int n )
   {
      return GenericFFT<T>::OptimizedLength( n );
   }
};
 
// ----------------------------------------------------------------------------
 
#undef m_dft
#undef m_rows
#undef m_cols
#undef m_monitor
#undef m_parallel
#undef m_maxProcessors
 
// ----------------------------------------------------------------------------
 
#ifndef __PCL_NO_FFT2D_INSTANTIATE
 
using FFFT2D = GenericFFT2D<float>;
 
using DFFT2D = GenericFFT2D<double>;
 
using FRealFFT2D = GenericRealFFT2D<float>;
 
using DRealFFT2D = GenericRealFFT2D<double>;
 
using FFT2D = FFFT2D;
 
using RealFFT2D = FRealFFT2D;
 
#endif // __PCL_NO_FFT2D_INSTANTIATE
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif   // __PCL_FFT2D_h
 
// ----------------------------------------------------------------------------
// EOF pcl/FFT2D.h - Released 2024-01-13T15:47:58Z