Spline Class Reference

#include <cubic_spline.h>

Collaboration diagram for Spline:

Public Types
enum	Type { LINEAR, CUBIC, CONSTRAINED_CUBIC }
enum	Extr { ZERO, CONSTANT, QUADRATIC }

Public Member Functions
template<typename Container >
void	set_points (const Container &x, const Container &y)
void	solve_coeffs_linear ()
void	solve_coeffs_cubic ()
void	solve_coeffs_constrained_cubic ()
template<typename Container >
void	constructAndReset (const Container &x, const Container &y)
Float	extrapolate_left (double x) const
Float	extrapolate_right (double x) const
Float	eval (Float x) const
Float	eval (Float x, int idx) const
void	print ()

Public Attributes
int	npoints
Type	splineType = CUBIC
Extr	extrapolate = CONSTANT

Protected Attributes
std::vector< Float >	m_a
std::vector< Float >	m_b
std::vector< Float >	m_c
std::vector< Float >	m_x
std::vector< Float >	m_y

Detailed Description

Definition at line 121 of file cubic_spline.h.

Member Enumeration Documentation

Extr

enum Spline::Extr

Enumerator
ZERO
CONSTANT
QUADRATIC

Definition at line 127 of file cubic_spline.h.

{ZERO, CONSTANT, QUADRATIC};

Type

enum Spline::Type

Enumerator
LINEAR
CUBIC
CONSTRAINED_CUBIC

Definition at line 125 of file cubic_spline.h.

{LINEAR, CUBIC, CONSTRAINED_CUBIC};

Member Function Documentation

constructAndReset()

template<typename Container >

void Spline::constructAndReset ( const Container &  x, const Container &  y  )

inline

Definition at line 246 of file cubic_spline.h.

                                                                  { // construct should be able to take any containers, as the points will be copied to the Spline's own storage anyway
        set_points(x,y);
    }

Here is the call graph for this function:

eval() [1/2]

Float Spline::eval ( Float  x ) const

inline

Definition at line 268 of file cubic_spline.h.

                                     {
        size_t n=m_x.size();

        if(x<m_x[0]) return extrapolate_left(x);
        else if(x>=m_x[n-1]) return extrapolate_right(x);
        else {
            // find the closest point m_x[idx] < x
//          std::vector<double>::const_iterator it;
//          it=std::lower_bound(m_x.begin(),m_x.end(),x);
//          int idx=std::max( int(it-m_x.begin())-1, 0);
            int it = my_lower_bound(x, m_x.data(), m_x.size());
            int idx=std::max(it-1,0); //std::max(it-1, 0);
//          int it_ub = my_upper_bound(x, m_x.data(), m_x.size());
//          int idx_ub=it_ub-1; //std::max(it_ub-1, 0);
//          std::cout << "x: " << x << " " << idx << " " << idx_ub << endl;
            
            return eval(x,idx);         
        }
    }

Here is the call graph for this function:

Here is the caller graph for this function:

eval() [2/2]

Float Spline::eval ( Float  x, int  idx  ) const

inline

Definition at line 288 of file cubic_spline.h.

                                             {
        double h=x-m_x[idx];
        return  ((m_a[idx]*h + m_b[idx])*h + m_c[idx])*h + m_y[idx];
    }

extrapolate_left()

Float Spline::extrapolate_left ( double  x ) const

inline

Definition at line 250 of file cubic_spline.h.

                                                 {
          double h=x-m_x[0];
          // extrapolation to the left
          if (extrapolate == ZERO) return 0;    // zero
          else if (extrapolate == CONSTANT) return m_y[0];  // constant
          else return ((m_b[0])*h + m_c[0])*h + m_y[0];  // quadratic
    }

Here is the caller graph for this function:

extrapolate_right()

Float Spline::extrapolate_right ( double  x ) const

inline

Definition at line 258 of file cubic_spline.h.

                                                  {
            size_t n=m_x.size();
            double h=x-m_x[n-1];
          // extrapolation to the right
            if (extrapolate == ZERO) return 0;  // zero
            else if (extrapolate == CONSTANT) return m_y[n-1]; // constant
            else return ((m_b[n-1])*h + m_c[n-1])*h + m_y[n-1];  // quadratic
    
    }

Here is the caller graph for this function:

print()

void Spline::print ( )

inline

Definition at line 293 of file cubic_spline.h.

                       {
        std::cout << "Spline: \n";
        for (int i=0; i<npoints; ++i){
            std::cout << std::right << std::setw(10) << std::setfill(' ') << m_x[i] << " " 
                      << std::right << std::setw(10) << std::setfill(' ') << m_y[i] << " "
                      << std::right << std::setw(10) << std::setfill(' ') << m_a[i] << " "
                      << std::right << std::setw(10) << std::setfill(' ') << m_b[i] << " "
                      << std::right << std::setw(10) << std::setfill(' ') << m_c[i] << " " 
                      << "\n";
        }
        std:: cout << std::endl;
    }

set_points()

template<typename Container >

void Spline::set_points ( const Container &  x, const Container &  y  )

inline

Definition at line 139 of file cubic_spline.h.

                                                           {    // construct should be able to take any containers, as the points will be copied to the Spline's own storage anyway
       assert(x.size() == y.size());
       m_x.assign(x.begin(), x.end());
       m_y.assign(y.begin(), y.end());
       
        npoints = x.size();
        m_a.resize(npoints);
        m_b.resize(npoints);
        m_c.resize(npoints);

       // Check that x is in increasing order
       for(int i=0; i<npoints-1; i++) {
          assert(m_x[i]<m_x[i+1]);
       }

       if( splineType == CUBIC) { // cubic spline interpolation
            solve_coeffs_cubic();
       } 
       else if (splineType == LINEAR) { // linear interpolation
            solve_coeffs_linear();
       }
       else if (splineType == CONSTRAINED_CUBIC){
            solve_coeffs_constrained_cubic();   
       }

    }

Here is the call graph for this function:

Here is the caller graph for this function:

solve_coeffs_constrained_cubic()

void Spline::solve_coeffs_constrained_cubic ( )

inline

Definition at line 219 of file cubic_spline.h.

                                                {
        int n = npoints;  

        std::vector <double> m(n);
        for (int i=1; i<n-1; ++i){
            double fplus  = (m_y[i+1]-m_y[i])/(m_x[i+1]-m_x[i]);
            double fminus = (m_y[i]-m_y[i-1])/(m_x[i]-m_x[i-1]);
            if (fplus*fminus < 0) m[i] = 0;
            else m[i] = 2/(1/fplus + 1/fminus);  // harmonic mean will tend towards smaller derivative
            //else m[i] = (fplus+fminus)/2; 
        }
        // Boundary conditions: zero curvature at x0 and xN-1 
        m[0] = 3*(m_y[1]-m_y[0])/(m_x[1]-m_x[0])/2 - m[1]/2;
        m[n-1] = 3*(m_y[n-1]-m_y[n-2])/(m_x[n-1]-m_x[n-2])/2 - m[n-2]/2;

        for (int i=0; i<n-1; ++i){
            m_c[i] = m[i];
            double hi = m_x[i+1]-m_x[i];
            m_b[i] = ( 3*(m_y[i+1]-m_y[i] - m[i]*hi) - hi*(m[i+1]-m[i]))/hi/hi;
            m_a[i] = (-2*(m_y[i+1]-m_y[i] - m[i]*hi) + hi*(m[i+1]-m[i]))/hi/hi/hi;
        }
        m_b[n-1] = 0;
        m_c[n-1] = m[n-1];
    }

Here is the caller graph for this function:

solve_coeffs_cubic()

void Spline::solve_coeffs_cubic ( )

inline

Definition at line 179 of file cubic_spline.h.

                                    {
        int n = npoints;  
          // setting up the matrix and right hand side of the equation system
          // for the parameters b[]
          //band_matrix A(n,1,1);
          std::vector <Float> l(n,0), m(n,0), u(n,0);   // lower, middle, and upper diagonals of thr matrix. These hold coefficients of b{i-1}, b{i}, b{i+1} respectively
          std::vector<double>  rhs(n);
          for(int i=1; i<n-1; i++) {
             l[i] =1.0/3.0*(m_x[i]-m_x[i-1]);
             m[i] =2.0/3.0*(m_x[i+1]-m_x[i-1]);
             u[i] =1.0/3.0*(m_x[i+1]-m_x[i]);
             rhs[i]=(m_y[i+1]-m_y[i])/(m_x[i+1]-m_x[i]) - (m_y[i]-m_y[i-1])/(m_x[i]-m_x[i-1]);
          }
          // boundary conditions, zero curvature b[0]=b[n-1]=0
          m[0]     = 2.0;
          u[0]     = 0.0;
          rhs[0]   = 0.0;   
          m[n-1]   = 2.0;
          l[n-1]   = 0.0;
          rhs[n-1] = 0.0;   

          // solve the equation system to obtain the parameters b[]
          thomas_solve(l.data(),m.data(),u.data(),rhs.data(),n);
            m_b = rhs;

          // calculate parameters a[] and c[] based on b[]
          for(int i=0; i<n-1; i++) {
             m_a[i]=1.0/3.0*(m_b[i+1]-m_b[i])/(m_x[i+1]-m_x[i]);
             m_c[i]=(m_y[i+1]-m_y[i])/(m_x[i+1]-m_x[i])
                    - 1.0/3.0*(2.0*m_b[i]+m_b[i+1])*(m_x[i+1]-m_x[i]);
          }
           // for the right boundary we define
           // f_{n-1}(x) = b*(x-x_{n-1})^2 + c*(x-x_{n-1}) + y_{n-1}
           double h=m_x[n-1]-m_x[n-2];
           // m_b[n-1] is determined by the boundary condition (turns out to be 0; = 0 for 0 curvature at xn-1)
           m_a[n-1] = 0.0;
           m_c[n-1]=3.0*m_a[n-2]*h*h+2.0*m_b[n-2]*h+m_c[n-2];   // = f'_{n-2}(x_{n-1})
    }

Here is the call graph for this function:

Here is the caller graph for this function:

solve_coeffs_linear()

void Spline::solve_coeffs_linear ( )

inline

Definition at line 167 of file cubic_spline.h.

                                     {
        int n = npoints;  
        for(int i=0; i<n-1; i++) {
             m_a[i]=0.0;
             m_b[i]=0.0;
             m_c[i]=(m_y[i+1]-m_y[i])/(m_x[i+1]-m_x[i]);
        }
        m_a[n-1] = m_b[n-1] = 0.0;
        m_c[n-1] = m_c[n-2];    // fn-1 is same as fn-2 - same line extends right for extrapolation
    }

Here is the caller graph for this function:

Member Data Documentation

extrapolate

Extr Spline::extrapolate = CONSTANT

Definition at line 128 of file cubic_spline.h.

m_a

std::vector<Float> Spline::m_a

protected

Definition at line 134 of file cubic_spline.h.

m_b

std::vector<Float> Spline::m_b

protected

Definition at line 134 of file cubic_spline.h.

m_c

std::vector<Float> Spline::m_c

protected

Definition at line 134 of file cubic_spline.h.

m_x

std::vector<Float> Spline::m_x

protected

Definition at line 135 of file cubic_spline.h.

m_y

std::vector<Float> Spline::m_y

protected

Definition at line 135 of file cubic_spline.h.

npoints

int Spline::npoints

Definition at line 123 of file cubic_spline.h.

splineType

Type Spline::splineType = CUBIC

Definition at line 126 of file cubic_spline.h.

---

The documentation for this class was generated from the following file:

• include/cubic_spline.h