This is an automated email from the git hooks/post-receive script. gert-guest pushed a commit to branch master in repository mia.
commit d47bc63994d9952c8f7aacea66b809d801a28043 Author: Gert Wollny <[email protected]> Date: Fri Jan 17 12:59:35 2014 +0100 Imported Upstream version 2.0.13 --- CMakeLists.txt | 6 ++-- ChangeLog | 15 ++++++++ mia/3d/fifof/label.cc | 12 +++---- mia/3d/interpolator.cxx | 34 ++++++++++++------ mia/3d/interpolator.hh | 5 +-- mia/3d/test_interpol.cc | 80 +++++++++++++++++++++++++++++++++++++++++++ mia/3d/transform.cc | 3 +- mia/core/factory.hh | 3 +- mia/core/iohandler.cxx | 5 +-- mia/core/splinekernel.cc | 16 ++++++++- mia/core/splinekernel.hh | 2 ++ mia/core/test_splinekernel.cc | 29 ++++++++++++++++ src/2dimageregistration.cc | 2 ++ src/3dgetslice.cc | 17 +++++---- 14 files changed, 196 insertions(+), 33 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 26a70bf..5268763 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -41,9 +41,9 @@ SET(VENDOR "Gert Wollny") SET(PACKAGE_NAME "mia") SET(MAJOR_VERSION 2) SET(MINOR_VERSION 0) -SET(MICRO_VERSION 12) -SET(INTERFACE_AGE 4) -SET(BINARY_AGE 4) +SET(MICRO_VERSION 13) +SET(INTERFACE_AGE 5) +SET(BINARY_AGE 5) SET(PACKAGE_VERSION "${MAJOR_VERSION}.${MINOR_VERSION}.${MICRO_VERSION}") SET(VERSION "${MAJOR_VERSION}.${MINOR_VERSION}") diff --git a/ChangeLog b/ChangeLog index 229385c..dc4f3d7 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,18 @@ +2.0.13 + + Fixes for tracked bugs: + + * #138 Handle compressed file extension properly when selecting the IO plugin. + * #143 Fix zero-fill boundary condition handling + * #144 Parallelize the prefiltering in the 3D interpolator + + Additional fixes: + + * 2dimageregistration: Exit if no cost functions are given + * 3dgetimage: add support for setting the minimum number of + digits in the output file names. + * factoryplugins: throw error if a factory plugin is not found + * make the basic interpolator operator in 3D thread save 2.0.12 diff --git a/mia/3d/fifof/label.cc b/mia/3d/fifof/label.cc index cff4502..01bb233 100644 --- a/mia/3d/fifof/label.cc +++ b/mia/3d/fifof/label.cc @@ -88,8 +88,8 @@ void C2DLabelStackFilter::grow( int x, int y, C2DBitImage& input, unsigned short void C2DLabelStackFilter::label_new_regions(C2DBitImage& input) { - C3DBitImage::iterator ii = input.begin(); - C3DUSImage::iterator usi = m_out_buffer.begin(); + auto ii = input.begin(); + auto usi = m_out_buffer.begin(); for (size_t y = 0; y < input.get_size().y; ++y) for (size_t x = 0; x < input.get_size().x; ++x, ++usi, ++ii) { @@ -112,7 +112,7 @@ void C2DLabelStackFilter::label_new_regions(C2DBitImage& input) void C2DLabelStackFilter::label(C2DBitImage& input) { // first grow all regions that are already labeled from the last slice - C3DUSImage::iterator usi = m_out_buffer.begin(); + auto usi = m_out_buffer.begin(); for (size_t y = 0; y < input.get_size().y; ++y) for (size_t x = 0; x < input.get_size().x; ++x, ++usi) { if ( *usi ) @@ -132,9 +132,9 @@ void C2DLabelStackFilter::new_label(C2DBitImage& input) void C2DLabelStackFilter::re_label(C2DBitImage& input) { - C3DUSImage::iterator usi = m_out_buffer.begin(); - C3DUSImage::iterator use = m_out_buffer.end(); - C2DBitImage::iterator ii = input.begin(); + auto usi = m_out_buffer.begin(); + auto use = m_out_buffer.end(); + auto ii = input.begin(); // maintain old labeling for new slice, and clean the input at // labeled positions diff --git a/mia/3d/interpolator.cxx b/mia/3d/interpolator.cxx index fba7e77..ac7ec39 100644 --- a/mia/3d/interpolator.cxx +++ b/mia/3d/interpolator.cxx @@ -20,6 +20,10 @@ #include <cmath> +#include <mia/core/threadedmsg.hh> +#include <tbb/parallel_for.h> +#include <tbb/blocked_range.h> + NS_MIA_BEGIN template <class T> @@ -117,38 +121,44 @@ void T3DConvoluteInterpolator<T>::prefilter(const T3DDatafield<T>& image) int cachYSize = image.get_size().y; int cachZSize = image.get_size().z; - { + + auto filter_x = [this, cachXSize, cachYSize, poles](const tbb::blocked_range<size_t>& range_z) { coeff_vector buffer(cachXSize); - for (int z = 0; z < cachZSize; z++){ + for (auto z = range_z.begin(); z != range_z.end() ; ++z){ for (int y = 0; y < cachYSize; y++) { m_coeff.get_data_line_x(y,z,buffer); m_xbc->filter_line(buffer, poles); m_coeff.put_data_line_x(y,z,buffer); } } - } + }; + parallel_for(tbb::blocked_range<size_t>(0, cachZSize, 1), filter_x); - { + + auto filter_y = [this, cachXSize, cachYSize, poles](const tbb::blocked_range<size_t>& range_z) { coeff_vector buffer(cachYSize); - for (int z = 0; z < cachZSize; z++){ + for (auto z = range_z.begin(); z != range_z.end() ; ++z){ for (int x = 0; x < cachXSize; x++) { m_coeff.get_data_line_y(x,z,buffer); m_ybc->filter_line(buffer, poles); m_coeff.put_data_line_y(x,z,buffer); } } - } + }; + parallel_for(tbb::blocked_range<size_t>(0, cachZSize, 1), filter_y); - { + + auto filter_z = [this, cachXSize, cachZSize, poles](const tbb::blocked_range<size_t>& range_y) { coeff_vector buffer(cachZSize); - for (int y = 0; y < cachYSize; y++){ + for (auto y = range_y.begin(); y != range_y.end() ; ++y){ for (int x = 0; x < cachXSize; x++) { m_coeff.get_data_line_z(x,y,buffer); m_zbc->filter_line(buffer, poles); m_coeff.put_data_line_z(x,y,buffer); } } - } + }; + parallel_for(tbb::blocked_range<size_t>(0, cachYSize, 1), filter_z); } @@ -207,7 +217,8 @@ struct add_3d<T3DDatafield< T >, 1> { const CSplineKernel::SCache& yc, const CSplineKernel::SCache& zc) { - return coeff(xc.index[0], yc.index[0], zc.index[0] ) ; + return xc.weights[0] * yc.weights[0] * zc.weights[0] * + coeff(xc.index[0], yc.index[0], zc.index[0] ) ; } }; @@ -267,7 +278,7 @@ T T3DConvoluteInterpolator<T>::operator () (const C3DFVector& x, CWeightCache& m_kernel->get_cached(x.y, cache.y); if (x.z != cache.z.x) - m_kernel->get_cached(x.z, cache.z); + m_kernel->get_cached(x.z, cache.z); U result = U(); // now we give the compiler a chance to optimize based on kernel size and data type. @@ -294,6 +305,7 @@ T T3DConvoluteInterpolator<T>::operator () (const C3DFVector& x, CWeightCache& template <typename T> T T3DConvoluteInterpolator<T>::operator () (const C3DFVector& x) const { + CScopedLock lock(m_cache_lock); typedef typename TCoeff3D::value_type U; // x will usually be the fastest changing index, therefore, it is of no use to use the cache diff --git a/mia/3d/interpolator.hh b/mia/3d/interpolator.hh index caa65f5..4b2f838 100644 --- a/mia/3d/interpolator.hh +++ b/mia/3d/interpolator.hh @@ -26,7 +26,7 @@ #include <mia/core/splinekernel.hh> #include <mia/core/boundary_conditions.hh> #include <mia/3d/image.hh> - +#include <tbb/mutex.h> NS_MIA_BEGIN @@ -154,7 +154,8 @@ private: T m_min; T m_max; - + + mutable tbb::mutex m_cache_lock; mutable CSplineKernel::SCache m_x_cache; mutable CSplineKernel::SCache m_y_cache; mutable CSplineKernel::SCache m_z_cache; diff --git a/mia/3d/test_interpol.cc b/mia/3d/test_interpol.cc index c7a746e..501781e 100644 --- a/mia/3d/test_interpol.cc +++ b/mia/3d/test_interpol.cc @@ -28,6 +28,7 @@ #include <mia/3d/interpolator.hh> #include <mia/core/msgstream.hh> +#include <mia/core/threadedmsg.hh> NS_MIA_USE using namespace std; @@ -197,6 +198,45 @@ struct FParallelInterpolator { }; +struct FParallelInterpolator2 { + T3DDatafield<float>& output; + const T3DConvoluteInterpolator<float>& src; + const C3DFVector& shift; + const T3DDatafield<float>& test_data; + + FParallelInterpolator2(T3DDatafield<float>& _output, + const T3DConvoluteInterpolator<float>& _src, + const C3DFVector& _shift, + const T3DDatafield<float>& _test_data + + ): + output(_output), + src(_src), + shift(_shift), + test_data(_test_data) + { + } + + void operator()( const tbb::blocked_range<int>& range ) const { + CThreadMsgStream thread_stream; + auto cache = src.create_cache(); + for (auto z = range.begin(); z != range.end(); ++z) + for (size_t y = 0; y < output.get_size().y; ++y) + for (size_t x = 0; x < output.get_size().x; ++x) { + C3DFVector loc(x,y,z); + output(x,y,z) = src(loc + shift, cache); + cvmsg() << loc << " weights:" + << " x = " << cache.x.weights << cache.x.index + << " y = " << cache.y.weights << cache.y.index + << " z = " << cache.z.weights << cache.z.index + << " res= " << output(x,y,z) + << " test= " << test_data(x,y,z) + << "\n"; + } + + } +}; + static void test_parallel_interpolator() { T3DDatafield<float> data(C3DBounds(10, 12, 11)); @@ -222,6 +262,44 @@ static void test_parallel_interpolator() } +static void test_parallel_interpolator_zerofill_shifted() +{ + C3DFVector shift(1,2,3); + + T3DDatafield<float> data(C3DBounds(9, 8, 7)); + T3DDatafield<float> test_data(C3DBounds(9, 8, 7)); + auto kernel = produce_spline_kernel(bspline0); + auto bc = produce_spline_boundary_condition("zero"); + + auto i = data.begin(); + for (size_t z = 0; z < data.get_size().z; ++z) + for (size_t y = 0; y < data.get_size().y; ++y) + for (size_t x = 0; x < data.get_size().x; ++x, ++i) { + *i = x + 10 * y + 100 * z; + if (x >= 1 && y >= 2 && z >= 3) + test_data(x - 1, y - 2, z - 3) = *i; + } + + T3DConvoluteInterpolator<float> src(data, kernel, *bc, *bc, *bc); + + T3DDatafield<float> output(data.get_size()); + FParallelInterpolator2 worker(output, src, shift, test_data); + + tbb::parallel_for(tbb::blocked_range<int>( 0, data.get_size().z), worker); + for (size_t z = 0; z < data.get_size().z; ++z) + for (size_t y = 0; y < data.get_size().y; ++y) + for (size_t x = 0; x < data.get_size().x; ++x) { + if (output(x,y,z) != test_data(x,y,z)) { + cvfail() << "FAIL:" << x<< "," << y << "," << z << ": " + << output(x,y,z) << " vs "<< test_data(x,y,z) << "\n"; + } + BOOST_CHECK_CLOSE(output(x,y,z), test_data(x,y,z), 0.01); + } + +} + + + static double omoms3(double x) { if (x < 0) @@ -336,4 +414,6 @@ void add_3dinterpol_tests( boost::unit_test::test_suite* suite) suite->add( BOOST_TEST_CASE( (&test_external_cache_interpolator))); suite->add( BOOST_TEST_CASE( (&test_parallel_interpolator))); + suite->add( BOOST_TEST_CASE( (&test_parallel_interpolator_zerofill_shifted))); + } diff --git a/mia/3d/transform.cc b/mia/3d/transform.cc index 6d7da39..f566d5c 100644 --- a/mia/3d/transform.cc +++ b/mia/3d/transform.cc @@ -302,9 +302,10 @@ void F3DTransformer<T>::operator() ( const tbb::blocked_range<int>& range ) cons CThreadMsgStream thread_stream; auto cache = interp.create_cache(); - auto r = result.begin_at(0,0,range.begin()); C3DBounds begin(0,0,range.begin()); C3DBounds end(result.get_size().x,result.get_size().y, range.end()); + + auto r = result.begin_at(0,0,range.begin()); cvdebug() << "range = " << begin << " - " << end << "\n"; diff --git a/mia/core/factory.hh b/mia/core/factory.hh index 9123ad4..c9d850b 100644 --- a/mia/core/factory.hh +++ b/mia/core/factory.hh @@ -272,7 +272,8 @@ typename I::Product *TFactoryPluginHandler<I>::produce_raw(const std::string& pa cvdebug() << "TFactoryPluginHandler<>::produce: Create plugin from '" << factory_name << "'\n"; auto factory = this->plugin(factory_name.c_str()); - DEBUG_ASSERT_RELEASE_THROW(factory, "A plug-in was not found but 'this->plugin' did not throw"); + if (!factory) + throw create_exception<std::invalid_argument>("Unable to find plugin for '", factory_name.c_str(), "'"); return factory->create(param_list.begin()->second,params.c_str()); } diff --git a/mia/core/iohandler.cxx b/mia/core/iohandler.cxx index 3ea832c..ac08d2c 100644 --- a/mia/core/iohandler.cxx +++ b/mia/core/iohandler.cxx @@ -57,9 +57,10 @@ TIOPluginHandler<I>::preferred_plugin_ptr(const std::string& fname) const std::string fsuffix = __bfs_get_extension(fpath); if (!fsuffix.empty()) { if (m_compress_sfx.find(fsuffix) != m_compress_sfx.end()) { + cvdebug() << "Got compression suffix '" << fsuffix << "\n"; // remove the last extension and get the one before bfs::path help(fpath.stem()); - fsuffix = __bfs_get_extension(fpath); + fsuffix = __bfs_get_extension(help); } }else fsuffix = fname; @@ -123,7 +124,7 @@ TIOPluginHandler<I>::preferred_plugin(const std::string& fname) const auto fsuffix = __bfs_get_extension(fpath); if (m_compress_sfx.find(fsuffix) != m_compress_sfx.end()) { bfs::path help(fpath.stem()); - fsuffix = __bfs_get_extension(fpath); + fsuffix = __bfs_get_extension(help); } cvdebug() << "Got suffix '" << fsuffix << "'\n"; auto p = m_suffixmap.find(fsuffix); diff --git a/mia/core/splinekernel.cc b/mia/core/splinekernel.cc index 205dae8..14013ae 100644 --- a/mia/core/splinekernel.cc +++ b/mia/core/splinekernel.cc @@ -123,8 +123,17 @@ void CSplineKernel::get_cached(double x, SCache& cache)const { int start_idx = get_start_idx_and_value_weights(x, cache.weights); cache.x = x; - if (start_idx == cache.start_idx) + + // here we need an additional field that stores whether the + // indices crossed a boundary + + // start index is same and inside the image, we can keep the new weights + // only zero-boundary conditions must check if weights were changed + if (start_idx == cache.start_idx && // cache.boundary_condition.is_zero() + start_idx >= 0 && + cache.start_idx <= cache.index_limit ) return; + cache.start_idx = start_idx; cache.is_flat = false; fill_index(start_idx, cache.index); @@ -237,6 +246,11 @@ int CSplineKernel::get_start_idx_and_value_weights(double x, VWeight& weights) c return result - (int)m_half_degree; } +int CSplineKernel::get_start_idx(double x) const +{ + return fastfloor(x + m_shift) - (int)m_half_degree; +} + int CSplineKernel::get_start_idx_and_derivative_weights(double x, VWeight& weights) const { const int result = fastfloor(x + m_shift); diff --git a/mia/core/splinekernel.hh b/mia/core/splinekernel.hh index fe3f8a1..e8c0ee1 100644 --- a/mia/core/splinekernel.hh +++ b/mia/core/splinekernel.hh @@ -262,6 +262,8 @@ protected: void add_pole(double x); private: + int get_start_idx(double x) const; + /** Helper function to fill the array index with consecutive values starting with i */ diff --git a/mia/core/test_splinekernel.cc b/mia/core/test_splinekernel.cc index 676233b..5ac1405 100644 --- a/mia/core/test_splinekernel.cc +++ b/mia/core/test_splinekernel.cc @@ -79,6 +79,35 @@ BOOST_FIXTURE_TEST_CASE( test_mirrored_big, TestSplineCacheFixture ) } +BOOST_AUTO_TEST_CASE( test_zero_cache_outside_never_flat ) +{ + auto bc = produce_spline_boundary_condition("zero"); + bc->set_width(10); + + CBSplineKernelMock skm; + + CSplineKernel::SCache cache(skm.size(), *bc, true); + skm.get_cached(9, cache); + + cvdebug() << "cache {" + << " x = " << cache.x + << " start_idx = " << cache.start_idx + << " is_flat = " << cache.is_flat + << " index = " << cache.index + << "\n"; + + skm.get_cached(9, cache); + BOOST_CHECK_EQUAL(cache.weights[0], 0); + BOOST_CHECK_EQUAL(cache.weights[1], 1); + BOOST_CHECK_EQUAL(cache.weights[2], 0); + + BOOST_CHECK_EQUAL(cache.index[0], 8); + BOOST_CHECK_EQUAL(cache.index[1], 9); + BOOST_CHECK_EQUAL(cache.index[2], 0); + +} + + CBSplineKernelMock::CBSplineKernelMock(): CSplineKernel(3, 0, ip_bspline3) { diff --git a/src/2dimageregistration.cc b/src/2dimageregistration.cc index 782664b..f9fb769 100644 --- a/src/2dimageregistration.cc +++ b/src/2dimageregistration.cc @@ -71,6 +71,8 @@ int do_main( int argc, char *argv[] ) auto cost_descrs = options.get_remaining(); + if (cost_descrs.empty()) + throw runtime_error("No registration criterion given"); C2DFullCostList costs; for (auto i = cost_descrs.begin(); i != cost_descrs.end(); ++i) diff --git a/src/3dgetslice.cc b/src/3dgetslice.cc index 9bb0fbe..57eb1bc 100644 --- a/src/3dgetslice.cc +++ b/src/3dgetslice.cc @@ -98,11 +98,12 @@ struct __dispatch<T, dir_yz> { template <EDirection s_dir> class TGetter : public TFilter<bool> { public: - TGetter(size_t start, size_t n, const string& fname, const string& type): + TGetter(size_t start, size_t n, const string& fname, const string& type, int digits): m_start(start), m_n(n), m_fname(fname), - m_type(type) + m_type(type), + m_digits(digits) { } @@ -115,7 +116,7 @@ public: for(size_t i = m_start; i < end; ++i) { P2DImage pimage(new T2DImage<T>(__dispatch<T, s_dir>::get_slice(i, image))); stringstream out_name; - out_name << m_fname << setw(4) << setfill('0') << i << "." << m_type; + out_name << m_fname << setw(m_digits) << setfill('0') << i << "." << m_type; retval &= save_image(out_name.str(), pimage); } return retval; @@ -125,6 +126,7 @@ private: size_t m_n; string m_fname; string m_type; + int m_digits; }; int do_main( int argc, char *argv[] ) @@ -136,6 +138,7 @@ int do_main( int argc, char *argv[] ) size_t start_slice = 0; size_t slice_number = 0; EDirection direction = dir_xy; + int digits = 4; const auto& imageio2d = C2DImageIOPluginHandler::instance(); @@ -149,6 +152,8 @@ int do_main( int argc, char *argv[] ) options.add(make_opt( out_type, imageio2d.get_set(), "type", 't', "output file type")); options.add(make_opt( start_slice, "start", 's',"start slice number")); options.add(make_opt( slice_number, "number", 'n', "number of slices (all=0)")); + options.add(make_opt( digits, "ndigits", 0, "minimum number of digits of the file name numbers")); + options.add(make_opt( direction, GDirectionmap, "dir", 'd', "slice direction (xy=axial, xz=coronal, yz=saggital)")); @@ -166,17 +171,17 @@ int do_main( int argc, char *argv[] ) switch (direction) { case dir_xy: result = mia::filter(TGetter<dir_xy>(start_slice, slice_number, - out_filename, out_suffix), + out_filename, out_suffix, digits), *in_image); break; case dir_xz: result = mia::filter(TGetter<dir_xz>(start_slice, slice_number, - out_filename, out_suffix), + out_filename, out_suffix, digits), *in_image); break; case dir_yz: result = mia::filter(TGetter<dir_yz>(start_slice, slice_number, - out_filename, out_suffix), + out_filename, out_suffix, digits), *in_image); break; default: -- Alioth's /git/debian-med/git-commit-notice on /srv/git.debian.org/git/debian-med/mia.git _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
