bspline.h

Go to the documentation of this file.

/* This software and supporting documentation are distributed by
 *     Institut Federatif de Recherche 49
 *     CEA/NeuroSpin, Batiment 145,
 *     91191 Gif-sur-Yvette cedex
 *     France
 *
 * This software is governed by the CeCILL-B license under
 * French law and abiding by the rules of distribution of free software.
 * You can  use, modify and/or redistribute the software under the
 * terms of the CeCILL-B license as circulated by CEA, CNRS
 * and INRIA at the following URL "http://www.cecill.info".
 *
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 *
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the
 * same conditions as regards security.
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-B license and that you accept its terms.
 */

#ifndef AIMSALGO_MATH_BSPLINE_H
#define AIMSALGO_MATH_BSPLINE_H

#include <aims/vector/vector.h>         // Point2dd
#include <cartobase/exception/assert.h>
#include <cmath>    // std::abs
#include <limits>   // std::numeric_limits
#include <vector>   // std::vector

namespace aims {

  //==========================================================================
  //
  //     B-SPLINE : ORDER AS TEMPLATE ARGUMENT
  //
  //==========================================================================

  template <unsigned Order>
  struct BSplineOrder
  {
    static double at( double x )
    {
      return spline(x);
    }

    static double spline( double x )
    {
      if( std::abs(x) > double(Order+1)/2 )
        return 0.;

      return ( dO()/2. + .5 - x ) / dO()
             * BSplineOrder<Order-1>::spline( x + .5 )
             + ( dO()/2. + .5 + x  ) / dO()
             * BSplineOrder<Order-1>::spline( x - .5 );
    }

    static double derivative( double x )
    {
      return BSplineOrder<Order-1>::spline( x + .5 )
             - BSplineOrder<Order-1>::spline( x - .5 );
    }

    static double derivative2( double x )
    {
      ASSERT( Order > 1 );
      return BSplineOrder<Order-1>::derivative( x + .5 )
             - BSplineOrder<Order-1>::derivative( x - .5 );
    }

    static unsigned order() { return Order; }

  private:
    static double dO()     { return (double)Order; }
  };

  template <>
  struct BSplineOrder<0>
  {
    static double spline( double value )
    {
      if( -.5 <= value && value < .5 )
        return 1.;
      else
        return 0.;
    }

    static unsigned order() { return 0; }
  };

  //==========================================================================
  //
  //     B-SPLINE : MULTI ORDER
  //
  //==========================================================================

  class BSpline
  {
  public:
    BSpline( unsigned order = 3, float scale = 1., bool shifted = false );
    BSpline( const BSpline & other );
    virtual ~BSpline();
    BSpline & operator= ( const BSpline & other );

    virtual
    double   operator() ( double x )                const;
    virtual
    double   at( double x )                     const;
    virtual
    double   derivative( double x, unsigned n = 1 ) const;
    virtual
    const Point2dd & support() const;

    virtual
    void reset( unsigned order = 3, float scale = 1., bool shifted = false );
    virtual
    void setOrder( unsigned order );
    virtual
    void setScale( float scale );
    virtual
    void setShifted( bool shifted );
    virtual
    void setCentered( bool centered );

    unsigned order()    const { return _order; }
    float    scale()    const { return _scale; }
    bool     shifted()  const { return _shift; }
    bool     centered() const { return !_shift; }

  protected:
    double dO()     const { return double( _order ); }
    virtual void setSupport();
    unsigned _order;
    float    _scale;
    bool     _shift;
    Point2dd _support;
  };

  //==========================================================================
  //
  //     B-SPLINE : FAST
  //
  //==========================================================================

  class FastBSpline: public BSpline
  {
  public:
    FastBSpline( unsigned order = 3, float scale = 1., bool shifted = false );
    FastBSpline( const FastBSpline & other );
    virtual ~FastBSpline();
    FastBSpline & operator= ( const FastBSpline & other );

    virtual
    double   at( double x )                     const;
    virtual
    double   derivative( double x, unsigned n = 1 ) const;

  };

  //==========================================================================
  //
  //     B-SPLINE : TABULATED VALUES
  //
  //==========================================================================

  class TabulBSpline: public BSpline
  {
  public:
    TabulBSpline( unsigned order = 3,
                  unsigned nder = 1,
                  float scale = 1.,
                  bool shifted = false,
                  size_t length = std::numeric_limits<unsigned short>::max() - 1 );
    TabulBSpline( const TabulBSpline & other );
    virtual ~TabulBSpline();
    TabulBSpline & operator= ( const TabulBSpline & other );

    virtual
    double   at( double x )                     const;
    virtual
    double   derivative( double x, unsigned n = 1 ) const;

    virtual
    void reset( unsigned order = 3,
                unsigned nder = 1,
                float scale = 1.,
                bool shifted = false,
                size_t length = std::numeric_limits<unsigned short>::max() - 1 );
    virtual
    void setOrder( unsigned order );
    virtual
    void setNbDer( unsigned nder );
    virtual
    void setLength( size_t length );

    virtual
    unsigned nder()     const;
    virtual
    size_t   length()   const;

  protected:
    size_t   index( double x )          const;
    bool     is_valid( size_t index ) const;
    void     setArray( unsigned nder, size_t length );
    std::vector<std::vector<double> > _values;
    double _tablescale;
    double _tableoffset;
  };

  //==========================================================================
  //
  //     DISCRETE B-SPLINE
  //
  //==========================================================================

  class DiscreteBSpline: public BSpline
  {
  public:
    DiscreteBSpline( unsigned order = 3, float scale = 1., bool shifted = false );
    DiscreteBSpline( const DiscreteBSpline & other );
    virtual ~DiscreteBSpline();
    DiscreteBSpline & operator= ( const DiscreteBSpline & other );

    virtual
    double operator() ( int x ) const;
    virtual
    double at( int x ) const;

    virtual
    void reset( unsigned order = 3, float scale = 1., bool shifted = false );
    virtual
    void setOrder( unsigned order );
    virtual
    void setScale( float scale );
    virtual
    void setShifted( bool shifted = true );
    virtual
    void setCentered( bool centered = true );

  protected:
    void setArray();
    void setSupport();
    std::vector<double> _values;
  };

} // namespace aims

#endif // AIMSALGO_MATH_BSPLINE_H