Interpolation.cpp

#include <iostream>

#include <vector/MultiArray.hpp>

#include <vector/CatmullSplineIpol.hpp>
#include <vector/CubicIpol.hpp>
#include <vector/FastCubicIpol.hpp>
#include <vector/LinearIpol.hpp>
#include <vector/NearestNeighborIpol.hpp>
#include <vector/Interpolate.hpp>

#include <time.h>

using namespace std;
using namespace Fiber;

template class  LinearIpol<double>;
template class  FastCubicIpol<double>;
template class  CubicIpol<double>;
template class  CatMullRomSpline<double>;

//template class MultiIndex<1>;
template class MultiIndex<2>;
template class MultiIndex<3>;
template class MultiIndex<4>;


int     MA2()
{
        enum { NX = 19, NY = 10 }; 

double  Data[ NX*NY ]; 

MultiArray<2, double> D(Data, MIndex(NX, NY) );

        D.set(0.0);

        D[ MIndex(2,5) ] = 9.0; 
        D[ MIndex(12,7) ] = 5.0;

        {for(int y=0; y<D.Size()[1]; y++)
        {
                for(int x=0; x<D.Size()[0]; x++)
                        cout << " " << D[y][x];
                cout << endl;
        }} 

//      cout << "------------------------" << endl; 
//      cout << D; 
//      cout << "------------------------" << endl; 


// Alternative indexing via MultiIndex
MultiIndex<2> I;
        I = 2,5;

double  value = D[ I ];

        cout << "Indexed Value = " << value << endl;

        // Note that indexing by consecuted [] operators is REVERSED
        // as to using a MultiIndex because the LAST index is running
        // fastest! 
        cout << "5th slice is " << D[5] << endl; 

        if (value != D[5][2])
                return 0;

        cout << "Value at 2,5 is " << value << endl;

FixedArray<double,2> point; point = 2.1, 5.25;

        cout << "------------------------" << endl;
Interpolate<2, double, LinearIpol<double>, double> LI(D, point);

double  V2 = LI.eval(); 
        if ( fabs(V2 -6.075) > 1E-10 )
                return 0;

        cout << "Ipol2D: value=" << V2 << endl; 

        for(point[0]=2.0; point[0]<3.0; point[0]+=0.1)
        {
                cout << "Ipol2D[ " << point << "] = " << LI.eval() << endl;
        } 

        cout << "2D Interpolation successfully demonstrated.\n" << endl;
        return 1;
}

static void pit(const Iterator<double>&D)
{
        for(index_t i=0; i<D.count(); i++)
                cout << D[i] << ' '; 
        cout << endl;
}

int     MA3()
{
        enum { NX=13, NY=10, NZ=23 }; 

        cout << "Demonstration of Interpolating within a 3D array" << endl;

double  Data[ NX*NY*NZ ]; 
MultiArray<3, double> D( Data, MIndex(NX,NY,NZ) );

        D.set(0.0);
MultiIndex<3> I;
        I = 5,6,7;

        D[ I ] = 100.0;

double  value = D[ I ];
        assert(value==100.0);

FixedArray<double,3> point;

Interpolate<3, double, LinearIpol<double>      , double> LI(D, point);
Interpolate<3, double, FastCubicIpol<double>   , double> LIfc(D, point);
Interpolate<3, double, CubicIpol<double>       , double> LIc(D, point);
Interpolate<3, double, CatMullRomSpline<double>, double> LIcs(D, point); 

        for(double fr=0; fr<1; fr+=0.1)
        {
                point = 5.1, 6.2, 7+fr;

                cout << "Interpolation @" << point 
                     << "\t Linear: "     << LI.eval() 
                     << "\t FCubic: "     << LIfc.eval()                     
                     << "\t Cubic: "      << LIc.eval()
                     << "\t Catmull: "    << LIcs.eval()
                     <<  endl;
        }

#if 0
{
cerr << "-------------------------------------------------" << endl;    
   clock_t start, finish;
   int N =1000000;
FixedArray<3, double> point; point = 5.1, 6.2, 7.6;   
Interpolate<3, double, NearestNeighborIpol<double>, double> LIn(D, point);
Interpolate<3, double, LinearIpol<double>         , double> LI(D, point);
Interpolate<3, double, FastCubicIpol<double>      , double> LIfc(D, point);
Interpolate<3, double, CubicIpol<double>          , double> LIc(D, point);
Interpolate<3, double, CatMullRomSpline<double>   , double> LIcs(D, point);
/*
Pentium IV, 1.6GHz, Intel Compiler, Optimization:
---------------------------
No     Ipol: 0.1 seconds
Linear Ipol: 1.1 seconds   11x    1x
FCubic Ipol: 1.1 seconds
Cubic  Ipol: 4.0 seconds   40x    3.9x
CSplineIpol: 4.6 seconds   46x    4.05x

  
Pentium IV, 1.6GHz, Intel Compiler, Debug mode:
----------------------------
No     Ipol: 2.4 seconds
Linear Ipol: 7.1 seconds      3x
FCubic Ipol: 7.3 seconds
Cubic  Ipol: 24.9 seconds    10x
CSplineIpol: 25.3 seconds

Centrino Duo, 1.66GHz, GCC 3.4.4, Optimize
----------------------------
No     Ipol: 0.1 seconds
Linear Ipol: 0.5 seconds
FCubic Ipol: 0.5 seconds
Cubic  Ipol: 2.3 seconds
CSplineIpol: 3.0 seconds
*/
        start = clock(); {for(int i=0; i<N; i++) LIn.eval();} finish = clock(); 
        fprintf(stderr, "No     Ipol: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC ); 

        start = clock(); {for(int i=0; i<N; i++) LI.eval();} finish = clock(); 
        fprintf(stderr, "Linear Ipol: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );

        start = clock(); {for(int i=0; i<N; i++) LIfc.eval();} finish = clock(); 
        fprintf(stderr, "FCubic Ipol: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );

        start = clock(); {for(int i=0; i<N; i++) LIc.eval();} finish = clock(); 
        fprintf(stderr, "Cubic  Ipol: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC );

        start = clock(); {for(int i=0; i<N; i++) LIcs.eval();} finish = clock(); 
        fprintf(stderr, "CSplineIpol: %2.1f seconds\n", (double)(finish - start) / CLOCKS_PER_SEC ); 
}
#endif

        return 1;
}


void plain()
{
MultiIndex<3> Dims;
        Dims = 10,10,10;

double  *d = new double[Dims.size()];   

MultiIndex<3> I;
        I = 5,5,5;

double  value = d[ I.linear(Dims) ];

        delete d; 
}

int InterpolationExample()
{
int     r=0; 
        r  = MA2();
        r += MA3();
        return r;
}

---