gradient2.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.                                   */
/*                                                                      */
/************************************************************************/
 
/// gradient.cc
///
/// If we create a b-spline over an array containing, at each grid point,
/// the sum of the grid point's coordinates, each 1D row, column, etc will
/// hold a linear gradient with first derivative == 1. If we use NATURAL
/// BCs, evaluating the spline with real coordinates anywhere inside the
/// defined range should produce precisely the sum of the coordinates.
/// This is a good test for both the precision of the evaluation and it's
/// correct functioning, particularly with higher-D arrays.
///
/// In this variant of the program, we use a vspline::domain functor.
/// This functor provides a handy way to access a b-spline with normalized
/// coordinates: instead of evaluating coordinates in the range of
/// [ 0 , core_shape - 1 ] (natural spline coordinates), we pass coordinates
/// in the range of [ 0 , 1 ] to the domain object, which is chained to
/// the evaluator and passes natural spline coordinates to it.
///
/// compile: clang++ -O3 -DUSE_VC -march=native -std=c++11 -pthread \
///          -o gradient2 gradient2.cc -lVc
/// (with Vc; use -DUSE_HWY and -lhwy for highway, or -std=c++17 and -DUSE_STDSIMD
/// for the std::simd backend)
/// or clang++ -O3 -march=native -std=c++11 -pthread -o gradient2 gradient2.cc
 
#include <random>
#include <iostream>
 
#include <vspline/vspline.h>
 
int main ( int argc , char * argv[] )
{
  typedef vspline::bspline < double , 3 > spline_type ;
  typedef typename spline_type::shape_type shape_type ;
  typedef typename spline_type::view_type view_type ;
  typedef typename spline_type::bcv_type bcv_type ;
  
  // let's have a knot point array with nicely odd shape
 
  shape_type core_shape = { 35 , 43 , 19 } ;
  
  // we have to use a longish call to the constructor since we want to pass
  // 0.0 to 'tolerance' and it's way down in the argument list, so we have to
  // explicitly pass a few arguments which usually take default values before
  // we have a chance to pass the tolerance
 
  spline_type bspl ( core_shape ,                    // shape of knot point array
                     3 ,                             // cubic b-spline
                     bcv_type ( vspline::NATURAL ) , // natural boundary conditions
                     0.0 ) ;                         // no tolerance
 
  // get a view to the bspline's core, to fill it with data
 
  view_type core = bspl.core ;
  
  // create the gradient in each dimension
 
  for ( int d = 0 ; d < bspl.dimension ; d++ )
  {
    for ( int c = 0 ; c < core_shape[d] ; c++ )
      core.bindAt ( d , c ) += c ;
  }
  
  // now prefilter the spline
 
  bspl.prefilter() ;
 
  // set up an evaluator
 
  typedef vigra::TinyVector < double , 3 > coordinate_type ;
  typedef vspline::evaluator < coordinate_type , double > evaluator_type ;
  
  evaluator_type inner_ev ( bspl ) ;
  
  // create the domain from the bspline:
  
  auto dom = vspline::domain < coordinate_type >
    ( coordinate_type(0) , coordinate_type(1) , bspl ) ;
  
  // chain domain and inner evaluator
  
  auto ev = dom + inner_ev ;
  
  // we want to bombard the evaluator with random in-range coordinates
  
  std::random_device rd;
  std::mt19937 gen(rd());
  // std::mt19937 gen(12345);   // fix starting value for reproducibility
 
  coordinate_type c ;
  
  // here comes our test, feed 100 random 3D coordinates and compare the
  // evaluator's result with the expected value, which is precisely the
  // sum of the coordinate's components
 
  for ( int times = 0 ; times < 100 ; times++ )
  {
    for ( int d = 0 ; d < bspl.dimension ; d++ )
    {
      // note the difference to the code in gardient.cc here:
      // our test coordinates are in the range of [ 0 , 1 ]
      c[d] = std::generate_canonical<double, 20>(gen) ;
    }
    double result ;
    ev.eval ( c , result ) ;
    
    // 'result' is the same as in gradient.cc, but we have to calculate
    // the expected value differently.
    
    double delta = result - sum ( c * ( core_shape - 1 ) ) ;
 
    std::cout << "eval(" << c << ") = "
              << result << " -> delta = " << delta << std::endl ;
  }
}