CatmullSplineIpol.hpp

00001 
00002 //
00003 // $Id: CatmullSplineIpol.hpp,v 1.1 2006/07/07 14:10:53 werner Exp $
00004 //
00005 // $Log: CatmullSplineIpol.hpp,v $
00006 // Revision 1.1  2006/07/07 14:10:53  werner
00007 // Intermediate state: multidimensional interpolation with on-demand
00008 // fractional loading appears to work, but still problem with vectorfield
00009 // interpolation.
00010 //
00011 // Revision 1.5  2004/07/09 20:15:26  werner
00012 // Added local alignment capability when interpolating vector fields and
00013 // option to compute the interpolated derivative.
00014 //
00015 // Revision 1.4  2004/07/04 17:24:03  werner
00016 // Recursive multidimensional Interpolation interface implemented.
00017 //
00018 // Revision 1.3  2004/06/18 23:29:46  werner
00019 // Intermediate version that supports the recursive interpolation scheme.
00020 // Need to beautify the interface, fix some indexing bug, and make it to
00021 // work with the ViewPtr.
00022 //
00023 // Revision 1.2  2004/06/15 22:45:31  werner
00024 // Enhanced the fixed array member functions.
00025 // Using better name for interpolation template parameters.
00026 //
00027 // Revision 1.1  2004/06/14 22:42:00  werner
00028 // Moved interpolation files from light++ to VecAl and Multindex from Fiber to VecAl.
00029 //
00030 // Revision 1.1  2002/07/11 16:11:52  werner
00031 // Support for generic spline interpolation and camera paths added.
00032 //
00033 // Created 2002/07/11 17:28:21  bzfbenge
00034 // Initial version
00035 //
00037 #ifndef __CATMULLROMSPLINE_HPP
00038 #define __CATMULLROMSPLINE_HPP "Created 11.07.2002 17:28:21 by bzfbenge"
00039 
00040 #include "IpolDelimiter.hpp"
00041 #include <assert.h>
00042 
00043 namespace Fiber
00044 {
00045 
00051 template <class T>
00052 class   CatMullRomSpline
00053 {
00054 
00055          double MatrixTerm(int i, double t, double t2, double t3)
00056          {
00057 static  const double Matrix [4][4] =  { { -1.0,  3.0, -3.0,  1.0 },
00058                                         {  2.0, -5.0,  4.0, -1.0 },
00059                                         { -1.0,    0,  1.0,    0 },
00060                                         {    0,  2.0,    0,    0 } };
00061 
00062                 return  Matrix[0][3-i] * t3 +
00063                         Matrix[1][3-i] * t2 +
00064                         Matrix[2][3-i] * t  +
00065                         Matrix[3][3-i]; 
00066          }
00067 
00068          double dMatrixTerm(int i, double t, double t2)
00069          {
00070 static  const double Matrix [4][4] =  { { -1.0,  3.0, -3.0,  1.0 },
00071                                         {  2.0, -5.0,  4.0, -1.0 },
00072                                         { -1.0,    0,  1.0,    0 },
00073                                         {    0,  2.0,    0,    0 } };
00074 
00075                 return  Matrix[0][3-i] * 3*t2 +
00076                         Matrix[1][3-i] * 2*t +
00077                         Matrix[2][3-i]; 
00078          }
00079 
00080 public: 
00081 
00082         template <class Storage1D, class Limiter>
00083 static  T interpolate(const Storage1D&Data, double t, index_t size, const Limiter&L)
00084         {
00085                 // TODO: Check for Samplings.size() < 2 !
00086                 assert(size > 2);
00087 
00088         int lastPoint = size - 1; // sub one because point numbering starts with 0 
00089 
00090                 // Catmull Rom cannot interpolate P0 - P1 and Pm-1 - Pm so we use 
00091                 // lin. interpolation here. 
00092 
00093                 // Do lin. interpol for P0-P1 and Pm-1 - Pm
00094                 if ( t>=0 && t<1 )
00095                 {
00096                 T       InterPol = Data[1] - Data[0]; 
00097                         return Data[0] +  t * InterPol;
00098                 } 
00099 
00100                 if ( t>lastPoint-1  && t <= lastPoint )
00101                 {
00102                 T       InterPol = Data[ lastPoint ] - Data[ lastPoint - 1 ]; 
00103                         return Data[ lastPoint - 1 ] +  ( t - ( lastPoint - 1 ) ) * InterPol;
00104                 }
00105 
00106         index_t index = index_t(t) + 2; 
00107                 t = t - floor(t); 
00108 
00109         T       sum; 
00110 
00111         double  t2 = t*t, 
00112                 t3 = t2*t;
00113 
00114 static  const double Matrix [4][4] =  { { -1.0,  3.0, -3.0,  1.0 },
00115                                         {  2.0, -5.0,  4.0, -1.0 },
00116                                         { -1.0,    0,  1.0,    0 },
00117                                         {    0,  2.0,    0,    0 } };
00118 
00119                for(int i=0; i<4; i++)
00120                {
00121                double MatrixTerm = Matrix[0][3-i] * t3 +
00122                                    Matrix[1][3-i] * t2 +
00123                                    Matrix[2][3-i] * t  +
00124                                    Matrix[3][3-i]; 
00125 
00126                        if (i==0)
00127                                sum  = MatrixTerm * Data[ index-i ]; 
00128                        else
00129                                sum += MatrixTerm * Data[ index-i ]; 
00130                } 
00131 
00132                return sum * 0.5;
00133         }
00134 };
00135 
00136 } /* namespace VecAl */ 
00137 
00138 #endif /* __CATMULLROMSPLINE_HPP */