[vspline] 18/72: polished mapping code now all the mappings resort to the raw mappings, reducing code redundancy. also some minor changes there, added more documentation.

Kay F. Jahnke kfj-guest at moszumanska.debian.org
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 79e725847b9f2c19f40b53cd4d618d71d6d6def2
Author: Kay F. Jahnke <kfjahnke at gmail.com>
Date:   Sat Dec 10 11:49:40 2016 +0100

    polished mapping code
    now all the mappings resort to the raw mappings, reducing code redundancy.
    also some minor changes there, added more documentation.
---
 common.h       |   6 +-
 eval.h         |  25 ++-
 interpolator.h |  10 +-
 mapping.h      | 628 +++++++++++++++++++--------------------------------------
 4 files changed, 245 insertions(+), 424 deletions(-)

diff --git a/common.h b/common.h
index d1f4986..74e4ff6 100644
--- a/common.h
+++ b/common.h
@@ -144,7 +144,11 @@ struct out_of_bounds
 /// upon, the code won't access memory outside the coefficient array. Using 1.0 for the
 /// remainder part is (with the widened brace on the right margin) a value which can't
 /// occur with coordinates inside the permitted range, but it's a value which will produce
-/// an in-range value when processed by the evaluation code.
+/// an in-range value when processed by the evaluation code. Finally, rv_out_of_bounds
+/// is the value assigned to an evaluation of a coordinate tagged with ic_out_of_bounds,
+/// rc_out_of_bounds. It's assigned if the evaluator object has template argument 'use_tag'
+/// set to true; if 'use_tag' is false, the tag values will be processed as if they were
+/// proper results, which may produce unwanted artifacts but is slightly faster.
 
 const long ic_out_of_bounds = 0 ;
 const double rc_out_of_bounds = 1.0 ;
diff --git a/eval.h b/eval.h
index f97cb8c..317d7c1 100644
--- a/eval.h
+++ b/eval.h
@@ -67,7 +67,7 @@
 
 #include "mapping.h"
 #include "bspline.h"
-#include "interpolator.h"
+// #include "interpolator.h" // currently unused
 
 namespace vspline {
 
@@ -466,7 +466,9 @@ template < int _dimension ,         ///< dimension of the spline
            class _ic_type ,         ///< singular integral coordinate, currently only int
            bool use_tag = false >   ///< flag switching tag checking on/off
 class evaluator
-: public interpolator < _dimension , _value_type , _rc_type , _ic_type >
+// we could inherit from class interpolator, since class evaluator satisfies the
+// interpolator interface, but there's no advantage in doing so.
+// : public interpolator < _dimension , _value_type , _rc_type , _ic_type >
 {
 public:
   
@@ -593,8 +595,8 @@ private:
   MultiArray < 1 , ic_type > component_offsets ;     ///< offsets in terms of ele_type, for vc op
   component_base_type component_base ;
   component_view_type component_view ;
-  nd_mapping_type mmap ;                 ///< mapping of real coordinates to spline coordinates
-  bspline_weight_functor_type wfd0 ;    ///< default weight functor: underived bspline
+  const nd_mapping_type mmap ;  ///< mapping of real coordinates to spline coordinates
+  bspline_weight_functor_type wfd0 ; ///< default weight functor: underived bspline
   const int spline_degree ;
   const int ORDER ;
   const int workspace_size ;
@@ -1323,6 +1325,21 @@ public:
     v_eval ( select , tune , result ) ;
   }
 
+//   void new_v_eval ( const nd_rc_v & input ,    // number of dimensions * coordinate vectors
+//                 mc_ele_v & result )  const // number of channels * value vectors
+//   {
+//     ic_v select ;
+//     nd_rc_v tune ;
+// 
+//     // map the incoming coordinates to the spline's range
+// 
+//     split_cordinate_v ( input , select , tune ) ;
+// 
+//     // delegate
+// 
+//     v_eval ( select , tune , result ) ;
+//   }
+
   /// This variant also operates on unsplit coordinates. Here again we require the calling function
   /// to pass pointers to vsize input and output values in contiguous memory. The data are
   /// (de)interleved in this routine, the actual evaluation is delegated to the variant working
diff --git a/interpolator.h b/interpolator.h
index caefe41..3aae68f 100644
--- a/interpolator.h
+++ b/interpolator.h
@@ -37,8 +37,13 @@
     In the context of vspline, this is used to present b-spline evaluators
     to the remap routines, but the interface can accomodate other
     interpolators as well.
-    
-    THIS Is A CONSTRUCTION SITE
+    This is meant to open up some of the remap routines to use alternative
+    interpolators, if their code isn't specific to b-splines. It might also
+    be possible to gain minimal speed improvements by handcoding stuff - I
+    handcoded nearest neighbour interpolation and the handcoded version seems
+    to run minimally faster (probably omits a futile multiplication with 1.0)
+    but it's not really a great deal.
+    Nevertheless, nice to have, no loss.
 */
 
 #ifndef VSPLINE_INTERPOLATOR_H
@@ -132,7 +137,6 @@ struct interpolator
 /// probably the simplest interpolation technique is 'nearest neighbour'
 /// interpolation. This serves as a model implementation of a class satisfying
 /// the interpolator interface:
-// TODO we're missing the treatment of values for out-of-range coordinates
 
 template < int _dimension ,
            class _value_type ,
diff --git a/mapping.h b/mapping.h
index f5c931a..dd2d70a 100644
--- a/mapping.h
+++ b/mapping.h
@@ -58,26 +58,11 @@
     to be found anywhere inside the defined range - it's instead the result of a mathematical
     operation. So the mapping in this case picks the value at the end of the defined range.
 
-    First we set up a few types to handle 'split coordinates'. A split coordinate consits
-    of an integral part and a fractional part. Typically, these coordinates are the result
-    of applying the modf operation to a real-valued coordinate. While commonly the 'splitting'
-    of the real-valued coordinates is performed during the spline evaluation, this can lead
-    to poorer performance, and if the splitting is done beforehand and the split coordinates
-    are reused, there is less arithmetic to perform - at the cost of higher memory use, since
-    split coordinates are bulkier.
-
-    The three classes we use in this context are all intended as n-dimensional objects.
-    The definition of 'split_type' as containing two n-dimensional components instead of
-    a TinyVector of 1D-split types is due to the fact that we need these separate components
-    in the evaluation: the first one determines the origin of the subarray of the braced spline
-    which will be processed into the result, and the second component constitutes a
-    'nd_fraction', which is needed to calculate the evaluator's weights. So we go with
-    the class of arrays approach instead an array of classs. A slight disadvantage here is
-    that 1D split coordinates are also formulated as manifolds (TinyVector < ... , 1  >).
-
+    First there is some collateral code to provide fmod and modf for vectors.
+    
     Next we define 'mappings'. Coordinates at which the spline is to be evaluated come in
-    originally as real values. These real values have to be split (therefore the types
-    above), but we also may want to perform the same mirroring or other boundary conditions
+    originally as real values. These real values have to be split, but we also may want to
+    perform the same mirroring or other boundary conditions
     on input coordinates, since we can only evaluate the spline within it's defined range.
     This is what mappings do. They apply the boundary conditions and then split the coordinate.
     Our approach to incoming coordinates outside the spline's defined range is this:
@@ -102,6 +87,13 @@
     access to memory would be more localized, since access within the same tile would be the
     likely case. I have some hopes for vigra's chunked arrays in this context.
 
+    Beyond the mappings which correspond to boundary conditions, there are 'mapping-only' mappings
+    Which don't try to mimick the application of a specific boundary condition, but provide some
+    sort of special treatment for coordinates outside the defined range. I provide code resulting
+    in an exception 'out_of_bounds' for out-of-bounds coordinates, code forcing the offending
+    vallues to the lowest/highest permissible value and code which 'tags' the result of the split
+    as out-of-bunds by setting the integral/remainder part to ic_out_of_bounds and rc_out_of_bounds,
+    respectively.
 */
 
 #ifndef VSPLINE_MAPPING_H
@@ -154,23 +146,16 @@ rc_v v_modf ( const rc_v& source ,
   }
 }
 
-/// for now, I'm taking the easy road to a vectorized fmod by using apply and a lambda expression
-/// TODO: handcoding this might increase performance
+/// ditto, for fmod.
 
 template <typename rc_v>
 rc_v v_fmod ( const rc_v& lhs ,
               const typename rc_v::EntryType rhs )
 {
-//   typedef Vc::SimdArray < int , rc_v::Size > ic_v ;
   typedef typename vector_traits < int , rc_v::Size > :: type ic_v ;
   
   ic_v iv ( lhs / rhs ) ;
   return lhs - rhs * rc_v ( iv ) ;
-//   return lhs.apply ( [&rhs] ( typename rc_v::EntryType lhs )
-//                      {
-//                        return fmod ( lhs , rhs ) ;
-//                      }
-//                    ) ;
 }
 
 #endif // USE_VC
@@ -179,7 +164,7 @@ rc_v v_fmod ( const rc_v& lhs ,
 /// processed by the solver, we use 'mapping' objects. These handle several tasks needed
 /// in this context:
 ///
-/// - they apply boundary conditions, like mirroring
+/// - they may apply boundary conditions, like mirroring
 /// - they split the resulting value into an integral and a real part
 /// - they take into account special requirements for odd and even splines
 ///
@@ -191,35 +176,26 @@ rc_v v_fmod ( const rc_v& lhs ,
 /// transformation routine on the coordinates first and then applying a minimal mapping
 /// like 'reject' or 'limit'.
 ///
-/// TODO: perform runtime measurements to see if the 'lumping together' really has a
-/// measurable performance advantage. If not, all mappings performing arithmetic could
-/// be abandoned, which would be cleaner design-wise.
-///
 /// I have made another design decision concerning the size of the bracing. Initially I was
 /// using a right brace which was as small as possible, but this forced me to check incoming
 /// coordinates to odd splines for v == M-1, because these then needed special treatment -
 /// the split coordinate had to be set to M-2, 1.0 instead of M-1, 0.0 which would have
 /// produced an out-of-bounds access. Now, for odd splines,  I use a right brace which is
-/// one slice larger, and accessing the spline at M-1, 0.0 is now safe, so the test can be omitted.
-/// I feel the sacrifice of one slice's worth of memory is worth the performance gain and increased
-/// code transparency.
+/// one larger, and accessing the spline at M-1, 0.0 is now safe, so the test can be omitted.
+/// I feel the sacrifice of one slice's worth of memory is worth the performance gain and
+/// increased code transparency.
 ///
 /// the simplest and fastest mapping is the 'raw' mapping. Here, only the split is performed,
 /// come what may. The user is responsible for the suitability of incoming coordinates, it is
-/// silently assumed that 0.0 <= v <= M-1; if this is not true, the results may be false or
+/// silently assumed that 0.0 <= v <= M-1. If this is not true, the results may be false or
 /// the program may crash.
-///
-/// Note how M is passed in but not used.
 
 template < typename ic_t , typename rc_t , int vsize = 1 >
 class odd_mapping_raw
 {
 public:
   
-  odd_mapping_raw ( int M )
-    { } ;
-    
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     rc_t fl_i ;
     fv = modf ( v , &fl_i ) ; // split v into integral and remainder from [0...1[
@@ -233,7 +209,7 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
     fv = v_modf ( v , &fl_i ) ; // split v into integral and remainder from [0...1[
@@ -243,19 +219,21 @@ public:
 #endif
 } ;
 
+/// the raw mapping for even splines first shifts the variable so that
+/// v == -0.5 'lands' at 0. Then it applies the modf function, resulting
+/// in an integral part and a remainder in [0,1], which is shifted back
+/// to range from [-0.5,0.5].
+
 template < typename ic_t , typename rc_t , int vsize = 1 >
 class even_mapping_raw
 {
 public:
 
-  even_mapping_raw ( int M )
-    { } ;
-  
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     rc_t fl_i ;
     fv = modf ( v + rc_t ( 0.5 ) , &fl_i ) - rc_t ( 0.5 ) ;
-    iv = ic_t ( fl_i ) ;     // set integer part from float representing it
+    iv = ic_t ( fl_i ) ;
   }
 
 #ifdef USE_VC
@@ -265,21 +243,21 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
     v += rc_t ( 0.5 ) ;
-    fv = v_modf ( v , &fl_i ) ; // split v into integral and remainder in [-0.5 ... 0.5]
+    fv = v_modf ( v , &fl_i ) ;
     fv -= rc_t ( 0.5 ) ;
-    iv = ic_v ( fl_i )  ;       // set integer part from float representing it
+    iv = ic_v ( fl_i )  ;
   }
 
 #endif
 } ;
 
-// All remaining mappings in this file are safe insofar as out-of-bounds incoming coordinates
-// will either result in an exception or be handled in a way specific to the mapping.
-
+/// All remaining mappings in this file are safe insofar as out-of-bounds incoming coordinates
+/// will either result in an exception or be handled in a way specific to the mapping.
+///
 /// the next mapping mode is LIMIT. Here any out-of-bounds coordinates are set to the
 /// nearest valid value.
 
@@ -294,7 +272,7 @@ public:
   : _ceiling ( M - 1 )
     { } ;
     
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     if ( v < 0.0 )
     {
@@ -302,15 +280,15 @@ public:
       fv = 0.0 ;
       return ;               // in this case we're done prematurely
     }
-    else if ( v > _ceiling )
+    
+    if ( v > _ceiling )
     {
-      iv = ic_t ( _ceiling ) ;
+      iv = ic_t ( _ceiling ) ; // above 'ceiling', pass ceiling, 0.
       fv = 0.0 ;
       return ;
-    }    
-    rc_t fl_i ;
-    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
+    }   
+    
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #ifdef USE_VC
@@ -320,20 +298,24 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
+    // set any out-of-bounds values in v to the closest endpoint value
+    
     v ( v > _ceiling ) = _ceiling ;
     
     v ( v < rc_t ( 0.0 ) ) = rc_t ( 0.0 ) ;
     
-    rc_v fl_i ;
-    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
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #endif
 } ;
 
+/// I use the same rejection criterion here as in odd splines, which is debatable:
+/// one might argue that coordinate values from -0.5 to 0.0 and M - 1.0 to M - 0.5
+/// can be processed meaningfully and hence should not be rejected.
+
 template < typename ic_t , typename rc_t , int vsize = 1 >
 class even_mapping_limit
 {
@@ -345,7 +327,7 @@ public:
   : _ceiling ( M - 1 )
     { } ;
   
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     if ( v > _ceiling )
       v = rc_t ( _ceiling ) ;
@@ -353,9 +335,7 @@ public:
     else if ( v < 0 )
       v = rc_t ( 0.0 ) ;
 
-    rc_t fl_i ;
-    fv = modf ( v + rc_t(0.5) , &fl_i ) - rc_t(0.5) ;
-    iv = ic_t ( fl_i ) ;
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #ifdef USE_VC
@@ -365,37 +345,21 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     v ( v > _ceiling ) = _ceiling ;
 
     v ( v < rc_t ( 0.0 ) ) = rc_t ( 0.0 ) ;
 
-    rc_v fl_i ;
-    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
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #endif
 } ;
 
 /// with mapping mode REJECT, out-of-bounds incoming coordinate values result in
-/// an exception (out_of_bounds). The calling code has to catch the exception if it
-/// wants to proceed. If the mapping is on vector types, some of the result may be
-/// valid. The out-of-bounds results can be recognized by the int part being set to
-/// 0 and the float part to rc_out_of_bounds, which is a value which can't be valid.
-/// Any out-of-bounds element triggers the exception.
-/// I tried out code which caught the exception in the evaluator and assigned some
-/// value (like 0 or 128) to the result for out-of-bounds coordinates. While, as long
-/// as the out-of-bounds exception didn't occur, the code was fast, if the exception
-/// was thrown and caught, things slowed down a great deal. So for this kind of
-/// processing the next mapping (TAG, TAGP) is better, since it has next to no
-/// run-time cost for out-of-bounds values. REJECT mapping, throwing exceptions, should
-/// be used if the calling code expects the exception won't occur at all and needs
-/// to do somethig on a higher level to deal with the problem.
-/// TODO: Therefor we might simply throw the exception and not bother with the masking.
+/// an exception (out_of_bounds). The exception is simply thrown, there is no attempt
+/// to pass information to the calling code beyond the fact of the offense.
 
 template < typename ic_t , typename rc_t , int vsize = 1 >
 class odd_mapping_reject
@@ -408,19 +372,14 @@ public:
   : _ceiling ( M - 1 )
     { } ;
     
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     if ( v < 0.0 || v > _ceiling )    // reject out-of-bounds values
     {
-      iv = ic_t ( ic_out_of_bounds ) ;
-      fv = rc_t ( rc_out_of_bounds ) ;
       throw out_of_bounds() ;
     }
 
-    // now we are sure that v is safely inside [ 0 : _ceiling ]
-    rc_t fl_i ;
-    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
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #ifdef USE_VC
@@ -430,26 +389,20 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
 
     auto too_large = ( v > rc_t ( _ceiling ) ) ;
     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 ) )
+    if ( any_of ( too_small | too_large ) )
     {
       // v has some out-of-bounds values
-      // 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 ) ;
       throw out_of_bounds() ;
     }
     
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #endif
@@ -470,19 +423,14 @@ public:
   : _ceiling ( M - 1 )
     { } ;
   
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     if ( v < 0.0 || v > _ceiling )    // reject out-of-bounds values
     {
-      iv = ic_t ( ic_out_of_bounds ) ;
-      fv = rc_t ( rc_out_of_bounds ) ;
       throw out_of_bounds() ;
     }
 
-    // split v into integral and remainder in [-0.5 ... 0.5]
-    rc_t fl_i ;
-    fv = modf ( v + rc_t(0.5) , &fl_i ) - rc_t(0.5) ;
-    iv = ic_t ( fl_i ) ;     // set integer part from float representing it
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #ifdef USE_VC
@@ -492,29 +440,20 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
     
     auto too_large = ( v > _ceiling ) ;
     auto too_small = ( v < rc_t ( 0.0 ) ) ;
-    auto mask = too_small | too_large ;
     
-    if ( any_of ( mask ) )
+    if ( any_of ( too_small | too_large ) )
     {
       // 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 ) ;
       throw out_of_bounds() ;
     }
-    // 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
+    
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #endif
@@ -548,19 +487,17 @@ public:
   : _ceiling ( M - 1 )
     { } ;
     
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv )  const
   {
-    if ( v < 0.0 || v > _ceiling )    // tag out-of-bounds values
+    if ( v < 0.0 || v > _ceiling )
     {
+      // tag out-of-bounds values and return
       iv = ic_t ( ic_out_of_bounds ) ;
       fv = rc_t ( rc_out_of_bounds ) ;
       return ;
     }
 
-    // now we are sure that v is safely inside [ 0 : _ceiling ]
-    rc_t fl_i ;
-    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
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #ifdef USE_VC
@@ -570,7 +507,7 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
 
@@ -578,12 +515,16 @@ 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
+    // since out-of-bounds values are acceptable with this mapping, we
+    // apply the raw mapping first, then overwrite those split values
+    // which resulted from out-of-bounds coordinates with the tag values
+    // ic_out_of_bounds and rc_out_of_bounds
+
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
 
     if ( any_of ( mask ) )
     {
-      // v has some out-of-bounds values
+      // v has some out-of-bounds values, tag them
       iv ( typename ic_v::mask_type ( mask ) ) = ic_t ( ic_out_of_bounds ) ;
       fv ( mask ) = rc_t ( rc_out_of_bounds ) ;
     }
@@ -593,6 +534,14 @@ public:
 #endif
 } ;
 
+/// r_tag mapping is a bit of an oddball. It is used for odd splines created
+/// with 'reflect' boundary condition. For this specific case, the range of
+/// acceptable coordinates goes from -0.5 to M - 0.5, so it's wider than
+/// ususal, and it uses a widened brace to enable accessing this range with
+/// an integral select value and a positive remainder. there is no even
+/// equivalent of this mapping, since even splines have smaller support and
+/// can handle the permitted range of coordinates without any trickery.
+
 template < typename ic_t , typename rc_t , int vsize = 1 >
 class odd_mapping_r_tag
 {
@@ -602,13 +551,13 @@ class odd_mapping_r_tag
 public:
   
   odd_mapping_r_tag ( int M )
-  : _ceiling ( M - 2 ) ,
-    _ceilx2 ( M + M - 4 )
+  : _ceiling ( M ) ,
+    _ceilx2 ( M + M )
     { } ;
   
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
-    v += rc_t ( 0.5 ) ;        // delete this to change back to origin at reflection
+    v += rc_t ( 0.5 ) ;
     rc_t fl_i ;
     if ( v < rc_t ( 0.0 ) || v > _ceiling )
     {
@@ -619,8 +568,8 @@ public:
     // 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
+
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
   
 #ifdef USE_VC
@@ -628,9 +577,7 @@ public:
   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 )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
     v += rc_t ( 0.5 ) ;
@@ -638,8 +585,9 @@ 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 from [0...1[
-    iv = ic_v ( fl_i )  ;      // set integer part from float representing it
+
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
+
     if ( any_of ( mask ) )
     {
       iv ( typename ic_v::mask_type ( mask ) ) = ic_t ( ic_out_of_bounds ) ;
@@ -650,7 +598,7 @@ public:
 #endif
 } ;
 
-/// I use the same rejectionion criterion here as in odd splines, which is debatable:
+/// I use the same rejection criterion here as in odd splines, which is debatable:
 /// one might argue that coordinate values from -0.5 to 0.0 and M - 1.0 to M - 0.5
 /// can be processed meaningfully and hence should not be tagged.
 
@@ -665,7 +613,7 @@ public:
   : _ceiling ( M - 1 )
     { } ;
   
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     if ( v < 0.0 || v > _ceiling )    // reject out-of-bounds values
     {
@@ -674,10 +622,7 @@ public:
       return ;
     }
 
-    // split v into integral and remainder in [-0.5 ... 0.5]
-    rc_t fl_i ;
-    fv = modf ( v + rc_t(0.5) , &fl_i ) - rc_t(0.5) ;
-    iv = ic_t ( fl_i ) ;     // set integer part from float representing it
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #ifdef USE_VC
@@ -685,9 +630,7 @@ public:
   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 )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
     
@@ -695,10 +638,7 @@ 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
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
 
     if ( any_of ( mask ) )
     {
@@ -711,11 +651,18 @@ public:
 #endif
 } ;
 
+/// Now we're done with mappings providing some sort of special treatment for
+/// out-of-bounds coordinates and move on to mappings which 'fold' out-of-bounds
+/// coordinates into the defined range. Here, the out-of-bounds coordinates are
+/// treated as acceptable input; the resulting split coordinates will not be
+/// distinguishable from values produced from coordinates which were already
+/// in-bounds in the first place.
+///
 /// the next mapping is for coefficients/data mirrored at the ends of the defined
 /// range. This is probably the most commonly used type and also the type which
-/// P. Thevenaz recommends. Like all mappings defined in this body of code, it comes
+/// P. Thevenaz recommends. Like most mappings defined in this body of code, it comes
 /// in two variants: one for odd splines and one for even ones. The handling is
-/// different for both cases: in an odd spline, the delta is in [0:1], in even
+/// different for both cases: for an odd spline, the delta is in [0:1], for even
 /// splines it's in [-0.5:+0.5].
 ///
 /// This is mirroring 'on the bounds':
@@ -736,11 +683,8 @@ public:
   : _ceilx2 ( M + M - 2 ) ,
     _ceiling ( M - 1 )
     { } ;
-    
-  // with odd splines we have to be careful not to access the coefficient matrix
-  // with iv == M - 1 and this results in extra safeguarding code.
-    
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     rc_t fl_i ;
     if ( v < 0.0 )              // apply mirror left boundary condition
@@ -753,9 +697,8 @@ public:
         v = _ceilx2 - v ;        // right border mirror
       }
     }
-    // now we are sure that v is safely inside [ 0 : _ceiling ]
-    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
+    
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #ifdef USE_VC
@@ -767,7 +710,7 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
     
@@ -779,13 +722,9 @@ public:
       v -= _ceiling ;                     // center
       v = abs ( v ) ;                     // map to half period
       v = _ceiling - v ;                  // flip
-      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
-      return ;
     }
-    // here we are sure that v is safely inside [ 0 : _ceiling ]
-    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
+    
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #endif
@@ -812,23 +751,21 @@ public:
     { } ;
   
   
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     rc_t fl_i ;
     if ( v < 0.0 )               // apply mirror left boundary condition
       v = -v ;
     if ( v > _ceiling )
     {
-      v = fmod ( v , _ceilx2 ) ; // apply right border mirror
+      v = fmod ( v , _ceilx2 ) ; // map to one full period
       if ( v > _ceiling )
       {
-        v = _ceilx2 - v ;        // no need to guard against too large v
+        v = _ceilx2 - v ;        // fold, no need to guard against too large v
       }
     }
-    // now v is <= _ceiling.
-    // split v into integral and remainder in [-0.5 ... 0.5]
-    fv = modf ( v + rc_t(0.5) , &fl_i ) - rc_t(0.5) ;
-    iv = ic_t ( fl_i ) ;     // set integer part from float representing it
+    
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #ifdef USE_VC
@@ -838,7 +775,7 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
     
@@ -850,11 +787,8 @@ public:
       v = abs ( v ) ;                   // map to half period
       v = _ceiling - v ;                // flip
     }
-    // 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
+    
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #endif
@@ -883,7 +817,7 @@ public:
   odd_mapping_periodic ( int M )
   : _ceiling ( M - 1 ) { } ;
   
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     rc_t fl_i ;
     if ( v < 0.0 )
@@ -894,8 +828,7 @@ public:
     {
       v = fmod ( v , _ceiling ) ; // force to first period (this also results in v <= _ceiling)
     }
-    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
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
   
 #ifdef USE_VC
@@ -905,16 +838,15 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
-    if ( any_of ( v < rc_t(0) ) || any_of ( v >= _ceiling ) )
+    if ( any_of ( ( v < rc_t(0) ) | ( v >= _ceiling ) ) )
     {
-      v = v_fmod ( v , _ceiling ) ;       // apply modulo (this also results in v <= _ceiling)
+      v = v_fmod ( v , _ceiling ) ;     // apply modulo (this also results in v <= _ceiling)
       v ( v < rc_t(0) ) += _ceiling ;   // shift values below zero up one period
     }
-    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
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #endif
@@ -930,7 +862,7 @@ public:
   even_mapping_periodic ( int M )
   : _ceiling ( M - 1 ) { } ;
   
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
     rc_t fl_i ;
     
@@ -942,9 +874,7 @@ public:
     {
       v = fmod ( v , _ceiling ) ; // force to first period
     }
-    // split v into integral and remainder in [-0.5 ... 0.5]
-    fv = modf ( v + rc_t(0.5) , &fl_i ) - rc_t(0.5) ;
-    iv = ic_t ( fl_i ) ;     // set integer part from float representing it
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #ifdef USE_VC
@@ -954,136 +884,15 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
     rc_v fl_i ;
-    if ( any_of ( v < rc_t(0) ) || any_of ( v >= _ceiling ) )
+    if ( any_of ( ( v < rc_t(0) ) | ( v >= _ceiling ) ) )
     {
-      v = v_fmod ( v , _ceiling ) ;     // apply modulo
+      v = v_fmod ( v , _ceiling ) ;   // apply modulo
       v ( v < rc_t(0) ) += _ceiling ; // shift values below zero up one period
     }
-    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
-} ;
-
-/// now the mapping for natural boundary conditions. Note that I use a generalization
-/// here: classically natural boundary conditions mean that all derivatives of the spline
-/// above the first are zero. Since the DSP approach to b-splines does not consider
-/// derivatives at the ends of the spline, there must be a different method to obtain the
-/// same result. I use point symmetry at the ends of the spline, which has the desired
-/// effect without using the derivatives explicitly.
-/// here we cannot obtain coordinates to data inside the defined range which would yield us
-/// values for the extrapolated signal (at least not in the adjoining area), since the
-/// values in the extrapolated signal are defined by
-/// f(x) - f(0) == f(0) - f(-x); f(x+n-1) - f(n-1) == f(n-1) - f (n-1-x)
-/// so there is an arithmetic operation involved which we can't represent by indexing.
-/// In other words, we can't 'fold' in the coordinates from outside the defined range
-/// as we could with periodic or mirror bcs.
-/// nevertheless we want to perform the usual transformations to real indices which
-/// correspond to data *inside* the defined range. What we do outside the defined range
-/// is arbitrary, so here I repeat the value at the bounds ad infinitum.
-/// TODO: this is silly, since it's not really specific to natural BCs.
-/// Might call it 'clamp' instead.
-
-template < typename ic_t , typename rc_t , int vsize = 1 >
-class odd_mapping_constant
-{
-  const rc_t _ceiling ;
-  
-public:
-  
-  odd_mapping_constant ( int M )
-  : _ceiling ( M - 1 ) { } ;
-  
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
-  {
-    rc_t fl_i ;
-    if ( v < 0.0 )
-    {
-      iv = 0 ;               // if v is below 0.0, we pass the value for v == 0.0
-      fv = 0.0 ;
-      return ;               // in this case we're done prematurely
-    }
-   else if ( v > _ceiling )
-    {
-      iv = ic_t ( _ceiling ) ;
-      fv = 0.0 ;
-      return ;
-    }
-    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 ( v < rc_t ( 0 ) ) = rc_t ( 0 ) ;
-    v ( v > _ceiling ) = _ceiling ;
-    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
-  }
-
-#endif
-} ;
-
-template < typename ic_t , typename rc_t , int vsize = 1 >
-class even_mapping_constant
-{
-  const rc_t _ceiling ;
-
-public:
-  
-  even_mapping_constant ( int M )
-  : _ceiling ( M - 1 ) { } ;
-  
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
-  {
-    if ( v < 0.0 )
-    {
-      iv = 0 ;               // if v is below 0.0, we pass the value for v == 0.0
-      fv = 0.0 ;
-      return ;               // in this case we're done prematurely
-    }
-   else if ( v > _ceiling )
-    {
-      iv = ic_t ( _ceiling ) ;
-      fv = 0.0 ;
-      return ;
-    }
-    // split v into integral and remainder in [-0.5 ... 0.5]
-    rc_t fl_i ;
-    fv = modf ( v + rc_t(0.5) , &fl_i ) - rc_t(0.5) ;
-    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 ( v < rc_t(0) ) = rc_t(0) ;
-    v ( v > _ceiling ) = _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
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #endif
@@ -1111,29 +920,16 @@ class odd_mapping_reflect
 public:
   
   odd_mapping_reflect ( int M )
-  : _ceiling ( M - 2 ) ,
-    _ceilx2 ( M + M - 4 )
+  : _ceiling ( M ) ,
+    _ceilx2 ( M + M )
     { } ;
   
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
-    v += rc_t ( 0.5 ) ;        // delete this to change back to origin at reflection
-    rc_t fl_i ;
-    if ( v < rc_t ( 0.0 ) )    // apply reflect left boundary condition
-      v = -v ;
-    if ( v > _ceiling )
-    {
-      v = fmod ( v , _ceilx2 ) ; // apply right border reflect
-      if ( v > _ceiling )
-      {
-        v = _ceilx2 - v ;
-      }
-    }
-    // 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
+    v += rc_t ( 0.5 ) ;        // shift coordinate to left point of reflection
+    odd_mapping_mirror<ic_t,rc_t,vsize> ( _ceiling + 1 ) ( v , iv , fv ) ; // mirror with wider range
+    v += rc_t ( 0.5 ) ;        // shift coordinate further to use with wider brace
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
   
 #ifdef USE_VC
@@ -1143,30 +939,24 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
-    rc_v fl_i ;
-    v += rc_t ( 0.5 ) ;
-    v = abs ( v ) ;
-    if ( any_of ( v > _ceiling ) )
-    {
-      v = v_fmod ( v , _ceilx2 ) ;        // map to one full period
-      v ( v > _ceiling ) = _ceilx2 - v ;
-    }
-    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
+    v += rc_t ( 0.5 ) ;        // shift coordinate to left point of reflection
+    odd_mapping_mirror<ic_t,rc_t,vsize> ( _ceiling + 1 ) ( v , iv , fv ) ; // mirror with wider range
+    v += rc_t ( 0.5 ) ;        // shift coordinate further to use with wider brace
+    odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #endif
 } ;
 
-/// mapping for REFLECT boundary conditions for even splines
-/// this mapping will map coordinates to [ -0.5 - M-0.5 ].
-/// for the even case, we don't need the widened brace, since it has smaller support,
-/// but we have to guard against the special case that v == M - 0.5, since this might
-/// be mapped to M, -0.5 which is outside the defined range. In this special case
-/// we use M-1, 0.5 instead which is equivalent and inside the range.
+///  mapping for REFLECT boundary conditions
+///
+///  for even splines this mapping will map coordinates to [ -0.5 - M-0.5 ]. for the 
+///  even case, we don't need the widened brace, since it has smaller support, but we 
+///  have to guard against the special case that v == M - 0.5, since this might be 
+///  mapped to M, -0.5 which is outside the defined range. In this special case we 
+///  use M-1, 0.5 instead which is equivalent and inside the range.
 
 template < typename ic_t , typename rc_t , int vsize = 1 >
 class even_mapping_reflect
@@ -1181,31 +971,12 @@ public:
     _ceilx2 ( M + M )
   { } ;
   
-  void operator() ( rc_t v , ic_t& iv , rc_t& fv )
+  void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
   {
-    v += rc_t ( 0.5 ) ;
-    rc_t fl_i ;
-    if ( v < 0.0 )              // apply reflect left boundary condition
-      v = -v ;
-    if ( v > _ceiling )
-    {
-      v = fmod ( v , _ceilx2 ) ; // apply right border reflect
-      if ( v > _ceiling )
-      {
-        v = _ceilx2 - v ;
-      }
-    }
-    if ( v >= _ceiling ) // guard against special == case; defensively use >= instead of ==
-    {
-      iv = _ceiling - 1 ;
-      fv = rc_t ( 0.5 ) ;
-      return ;
-    }
-    // now we have v in [ 0 | _ceiling ]
-    // which corresponds to spline coordinates [ -0.5 | _ceiling + 0.5 [
-    // split v into integral and remainder in [-0.5 ... 0.5]
-    fv = modf ( v , &fl_i ) - rc_t(0.5) ;
-    iv = ic_t ( fl_i ) ;     // set integer part from float representing it
+    v += rc_t ( 0.5 ) ;                 // shift coordinate origin to left point of reflection
+    even_mapping_mirror<ic_t,rc_t,vsize> ( _ceiling + 1 ) ( v , iv , fv ) ; // mirror with wider range
+    v -= rc_t ( 0.5 ) ;                 // shift coordinate back
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #ifdef USE_VC
@@ -1215,25 +986,12 @@ public:
 
   /// this is the operator() for vectorized operation
 
-  void operator() ( rc_v v , ic_v & iv , rc_v & fv )
+  void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
   {
-    rc_v fl_i ;
-    v += rc_t ( 0.5 ) ;
-    v = abs ( v ) ;
-    if ( any_of ( v > _ceiling ) )
-    {
-      v = v_fmod ( v , _ceilx2 ) ;        // map to one full period
-      v ( v > _ceiling ) = _ceilx2 - v ;
-    }
-    fv = v_modf ( v , &fl_i ) - rc_t(0.5) ;
-    iv = ic_v ( fl_i )  ;      // set integer part from float representing it
-    auto mask = ( v >= _ceiling ) ;  // guard against special == case; defensively use >= instead of ==
-    if ( any_of ( mask ) )
-    {
-      // mask cast is safe as sizes match
-      iv ( typename ic_v::mask_type ( mask ) ) = _ceiling - 1 ;
-      fv ( mask ) = rc_t ( 0.5 ) ;
-    }
+    v += rc_t ( 0.5 ) ;                 // shift coordinate origin to left point of reflection
+    even_mapping_mirror<ic_t,rc_t,vsize> ( _ceiling + 1 ) ( v , iv , fv ) ; // mirror with wider range
+    v -= rc_t ( 0.5 ) ;                 // shift coordinate back
+    even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
   }
 
 #endif
@@ -1246,56 +1004,94 @@ template < typename ic_t ,
            typename rc_t ,
            int vsize ,
            typename map_func >
-map_func pick_mapping ( bc_code bc , int spline_degree , int M )
+map_func pick_mapping ( bc_code bc ,        ///< mapping mode
+                        int spline_degree , ///< degree of the spline
+                        int M )             ///< shape along the axis to process
 {
+   // spline_degree is only needed to tell odd from even splines:
   if ( spline_degree & 1 )
   {
     switch ( bc )
     {
       case RAW :
       {
-        return odd_mapping_raw < ic_t , rc_t , vsize > ( M ) ;
+        // raw mapping simply splits the coordinate with no checks
+        // whatsoever. So it's constructor doesn't take an argument.
+        odd_mapping_raw < ic_t , rc_t , vsize > () ;
         break ;
       }
       case LIMIT :
       {
+        // limit mapping returns the nearest boundary value for
+        // out-of-bound coordinates.
         return odd_mapping_limit < ic_t , rc_t , vsize > ( M ) ;
         break ;
       }
       case LIMITP :
       {
+        // like limit mapping, but for periodic cases, where the defined
+        // range is one unit step larger than the shape along this axis
         return odd_mapping_limit < ic_t , rc_t , vsize > ( M + 1 ) ;
         break ;
       }
       case REJECT :
       {
+        // any out-of-bounds coordinate results in an
+        // 'out_of_bounds' exception
         return odd_mapping_reject < ic_t , rc_t , vsize > ( M ) ;
         break ;
       }
       case REJECTP :
       {
+        // like reject, but for periodic cases, where the defined
+        // range is one unit step larger than the shape along this axis
         return odd_mapping_reject < ic_t , rc_t , vsize > ( M + 1 ) ;
         break ;
       }
       case TAG :
       {
+        // with tag mapping, out-of-bounds incoming coordinates result in
+        // split coordinates with the special value pair rc_out_of_bounds,
+        // ic_out_of_bounds. This value combination is chosen so that it
+        // will not cause an out-of-bounds memory access if the calling
+        // code fails to react to these special values, but at the same
+        // time it can't occur for acceptable coordinates.
         return odd_mapping_tag < ic_t , rc_t , vsize > ( M ) ;
         break ;
       }
       case TAGP :
       {
+        // like tag, but for periodic cases, where the defined
+        // range is one unit step larger than the shape along this axis
         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 ) ;
+        // special treatment for odd splines with reflect boundary condition
+        return odd_mapping_r_tag < ic_t , rc_t , vsize > ( M ) ;
         break ;
       }
       case CONSTANT :
       case NATURAL :
       {
-        return odd_mapping_constant < ic_t , rc_t , vsize > ( M ) ;
+        /// now the mapping for natural boundary conditions. Note that I use a generalization
+        /// here: classically natural boundary conditions mean that all derivatives of the spline
+        /// above the first are zero. Since the DSP approach to b-splines does not consider
+        /// derivatives at the ends of the spline, there must be a different method to obtain the
+        /// same result. I use point symmetry at the ends of the spline, which has the desired
+        /// effect without using the derivatives explicitly.
+        /// here we cannot obtain coordinates to data inside the defined range which would yield us
+        /// values for the extrapolated signal (at least not in the adjoining area), since the
+        /// values in the extrapolated signal are defined by
+        /// f(x) - f(0) == f(0) - f(-x); f(x+n-1) - f(n-1) == f(n-1) - f (n-1-x)
+        /// so there is an arithmetic operation involved which we can't represent by indexing.
+        /// In other words, we can't 'fold' in the coordinates from outside the defined range
+        /// as we could with periodic or mirror bcs.
+        /// nevertheless we want to perform the usual transformations to real indices which
+        /// correspond to data *inside* the defined range. What we do outside the defined range
+        /// is arbitrary, so here I repeat the value at the bounds ad infinitum.
+        return odd_mapping_limit < ic_t , rc_t , vsize > ( M ) ;
         break ;
       }
       case MIRROR :
@@ -1306,7 +1102,7 @@ map_func pick_mapping ( bc_code bc , int spline_degree , int M )
       case REFLECT :
       case SPHERICAL :
       {
-        return odd_mapping_reflect < ic_t , rc_t , vsize > ( M + 2 ) ;
+        return odd_mapping_reflect < ic_t , rc_t , vsize > ( M ) ;
         break ;
       }
       case PERIODIC :
@@ -1327,7 +1123,7 @@ map_func pick_mapping ( bc_code bc , int spline_degree , int M )
     {
       case RAW :
       {
-        return even_mapping_raw < ic_t , rc_t , vsize > ( M ) ;
+        even_mapping_raw < ic_t , rc_t , vsize > () ;
         break ;
       }
       case LIMIT :
@@ -1368,7 +1164,7 @@ map_func pick_mapping ( bc_code bc , int spline_degree , int M )
       case NATURAL :
       case CONSTANT :
       {
-        return even_mapping_constant < ic_t , rc_t , vsize > ( M ) ;
+        return even_mapping_limit < ic_t , rc_t , vsize > ( M ) ;
         break ;
       }
       case MIRROR :
@@ -1397,8 +1193,8 @@ map_func pick_mapping ( bc_code bc , int spline_degree , int M )
 }
 
 /// struct mapping contains a std::function for singular mappings, and, if USE_VC
-/// is defined, a std::function for vectorized mappings. The specific function(s) used
-/// are determined by class pick_mapping.
+/// is defined, a std::function for vectorized mappings. The specific function(s)
+/// used are determined by class pick_mapping.
 
 template < typename ic_t , typename rc_t , int vsize = 1 >
 struct mapping
@@ -1491,7 +1287,7 @@ struct nd_mapping
 //   } ;
   
   /// apply the mapping along axis 'axis' to coordinate v, resulting in the setting
-  /// of the integral part iv and the fraczional part fv
+  /// of the integral part iv and the fractional part fv
   
   void operator() ( rc_t v , ic_t& iv , rc_t& fv , const int & axis ) const
   {

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



More information about the debian-science-commits mailing list