Bspline.hpp

00001 // This file is part of the imaging2 class library.
00002 //
00003 // University of Innsbruck, Infmath Imaging, 2009.
00004 // http://infmath.uibk.ac.at
00005 //
00006 // All rights reserved.
00007 
00008 
00009 #ifndef BSPLINE_H
00010 #define BSPLINE_H
00011 
00012 #include <core/imaging2.hpp>
00013 
00014 namespace imaging
00015 {
00023   template <class DATA_t>
00024   class Bspline
00025   {
00026     size_t _spline_order;
00027     ublas::vector<DATA_t> _coefficients;
00028     ublas::vector<float_t> _knots;
00029 
00030     size_t determine_interval(float_t t) const
00031     {
00032       size_t i = 0;
00033 
00034       while( true )
00035       {
00036         if(_knots.size() == i)
00037           break;
00038 
00039         if(t < _knots(i))
00040           break;
00041 
00042         i++;
00043       }
00044       i = i - 1;
00045 
00046       if(i + _spline_order - 1 >= _knots.size() - 1 )
00047         i = _knots.size() - _spline_order - 1;
00048 
00049       return i;
00050     }
00051 
00052     bool is_in_spline_support(float_t t) const
00053     {
00054       return ( t >= first_knot() && t <= last_knot() );
00055     }
00056 
00057   public:
00058     Bspline() : _knots(), _coefficients() {}
00059 
00061     Bspline(size_t spline_order, size_t n_coefficients) : _spline_order(spline_order), _coefficients(n_coefficients + 1), _knots(n_coefficients + spline_order + 1) {}
00062     
00064     Bspline(const Bspline & source) : _spline_order(source._spline_order), _knots(source._knots), _coefficients(source._coefficients) {}
00065 
00067     const Bspline & operator=(const Bspline & rhs)
00068     { _spline_order = rhs._spline_order; _knots = rhs._knots; _coefficients = rhs._coefficients; return *this; }
00069 
00071     size_t spline_order() const { return _spline_order; }
00072     
00074     float_t knot(size_t i) const { return _knots(i); }
00075     
00077     float_t first_knot() const { return _knots(_spline_order - 1); }
00078     
00080     float_t last_knot() const { return _knots(n_knots() - _spline_order); }
00081     
00083     void set_knot(size_t i, float_t value)
00084     {
00085 #ifdef DEBUG
00086       if(i >= n_knots())
00087         throw Exception("Knot index out of bounds in Bspline::set_knot().");
00088 #endif
00089 
00090       _knots(i) = value;
00091 
00092       if(i == n_knots() - 1)
00093         _knots(n_knots()) = value;
00094     }
00095 
00097     const DATA_t & coefficient(size_t i) const { return _coefficients[i]; }
00098     
00100     void set_coefficient(size_t i, const DATA_t & value)
00101     {
00102 #ifdef DEBUG
00103       if(i >= n_coefficients())
00104         throw Exception("Coefficient index out of bounds in Bspline::set_coefficient().");
00105 #endif
00106 
00107       _coefficients(i) = value;
00108     }
00109 
00111     void resize(size_t spline_order, size_t n_coefficients)
00112     {
00113       _spline_order = spline_order;
00114       _coefficients.resize(n_coefficients + 1);
00115       _knots.resize(n_coefficients + _spline_order + 1);
00116       _coefficients(_coefficients.size() - 1) = DATA_t(0.0);
00117     }
00118 
00120     size_t n_knots() const { return _knots.size() - 1; }
00121     
00123     size_t n_coefficients() const { return _coefficients.size() - 1; }
00124 
00126     DATA_t operator()(float_t t, DATA_t & first_derivative) const
00127       { DATA_t temp; return operator()(t, first_derivative, temp); }
00128       
00130     DATA_t operator()(float_t t) const
00131       { DATA_t temp_1, temp_2; return operator()(t, temp_1, temp_2); }
00132 
00134     void regular_knots(float_t x_0, float_t x_1)
00135     {
00136       if(n_coefficients() < _spline_order)
00137         throw Exception("Too little coefficients for Bspline::init_regular_knots().");
00138 
00139 
00140       float_t x = x_0;
00141       size_t i;
00142 
00143       for(i = 0; i < _spline_order; ++i)
00144         set_knot(i, x_0);
00145 
00146       for( ; i < n_knots() - _spline_order; ++i)
00147       {
00148         x += (x_1 - x_0) / float_t(n_knots() - 2 * _spline_order + 1);
00149         set_knot(i, x);
00150       }
00151 
00152       for(; i < n_knots(); ++i)
00153         set_knot(i, x_1);
00154     }
00155 
00157     DATA_t operator()(float_t t, DATA_t & first_derivative, DATA_t & second_derivative) const
00158     {
00159       size_t k = _spline_order;
00160 
00161       if ( ! is_in_spline_support(t) )
00162       {
00163         first_derivative = DATA_t(0.0);
00164         second_derivative = DATA_t(0.0);
00165         return DATA_t(0.0);
00166       }
00167 
00168       size_t i = determine_interval(t);
00169 
00170       if(t < _knots(_spline_order - 1))
00171       {
00172         first_derivative = DATA_t(0.0);
00173         second_derivative = DATA_t(0.0);
00174         return DATA_t(0.0);
00175       }
00176 
00177       size_t i_offset = i - k + 1;
00178       ublas::matrix<float_t> basis_splines(k + 1, k);
00179 
00180       typename ublas::matrix<float_t>::iterator1 iter1;
00181       typename ublas::matrix<float_t>::iterator2 iter2;
00182 
00183       for(iter1 = basis_splines.begin1(); iter1 != basis_splines.end1(); ++iter1)
00184         for(iter2 = iter1.begin(); iter2 != iter1.end(); ++iter2)
00185           *iter2 = 0.0;
00186 
00187       basis_splines(i - i_offset, 0) = 1;
00188 
00189       for(size_t j = 1; j < k; ++j)
00190         for(size_t ii = i - j; ii <= i; ++ii)
00191         {
00192           float_t n_1, n_2;
00193 
00194           if(_knots(ii + j) == _knots(ii))
00195             n_1 = 0;
00196           else
00197             n_1 = (t - _knots(ii)) /
00198                   (_knots(ii + j) - _knots(ii));
00199 
00200           if(_knots(ii + j + 1) == _knots(ii + 1))
00201             n_2 = 1;
00202           else
00203             n_2 = (_knots(ii + j + 1) - t) /
00204                   (_knots(ii + j + 1) - _knots(ii + 1));
00205 
00206           basis_splines(ii - i_offset, j) =
00207             basis_splines(ii - i_offset, j - 1) * n_1 +
00208             basis_splines(ii + 1 - i_offset, j - 1) * n_2;
00209         }
00210 
00211 
00212       DATA_t value = DATA_t(0.0);
00213       for (size_t ii = i - k + 1; ii <= i; ++ii)
00214         value += basis_splines(ii - i_offset, k - 1) * _coefficients(ii);
00215 
00216 
00217       ublas::vector<DATA_t> coefficient_differences(k); // corresponds to A^(1) in (13) in de Boor's paper
00218       size_t offset = i - k + 2; // offset between the 0 index of coefficient_differences and its real index.
00219 
00220 
00221       first_derivative = DATA_t(0.0);
00222       for (size_t ii = i - k + 2; ii <= i; ++ii)
00223       {
00224         if(_knots(ii + k - 1) != _knots(ii))
00225           coefficient_differences(ii - offset) = (_coefficients(ii) - _coefficients(ii - 1) ) /
00226                                                  (_knots(ii + k - 1) - _knots(ii) );
00227         else
00228           coefficient_differences(ii - offset) = DATA_t(0.0);
00229         first_derivative += basis_splines(ii - i_offset, k - 2)
00230                             * coefficient_differences(ii - offset);
00231       }
00232       first_derivative *= k - 1;
00233 
00234       second_derivative = DATA_t(0.0);
00235       for (size_t ii = i - k + 3; ii <= i; ++ii)
00236         if(_knots(ii + k - 2) != _knots(ii))
00237           second_derivative +=
00238             basis_splines(ii - i_offset, k - 3) /
00239             (_knots(ii + k - 2) - _knots(ii) )
00240             * (coefficient_differences(ii - offset) - coefficient_differences(ii - 1 - offset));
00241       second_derivative *= (k - 1) * (k - 2);
00242 
00243 
00244       return value;
00245 
00246     }
00247     
00249     bool is_in_basis_spline_support(size_t i, float_t t) const
00250     {
00251       return t >= knot(i) && t <= knot(i + spline_order());
00252     }
00253   };
00254 }
00255 
00256 #endif

Generated on Tue Feb 10 10:01:31 2009 for imaging2 by  doxygen 1.5.5