gradient.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 gradient.cc
///
/// \brief evaluating a specific spline, derivatives, precision
///
/// 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.
///
/// compile: clang++ -O3 -DUSE_VC -march=native -std=c++11 -pthread
///          -o gradient gradient.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 gradient gradient.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. This is more to make sure that the prefilter
  // does not do anything wrong - for the given signal, using 'brace()' would
  // be sufficient, because the prefilter on a linear gradient should not have
  // an effect on the coefficients (due to symmetry)
 
  bspl.prefilter() ;
 
  // set up the evaluator type
 
  typedef vigra::TinyVector < double , 3 > coordinate_type ;
  typedef vspline::evaluator < coordinate_type , double > evaluator_type ;
  
  // we also want to verify the derivative along each axis
  
  typedef typename evaluator_type::derivative_spec_type deriv_t ;
  deriv_t dsx , dsy , dsz ;
 
  dsx[0] = 1 ; // first derivative along axis 0
  dsy[1] = 1 ; // first derivative along axis 1
  dsz[2] = 1 ; // first derivative along axis 2
  
  // set up the evaluator for the underived result
  
  evaluator_type ev ( bspl ) ;
  
  // and evaluators for the three first derivatives
  
  evaluator_type ev_dx ( bspl , dsx ) ;
  evaluator_type ev_dy ( bspl , dsy ) ;
  evaluator_type ev_dz ( bspl , dsz ) ;
  
  // 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. The printout of the derivatives
  // is boring: it's always 1. But this assures us that the b-spline is
  // perfectly plane, even off the grid points. Towards the surface of the
  // volume, the derivative remains at 1.0 because we're using NATURAL
  // boundary conditions.
 
  for ( int times = 0 ; times < 100 ; times++ )
  {
    for ( int d = 0 ; d < bspl.dimension ; d++ )
      c[d] = ( core_shape[d] - 1 ) * std::generate_canonical<double, 20>(gen) ;
    double result ;
    ev.eval ( c , result ) ;
    double delta = result - sum ( c ) ;
 
    std::cout << "eval(" << c << ") = " << result
              << " -> delta = " << delta << std::endl ;
    
    ev_dx.eval ( c , result ) ;
    std::cout << "dx: " << result ;
    
    ev_dy.eval ( c , result ) ;
    std::cout << " dy: " << result ;
    
    ev_dz.eval ( c , result ) ;
    std::cout << " dz: " << result << std::endl ;
  }
}