[vspline] 52/72: wrote restore function using grid_eval, changes to roundtrip.cc
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 64234aa43aaa4ab79e37aa7e3061014357eaf480
Author: Kay F. Jahnke <kfjahnke at gmail.com>
Date: Tue May 2 18:45:02 2017 +0200
wrote restore function using grid_eval, changes to roundtrip.cc
---
README.rst | 18 ++++++++-------
basis.h | 13 ++++++-----
example/roundtrip.cc | 63 +++++++++++++++++++++++++++++++++++++++++++++-------
prefilter.h | 58 -----------------------------------------------
remap.h | 57 +++++++++++++++++++++++++++++++++++++++++++++++
vspline.doxy | 2 +-
6 files changed, 131 insertions(+), 80 deletions(-)
diff --git a/README.rst b/README.rst
index feaf541..af28cff 100644
--- a/README.rst
+++ b/README.rst
@@ -17,7 +17,7 @@ Fellow, IEEE.
While there are several freely available packets of B-spline code, I failed to find
one which is comprehensive, efficient and generic at once. vspline attempts to be
-all that, making use of generic programming in C++, and of common, but often underused
+all that, making use of generic programming in C++11, and of common, but often underused
hardware features in modern processors. Overall, there is an emphasis on speed, even
if this makes the code more complex. I tried to eke as much performance out of the
hardware at my disposal as possible, only compromising when the other design goals
@@ -66,19 +66,21 @@ On the evaluation side I provide
- evaluation of the spline at point locations in the defined range
- evaluation of the spline's derivatives
- mapping of arbitrary coordinates into the defined range
-- evaluation of nD arrays of coordinates (generalized remap() function)
-- coordinate-transformation based remap function
+- evaluation of nD arrays of coordinates ('remap' function)
+- coordinate-fed remap function ('index_remap')
+- functor-based remap, aka 'transform' function
+- functor-based 'apply' function
-The code at the very core of my B-spline coefficient generation code evolved from the code by Philippe Thévenaz which he published here_, with some of the boundary condition treatment code derived from formulae which Philippe Thévenaz communicated to me. Next I needed code to handle multidimensional arrays in a generic fashion in C++. I chose to use Vigra_. Since my work has a focus on signal (and, more specifically image) processing, it's an excellent tool, as it provides a large body o [...]
+The code at the very core of my B-spline coefficient generation code evolved from the code by Philippe Thévenaz which he published here_, with some of the boundary condition treatment code derived from formulae which Philippe Thévenaz communicated to me. Next I needed code to handle multidimensional arrays in a generic fashion in C++. I chose to use Vigra_. Since my work has a focus on signal (and, more specifically image) processing, it's an excellent choice, as it provides a large body [...]
-I did all my programming on a Kubuntu_ system, running on an intel(R) Core (TM) i5-4570 CPU, and used GNU gcc_ 5.4.0 and clang_ 3.8.0 to compile the code in C++11 dialect. While I am confident that the code runs on other CPUs, I have not tested it with other compilers or operating systems (yet).
+I did all my programming on a Kubuntu_ system, running on an intel(R) Core (TM) i5-4570 CPU, and used GNU g++_ and clang++_ to compile the code in C++11 dialect. While I am confident that the code runs on other CPUs, I have not tested it with other compilers or operating systems (yet).
.. _here: http://bigwww.epfl.ch/thevenaz/interpolation/
.. _Vigra: http://ukoethe.github.io/vigra/
.. _Vc: https://compeng.uni-frankfurt.de/index.php?id=vc
.. _Kubuntu: http://kubuntu.org/
-.. _gcc: https://gcc.gnu.org/
-.. _clang: http://http://clang.llvm.org/
+.. _g++: https://gcc.gnu.org/
+.. _clang++: http://http://clang.llvm.org/
.. [CIT2000] Interpolation Revisited by Philippe Thévenaz, Member,IEEE, Thierry Blu, Member, IEEE, and Michael Unser, Fellow, IEEE in IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 19, NO. 7, JULY 2000, available online here_
@@ -132,7 +134,7 @@ Using double precision arithmetics, vectorization doesn't help so much, and pref
History
----------
-Some years ago, I needed uniform B-splines for a project in python. I looked for code in C which I could easily wrap with cffi_, as I intended to use it with pypy_, and found K. P. Esler's libeinspline_. I proceeded to code the wrapper, which also contained a layer to process Numpy_ nD-arrays, but at the time I did not touch the C code in libeinspline. The result of my efforts is still available from the repository_ I set up at the time. I did not use the code much and occupied myself wi [...]
+Some years ago, I needed uniform B-splines for a project in python. I looked for code in C which I could easily wrap with cffi_, as I intended to use it with pypy_, and found K. P. Esler's libeinspline_. I proceeded to code the wrapper, which also contained a layer to process Numpy_ nD-arrays, but at the time I did not touch the C code in libeinspline. The result of my efforts is still available from the repository_ I set up at the time. I did not use the code much and occupied myself wi [...]
.. _cffi: https://cffi.readthedocs.org/en/latest/
.. _pypy: http://pypy.org/
diff --git a/basis.h b/basis.h
index b1c84df..7573d45 100644
--- a/basis.h
+++ b/basis.h
@@ -163,11 +163,14 @@ real_type cdb_bspline_basis ( int x , int degree , int derivative = 0 )
}
/// see bspline_basis() below!
-///
-/// this helper routine works with the doubled value of x, so it can capture calls equivalent
-/// to basis ( x + .5 ) or basis ( x - .5 ) as basis2 ( x + 1 ) and basis2 ( x - 1 )
-/// having precalculated the basis function at .5 steps, we can therefore avoid
-/// using the general recursion formula. This is a big time-saver for high degrees.
+/// this helper routine works with the doubled value of x, so it can serve for calls
+/// equivalent to basis ( x + .5 ) or basis ( x - .5 ) as basis2 ( x + 1 ) and
+/// basis2 ( x - 1 ) having precalculated the basis function at .5 steps, we can
+/// therefore avoid using the general recursion formula. This is a big time-saver
+/// for high degrees. Note, though, that calculating the basis function for a
+/// spline's derivatives still needs recursion, with two branches per level.
+/// So calculating the basis function's value for high derivatives still consumes
+/// a fair amount of time.
template < class real_type >
real_type bspline_basis_2 ( int x2 , int degree , int derivative )
diff --git a/example/roundtrip.cc b/example/roundtrip.cc
index e18b29b..fc79c4c 100644
--- a/example/roundtrip.cc
+++ b/example/roundtrip.cc
@@ -152,16 +152,45 @@ void run_test ( view_type & data ,
// get a view to the core coefficients (those which aren't part of the brace)
view_type cfview = bsp.core ;
- // get the coordinate type from the evaluator
+ // set the coordinate type
typedef vigra::TinyVector < rc_type , 2 > coordinate_type ;
- // get evaluator type
+ // set the evaluator type
typedef vspline::evaluator<coordinate_type,pixel_type> eval_type ;
- // create the evaluator
- eval_type ev ( bsp ) ;
+ // create the evaluator. We create an evaluator using a mapping which corresponds
+ // to the boundary conditions we're using, instead of the default 'REJECT' mapping.
+ // For natural boundary conditions, there is no corresponding mapping and we use
+ // MAP_LIMIT instead.
+
+ vspline::map_code mc ;
+ switch ( bc )
+ {
+ case vspline::MIRROR:
+ mc = vspline::MAP_MIRROR ;
+ break ;
+ case vspline::NATURAL:
+ mc = vspline::MAP_LIMIT ;
+ break ;
+ case vspline::REFLECT:
+ mc = vspline::MAP_REFLECT ;
+ break ;
+ case vspline::PERIODIC:
+ mc = vspline::MAP_PERIODIC ;
+ break ;
+ default:
+ mc = vspline::MAP_REJECT ;
+ break ;
+ }
+
+ // create the evaluator for the b-spline, using plain evaluation (no derivatives)
+ // and the same mapping mode for both axes:
+
+ eval_type ev ( bsp , // spline
+ { 0 , 0 } , // (no) derivatives
+ { mc , mc } ) ; // mapping code as per switch above
- // get type for coordinate array
+ // type for coordinate array
typedef vigra::MultiArray<2, coordinate_type> coordinate_array ;
int Tx = Nx ;
@@ -253,11 +282,29 @@ void run_test ( view_type & data ,
<< " ms" << endl ;
#endif
+ cout << "difference original data/restored data:" << endl ;
check_diff<view_type> ( target , data ) ;
+ // fifth test: use 'restore' which internally delegates to grid_eval. This is
+ // usually slightly faster than the previous way to restore the original data,
+ // but otherwise makes no difference.
+
+#ifdef PRINT_ELAPSED
+ start = std::chrono::system_clock::now();
+#endif
+
+ for ( int times = 0 ; times < TIMES ; times++ )
+ vspline::restore < 2 , pixel_type , rc_type > ( bsp , target ) ;
+
+#ifdef PRINT_ELAPSED
+ end = std::chrono::system_clock::now();
+ cout << "avg " << TIMES << " x restore original data: .......... "
+ << std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count() / float(TIMES)
+ << " ms" << endl ;
+#endif
+
cout << "difference original data/restored data:" << endl ;
- vspline::restore_from_braced<view_type> ( bsp.coeffs , bsp.coeffs , DEGREE ) ;
- check_diff<view_type> ( cfview , data ) ;
+ check_diff<view_type> ( data , target ) ;
cout << endl ;
}
@@ -310,7 +357,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 < 4 ; spline_degree++ )
+ for ( int spline_degree = 0 ; spline_degree < 6 ; spline_degree++ )
{
#ifdef USE_VC
cout << "testing bc code " << vspline::bc_name[bc]
diff --git a/prefilter.h b/prefilter.h
index 259fea5..2cdc7ca 100644
--- a/prefilter.h
+++ b/prefilter.h
@@ -122,12 +122,6 @@
#include <array>
#include <assert.h>
-#include <vigra/multi_array.hxx>
-#include <vigra/multi_iterator.hxx>
-#include <vigra/navigator.hxx>
-#include <vigra/bordertreatment.hxx>
-#include <vigra/multi_convolution.hxx>
-
#include "common.h"
#include "filter.h"
#include "basis.h"
@@ -215,58 +209,6 @@ void solve ( input_array_type & input ,
njobs ) ;
}
-/// An interlude: restoration of the original knot point data from the spline coefficients.
-/// This is easily done by a simple convolution with the values of the basis function
-/// taken at discrete points inside the defined range.
-/// The function can take arbitrary spline degrees as template argument.
-/// there are two functions to restore the original data from a spline: the first one takes a
-/// braced spline and convolves it with BORDER_TREATMENT_AVOID. This way the explicit
-/// border treatment manifest in the brace is used, and splines with arbitrary border conditions
-/// can be verified.
-/// TODO: bit rough and ready - restoration from braced should really only produce the inner
-/// part, not the area covered by the brace. And it's neither multithreaded nor vectorized...
-
-template < class array >
-void restore_from_braced ( array &spline , array& target , int SplineDegree )
-{
- typedef typename array::value_type value_type ;
- const int half_ext = SplineDegree / 2 ;
- vigra::Kernel1D<value_type> spline_kernel ;
- spline_kernel.initExplicitly(-half_ext, half_ext) ;
- for ( int i = -half_ext ; i <= half_ext ; i++ )
- spline_kernel[i] = bspline_basis<double> ( i , SplineDegree ) ;
-// cout << "using kernel" ;
-// for ( int i = -half_ext ; i <= half_ext ; i++ )
-// cout << " " << spline_kernel[i] ;
-// cout << endl ;
- spline_kernel.setBorderTreatment ( BORDER_TREATMENT_AVOID ) ;
- separableConvolveMultiArray(spline, target, spline_kernel);
-}
-
-/// the second function takes a border treatment mode from vigra's collection and works
-/// on an unbraced spline. Some border treatment modes implemented here aren't available
-/// from vigra or aren't applicable
-/// - vigra vspline
-/// - BORDER_TREATMENT_AVOID used with braced splines
-/// - BORDER_TREATMENT_CLIP not sure
-/// - BORDER_TREATMENT_REPEAT nearest is flat or REFLECT, same effect in cubic splines
-/// - BORDER_TREATMENT_REFLECT MIRROR
-/// - BORDER_TREATMENT_WRAP PERIODIC
-/// - BORDER_TREATMENT_ZEROPAD -
-
-template < class array >
-void restore_from_unbraced ( array &spline , array& target , BorderTreatmentMode btm , int SplineDegree )
-{
- typedef typename array::value_type value_type ;
- const int half_ext = SplineDegree / 2 ;
- vigra::Kernel1D<value_type> spline_kernel ;
- spline_kernel.initExplicitly(-half_ext, half_ext) ;
- for ( int i = -half_ext ; i <= half_ext ; i++ )
- spline_kernel[i] = bspline_basis<double> ( i , SplineDegree ) ;
- spline_kernel.setBorderTreatment ( btm ) ;
- separableConvolveMultiArray(spline, target, spline_kernel);
-}
-
} ; // namespace vspline
#endif // VSPLINE_PREFILTER_H
diff --git a/remap.h b/remap.h
index 2c41975..3570645 100644
--- a/remap.h
+++ b/remap.h
@@ -102,6 +102,9 @@
/// While the code presented here is quite involved and there are several types and routines
/// the use(fulness) of which isn't immediately apparent, most use cases will be able to get
/// by using only remap() or index_remap().
+///
+/// Finally, this file also has a routine to restore the original knot point data from a
+/// bspline object. This is done very efficiently using grid_eval.
#ifndef VSPLINE_REMAP_H
#define VSPLINE_REMAP_H
@@ -1070,6 +1073,8 @@ void st_grid_eval ( shape_range_type < target_type::actual_dimension > range ,
/// to the prefiltering routine. Of course any other way of smoothing can
/// be used just the same, like a Burt filter or Gaussian smoothing.
+// TODO template args are too many and too unspecific, especially weight_type
+
template < typename evaluator_type , // type offering eval()
typename target_type ,
typename weight_type , // singular weight data type
@@ -1084,6 +1089,58 @@ void grid_eval ( rc_type ** const grid_coordinate ,
ncores * 8 , range , grid_coordinate , &itp , &result ) ;
}
+/// grid_eval allows us to code a function to restore the original knot point
+/// date from a bspline. We simply fill in the discrete coordinates into the
+/// grid coordinate vectors and call grid_eval with them.
+/// note that this routine can't operate in-place, so you can't overwrite
+/// a bspline object's core with the restored knot point data, you have to
+/// provide a separate target array.
+/// This routine is potentially faster than running an index_remap with
+/// the same target, due to the precalculated weight components.
+
+template < int dimension , typename value_type , typename rc_type = float >
+void restore ( const vspline::bspline < value_type , dimension > & bspl ,
+ vigra::MultiArrayView < dimension , value_type > & target )
+{
+ if ( target.shape() != bspl.core.shape() )
+ throw shape_mismatch
+ ( "restore: spline's core shape and target array shape must match" ) ;
+
+ typedef vigra::TinyVector < rc_type , dimension > coordinate_type ;
+ typedef vigra::MultiArrayView < dimension , value_type > target_type ;
+ typedef typename vigra::ExpandElementResult < value_type > :: type weight_type ;
+
+ // set up the coordinate component vectors
+ rc_type * p_ruler [ dimension ] ;
+ for ( int d = 0 ; d < dimension ; d++ )
+ {
+ p_ruler[d] = new rc_type [ target.shape ( d ) ] ;
+ for ( int i = 0 ; i < target.shape ( d ) ; i++ )
+ p_ruler[d][i] = rc_type(i) ;
+ }
+
+ // we specialize by 'evenness' to get evaluators using raw mapping for
+ // maximum speed.
+
+ if ( bspl.spline_degree & 1 )
+ {
+ typedef vspline::evaluator < coordinate_type , value_type , -1 , 0 > ev_type ;
+ ev_type ev ( bspl ) ;
+ vspline::grid_eval < ev_type , target_type , weight_type , rc_type >
+ ( p_ruler , ev , target ) ;
+ }
+ else
+ {
+ typedef vspline::evaluator < coordinate_type , value_type , -1 , 1 > ev_type ;
+ ev_type ev ( bspl ) ;
+ vspline::grid_eval < ev_type , target_type , weight_type , rc_type >
+ ( p_ruler , ev , target ) ;
+ }
+
+ for ( int d = 0 ; d < dimension ; d++ )
+ delete[] p_ruler[d] ;
+}
+
} ; // end of namespace vspline
#endif // VSPLINE_REMAP_H
diff --git a/vspline.doxy b/vspline.doxy
index 7412ac1..7767592 100644
--- a/vspline.doxy
+++ b/vspline.doxy
@@ -38,7 +38,7 @@ PROJECT_NAME = "vspline"
# could be handy for archiving the generated documentation or if some version
# control system is used.
-PROJECT_NUMBER = 16
+PROJECT_NUMBER = 17
# Using the PROJECT_BRIEF tag one can provide an optional one line description
# for a project that appears at the top of each page and should give viewer a
--
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