CatmullSplineIpol.hpp

 
//
// $Id: CatmullSplineIpol.hpp,v 1.1 2006/07/07 14:10:53 werner Exp $
//
// $Log: CatmullSplineIpol.hpp,v $
// Revision 1.1  2006/07/07 14:10:53  werner
// Intermediate state: multidimensional interpolation with on-demand
// fractional loading appears to work, but still problem with vectorfield
// interpolation.
//
// Revision 1.5  2004/07/09 20:15:26  werner
// Added local alignment capability when interpolating vector fields and
// option to compute the interpolated derivative.
//
// Revision 1.4  2004/07/04 17:24:03  werner
// Recursive multidimensional Interpolation interface implemented.
//
// Revision 1.3  2004/06/18 23:29:46  werner
// Intermediate version that supports the recursive interpolation scheme.
// Need to beautify the interface, fix some indexing bug, and make it to
// work with the ViewPtr.
//
// Revision 1.2  2004/06/15 22:45:31  werner
// Enhanced the fixed array member functions.
// Using better name for interpolation template parameters.
//
// Revision 1.1  2004/06/14 22:42:00  werner
// Moved interpolation files from light++ to VecAl and Multindex from Fiber to VecAl.
//
// Revision 1.1  2002/07/11 16:11:52  werner
// Support for generic spline interpolation and camera paths added.
//
// Created 2002/07/11 17:28:21  bzfbenge
// Initial version
//
#ifndef __CATMULLROMSPLINE_HPP
#define __CATMULLROMSPLINE_HPP "Created 11.07.2002 17:28:21 by bzfbenge"
 
#include "IpolDelimiter.hpp"
#include <assert.h>
 
namespace Fiber
{
 
template <class T>
class   CatMullRomSpline
{
 
double MatrixTerm(int i, double t, double t2, double t3)
{
static  const double Matrix [4][4] =
                  { { -1.0,  3.0, -3.0,  1.0 },
                                        {  2.0, -5.0,  4.0, -1.0 },
                                        { -1.0,    0,  1.0,    0 },
                                        {    0,  2.0,    0,    0 } };
 
                return  Matrix[0][3-i] * t3 +
                Matrix[1][3-i] * t2 +
                Matrix[2][3-i] * t  +
                Matrix[3][3-i];
}
 
double dMatrixTerm(int i, double t, double t2)
{
static  const double Matrix [4][4] =
                  { { -1.0,  3.0, -3.0,  1.0 },
                                        {  2.0, -5.0,  4.0, -1.0 },
                                        { -1.0,    0,  1.0,    0 },
                                        {    0,  2.0,    0,    0 } };
 
                return  Matrix[0][3-i] * 3*t2 +
                Matrix[1][3-i] * 2*t +
                Matrix[2][3-i];
}
 
public:
 
template <class Storage1D, class Limiter>
static  T interpolate(const Storage1D&Data, double t, index_t size, const Limiter&L)
{
    /* there is something to do */
                // TODO: Check for Samplings.size() < 2 !
                assert(size > 2);
 
int lastPoint = size - 1; // point numbering starts with 0
 
/*
    CatmullRomSpline cannot interpolate the first and the last point,
    so LinearInterpolation is used there
 
    (mn)    using the first/last point twice would also help and so
            no LinearInterpolation is needed
*/
                if ( t>=0 && t<1 )
                {
                T       InterPol = Data[1] - Data[0];
                        return Data[0] +  t * InterPol;
                }
 
                if ( t>lastPoint-1  && t <= lastPoint )
                {
                T       InterPol = Data[ lastPoint ] - Data[ lastPoint - 1 ];
                        return Data[ lastPoint - 1 ] +  ( t - ( lastPoint - 1 ) ) * InterPol;
                }
 
index_t index = index_t(t) + 2;
                t = t - floor(t);
 
T       sum;
 
double  t2 = t*t,
                t3 = t2*t;
 
static  const double Matrix [4][4] =
                  { { -1.0,  3.0, -3.0,  1.0 },
                                        {  2.0, -5.0,  4.0, -1.0 },
                                        { -1.0,    0,  1.0,    0 },
                                        {    0,  2.0,    0,    0 } };
 
        for(int i=0; i<4; i++)
        {
        double MatrixTerm = Matrix[0][3-i] * t3 +
                            Matrix[1][3-i] * t2 +
                            Matrix[2][3-i] * t  +
                            Matrix[3][3-i];
 
                if (i==0)
                        sum  = MatrixTerm * Data[ index-i ];
                else
                        sum += MatrixTerm * Data[ index-i ];
        }
 
        return sum * 0.5;
}
};  // class CatMullRomSpline
 
} /* namespace VecAl */
 
#endif /* __CATMULLROMSPLINE_HPP */