[vspline] 51/72: separated boundary conditions and mapping codes I was unhappy with both voundary conditions and mappings sharing the same enum and being conceptually linked. Now there is a separate set of mapping codes (in common.h) and the mapping code in mapping.h is based on them. When creating evaluators from bspline objects, the old automatism of picking a mapping which 'echoed' the spline's boundary conditions was removed, the default now is to use MAP_REJECT, which throws vspline::out_of_bounds for out-of-bounds coordinates - but any of the other mappings can be used as well by simply passing in the desired mapping code(s) or an mmap object. The coundary conditions in bspline objects are now const, and the old mechanism of changing the spline's bc's in order to get different mapping behaviour when automatically creating an evaluator from a bspline object can no longer be used, which is a relief, since it simplifies matters, especially when shifting the spline; the saving of initial boundary conditions could be removed. The code in map.h to produce mapping objects which are fully templated and hence potentially faster than the std::function-based mappings in mapping.h has been taken into the code base but is currently unused.
Kay F. Jahnke
kfj-guest at moszumanska.debian.org
Sun Jul 2 09:02:42 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 5002ac10907cd8bddd0bb05ccd1c613f496a38a7
Author: Kay F. Jahnke <kfjahnke at gmail.com>
Date: Fri Apr 28 08:17:52 2017 +0200
separated boundary conditions and mapping codes
I was unhappy with both voundary conditions and mappings sharing the same
enum and being conceptually linked. Now there is a separate set of mapping
codes (in common.h) and the mapping code in mapping.h is based on them.
When creating evaluators from bspline objects, the old automatism of
picking a mapping which 'echoed' the spline's boundary conditions was
removed, the default now is to use MAP_REJECT, which throws
vspline::out_of_bounds for out-of-bounds coordinates - but any of the
other mappings can be used as well by simply passing in the desired
mapping code(s) or an mmap object.
The coundary conditions in bspline objects are now const, and the old
mechanism of changing the spline's bc's in order to get different mapping
behaviour when automatically creating an evaluator from a bspline object
can no longer be used, which is a relief, since it simplifies matters,
especially when shifting the spline; the saving of initial boundary conditions
could be removed.
The code in map.h to produce mapping objects which are fully templated
and hence potentially faster than the std::function-based mappings in
mapping.h has been taken into the code base but is currently unused.
---
brace.h | 15 +-
bspline.h | 60 ++---
common.h | 68 ++---
eval.h | 28 +-
example/eval.cc | 54 +++-
example/gsm.cc | 9 +-
example/roundtrip.cc | 2 +-
example/splinus.cc | 2 +-
filter.h | 78 +++---
map.h | 185 +++++++++++--
mapping.h | 726 ++++++++++++++++-----------------------------------
remap.h | 1 +
12 files changed, 578 insertions(+), 650 deletions(-)
diff --git a/brace.h b/brace.h
index f6997cd..aab113a 100644
--- a/brace.h
+++ b/brace.h
@@ -177,9 +177,9 @@ struct bracer
static int left_brace_size ( int spline_degree , bc_code bc )
{
- if ( bc == REFLECT || bc == SPHERICAL || bc == RTAG )
- return ( spline_degree + 1 ) / 2 ;
- else
+// if ( bc == REFLECT || bc == SPHERICAL )
+// return ( spline_degree + 1 ) / 2 ;
+// else
return spline_degree / 2 ;
}
@@ -373,8 +373,6 @@ struct bracer
}
case NATURAL :
case MIRROR :
- case TAG :
- case TAGP :
{
ls = l0 + 2 ;
rs = r0 - 2 ;
@@ -383,7 +381,6 @@ struct bracer
case CONSTANT :
case SPHERICAL :
case REFLECT :
- case RTAG :
{
ls = l0 + 1 ;
rs = r0 - 1 ;
@@ -415,9 +412,6 @@ struct bracer
case PERIODIC :
case MIRROR :
case REFLECT :
- case TAG :
- case RTAG :
- case TAGP :
{
// with these three bracing modes, we simply copy from source to target
a.bindAt ( axis , lt ) = a.bindAt ( axis , ls ) ;
@@ -467,9 +461,6 @@ struct bracer
case PERIODIC :
case MIRROR :
case REFLECT :
- case TAG :
- case RTAG :
- case TAGP :
{
// with these three bracing modes, we simply copy from source to target
a.bindAt ( axis , rt ) = a.bindAt ( axis , rs ) ;
diff --git a/bspline.h b/bspline.h
index ac18706..cf013e4 100644
--- a/bspline.h
+++ b/bspline.h
@@ -203,15 +203,15 @@ private:
array_type _coeffs ;
prefilter_strategy strategy ;
- bcv_type boundary_conditions ; ///< original boundary conditions, see common.h
public:
+ const bcv_type bcv ; ///< coundary conditions, see common.h
+
view_type container ; ///< view to container array
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 ; ///< 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
@@ -237,7 +237,7 @@ public:
container_shape = braced_shape = core_shape ;
left_brace = right_brace = left_frame = right_frame = shape_type() ;
braced = false ;
- break ;
+ return ;
case BRACED:
// again an implicit prefiltering scheme will be used, but here we add
// a 'brace' to the core data, which makes the resulting bspline object
@@ -249,8 +249,6 @@ public:
right_brace = bracer<view_type>::right_corner ( bcv , spline_degree ) ;
left_frame = left_brace + headroom ;
right_frame = right_brace + headroom ;
- container_shape = core_shape + left_frame + right_frame ;
- braced = true ;
break ;
case EXPLICIT:
// here we prepare for an explicit extrapolation. This requires additional
@@ -264,10 +262,25 @@ public:
right_brace = bracer<view_type>::right_corner ( bcv , spline_degree ) ;
left_frame = left_brace + horizon + headroom ;
right_frame = right_brace + horizon + headroom ;
- container_shape = core_shape + left_frame + right_frame ;
- braced = true ;
break ;
}
+ braced = true ;
+
+ // for odd splines with REFLECT boundary conditions we increase the size of
+ // the left frame, so that access to coordinates in [-0.5:0] will not result
+ // in accessing memory outside the coefficient array. Access to [M-1:M-0.5]
+ // is safe with the right brace size the bracer has returned.
+
+ if ( spline_degree & 1 )
+ {
+ for ( int d = 0 ; d < dimension ; d++ )
+ {
+ if ( bcv[d] == REFLECT || bcv[d] == SPHERICAL )
+ left_frame[d]++ ;
+ }
+ }
+
+ container_shape = core_shape + left_frame + right_frame ;
}
/// this method calculates the size of container needed by a bspline object with
@@ -345,7 +358,6 @@ public:
: core_shape ( _core_shape ) ,
spline_degree ( _spline_degree ) ,
bcv ( _bcv ) ,
- boundary_conditions ( _bcv ) ,
smoothing ( _smoothing ) ,
strategy ( _strategy )
{
@@ -374,6 +386,9 @@ public:
else
horizon = _horizon ; // whatever the user specifies
+// if ( _strategy == BRACED || _strategy == EXPLICIT )
+// headroom ++ ;
+
// first, calculate all the various shapes and sizes used internally
setup_metrics ( headroom ) ;
@@ -415,7 +430,6 @@ public:
bspline ( const bspline& other )
: strategy ( other.strategy ) ,
- boundary_conditions ( other.boundary_conditions ) ,
container ( other.container ) ,
coeffs ( other.coeffs ) ,
core ( other.core ) ,
@@ -493,11 +507,6 @@ public:
void prefilter ( int njobs = default_njobs , ///< intended number of jobs 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
@@ -620,7 +629,7 @@ public:
bspline ( long _core_shape , ///< shape of knot point data
int _spline_degree = 3 , ///< spline degree with reasonable default
- bc_code _bcv = bcv_type ( MIRROR ) , ///< boundary conditions and common default
+ bc_code _bc = MIRROR , ///< boundary conditions and common default
prefilter_strategy _strategy = BRACED , ///< default strategy is the 'implicit' scheme
int _horizon = -1 , ///< width of frame for explicit scheme
view_type _space = view_type() , ///< coefficient storage to 'adopt'
@@ -630,7 +639,7 @@ public:
)
:bspline ( TinyVector < long , 1 > ( _core_shape ) ,
_spline_degree ,
- bcv_type ( _bcv ) ,
+ bcv_type ( _bc ) ,
_strategy ,
_horizon ,
_space ,
@@ -658,14 +667,6 @@ public:
/// This is a quick-shot solution to the problem of scaled-down interpolated
/// results; it may be better in some situations to shift the spline up and
/// evaluate than to apply smoothing to the source data.
- /// 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.
// TODO consider moving the concept of shifting to class evaluator
@@ -675,9 +676,6 @@ 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 ) ;
@@ -701,7 +699,6 @@ public:
d = 0 ;
}
- bcv = bcv_backup ; // restore saved mapping mode
return d ;
}
@@ -711,14 +708,11 @@ public:
{
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:................ " ;
+ osr << "boundary conditions:......... " ;
for ( auto bc : bsp.bcv )
osr << " " << bc_name [ bc ] ;
osr << endl ;
+ 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 ;
diff --git a/common.h b/common.h
index 591bf2c..899bf34 100644
--- a/common.h
+++ b/common.h
@@ -188,27 +188,31 @@ const double rv_out_of_bounds = 0.0 ;
/// the core array. Finally, mapping of coordinates into the defined range uses some of the
/// same codes.
-typedef enum { PERIODIC, ///< periodic boundary conditions / periodic mapping
- NATURAL, ///< natural boundary conditions / uses constant mapping
- MIRROR , ///< mirror on the bounds, both bc and mapping
- REFLECT , ///< reflect, so that f(-1) == f(0), both bc and mapping
- SAFE_BEYOND , ///< assume the input extends beyond [0...M]
- CONSTANT , ///< used for framing, with explicit prefilter scheme, and for mapping
- ZEROPAD , ///< used for boundary condition, bracing
- IGNORE , ///< used for boundary condition, bracing
- IDENTITY , ///< used as solver argument, mostly internal use
- GUESS , ///< used with EXPLICIT scheme to keep margin errors low
- SPHERICAL , ///< intended use for spherical panoramas, needs work
- REJECT , ///< mapping mode, throws out_of_bounds for out-of-bounds coordinates
- REJECTP , ///< mapping mode, like REJECT, but ceiling 1 higher
- LIMIT , ///< mapping mode, maps out-of-bounds coordinates to left/right boundary
- 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
+typedef enum {
+ MIRROR , ///< mirror on the bounds
+ PERIODIC, ///< periodic boundary conditions
+ REFLECT , ///< reflect, so that f(-1) == f(0)
+ NATURAL, ///< natural boundary conditions
+ CONSTANT , ///< used for framing, with explicit prefilter scheme
+ ZEROPAD , ///< used for boundary condition, bracing
+ IGNORE , ///< used for boundary condition, bracing
+ IDENTITY , ///< used as solver argument, mostly internal use
+ GUESS , ///< used with EXPLICIT scheme to keep margin errors low
+ SPHERICAL , ///< use for spherical panoramas y axis
+// SAFE_BEYOND , ///< assume the input extends beyond [0...M]
} bc_code;
+typedef enum {
+ MAP_MIRROR , ///< mirror on the bounds
+ MAP_PERIODIC, ///< periodic mapping
+ MAP_REFLECT , ///< mirror between bounds
+ MAP_LIMIT , ///< maps out-of-bounds coordinates to bounds
+ MAP_REJECT , ///< throw out_of_bounds for out-of-bounds coordinates
+ MAP_RAW , ///< process coordinates unchecked, may crash!
+ MAP_CONSTANT , ///< assign constant value to out-of-bounds coordinates
+ MAP_FUNCTION ///< applies std::function to coordinates
+} map_code;
+
/// This enumeration is used by the convenience class 'bspline' to determine the prefiltering
/// scheme to be used.
@@ -222,25 +226,29 @@ typedef enum { UNBRACED , ///< implicit scheme, no bracing applied
const std::string bc_name[] =
{
- "PERIODIC",
- "NATURAL",
"MIRROR" ,
+ "PERIODIC",
"REFLECT" ,
- "SAFE_BEYOND" ,
+ "NATURAL",
"CONSTANT" ,
"ZEROPAD" ,
"IGNORE" ,
"IDENTITY" ,
"GUESS" ,
"SPHERICAL" ,
- "REJECT" ,
- "REJECTP" ,
- "LIMIT" ,
- "LIMITP" ,
- "TAG" ,
- "TAGP" ,
- "RTAG" ,
- "RAW"
+// "SAFE_BEYOND" , // currently unused
+} ;
+
+const std::string mc_name[] =
+{
+ "MAP_MIRROR" ,
+ "MAP_PERIODIC" ,
+ "MAP_REFLECT" ,
+ "MAP_LIMIT" ,
+ "MAP_REJECT" ,
+ "MAP_RAW" ,
+ "MAP_CONSTANT" ,
+ "MAP_FUNCTION"
} ;
} ; // end of namespace vspline
diff --git a/eval.h b/eval.h
index 5cf7915..4ea819d 100644
--- a/eval.h
+++ b/eval.h
@@ -635,17 +635,21 @@ public:
}
} ;
- /// simplified constructor from a bspline object
+ /// simplified constructors from a bspline object
/// when using the higher-level interface to vspline's facilities - via class bspline -
/// class bspline objects already have a fair amount of metadata which are perfectly
/// sufficient to create a suitable evaluator object. If evaluation of the spline's value
/// (no derivative spec) is wanted, all this constructor takes is a const reference
/// to a bspline object, which is about as simple as it can get:
+ typedef vigra::TinyVector < vspline::map_code , dimension > mcv_type ;
+
evaluator ( const bspline < value_type , dimension > & bspl ,
- derivative_spec_type _derivative = derivative_spec_type ( 0 ) )
+ derivative_spec_type _derivative ,
+ mcv_type mcv
+ )
: evaluator ( bspl.coeffs ,
- nd_mapping_type ( bspl.bcv , bspl.spline_degree , bspl.core_shape ) ,
+ nd_mapping_type ( bspl.bcv , mcv , bspl.spline_degree , bspl.core_shape ) ,
bspl.spline_degree ,
_derivative )
{
@@ -653,6 +657,24 @@ public:
throw not_supported ( "for spline degree > 1: evaluation needs braced coefficients" ) ;
} ;
+ /// constructor variant accepting a common derivative specification and mapping code
+ /// for all axes of the spline
+
+ evaluator ( const bspline < value_type , dimension > & bspl ,
+ int _derivative = 0 ,
+ vspline::map_code mc = vspline::MAP_REJECT
+ )
+ : evaluator ( bspl , derivative_spec_type ( _derivative ) , mcv_type ( mc ) )
+ { } ;
+
+ /// constructor variant for derivative-0 splines accepting a common mapping code
+ /// for all axes of the spline
+
+ evaluator ( const bspline < value_type , dimension > & bspl ,
+ vspline::map_code mc )
+ : evaluator ( bspl , derivative_spec_type ( 0 ) , mcv_type ( mc ) )
+ { } ;
+
const nd_mapping_type& get_mapping() const
{
return mmap ;
diff --git a/example/eval.cc b/example/eval.cc
index 393505a..1df5ad1 100644
--- a/example/eval.cc
+++ b/example/eval.cc
@@ -61,9 +61,9 @@ int main ( int argc , char * argv[] )
{
cout << "choose boundary condition" << endl ;
cout << "1) MIRROR" << endl ;
- cout << "2) REFLECT" << endl ;
- cout << "3) NATURAL" << endl ;
- cout << "4) PERIODIC" << endl ;
+ cout << "2) PERIODIC" << endl ;
+ cout << "3) REFLECT" << endl ;
+ cout << "4) NATURAL" << endl ;
cin >> bci ;
}
@@ -73,18 +73,58 @@ int main ( int argc , char * argv[] )
bc = MIRROR ;
break ;
case 2 :
- bc = REFLECT ;
+ bc = PERIODIC ;
break ;
case 3 :
- bc = NATURAL ;
+ bc = REFLECT ;
break ;
case 4 :
- bc = PERIODIC ;
+ bc = NATURAL ;
break ;
}
// put the BC code into a TinyVector
TinyVector < bc_code , 1 > bcv ( bc ) ;
+ int mci = -1 ;
+ map_code mc ;
+
+ while ( mci < 1 || mci > 6 )
+ {
+ cout << "choose mapping mode" << endl ;
+ cout << "1) MAP_MIRROR" << endl ;
+ cout << "2) MAP_PERIODIC" << endl ;
+ cout << "3) MAP_REFLECT" << endl ;
+ cout << "4) MAP_LIMIT" << endl ;
+ cout << "5) MAP_REJECT" << endl ;
+ cout << "6) MAP_RAW" << endl ;
+ cin >> mci ;
+ }
+
+ switch ( mci )
+ {
+ case 1 :
+ mc = MAP_MIRROR ;
+ break ;
+ case 2 :
+ mc = MAP_PERIODIC ;
+ break ;
+ case 3 :
+ mc = MAP_REFLECT ;
+ break ;
+ case 4 :
+ mc = MAP_LIMIT ;
+ break ;
+ case 5 :
+ mc = MAP_REJECT ;
+ break ;
+ case 6 :
+ mc = MAP_RAW ;
+ break ;
+ }
+ // put the mapping code into a TinyVector
+ TinyVector < map_code , 1 > mcv ( mc ) ;
+
+ TinyVector < int , 1 > deriv_spec ( 0 ) ;
// obtain knot point values
double v ;
@@ -117,7 +157,7 @@ int main ( int argc , char * argv[] )
// fix the type for the evaluator and create it
typedef evaluator < double , double > eval_type ;
- eval_type ev ( bsp ) ;
+ eval_type ev ( bsp , deriv_spec , mcv ) ;
auto map = ev.get_mapping() ;
int ic ;
double rc ;
diff --git a/example/gsm.cc b/example/gsm.cc
index 62e5e82..79a1c84 100644
--- a/example/gsm.cc
+++ b/example/gsm.cc
@@ -90,13 +90,8 @@ int main ( int argc , char * argv[] )
// with the desired derivative degree for each dimension. Here we want the
// first derivative in x and y direction:
- vigra::TinyVector < float , 2 > dx1_spec { 1 , 0 } ;
- vigra::TinyVector < float , 2 > dy1_spec { 0 , 1 } ;
-
- // now we can construct the evaluators
-
- ev_type xev ( bspl , dx1_spec ) ;
- ev_type yev ( bspl , dy1_spec ) ;
+ ev_type xev ( bspl , { 1 , 0 } ) ;
+ ev_type yev ( bspl , { 0 , 1 } ) ;
// this is where the result should go:
target_type target ( imageInfo.shape() ) ;
diff --git a/example/roundtrip.cc b/example/roundtrip.cc
index 5440e95..e18b29b 100644
--- a/example/roundtrip.cc
+++ b/example/roundtrip.cc
@@ -310,7 +310,7 @@ void process_image ( char * name )
for ( int b = 0 ; b < 4 ; b++ )
{
vspline::bc_code bc = bcs[b] ;
- for ( int spline_degree = 0 ; spline_degree < 8 ; spline_degree++ )
+ for ( int spline_degree = 0 ; spline_degree < 4 ; spline_degree++ )
{
#ifdef USE_VC
cout << "testing bc code " << vspline::bc_name[bc]
diff --git a/example/splinus.cc b/example/splinus.cc
index b84abf2..555eaf2 100644
--- a/example/splinus.cc
+++ b/example/splinus.cc
@@ -79,7 +79,7 @@ int main ( int argc , char * argv[] )
vspline::evaluator < double , // evaluator taking double coordinates
double // and producing double results
- > ev ( bsp ) ; // from the bspline object we just made
+ > ev ( bsp , 0 ) ; // from the bspline object we just made
while ( true )
{
diff --git a/filter.h b/filter.h
index b0f7e62..5f45b89 100644
--- a/filter.h
+++ b/filter.h
@@ -461,44 +461,44 @@ value_type iacc_periodic ( out_iter c , int k )
return Sum ;
}
-/// icc_safe_beyond and iacc_safe_beyond are for situations where c[0] ... c[M-1]
-/// merely define a section of a larger array of data, so that it's safe to access
-/// horizon[k] values beyond either end of c[0] ... c[M-1].
-/// currently unused, because it doesn't work with buffering
-
-template < class IT >
-value_type icc_safe_beyond ( IT c , int k )
-{
- value_type p = value_type ( pole[k] ) ;
- int z = horizon[k] ;
- value_type icc = c [ - z ] ;
-
- for ( int i = - z + 1 ; i <= 0 ; i++ )
- {
- icc += p * icc + c[i] ;
- }
-
- return icc ;
-}
-
-value_type iacc_safe_beyond ( out_iter c , int k )
-{
- value_type p = value_type ( pole[k] ) ;
- int z = horizon[k] ;
- value_type iacc = c [ M - 1 + z ] ;
-
- for ( int i = z - 1 ; i >= 0 ; i-- )
- iacc = p * ( iacc - c [ M - 1 + i ] ) ;
- return iacc ;
-}
+// /// icc_safe_beyond and iacc_safe_beyond are for situations where c[0] ... c[M-1]
+// /// merely define a section of a larger array of data, so that it's safe to access
+// /// horizon[k] values beyond either end of c[0] ... c[M-1].
+// /// currently unused, because it doesn't work with buffering
+//
+// template < class IT >
+// value_type icc_safe_beyond ( IT c , int k )
+// {
+// value_type p = value_type ( pole[k] ) ;
+// int z = horizon[k] ;
+// value_type icc = c [ - z ] ;
+//
+// for ( int i = - z + 1 ; i <= 0 ; i++ )
+// {
+// icc += p * icc + c[i] ;
+// }
+//
+// return icc ;
+// }
+//
+// value_type iacc_safe_beyond ( out_iter c , int k )
+// {
+// value_type p = value_type ( pole[k] ) ;
+// int z = horizon[k] ;
+// value_type iacc = c [ M - 1 + z ] ;
+//
+// for ( int i = z - 1 ; i >= 0 ; i-- )
+// iacc = p * ( iacc - c [ M - 1 + i ] ) ;
+// return iacc ;
+// }
// guess the initial coefficient. This tries to minimize the effect
// of starting out with a hard discontinuity as it occurs with zero-padding,
// while at the same time requiring little arithmetic effort
-// TODO needs some thought
// for the forward filter, we guess an extrapolation of the signal to the left
// repeating c[0] indefinitely, which is cheap to compute:
+
template < class IT >
value_type icc_guess ( IT c , int k )
{
@@ -506,6 +506,7 @@ value_type icc_guess ( IT c , int k )
}
// for the backward filter, we assume mirror BC, which is also cheap to compute:
+
value_type iacc_guess ( out_iter c , int k )
{
return iacc_mirror ( c , k ) ;
@@ -767,12 +768,6 @@ filter ( int _M , ///< number of input/output elements (DataLength
_p_icc2 = & filter_type::icc_reflect<out_iter> ;
_p_iacc = & filter_type::iacc_reflect ;
}
- else if ( bc == SAFE_BEYOND )
- {
- _p_icc1 = & filter_type::icc_safe_beyond<in_iter> ;
- _p_icc2 = & filter_type::icc_safe_beyond<out_iter> ;
- _p_iacc = & filter_type::iacc_safe_beyond ;
- }
else if ( bc == ZEROPAD || bc == IGNORE )
{
_p_icc1 = & filter_type::icc_identity<in_iter> ;
@@ -783,12 +778,19 @@ filter ( int _M , ///< number of input/output elements (DataLength
{
_p_solve = & filter_type::solve_identity ;
}
- else if ( bc == GUESS ) // currently unused
+ else if ( bc == GUESS )
{
_p_icc1 = & filter_type::icc_guess<in_iter> ;
_p_icc2 = & filter_type::icc_guess<out_iter> ;
_p_iacc = & filter_type::iacc_guess ;
}
+// /// currently unused, because it doesn't work with buffering
+// else if ( bc == SAFE_BEYOND )
+// {
+// _p_icc1 = & filter_type::icc_safe_beyond<in_iter> ;
+// _p_icc2 = & filter_type::icc_safe_beyond<out_iter> ;
+// _p_iacc = & filter_type::iacc_safe_beyond ;
+// }
else
{
std::cout << "boundary condition " << bc << " not supported by vspline::filter" << std::endl ;
diff --git a/map.h b/map.h
index b403d77..b5047e6 100644
--- a/map.h
+++ b/map.h
@@ -55,7 +55,7 @@
the gate_types and applies them to the components in turn.
The final mapper object is a functor which converts an arbitrary incoming
- coordinate into an treated coordinate (or, for REJECT mode, throws an
+ coordinate into a 'treated' coordinate (or, for REJECT mode, may throw an
out_of_bounds exception). If the treatment succeeded to produce a
suitable in-range coordinate, accessing the evaluator with it is safe.
@@ -64,6 +64,19 @@
'unmirror' on the bounds or to 'unperiodize'. Finally there is a version
taking two std::functions, one for single values and one for simdized types,
which applies the appropriate functor to incoming coordinates.
+
+ This code more or less duplicates the capabilities of the 'mmap' object
+ class evaluator uses: this object can be constructed with the desired mapping
+ modes at run-time and is therefore more flexible, and at times it's code is less
+ involved, since it is less general. But using mapping ocde fixed at compile-time,
+ as is done in this body of code, is potentially faster.
+
+ On my system I noticed little difference when using clang++, but compiled with
+ g++ using the templated code in this file performed noticeably (ca. 10%) faster.
+
+ gate_type objects are derived from vspline::unary_functor, so they fit in well
+ with other code in vspline and can easily be wrapped in other unary_functor
+ objects, or used stand-alone.
*/
#ifndef VSPLINE_MAP_H
@@ -79,7 +92,7 @@ namespace vspline
/// actually do any work.
template < typename rc_type ,
- vspline::bc_code mode ,
+ vspline::map_code mode ,
int vsize >
struct gate_type
{
@@ -90,11 +103,13 @@ struct gate_type
/// the most basic gate passes it's input to it's output unmodified.
/// Like all gate_types, this gate_type inherits from a vspline::unary_functor,
/// with input and output types equal. this specialization needs no
-/// constructor, since it's independent of a range.
+/// constructor, since it's independent of a range, but I have included a
+/// constructor with two default arguments to make it constructor-signature
+/// compatible with the other ones.
template < typename _rc_type ,
int _vsize >
-struct gate_type < _rc_type , vspline::RAW , _vsize >
+struct gate_type < _rc_type , vspline::MAP_RAW , _vsize >
: public vspline::unary_functor < _rc_type , _rc_type , _vsize >
{
typedef _rc_type rc_type ;
@@ -107,6 +122,10 @@ struct gate_type < _rc_type , vspline::RAW , _vsize >
using_unary_functor_types(base_type) ;
+ gate_type ( rc_type _lower = 0 ,
+ rc_type _upper = 0 )
+ {} ;
+
// implement vspline::unary_functor's two pure virtual methods:
void eval ( const in_type & c ,
@@ -114,20 +133,24 @@ struct gate_type < _rc_type , vspline::RAW , _vsize >
{
result = c ;
}
-
+
+#ifdef USE_VC
+
void eval ( const in_v & c ,
out_v & result ) const
{
result = c ;
}
+#endif
+
} ;
/// gate with REJECT throws vspline::out_of_bounds for invalid coordinates
template < typename _rc_type ,
int _vsize >
-struct gate_type < _rc_type , vspline::REJECT , _vsize >
+struct gate_type < _rc_type , vspline::MAP_REJECT , _vsize >
: public vspline::unary_functor < _rc_type , _rc_type , _vsize >
{
typedef _rc_type rc_type ;
@@ -153,6 +176,8 @@ struct gate_type < _rc_type , vspline::REJECT , _vsize >
result = c ;
}
+#ifdef USE_VC
+
void eval ( const in_v & c ,
out_v & result ) const
{
@@ -161,13 +186,15 @@ struct gate_type < _rc_type , vspline::REJECT , _vsize >
result = c ;
}
+#endif
+
} ;
/// gate with LIMIT clamps out-of-bounds values
template < typename _rc_type ,
int _vsize >
-struct gate_type < _rc_type , vspline::LIMIT , _vsize >
+struct gate_type < _rc_type , vspline::MAP_LIMIT , _vsize >
: public vspline::unary_functor < _rc_type , _rc_type , _vsize >
{
typedef _rc_type rc_type ;
@@ -194,6 +221,8 @@ struct gate_type < _rc_type , vspline::LIMIT , _vsize >
result = c ;
}
+#ifdef USE_VC
+
void eval ( const in_v & c ,
out_v & result ) const
{
@@ -203,13 +232,17 @@ struct gate_type < _rc_type , vspline::LIMIT , _vsize >
result[0] = cc ;
}
+#endif
+
} ;
-/// constant gate assigns 'fix' to out-of-bounds coordinates
+/// constant gate assigns 'fix' to out-of-bounds coordinates; by default
+/// it will return 0 for out-of-bounds input, but the value can be chosen
+/// by passing a third argument to the constructor.
template < typename _rc_type ,
int _vsize >
-struct gate_type < _rc_type , vspline::CONSTANT , _vsize >
+struct gate_type < _rc_type , vspline::MAP_CONSTANT , _vsize >
: public vspline::unary_functor < _rc_type , _rc_type , _vsize >
{
typedef _rc_type rc_type ;
@@ -222,7 +255,7 @@ struct gate_type < _rc_type , vspline::CONSTANT , _vsize >
gate_type ( rc_type _lower ,
rc_type _upper ,
- rc_type _fix
+ rc_type _fix = rc_type(0)
)
: lower ( _lower ) ,
upper ( _upper ) ,
@@ -240,6 +273,8 @@ struct gate_type < _rc_type , vspline::CONSTANT , _vsize >
result = c ;
}
+#ifdef USE_VC
+
void eval ( const in_v & c ,
out_v & result ) const
{
@@ -249,16 +284,18 @@ struct gate_type < _rc_type , vspline::CONSTANT , _vsize >
result[0] = cc ;
}
+#endif
+
} ;
-/// gate with mirror moves coordinates into the range. It places the
+/// gate with mirror 'folds' coordinates into the range. It places the
/// result value so that mirroring it on both lower and upper (which
/// produces an infinite number of mirror images) will produce one
/// mirror image coinciding with the input.
template < typename _rc_type ,
int _vsize >
-struct gate_type < _rc_type , vspline::MIRROR , _vsize >
+struct gate_type < _rc_type , vspline::MAP_MIRROR , _vsize >
: public vspline::unary_functor < _rc_type , _rc_type , _vsize >
{
typedef _rc_type rc_type ;
@@ -291,6 +328,8 @@ struct gate_type < _rc_type , vspline::MIRROR , _vsize >
result = cc ;
}
+#ifdef USE_VC
+
void eval ( const in_v & c ,
out_v & result ) const
{
@@ -312,11 +351,81 @@ struct gate_type < _rc_type , vspline::MIRROR , _vsize >
result[0] = cc + lower ;
}
+#endif
+
} ;
+// same as MIRROR:
+
template < typename _rc_type ,
int _vsize >
-struct gate_type < _rc_type , vspline::PERIODIC , _vsize >
+struct gate_type < _rc_type , vspline::MAP_REFLECT , _vsize >
+: public vspline::unary_functor < _rc_type , _rc_type , _vsize >
+{
+ typedef _rc_type rc_type ;
+ typedef vspline::unary_functor < rc_type , rc_type , _vsize > base_type ;
+ using_unary_functor_types(base_type) ;
+
+ const rc_type lower ;
+ const rc_type upper ;
+
+ gate_type ( rc_type _lower ,
+ rc_type _upper )
+ : lower ( _lower ) ,
+ upper ( _upper )
+ { } ;
+
+ void eval ( const in_type & c ,
+ out_type & result ) const
+ {
+ in_type cc = c ;
+ if ( cc < lower )
+ cc = 2 * lower - cc ;
+ if ( cc > upper )
+ {
+ cc -= lower ;
+ cc = std::fmod ( cc , 2 * ( upper - lower ) ) ;
+ cc += lower ;
+ if ( cc > upper )
+ cc = 2 * upper - cc ;
+ }
+ result = cc ;
+ }
+
+#ifdef USE_VC
+
+ void eval ( const in_v & c ,
+ out_v & result ) const
+ {
+ out_ele_v cc ;
+
+ cc = c[0] - lower ;
+ auto w = upper - lower ;
+
+ cc = abs ( cc ) ; // left mirror, v is now >= 0
+
+ if ( any_of ( cc > w ) )
+ {
+ cc = v_fmod ( cc , 2 * w ) ; // map to one full period
+ cc -= w ; // ccenter
+ cc = abs ( cc ) ; // map to half period
+ cc = w - cc ; // flip
+ }
+
+ result[0] = cc + lower ;
+ }
+
+#endif
+
+} ;
+
+/// the periodic mapping also folds the incoming value into the allowed range.
+/// The resulting value will be ( N * period ) from the input value and inside
+/// the range, period being upper - lower.
+
+template < typename _rc_type ,
+ int _vsize >
+struct gate_type < _rc_type , vspline::MAP_PERIODIC , _vsize >
: public vspline::unary_functor < _rc_type , _rc_type , _vsize >
{
typedef _rc_type rc_type ;
@@ -345,6 +454,8 @@ struct gate_type < _rc_type , vspline::PERIODIC , _vsize >
result = cc + lower ;
}
+#ifdef USE_VC
+
void eval ( const in_v & c ,
out_v & result ) const
{
@@ -362,15 +473,18 @@ struct gate_type < _rc_type , vspline::PERIODIC , _vsize >
result[0] = cc + lower ;
}
+#endif
+
} ;
-/// currently abusing GUESS to delegate to a pair of std::functions
+/// mode MAP_FUNCTION is used to delegate to a pair of std::functions
/// this allows user code to introduce arbitrary gate functions when
-/// the provided ready-made ones aren't sufficient.
+/// the provided ready-made ones aren't sufficient. This is a mode which
+/// evaluator's mmap doesn't use.
template < typename _rc_type ,
int _vsize >
-struct gate_type < _rc_type , vspline::GUESS , _vsize >
+struct gate_type < _rc_type , vspline::MAP_FUNCTION , _vsize >
: public vspline::unary_functor < _rc_type , _rc_type , _vsize >
{
typedef _rc_type rc_type ;
@@ -378,27 +492,37 @@ struct gate_type < _rc_type , vspline::GUESS , _vsize >
using_unary_functor_types(base_type) ;
typedef std::function < out_type ( in_type ) > func_type ;
- typedef std::function < out_ele_v ( in_ele_v ) > vfunc_type ;
-
const func_type f ;
+
+#ifdef USE_VC
+
+ typedef std::function < out_ele_v ( in_ele_v ) > vfunc_type ;
const vfunc_type vf ;
gate_type ( func_type _f , vfunc_type _vf )
: f ( _f ) , vf ( _vf )
{ } ;
- void eval ( const in_type & c ,
- out_type & result ) const
- {
- result = f ( c ) ;
- }
-
void eval ( const in_v & c ,
out_v & result ) const
{
result[0] = vf ( c[0] ) ;
}
+#else
+
+ gate_type ( func_type _f )
+ : f ( _f )
+ { } ;
+
+#endif
+
+ void eval ( const in_type & c ,
+ out_type & result ) const
+ {
+ result = f ( c ) ;
+ }
+
} ;
// /// struct mapper1 applies a gate to it's incoming coordinate,
@@ -521,6 +645,19 @@ struct mapper
}
} ;
+/// struct raw_mapper just passes it's input through to it's output.
+/// This is a shorter version of using a set of gate_types with RAW mapping.
+
+struct raw_mapper
+{
+ template < typename nd_rc_type >
+ void operator() ( const nd_rc_type & in ,
+ nd_rc_type & out )
+ {
+ out = in ;
+ }
+} ;
+
} ; // namespace vspline
#endif // #ifndef VSPLINE_MAP_H
diff --git a/mapping.h b/mapping.h
index aeb3e7b..07ae82b 100644
--- a/mapping.h
+++ b/mapping.h
@@ -41,67 +41,68 @@
\brief code to handle the processing of incoming real coordinates
- One might argue that incoming coordinates should be required to be inside
- the defined range of the coefficients. While this may or may not be enforced
- or checked with the code given here, I also provide 'mapping' routines which
- correspond to some common boundary conditions used for a spline. With these
- mappings, incoming coordinates are 'folded into' the defined range of the spline
- by applying the relevant transformation, like mirroring or applying a modulo
- operation for periodic BCs.
-
- While the mappings are central to the evaluation process, there isn't usually a
- need to directly handle them: If a bspline object is used to contain a coefficient
- array and it's metadata, the appropriate mapping can be derived from it, and the
- mapping doesn't have to be made explicit.
-
- Note how we use the same enumeration from common.h to codify boundary
- conditions and bracings, because often they correspond, and, more specifically,
- those boundary condition codes used to create coefficient arrays (bspline objects)
- with implicit prefiltering schemes have corresponding mappings.
-
- for BC codes MIRROR, REFLECT and PERIODIC this means that incoming coordinates will
- be mapped so that they 'land' at a coordinate inside the defined range. But note that
- for NATURAL BCs, there is no coordinate inside the defined range that would produce
- the extrapolated value used for natural boundary conditions, since this value is not
- 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 there is some collateral code to provide fmod and modf for vectors.
+ mapping.h provides code to make sure that incoming coordinates are inside the
+ defined range, or that out-of-bounds coordinates are noticed and rejected.
- 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, 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:
- If the boundary conditions lend themselves to it, such coordinates are mapped onto
- coordinates inside the defined range ('folded in', as I say). If the boundary conditions
- don't work this way (like with 'natural' boundary conditions, where the extrapolated signal
- is the result of an arithmetic operation), we do something arbitrary to the outlying
- coordinate (like clamp it to the neares extremal value) in the conviction that a user
- requesting such boundary conditions is aware of the problem.
-
- vectorized operation with long ints is not yet available, so for now I'll have to limit
- the vectorized code to 32bit ints. This will still suffice for most cases, but if the arrays
- become really large it'll eventually break.
+ It aso provides code for the 'splitting' of coordinates, which separates
+ coordinates into an integral part and a small remainder, which is what is
+ needed for the evaluation of a b-spline.
+
+ vspline::evaluator objects created with the template argument 'evenness'
+ set to 0 or 1 directly use a 'raw' mapping, which splits incoming coordinates
+ unchecked, leaving the user to ensure that all coordinates are in range.
+ This is achieved by a direct call to a functor of type odd_mapping_raw or
+ even_mapping_raw (see below).
+
+ If vspline::evaluator objects are created with 'evenness' set to any other
+ value - like -1, which is the default - they either accept a parameter
+ 'mmap' which constitutes a functor which is applied to incoming coordinates,
+ or (when constructed from a b-spline object) a vigra::TinyVector of
+ enum map_code (see common.h). In both cases, an nd_mapping object (see below)
+ is used; in the first case this object is passed in, in the second it is
+ built to the specification found in the TinyVector of map_codes.
+
+ The nd_mapping object performs it's task using one or several of the mapping
+ objects in this header. The mapping(s) are incorporated into the nd_mapping
+ object as std::functions, which is possible since they all share the same
+ signature, and allows to construct a uniform type even if the mappings differ.
+ The downside of this method is that applying the nd_mapping is slightly slower
+ than producing distinct mappings templated by the mapping mode (see the code in
+ map.h which does just that)
+
+ So what is gained from using the nd_mapping object? Since it offers a single type
+ for both odd and even splines with arbitrary combinations of mapping modes, it
+ is well suited to provide code for situations where the desired mapping has to
+ be provided - at least partly - automatically. If maximum performance is needed
+ for time-critical applications, evaluation with 'evenness' set to 0 or 1 and
+ some user-provided coordinate treatment can be used, otherwise some default
+ or simply expressed run-time selectable parameter can be used to set up coordinate
+ treatment which is almost as efficient. In fact my trials suggest that when using
+ clang++, using nd_mapping object is pretty much as fast as using handcrafted
+ code - only with g++ I notice a significant but slight (~10%) speed difference.
+
+ Another nice feature of using nd_mapping objects is that they determine the
+ bounds which are adequate to a b-spline with given boundary conditions automatically.
+ The constructor of class nd_mapping takes both a TinyVector of boundary condition
+ codes and a shape and sets the bounds of the mappings accordingly. this is important
+ since, for example, periodic splines have a wider defined range than nonperiodic
+ ones.
+
+ So, all in all, nd_mappings are helpful, especially when things are kept simple
+ and evaluators are constructed straight from bsplne objects. constructing an
+ evaluator can be as simple as
+
+ evaluator_type ev ( bspline ) ;
+
+ which will pick MAP_REJECT mapping mode by default, throwing an exception for
+ out-of-bounds coordinates, or, slightly more involved
+
+ evaluator_type ev ( bspline , vspline::MAP_MIRROR ) ;
+
+ To use coordinate folding into the range which works as if the coordinates were mirrored
+ on the bounds until they fall into the range, effectively mapping an infinite range of
+ coordinates into the spline's defined range.
- I hope that by the time this becomes an issue, 64bit int SIMD operation becomes available,
- until then, if it has to be used, only the non-vectorized variant can be used for long ints,
- or (slow) workaround code for vectors of long values has to be provided (which I did out of
- curiosity: it works, but it makes the vectorized code slower than the scalar)
- Anyway I'd advocate splitting very large coefficient arrays into tiles. Since the support
- of the evaluation routine is only as large as the spline's order, duplicated strips of
- coefficients in the margins would take up a comparatively small part of the tile and the
- 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
@@ -119,42 +120,42 @@ using namespace vigra ;
#ifdef USE_VC
-/// since Vc doesn't offer a vectorized modf function, we have to code it. We make the effort not
-/// to simply iterate over the vector's components but write 'proper' vector code to make this
-/// operation as efficient as possible.
-
-template <typename rc_v>
-rc_v v_modf ( const rc_v& source ,
- rc_v * const iptr )
-{
- typedef typename rc_v::mask_type mask_type ;
- typedef typename rc_v::EntryType single_type ;
-
- mask_type negative ( source < single_type(0) ) ; // Vc::isnegative ( source ) ;
- rc_v help ( source ) ;
-
- // we treat the case that any of the incoming vector components is negative separately
- // to avoid the corresponding code altogether in the - hopefully - most common case
- // where none of the incoming values are in fact negative.
-
- if ( any_of ( negative ) )
- {
- help(negative) = -source ;
- (*iptr) = floor ( help ) ;
- help -= (*iptr) ;
- (*iptr)(negative) = -(*iptr) ;
- help(negative) = -help ;
- return help ;
- }
- else
- {
- // for all positive components, the operation is trivial:
- (*iptr) = floor ( source ) ;
- return ( help - *iptr ) ;
- }
-}
+// /// since Vc doesn't offer a vectorized modf function, we have to code it. We make the effort not
+// /// to simply iterate over the vector's components but write 'proper' vector code to make this
+// /// operation as efficient as possible.
+//
+// template <typename rc_v>
+// rc_v v_modf ( const rc_v& source ,
+// rc_v * const iptr )
+// {
+// typedef typename rc_v::mask_type mask_type ;
+// typedef typename rc_v::EntryType single_type ;
+//
+// mask_type negative ( source < single_type(0) ) ; // Vc::isnegative ( source ) ;
+// rc_v help ( source ) ;
+//
+// // we treat the case that any of the incoming vector components is negative separately
+// // to avoid the corresponding code altogether in the - hopefully - most common case
+// // where none of the incoming values are in fact negative.
+//
+// if ( any_of ( negative ) )
+// {
+// help(negative) = -source ;
+// (*iptr) = floor ( help ) ;
+// help -= (*iptr) ;
+// (*iptr)(negative) = -(*iptr) ;
+// help(negative) = -help ;
+// return help ;
+// }
+// else
+// {
+// // for all positive components, the operation is trivial:
+// (*iptr) = floor ( source ) ;
+// return ( help - *iptr ) ;
+// }
+// }
-/// ditto, for fmod.
+/// vectorized fmod function
template <typename rc_v>
rc_v v_fmod ( const rc_v& lhs ,
@@ -164,10 +165,6 @@ rc_v v_fmod ( const rc_v& lhs ,
ic_v iv ( lhs / rhs ) ;
return lhs - rhs * rc_v ( iv ) ;
-// or, to always get a positive result:
-// auto result = lhs - rhs * rc_v ( iv ) ;
-// result ( result < 0 ) = result + rhs ;
-// return result ;
}
#endif // USE_VC
@@ -199,8 +196,13 @@ rc_v v_fmod ( const rc_v& lhs ,
///
/// 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
-/// the program may crash.
+/// silently assumed that the coordinates are in range. If this is not true, the results may
+/// be false or the program may even crash.
+///
+/// Contrary to my initial implementation which used modf() to perform the split, I now use
+/// floor(). This ensures that the split is performed correctly even for negative v, which
+/// can occur for example with REFLECT boundary conditions and can be in-range, since
+/// splines done with REFLECT have wider headroom.
template < typename ic_t , typename rc_t , int vsize = 1 >
class odd_mapping_raw
@@ -209,9 +211,12 @@ public:
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[
- iv = ic_t ( fl_i ) ; // set integer part from float representing it
+ rc_t fl_i = std::floor ( v ) ;
+ fv = v - fl_i ;
+ iv = ic_t ( fl_i ) ;
+// 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
}
#ifdef USE_VC
@@ -223,18 +228,21 @@ public:
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[
- iv = ic_v ( fl_i ) ; // set integer part from float representing it
+ rc_v fl_i = floor ( v ) ;
+ fv = v - fl_i ;
+ iv = ic_v ( fl_i ) ;
+// 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
}
#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].
+/// the raw mapping for even splines produces the integral part by rounding
+/// the incoming coordinate. The fractional part is obtained by subtracting
+/// the rounded value from the incoming coordinate, resulting in a fractional
+/// part in [-0.5,0.5], which is precisely what is needed for even splines.
template < typename ic_t , typename rc_t , int vsize = 1 >
class even_mapping_raw
@@ -243,9 +251,13 @@ public:
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 ) ;
+ rc_t fl_i = std::round ( v ) ;
+ fv = v - fl_i ;
iv = ic_t ( fl_i ) ;
+// previously I used:
+// rc_t fl_i ;
+// fv = modf ( v + rc_t ( 0.5 ) , &fl_i ) - rc_t ( 0.5 ) ;
+// iv = ic_t ( fl_i ) ;
}
#ifdef USE_VC
@@ -257,11 +269,15 @@ public:
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 ) ;
- fv -= rc_t ( 0.5 ) ;
- iv = ic_v ( fl_i ) ;
+ rc_v fl_i = round ( v ) ;
+ fv = v - fl_i ;
+ iv = ic_v ( fl_i ) ;
+// previously I used:
+// rc_v fl_i ;
+// v += rc_t ( 0.5 ) ;
+// fv = v_modf ( v , &fl_i ) ;
+// fv -= rc_t ( 0.5 ) ;
+// iv = ic_v ( fl_i ) ;
}
#endif
@@ -471,198 +487,6 @@ public:
#endif
} ;
-/// with mapping mode TAG, out-of-bounds incoming coordinate values don't result in
-/// an exception, but offending values are 'tagged' by the remainder part being set to
-/// rc_out_of_bounds, which is a value which can't be valid. The calling code is
-/// responsible for taking appropriate action.
-/// If the calling code fails to treat the 'impossible' value appropriately, this
-/// will result in a wrong result but not in an access to an invalid memory location,
-/// which might be the result of passing some impossible value for the int part.
-
-// TODO: there is a fundamental problem here: if the data are, say, periodic, we need
-// a different ceiling to when they are, say, mirrored. There have to be different
-// mappings for these scenarios, because if we use the ceiling for mirror data (M-1)
-// we have a 1-pixel gap when the data are in fact periodic, but if we use M as the
-// ceiling, this would let out-of-bounds acesses at the right margin pass for mirror
-// data. And this applies to REJECT and RAW as well.
-// For now I provide the right behaviour via a variant bc code 'TAGP' which feeds
-// a value of M to the constructor which is 1 larger than that for plain TAG.
-
-template < typename ic_t , typename rc_t , int vsize = 1 >
-class odd_mapping_tag
-{
- const rc_t _ceiling ;
-
-public:
-
- odd_mapping_tag ( int M )
- : _ceiling ( M - 1 )
- { } ;
-
- void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
- {
- 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 ;
- }
-
- odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
- }
-
-#ifdef USE_VC
-
- typedef typename vector_traits < ic_t , vsize > :: ele_v ic_v ;
- typedef typename vector_traits < rc_t , vsize > :: ele_v rc_v ;
-
- /// this is the operator() for vectorized operation
-
- 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 ;
-
- // 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, tag them
- iv ( typename ic_v::mask_type ( mask ) ) = ic_t ( ic_out_of_bounds ) ;
- fv ( mask ) = rc_t ( rc_out_of_bounds ) ;
- }
-
- }
-
-#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 cefficient array 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
-{
- const rc_t _ceiling ;
- const rc_t _ceilx2 ;
-
-public:
-
- odd_mapping_r_tag ( int M )
- : _ceiling ( M ) ,
- _ceilx2 ( M + M )
- { } ;
-
- void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
- {
- v += rc_t ( 0.5 ) ;
- 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 ;
-
- odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
- }
-
-#ifdef USE_VC
-
- typedef typename vector_traits < ic_t , vsize > :: ele_v ic_v ;
- typedef typename vector_traits < rc_t , vsize > :: ele_v rc_v ;
-
- void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
- {
- 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) ;
-
- 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 ) ;
- fv ( mask ) = rc_t ( rc_out_of_bounds ) ;
- }
- }
-
-#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 tagged.
-
-template < typename ic_t , typename rc_t , int vsize = 1 >
-class even_mapping_tag
-{
- const rc_t _ceiling ;
-
-public:
-
- even_mapping_tag ( int M )
- : _ceiling ( M - 1 )
- { } ;
-
- 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 ) ;
- return ;
- }
-
- even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
- }
-
-#ifdef USE_VC
-
- typedef typename vector_traits < ic_t , vsize > :: ele_v ic_v ;
- typedef typename vector_traits < rc_t , vsize > :: ele_v rc_v ;
-
- 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 ;
-
- even_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
-
- if ( any_of ( mask ) )
- {
- // v has some 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 ) ;
- }
- }
-
-#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
@@ -913,7 +737,7 @@ public:
/// Mapping for REFLECT boundary conditions for odd splines
///
/// This mapping will map coordinates to [ -0.5 - M-0.5 ] and needs a spline with
-/// a wider bracing. If the spline is created via a bspline object with REFLECT BCs,
+/// a wider left frame. If the spline is created via a bspline object with REFLECT BCs,
/// this will be done automatically. The widened range is due to the point of reflection,
/// which is half a unit spacing 'further out' than for MIRROR BCs.
/// In my initial implementation I was using a coordinate shift by 0.5 to put the
@@ -922,29 +746,35 @@ public:
/// consistency and easier automatic testing - this way, the restoration with grid
/// coordinates, which is used in the roundtrip test, will produce the expected result.
/// So, the left point of reflection as at -0.5, the right point of reflection at M - 0.5.
-
-// TODO I think this is wrong, the initial application of odd_mapping_mirror
-// can't have an effect
+/// This mapping is slower than the other ones since it needs more arithmetic.
template < typename ic_t , typename rc_t , int vsize = 1 >
class odd_mapping_reflect
{
- const rc_t _ceiling ;
- const rc_t _ceilx2 ;
+ const rc_t upper ;
+ const rc_t lower ;
public:
odd_mapping_reflect ( int M )
- : _ceiling ( M ) ,
- _ceilx2 ( M + M )
+ : upper ( M - rc_t(0.5) ) ,
+ lower ( rc_t(-0.5) )
{ } ;
void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
{
- 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 ) ;
+ auto cc = v ;
+ if ( cc < lower )
+ cc = 2 * lower - cc ;
+ if ( cc > upper )
+ {
+ cc -= lower ;
+ cc = std::fmod ( cc , 2 * ( upper - lower ) ) ;
+ cc += lower ;
+ if ( cc > upper )
+ cc = 2 * upper - cc ;
+ }
+ odd_mapping_raw<ic_t,rc_t,vsize>() ( cc , iv , fv ) ;
}
#ifdef USE_VC
@@ -956,41 +786,21 @@ public:
void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
{
- 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
-} ;
-
-template < typename ic_t , typename rc_t , int vsize = 1 >
-class odd_mapping_reflect_raw
-{
-
-public:
-
- odd_mapping_reflect_raw ()
- { } ;
-
- void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
- {
- v += rc_t ( 1.0 ) ; // shift coordinate further to use with wider brace
- odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
- }
-
-#ifdef USE_VC
-
- typedef typename vector_traits < ic_t , vsize > :: ele_v ic_v ;
- typedef typename vector_traits < rc_t , vsize > :: ele_v rc_v ;
+ auto cc = v - lower ;
+ auto w = upper - lower ;
- /// this is the operator() for vectorized operation
+ cc = abs ( cc ) ; // left mirror, v is now >= 0
- void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
- {
- v += rc_t ( 1.0 ) ; // shift coordinate further to use with wider brace
- odd_mapping_raw<ic_t,rc_t,vsize>() ( v , iv , fv ) ;
+ if ( any_of ( cc > w ) )
+ {
+ cc = v_fmod ( cc , 2 * w ) ; // map to one full period
+ cc -= w ; // ccenter
+ cc = abs ( cc ) ; // map to half period
+ cc = w - cc ; // flip
+ }
+
+ cc += lower ;
+ odd_mapping_raw<ic_t,rc_t,vsize>() ( cc , iv , fv ) ;
}
#endif
@@ -999,30 +809,35 @@ public:
/// 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.
+/// even case, we don't need the widened frame, since it has smaller support.
template < typename ic_t , typename rc_t , int vsize = 1 >
class even_mapping_reflect
{
- const rc_t _ceiling ;
- const rc_t _ceilx2 ;
+ const rc_t upper ;
+ const rc_t lower ;
public:
even_mapping_reflect ( int M )
- : _ceiling ( M ) ,
- _ceilx2 ( M + M )
- { } ;
+ : upper ( M - rc_t(0.5) ) ,
+ lower ( rc_t(-0.5) )
+ { } ;
void operator() ( rc_t v , ic_t& iv , rc_t& fv ) const
{
- 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 ) ;
+ auto cc = v ;
+ if ( cc < lower )
+ cc = 2 * lower - cc ;
+ if ( cc > upper )
+ {
+ cc -= lower ;
+ cc = std::fmod ( cc , 2 * ( upper - lower ) ) ;
+ cc += lower ;
+ if ( cc > upper )
+ cc = 2 * upper - cc ;
+ }
+ even_mapping_raw<ic_t,rc_t,vsize>() ( cc , iv , fv ) ;
}
#ifdef USE_VC
@@ -1034,10 +849,21 @@ public:
void operator() ( rc_v v , ic_v & iv , rc_v & fv ) const
{
- 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 ) ;
+ auto cc = v - lower ;
+ auto w = upper - lower ;
+
+ cc = abs ( cc ) ; // left mirror, v is now >= 0
+
+ if ( any_of ( cc > w ) )
+ {
+ cc = v_fmod ( cc , 2 * w ) ; // map to one full period
+ cc -= w ; // ccenter
+ cc = abs ( cc ) ; // map to half period
+ cc = w - cc ; // flip
+ }
+
+ cc += lower ;
+ even_mapping_raw<ic_t,rc_t,vsize>() ( cc , iv , fv ) ;
}
#endif
@@ -1050,190 +876,101 @@ template < typename ic_t ,
typename rc_t ,
int vsize ,
typename map_func >
-map_func pick_mapping ( bc_code bc , ///< mapping mode
+map_func pick_mapping ( bc_code bc , ///< spline's boundary condition
+ map_code mc , ///< 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:
+ // periodic splines are okay to evaluate in [0...M]
+ if ( bc == PERIODIC )
+ M++ ;
+ // spline_degree is only needed to tell odd from even splines:
if ( spline_degree & 1 )
{
- switch ( bc )
+ switch ( mc )
{
- case RAW :
+ case MAP_RAW :
{
// raw mapping simply splits the coordinate with no checks
// whatsoever. So it's constructor doesn't take an argument.
return odd_mapping_raw < ic_t , rc_t , vsize > () ;
break ;
}
- case LIMIT :
+ case MAP_LIMIT :
{
// limit mapping returns the nearest boundary value for
- // out-of-bound coordinates.
+ // out-of-bound coordinates. aka 'clamp'
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 :
+ case MAP_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 :
- {
- // 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 :
- {
- /// 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 clamp value at the bounds.
- return odd_mapping_limit < ic_t , rc_t , vsize > ( M ) ;
- break ;
- }
- case MIRROR :
+ case MAP_MIRROR :
{
return odd_mapping_mirror < ic_t , rc_t , vsize > ( M ) ;
break ;
}
- case REFLECT :
- case SPHERICAL :
+ case MAP_REFLECT :
{
return odd_mapping_reflect < ic_t , rc_t , vsize > ( M ) ;
break ;
}
- case PERIODIC :
+ case MAP_PERIODIC :
{
- return odd_mapping_periodic < ic_t , rc_t , vsize > ( M + 1 ) ;
+ return odd_mapping_periodic < ic_t , rc_t , vsize > ( M ) ;
break ;
}
default:
{
- std::cerr << "unsupported BC code " << bc << std::endl ;
- throw not_supported ( "mapping for this BC code is not supported" ) ;
+ std::cerr << "unsupported mapping code " << mc << std::endl ;
+ throw not_supported ( "automatic mapping for this mapping code is not supported" ) ;
break ;
}
}
}
else
{
- switch ( bc )
+ switch ( mc )
{
- case RAW :
+ case MAP_RAW :
{
return even_mapping_raw < ic_t , rc_t , vsize > () ;
break ;
}
- case LIMIT :
+ case MAP_LIMIT :
{
return even_mapping_limit < ic_t , rc_t , vsize > ( M ) ;
break ;
}
- case LIMITP :
- {
- return even_mapping_limit < ic_t , rc_t , vsize > ( M + 1 ) ;
- break ;
- }
- case REJECT :
+ case MAP_REJECT :
{
return even_mapping_reject < ic_t , rc_t , vsize > ( M ) ;
break ;
}
- case REJECTP :
- {
- return even_mapping_reject < ic_t , rc_t , vsize > ( M + 1 ) ;
- break ;
- }
- case TAG :
- {
- return even_mapping_tag < ic_t , rc_t , vsize > ( M ) ;
- break ;
- }
- case TAGP :
- {
- 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 :
- {
- return even_mapping_limit < ic_t , rc_t , vsize > ( M ) ;
- break ;
- }
- case MIRROR :
+ case MAP_MIRROR :
{
return even_mapping_mirror < ic_t , rc_t , vsize > ( M ) ;
break ;
}
- case REFLECT :
- case SPHERICAL :
+ case MAP_REFLECT :
{
return even_mapping_reflect < ic_t , rc_t , vsize > ( M ) ;
break ;
}
- case PERIODIC :
+ case MAP_PERIODIC :
{
- return even_mapping_periodic < ic_t , rc_t , vsize > ( M + 1 ) ;
+ return even_mapping_periodic < ic_t , rc_t , vsize > ( M ) ;
break ;
}
default:
{
- std::cerr << "unsupported BC code " << bc << std::endl ;
- throw not_supported ( "mapping for this BC code is not supported" ) ;
+ std::cerr << "unsupported mapping code " << mc << std::endl ;
+ throw not_supported ( "automatic mapping for this mapping code is not supported" ) ;
break ;
}
}
@@ -1257,11 +994,11 @@ struct mapping
map_func_v map_v ;
#endif
- mapping ( bc_code bc , int spline_degree , int M )
+ mapping ( bc_code bc , map_code mc , int spline_degree , int M )
{
- map = pick_mapping < ic_t , rc_t , vsize , map_func > ( bc , spline_degree , M ) ;
+ map = pick_mapping < ic_t , rc_t , vsize , map_func > ( bc , mc , spline_degree , M ) ;
#ifdef USE_VC
- map_v = pick_mapping < ic_t , rc_t , vsize , map_func_v > ( bc , spline_degree , M ) ;
+ map_v = pick_mapping < ic_t , rc_t , vsize , map_func_v > ( bc , mc , spline_degree , M ) ;
#endif
}
@@ -1274,9 +1011,6 @@ struct mapping
/// to which the mapping is applied is determined by an additional parameter,
/// axis.
/// The additional routines (without axis parameter) iterate over all axes.
-/// Since in this object we keep base class pointers to the actual mappings,
-/// we need additional code to assure proper deletion of the mapping objects,
-/// which are created by new. Hence the copy constructor and assignment operator.
template < typename ic_t , typename rc_t , int dimension , int vsize = 1 >
struct nd_mapping
@@ -1285,6 +1019,7 @@ struct nd_mapping
typedef typename vigra::TinyVector < mapping_type , dimension > nd_mapping_type ;
typedef typename vigra::TinyVector < bc_code , dimension > bcv_type ;
+ typedef typename vigra::TinyVector < map_code , dimension > mcv_type ;
typedef typename vigra::MultiArrayShape < dimension > :: type shape_type ;
typedef vigra::TinyVector < rc_t , dimension > nd_rc_t ;
@@ -1302,6 +1037,7 @@ struct nd_mapping
nd_mapping_type map ; // container for the mappings
bcv_type bcv ;
+ mcv_type mcv ;
int spline_degree ;
nd_ic_t shape ;
@@ -1310,16 +1046,18 @@ struct nd_mapping
/// other constructors delegate to it.
nd_mapping ( const bcv_type& _bcv ,
+ const mcv_type& _mcv ,
const int& _spline_degree ,
const shape_type& _shape )
: bcv ( _bcv ) ,
+ mcv ( _mcv ) ,
spline_degree ( _spline_degree ) ,
shape ( _shape )
{
// here vigra's shape_type, which uses long ints, is explicitly cast down to int.
// TODO is this a problem?
for ( int d = 0 ; d < dimension ; d++ )
- map [ d ] = mapping_type ( bcv[d] , spline_degree , ic_t ( shape[d] ) ) ;
+ map [ d ] = mapping_type ( bcv[d] , mcv[d] , spline_degree , ic_t ( shape[d] ) ) ;
}
/// apply the mapping along axis 'axis' to coordinate v, resulting in the setting
diff --git a/remap.h b/remap.h
index 4346966..2c41975 100644
--- a/remap.h
+++ b/remap.h
@@ -594,6 +594,7 @@ int remap ( const MultiArrayView < coordinate_traits < coordinate_type > :: dime
// create an evaluator over the bspline
typedef evaluator < coordinate_type , value_type > evaluator_type ;
+
evaluator_type ev ( bsp ) ;
// and call the other remap variant,
--
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