GaussianFilter.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.8.5
// ----------------------------------------------------------------------------
// pcl/GaussianFilter.h - Released 2024-12-28T16:53:48Z
// ----------------------------------------------------------------------------
// 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_GaussianFilter_h
#define __PCL_GaussianFilter_h
 
 
#include <pcl/Defs.h>
 
#include <pcl/KernelFilter.h>
#include <pcl/Math.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
class PCL_CLASS GaussianFilter : public KernelFilter
{
public:
 
   GaussianFilter() = default;
 
   GaussianFilter( float sigma, float epsilon = 0.01, const String& name = String() )
   {
      Initialize( sigma, epsilon, 1, 0 );
      Rename( name );
   }
 
   GaussianFilter( float sigma, float epsilon, float rho, float theta = 0, const String& name = String() )
   {
      Initialize( sigma, epsilon, rho, theta );
      Rename( name );
   }
 
   GaussianFilter( int n, float epsilon = 0.01, const String& name = String() )
   {
      Initialize( n, epsilon, 1, 0 );
      Rename( name );
   }
 
   GaussianFilter( int n, float epsilon, float rho, float theta = 0, const String& name = String() )
   {
      Initialize( n, epsilon, rho, theta );
      Rename( name );
   }
 
   GaussianFilter( const GaussianFilter& ) = default;
 
   GaussianFilter( GaussianFilter&& ) = default;
 
   KernelFilter* Clone() const override
   {
      return new GaussianFilter( *this );
   }
 
   SeparableFilter AsSeparableFilter( float tolerance = __PCL_DEFAULT_FILTER_SEPARABILITY_TOLERANCE ) const override
   {
      if ( m_rho == 1 )
      {
         FVector v = coefficients.RowVector( Size()>>1 );
         return SeparableFilter( v, v );
      }
      return SeparableFilter();
   }
 
   bool IsSeparable() const override
   {
      return m_rho == 1;
   }
 
   GaussianFilter& operator =( const GaussianFilter& x )
   {
      (void)KernelFilter::operator =( x );
      m_sigma   = x.m_sigma;
      m_epsilon = x.m_epsilon;
      m_rho     = x.m_rho;
      m_theta   = x.m_theta;
      return *this;
   }
 
   GaussianFilter& operator =( GaussianFilter&& x )
   {
      (void)KernelFilter::operator =( std::move( x ) );
      m_sigma   = x.m_sigma;
      m_epsilon = x.m_epsilon;
      m_rho     = x.m_rho;
      m_theta   = x.m_theta;
      return *this;
   }
 
   float SigmaX() const
   {
      return m_sigma;
   }
 
   float SigmaY() const
   {
      return m_rho*m_sigma;
   }
 
   float Sigma() const
   {
      return SigmaX();
   }
 
   float Truncation() const
   {
      return m_epsilon;
   }
 
   float AspectRatio() const
   {
      return m_rho;
   }
 
   float RotationAngle() const
   {
      return m_theta;
   }
 
   double FWHMx() const
   {
      return 2.3548200450309493 * m_sigma;
   }
 
   double FWHMy() const
   {
      return m_rho * FWHMx();
   }
 
   double FWHM() const
   {
      return FWHMx();
   }
 
   void Set( float sigma, float epsilon, float rho, float theta )
   {
      Initialize( sigma, epsilon, rho, theta );
   }
 
   void Set( float sigma, float epsilon, float rho )
   {
      Initialize( sigma, epsilon, rho, m_theta );
   }
 
   void Set( float sigma, float epsilon )
   {
      Initialize( sigma, epsilon, m_rho, m_theta );
   }
 
   void Set( float sigma )
   {
      Initialize( sigma, m_epsilon, m_rho, m_theta );
   }
 
   void SetSigma( float sigma )
   {
      Set( sigma );
   }
 
   void SetTruncation( float epsilon )
   {
      Set( m_sigma, epsilon );
   }
 
   void SetAspectRatio( float rho )
   {
      Set( m_sigma, m_epsilon, rho );
   }
 
   void SetRotationAngle( float theta )
   {
      Set( m_sigma, m_epsilon, m_rho, theta );
   }
 
   void Resize( int n ) override
   {
      Initialize( n, m_epsilon, m_rho, m_theta );
   }
 
   friend void Swap( GaussianFilter& x1, GaussianFilter& x2 )
   {
      pcl::Swap( static_cast<KernelFilter&>( x1 ), static_cast<KernelFilter&>( x2 ) );
      pcl::Swap( x1.m_sigma,   x2.m_sigma );
      pcl::Swap( x1.m_rho,     x2.m_rho );
      pcl::Swap( x1.m_theta,   x2.m_theta );
      pcl::Swap( x1.m_epsilon, x2.m_epsilon );
   }
 
private:
 
   float m_sigma = 2;       // standard deviation, horizontal axis
   float m_rho = 1;         // vertical:horizontal axes ratio
   float m_theta = 0;       // rotation angle in radians, [0,+pi]
   float m_epsilon = 0.01F; // maximum truncation error in sigma units
 
   void Initialize( float s, float e, float r, float a )
   {
      PCL_PRECONDITION( s > 0 )
      PCL_PRECONDITION( e > 0 )
      PCL_PRECONDITION( r >= 0 && r <= 1 )
      PCL_PRECONDITION( a >= 0 && a <= Const<float>::pi() )
      m_sigma = Abs( s );
      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( -2*Ln( m_epsilon ) ) ) ) << 1) );
      Rebuild();
   }
 
   void Initialize( int n, float e, float r, float a )
   {
      PCL_PRECONDITION( n == 0 || n >= 3 && (n & 1) != 0 )
      PCL_PRECONDITION( e > 0 )
      PCL_PRECONDITION( r >= 0 && r <= 1 )
      PCL_PRECONDITION( a >= 0 && a <= Const<float>::pi() )
      KernelFilter::Resize( n );
      m_epsilon = Abs( e );
      m_sigma = (Size() >> 1)/Sqrt( -2*Ln( m_epsilon ) );
      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* __restrict__ h = *coefficients;
      double sx = m_sigma;
      double sy = m_rho * sx;
 
      if ( m_theta == 0 || m_rho == 1 )
      {
         double twosx2 = 2*sx*sx;
         double twosy2 = 2*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( Exp( -(dx*dx/twosx2 + dy*dy/twosy2) ) );
      }
      else
      {
         double st, ct; SinCos( double( m_theta ), st, ct );
         double sct    = st*ct;
         double st2    = st*st;
         double ct2    = ct*ct;
         double twosx2 = 2*sx*sx;
         double twosy2 = 2*sy*sy;
         double p1     = ct2/twosx2 + st2/twosy2;
         double p2     = sct/twosy2 - sct/twosx2;
         double p3     = st2/twosx2 + ct2/twosy2;
         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( Exp( -(p1*dx*dx + twop2dy*dx + p3dy2) ) );
         }
      }
   }
};
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif   // __PCL_GaussianFilter_h
 
// ----------------------------------------------------------------------------
// EOF pcl/GaussianFilter.h - Released 2024-12-28T16:53:48Z