CubicSplineInterpolation.h

Go to the documentation of this file.

//     ____   ______ __
//    / __ \ / ____// /
//   / /_/ // /    / /
//  / ____// /___ / /___   PixInsight Class Library
// /_/     \____//_____/   PCL 2.7.0
// ----------------------------------------------------------------------------
// pcl/CubicSplineInterpolation.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_CubicSplineInterpolation_h
#define __PCL_CubicSplineInterpolation_h
 
 
#include <pcl/Defs.h>
#include <pcl/Diagnostics.h>
 
#include <pcl/UnidimensionalInterpolation.h>
 
namespace pcl
{
 
// ----------------------------------------------------------------------------
 
#define m_x this->m_x
#define m_y this->m_y
 
template <typename T = double>
class PCL_CLASS CubicSplineInterpolation : public UnidimensionalInterpolation<T>
{
public:
 
   using vector_type = typename UnidimensionalInterpolation<T>::vector_type;
 
   CubicSplineInterpolation() = default;
 
   CubicSplineInterpolation( const CubicSplineInterpolation& ) = default;
 
   CubicSplineInterpolation( CubicSplineInterpolation&& ) = default;
 
   ~CubicSplineInterpolation() override
   {
   }
 
   void GetBoundaryConditions( double& y1, double& yn ) const
   {
      y1 = m_dy1;
      yn = m_dyn;
   }
 
   void SetBoundaryConditions( double y1, double yn )
   {
      Clear();
      m_dy1 = y1;
      m_dyn = yn;
   }
 
   void Initialize( const vector_type& x, const vector_type& y ) override
   {
      if ( y.Length() < 2 )
         throw Error( "CubicSplineInterpolation::Initialize(): Less than two data points specified." );
 
      try
      {
         Clear();
         UnidimensionalInterpolation<T>::Initialize( x, y );
 
         int n = this->Length();
         m_dy2 = DVector( n );
         m_current = -1; // prepare for 1st interpolation
         DVector w( n ); // working vector
 
         if ( m_x )
         {
            /*
             * Cubic splines with explicit x[i] for i = {0,...,n-1}.
             */
            if ( m_dy1 == 0 && m_dyn == 0 )
            {
               /*
                * Natural cubic spline.
                */
               m_dy2[0] = m_dy2[n-1] = w[0] = 0;
 
               for ( int i = 1; i < n-1; ++i )
               {
                  double s = (double( m_x[i] ) - double( m_x[i-1] )) / (double( m_x[i+1] ) - double( m_x[i-1] ));
                  double p = s*m_dy2[i-1] + 2;
                  m_dy2[i] = (s - 1)/p;
                  w[i] = (double( m_y[i+1] ) - double( m_y[i  ] )) / (double( m_x[i+1] ) - double( m_x[i  ] ))
                       - (double( m_y[i  ] ) - double( m_y[i-1] )) / (double( m_x[i  ] ) - double( m_x[i-1] ));
                  w[i] = (6*w[i]/(double( m_x[i+1] ) - double( m_x[i-1] )) - s*w[i-1])/p;
               }
 
               for ( int i = n-2; i > 0; --i ) // N.B. w[0] is not defined
                  m_dy2[i] = m_dy2[i]*m_dy2[i+1] + w[i];
            }
            else
            {
               /*
                * Cubic spline with prescribed end point derivatives.
                */
               w[0] = 3/(double( m_x[1] ) - double( m_x[0] ))
                    * ((double( m_y[1] ) - double( m_y[0] ))/(double( m_x[1] ) - double( m_x[0] )) - m_dy1);
 
               m_dy2[0] = -0.5;
 
               for ( int i = 1; i < n-1; ++i )
               {
                  double s = (double( m_x[i] ) - double( m_x[i-1] )) / (double( m_x[i+1] ) - double( m_x[i-1] ));
                  double p = s*m_dy2[i-1] + 2;
                  m_dy2[i] = (s - 1)/p;
                  w[i] = (double( m_y[i+1] ) - double( m_y[i  ] )) / (double( m_x[i+1] ) - double( m_x[i  ] ))
                       - (double( m_y[i  ] ) - double( m_y[i-1] )) / (double( m_x[i  ] ) - double( m_x[i-1] ));
                  w[i] = (6*w[i]/(double( m_x[i+1] ) - double( m_x[i-1] )) - s*w[i-1])/p;
               }
 
               m_dy2[n-1] = (3/(double( m_x[n-1] ) - double( m_x[n-2] ))
                             * (m_dyn - (double( m_y[n-1] ) - double( m_y[n-2] ))/(double( m_x[n-1] ) - double( m_x[n-2] ))) - w[n-2]/2)
                          / (1 + m_dy2[n-2]/2);
 
               for ( int i = n-2; i >= 0; --i )
                  m_dy2[i] = m_dy2[i]*m_dy2[i+1] + w[i];
            }
         }
         else
         {
            /*
             * Natural cubic spline with
             * implicit x[i] = i for i = {0,1,...,n-1}.
             */
            m_dy2[0] = m_dy2[n-1] = w[0] = 0;
 
            for ( int i = 1; i < n-1; ++i )
            {
               double p = m_dy2[i-1]/2 + 2;
               m_dy2[i] = -0.5/p;
               w[i] = double( m_y[i+1] ) - 2*double( m_y[i] ) + double( m_y[i-1] );
               w[i] = (3*w[i] - w[i-1]/2)/p;
            }
 
            for ( int i = n-2; i > 0; --i ) // N.B. w[0] is not defined
               m_dy2[i] = m_dy2[i]*m_dy2[i+1] + w[i];
         }
      }
      catch ( ... )
      {
         Clear();
         throw;
      }
   }
 
   double operator()( double x ) const override
   {
      PCL_PRECONDITION( IsValid() )
 
      int n = this->Length();
 
      if ( m_x )
      {
         /*
          * Cubic spline with explicit x[i] for i = {0,...,n-1}.
          */
 
         /*
          * Bracket the evaluation point x by binary search of the closest
          * pair of data points, if needed. m_current < 0 signals first-time
          * evaluation since initialization.
          */
         int j0 = m_current, j1;
         if ( j0 < 0 || x < m_x[j0] || m_x[j0+1] < x )
            for ( j0 = 0, j1 = n-1; j1-j0 > 1; )
            {
               int m = (j0 + j1) >> 1;
               if ( x < m_x[m] )
                  j1 = m;
               else
                  j0 = m;
            }
         else
            j1 = j0 + 1;
         m_current = j0;
 
         /*
          * Distance h between the closest neighbors. Will be zero (or
          * insignificant) if two or more x values are equal with respect to
          * the machine epsilon for type T.
          */
         double h = double( m_x[j1] ) - double( m_x[j0] );
         if ( 1 + h == 1 )
            return 0.5*(double( m_y[j0] ) + double( m_y[j1] ));
 
         double a = (double( m_x[j1] ) - x)/h;
         double b = (x - double( m_x[j0] ))/h;
         return a*double( m_y[j0] )
              + b*double( m_y[j1] )
              + ((a*a*a - a)*m_dy2[j0] + (b*b*b - b)*m_dy2[j1])*h*h/6;
      }
      else
      {
         /*
          * Natural cubic spline with implicit x[i] = i for i = {0,1,...,n-1}.
          */
         int j0 = pcl::Range( pcl::TruncInt( x ), 0, n-1 );
         int j1 = pcl::Min( n-1, j0+1 );
         double a = j1 - x;
         double b = x - j0;
         return a*double( m_y[j0] )
              + b*double( m_y[j1] )
              + ((a*a*a - a)*m_dy2[j0] + (b*b*b - b)*m_dy2[j1])/6;
      }
   }
 
   void Clear() override
   {
      UnidimensionalInterpolation<T>::Clear();
      m_dy2.Clear();
   }
 
   bool IsValid() const override
   {
      return m_dy2;
   }
 
private:
 
           double  m_dy1 = 0;      // 1st derivative of spline at the first data point
           double  m_dyn = 0;      // 1st derivative of spline at the last data point
           DVector m_dy2;          // second derivatives of the interpolating function at x[i]
   mutable int     m_current = -1; // index of the current interpolation segment
};
 
#undef m_x
#undef m_y
 
// ----------------------------------------------------------------------------
 
} // pcl
 
#endif  // __PCL_CubicSplineInterpolation_h
 
// ----------------------------------------------------------------------------
// EOF pcl/CubicSplineInterpolation.h - Released 2024-06-18T15:48:54Z