vspline-dev 0.3.1-1 / usr / include / vspline / basis.h

/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2015 - 2017 by Kay F. Jahnke                    */
/*                                                                      */
/*    The git repository for this software is at                        */
/*                                                                      */
/*    https://bitbucket.org/kfj/vspline                                 */
/*                                                                      */
/*    Please direct questions, bug reports, and contributions to        */
/*                                                                      */
/*    kfjahnke+vspline@gmail.com                                        */
/*                                                                      */
/*    Permission is hereby granted, free of charge, to any person       */
/*    obtaining a copy of this software and associated documentation    */
/*    files (the "Software"), to deal in the Software without           */
/*    restriction, including without limitation the rights to use,      */
/*    copy, modify, merge, publish, distribute, sublicense, and/or      */
/*    sell copies of the Software, and to permit persons to whom the    */
/*    Software is furnished to do so, subject to the following          */
/*    conditions:                                                       */
/*                                                                      */
/*    The above copyright notice and this permission notice shall be    */
/*    included in all copies or substantial portions of the             */
/*    Software.                                                         */
/*                                                                      */
/*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
/*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
/*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
/*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
/*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
/*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
/*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
/*    OTHER DEALINGS IN THE SOFTWARE.                                   */
/*                                                                      */
/************************************************************************/

/*! \file basis.h

    \brief Code to calculate the value B-spline basis function
    and it's derivatives.

    There are several variants in here. First, there is a perfectly general
    routine, using the Cox-de Boor recursion. While this is 'nice to have',
    vspline does not actually use it (except as a reference in unit testing).

    vspline only needs evaluation of the B-spline basis function at multiples
    of 0.5. With these values it can construct it's evaluators which in turn
    are capable of evaluating the spline at real coordinates.
    
    So next is a specialized routine using an adapted version of the recursion
    to calculate the basis function's value for integral operands. This isn't
    used in vspline either - instead vspline uses a third version which abbreviates
    the recursion by relying on precomputed values for the basis function with
    derivative 0, which the recursion reaches after as many levels as the
    requested derivative, so seldom deeper than 2. That makes it very fast.

    For comparison there is also a routine calculating an approximation of the
    basis function's value (only derivative 0) by means of a gaussian. This
    routine isn't currently used in vspline.

    for a discussion of the b-spline basis function, have a look at
    http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-basis.html
*/

#ifndef VSPLINE_BASIS_H
#define VSPLINE_BASIS_H

// poles.h has precomputed basis function values sampled at n * 1/2.
// These values were calculated with the general routine (gen_bspline_basis)
// using long doubles for maximum precision. If desired, the process can be
// repeated by running prefilter_poles.cc, which first calculates the basis
// function values and then the prefilter polesbased on the basis function
// values.

#include <vspline/poles.h>

namespace vspline {

/// Implementation of the Cox-de Boor recursion formula to calculate
/// the value of the bspline basis function. This code is taken from vigra
/// but modified to take the spline degree as a parameter.
///
/// This code is quite expensive for higer spline orders
/// because the routine calls itself twice recursively, so the performance is
/// N*N with the spline's degree. Luckily there are ways around using this routine
/// at all - whenever we need the b-spline basis function value in vspline, it is
/// at multiples of 1/2, and poles.h has precomputed values for all spline
/// degrees covered by vspline. I leave the code in here for reference purposes,
/// and also for bootstrapping when precalculating the basis function values
/// and prefilter poles from scratch, like in prefiler_poles.cc, but even there
/// I now use an alternative. Still gen_bspline_basis yields the value of the
/// b-spline basis function at arbitrary x, so it should not miss in a b-spline
/// library.

template < class real_type >
real_type gen_bspline_basis ( real_type x , int degree , int derivative )
{
  if ( degree == 0 )
  {
    if ( derivative == 0 )
        return ( x < real_type(0.5) && real_type(-0.5) <= x )
               ? real_type(1.0)
               : real_type(0.0) ;
    else
        return real_type(0.0);
  }
  if ( derivative == 0 )
  {
    real_type n12 = real_type((degree + 1.0) / 2.0);
    return (     ( n12 + x )
                * gen_bspline_basis<real_type> ( x + real_type(0.5) , degree - 1 , 0 )
              +   ( n12 - x )
                * gen_bspline_basis<real_type> ( x - real_type(0.5) , degree - 1 , 0 )
            )
            / degree;
  }
  else
  {
    --derivative;
    return   gen_bspline_basis<real_type> ( x + real_type(0.5) , degree - 1 , derivative )
           - gen_bspline_basis<real_type> ( x - real_type(0.5) , degree - 1 , derivative ) ;
  }
}

/// this routine is a helper routine to cdb_bspline_basis (below), the
/// modified Cox-de Boor recursion formula to calculate the b-spline basis function
/// for integral operands, operating in int as long as possible. This is achieved by
/// working with 'x2', the doubled x value. Since in the 'real' recursion, the next
/// iteration is called with x +/- 1/2, we can call the 'doubled' version with x +/- 1.
/// This routine recurses 'all the way down to degree 0, So the result is, disregarding
/// arithmetic errors, the same as the result obtained with the general routine.
/// This routine is used in prefiler_poles.cc to obtain the b-spline basis function
/// values at half unit steps, the result is identical to using gen_bspline_basis
/// with long float arguments.

template < class real_type >
real_type cdb_bspline_basis_2 ( int x2 , int degree , int derivative )
{
  if ( degree == 0 )
  {
    if ( derivative == 0 )
        return ( x2 < 1 && -1 <= x2 )
               ? real_type(1.0)
               : real_type(0.0) ;
    else
        return real_type(0.0);
  }
  if ( derivative == 0 )
  {
    int n122 = degree + 1 ;
    return (     ( n122 + x2 )
                * cdb_bspline_basis_2<real_type> ( x2 + 1 , degree - 1 , 0 )
              +   ( n122 - x2 )
                * cdb_bspline_basis_2<real_type> ( x2 - 1 , degree - 1 , 0 )
            )
            / ( 2 * degree ) ;
  }
  else
  {
    --derivative;
    return   cdb_bspline_basis_2<real_type> ( x2 + 1 , degree - 1 , derivative )
           - cdb_bspline_basis_2<real_type> ( x2 - 1 , degree - 1 , derivative ) ;
  }
}

/// modified Cox-de Boor recursion formula to calculate the b-spline basis function
/// for integral operands, delegates to the 'doubled' routine above

template < class real_type >
real_type cdb_bspline_basis ( int x , int degree , int derivative = 0 )
{
  return cdb_bspline_basis_2<real_type> ( x + x , degree , derivative ) ;
}

/// see bspline_basis() below!
/// this helper routine works with the doubled value of x, so it can serve for calls
/// equivalent to basis ( x + .5 ) or basis ( x - .5 ) as basis2 ( 2 * x + 1 ) and
/// basis2 ( 2 * x - 1 ). Having precalculated the basis function at .5 steps, we can
/// therefore avoid using the general recursion formula. This is a big time-saver
/// for high degrees. Note, though, that calculating the basis function for a
/// spline's derivatives still needs recursion, with two branches per level.
/// So calculating the basis function's value for high derivatives still consumes
/// a fair amount of time.

template < class real_type >
real_type bspline_basis_2 ( int x2 , int degree , int derivative )
{
  if ( degree == 0 )
  {
    if ( derivative == 0 )
        return ( x2 < 1 && -1 <= x2 )
               ? real_type(1.0)
               : real_type(0.0) ;
    else
        return real_type(0.0);
  }
  if ( derivative == 0 )
  {
    if ( abs ( x2 ) > degree )
      return real_type ( 0 ) ;
    // for derivative 0 we have precomputed values:
    const long double * pk
      = vspline_constants::precomputed_basis_function_values [ degree ] ;
    return pk [ abs ( x2 ) ] ;
  }
  else
  {
    --derivative;
    return   bspline_basis_2<real_type> ( x2 + 1 , degree - 1 , derivative )
           - bspline_basis_2<real_type> ( x2 - 1 , degree - 1 , derivative ) ;
  }
}

/// bspline_basis produces the value of the b-spline basis function for
/// integral operands, the given degree 'degree' and the desired derivative.
/// It turns out that this is all we ever need inside vspline, the calculation
/// of the basis function at arbitrary points is performed via the matrix
/// multiplication in the weight generating functor, and this functor sets
/// it's internal matrix up with bspline basis function values at integral
/// locations.
///
/// bspline_basis delegates to bspline_basis_2 above, which picks precomputed
/// values as soon as derivative becomes 0. This abbreviates the recursion
/// a lot, since usually the derivative requested is 0 or a small integer.
/// all internal calculations in vspline accessing b-spline basis function
/// values are currently using this routine, not the general routine.

template < class real_type >
real_type bspline_basis ( int x , int degree , int derivative = 0 )
{
  return bspline_basis_2<real_type> ( x + x , degree , derivative ) ;
}

/// Gaussian approximation to B-spline basis function. This routine
/// approximates the basis function of degree spline_degree for real x.
/// I checked for all degrees up to 20. The partition of unity quality of the
/// resulting reconstruction filter is okay for larger degrees, the cumulated
/// error over the covered interval is quite low. Still, as the basis function
/// is never actually evaluated in vspline (whenever it's needed, it is needed
/// at n * 1/2 and we have precomputed values for that) there is not much point
/// in having this function around. I leave the code in for now.

template < typename real_type >
real_type gaussian_bspline_basis_approximation ( real_type x , int degree )
{
  real_type sigma = ( degree + 1 ) / 12.0 ;
  return   real_type(1.0)
         / sqrt ( real_type(2.0 * M_PI) * sigma )
         * exp ( - ( x * x ) / ( real_type(2.0) * sigma ) ) ;
}

} ; // end of namespace vspline

#endif // #define VSPLINE_BASIS_H