vspline 1.1.0
Generic C++11 Code for Uniform B-Splines
Classes | Namespaces | Typedefs | Functions
transform.h File Reference

set of generic remap, transform and apply functions More...

#include "multithread.h"
#include "eval.h"
#include "poles.h"
#include "convolve.h"
#include "wielding.h"

Go to the source code of this file.

Classes

struct  vspline::detail::grev_generator< ET >
 we need a 'generator functor' to implement grid_eval using the code in wielding.h. this functor precalculates the b-spline evaluation weights corresponding to the coordinates in the grid and stores them in vectorized format, to speed up their use as much as possible. More...
 
struct  vspline::detail::grid_eval_functor< _inner_type, ic_type >
 grid_eval_functor is used for 'generalized' grid evaluation, where the functor passed in is not a bspline evaluator. For this general case, we can't precalculate anything to speed things up, all we need to do is pick the grid values and put them together to form arguments for calls to the functor. More...
 

Namespaces

namespace  vspline
 
namespace  vspline::detail
 

Typedefs

template<unsigned int dimension>
using vspline::bcv_type = vigra::TinyVector< vspline::bc_code, dimension >
 a type for a set of boundary condition codes, one per axis More...
 
template<unsigned int dimension, typename rc_ele_type = float>
using vspline::grid_spec = vigra::TinyVector< vigra::MultiArrayView< 1, rc_ele_type >, dimension >
 

Functions

template<typename unary_functor_type , unsigned int dimension>
void vspline::transform (const unary_functor_type &functor, const vigra::MultiArrayView< dimension, typename unary_functor_type::in_type > &input, vigra::MultiArrayView< dimension, typename unary_functor_type::out_type > &output, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 implementation of two-array transform using wielding::coupled_wield. More...
 
template<typename sink_functor_type , unsigned int dimension>
void vspline::reduce (const sink_functor_type &functor, const vigra::MultiArrayView< dimension, typename sink_functor_type::in_type > &input, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 reduce function processing by value. sink_functor_type is a functor taking values (which are loaded from an array). The reduce function, by delegation to wielding::value_reduce, performs the feeding part by parcelling the data, the functor receives either single in_type or vectorized data. The functor has to be a construct capable of accumulating data with several threads in parallel. A typical example would be a functor cumulating to some internal structure and 'offloading' to a mutex-protected pooling target during it's destructor. The wielding code creates a copy of the functor which is passed here for each worker thread and guarantees that this copy's destructor is executed before the call to reduce returns, so the user can be assured of having collected the partial results of all worker threads in the pooling target. The functor needs a copy constructor to communicate the pooling target to the copies made for the individual threads, see the index-based reduce function below this one for an illustration of the simple and elegant mechanism. Note that sink_functor will normally be a type inheriting from vspline::sink_functor, a close relative of vspline::unary_functor and also defined in unary_functor.h. This inheritance does not produce any functionality, just the standard vspline functor argument type system, which this code relies on. Of course the sink functor can be coded 'manually' as well - inheriting form vspline::sink_functor is merely convenient. More...
 
template<class unary_functor_type >
void vspline::transform (const unary_functor_type &functor, vigra::MultiArrayView< unary_functor_type::dim_in, typename unary_functor_type::out_type > &output, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 reduce function processing by coordinate. sink_functor_type is a functor receiving integral coordinates, with no fixed meaning, but usually pertaining to an array. The reduce function feeds all coordinates encompassed by the shape, either in vectorized or single form. Here's an example of a functor for this function: More...
 
template<typename sink_functor_type >
void vspline::reduce (const sink_functor_type &functor, typename std::conditional<(sink_functor_type::dim_in==1), int, vigra::TinyVector< int, sink_functor_type::dim_in > > ::type shape, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 while the function above fills an array by calling a functor with the target coordinates, this function accumulates the results by means of a sink_functor. The effect is as if an index-based transform was executed first, filling an intermediate array, followed by a call to reduce taking the intermediate array as input. The combination of both operations 'cuts out' the intermediate. Since this function does not have a target array, the shape has to be communicated explicitly. The functor will be invoked for every possible coordinate within the shape's scope. This is done by the function 'index_reduce' in wielding.h the type of the argument 'shape' is derived from the sink_functor's 'dim_in' value, to make the code 'respect the singular': if the data are 1D, reduce will accept a single integer as 'shape', whereas nD data require a TinyVector of int. For suitable sink_functors, see the remarks and example given above in the comments for the array-base overload of 'reduce'. More...
 
template<typename unary_functor_type , unsigned int dimension>
void vspline::apply (const unary_functor_type &ev, vigra::MultiArrayView< dimension, typename unary_functor_type::out_type > &output, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 we code 'apply' as a special variant of 'transform' where the output is also used as input, so the effect is to feed the unary functor each 'output' value in turn, let it process it and store the result back to the same location. While this looks like a rather roundabout way of performing an apply, it has the advantage of using the same type of functor (namely one with const input and writable output), rather than a different functor type which modifies it's argument in-place. While, at this level, using such a functor looks like a feasible idea, It would require specialized code 'further down the line' when complex functors are built with vspline's functional programming tools: the 'apply-capable' functors would need to read the output values first and write them back after anyway, resulting in the same sequence of loads and stores as we get with the current 'fake apply' implementation. From the direct delegation to the two-array overload of 'transform', you can see that this is merely syntactic sugar. More...
 
template<typename original_type , typename result_type , typename coordinate_type , unsigned int cf_dimension, unsigned int trg_dimension, int bcv_dimension = cf_dimension>
void vspline::remap (const vigra::MultiArrayView< cf_dimension, original_type > &input, const vigra::MultiArrayView< trg_dimension, coordinate_type > &coordinates, vigra::MultiArrayView< trg_dimension, result_type > &output, bcv_type< bcv_dimension > bcv=bcv_type< bcv_dimension >(MIRROR), int degree=3, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 Implementation of 'classic' remap, which directly takes an array of values and remaps it, internally creating a b-spline of given order just for the purpose. This is used for one-shot remaps where the spline isn't reused, and specific to b-splines, since the functor used is a b-spline evaluator. The spline defaults to a cubic b-spline with mirroring on the bounds. More...
 
template<typename functor_type >
void vspline::detail::transform (std::true_type, grid_spec< functor_type::dim_in, typename functor_type::in_ele_type > &grid, const functor_type &functor, vigra::MultiArrayView< functor_type::dim_in, typename functor_type::out_type > &result, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 given the generator functor 'grev_generator' (above), performing the grid_eval becomes trivial: construct the generator and pass it to the wielding code. The signature looks complex, but the types needed can all be inferred by ATD - callers just need to pass a grid_spec object, a b-spline evaluator and a target array to accept the result - plus, optionally, a pointer to a cancellation atomic, which, if given, may be set by the calling code to gracefully abort the calculation at the next convenient occasion. Note that grid_eval is specific to b-spline evaluation and will not work with any evaluator type which is not a 'plain', 'raw' b-spline evaluator. So, even functors gained by calling vspline's factory functions (make_evaluator, make_safe_evaluator) will not work here, since they aren't of type vspline::evaluator (they are in fact grok_type). To use arbitrary functors on a grid, use gen_grid_eval, below. More...
 
template<typename functor_type >
void vspline::detail::transform (std::false_type, grid_spec< functor_type::dim_in, typename functor_type::in_ele_type > &grid, const functor_type &functor, vigra::MultiArrayView< functor_type::dim_in, typename functor_type::out_type > &result, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 generalized grid evaluation. The production of result values from input values is done by an instance of grid_eval_functor, see above. The template argument, ev_type, has to be a functor (usually this will be a vspline::unary_functor). If the functor's in_type has dim_in components, grid_spec must also contain dim_in 1D arrays, since ev's input is put together by picking a value from each of the arrays grid_spec points to. The result obviously has to have as many dimensions. More...
 
template<typename functor_type >
void vspline::transform (grid_spec< functor_type::dim_in, typename functor_type::in_ele_type > &grid, const functor_type &functor, vigra::MultiArrayView< functor_type::dim_in, typename functor_type::out_type > &result, int njobs, vspline::atomic< bool > *p_cancel, std::false_type)
 This overload of 'transform' takes a vspline::grid_spec as it's first parameter. This object defines a grid where each grid point is made up from components taken from the per-axis grid values. Let's say the grid_spec holds. More...
 
template<typename functor_type >
void vspline::transform (grid_spec< 1, typename functor_type::in_ele_type > &grid, const functor_type &functor, vigra::MultiArrayView< 1, typename functor_type::out_type > &result, int njobs, vspline::atomic< bool > *p_cancel, std::true_type)
 overload of grid-based transform for 1D grids. Obviously, this is a case for using 'transform' directly taking the single set of per-axis grid values as input, and this overload is here only for completeness' sake - 'normal' user code will use an array-based transform straight away. More...
 
template<typename functor_type >
void vspline::transform (grid_spec< functor_type::dim_in, typename functor_type::in_ele_type > &grid, const functor_type &functor, vigra::MultiArrayView< functor_type::dim_in, typename functor_type::out_type > &result, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 
template<typename functor_type >
void vspline::grid_eval (grid_spec< functor_type::dim_in, typename functor_type::in_ele_type > &grid, const functor_type &functor, vigra::MultiArrayView< functor_type::dim_in, typename functor_type::out_type > &result, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 for backward compatibility, we provide 'grid_eval'. This is deprecated and will eventually go away, new code should use 'transform' above. More...
 
template<typename functor_type >
void vspline::gen_grid_eval (grid_spec< functor_type::dim_in, typename functor_type::in_ele_type > &grid, const functor_type &functor, vigra::MultiArrayView< functor_type::dim_in, typename functor_type::out_type > &result, int njobs=vspline::default_njobs, vspline::atomic< bool > *p_cancel=0)
 for backward compatibility, we provide 'gen_grid_eval'. This is deprecated and will eventually go away, new code should use 'transform' above. More...
 
template<unsigned int dimension, typename value_type , typename math_ele_type = typename vigra::NumericTraits < typename vigra::ExpandElementResult < value_type > :: type > :: RealPromote, size_t vsize = vspline::vector_traits < math_ele_type > :: size>
void vspline::restore (const vspline::bspline< value_type, dimension > &bspl, vigra::MultiArrayView< dimension, value_type > &target)
 restore restores the original data from the b-spline coefficients. This is done efficiently using a separable convolution, the kernel is simply a unit-spaced sampling of the basis function. Since the filter uses internal buffering, using this routine in-place is safe - meaning that 'target' may be bspl.core itself. math_type, the data type for performing the actual maths on the buffered data, and the type the data are converted to when they are placed into the buffer, can be chosen, but I could not detect any real benefits from using anything but the default, which is to leave the data in their 'native' type. More...
 
template<int dimension, typename value_type , typename math_ele_type = typename vigra::NumericTraits < typename vigra::ExpandElementResult < value_type > :: type > :: RealPromote, size_t vsize = vspline::vector_traits<math_ele_type>::size>
void vspline::restore (vspline::bspline< value_type, dimension > &bspl)
 overload of 'restore' writing the result of the operation back to the array which is passed in. This looks like an in-place operation, but the data are in fact moved to a buffer stripe by stripe, then the arithmetic is done on the buffer and finally the buffer is written back. This is repeated for each dimension of the array. More...
 

Detailed Description

set of generic remap, transform and apply functions

My foremost reason to have efficient B-spline processing is the formulation of generic remap-like functions. remap() is a function which takes an array of real-valued nD coordinates and an interpolator over a source array. Now each of the real-valued coordinates is fed into the interpolator in turn, yielding a value, which is placed in the output array at the same place the coordinate occupies in the coordinate array. To put it concisely, if we have

remap defines the operation

t[j] = i(c[j]) for all j

Now we widen the concept of remapping to a 'transform' function. Instead of limiting the process to the use of an 'interpolator', we use an arbitrary unary functor transforming incoming values to outgoing values, where the type of the incoming and outgoing values is determined by the functor. If the functor actually is an interpolator, we have a 'true' remap transforming coordinates into values, but this is merely a special case. So here we have:

transform performs the operation

t[j] = f(c[j]) for all j

remaps/transforms to other-dimensional objects are supported. This makes it possible to, for example, remap from a volume to a 2D image, using a 2D coordinate array containing 3D coordinates ('slicing' a volume)

There is also a variant of this transform function in this file, which doesn't take an input array. Instead, for every target location, the location's discrete coordinates are passed to the unary_functor type object. This way, transformation-based remaps can be implemented easily: the user code just has to provide a suitable functor to yield values for coordinates. This functor will internally take the discrete incoming coordinates (into the target array) and take it from there, eventually producing values of the target array's value_type.

Here we have:

'index-based' transform performs the operation

t[j] = f(j) for all j

This file also has code to evaluate a b-spline at coordinates in a grid, which can be used for scaling, and for separable geometric transformations.

Finally there is a function to restore the original data from a b-spline to the precision possible with the given data type and degree of the spline. This is done with a separable convolution, using a unit-stepped sampling of the basis function as the convolution kernel along every axis.

Let me reiterate the strategy used to perform the transforms and remaps in this file. The approach is functional: A 'processing chain' is set up and encoded as a functor providing two evaluation functions: one for 'single' data and one for vectorized data. This functor is applied to the data by 'wielding' code, which partitions the data into several jobs to be performed by individual worker threads, invokes the worker threads, and, once in the worker thread, feeds the data to the functor in turn, using hardware vectorization if possible. So while at the user code level a single call to some 'transform' or 'remap' routine is issued, passing in arrays of data and functors, all the 'wielding' is done automatically without any need for the user to even be aware of it, while still using highly efficient vector code with a tread pool, potentially speeding up the operation greatly as compared to single-threaded, unvectorized operation, of course depending on the presence of several cores and vector units. On my system (Haswell 4-core i5 with AVX2) the speedup is about one order of magnitude. The only drawback is the production of a hardware-specific binary if vectorization is used. Both Vc use and multithreading can be easily activated/deactivated by #define switches, providing a clean fallback solution to produce code suitable for any target, even simple single-core machines with no vector units. Vectorization will be used if possible - either explicit vectorization (by defining USE_VC) - or autovectorization per default. Defining VSPLINE_SINGLETHREAD will disable multithreading. The code accessing multithreading and/or Vc use is #ifdeffed, so if these features are disabled, their code 'disappears' and the relevant headers are not included, nor do the corresponding libraries have to be present.

Definition in file transform.h.