interpolnd.h

Go to the documentation of this file.

#ifndef interpolnd_h
#define interpolnd_h

/*
________________________________________________________________________

 (C) dGB Beheer B.V.; (LICENSE) http://opendtect.org/OpendTect_license.txt
 Author:        Bert
 Date:          Apr 2010
 RCS:           $Id$
________________________________________________________________________

*/

#include "math2.h"


namespace Interpolate
{

template <class T>
inline T linearRegND( int N, const T* v, const T* pos )
{
    if ( N == 0 )
        return v[0];
    else if ( N == 1 )
        return v[0] * (1-pos[0]) + v[1] * pos[0];
    else if ( N == 2 )
    {
        T a[4];
        a[0] = v[0];
        a[1] = v[1] - v[0];
        a[2] = v[2] - v[0];
        a[3] = v[3] + v[0] - v[1] - v[2];
        return a[0] + a[1]*pos[0] + a[2]*pos[1] + a[3]*pos[0]*pos[1];
    }
    else
    {
        const int lowerN = N-1;
        const od_int64 nlowerpts = Math::IntPowerOf( ((od_int64)2), lowerN );
        float* lowerv = new float [nlowerpts];
        const float lastpos = pos[lowerN];
        for ( od_int64 idx=0; idx<nlowerpts; idx++ )
        {
            const float v0 = v[idx];
            const float v1 = v[idx+nlowerpts];
            lowerv[idx] = (1-lastpos) * v0 + lastpos * v1;
        }
        const float res = linearRegND( lowerN, lowerv, pos );
        delete [] lowerv;
        return res;
    }
}


} // namespace Interpolate

#endif