int_spline.cc

Go to the documentation of this file.

/************************************************************************/
/*                                                                      */
/*    vspline - a set of generic tools for creation and evaluation      */
/*              of uniform b-splines                                    */
/*                                                                      */
/*            Copyright 2015 - 2023 by Kay F. Jahnke                    */
/*                                                                      */
/*    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 int_spline.cc
 
    \brief using a b-spline with integer coefficients
  
    vspline can use integral data types as b-spline coefficients,
    but doing so naively will produce imprecise results. This example
    demonstrates how to do it right. What needs to be done is using
    'boosted', or amplified, coefficients, so that the range of
    possible coefficients is 'spread out' to the range of the
    integral coefficient type. After evaluation, the 'raw' result
    value has to be attenuated back to the range of original values.
    
    You can pass the 'boost' factor on the command line to see how
    different values influence precision of the results. You'll also
    notice that, when passing too high 'boost' values, the results will
    turn wrong due to the int coefficients overflowing.
    
    Prefiltering will create a signal with higher peak values, so
    the maximal boost factor has to be chosen to take this into
    account, see 'max_boost' below.
*/
 
#include <iostream>
#include <random>
#include <limits>
#include <vspline/vspline.h>
 
/// given cf_type, the type for a b-spline coefficient, min and max
/// as the signal's minimum and maximum, and the spline's degree,
/// max_boost calculates the maximal safe boost to apply during
/// prefiltering so that the signal won't be damaged. This is
/// conservative, assuming the worst case (the value oscillating
/// between min and max). Instead of doing the 'proper' math, we
/// take the easy route: build a periodic b-spline over min and max
/// as the sole knot points, prefilter, get the coefficient with
/// the highest absolute value. Then divide the maximal value
/// of cf_type by this coefficient's absolute value.
 
// TODO might put this function into the library code
 
using vspline::xlf_type ;
 
template < typename cf_type >
static xlf_type max_boost ( xlf_type min , xlf_type max ,
                            const int & spline_degree )
{
  vspline::bspline < xlf_type , 1 > sample
    ( 2 , spline_degree , vspline::PERIODIC ) ;
    
  sample.core[0] = min ;
  sample.core[1] = max ;
  
  sample.prefilter() ;
  
  xlf_type d = std::max ( std::abs ( sample.core[0] ) ,
                          std::abs ( sample.core[1] ) ) ;
  
  if ( d != 0 )
    return std::numeric_limits<cf_type>::max() / d ;
 
  return 1 ;
}
 
int main ( int argc , char * argv[] )
{
  double min_x = -1.0 ;
  double max_x = 1.0 ;
  
  // per default, apply maximum boost
 
  double boost = max_boost<int> ( min_x , max_x , 3 ) ;
 
  // user has passed a different value on the CL; use that instead
 
  if ( argc > 1 )
    boost = std::atof ( argv[1] ) ;
  
  assert ( boost != 0.0 ) ;
  
  // we'll create a spline over two know points only:
 
  vigra::MultiArray < 1 , double > original ( 2 ) ;
  original[0] = 1 ;
  original[1] = -1 ;
  
  // create a bspline object holding int coefficients
  
  vspline::bspline < int , 1 > ibspl ( 2 , 3 , vspline::PERIODIC ) ;
  
  // prefilter with 'boost'
 
  ibspl.prefilter ( original , boost ) ;
 
  // create a reference b-spline with double coefficients over the
  // same data, using no 'boost'.
 
  vspline::bspline < double , 1 > bspl ( 2 , 3 , vspline::PERIODIC ) ;
  bspl.prefilter ( original ) ;
  
  // to compare the results of using the 'boosted' int spline with the
  // 'ordinary' b-spline with double coefficients, we create evaluators
  // of both types. the first evaluator, which takes int coefficients,
  // is piped to an amplify_type object, which multiplies with a factor
  // of 1/boost, putting the signal back in range.
  
  // we need these template arguments, in sequence:
  // evaluator taking float coordinates and yielding double, aggregating
  // into SIMD types with 8 elements, using general b-spline evaluation,
  // doing maths in double precision and using int spline coefficients:
  
  typedef vspline::evaluator
          < float , double , 8 , -1 , double , int > iev_t ;
  
  // create such an evaluator from 'ibspl', the b-spline holding
  // (boosted) int coefficients
          
  iev_t _iev ( ibspl ) ;
  
  // to 'undo' the boost we use an amplify_type object with a
  // factor which is the reciprocal of 'boost'
  
  auto attenuate = vspline::amplify_type
                   < double , double , double , 8 > ( 1.0 / boost ) ;
  
  // chain 'attenuate' to the evaluator, so that the evaluator's
  // output is multiplied with the factor in 'attenuate'
                   
  auto iev = _iev + attenuate ;
  
  // the equivalent evaluator using the spline 'bspl' which holds
  // double coefficients with no boost, first it's type, then the object
  
  typedef vspline::evaluator < float , double , 8 > ev_t ;
  
  ev_t ev ( bspl ) ;
  
  // now we call both evaluators at 0 and 1, the extrema,
  // expecting to see the results 1 and -1, respectively.
  // If boost is too low, we may find quantization errors,
  // if it is too high, the signal will be damaged.
  
  std::cout << "iev ( 0 ) = " << iev ( 0.0f ) << std::endl ;
  std::cout << " ev ( 0 ) = " <<  ev ( 0.0f ) << std::endl ;
  std::cout << "iev ( 1 ) = " << iev ( 1.0f ) << std::endl ;
  std::cout << " ev ( 1 ) = " <<  ev ( 1.0f ) << std::endl ;
  
  // just to doublecheck: vectorized operation with 'iev'
  
  auto iv = vspline::simdized_type < float , 8 > ( 0 ) ;
  std::cout << iv << " -> " << iev ( iv ) << std::endl ;
  iv = vspline::simdized_type < float , 8 > ( 1 ) ;
  std::cout << iv << " -> " << iev ( iv ) << std::endl ;
  
  std::cout << "used boost = " << boost << std::endl ;
}