Bspline.hpp

// This file is part of the imaging2 class library.
//
// University of Innsbruck, Infmath Imaging, 2009.
// http://infmath.uibk.ac.at
//
// All rights reserved.


#ifndef BSPLINE_H
#define BSPLINE_H

#include <core/imaging2.hpp>

namespace imaging
{
  template <class DATA_t>
  class Bspline
  {
    size_t _spline_order;
    ublas::vector<DATA_t> _coefficients;
    ublas::vector<float_t> _knots;

    size_t determine_interval(float_t t) const
    {
      size_t i = 0;

      while( true )
      {
        if(_knots.size() == i)
          break;

        if(t < _knots(i))
          break;

        i++;
      }
      i = i - 1;

      if(i + _spline_order - 1 >= _knots.size() - 1 )
        i = _knots.size() - _spline_order - 1;

      return i;
    }

    bool is_in_spline_support(float_t t) const
    {
      return ( t >= first_knot() && t <= last_knot() );
    }

  public:
    Bspline() : _knots(), _coefficients() {}

    Bspline(size_t spline_order, size_t n_coefficients) : _spline_order(spline_order), _coefficients(n_coefficients + 1), _knots(n_coefficients + spline_order + 1) {}
    
    Bspline(const Bspline & source) : _spline_order(source._spline_order), _knots(source._knots), _coefficients(source._coefficients) {}

    const Bspline & operator=(const Bspline & rhs)
    { _spline_order = rhs._spline_order; _knots = rhs._knots; _coefficients = rhs._coefficients; return *this; }

    size_t spline_order() const { return _spline_order; }
    
    float_t knot(size_t i) const { return _knots(i); }
    
    float_t first_knot() const { return _knots(_spline_order - 1); }
    
    float_t last_knot() const { return _knots(n_knots() - _spline_order); }
    
    void set_knot(size_t i, float_t value)
    {
#ifdef DEBUG
      if(i >= n_knots())
        throw Exception("Knot index out of bounds in Bspline::set_knot().");
#endif

      _knots(i) = value;

      if(i == n_knots() - 1)
        _knots(n_knots()) = value;
    }

    const DATA_t & coefficient(size_t i) const { return _coefficients[i]; }
    
    void set_coefficient(size_t i, const DATA_t & value)
    {
#ifdef DEBUG
      if(i >= n_coefficients())
        throw Exception("Coefficient index out of bounds in Bspline::set_coefficient().");
#endif

      _coefficients(i) = value;
    }

    void resize(size_t spline_order, size_t n_coefficients)
    {
      _spline_order = spline_order;
      _coefficients.resize(n_coefficients + 1);
      _knots.resize(n_coefficients + _spline_order + 1);
      _coefficients(_coefficients.size() - 1) = DATA_t(0.0);
    }

    size_t n_knots() const { return _knots.size() - 1; }
    
    size_t n_coefficients() const { return _coefficients.size() - 1; }

    DATA_t operator()(float_t t, DATA_t & first_derivative) const
      { DATA_t temp; return operator()(t, first_derivative, temp); }
      
    DATA_t operator()(float_t t) const
      { DATA_t temp_1, temp_2; return operator()(t, temp_1, temp_2); }

    void regular_knots(float_t x_0, float_t x_1)
    {
      if(n_coefficients() < _spline_order)
        throw Exception("Too little coefficients for Bspline::init_regular_knots().");


      float_t x = x_0;
      size_t i;

      for(i = 0; i < _spline_order; ++i)
        set_knot(i, x_0);

      for( ; i < n_knots() - _spline_order; ++i)
      {
        x += (x_1 - x_0) / float_t(n_knots() - 2 * _spline_order + 1);
        set_knot(i, x);
      }

      for(; i < n_knots(); ++i)
        set_knot(i, x_1);
    }

    DATA_t operator()(float_t t, DATA_t & first_derivative, DATA_t & second_derivative) const
    {
      size_t k = _spline_order;

      if ( ! is_in_spline_support(t) )
      {
        first_derivative = DATA_t(0.0);
        second_derivative = DATA_t(0.0);
        return DATA_t(0.0);
      }

      size_t i = determine_interval(t);

      if(t < _knots(_spline_order - 1))
      {
        first_derivative = DATA_t(0.0);
        second_derivative = DATA_t(0.0);
        return DATA_t(0.0);
      }

      size_t i_offset = i - k + 1;
      ublas::matrix<float_t> basis_splines(k + 1, k);

      typename ublas::matrix<float_t>::iterator1 iter1;
      typename ublas::matrix<float_t>::iterator2 iter2;

      for(iter1 = basis_splines.begin1(); iter1 != basis_splines.end1(); ++iter1)
        for(iter2 = iter1.begin(); iter2 != iter1.end(); ++iter2)
          *iter2 = 0.0;

      basis_splines(i - i_offset, 0) = 1;

      for(size_t j = 1; j < k; ++j)
        for(size_t ii = i - j; ii <= i; ++ii)
        {
          float_t n_1, n_2;

          if(_knots(ii + j) == _knots(ii))
            n_1 = 0;
          else
            n_1 = (t - _knots(ii)) /
                  (_knots(ii + j) - _knots(ii));

          if(_knots(ii + j + 1) == _knots(ii + 1))
            n_2 = 1;
          else
            n_2 = (_knots(ii + j + 1) - t) /
                  (_knots(ii + j + 1) - _knots(ii + 1));

          basis_splines(ii - i_offset, j) =
            basis_splines(ii - i_offset, j - 1) * n_1 +
            basis_splines(ii + 1 - i_offset, j - 1) * n_2;
        }


      DATA_t value = DATA_t(0.0);
      for (size_t ii = i - k + 1; ii <= i; ++ii)
        value += basis_splines(ii - i_offset, k - 1) * _coefficients(ii);


      ublas::vector<DATA_t> coefficient_differences(k); // corresponds to A^(1) in (13) in de Boor's paper
      size_t offset = i - k + 2; // offset between the 0 index of coefficient_differences and its real index.


      first_derivative = DATA_t(0.0);
      for (size_t ii = i - k + 2; ii <= i; ++ii)
      {
        if(_knots(ii + k - 1) != _knots(ii))
          coefficient_differences(ii - offset) = (_coefficients(ii) - _coefficients(ii - 1) ) /
                                                 (_knots(ii + k - 1) - _knots(ii) );
        else
          coefficient_differences(ii - offset) = DATA_t(0.0);
        first_derivative += basis_splines(ii - i_offset, k - 2)
                            * coefficient_differences(ii - offset);
      }
      first_derivative *= k - 1;

      second_derivative = DATA_t(0.0);
      for (size_t ii = i - k + 3; ii <= i; ++ii)
        if(_knots(ii + k - 2) != _knots(ii))
          second_derivative +=
            basis_splines(ii - i_offset, k - 3) /
            (_knots(ii + k - 2) - _knots(ii) )
            * (coefficient_differences(ii - offset) - coefficient_differences(ii - 1 - offset));
      second_derivative *= (k - 1) * (k - 2);


      return value;

    }
    
    bool is_in_basis_spline_support(size_t i, float_t t) const
    {
      return t >= knot(i) && t <= knot(i + spline_order());
    }
  };
}

#endif

---