[vspline] 34/72: removed tf_remap I have done away with tf_remap in favour of only providing index_remap. index_remap can be used to the same effect, but burdens the implementation of a transformation object on the calling code. Yet wrapping an evaluator in an outer object which can do the necessary transformation isn't hard, and doing so is more efficient than using the slightly awkward vspline::transform object (though it's auto-broadcast method seemed like a good idea at the time...) removing tf_remap and vspline::transform made the code in remap.h a good bit slimmer, so I'm quite happy with it now. In pv, I have also removed all use of the mmap object in class evaluator, but leaving the code in mapping.h for the time being does no harm, though I may eventually remove automatic mapping for all BC codes other than those which can also be used as boundary conditions.
Kay F. Jahnke
kfj-guest at moszumanska.debian.org
Sun Jul 2 09:02:40 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 438354cb6910f60f25d0f4aff89a0658eff69e29
Author: Kay F. Jahnke <kfjahnke at gmail.com>
Date: Mon Feb 13 11:37:17 2017 +0100
removed tf_remap
I have done away with tf_remap in favour of only providing index_remap.
index_remap can be used to the same effect, but burdens the implementation
of a transformation object on the calling code. Yet wrapping an evaluator
in an outer object which can do the necessary transformation isn't hard,
and doing so is more efficient than using the slightly awkward
vspline::transform object (though it's auto-broadcast method seemed like
a good idea at the time...)
removing tf_remap and vspline::transform made the code in remap.h a good
bit slimmer, so I'm quite happy with it now.
In pv, I have also removed all use of the mmap object in class evaluator,
but leaving the code in mapping.h for the time being does no harm, though
I may eventually remove automatic mapping for all BC codes other than
those which can also be used as boundary conditions.
---
eval.h | 56 +++++---
interpolator.h | 22 +++
mapping.h | 4 +-
remap.h | 439 +++++----------------------------------------------------
4 files changed, 100 insertions(+), 421 deletions(-)
diff --git a/eval.h b/eval.h
index 770a84e..63a27d8 100644
--- a/eval.h
+++ b/eval.h
@@ -448,6 +448,8 @@ struct alt_bspline_derivative_weight_functor
/// as possible. The strategy is to have all dependencies of the evaluation except for the actual
/// coordinates taken care of by the constructor - and immutable for the evaluator's lifetime.
/// The resulting object has no state which is modified after construction, making it thread-safe.
+/// It also constitutes a 'pure' function in a functional-programming sense, because it has
+/// no mutable state and no side-effects.
/// The eval() variants take only coordinates (in their various guises) as arguments, relying
/// on the state of the evaluator fixed at it's construction. The eval() overloads also
/// form a hierarchy, as evaluation progresses from accepting unsplit real coordinates to
@@ -459,12 +461,29 @@ struct alt_bspline_derivative_weight_functor
/// and finally subjects the resultant values to some postprocessing. All these processing
/// steps can be bound into a single functor, and the calling code can be reduced to polling
/// this functor until it has obtained the desired number of output values.
+/// While the 'unspecialized' evaluator will try and do 'the right thing' by using general
+/// purpose code fit for all eventualities, for time-critical operation there are two specializations
+/// which can be used to make the code faster:
+///
+/// - template argument 'specialize' can be set to 0 to forcibly use (more efficient) nearest
+/// neighbour interpolation, which has the same effect as simply running with degree 0 but avoids
+/// code which isn't needed for nearest neighbour interpolation (like the application of weights,
+/// which is futile under the circumstances, the weight always being 1.0).
+/// specialize can also be set to 1 for explicit bilinear interpolation. Any value but 0 or 1 will
+/// result in the general-purpose code being used.
+///
+/// - template argument 'evenness' can be passed in as 1 to enforce raw-mode coordinate splitting
+/// for an even b-spline. Since this is done 'hardwired' it's potentially slightly faster than
+/// using the general-purpose mapping via an mmap object. Passing 0 for 'evenness' does the equivalent
+/// for odd splines, enforcing a raw mapping. Any other value for 'evenness' will result in the use of
+/// the mmap object passed in at the evaluator's creation, which provides cordinate handling services
+/// with the routines defined in mapping.h.
template < int _dimension , ///< dimension of the spline
class _value_type , ///< type of spline coefficients/result (pixels etc.)
class _rc_type , ///< singular real coordinate, float or double
class _ic_type , ///< singular integral coordinate, currently only int
- bool use_tag = false , ///< flag switching tag checking on/off
+// bool use_tag = false , ///< flag switching tag checking on/off
int specialize = -1 ,
int evenness = -1 >
class evaluator
@@ -626,6 +645,7 @@ public:
/// (in terms of member data) is fixed at construction time and remains constant afterwards.
// TODO: It would be desirable to make this explicit for all member data, but some of them need
// more involved initialization than the one-liners after the colon... so I'm a bit sloppy here.
+ // in the calling code, it's fine to use a const evaluator, that sorts it as well.
evaluator ( const array_type& _coefficients ,
nd_mapping_type _mmap ,
@@ -1467,22 +1487,23 @@ public:
v_eval ( select , weight , result ) ;
}
- if ( use_tag )
- {
- // this bit of code is optional, it's only needed if mapping mode TAG is used.
- // If use_tag is false and mapping mode TAG is used, out-of-range
- // coordinates produce undefined results but the program will not crash.
-
- for ( int d = 0 ; d < dimension ; d++ )
- {
- auto mask = ( tune[d] == rc_type ( rc_out_of_bounds ) ) ; // test for oob value
- if ( any_of ( mask ) )
- {
- for ( int ch = 0 ; ch < channels ; ch++ ) // coordinate is oob, set to 0
- result[ch] ( mask ) = ele_type ( rv_out_of_bounds ) ;
- }
- }
- }
+// use_tag is no longer used.
+// if ( use_tag )
+// {
+// // this bit of code is optional, it's only needed if mapping mode TAG is used.
+// // If use_tag is false and mapping mode TAG is used, out-of-range
+// // coordinates produce undefined results but the program will not crash.
+//
+// for ( int d = 0 ; d < dimension ; d++ )
+// {
+// auto mask = ( tune[d] == rc_type ( rc_out_of_bounds ) ) ; // test for oob value
+// if ( any_of ( mask ) )
+// {
+// for ( int ch = 0 ; ch < channels ; ch++ ) // coordinate is oob, set to 0
+// result[ch] ( mask ) = ele_type ( rv_out_of_bounds ) ;
+// }
+// }
+// }
}
/// here we transform incoming nD coordinates into offsets into the coefficient
@@ -1669,6 +1690,7 @@ public:
else
{
// no specialization; fall back to calling mmap
+ std::cerr << "warning: using fallback to old mapping code" << std::endl ;
mmap ( input , select , tune ) ;
}
diff --git a/interpolator.h b/interpolator.h
index 513043b..34ec708 100644
--- a/interpolator.h
+++ b/interpolator.h
@@ -164,6 +164,28 @@ struct interpolator
}
} ;
+ /// v_eval taking a reference to a simdized coordinate and a pointer
+ /// to memory accommodating the result, expecting this memory to be
+ /// contiguous.
+
+ void v_eval ( const nd_rc_v & cv , // simdized coordinate
+ value_type * result ) const // -> vsize result values
+ {
+ mc_ele_v ev ;
+ v_eval ( cv , ev ) ;
+ if ( channels == 1 )
+ {
+ ev[0].store ( (ele_type*) result ) ;
+ }
+ else
+ {
+ for ( int ch = 0 ; ch < channels ; ch++ )
+ {
+ ev[ch].scatter ( (ele_type*) result , ic_v::IndexesFromZero() * channels + ch ) ;
+ }
+ }
+ } ;
+
#else
enum { vsize = 1 } ;
diff --git a/mapping.h b/mapping.h
index 723c423..e77509f 100644
--- a/mapping.h
+++ b/mapping.h
@@ -541,8 +541,8 @@ public:
/// 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
+/// 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.
diff --git a/remap.h b/remap.h
index 90aca7d..65ed1de 100644
--- a/remap.h
+++ b/remap.h
@@ -73,19 +73,14 @@
/// data access modes here.
///
/// There is also a second set of remap functions in this file, which don't take a
-/// 'warp' array but a functor instead. This functor receives, for every location in the
-/// output, the corresponding (discrete) coordinates, and returns (real) coordinates pointing
-/// into the input array. This variant exists, because a common situation in image transformations
-/// is that the geometric transformation from one image to another is defined by precisely
-/// such a functor, and the creation of the warp array only makes sense if the transformation
-/// is performed several times. With 'cheap' transformations the difference isn't large,
-/// but if the transformation is very complex, reusing it can save a lot of processing.
-/// What I have in mind here is image transformations as used in panorama stitching,
-/// lens correction, etc.
-///
-/// In the 'example' folder, there is a program called pano_extract.cc, which demonstrates
-/// the use of a transformation-based remap. There's also pv.cc which uses remap extensively,
-/// and, being a panorama viewer, gives instant visual feedback showing the results as images.
+/// 'warp' array. Instead, for every target location, the location's discrete coordinates
+/// are passed to the interpolator_type object. This way, transformation-based remaps
+/// can be implemented easily: the user code just has to provide a suitable 'interpolator'
+/// to yield values for coordinates. This interpolator will internally take the discrete
+/// incoming coordinates (into the target array) and transform them as required, internally
+/// producing coordinates suitable for the 'actual' interpolation using a b-spline or some
+/// other object capable of producing values for real coordinates. The routine offering
+/// this service is called index_remap.
///
/// This file also has code to evaluate a b-spline at positions in a mesh grid, which can be used
/// for scaling, and for separable geometric transformations.
@@ -112,167 +107,6 @@ using namespace vigra ;
template < int dimension >
using bcv_type = vigra::TinyVector < bc_code , dimension > ;
-// Next we have some collateral code to get ready for the transformation-based remap
-
-/// declaration of a generalized coordinate transformation. In vspline, we're using
-/// a 'reverse transformation': For every coordinate ct into the target array (the result
-/// of the 'remap' procedure, the final picture) we calculate a corresponding coordinate cs
-/// into the source data/spline coefficients, where we evaluate the interpolator to yield
-/// the value which has to be written to the target array at ct.
-///
-/// target [ ct ] = interpolator ( cs )
-///
-/// target [ ct ] = interpolator ( transform ( ct ) )
-///
-/// the coordinates ct and cs can be of different dimensionality (dim_target, dim_source)
-/// but share their component type, so that
-///
-/// std::is_same < decltype ( cs[0] ) , decltype ( ct[0] ) > :: value == true
-
-template < class component_type , // 1D coordinate type, single value or SimdArray
- int dim_target , // number of dimensions of target array, incoming coordinate
- int dim_source > // number of dimensions of coefficient array, resultant coordinate
-using transform
-= std::function < void ( const vigra::TinyVector < component_type , dim_target > & ,
- vigra::TinyVector < component_type , dim_source > & ) > ;
-
-/// class transformation is a (multi-) functor, handling coordinate transformations
-/// this class simplifies pulling in transformations by automatically providing a brodcasting
-/// method to apply a single-value transformation to every element of an incoming SimdArray,
-/// if there is no vectorized code at hand. This is less efficient, but often just enough,
-/// especially if the transformation is simple.
-///
-/// struct transformation is constructed either with a single coordinate transformation
-/// functor, or, if USE_VC is defined, optionally with a vectorized coordinate transformation
-/// functor as second constructor argument. The functors have to satisfy the two type aliases
-/// transform and vtransform (see above). For illustration, consider this:
-///
-/// tf_identity defines a transformation function template leaving the coordinate unchanged
-///
-/// template < class coordinate_type >
-/// void tf_identity ( const TinyVector < coordinate_type , 2 > & c_in ,
-/// TinyVector < coordinate_type , 2 > & c_out )
-/// {
-/// c_out = c_in ;
-/// }
-///
-/// now we can create a vspline::transformation object, passing in two specializations
-/// of the transformation function, one for unvectorized and one for vectorized coordinates
-///
-/// vspline::transformation < float , 2 , 2 >
-/// tf ( tf_identity<float> , vtf_identity<Vc::float_v> ) ;
-///
-/// we could also using the single-argument constructor to the same effect, not using
-/// vectorized transformations:
-///
-/// vspline::transformation < float , 2 , 2 >
-/// tf ( tf_identity<float> ) ;
-
-template < class rc_type , /// elementary coordinate type
- int dim_target , /// dimension of incoming coordinates (== dim. of target array)
- int dim_source , /// dimension of outgoing coordinates (== dim. of coefficient array)
- int vsize = 1 /// number of elements in simdized coordinate
- >
-struct transformation
-{
- typedef transform < rc_type , dim_target , dim_source > transform_type ;
-
- transform_type tf ; /// functor for single-value coordinate transformation
-
- typedef vigra::TinyVector < rc_type , dim_target > target_nd_rc_type ;
- typedef vigra::TinyVector < rc_type , dim_source > source_nd_rc_type ;
-
-#ifdef USE_VC
-
- typedef typename vector_traits < rc_type , vsize > :: type rc_v ;
-
- typedef vigra::TinyVector < rc_v , dim_target > target_nd_rc_v ;
- typedef vigra::TinyVector < rc_v , dim_source > source_nd_rc_v ;
-
- typedef transform < rc_v , dim_target , dim_source > vtransform_type ;
-
- vtransform_type vtf ; /// functor for vectorized operation
-
- /// broadcast calls the single-element transform repeatedly to emulate vectorized operation.
- /// This is used to build up to a vectorized transformation: usually the unvectorized code
- /// is easier to provide or even already present. By broadcasting it, we get a functioning
- /// prototype and a reference to test emerging vector code against. Once a vectorized version
- /// of the transformation has been written and is passed into class transformation's constructor,
- /// bradcasting is no longer used.
-
- struct broadcast
- {
- transform_type tf ; /// functor for single-value coordinate transformation
-
- broadcast ( transform_type _tf )
- : tf ( _tf )
- { } ;
-
- void operator() ( const target_nd_rc_v & c_target_v ,
- source_nd_rc_v & c_source_v )
- {
- target_nd_rc_type c_target ;
- source_nd_rc_type c_source ;
-
- for ( int e = 0 ; e < vsize ; e++ )
- {
- for ( int d = 0 ; d < dim_source ; d++ )
- c_target[d] = c_target_v[d][e] ;
- tf ( c_target , c_source ) ;
- for ( int d = 0 ; d < dim_target ; d++ )
- c_source_v[d][e] = c_source[d] ;
- }
- }
- } ;
-
-public:
-
- /// the constructor takes tranformation functors for both unvectorized
- /// and vectorized coordinate transformations. If only the unvectorized
- /// transform is passed, it is automatically broadcasted.
-
- transformation ( transform_type _tf ,
- vtransform_type _vtf = 0 )
- : tf ( _tf ) ,
- vtf ( _vtf ? _vtf : broadcast ( _tf ) )
- { } ;
-
-#else
-
-public:
-
- /// if USE_VC is not defined, we only have this single constructor:
-
- transformation ( transform_type _tf )
- : tf ( _tf )
- { } ;
-
-#endif
-
-public :
-
- /// class transformation's operator() delegates to the functor passed in at construction.
- /// it transforms a target coordinate, c_target, into a source coordinate c_source
-
- void operator() ( const target_nd_rc_type& c_target ,
- source_nd_rc_type& c_source ) const
- {
- tf ( c_target , c_source ) ;
- }
-
-#ifdef USE_VC
-
- /// if USE_VC is defined, we overload operator() for vecorized arguments
-
- void operator() ( const target_nd_rc_v& c_target ,
- source_nd_rc_v& c_source ) const
- {
- vtf ( c_target , c_source ) ;
- }
-
-#endif
-} ;
-
/// struct coordinate_iterator_v provides the vectorized equivalent of vigra's
/// MultiCoordinateIterator. The iterator itself is not very elaborate, it only
/// supports dereferencing and preincrement, but for the intended purpose that's
@@ -860,21 +694,19 @@ int remap ( const MultiArrayView < coordinate_traits < coordinate_type > :: dime
return 0 ;
}
-/// struct transform_generator provides target values by performing a sequence of steps:
-///
-/// - start with a discrete coordinate into the target array
-///
-/// - transform it into a real-valued coordinate suitable for use with an interpolator
-/// by means of a vspline::tranformation, which contains the desired geometric transformation
-///
-/// - feed the real-valued coordinate to the interpolator, yielding the result
-///
-/// transform_generator satisfies the interface needed by the fill type of routines,
-/// so it can be fed to fill(), just like a warp_generator. Doing so results in a
-/// 'transform-based remap', see tf_remap below.
+/// index_generator is like warp_generator, but instead of picking coordinates from
+/// an array, the corresponding discrete coordinates are passed to the interpolator.
+/// Instead the interpolator is directly fed discrete target coordinates, and it's up
+/// This is a convenience
+/// routine, because it delivers the logic of both coordinate generation and value
+/// storage, in several threads (due to the use of fill) - and leaves the interesting
+/// work to be done in the 'interpolator'.
+/// In a way calling it 'interpolator' is a bit off since it does much more than merely
+/// interpolating at a given position: it performs the entire processing chain from
+/// discrete target coordinate to final result.
template < int dim_target , typename interpolator_type >
-struct transform_generator
+struct index_generator
{
typedef typename interpolator_type::value_type value_type ;
@@ -882,7 +714,6 @@ struct transform_generator
#ifdef USE_VC
-// typedef typename vigra::ExpandElementResult < value_type > :: type ele_type ;
enum { vsize = interpolator_type :: vsize } ;
#else
@@ -895,22 +726,17 @@ struct transform_generator
typedef typename interpolator_type::rc_type rc_type ;
typedef typename interpolator_type::nd_ic_type nd_ic_type ;
typedef typename interpolator_type::nd_rc_type nd_rc_type ;
-
- typedef transformation < rc_type , dim_target , dim_source , vsize > transform_type ;
const interpolator_type & itp ;
- const transform_type & transform ;
const shape_range_type < dim_target > & range ;
vigra::MultiCoordinateIterator < dim_source > mci ;
- transform_generator
- ( const transform_type & _transform ,
- const interpolator_type & _itp ,
+ index_generator
+ ( const interpolator_type & _itp ,
const shape_range_type < dim_target > & _range
)
- : transform ( _transform ) ,
- itp ( _itp ) ,
+ : itp ( _itp ) ,
range ( _range ) ,
mci ( _range[1] - _range[0] )
#ifdef USE_VC
@@ -921,38 +747,30 @@ struct transform_generator
void operator() ( value_type & target )
{
- typename transform_type::target_nd_rc_type trg_coordinate = * ( mci + range[0] ) ;
- typename transform_type::source_nd_rc_type src_coordinate ;
- transform ( trg_coordinate , src_coordinate ) ;
- itp.eval ( src_coordinate , target ) ;
+ itp.eval ( * ( mci + range[0] ) , target ) ;
++mci ;
}
#ifdef USE_VC
- typedef typename interpolator_type::rc_v rc_v ;
coordinate_iterator_v < rc_type , ic_type , dim_source , vsize > civ ;
void operator() ( value_type * target )
{
- typedef typename transform_type::target_nd_rc_v trg_coordinate_type ;
- trg_coordinate_type trg_coordinate = trg_coordinate_type ( *civ ) ;
- typename transform_type::source_nd_rc_v src_coordinate ;
- transform ( trg_coordinate , src_coordinate ) ;
- itp.v_eval ( src_coordinate , target ) ;
+ itp.v_eval ( *civ , target ) ;
++civ ;
}
#endif
- transform_generator < dim_target , interpolator_type >
+ index_generator < dim_target , interpolator_type >
subrange ( const shape_range_type < dim_target > & range ) const
{
- return transform_generator < dim_target , interpolator_type >
- ( transform , itp , range ) ;
+ return index_generator < dim_target , interpolator_type >
+ ( itp , range ) ;
}
- transform_generator < dim_target , interpolator_type >
+ index_generator < dim_target , interpolator_type >
bindOuter ( const int & c ) const
{
nd_ic_type slice_start = range[0] , slice_end = range[1] ;
@@ -960,211 +778,28 @@ struct transform_generator
slice_start [ dim_target - 1 ] = c ;
slice_end [ dim_target - 1 ] = c ;
- return transform_generator < dim_target , interpolator_type >
- ( transform,
- itp ,
- shape_range_type < dim_target > ( slice_start , slice_end ) ) ;
+ return index_generator < dim_target , interpolator_type >
+ ( itp , shape_range_type < dim_target > ( slice_start , slice_end ) ) ;
}
} ;
-/// implementation of tf_remap() by delegation to the more general fill() routine, passing
-/// in the coordinate transform and the interpolator via a transform_generator object which
-/// handles the full chain of operations from picking the right discrete target coordinates
-/// to producing result values.
-/// Having the generator object makes the implementation of the remap routine trivial:
-/// just create the generator and pass it to fill().
-///
-/// this function provides a generalized geometric transformation
-/// In image transformations, often the source location at which interpolation
-/// of the source image is required is defined by a transformation of the
-/// target coordinate. There are (at least) two approaches to handle this.
-/// The first is what is implemented in the previous implementation of remap: we
-/// calculate a 'warp' array full of transformed coordinates and process it en bloc.
-/// Alternatively, with a coordinate transformation at hand, we can write a remap
-/// variant which performs the coordinate transformation for each target coordinate
-/// as it goes along, saving the construction of the warp array.
-///
-/// The transformation is passed in by using a 'transformation' object, which is a
-/// wrapper around the actual transformation routine(s). Note that this routine is
-/// also suitable for any type of interpolator satisfying the interface defined in
-/// class interpolator in interpolator.h.
-
template < class interpolator_type , // type satisfying the interface in class interpolator
int dim_target > // dimension of target array
-void tf_remap ( const interpolator_type & ev ,
- transformation < typename interpolator_type::rc_type ,
- dim_target ,
- interpolator_type::dimension ,
- interpolator_type::vsize > tf ,
- MultiArrayView < dim_target ,
- typename interpolator_type::value_type > & output ,
- bool use_vc = true
- )
+void index_remap ( const interpolator_type & ev ,
+ MultiArrayView < dim_target ,
+ typename interpolator_type::value_type > & output ,
+ bool use_vc = true
+ )
{
typedef typename interpolator_type::value_type value_type ;
typedef TinyVector < int , dim_target > nd_ic_type ;
- typedef transform_generator < dim_target , interpolator_type > gen_t ;
+ typedef index_generator < dim_target , interpolator_type > gen_t ;
shape_range_type < dim_target > range ( nd_ic_type() , output.shape() ) ;
- gen_t gen ( tf , ev , range ) ;
+ gen_t gen ( ev , range ) ;
fill < gen_t , dim_target > ( gen , output , use_vc ) ;
}
-/// this highest-level transform-based remap takes an input array and creates
-/// a b-spline over it internally. Then it calles the previous routine, which
-/// takes an interpolator as it's first parameter. Like with warp-array-based remap,
-/// this is for on-the-fly remapping where the b-spline won't be reused.
-
-template < typename input_value_type , // type of values in input
- typename output_value_type , // type of values in output
- typename rc_type , // elementary type of coordinates
- int dim_target , // dimension of incoming coordinates (== dim. of target array)
- int dim_source > // dimension of outgoing coordinates (== dim. of coefficient array)
-void tf_remap ( MultiArrayView < dim_source ,
- input_value_type > input ,
- transformation < rc_type ,
- dim_target ,
- dim_source ,
- vector_traits < typename vigra::ExpandElementResult < input_value_type > :: type > :: vsize > tf ,
- MultiArrayView < dim_target , output_value_type > & output ,
- bcv_type < dim_source > bcv = bcv_type < dim_source > ( MIRROR ) ,
- int degree = 3 ,
- bool use_vc = true ,
- int nthreads = ncores
- )
-{
- // create the bspline object
- bspline < input_value_type , dim_source > bsp ( input.shape() , degree , bcv ) ;
- // copy in the data
- bsp.prefilter ( use_vc , nthreads * 8 , input ) ;
- // create an evaluator for this spline
- typedef evaluator < dim_source , input_value_type , rc_type , int > evaluator_type ;
- evaluator_type ev ( bsp ) ;
- tf_remap < evaluator_type , dim_target >
- ( ev , tf , output , use_vc , nthreads ) ;
-}
-
-/// this function creates a warp array instead of evaluating a spline. The warp array can
-/// then be used as input to a warp-array-based remap. This way, if the transformation
-/// is computationally expensive and it's resulting transformed coordinates can be reused
-/// several times, the coordinate transformation can be separated from the remap procedure.
-/// This is the single-threaded routine; further below we will use divide_and_conquer
-/// to apply this routine to several chunks of data in parallel.
-
-// TODO rewrite using fill()
-
-template < class value_type , // like, float, double, complex<>, pixels, TinyVectors etc.
- class rc_type , // float or double, for coordinates
- int dim_in , // number of dimensions of input array (and of 'warp' coordinates)
- int dim_out > // number of dimensions of output array
-void st_make_warp_array ( shape_range_type < dim_out > range ,
- transformation < rc_type ,
- dim_in ,
- dim_out ,
- vector_traits < typename vigra::ExpandElementResult < value_type > :: type > :: vsize > tf ,
- MultiArrayView < dim_out ,
- vigra::TinyVector < rc_type , dim_in > >
- _warp_array ,
- bool use_vc = true
- )
-{
- auto warp_array = _warp_array.subarray ( range[0] , range[1] ) ;
- auto offset = range[0] ;
-
- /// incoming coordinates (into output, which has dim_out dims)
- /// generated by vigra's CoupledIterator or, in the vectorized code,
- /// by coordinate_iterator_v
-
- typedef vigra::TinyVector < rc_type , dim_out > nd_rc_in ;
-
- /// transformed coordinates (into input, which has dim_in dims)
- /// populating the resulting warp array
-
- typedef vigra::TinyVector < rc_type , dim_in > nd_rc_out ;
-
- nd_rc_in cs_in ;
- nd_rc_out cs_out ;
-
- auto target_it = createCoupledIterator ( warp_array ) ;
- int leftover = warp_array.elementCount() ;
-
-#ifdef USE_VC
-
- typedef typename ExpandElementResult<value_type>::type ele_type ;
- typedef typename vector_traits < ele_type > :: type ele_v ;
- const int vsize = ele_v::Size ;
- typedef typename vector_traits < rc_type , vsize > :: type rc_v ;
-
- /// vecorized incoming coordinates (into output, which has dim_out dims)
- typedef vigra::TinyVector < rc_v , dim_out > nd_rc_in_v ;
-
- /// vectorized transformed coordinates (into input, which has dim_in dims)
- typedef vigra::TinyVector < rc_v , dim_in > nd_rc_out_v ;
-
- if ( use_vc )
- {
- int aggregates = warp_array.elementCount() / vsize ; // number of full vectors
- leftover = warp_array.elementCount() - aggregates * vsize ; // any leftover single values
- nd_rc_out_v target_buffer ; // vectorized warp coordinates
-
- // use utility class coordinate_iterator_v to produce vectorized target coordinates:
- coordinate_iterator_v < rc_type , int , dim_out , vsize >
- civ ( offset , offset + warp_array.shape() ) ;
-
- for ( int a = 0 ; a < aggregates ; a++ )
- {
- tf ( *civ , target_buffer ) ; // perform the coordinate transformation
- ++civ ;
- // interleave target_buffer into warp array (TODO: maybe use scatter?)
- for ( int e = 0 ; e < vsize ; e++ )
- {
- for ( int d = 0 ; d < dim_in ; d++ )
- target_it.template get<1>()[d] = target_buffer[d][e] ;
- ++target_it ;
- }
- }
- }
-
-#endif // USE_VC
-
- // process leftovers, if any - if vc isn't used, this loop does all the processing
- while ( leftover-- )
- {
- // process leftovers
- cs_in = target_it.template get<0>() + offset ;
- tf ( cs_in , cs_out ) ;
- target_it.template get<1>() = cs_out ;
- ++target_it ;
- }
-}
-
-/// multithreaded code to make a warp array using a coordinate transformation.
-
-template < class value_type , // like, float, double, complex<>, pixels, TinyVectors etc.
- class rc_type , // float or double for coordinates
- int dim_in , // number of dimensions of input array (and of 'warp' coordinates)
- int dim_out > // number of dimensions of output array
-void make_warp_array ( transformation < rc_type ,
- dim_in ,
- dim_out ,
- vector_traits < typename vigra::ExpandElementResult < value_type > :: type > :: vsize > tf ,
- MultiArrayView < dim_out ,
- vigra::TinyVector < rc_type , dim_in > >
- & warp_array ,
- bool use_vc = true
- )
-{
- int njobs = warp_array.shape ( dim_out - 1 ) / 4 ;
- if ( njobs < 16 )
- njobs = 16 ;
-
- shape_range_type < dim_out > range ( shape_type < dim_out > () , warp_array.shape() ) ;
-
- multithread ( st_make_warp_array < value_type , rc_type , dim_in , dim_out > ,
- njobs , range , tf , warp_array , use_vc ) ;
-
-}
-
namespace detail // workhorse code for grid_eval
{
// in grid_weight, for every dimension we have a set of ORDER weights
--
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