polish.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                    */
/*                                                                      */
/*    The git repository for this software is at                        */
/*                                                                      */
/*    https://bitbucket.org/kfj/vspline                                 */
/*                                                                      */
/*    Please direct questions, bug reports, and contributions to        */
/*                                                                      */
/*    kfjahnke+vspline@gmail.com                                        */
/*                                                                      */
/*    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 polish.cc
 
    \brief 'polish' a b-spline several times, to see if it's precision
          can be improved. It turns out that, starting out with zero
          tolerance for the coefficient generation, there is very little
          room for improvement.
*/
    
#include <vigra/multi_array.hxx>
#include <vigra/accumulator.hxx>
#include <vigra/multi_math.hxx>
#include <iostream>
#include <typeinfo>
#include <random>
#include <ctime>
#include <chrono>
 
#include <vspline/vspline.h>
 
bool verbose = true ; // false ;
 
// 'condense' aggregate types (TinyVectors etc.) into a single value
 
template < typename T >
double condense ( const T & t , std::true_type )
{
  return std::abs ( t ) ;
}
 
template < typename T >
double condense ( const T & t , std::false_type )
{
  return sqrt ( sum ( t * t ) ) / t.size() ;
}
 
template < typename T >
double condense ( const std::complex<T> & t , std::false_type )
{
  return std::abs ( t ) ;
}
 
template < class T >
using is_singular = typename
  std::conditional
  < std::is_fundamental < T > :: value ,
    std::true_type ,
    std::false_type
  > :: type ;
 
template < typename T >
double condense ( const T & t )
{
  return condense ( t , is_singular<T>() ) ;
}
 
template < int dim , typename T >
double check_diff ( vigra::MultiArrayView < dim , T > & reference ,
                    vigra::MultiArrayView < dim , T > & candidate )
{
  using namespace vigra::multi_math ;
  using namespace vigra::acc;
  
  assert ( reference.shape() == candidate.shape() ) ;
  
  vigra::MultiArray < 1 , double >
    error_array ( vigra::Shape1(reference.size() ) ) ;
    
  for ( int i = 0 ; i < reference.size() ; i++ )
  {
    auto help = candidate[i] - reference[i] ;
//     std::cerr << reference[i] << " <> " << candidate[i]
//               << " CFD " << help << std::endl ;
    error_array [ i ] = condense ( help ) ;
  }
 
  AccumulatorChain < double , Select < Mean, Maximum > > ac ;
  extractFeatures ( error_array.begin() , error_array.end() , ac ) ;
  double mean = get<Mean>(ac) ;
  double max = get<Maximum>(ac) ;
  if ( verbose )
  {
    std::cout << "delta Mean:    "
              << mean << std::endl;
    std::cout << "delta Maximum: "
              << max << std::endl;
  }
  return max ;
}
 
template < int dim , typename T > void
polish_test ( vigra::TinyVector < int , dim > shape ,
              vspline::bc_code bc ,
              int spline_degree )
{
  typedef vigra::MultiArray < dim , T > array_type ;
  typedef vspline::bspline < T , dim > spline_type ;
  
  array_type arr ( shape ) ;
  vigra::TinyVector < vspline::bc_code , dim > bcv { bc } ;
  spline_type bsp  ( shape , spline_degree , bcv , 0.0 ) ;
 
  if ( verbose )
    std::cout << "created b-spline:" << std::endl << bsp << std::endl ;
  
  std::random_device rd ;
  std::mt19937 gen ( rd() ) ;
//   gen.seed ( 765 ) ;              // level playing field
  std::uniform_real_distribution<> dis ( -1.0 , 1.0 ) ;
  for ( auto & e : arr )
    e = dis ( gen ) ;     // fill array with random data
  
  array_type reference = arr ; // hold a copy of these data
  
  bsp.prefilter ( arr ) ; // suck into b-spline via prefilter
  
  vspline::restore < dim , T > ( bsp , arr ) ; // restore back to arr
 
  if ( verbose )
    std::cout << "after restoration of original data:" << std::endl ;
  
  double emax = check_diff < dim , T > ( reference , arr ) ;
  
  if ( verbose )
  {
    // print a summary, can use '| grep CF' on cout
    std::cout
    << typeid(T()).name()
    << " CF "
//     << vspline::pfs_name [ pfs ][0]
    << " D " << dim
    << " " << arr.shape()
    << " BC " << vspline::bc_name[bc]
    << " DG " << spline_degree
    << " TOL " << bsp.tolerance
    << " EMAX "
    << emax << std::endl ;
  }
  
  for ( int times = 1 ; times < 5 ; times++ )
  {
    // difference original / restored
    arr -= reference ;
    // create another bspline
    spline_type polish_spl ( arr.shape() , spline_degree , bcv , 0.0 ) ;
    // prefilter it, sucking in the difference
    polish_spl.prefilter ( arr ) ;
    // and layer the resulting coeffs on top of the previous set
    bsp.core -= polish_spl.core ;
    // brace the 'augmented' spline
    bsp.brace() ;
    // and check it's quality
    vspline::restore < dim , T > ( bsp , arr ) ;
 
    if ( verbose )
      std::cout << "after polishing run " << times << std::endl ;
 
    double emax2 = check_diff < dim , T > ( reference , arr ) ;
    
    if ( verbose )
    {
      // print a summary, can use '| grep CF' on cout
      std::cout
      << typeid(T()).name()
      << " CF "
  //     << vspline::pfs_name [ pfs ][0]
      << " D " << dim
      << " " << arr.shape()
      << " BC " << vspline::bc_name[bc]
      << " DG " << spline_degree
      << " TOL " << bsp.tolerance
      << " EMAX "
      << emax2 << std::endl ;
    }
  }
}
 
int main ( int argc , char * argv[] )
{
  if ( argc < 2 )
  {
    std::cerr << "pass the spline degree as parameter" << std::endl ;
    exit ( -1 ) ;
  }
  
  int degree = std::atoi ( argv[1] ) ;
  
  polish_test < 2 , float >
    ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
    
  polish_test < 2 , double >
    ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
    
  polish_test < 2 , long double >
    ( { 1000 , 1000 } ,vspline::PERIODIC , degree ) ;
}