This is the new pcl/AkimaInterpolation.h header that will be released in the next version of the PCL framework. It implements Akima subspline interpolation, which is used in the new CurvesTransform and ColorSaturation processes and interfaces (see
this thread on the Announcements section.
Since the AkimaInterpolation class is completely implemented in this header, it can be readily incorporated into existing installations and will work without problems.
// ----------------------------------------------------------------------------
// PixInsight Class Library - PCL 01.00.15.0125
// Copyright (c) 2003-2007 Pleiades Software. All Rights Reserved.
// ----------------------------------------------------------------------------
// pcl/AkimaInterpolation.h - Released 2007/08/18 18:53:15 UTC
// ----------------------------------------------------------------------------
// This file is part of the PCL. The PCL is an ISO C++ API for development of
// PixInsight modules. The PCL is provided in the hope that it will be useful,
// but it is distributed "AS IS" WITHOUT ANY WARRANTY; without even the
// warranties of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// ----------------------------------------------------------------------------
#ifndef __PCL_AkimaInterpolation_h
#define __PCL_AkimaInterpolation_h
/// \file pcl/AkimaInterpolation.h
#ifndef __PCL_Defs_h
#include <pcl/Defs.h>
#endif
#ifndef __PCL_Diagnostics_h
#include <pcl/Diagnostics.h>
#endif
#ifndef __PCL_UnidimensionalInterpolation_h
#include <pcl/UnidimensionalInterpolation.h>
#endif
#ifndef __PCL_Math_h
#include <pcl/Math.h>
#endif
#ifndef __PCL_Search_h
#include <pcl/Search.h>
#endif
#ifndef __PCL_StatusMonitor_h
#include <pcl/StatusMonitor.h>
#endif
namespace pcl
{
// ----------------------------------------------------------------------------
#define n this->n
#define xv this->xv
#define yv this->yv
// ----------------------------------------------------------------------------
/*! \class AkimaInterpolation
\brief Akima subspline interpolation algorithm
Reference:
Hiroshi Akima, A new method of interpolation and smooth curve fitting
based on local procedures, Journal of the ACM, Vol. 17, No. 4, October 1970,
pages 589-602.
Our implementation is based on the book Numerical Algorithms with C, by
G. Engeln-Müllges and F. Uhlig (Springer, 1996), section 13.1. We properly
represent corners when a data point lies between two adjacent straight lines
with different slopes. This means that our implementation does not impose
continuous differentiability, which deviates from the original work by
Akima. Supporting the accurate representation of corners has several
practical advantages in our opinion; one of them is the enhanced flexibility
for the application of Akima interpolation to graphical representations of
curves given by a set of prescribed x,y data points.
*/
template <class T>
class AkimaInterpolation : public UnidimensionalInterpolation<T>
{
public:
/*!
Constructs an %AkimaInterpolation object.
*/
AkimaInterpolation() : UnidimensionalInterpolation<T>(), b( 0 ), c( 0 ), d( 0 )
{
}
/*!
Destroys an %AkimaInterpolation object.
*/
virtual ~AkimaInterpolation()
{
}
/*!
Initializes a new interpolation.
\param xData Vector of x-values:
<ul>
<li>If xData != 0: Must be a set of monotonically increasing,
distinct values: x[0] < x[1] < ... < x[n-1].</li>
<li>If xData == 0: The interpolation will use implicit
x[i] = i for i = {0,1,2,3,...,n-1}.</li>
</ul>
\param yData Vector of function values for i = {0,1,2,3,...,n-1}.
\param nelem Number of data elements. Must be >= 5. This is because
Akima interpolation requires at least 4 subintervals.
\param monitor Optional address of a status monitoring object. Derived
classes may use this object to report status
information during initialization.
*/
virtual void Initialize( const T* xData, const T* yData, size_type nelem, StatusMonitor* monitor = 0 )
{
PCL_PRECONDITION( nelem >= 5 )
if ( nelem < 5 )
throw Error( "Need five or more data items for AkimaInterpolation" );
Clear();
UnidimensionalInterpolation<T>::Initialize( xData, yData, nelem, monitor );
size_type N = n-1; // Number of subintervals
T* m0 = 0; // Chordal slopes
T* tL = 0; // Akima left-hand slopes to support corners
try
{
b = new T[ N ];
c = new T[ N ];
d = new T[ N ];
m0 = new T[ N + 4 ]; // Room for 4 additional prescribed slopes
T* m = m0 + 2; // Allow negative subscripts m[-1] and m[-2]
tL = new T[ n ];
// Calculate chordal slopes for each subinterval
for ( size_type i = 0; i < N; ++i )
m[i] = (yv[i+1] - yv[i])/(xv[i+1] - xv[i]);
// Prescribed slopes at ending locations
m[-2 ] = 3*m[ 0] - 2*m[ 1];
m[-1 ] = 2*m[ 0] - m[ 1];
m[ N ] = 2*m[N-1] - m[N-2];
m[N+1] = 3*m[N-1] - 2*m[N-2];
// Akima left-hand and right-hand slopes
// Right-hand slopes are just the interpolation coefficients bi.
for ( int i = 0; size_type( i ) < n; ++i )
{
T f = Abs( m[i-1] - m[i-2] );
T e = Abs( m[i+1] - m[i] ) + f;
if ( 1 + e != 1 )
{
tL[i] = m[i-1] + f*(m[i] - m[i-1])/e;
if ( size_type( i ) < N )
b[i] = tL[i];
}
else
{
tL[i] = m[i-1];
if ( size_type( i ) < N )
b[i] = m[i];
}
if ( monitor != 0 )
++*monitor;
}
// Interpolation coefficients b[i], c[i], d[i]. a[i] coefficients are
// the yv[i] ordinate values.
for ( size_type i = 0; i < N; ++i )
{
c[i] = 3*m[i] - 2*b[i] - tL[i+1];
d[i] = b[i] + tL[i+1] - 2*m[i];
if ( xv != 0 )
{
T h = xv[i+1] - xv[i];
c[i] /= h;
d[i] /= h*h;
}
}
delete [] tL, tL = 0;
delete [] m0, m0 = 0;
}
catch ( ... )
{
if ( tL != 0 )
delete [] tL;
if ( m0 != 0 )
delete [] m0;
Clear();
throw;
}
}
/*!
Returns an interpolated function value at \a x location.
*/
virtual double operator()( double x ) const
{
PCL_PRECONDITION( b != 0 && c != 0 && d != 0 )
// Find the subinterval i such that xv[i] <= x < xv[i+1]
// Find the distance dx = x - xv[i], or dx = x - i for implicit x = {0,1,2,3,...n-1}
size_type i;
T dx;
if ( xv != 0 )
{
if ( x <= *xv )
return *yv;
if ( x >= xv[n-1] )
return yv[n-1];
i = InsertionPoint( xv, xv+n, x ) - xv - 1; // -1 because i is a *subinterval* index
dx = x - xv[i];
}
else
{
if ( x <= 0 )
return *yv;
if ( x >= n-1 )
return yv[n-1];
i = size_type( x );
dx = x - i;
}
// Use a Horner scheme to calculate b[i]*dx + c[i]*dx^2 + d[i]*dx^3.
return yv[i] + dx*(b[i] + dx*(c[i] + dx*d[i]));
}
/*!
Frees internal data structures in this AkimaInterpolation object.
*/
virtual void Clear()
{
if ( b != 0 )
delete [] b, b = 0;
if ( c != 0 )
delete [] c, c = 0;
if ( d != 0 )
delete [] d, d = 0;
UnidimensionalInterpolation<T>::Clear();
}
// -------------------------------------------------------------------------
protected:
// Interpolating coefficients for each subinterval.
// The coefficients for dx^0 are the input ordinate values in the yv vector.
T* b; // coefficients for dx^1
T* c; // coefficients for dx^2
T* d; // coefficients for dx^3
};
// ----------------------------------------------------------------------------
#undef n
#undef xv
#undef yv
// ----------------------------------------------------------------------------
} // pcl
#endif // __PCL_AkimaInterpolation_h
// ----------------------------------------------------------------------------
// EOF pcl/AkimaInterpolation.h - Released 2007/08/18 18:53:15 UTC
