vspline-dev 0.3.1-1 / usr / include / vspline / bspline.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 bspline.h

    \brief defines class bspline

  class bspline is an object to contain a b-spline's coefficients and some
  metadata in one handy package. It also provides easy access to b-spline
  prefiltering. The idea is that user code establishes a bspline object
  representing the data at hand and then proceeds to create 'evaluators'
  to evaluate the spline. You may be reminded of SciPy's bisplrep object,
  and I admit that SciPy's bspline code has been one of my inspirations.

  It attempts to do 'the right thing' by automatically creating suitable helper
  objects and parametrization so that the spline does what it's supposed to do.
  Most users will not need anything else, and using class bspline is quite
  straightforward. It's quite possible to have a b-spline up and running with
  a few lines of code without even having to make choices concerning it's
  parametrization, since there are sensible defaults for everything. At the same
  time, pretty much everything *can* be parametrized even at this level.
  bspline objects can be used without any knowledge of their internals,
  e.g. as parameters to the remap functions.
  
  Note that class bspline does not provide evaluation of the spline. To evaluate,
  objects of class evaluator (see eval.h) are used, which construct from a bspline
  object with additional parameters, like, whether to calculate the spline's
  value or it's derivative(s) or whether to use optimizations for special cases.
  
  While using 'raw' coefficient arrays with an evaluation scheme which applies
  boundary conditions is feasible and most memory-efficient, it's not so well
  suited for very fast evaluation, since the boundary treatment needs conditionals,
  and 'breaks' uniform access, which is especially detrimental when using
  vectorization. So vspline uses coefficient arrays with a few extra coefficients
  'framing' or 'bracing' the 'core' coefficients. Since evaluation of the spline
  looks at a certain small section of coefficients (the evaluator's 'support'),
  the bracing is chosen so that this lookup will always succeed without having to
  consider boundary conditions: the brace is set up to make the boundary conditions
  explicit, and the evaluation can proceed blindly without bounds checking. With
  large coefficient arrays, the little extra space needed for the brace becomes
  negligible, but the code for evaluation becomes much faster and simpler.
  In effect, 'bracing' is taken a little bit further than merely providing
  enough extra coefficients to cover the support: additional coefficients are
  produced to allow for the spline to be evaluated without bounds checking
  
  - at the lower and upper limit of the spline's defined range
  
  - and even slightly beyond those limits
  
  This makes the code more robust: being very strict about the ends of the
  defined range can easily result in quantization errors producing out-of-bounds
  access to the coefficient array, so the slightly wider brace acts as a safeguard.

  class bspline handles two views to the coefficients it operates on, these
  are realized as vigra::MultiArrayViews, and they share the same storage:

  - the 'core', which is a view to an array of data precisely the same shape as
    the knot point data over which the spline is calculated.

  - 'container', which surrounds the view above with an additional frame of coefficients
    used for 'bracing', for additional coefficints needed for using the 'explicit'
    scheme of extrapolation before prefiltering, and as extra 'headroom' if 'shifting'
    the spline is intended. 

  Using class bspline, there is a choice of 'strategy'. The simplest strategy is
  'UNBRACED'. With this strategy, after putting the knot point data into the bspline's
  'core' area and calling prefilter(), the core area will contain the b-spline coefficients.
  The resulting b-spline object can't be safely evaluated with the code in eval.h.
  this mode of operation is intended for users who want to do their own processing of the
  coefficients and don't need the code in eval.h. prefiltering is done using an implicit
  scheme as far as the boundary conditions are concerned.
  
  The 'standard' strategy is 'BRACED'. Here, after prefiltering, the container
  in the bspline object will contain the b-spline coefficients, surrounded by a 'brace'
  of coefficients which allows code in eval.h to process them without special treatment
  for values near the border (the brace covers what support is needed by marginal
  coefficients). Again, an implicit scheme is used.
  
  The third strategy, 'EXPLICIT', extrapolates the knot point data in the 'core' area
  sufficiently far to suppress margin effects when the prefiltering is performed without
  initial coefficient calculation. If the 'frame' of extrapolated data is large enough,
  the result is just the same. The inner part of the frame is taken as the brace, so no
  bracing needs to be performed explicitly. The resulting b-spline object will work with
  vspline's evaluation code. Note that the explicit scheme uses 'GUESS' boundary conditions
  on the (framed) array, which tries to minimize margin effects further.

  Also note that the additional memory needed for the 'frame' will be held throughout the bspline
  object's life, the only way to 'shrink' the coefficient array to the size of the braced or core
  coefficients is by copying them out to a smaller array.

  The fourth strategy, 'MANUAL', is identical to 'EXPLICIT', except that automatic extrapolation
  of the core data to the frame is not performed. Instead, this strategy relies on the user to
  fill the frame with extrapolated data. This is to allow for the user to apply custom
  extrapolation schemes. The procedure would be to create the bspline object, fill the core,
  apply the extrapolation, then call prefilter.

  Probably the most common scenario is that the source data for the spline are available from
  someplace like a file. Instead of reading the file's contents into memory first and passing
  the memory to class bspline, there is a more efficient way: a bspline object is set up
  first, with the specification of the size of the incoming data and the intended mode of
  operation. The bspline object allocates the memory it will need for the purpose, but
  doesn't do anything else. The 'empty' bspline object is then 'filled' by the user
  by putting data into it's 'core' area. Subsequently, prefilter() is called, which converts
  the data to b-spline coefficients. This way, only one block of memory is used throughout,
  the initial data are overwritten by the coefficients, operation is in-place and most efficient.

  If this pattern can't be followed, there are alternatives:

  - if a view to an array at least the size of the container array is passed into bspline's
    constructor, this view is 'adopted' and all operations will use the data it refers to.
    The caller is responsible for keeping these data alive while the bspline object exists,
    and relinquishes control over the data, which may be changed by the bspline object.

  - if data are passed to prefilter(), they will be taken as containing the knot point data,
    rather than expecting the knot point data to be in the bspline oject's memory already.
    This can also be used to reuse a bspline object with new data. The data passed in will
    not be modified. This is most efficient when using an implicit scheme; when used together
    with EXPLICIT, the data are (automatically) copied into the core area before prefiltering,
    which is unnecessary with the implicit schemes - they can 'pull in' data in the course
    of their operation.

  While there is no explicit code to create a 'smoothing spline' - a b-spline evaluating
  the source data without prefiltering them - this can be achieved simply by creating a b-spline
  object with spline degree 0 and 'shifting' it to the desired degree for evaluation. Note that
  you'll need the EXPLICIT strategy for the purpose, or specify extra 'headroom', because
  otherwise the spline won't have enough 'headroom' for shifting.
  
  If stronger smoothing is needed, this can be achieved with the code in filter.h, passing in
  appropriate pole values. a single-pole filter with a positive pole in ] 0 , 1 [ will do the
  trick - the larger the pole, the stronger the smoothing. Note that smoothing with large pole
  values will need a large 'horizon' as well to handle the margins properly. This is what's
  used when class bspline's creator is called with a 'smoothing' paramater greater than 0.

  With shifting, you can also create a 'poor man's pyramid'. While using no additional storage,
  you can extract smoothed data from the spline by shifting it up. This only goes so far, though,
  because even a degree-20 b-spline reconstruction kernel's equivalent gaussian doesn't have a very
  large standard deviation, and evaluation times become very long. From the gaussian approximation
  of the b-spline basis function, it can be seen that the equivalent gaussian's standard deviation
  is ( degree + 1 ) / 12.0, so a quintic spline will have a standard deviation of 0.5 only.
*/

#ifndef VSPLINE_BSPLINE_H
#define VSPLINE_BSPLINE_H

#include "prefilter.h"
#include "brace.h"

namespace vspline {

/// This enumeration is used to determine the prefiltering scheme to be used.

typedef enum { UNBRACED , ///< implicit scheme, no bracing applied
               BRACED ,   ///< implicit scheme, bracing will be applied
               EXPLICIT , ///< explicit scheme, frame with extrapolated signal, brace
               MANUAL     ///< like explicit, but don't frame before filtering
} prefilter_strategy  ;

const std::string pfs_name[] =
{
  "UNBRACED" ,
  "BRACED  " ,
  "EXPLICIT" ,
  "MANUAL  "
} ;

/// struct bspline is a convenience class which bundles a coefficient array (and it's creation)
/// with a set of metadata describing the parameters used to create the coefficients and the
/// resulting data. I have chosen to implement class bspline so that there is only a minimal
/// set of template arguments, namely the spline's data type (like pixels etc.) and it's dimension.
/// All other parameters relevant to the spline's creation are passed in at construction time.
/// This way, if explicit specialization becomes necessary (like, to interface to code which
/// can't use templates) the number of specializations remains manageable. This design decision
/// pertains specifically to the spline's degree, which could also be implemented as a template
/// argument, allowing for some optimization by making some members static. Yet going down this
/// path requires explicit specialization for every spline degree used and the performance gain
/// I found doing so was hardly measurable, while automatic testing became difficult and compilation
/// times grew.
///
/// I chose making bspline a struct for now, but messing with the data inside is probably
/// not a good idea.
///
/// class bspline may or may not 'own' the coefficient data it refers to - this
/// depends on the specific initialization used, but is handled privately by
/// class b-spline, using a shared_ptr to the data if they are owned, which
/// makes bspline objects trivially copyable.
///
/// There are two views in class bspline: the view 'core' corresponds with the
/// knot point values (aka the original data the spline was constructed over).
/// The view 'container' contains 'core' and an additional 'frame' around it
/// which allows for evaluating in-bounds coordinates without bounds checking.
/// For splines constructed with EXPLICIT prefilter strategy and/or additional
/// 'headroom' this frame may be even wider. We have a simple relationship:
/// core.shape() + left_frame + right_frame == container.shape()

template < class _value_type , int _dimension >
struct bspline
{
  typedef _value_type value_type ;
  /// pull the template arg into an enum
  enum { dimension = _dimension } ;
  /// if the coefficients are owned, this array holds the data
  typedef vigra::MultiArray < dimension , value_type > array_type ;
  /// data are read and written to vigra MultiArrayViews
  typedef vigra::MultiArrayView < dimension , value_type > view_type ;
  /// multidimensional index type
  typedef typename view_type::difference_type shape_type ;
  /// nD type for one boundary condition per axis
  typedef typename vigra::TinyVector < bc_code , dimension > bcv_type ;

  /// elementary type of value_type, like float or double
  typedef typename ExpandElementResult < value_type >::type real_type ;
  enum { channels = ExpandElementResult < value_type >::size } ;

  typedef bspline < real_type , dimension > channel_view_type ;
  
  // for internal calculations in the filter, we use the elementary type of value_type.
  // Note how in class bspline we make very specific choices about the
  // source data type, the target data type and the type used for arithmetics: we use
  // the same value_type for source and target array and perform the arithmetics with
  // this value_type's elementary type. The filter itself is much more flexible; all of
  // the three types can be different, the only requirements are that the value_types
  // must be vigra element-expandable types with an elementary type that can be cast
  // to and from math_type, and math_type must be a real data type, with the additional
  // requirement that it can be vectorized by VC if Vc is used.

  typedef real_type math_type ; // arbitrary, can use float or double here.

private:

  // _p_coeffs points to a vigra::MultiArray, which is either default-initialized
  // and contains no data, or holds data which is used by 'container'. Using a
  // std::shared_ptr here has the pleasant side-effect that class bspline objects 
  // can use the default copy and assignment operators.
  
  std::shared_ptr < array_type > _p_coeffs ;
  prefilter_strategy strategy ;

public:

  const bcv_type bcv ;      ///< boundary conditions, see common.h
  int spline_degree ;       ///< degree of the spline (3 == cubic spline)

  double tolerance ;        ///< acceptable error
  double smoothing ;        ///< E ] 0 : 1 [; apply smoothing to data before prefiltering
  bool braced ;             ///< whether coefficient array is 'braced' or not
  int horizon ;             ///< additional frame size for explicit scheme
  
  view_type container ;     ///< view to container array
  view_type core ;          ///< view to the core part of the coefficient array
  
  shape_type left_frame ;   ///< total width(s) of the left handside frame
  shape_type right_frame ;  ///< total width(s) of the right handside frame

  /// lower_limit returns the lower bound of the spline's defined range.
  /// This is usually 0.0, but with REFLECT boundary condition it's -0.5,
  /// the lower point of reflection. The lowest coordinate at which the
  /// spline can be accessed may be lower: even splines have wider support,
  /// and splines with extra headroom add even more room to manoevre.
  
  // TODO might introduce code to provide the 'technical limits'

  long double lower_limit ( const int & axis ) const
  {
    long double limit = 0.0L ;
    
    if ( bcv [ axis ] == vspline::REFLECT )
      limit = -0.5L ;
    else if ( bcv [ axis ] == vspline::SPHERICAL )
      limit = -0.5L ;
    
    return limit ;
  }
  
  vigra::TinyVector < long double , dimension > lower_limit() const
  {
    vigra::TinyVector < double , dimension > limit ;
    for ( int d = 0 ; d < dimension ; d++ )
      limit[d] = lower_limit ( d ) ;
    return limit ;
  }
  
  /// upper_limit returns the upper bound of the spline's defined range.
  /// This is normally M - 1 if the shape for this axis is M. Splines with
  /// REFLECT boundary condition use M - 0.5, the upper point of reflection,
  /// and periodic splines use M. The highest coordinate at which the spline
  /// may be accessed safely may be higher.
  
  long double upper_limit ( const int & axis ) const
  {
    long double limit = core.shape ( axis ) - 1 ;
    
    if ( bcv [ axis ] == vspline::REFLECT )
      limit += 0.5L ;
    else if ( bcv [ axis ] == vspline::SPHERICAL )
      limit += 0.5L ;
    else if ( bcv [ axis ] == vspline::PERIODIC )
      limit += 1.0L ;
    
    return limit ;
  }
  
  vigra::TinyVector < long double , dimension > upper_limit() const
  {
    vigra::TinyVector < long double , dimension > limit ;
    for ( int d = 0 ; d < dimension ; d++ )
      limit[d] = upper_limit ( d ) ;
    return limit ;
  }
  
  /// get_left_brace_size and get_right_brace_size calculate the size of
  /// the brace vspline puts around the 'core' coefficients to allow evaluation
  /// inside the defined range (and even slightly beyond) without bounds
  /// checking. These routines are static to allow user code to establish
  /// vspline's bracing requirements without creating a bspline object.
  /// user code might use this information to generate coefficient arrays
  /// suitable for use with vspline evaluation code, sidestepping use of
  /// a bspline object.

  static shape_type get_left_brace_size ( int spline_degree , bcv_type bcv )
  {
    int support = spline_degree / 2 ;
    // we start out with left_brace as large as the support
    // of the reconstruction kernel
    shape_type left_brace ( support ) ;
    // for some situations, we extend the array further along a specific axis
    for ( int d = 0 ; d < dimension ; d++ )
    {
      // If the spline is done with REFLECT or SPHERICAL
      // boundary conditions, the lower and upper limits are
      // between bounds.
      // the lower limit in this case is -0.5. When using
      // floor() or round() on this value, we receive -1,
      // so we need to extend the left brace.
      // if rounding could be done so that -0.5 is rounded
      // towards zero, this brace increase could be omitted
      // for even splines, but this would also bring operation
      // close to dangerous terrain: an infinitesimal undershoot
      // would already produce an out-of-bounds access.
      
      if ( bcv[d] == vspline::REFLECT || bcv[d] == vspline::SPHERICAL )
      {
        left_brace[d] ++ ;
      }
      
      // for other boundary conditions, the lower limit is 0.
      // for odd splines,
      // as long as evaluation is at positions >= 0 this is
      // okay, but as soon as evaluation is tried with a value
      // even just below 0, we'd have an out-of-bounds access,
      // with potential memory fault. Rather than requiring
      // evaluation to never undershoot, we err on the side
      // of caution and extend the left brace, so that
      // quantization errors won't result in a crash.
      // This is debatable and could be omitted, if it can
      // be made certain that evaluation will never be tried
      // at values even infinitesimally below zero.
      // for even splines, this problem does not exist, since
      // coordinate splitting is done with std::round.
      
      else if ( spline_degree & 1 )
      {
        left_brace[d]++ ;
      }
    }
    return left_brace ;
  }
  
  static shape_type get_right_brace_size ( int spline_degree , bcv_type bcv )
  {
    int support = spline_degree / 2 ;
    // we start out with right_brace as large as the support
    // of the reconstruction kernel
    shape_type right_brace ( support ) ;
    // for some situations, we extend the array further along a specific axis
    for ( int d = 0 ; d < dimension ; d++ )
    {
      // If the spline is done with REFLECT or SPHERICAL
      // boundary conditions, the lower and upper limits are
      // between bounds.
      // So the upper limit is Z + 0.5 where Z is integer.
      // using floor on this value lands at Z, which fine,
      // but using round() (as is done for even splines)
      // lands at Z+1, so for this case we need to extend
      // the right brace. If we could use a rounding mode
      // rounding towards zero, we could omit this extension,
      // but we'd also be cutting it very fine.
      
      if ( bcv[d] == vspline::REFLECT || bcv[d] == vspline::SPHERICAL )
      {
        if ( ! ( spline_degree & 1 ) )
          right_brace[d] ++ ;
      }
      
      // The upper limit is M-1 for most splines, and M-1+0.5 for
      // splines with REFLECT and SPHERICAL BCs. When accessing
      // the spline at this value, we'd be out of bounds.
      // For splines done with REFLECT and SPHERICAL BCs, we have
      // to extend the right brace to allow access to coordinates
      // in [M-1,M-1+0.5], there is no other option.
      // For other splines, We could require the evaluation code
      // to check and split incoming values of M-1 to M-2, 1.0, but this
      // would require additional inner-loop code. So we add another
      // coefficient on the upper side for these as well.
      // This is debatable, but with the current implementation of the
      // evaluation it's necessary.

      // So, erring on the side of caution, we add the extra coefficient
      // for all odd splines.

      if ( spline_degree & 1 )
      {
        right_brace[d]++ ;
      }
      
      // periodic splines need an extra coefficient on the upper
      // side, to allow evaluation in [M-1,M]. This interval is
      // implicit in the original data since the value at M is
      // equal to the value at 0, but if we want to process without
      // bounds checking and index manipulations, we must provide
      // an extra coefficient.
      
      if ( bcv[d] == vspline::PERIODIC )
      {
        right_brace[d]++ ;
      }
    }
    return right_brace ;
  }
  
  /// construct a bspline object with appropriate storage space to contain and process an array
  /// of knot point data with shape core.shape. Depending on the strategy chosen and the other
  /// parameters passed, more space than core.shape may be allocated. Once the bspline object
  /// is ready, usually it is filled with the knot point data and then the prefiltering needs
  /// to be done. This sequence assures that the knot point data are present in memory only once,
  /// the prefiltering is done in-place. So the user can create the bspline, fill in data (like,
  /// from a file), prefilter, and then evaluate.
  ///
  /// alternatively, if the knot point data are already manifest elsewhere, they can be passed
  /// to prefilter(). With this mode of operation, they are 'pulled in' during prefiltering.
  ///
  /// It's possible to pass in a view to an array providing space for the coefficients,
  /// or even the coefficients themselves. This is done via the parameter _space. This has
  /// to be an array of the same or larger shape than the container array would end up having
  /// given all the other parameters. This view is then 'adopted' and subsequent processing
  /// will operate on it's data.
  ///
  /// with the EXPLICIT scheme, the horizon is set by default to a value which is
  /// deemed to be 'sufficiently large' to keep the error 'low enough'. the expression
  /// used here produces a frame which is roughly the size needed to make any margin
  /// effects vanish by the time the prefilter hits the core, but it's a bit 'rule of thumb'.
  ///
  /// The additional parameter 'headroom' is used to make the 'frame' even wider. This is
  /// needed if the spline is to be 'shifted' up (evaluated as if it had been prefiltered
  /// with a higher-degree prefilter) - see shift(). If extreme precision is not an issue,
  /// when using the explicit extrapolation scheme, extra headroom for shifting is not
  /// necessary, since the data in the 'frame' are 'good enough' for the purpose.
  ///
  /// While bspline objects allow very specific parametrization, most use cases won't use
  /// parameters beyond the first few. The only mandatory parameter is, obviously, the
  /// shape of the knot point data, the original data which the spline is built over.
  /// This shape 'returns' as the bspline object's 'core' shape. If this is the only
  /// parameter passed to the constructor, the resulting bspline object will be a
  /// cubic b-spline with mirror boundary conditions, generated with an implicit
  /// extrapolation scheme to a 'good' quality, no smoothing, and allocating it's own
  /// storage for the coefficients, and the resuling bspline object will be suitable for
  /// use with vspline's evaluation code.
  
  // TODO: when bracing/framing is applied, we might widen the array size to a
  // multiple of the Vc:Vector's Size for the given data type to have better-aligned
  // access. This may or may not help, has to be tested. We might also want to position
  // the origin of the brace to an aligned position, since evaluation is based there.
  
  // TODO: while the coice to keep the value_types and math_type closely related makes
  // for simple code, with the more flexible formulation of the prefiltering code we might
  // widen class bspline's scope to accept input of other types and/or use a different
  // math_type.

  bspline ( shape_type _core_shape ,                ///< shape of knot point data
            int _spline_degree = 3 ,                ///< spline degree with reasonable default
            bcv_type _bcv = bcv_type ( MIRROR ) ,   ///< boundary conditions and common default
            prefilter_strategy _strategy = BRACED , ///< default strategy is the 'implicit' scheme
            int _horizon = -1 ,                 ///< width of frame for explicit scheme (-1: auto)
            double _tolerance = -1.0 ,          ///< acceptable error (relative to unit pulse) (-1: auto)
            double _smoothing = 0.0 ,               ///< apply smoothing to data before prefiltering
            int headroom = 0 ,                      ///< additional headroom, for 'shifting'
            view_type _space = view_type()          ///< coefficient storage to 'adopt'
          )
  : spline_degree ( _spline_degree ) ,
    bcv ( _bcv ) ,
    smoothing ( _smoothing ) ,
    strategy ( _strategy ) ,
    tolerance ( _tolerance )
  {
    if ( _tolerance < 0.0 )
    {
      // heuristic: 'reasonable' defaults
      // TODO derive from <limits>
      if ( std::is_same < real_type , float > :: value )
        tolerance = .000001 ;
      else if ( std::is_same < real_type , double > :: value )
        tolerance = .0000000000001 ;
      else
        tolerance = 0.0000000000000000001 ;
    }

    // heuristic: horizon for reasonable precision - we assume that no one in their right
    // minds would want a negative horizon ;)

    real_type max_pole = .1 ;
    if ( spline_degree > 1 )
      max_pole = fabs ( vspline_constants::precomputed_poles [ spline_degree ] [ 0 ] ) ;
    if ( smoothing > max_pole )
      max_pole = smoothing ;

    if ( _horizon < 0 )
      horizon = ceil ( log ( tolerance ) / log ( max_pole ) ) ; // TODO what if tolerance == 0.0?
    else
      horizon = _horizon ; // whatever the user specifies

    // first, calculate the shapes and sizes used internally
    
    shape_type container_shape = _core_shape ;
    left_frame = right_frame = shape_type() ;
    
    switch ( strategy )
    {
      case UNBRACED:
      case MANUAL:
        
        braced = false ;
        break ;

      case EXPLICIT:
        
        // for EXPLICIT extrapolation, we extend the coefficients
        // by 'horizon' so the filter can 'run up to precision'.
        // apart from that, everything is the same as for BRACED
        // prefilter strategy, so we fall through to that case.
        
        left_frame = horizon ;
        right_frame = horizon ;
        
      case BRACED:
        
        // next we add the size of the 'brace' which is needed to
        // perform evaluation without bounds checking, and any
        // additional headroom the caller has requested.

        braced = true ;
        
        left_frame += get_left_brace_size ( spline_degree , bcv ) ;
        left_frame += headroom ;

        right_frame += get_right_brace_size ( spline_degree , bcv ) ;
        right_frame += headroom ;
        
        // the container shape is the core shape plus the frame

        container_shape += left_frame + right_frame ;
        
        break ;
    }

    // now either adopt external memory or allocate memory for the coefficients
    
    if ( _space.hasData() )
    {
      // caller has provided space for the coefficient array. This space
      // has to be at least as large as the container_shape we have
      // determined to make sure it's compatible with the other parameters.
      // With the array having been provided by the caller, it's the caller's
      // responsibility to keep the data alive as long as the bspline object
      // is used to access them.

      if ( ! ( allGreaterEqual ( _space.shape() , container_shape ) ) )
        throw shape_mismatch ( "the intended container shape does not fit into the shape of the storage space passed in" ) ;
      
      // if the shape matches, we adopt the data in _space.
      // We take a view to the container_shape-sized subarray only.
      
      container = _space.subarray ( shape_type() , container_shape ) ;
      
      // _p_coeffs is made to point to a default-constructed MultiArray,
      // which holds no data.
      
      _p_coeffs = std::make_shared < array_type >() ;
    }
    else
    {
      // _space was default-constructed and has no data.
      // in this case we allocate a container array and hold a shared_ptr
      // to it. so we can copy bspline objects without having to worry about
      // dangling pointers, or who destroys the array.
      
      _p_coeffs = std::make_shared < array_type > ( container_shape ) ;
      
      // container is made to refer to a view to this array.
      
      container = *_p_coeffs ;
    }
    // finally we set the view to the core area
    
    core = container.subarray ( left_frame , left_frame + _core_shape ) ;
  } ;

  /// get a bspline object for a single channel of the data. This is lightweight
  /// and requires the viewed data to remain present as long as the channel view is used.
  /// the channel view inherits all metrics from it's parent, only the MultiArrayViews
  /// to the data are different.
  
  channel_view_type get_channel_view ( const int & channel )
  {
    assert ( channel < channels ) ;
    
    real_type * base = (real_type*) ( container.data() ) ;
    base += channel ;
    auto stride = container.stride() ;
    stride *= channels ;
    
    MultiArrayView < dimension , real_type >
      channel_container ( container.shape() , stride , base ) ;

    return channel_view_type ( core.shape() , 
                               spline_degree ,
                               bcv ,
                               strategy ,
                               horizon ,
                               tolerance ,
                               smoothing ,
                               0 ,
                               channel_container // coefficient storage to 'adopt'
                             ) ;
  } ;

  /// prefilter converts the knot point data in the 'core' area into b-spline
  /// coefficients. Depending on the strategy chosen in the b-spline object's
  /// constructor, bracing/framing may be applied. Even if the degree of the
  /// spline is zero or one, prefilter() should be called because it also
  /// performs the bracing, if needed.
  ///
  /// If data are passed in, they have to have precisely the shape
  /// we have set up in core (_core_shape passed into the constructor).
  /// These data will then be used in place of any data present in the
  /// bspline object to calculate the coefficients. They won't be looked at
  /// after prefilter() terminates, so it's safe to pass in some MultiArrayView
  /// which is destroyed after the call to prefilter() returns.

  void prefilter ( view_type data = view_type() ) ///< view to knot point data to use instead of 'core'
  {
    if ( data.hasData() )
    {
      // if the user has passed in data, they have to have precisely the shape
      // we have set up in core (_core_shape passed into the constructor).
      // This can have surprising effects if the container array isn't owned by the
      // spline but constitutes a view to data kept elsewhere (by passing _space the
      // to constructor): the data held by whatever constructed the bspline object
      // will be overwritten with the (prefiltered) data passed in via 'data'.
      
      if ( data.shape() != core.shape() )
        throw shape_mismatch
         ( "when passing data to prefilter, they have to have precisely the core's shape" ) ;
      if ( strategy == EXPLICIT )
      {
        // the explicit scheme requires the data and frame to be together in the
        // containing array, so we have to copy the data into the core.
        core = data ;
      }
      // the other strategies can move the data from 'data' into the spline's memory
      // during coefficient generation, so we needn't copy them in first.
    }
    else
    {
      // otherwise, we assume data are already in 'core' and we operate in-place
      // note, again, the semantics of the assignment here: since 'data' has no data,
      // the assignment results in 'adopting' the data in core rather than copying them
      data = core ;
    }

    // per default the output will be braced. This does require the output
    // array to be sufficiently larger than the input. class bracer handles
    // filling of the brace/frame

    bracer<view_type> br ;

    // for the explicit scheme, we use boundary condition 'guess' which tries to
    // provide a good guess for the initial coefficients with a small computational
    // cost. using zero-padding instead introduces a sharp discontinuity at the
    // margins, which we want to avoid.

    bcv_type explicit_bcv ( GUESS ) ;

    switch ( strategy )
    {
      case UNBRACED:
        // only call the prefilter, don't do any bracing. If necessary, bracing can be
        // applied later by a call to brace() - provided the bspline object has space
        // for the brace.
        solve < view_type , view_type , math_type >
              ( data ,
                core ,
                bcv ,
                spline_degree ,
                tolerance ,
                smoothing
              ) ;
        break ;
      case BRACED:
        // prefilter first, passing in BC codes to pick out the appropriate functions to
        // calculate the initial causal and anticausal coefficient, then brace result.
        // note how, just as in brace(), the whole frame is filled, which may be more
        // than is strictly needed by the evaluator.
        solve < view_type , view_type , math_type >
              ( data ,
                core ,
                bcv ,
                spline_degree ,
                tolerance ,
                smoothing
              ) ;
        // fill the frame
        for ( int d = 0 ; d < dimension ; d++ )
          br.apply ( container , bcv[d] , left_frame[d] , right_frame[d] , d ) ;
        break ;
      case EXPLICIT:
        // first fill frame using BC codes passed in, then solve with BC code GUESS
        // this automatically fills the brace as well, since it's part of the frame.
        // TODO: the values in the frame will not come out precisely the same as they
        // would by filling the brace after the coefficients have been calculated.
        // The difference will be larger towards the margin of the frame; we assume
        // that due to the small support of the evaluation the differences near the
        // margin of the core data will be negligible, having picked a sufficiently
        // large frame size. This is debatable. If it's a problem, a call to brace()
        // after prefilter() will brace again, now with coefficients from the core.
        for ( int d = 0 ; d < dimension ; d++ )
          br.apply ( container , bcv[d] , left_frame[d] , right_frame[d] , d ) ;
        solve < view_type , view_type , math_type >
              ( container ,
                container ,
                explicit_bcv ,
                spline_degree ,
                tolerance ,
                smoothing
              ) ;
        break ;
      case MANUAL:
        // like EXPLICIT, but don't apply a frame, assume a frame was applied
        // by external code. process whole container with GUESS BC. For cases
        // where the frame can't be constructed by applying any of the stock bracing
        // modes. Note that if any data were passed into this routine, in this case
        // they will be silently ignored (makes no sense overwriting the core after
        // having manually framed it in some way)
        solve < view_type , view_type , math_type >
              ( container ,
                container ,
                explicit_bcv ,
                spline_degree ,
                tolerance ,
                smoothing
              ) ;
        break ;
    }
  }

  /// if the spline coefficients are already known, they obviously don't need to be
  /// prefiltered. But in order to be used by vspline's evaluation code, they need to
  /// be 'braced' - the 'core' coefficients have to be surrounded by more coeffcients
  /// covering the support the evaluator needs to operate without bounds checking
  /// inside the spline's defined range. brace() performs this operation. brace()
  /// assumes the bspline object has been set up with the desired initial parameters,
  /// so that the boundary conditions and metrics are already known and storage is
  /// available. If brace() is called with an empty view (or without parameters),
  /// it assumes the coefficients are in the spline's core already and simply
  /// fills in the 'empty' space around them. If data are passed to brace(), they
  /// have to be the same size as the spline's core and are copied into the core
  /// before the bracing is applied - the same behaviour as prefilter exhibits.

  void brace ( view_type data = view_type() ) ///< view to knot point data to use instead of 'core'
  {
    if ( data.hasData() )
    {
      // if the user has passed in data, they have to have precisely the shape
      // we have set up in core

      if ( data.shape() != core.shape() )
        throw shape_mismatch
         ( "when passing data to brace(), they have to have precisely the core's shape" ) ;
      
      // we copy the data into the core
      core = data ;
    }

    // we use class bracer to do the work, filling the brace for all axes in turn

    bracer<view_type> br ;

    for ( int d = 0 ; d < dimension ; d++ )
      br.apply ( container , bcv[d] , left_frame[d] , right_frame[d] , d ) ;
  }

  /// overloaded constructor for 1D splines. This is useful because if we don't
  /// provide it, the caller would have to pass TinyVector < T , 1 > instead of T
  /// for the shape and the boundary condition.
  
  bspline ( long _core_shape ,                      ///< shape of knot point data
            int _spline_degree = 3 ,                ///< spline degree with reasonable default
            bc_code _bc = MIRROR ,                  ///< boundary conditions and common default
            prefilter_strategy _strategy = BRACED , ///< default strategy is the 'implicit' scheme
            int _horizon = -1 ,                     ///< width of frame for explicit scheme
            double _tolerance = -1.0 ,              ///< acceptable error (relative to unit pulse)
            double _smoothing = 0.0 ,               ///< apply smoothing to data before prefiltering
            int headroom = 0 ,                      ///< additional headroom, for 'shifting'
            view_type _space = view_type()          ///< coefficient storage to 'adopt'
          )
  :bspline ( TinyVector < long , 1 > ( _core_shape ) ,
             _spline_degree ,
             bcv_type ( _bc ) ,
             _strategy ,
             _horizon ,
             _tolerance ,
             _smoothing ,
             headroom ,
             _space
           )
  {
    static_assert ( _dimension == 1 , "bspline: 1D constructor only usable for 1D splines" ) ;
  } ;

  /// shift will change the interpretation of the data in a bspline object.
  /// d is taken as a difference to add to the current spline degree. The coefficients
  /// remain the same, but creating an evaluator from the shifted spline will make
  /// the evaluator produce data *as if* the coefficients were those of a spline
  /// of the changed order. Shifting with positive d will efectively blur the
  /// interpolated signal, shifting with negative d will sharpen it.
  /// For shifting to work, the spline has to have enough 'headroom', meaning that
  /// spline_degree + d, the new spline degree, has to be greater or equal to 0
  /// and smaller than the largest supported spline degree (lower twenties) -
  /// and, additionally, there has to bee a wide-enough brace to allow evaluation
  /// with the wider kernel of the higher-degree spline's reconstruction filter.
  /// So if a spline is set up with degree 0 and shifted to degree 5, it has to be
  /// constructed with an additional headroom of 3 (see the constructor).
  ///
  /// shiftable() is called with a desired change of spline_degree. If it
  /// returns true, interpreting the data in the container array as coefficients
  /// of a spline with the changed degree is safe. If not, the frame size is
  /// not sufficient or the resulting degree is invalid ans shiftable()
  /// returns false. Note how the decision is merely technical: if the new
  /// degree is okay and the *frame* is large enough, the shift will be
  /// considered permissible, even if the frame only holds data used for
  /// the EXPLICIT extrapolation scheme, which aren't really 'brace quality'.
  
  bool shiftable ( int d ) const
  {
    int new_degree = spline_degree + d ;
    if ( new_degree < 0 || new_degree > 24 )
      return false ;

    shape_type new_left_brace = get_left_brace_size ( new_degree , bcv ) ;
    shape_type new_right_brace = get_right_brace_size ( new_degree , bcv ) ;
    if (    allLessEqual ( new_left_brace , left_frame )
         && allLessEqual ( new_right_brace , right_frame ) )
    {
      return true ;
    }

    return false ;
  }
  
  /// shift() actually changes the interpretation of the data. The data
  /// will be taken to be coefficients of a spline with degree
  /// spline_degree + d, and the original degree is lost. This operation
  /// is only performed if it is technically safe (see shiftable()).
  /// If the shift was performed successfully, this function returns true,
  /// false otherwise.

  bool shift ( int d )
  {
    if ( shiftable ( d ) )
    {
      spline_degree += d ;
      return true ;
    }
    return false ;
  }

  /// helper function to << a bspline object to an ostream

  friend ostream& operator<< ( ostream& osr , const bspline& bsp )
  {
    osr << "dimension:................... " << bsp.dimension << endl ;
    osr << "degree:...................... " << bsp.spline_degree << endl ;
    osr << "boundary conditions:......... " ;
    for ( auto bc : bsp.bcv )
      osr << " " << bc_name [ bc ] ;
    osr << endl ;
    osr << endl ;
    osr << "shape of container array:.... " << bsp.container.shape() << endl ;
    osr << "shape of core:............... " << bsp.core.shape() << endl ;
    osr << "braced:...................... " << ( bsp.braced ? std::string ( "yes" ) : std::string ( "no" ) ) << endl ;
    osr << "horizon:..................... " << bsp.horizon << endl ;
    osr << "left frame:.................. " << bsp.left_frame << endl ;
    osr << "right frame:................. " << bsp.right_frame << endl ;
    osr << ( bsp._p_coeffs->hasData()
             ? "bspline object owns data"
             : "data are owned externally" ) << endl ;
    osr << "container base adress:....... " << bsp.container.data() << endl ;
    osr << "core base adress:............ " << bsp.core.data() << endl ;
    return osr ;
  }

} ;

/// using declaration for a coordinate suitable for bspline, given
/// elementary type rc_type. This produces the elementary type itself
/// if the spline is 1D, a TinyVector of rc_type otherwise.

template < class spline_type , typename rc_type >
using bspl_coordinate_type = typename std::conditional
                       < spline_type::dimension == 1 ,
                         rc_type ,
                         vigra::TinyVector
                         < rc_type , spline_type::dimension >
                       > :: type ;

/// using declaration for a bspline's value type

template < class spline_type >
using bspl_value_type = typename spline_type::value_type ;

} ; // end of namespace vspline

#endif // VSPLINE_BSPLINE_H