MoffatFilter.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/MoffatFilter.h - Released 2024-06-18T15:48:54Z
// ----------------------------------------------------------------------------
// 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_MoffatFilter_h
#define __PCL_MoffatFilter_h
 
 
#include <pcl/Defs.h>
 
#include <pcl/KernelFilter.h>
#include <pcl/Math.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS MoffatFilter : public KernelFilter
{
public:
 
   MoffatFilter() = default;
 
   MoffatFilter( float sigma, float beta = 4, float epsilon = 0.01, const String& name = String() )
   {
      Initialize( sigma, beta, epsilon, 1, 0 );
      Rename( name );
   }
 
   MoffatFilter( float sigma, float beta, float epsilon, float rho, float theta = 0, const String& name = String() )
   {
      Initialize( sigma, beta, epsilon, rho, theta );
      Rename( name );
   }
 
   MoffatFilter( int n, float beta = 4, float epsilon = 0.01, const String& name = String() )
   {
      Initialize( n, beta, epsilon, 1, 0 );
      Rename( name );
   }
 
   MoffatFilter( int n, float beta, float epsilon, float rho, float theta = 0, const String& name = String() )
   {
      Initialize( n, beta, epsilon, rho, theta );
      Rename( name );
   }
 
   MoffatFilter( const MoffatFilter& ) = default;
 
   MoffatFilter( MoffatFilter&& ) = default;
 
   KernelFilter* Clone() const override
   {
      return new MoffatFilter( *this );
   }
 
   SeparableFilter AsSeparableFilter( float tolerance = __PCL_DEFAULT_FILTER_SEPARABILITY_TOLERANCE ) const override
   {
      return SeparableFilter();
   }
 
   bool IsSeparable() const override
   {
      return false;
   }
 
   MoffatFilter& operator =( const MoffatFilter& ) = default;
 
   MoffatFilter& operator =( MoffatFilter&& ) = default;
 
   float SigmaX() const
   {
      return m_sigma;
   }
 
   float SigmaY() const
   {
      return m_rho*m_sigma;
   }
 
   float Sigma() const
   {
      return SigmaX();
   }
 
   float Beta() const
   {
      return m_beta;
   }
 
   float Truncation() const
   {
      return m_epsilon;
   }
 
   float AspectRatio() const
   {
      return m_rho;
   }
 
   float RotationAngle() const
   {
      return m_theta;
   }
 
   double FWHMx() const
   {
      return 2 * Sqrt( Pow2( 1/m_beta ) - 1 ) * m_sigma;
 
   }
 
   double FWHMy() const
   {
      return m_rho * FWHMx();
 
   }
 
   double FWHM() const
   {
      return FWHMx();
   }
 
   void Set( float sigma, float beta, float epsilon, float rho, float theta )
   {
      Initialize( sigma, beta, epsilon, rho, theta );
   }
 
   void Set( float sigma, float beta, float epsilon, float rho )
   {
      Initialize( sigma, beta, epsilon, rho, m_theta );
   }
 
   void Set( float sigma, float beta, float epsilon )
   {
      Initialize( sigma, beta, epsilon, m_rho, m_theta );
   }
 
   void Set( float sigma, float beta )
   {
      Initialize( sigma, beta, m_epsilon, m_rho, m_theta );
   }
 
   void Set( float sigma )
   {
      Initialize( sigma, m_beta, m_epsilon, m_rho, m_theta );
   }
 
   void SetSigma( float sigma )
   {
      Set( sigma );
   }
 
   void SetBeta( float beta )
   {
      Set( m_sigma, beta );
   }
 
   void SetTruncation( float epsilon )
   {
      Set( m_sigma, m_beta, epsilon );
   }
 
   void SetAspectRatio( float rho )
   {
      Set( m_sigma, m_beta, m_epsilon, rho );
   }
 
   void SetRotationAngle( float theta )
   {
      Set( m_sigma, m_beta, m_epsilon, m_rho, theta );
   }
 
   void Resize( int n ) override
   {
      Initialize( n, m_beta, m_epsilon, m_rho, m_theta );
   }
 
private:
 
   float m_sigma = 2.0F;    // standard deviation, horizontal axis
   float m_beta = 4.0F;     // beta exponent
   float m_rho = 1.0F;      // vertical:horizontal axis ratio
   float m_theta = 0.0F;    // rotation angle in radians, [0,+pi]
   float m_epsilon = 0.01F; // maximum truncation error in sigma units
 
   void Initialize( float s, float b, float e, float r, float a )
   {
      PCL_PRECONDITION( s > 0 )
      PCL_PRECONDITION( b > 0 )
      PCL_PRECONDITION( e > 0 )
      PCL_PRECONDITION( r >= 0 && r <= 1 )
      PCL_PRECONDITION( a >= 0 && a <= Const<float>::pi() )
      m_sigma = Abs( s );
      m_beta = Abs( b );
      m_epsilon = Abs( e );
      m_rho = Range( r, 0.0F, 1.0F );
      m_theta = Range( a, 0.0F, Const<float>::pi() );
      KernelFilter::Resize( 1 + (Max( 1, RoundInt( m_sigma*Sqrt( Pow( m_epsilon, -1/m_beta ) - 1 ) ) ) << 1) );
      Rebuild();
   }
 
   void Initialize( int n, float b, float e, float r, float a )
   {
      PCL_PRECONDITION( n == 0 || n >= 3 && (n & 1) != 0 )
      PCL_PRECONDITION( b > 0 )
      PCL_PRECONDITION( e > 0 )
      PCL_PRECONDITION( r >= 0 && r <= 1 )
      PCL_PRECONDITION( a >= 0 && a <= Const<float>::pi() )
      KernelFilter::Resize( n );
      m_beta = Abs( b );
      m_epsilon = Abs( e );
      m_sigma = (Size() >> 1)/Sqrt( Pow( m_epsilon, -1/m_beta ) - 1 );
      m_rho = Range( r, 0.0F, 1.0F );
      m_theta = Range( a, 0.0F, Const<float>::pi() );
      Rebuild();
   }
 
   void Rebuild()
   {
      int size = Size();
      if ( size == 0 )
         return;
 
      float* h = *coefficients;
      double sx = m_sigma;
      double sy = m_rho * sx;
      double mb = -m_beta;
 
      if ( m_theta == 0 || m_rho == 1 )
      {
         double sx2 = sx*sx;
         double sy2 = sy*sy;
         for ( int n2 = size >> 1, dy = -n2; dy <= n2; ++dy )
            if ( dy > 0 )
               for ( int dx = 0; dx < size; ++dx, ++h )
                  *h = *(h - ((dy+dy)*size));
            else
               for ( int dx = -n2; dx <= n2; ++dx, ++h )
                  *h = (dx > 0) ? *(h - (dx+dx)) : float( Pow( 1 + dx*dx/sx2 + dy*dy/sy2, mb ) );
      }
      else
      {
         double st, ct; SinCos( double( m_theta ), st, ct );
         double sct = st*ct;
         double st2 = st*st;
         double ct2 = ct*ct;
         double sx2 = sx*sx;
         double sy2 = sy*sy;
         double p1  = ct2/sx2 + st2/sy2;
         double p2  = sct/sy2 - sct/sx2;
         double p3  = st2/sx2 + ct2/sy2;
         for ( int n2 = size >> 1, dy = -n2; dy <= n2; ++dy )
         {
            double twop2dy = 2*p2*dy;
            double p3dy2   = p3*dy*dy;
            for ( int dx = -n2; dx <= n2; ++dx, ++h )
               *h = float( Pow( 1 + p1*dx*dx + twop2dy*dx + p3dy2, mb ) );
         }
      }
   }
};
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif   // __PCL_MoffatFilter_h
 
// ----------------------------------------------------------------------------
// EOF pcl/MoffatFilter.h - Released 2024-06-18T15:48:54Z