ChebyshevFit.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/ChebyshevFit.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_ChebyshevFit_h
#define __PCL_ChebyshevFit_h
 
 
#include <pcl/Defs.h>
 
#include <pcl/Constants.h>
#include <pcl/ErrorHandler.h>
#include <pcl/Exception.h>
#include <pcl/MultiVector.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
template <typename Tx, typename Ty> class GenericScalarChebyshevFit;
 
template <typename Tx, typename Ty>
class GenericChebyshevFit
{
public:
 
   using coefficients = GenericVector<Ty>;
 
   using coefficient_series = GenericMultiVector<Ty>;
 
   using function_value = GenericVector<Ty>;
 
   using function_values = GenericMultiVector<Ty>;
 
   template <class F>
   GenericChebyshevFit( F f, Tx x1, Tx x2, int N, int n )
      : dx( Abs( x2 - x1 ) )
      , x0( (x1 + x2)/2 )
      , m( Max( 2, n ), Max( 1, N ) )
   {
      PCL_PRECONDITION( N > 0 )
      PCL_PRECONDITION( n > 1 )
      if ( 1 + dx == 1 )
         throw Error( "GenericChebyshevFit: Empty or insignificant function evaluation interval." );
      N = m.Length();
      n = m[0];
      Tx dx2 = dx/2;
      function_values y( n, N );
      for ( int j = 0; j < n; ++j )
         y[j] = f( x0 + Cos( Const<Ty>::pi()*(j + 0.5)/n )*dx2 );
      c = coefficient_series( N, n );
      Ty k = 2.0/n;
      for ( int i = 0; i < N; ++i )
         for ( int j = 0; j < n; ++j )
         {
            Ty s = Ty( 0 );
            for ( int k = 0; k < n; ++k )
               s += y[k][i] * Cos( Const<Ty>::pi()*j*(k + 0.5)/n );
            c[i][j] = k*s;
         }
   }
 
   GenericChebyshevFit( const coefficient_series& ck, Tx x1, Tx x2 )
      : dx( Abs( x2 - x1 ) )
      , x0( (x1 + x2)/2 )
      , c( ck )
      , m( int( ck.Length() ) )
   {
      if ( 1 + dx == 1 )
         throw Error( "GenericChebyshevFit: Empty or insignificant function evaluation interval." );
      if ( c.IsEmpty() )
         throw Error( "GenericChebyshevFit: Empty polynomial expansion." );
      for ( int i = 0; i < m.Length(); ++i )
         if ( (m[i] = c[i].Length()) < 1 )
            throw Error( "GenericChebyshevFit: Invalid coefficient series for dimension " + String( i ) + '.' );
   }
 
   GenericChebyshevFit() = default;
 
   GenericChebyshevFit( const GenericChebyshevFit& ) = default;
 
   GenericChebyshevFit( GenericChebyshevFit&& ) = default;
 
   GenericChebyshevFit& operator =( const GenericChebyshevFit& ) = default;
 
   GenericChebyshevFit& operator =( GenericChebyshevFit&& ) = default;
 
   bool IsValid() const
   {
      return !c.IsEmpty();
   }
 
   Tx LowerBound() const
   {
      return x0 - dx/2;
   }
 
   Tx UpperBound() const
   {
      return x0 + dx/2;
   }
 
   int NumberOfComponents() const
   {
      return m.Length();
   }
 
   int Length( int i = 0 ) const
   {
      PCL_PRECONDITION( i >= 0 && i < NumberOfComponents() )
      return c[i].Length();
   }
 
   int TruncatedLength( int i = -1 ) const
   {
      PCL_PRECONDITION( i < 0 || i < NumberOfComponents() )
      if ( i < 0 )
         return m.MaxComponent();
      return m[i];
   }
 
   int NumberOfCoefficients() const
   {
      int N = 0;
      for ( int i = 0; i < NumberOfComponents(); ++i )
         N += c[i].Length();
      return N;
   }
 
   int NumberOfTruncatedCoefficients() const
   {
      return m.Sum();
   }
 
   bool IsTruncated( int i = -1 ) const
   {
      PCL_PRECONDITION( i < 0 || i < NumberOfComponents() )
      if ( i < 0 )
         return m.MaxComponent() < Length();
      return m[i] < c[i].Length();
   }
 
   Ty TruncationError( int i = -1 ) const
   {
      PCL_PRECONDITION( i < 0 || i < NumberOfComponents() )
      if ( i < 0 )
      {
         int N = NumberOfComponents();
         function_value e( N );
         for ( int j = 0; j < N; ++j )
            e[j] = TruncationError( j );
         return e.MaxComponent();
      }
      if ( m[i] < c[i].Length() )
      {
         Ty e = Ty( 0 );
         for ( int j = c[i].Length(); --j >= m[i]; )
            e += Abs( c[i][j] );
         return e;
      }
      return Abs( c[i][c[i].Length()-1] );
   }
 
   bool Truncate( Ty e, int mmin = 2 )
   {
      e = Abs( e );
      mmin = Max( 2, mmin );
      int N = NumberOfComponents();
      int tc = 0;
      for ( int i = 0; i < N; ++i )
      {
         Ty s = Ty( 0 );
         for ( m[i] = c[i].Length(); m[i] > mmin; --m[i] )
            if ( (s += Abs( c[i][m[i]-1] )) >= e )
               break;
         if ( m[i] < c[i].Length() )
            ++tc;
      }
      return tc == N;
   }
 
   bool Truncate( const coefficients& eps, int mmin = 2 )
   {
      mmin = Max( 2, mmin );
      int N = NumberOfComponents();
      int tc = 0;
      for ( int i = 0; i < N; ++i )
      {
         Ty e = Abs( eps[i] );
         Ty s = Ty( 0 );
         for ( m[i] = c[i].Length(); m[i] > mmin; --m[i] )
            if ( (s += Abs( c[i][m[i]-1] )) >= e )
               break;
         if ( m[i] < c[i].Length() )
            ++tc;
      }
      return tc == N;
   }
 
   const coefficient_series& Coefficients() const
   {
      return c;
   }
 
   function_value Evaluate( Tx x ) const
   {
      PCL_PRECONDITION( x >= LowerBound() )
      PCL_PRECONDITION( x <= UpperBound() )
      const Ty y0 = Ty( 2*(x - x0)/dx );
      const Ty y2 = 2*y0;
      function_value y( NumberOfComponents() );
      for ( int i = 0; i < y.Length(); ++i )
      {
         Ty d0 = Ty( 0 );
         Ty d1 = Ty( 0 );
         const Ty* k = c[i].At( m[i] );
         for ( int j = m[i]; --j > 0; )
         {
            Ty d = d1;
            d1 = y2*d1 - d0 + *--k;
            d0 = d;
         }
         y[i] = y0*d1 - d0 + *--k/2;
      }
      return y;
   }
 
   function_value operator ()( Tx x ) const
   {
      return Evaluate( x );
   }
 
   GenericChebyshevFit Derivative() const
   {
      int N = NumberOfComponents();
      GenericChebyshevFit ch1;
      ch1.dx = dx;
      ch1.x0 = x0;
      ch1.c = coefficient_series( N );
      ch1.m = IVector( N );
      for ( int i = 0; i < N; ++i )
      {
         int n = Max( 1, c[i].Length()-1 );
         ch1.c[i] = coefficients( n );
         ch1.m[i] = n;
         if ( n > 1 )
         {
            ch1.c[i][n-1] = 2*n*c[i][n];
            ch1.c[i][n-2] = 2*(n-1)*c[i][n-1];
            for ( int j = n-3; j >= 0; --j )
               ch1.c[i][j] = ch1.c[i][j+2] + 2*(j+1)*c[i][j+1];
            ch1.c[i] *= Ty( 2 )/dx;
         }
         else
            ch1.c[i] = Ty( 0 );
      }
      return ch1;
   }
 
   GenericChebyshevFit Integral() const
   {
      int N = NumberOfComponents();
      GenericChebyshevFit ch;
      ch.dx = dx;
      ch.x0 = x0;
      ch.c = coefficient_series( N );
      ch.m = IVector( N );
      for ( int i = 0; i < N; ++i )
      {
         int n = c[i].Length();
         ch.c[i] = coefficients( n );
         ch.m[i] = n;
         if ( n > 1 )
         {
            const Ty k = Ty( dx )/4;
            Ty s = Ty( 0 );
            int f = 1;
            for ( int j = 1; j < n-1; ++j, f = -f )
               s += f*(ch.c[i][j] = k*(c[i][j-1] - c[i][j+1])/j);
            ch.c[i][0] = 2*(s + f*(ch.c[i][n-1] = k*c[i][n-2]/(n-1)));
         }
         else
            ch.c[i] = Ty( 0 );
      }
      return ch;
   }
 
private:
 
   Tx                 dx; // x2 - x1
   Tx                 x0; // (x1 + x2)/2
   coefficient_series c;  // Chebyshev polynomial coefficients
   IVector            m;  // length of the truncated coefficient series
 
   friend class GenericScalarChebyshevFit<Tx,Ty>;
};
 
// ----------------------------------------------------------------------------
 
template <typename Tx, typename Ty>
class GenericScalarChebyshevFit : public GenericChebyshevFit<Tx,Ty>
{
public:
 
   template <class F>
   GenericScalarChebyshevFit( F f, Tx x1, Tx x2, int n )
      : GenericChebyshevFit<Tx,Ty>( [f]( Tx x ){ return GenericVector<Ty>( f( x ), 1 ); }, x1, x2, 1, n )
   {
   }
 
   GenericScalarChebyshevFit() = default;
 
   GenericScalarChebyshevFit( const GenericScalarChebyshevFit& ) = default;
 
   GenericScalarChebyshevFit( GenericScalarChebyshevFit&& ) = default;
 
   GenericScalarChebyshevFit& operator =( const GenericScalarChebyshevFit& ) = default;
 
   GenericScalarChebyshevFit& operator =( GenericScalarChebyshevFit&& ) = default;
 
   int TruncatedLength() const
   {
      return GenericChebyshevFit<Tx,Ty>::TruncatedLength( 0 );
   }
 
   bool IsTruncated() const
   {
      return GenericChebyshevFit<Tx,Ty>::IsTruncated( 0 );
   }
 
   Ty TruncationError() const
   {
      return GenericChebyshevFit<Tx,Ty>::TruncationError( 0 );
   }
 
   Ty Evaluate( Tx x ) const
   {
      return GenericChebyshevFit<Tx,Ty>::Evaluate( x )[0];
   }
 
   Ty operator ()( Tx x ) const
   {
      return Evaluate( x );
   }
 
   GenericScalarChebyshevFit Derivative() const
   {
      return GenericChebyshevFit<Tx,Ty>::Derivative();
   }
 
   GenericScalarChebyshevFit Integral() const
   {
      return GenericChebyshevFit<Tx,Ty>::Integral();
   }
 
private:
 
   GenericScalarChebyshevFit( const GenericChebyshevFit<Tx,Ty>& T )
      : GenericChebyshevFit<Tx,Ty>()
   {
      this->dx = T.dx;
      this->x0 = T.x0;
      this->c = typename GenericChebyshevFit<Tx,Ty>::coefficient_series( 1 );
      this->c[0] = T.c[0];
      this->m = IVector( 1 );
      this->m[0] = T.m[0];
   }
};
 
// ----------------------------------------------------------------------------
 
#ifndef __PCL_NO_CHEBYSHEV_FIT_INSTANTIATE
 
using F32ChebyshevFit = GenericChebyshevFit<float, float>;
 
using FChebyshevFit = F32ChebyshevFit;
 
using F64ChebyshevFit = GenericChebyshevFit<double, double>;
 
using DChebyshevFit = F64ChebyshevFit;
 
using ChebyshevFit = DChebyshevFit;
 
#ifndef _MSC_VER
 
using F80ChebyshevFit = GenericChebyshevFit<long double, long double>;
 
using LDChebyshevFit = F80ChebyshevFit;
 
#endif   // !_MSC_VER
 
using F32ScalarChebyshevFit = GenericScalarChebyshevFit<float, float>;
 
using FScalarChebyshevFit = F32ScalarChebyshevFit;
 
using F64ScalarChebyshevFit = GenericScalarChebyshevFit<double, double>;
 
using DScalarChebyshevFit = F64ScalarChebyshevFit;
 
using ScalarChebyshevFit = DScalarChebyshevFit;
 
#ifndef _MSC_VER
 
using F80ScalarChebyshevFit = GenericScalarChebyshevFit<long double, long double>;
 
using LDScalarChebyshevFit = F80ScalarChebyshevFit;
 
#endif   // !_MSC_VER
 
#endif   // !__PCL_NO_CHEBYSHEV_FIT_INSTANTIATE
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif  // __PCL_ChebyshevFit_h
 
// ----------------------------------------------------------------------------
// EOF pcl/ChebyshevFit.h - Released 2024-06-18T15:48:54Z