CubicIpol.hpp

 
//
// $Id: CubicIpol.hpp,v 1.3 2007/02/22 19:55:34 werner Exp $
//
// $Log: CubicIpol.hpp,v $
// Revision 1.3  2007/02/22 19:55:34  werner
// adjusted fixed array interface
//
// Revision 1.2  2006/07/30 20:43:42  werner
// Partial Derivative of multidimensional arrays enabled.
//
// 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.9  2004/08/17 16:59:52  werner
// Towards compilation with gcc 3.3.3
//
// Revision 1.8  2004/08/12 16:08:45  werner
// Various improvements for implementing the tangential space
//
// Revision 1.7  2004/07/14 15:45:14  werner
// Outsourcing of quadratic matrices into its own header file.
// Added great gaus decomposition from Marcus Weber that allows
// multiple right-hand sides.
//
// Revision 1.6  2004/07/12 10:29:24  werner
// Computing Lower Christoffels in demo application.
//
// 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.
//
#ifndef __IPOL_CUBIC_HPP
#define __IPOL_CUBIC_HPP "Created 11.06.2004 22:28:21 by bzfbenge"
 
#include "Index.hpp"
#include <eagle/vecmath.hpp>
#include <assert.h>
 
namespace Fiber
{
 
template <class T>
class   CubicIpol
{
public:
        typedef T value_type;
 
        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);
 
        //      printf("CubicIpol[%lg]\n", t);
 
                if (t<0     ) return L.limit(Data[0]     );
                if (t>=size ) return L.limit(Data[size-1]);
 
        const double derivation_factor = 2;
 
        index_t i = index_t(t);
        double  x = t - floor(t);
        double  x2 = x*x,
                x3 = x2*x,
 
                // Hermite'sche Polynome
                h00 =  2*x3 - 3*x2 + 1,
//              h10 =  1 - h00,
                h01 = x3 - 2*x2 + x,
                h11 = x3 - x2;
 
                // Boundary conditions.
        index_t i_1 = i-1, i1 = i+1, i2 = i+2;
                if (i_1<0 || i_1>=size) i_1 = 0;
                if (i < 0   )           i  = 0;
                if (i >=size)           i  = size-1;
                if (i1<0)               i1 = 0;
                if (i2<0)               i2 = 0;
                if (i1>=size)           i1 = size-1;
                if (i2>=size)           i2 = size-1;
 
        //      printf("CubicIpol[%d,%d,%d,%d]\n", i_1, i, i1, i2);
 
/*
                // Estimate the Derivatives at A and B
        T       d0 = ( Data[i+1]  - Data[i-1] ) * (1./derivation_factor);
                d1 = ( Data[i+2] -  Data[i  ] ) * (1./derivation_factor);
 
                //         < linear term >             < cubic term >
                return Data[i]*h00 + Data[i+1]*h10 + d0 * h01 + d1 * h11;
 
        Reordering to reduce the number of data array evaluations yields:
*/
 
//              ( (L.limit( Data[i  ] ))*(h00 - (1./derivation_factor) * h11 )
//              + (L.limit( Data[i1 ] ))*(h10 + (1./derivation_factor) * h01 )
//              - (L.limit( Data[i_1] ))*       (1./derivation_factor) * h01
//              + (L.limit( Data[i2 ] ))*       (1./derivation_factor) * h11 );
 
// this version also works for Eagle::PhysicalSpace::point
                return T(          L.limit(Data[i1])
                                + (L.limit(Data[i])-L.limit(Data[i1]))*h00
                                + (L.limit(Data[i1])-L.limit(Data[i_1]))*h01*(1./derivation_factor)
                                + (L.limit(Data[i2])-L.limit(Data[i]))*h11*(1./derivation_factor));
        }
 
        template <class Storage1D, class Limiter>
static  T derivative(const Storage1D&Data, double t, index_t size, const Limiter&L)
        {
        index_t i = index_t(t);
        double  x = t - floor(t),
                x2 = x*x,
 
                // Ableitungen der Hermite'schen Polynome
                h00 =  6*x*(x - 1) ,
                h10 =    - h00,
                h01 = 3*x2 - 4*x + 1,
                h11 = 3*x2 - 2*x;
 
        //      printf("Deriv CubicIpol[%lg] ->", t);
 
        const double derivation_factor = 2;
 
                // Boundary conditions.
        index_t i_1 = i-1, i1 = i+1, i2 = i+2;
                if (i_1<0 || i_1>=size) i_1 = 0;
                if (i < 0   )           i  = 0;
                if (i >=size)           i  = size-1;
                if (i1<0)               i1 = 0;
                if (i2<0)               i2 = 0;
                if (i1>=size)           i1 = size-1;
                if (i2>=size)           i2 = size-1;
 
        //      printf("Indices: [%d,%d,%d,%d]", i_1, i, i1, i2);
        //      printf("Values:  [%lg %lg %lg %lg]\n", Data[i_1], Data[i], Data[i1], Data[i2]);
 
        T       Di_1 = L.limit(Data[i_1]),
                Di   = L.limit(Data[i ]),
                Di1  = L.limit(Data[i1]),
                Di2  = L.limit(Data[i2]);
 
        //      printf("SValues:  [%lg %lg %lg %lg] %s\n",
        //              double(Di_1), double(Di), double(Di1), double(Di2), typeid(Data).name() );
 
                // Estimate the Derivatives at A and B
        T       d0 = ( Di1 - Di_1 ) * (1./derivation_factor),
                d1 = ( Di2 - Di   ) * (1./derivation_factor);
 
                return Di*h00 + Di1*h10 + d0 * h01 + d1 * h11;
        }
};
 
} /* namespace VecAl */
 
#endif /* __IPOL_CUBIC_HPP */