splinus.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 splinus.cc
///
/// \brief compare a periodic b-spline with a sine
/// 
/// This is a simple example using a periodic b-spline
/// over just two values: 1 and -1. This spline is used to approximate
/// a sine function. You pass the spline's desired degree on the command
/// line. Next you enter a number (interpreted as degrees) and the program
/// will output the sine and the 'splinus' of the given angle.
/// As you can see when playing with higher degrees, the higher the spline's
/// degree, the closer the match with the sine. So apart from serving as a
/// very simple demonstration of using a 1D periodic b-spline, it teaches us
/// that a periodic b-spline can approximate a sine.
/// To show off, we use long double as the spline's data type.
/// This program also shows a bit of functional programming magic, putting
/// together the 'splinus' functor from several of vspline's functional
/// bulding blocks.
///
/// compile with: clang++ -pthread -O3 -std=c++11 splinus.cc -o splinus
 
#include <iostream>
#include <assert.h>
#include <vspline/vspline.h>
 
int main ( int argc , char * argv[] )
{
  if ( argc < 2 )
  {
    std::cerr << "please pass the spline degree on the command line"
              << std::endl ;
    exit ( 1 ) ;
  }
  
  int degree = std::atoi ( argv[1] ) ;
  
  if ( degree < 4 )
  {
    std::cout << "rising degree to 4, the minimum for this program" << std::endl ;
    degree = 4 ;
  }
  
  assert ( degree >= 4 && degree <= vspline_constants::max_degree ) ;
  
  // create the bspline object
  
  typedef vspline::bspline < long double , 1 > spline_type ;
  
  spline_type bsp ( 2 ,                   // two values
                    degree ,              // degree as per command line
                    vspline::PERIODIC ,   // periodic boundary conditions
                    0.0 ) ;               // no tolerance
          
  // the bspline object's 'core' is a MultiArrayView to the knot point
  // data, which we set one by one for this simple example:
  
  bsp.core[0] = 1.0L ;
  bsp.core[1] = -1.0L ;
  
  // now we prefilter the data
  
  bsp.prefilter() ;
  
  // we build 'splinus' as a functional construct. Inside the brace,
  // we 'chain' several vspline::unary_functors:
  // - a 'domain' which scales and shifts input to the spline's range.
  //   with the 'incoming' range of [ 90 , 450 ] and the spline's
  //   range of [ 0 , 2 ], the translation of incoming coordinates is:
  //   x' = 0 + 2 * ( x - 90 ) / ( 450 - 90 )
  // - a 'safe' evaluator for the spline. since the spline has been
  //   built with PERIODIC boundary conditions, this evaluator will
  //   map incoming coordinates into the first period, [0,2].
 
  auto splinus =
       (   vspline::domain ( 90.0L , 450.0L , bsp )
         + vspline::make_safe_evaluator < spline_type , long double > ( bsp ) ) ;
 
  // alternatively we can use this construct. This will work just about
  // the same, but has a potential flaw: If arithmetic imprecision should
  // land the output of the periodic gate just slightly below 90.0, the
  // domain may produce a value just below 0.0, resulting in a slightly
  // out-of-bounds access. So the construct above is preferable.
  // Just to demonstrate that vspline::grok also produces an object that
  // can be used with function call syntax, we use vspline::grok here.
  // Note how 'grokking' the chain of functors produces a simply typed
  // object, rather than the complexly typed result of the chaining
  // operation inside the brace:
      
  vspline::grok_type < long double , long double > splinus2
  = vspline::grok
    (   vspline::periodic ( 90.0L , 450.0L )
      + vspline::domain ( 90.0L , 450.0L , bsp )
      + vspline::evaluator < long double , long double > ( bsp ) ) ;
 
  // we throw derivatives in the mix. If our spline models a sine,
  // it's derivatives should model the sine's derivatives, cos etc.
  // note how this is obscured by the higher 'steepness' of the spline,
  // which is over [ 0 , 2 ] , not [ 0 , 2 * pi ]. Hence the derivatives
  // come out amplified with a power of pi, which we compensate for.
  
  vigra::TinyVector < int , 1 > derivative ;
  derivative[0] = 1 ;
  vspline::evaluator < long double , long double > evd1 ( bsp , derivative ) ; 
  derivative[0] = 2 ;
  vspline::evaluator < long double , long double > evd2 ( bsp , derivative ) ; 
  derivative[0] = 3 ;
  vspline::evaluator < long double , long double > evd3 ( bsp , derivative ) ; 
 
  auto splinus_d1
  = (   vspline::domain ( 90.0L , 450.0L , bsp )
      + vspline::periodic ( 0.0L , 2.0L )
      + evd1 ) ;
  
  auto splinus_d2
  = (   vspline::domain ( 90.0L , 450.0L , bsp )
      + vspline::periodic ( 0.0L , 2.0L )
      + evd2 ) ;
  
  auto splinus_d3
  = (   vspline::domain ( 90.0L , 450.0L , bsp )
      + vspline::periodic ( 0.0L , 2.0L )
      + evd3 ) ;
  
  // now we let the user enter arbitrary angles, calculate the sine
  // and the 'splinus' and print the result and difference:
 
  while ( std::cin.good() )
  {
    std::cout << " angle in degrees > " ;
    long double x ;
    std::cin >> x ;                       // get an angle
 
    if ( std::cin.eof() )
    {
      std::cout << std::endl ;
      break ;
    }
 
    long double xs = x * M_PI / 180.0L ;  // note: sin() uses radians 
    
    // finally we can produce both results. Note how we can use periodic_ev,
    // the combination of gate and evaluator, like an ordinary function.
 
    std::cout << "sin(" << x << ") = "
              << sin ( xs ) << std::endl
              << "cos(" << x << ") = "
              << cos ( xs ) << std::endl
              << "splinus(" << x << ") = "
              << splinus ( x ) << std::endl
              << "splinus2(" << x << ") = "
              << splinus2 ( x ) << std::endl
              << "difference sin/splinus: "
              << sin ( xs ) - splinus ( x ) << std::endl
              << "difference sin/splinus2: "
              << sin ( xs ) - splinus2 ( x ) << std::endl
              << "difference splinus/splinus2: "
              << splinus2 ( x ) - splinus ( x ) << std::endl
              << "derivatives: " << splinus_d1 ( x ) / M_PI
              << " " << splinus_d2 ( x ) / ( M_PI * M_PI )
              << " " << splinus_d3 ( x ) / ( M_PI * M_PI * M_PI ) << std::endl ;
 
  }
}