verify.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 verify.cc
 
    \brief verify bspline interpolation against polynomial
    
    A b-spline is a piecewise polynomial function. Therefore, it should
    model a polynomial of the same degree precisely. This program tests
    this assumption.
    
    While the test should hold throughout, we have to limit it to
    'reasonable' degrees, because we have to create the spline over
    a range sufficiently large to make margin errors disappear,
    so even if we want to, say, only look at the spline's values
    between 0 and 1, we have to have a few ten or even 100 values
    to the left and right of this interval, because the polynomial
    does not exhibit convenient features like periodicity or
    mirroring on the bounds. But since a polynomial, outside
    [-1,1], grows with the power of it's degree, the values around
    the test interval become very large very quickly with high
    degrees. We can't reasonable expect to calculate a meaningful
    spline over such data. The test shows how the measured fidelity
    degrades with higher degrees due to this effect.
    
    still , with 'reasonable' degrees, we see that the spline fits the
    signal very well indeed, demonstrating that the spline can
    faithfully represent a polynomial of equal degree.
*/
 
#include <iostream>
#include <typeinfo>
#include <random>
 
#include <vigra/multi_array.hxx>
#include <vigra/accumulator.hxx>
#include <vigra/multi_math.hxx>
 
#include <vspline/vspline.h>
 
template < class dtype >
struct random_polynomial
{
  int degree ;
  std::vector < dtype > coefficient ;
 
  // set up a polynomial with random coefficients in [0,1]
 
  random_polynomial ( int _degree )
  : degree ( _degree ) ,
    coefficient ( _degree + 1 )
  {
    std::random_device rd ;
    std::mt19937 gen ( rd() ) ;
    std::uniform_real_distribution<> dis ( -1.0 , 1.0 ) ;
    for ( auto & e : coefficient )
      e = dis ( gen ) ;
  }
  
  // evaluate the polynomial at x
 
  dtype operator() ( dtype x )
  {
    dtype power = 1 ;
    dtype result = 0 ;
    for ( auto e : coefficient )
    {
      result += e * power ;
      power *= x ;
    }
    return result ;
  }
} ;
 
template < class dtype >
void polynominal_test ( int degree , const char * type_name )
{
  // this is the function we want to model:
  
  random_polynomial < long double > rp ( degree ) ;
  
  // we evaluate this function in the range [-200,200[
  // note that for high degrees, the signal will grow very
  // large outside [-1,1], 'spoiling' the test
  
  vigra::MultiArray < 1 , dtype > data ( vigra::Shape1 ( 400 ) ) ;
 
  for ( int i = 0 ; i < 400 ; i++ )
  {
    dtype x = ( i - 200 ) ;
    data[i] = rp ( x ) ;
  }
  
  // we create a b-spline over these data
  
  vspline::bspline < dtype , 1 > bspl ( 400 , degree , vspline::NATURAL , 0.0 ) ;
  bspl.prefilter ( data ) ;
  
  auto ev = vspline::make_evaluator < decltype(bspl), dtype > ( bspl ) ;
 
  // to test the spline against the polynomial, we generate random
  // arguments in [-2,2]
 
  std::random_device rd ;
  std::mt19937 gen ( rd() ) ;
  std::uniform_real_distribution<long double> dis ( -2.0 , 2.0 ) ;
  
  long double signal = 0 ;
  long double spline = 0 ;
  long double noise = 0 ;
  long double error ;
  long double max_error = 0 ;
  
  // now we evaluate the spline and the polynomial at equal arguments
  // and do the statistics
 
  for ( int i = 0 ; i < 10000 ; i++ )
  {
    long double x = dis ( gen ) ;
    long double p = rp ( dtype ( x ) ) ;
    // note how we have to translate x to spline coordinates here
    long double s = ev ( dtype ( x + 200 ) ) ;
    error = std::fabs ( p - s ) ;
    if ( error > max_error )
      max_error = error ;
    signal += std::fabs ( p ) ;
    spline += std::fabs ( s ) ;
    noise += error ;
  }
  
  // note: with optimized code, this does not work:
 
  if ( std::isnan ( noise ) || std::isnan ( noise ) )
  {
    std::cout << type_name << " aborted due to numeric overflow" << std::endl ;
    return  ;
  }
  
  long double mean_error = noise / 10000.0L ;
  
  // finally we echo the results of the test
 
  std::cout << type_name
            << " Mean error: "
            << mean_error
            << " Maximum error: "
            << max_error
            << " SNR "
            << int ( 20 * std::log10 ( signal / noise ) )
            << "dB"
            << std::endl ;
}
 
int main ( int argc , char * argv[] )
{
  std::cout << std::fixed << std::showpos << std::showpoint
            << std::setprecision(18) ;
            
  for ( int degree = 1 ; degree < 15 ; degree++ )
  {
    std::cout << "testing spline against polynomial, degree " << degree << std::endl ;
    polynominal_test < float >       ( degree , "using float........" ) ;
    polynominal_test < double >      ( degree , "using double......." ) ;
    polynominal_test < long double > ( degree , "using long double.." ) ;
  }
}