Sun Jul 2 09:02:38 UTC 2017

This is an automated email from the git hooks/post-receive script.

kfj-guest pushed a commit to branch master
in repository vspline.

commit 3cde5edbc8759ba44e4947108d51b02bbe54e74b
Author: Kay F. Jahnke <kfjahnke at gmail.com>
Date:   Tue Dec 6 08:55:41 2016 +0100

    changed shifting code, smaller changes to pv
    I now save the initial boundary conditions used for prefiltering and revert to
    them during 'shifting' of the spline, because the shifting produced erroneous
    results with REFLECT bcs when the spline was later run with RTAG maping mode.
    pv now runs fine with all combinations of pan-mode and HQ splines and with
    full and partial sphericals - it takes a few extra command line arguments
    to fix the degrees for same, and has the arguments for partials.
---
 bspline.h | 97 ++++++++++++++++++++++++++++++++++++---------------------------
 common.h  |  2 ++
 eval.h    |  2 ++
 mapping.h | 88 ++++++++++++++++++++++++++++++++++++++++++++++++---------
 4 files changed, 133 insertions(+), 56 deletions(-)

diff --git a/bspline.h b/bspline.h
index d88e856..048b800 100644
--- a/bspline.h
+++ b/bspline.h
@@ -189,6 +189,7 @@ private:
 
   array_type _coeffs ;
   prefilter_strategy strategy ;
+  bcv_type boundary_conditions ;   ///< original boundary conditions, see common.h
 
 public:
 
@@ -196,7 +197,7 @@ public:
   view_type coeffs ;        ///< view to the braced coefficient array
   view_type core ;          ///< view to the core part of the coefficient array
   int spline_degree ;       ///< degree of the spline (3 == cubic spline)
-  bcv_type bcv ;            ///< boundary condition, see common.h
+  bcv_type bcv ;            ///< mapping mode, see common.h
   double tolerance ;        ///< acceptable error
   double smoothing ;        ///< E ] 0 : 1 [; apply smoothing to data before prefiltering
   bool braced ;             ///< whether coefficient array is 'braced' or not
@@ -316,6 +317,7 @@ public:
   : core_shape ( _core_shape ) ,
     spline_degree ( _spline_degree ) ,
     bcv ( _bcv ) ,
+    boundary_conditions ( _bcv ) ,
     smoothing ( _smoothing ) ,
     strategy ( _strategy )
   {
@@ -393,6 +395,11 @@ public:
                    int nthreads = ncores , ///< number of threads to use
                    view_type data = view_type() ) ///< view to knot point data to use instead of 'core'
   {
+    // if the user should have modified 'bcv' since the bspline object's creation,
+    // we want to have a record of the boundary conditions used for prefiltering:
+
+    boundary_conditions = bcv ;
+
     if ( data.hasData() )
     {
       // if the user has passed in data, they have to have precisely the shape
@@ -419,22 +426,6 @@ public:
       data = core ;
     }
 
-//     // we call the solver via a function pointer
-//     p_solve solve ;
-// 
-//     // for simplicity's sake, if USE_VC isn't defined we use solve_vigra
-//     // always and simply ignore the flag use_vc.
-// 
-// #ifdef USE_VC
-// //     we have two variants, one is using Vc and the other doesn't
-//     if ( use_vc )
-//       solve = & solve_vc < view_type , view_type > ;
-//     else
-//       solve = & solve_vigra < view_type , view_type > ;
-// #else
-//     solve = & solve_vigra < view_type , view_type > ;
-// #endif
-
     // per default the output will be braced. This does require the output
     // array to be sufficiently larger than the input; class bracer has code
     // to provide the right sizes
@@ -551,9 +542,16 @@ public:
   /// spline_degree + d, the new spline degree, has to be greater or equal to 0
   /// and smaller than the largest supported spline degree (lower twenties)
   /// This is a quick-shot solution to the problem of scaled-down interpolated
-  /// results; it may be faster in some situations to shift the spline up and
+  /// results; it may be better in some situations to shift the spline up and
   /// evaluate than to apply smoothing to the source data.
-  /// TODO: look into the transfer function
+  /// While the spline is 'shifted', it will revert to the 'original' boundary
+  /// conditions (those which were used for prefiltering or at creation of the
+  /// bspline object) during the operation. Normally, after prefiltering, the
+  /// only use of 'bcv' is to contain a mapping mode, which, per default, is the
+  /// same as the boundary conditions (so that evaluators can be created automatically).
+  /// But if bcv is changed to another mapping mode, the shift might not work out,
+  /// the point in case being REFLECT bondary conditions with their widened brace
+  /// for odd spline degrees.
 
   int shift ( int d )
   {
@@ -561,21 +559,32 @@ public:
     if ( new_degree < 0 || new_degree > 24 )
       return 0 ;
 
+    bcv_type bcv_backup ( bcv ) ; // save current mapping mode 
+    bcv = boundary_conditions ;   // revert to 'original' boundary conditions
+
     bracer<view_type> br ;
     shape_type new_left_brace = br.left_corner ( bcv , new_degree ) ;
     shape_type new_right_brace = br.right_corner ( bcv , new_degree ) ;
-    if ( ! (    allLessEqual ( new_left_brace , left_frame )
-             && allLessEqual ( new_right_brace , right_frame ) ) )
-      return 0 ;
-
-    spline_degree = new_degree ;
-    left_brace = new_left_brace ;
-    right_brace = new_right_brace ;
-    braced_shape = core_shape + left_brace + right_brace ;
+    if (    allLessEqual ( new_left_brace , left_frame )
+             && allLessEqual ( new_right_brace , right_frame ) )
+    {
+      // perform the shift
+      spline_degree = new_degree ;
+      left_brace = new_left_brace ;
+      right_brace = new_right_brace ;
+      braced_shape = core_shape + left_brace + right_brace ;
+
+      shape_type coefs_offset = left_frame - new_left_brace ;
+      coeffs.reset() ;
+      coeffs = container.subarray ( coefs_offset , coefs_offset + braced_shape ) ;
+    }
+    else
+    {
+      // can't shift
+      d = 0 ;
+    }
 
-    shape_type coefs_offset = left_frame - new_left_brace ;
-    coeffs.reset() ;
-    coeffs = container.subarray ( coefs_offset , coefs_offset + braced_shape ) ;
+    bcv = bcv_backup ; // restore saved mapping mode
     return d ;
   }
 
@@ -583,20 +592,24 @@ public:
 
   friend ostream& operator<< ( ostream& osr , const bspline& bsp )
   {
-    osr << "dimension:.................. " << bsp.dimension << endl ;
-    osr << "degree:..................... " << bsp.spline_degree << endl ;
-    osr << "boundary conditions:........" ;
+    osr << "dimension:................... " << bsp.dimension << endl ;
+    osr << "degree:...................... " << bsp.spline_degree << endl ;
+    osr << "original boundary conditions: " ;
+    for ( auto bc : bsp.boundary_conditions )
+      osr << " " << bc_name [ bc ] ;
+    osr << endl ;
+    osr << "mapping mode:................ " ;
     for ( auto bc : bsp.bcv )
       osr << " " << bc_name [ bc ] ;
     osr << endl ;
-    osr << "shape of container array:... " << bsp.container.shape() << endl ;
-    osr << "shape of braced coefficients " << bsp.coeffs.shape() << endl ;
-    osr << "shape of core:.............. " << bsp.core.shape() << endl ;
-    osr << "braced:..................... " << ( bsp.braced ? std::string ( "yes" ) : std::string ( "no" ) ) << endl ;
-    osr << "left brace:................. " << bsp.left_brace << endl ;
-    osr << "right brace:................ " << bsp.right_brace << endl ;
-    osr << "left frame:................. " << bsp.left_frame << endl ;
-    osr << "right frame:................ " << bsp.right_frame << endl ;
+    osr << "shape of container array:.... " << bsp.container.shape() << endl ;
+    osr << "shape of braced coefficients: " << bsp.coeffs.shape() << endl ;
+    osr << "shape of core:............... " << bsp.core.shape() << endl ;
+    osr << "braced:...................... " << ( bsp.braced ? std::string ( "yes" ) : std::string ( "no" ) ) << endl ;
+    osr << "left brace:.................. " << bsp.left_brace << endl ;
+    osr << "right brace:................. " << bsp.right_brace << endl ;
+    osr << "left frame:.................. " << bsp.left_frame << endl ;
+    osr << "right frame:................. " << bsp.right_frame << endl ;
     osr << ( bsp._coeffs.hasData() ? "bspline object owns data" : "data are owned externally" ) << endl ;
     return osr ;
   }
@@ -605,4 +618,4 @@ public:
 
 } ; // end of namespace vspline
 
-#endif // VSPLINE_BSPLINE_H
\ No newline at end of file
+#endif // VSPLINE_BSPLINE_H
diff --git a/common.h b/common.h
index 748e0ed..d1f4986 100644
--- a/common.h
+++ b/common.h
@@ -176,6 +176,7 @@ typedef enum { PERIODIC,   ///< periodic boundary conditions / periodic mapping
                LIMITP ,    ///< mapping mode, like LIMIT, but ceiling 1 higher
                TAG ,       ///< mapping mode, out-of-bounds coordinates produce delta -2.0
                TAGP ,      ///< mapping mode, like TAG, but ceiling 1 higher
+               RTAG ,      ///< mapping mode, like TAG, but for widened brace
                RAW   ,     ///< mapping mode, processes coordinates unchecked, may crash the program
 } bc_code;
 
@@ -209,6 +210,7 @@ std::vector < std::string > bc_name =
   "LIMITP" ,
   "TAG" ,
   "TAGP" ,
+  "RTAG" ,
   "RAW"
 } ;
 
diff --git a/eval.h b/eval.h
index 9874e72..6295651 100644
--- a/eval.h
+++ b/eval.h
@@ -67,6 +67,7 @@
 
 #include "mapping.h"
 #include "bspline.h"
+#include "interpolator.h"
 
 namespace vspline {
 
@@ -465,6 +466,7 @@ template < int dimension ,         ///< dimension of the spline
            class ic_type = int >   ///< singular integral coordinate, currently only int
 
 class evaluator
+: public interpolator < dimension , value_type , rc_type , ic_type >
 {
 public:
   
diff --git a/mapping.h b/mapping.h
index 746cb07..439f9e4 100644
--- a/mapping.h
+++ b/mapping.h
@@ -596,19 +596,73 @@ public:
     auto too_small = ( v < rc_t ( 0.0 ) ) ;
     auto mask = too_small | too_large ;
     
+    fv = v_modf ( v , &fl_i ) ; // split v into integral and remainder from [0...1[
+    iv = ic_v ( fl_i )  ;       // set integer part from float representing it
+
     if ( any_of ( mask ) )
     {
       // v has some out-of-bounds values
-      fv = v_modf ( v , &fl_i ) ;   // split v into integral and remainder from [0...1[
-      iv = ic_v ( fl_i )  ;         // set integer part from float representing it
-      // set fv to ic_out_of_bounds , rc_out_of_bounds at out-of-bounds values.
       iv ( typename ic_v::mask_type ( mask ) ) = ic_t ( ic_out_of_bounds ) ;
       fv ( mask ) = rc_t ( rc_out_of_bounds ) ;
-      return ;
     }
     
+  }
+
+#endif
+} ;
+
+template < typename ic_t , typename rc_t , int vsize = 1 >
+class odd_mapping_r_tag
+{
+  const rc_t _ceiling ;
+  const rc_t _ceilx2 ;
+  
+public:
+  
+  odd_mapping_r_tag ( int M )
+  : _ceiling ( M - 2 ) ,
+    _ceilx2 ( M + M - 4 )
+    { } ;
+  
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  {
+    v += rc_t ( 0.5 ) ;        // delete this to change back to origin at reflection
+    rc_t fl_i ;
+    if ( v < rc_t ( 0.0 ) || v > _ceiling )
+    {
+      iv = ic_t ( ic_out_of_bounds ) ;
+      fv = rc_t ( rc_out_of_bounds ) ;
+      return ;
+    }
+    // now we have v in [ 0 | _ceiling ]
+    // which corresponds to spline coordinates [ 0.5 , _ceiling + 0.5 ]
+    v += 0.5 ;
+    fv = modf ( v , &fl_i ) ; // split v into integral and remainder from [0...1[
+    iv = ic_t ( fl_i ) ;     // set integer part from float representing it
+  }
+  
+#ifdef USE_VC
+
+  typedef typename vector_traits < ic_t , vsize > :: type ic_v ;
+  typedef typename vector_traits < rc_t , vsize > :: type rc_v ;
+
+  /// this is the operator() for vectorized operation
+
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  {
+    rc_v fl_i ;
+    v += rc_t ( 0.5 ) ;
+    auto too_large = ( v > rc_t ( _ceiling ) ) ;
+    auto too_small = ( v < rc_t ( 0.0 ) ) ;
+    auto mask = too_small | too_large ;
+    v += rc_t(0.5) ;
     fv = v_modf ( v , &fl_i ) ; // split v into integral and remainder from [0...1[
-    iv = ic_v ( fl_i )  ;       // set integer part from float representing it
+    iv = ic_v ( fl_i )  ;      // set integer part from float representing it
+    if ( any_of ( mask ) )
+    {
+      iv ( typename ic_v::mask_type ( mask ) ) = ic_t ( ic_out_of_bounds ) ;
+      fv ( mask ) = rc_t ( rc_out_of_bounds ) ;
+    }
   }
 
 #endif
@@ -659,21 +713,17 @@ public:
     auto too_small = ( v < rc_t ( 0.0 ) ) ;
     auto mask = too_small | too_large ;
     
+    v += rc_t(0.5) ;
+    fv = v_modf ( v , &fl_i ) ; // split v into integral and remainder in [-0.5, 0.5]
+    fv -= rc_t(0.5) ;
+    iv = ic_v ( fl_i )  ;       // set integer part from float representing it
+
     if ( any_of ( mask ) )
     {
       // v has some out-of-bounds values
-      fv = v_modf ( v , &fl_i ) ;    // split v into integral and remainder from [0...1[
-      iv = ic_v ( fl_i )  ;         // set integer part from float representing it
-      // set iv to -1 at out-of-bounds values. mask cast is safe as sizes match.
       iv ( typename ic_v::mask_type ( mask ) ) = ic_t ( ic_out_of_bounds ) ;
       fv ( mask ) = rc_t ( rc_out_of_bounds ) ;
-      return ;
     }
-    // now we are sure that v is safely inside [ 0 : _ceiling [
-    v += rc_t(0.5) ;
-    fv = v_modf ( v , &fl_i ) ;           // split v into integral and remainder in [-0.5 ... 0.5]
-    fv -= rc_t(0.5) ;
-    iv = ic_v ( fl_i )  ;              // set integer part from float representing it
   }
 
 #endif
@@ -1255,6 +1305,11 @@ map_func pick_mapping ( bc_code bc , int spline_degree , int M )
         return odd_mapping_tag < ic_t , rc_t , vsize > ( M + 1 ) ;
         break ;
       }
+      case RTAG :
+      {
+        return odd_mapping_r_tag < ic_t , rc_t , vsize > ( M + 2 ) ;
+        break ;
+      }
       case CONSTANT :
       case NATURAL :
       {
@@ -1323,6 +1378,11 @@ map_func pick_mapping ( bc_code bc , int spline_degree , int M )
         return even_mapping_tag < ic_t , rc_t , vsize > ( M + 1 ) ;
         break ;
       }
+      case RTAG :
+      {
+        return even_mapping_tag < ic_t , rc_t , vsize > ( M ) ;
+        break ;
+      }
       case NATURAL :
       case CONSTANT :
       {

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/vspline.git

