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)
        {
                // TODO: Check for Samplings.size() < 2 !
                assert(size > 2);

        int lastPoint = size - 1; // sub one because point numbering starts with 0 

                // Catmull Rom cannot interpolate P0 - P1 and Pm-1 - Pm so we use 
                // lin. interpolation here. 

                // Do lin. interpol for P0-P1 and Pm-1 - Pm
                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;
        }
};

} /* namespace VecAl */ 

#endif /* __CATMULLROMSPLINE_HPP */