[Olena-patches] 3305: Make distance-geodesic-based routines have fastest versions.

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Thierry Geraud &lt;thierry.geraud(a)lrde.epita.fr&gt; Make distance-geodesic-based routines have fastest versions. * mln/canvas/distance_geodesic.hh (distance_geodesic_tests, distance_geodesic_dispatch): New. (distance_geodesic_fastest): New. (distance_geodesic): Call dispatch. * mln/transform/internal/influence_zone_functor.hh (init_p): New; fix missing. * mln/transform/internal/influence_zone_functor.hh, * mln/transform/internal/closest_point_functor.hh, * mln/transform/internal/distance_functor.hh, (P): Fix. (init_, init_p_, process_): New. (inqueue_p_wrt_input_p_, inqueue_p_wrt_input_n_): New. * mln/transform/closest_point_geodesic.hh: New. * mln/transform/all.hh: Update. * tests/transform/bench_closest_point_geodesic.cc: New. * tests/transform/closest_point_geodesic.cc: New. * tests/transform/Makefile.am: Update. mln/canvas/distance_geodesic.hh | 228 ++++++++++++++++++++++- mln/transform/all.hh | 1 mln/transform/closest_point_geodesic.hh | 81 ++++++++ mln/transform/internal/closest_point_functor.hh | 10 - mln/transform/internal/distance_functor.hh | 9 mln/transform/internal/influence_zone_functor.hh | 9 tests/transform/Makefile.am | 4 tests/transform/bench_closest_point_geodesic.cc | 63 ++++++ tests/transform/closest_point_geodesic.cc | 54 +++++ 9 files changed, 448 insertions(+), 11 deletions(-) Index: mln/transform/closest_point_geodesic.hh --- mln/transform/closest_point_geodesic.hh (revision 0) +++ mln/transform/closest_point_geodesic.hh (revision 0) @@ -0,0 +1,81 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +#ifndef MLN_TRANSFORM_CLOSEST_POINT_GEODESIC_HH +# define MLN_TRANSFORM_CLOSEST_POINT_GEODESIC_HH + +/// \file mln/transform/closest_point_geodesic.hh +/// +/// Geodesic closest point transform. +/// +/// \todo Add a version to retrieve both distance and closest point +/// maps. + +# include <mln/canvas/distance_geodesic.hh> +# include <mln/transform/internal/closest_point_functor.hh> + + + +namespace mln +{ + + namespace transform + { + + /// Discrete geodesic distance transform. + template <typename I, typename N, typename D> + mln_ch_value(I, mln_psite(I)) + closest_point_geodesic(const Image<I>& input, const Neighborhood<N>& nbh, D max); + + +# ifndef MLN_INCLUDE_ONLY + + template <typename I, typename N, typename D> + inline + mln_ch_value(I, mln_psite(I)) + closest_point_geodesic(const Image<I>& input, const Neighborhood<N>& nbh, D max) + { + trace::entering("transform::closest_point_geodesic"); + + mln_precondition(exact(input).is_valid()); + mln_precondition(exact(nbh).is_valid()); + + internal::closest_point_functor<I> f; + (void) mln::canvas::distance_geodesic(input, nbh, max, f); + + trace::exiting("transform::closest_point_geodesic"); + return f.cp_ima; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::transform + +} // end of namespace mln + + +#endif // ! MLN_TRANSFORM_CLOSEST_POINT_GEODESIC_HH Index: mln/transform/all.hh --- mln/transform/all.hh (revision 3304) +++ mln/transform/all.hh (working copy) @@ -42,6 +42,7 @@ } // end of namespace mln +# include <mln/transform/closest_point_geodesic.hh> # include <mln/transform/distance_front.hh> # include <mln/transform/distance_geodesic.hh> # include <mln/transform/influence_zone_front.hh> Index: mln/transform/internal/influence_zone_functor.hh --- mln/transform/internal/influence_zone_functor.hh (revision 3304) +++ mln/transform/internal/influence_zone_functor.hh (working copy) @@ -49,14 +49,21 @@ struct influence_zone_functor { typedef mln_value(I) V; - typedef mln_site(I) P; + typedef mln_psite(I) P; mln_concrete(I) output; void init(const I& input); bool inqueue_p_wrt_input_p(const V& input_p); bool inqueue_p_wrt_input_n(const V& input_n); + void init_p(const P&) {} // FIXME: move def below. void process(const P& p, const P& n); + + void init_(const I& input) { output = duplicate(input); } + bool inqueue_p_wrt_input_p_(const V& input_p) { return input_p != 0u; } + bool inqueue_p_wrt_input_n_(const V& input_n) { return input_n == 0u; } + void init_p_(unsigned) {} + void process_(unsigned p, unsigned n) { output.element(n) = output.element(p); } }; Index: mln/transform/internal/closest_point_functor.hh --- mln/transform/internal/closest_point_functor.hh (revision 3304) +++ mln/transform/internal/closest_point_functor.hh (working copy) @@ -48,7 +48,9 @@ struct closest_point_functor { typedef mln_value(I) V; - typedef mln_site(I) P; + typedef mln_psite(I) P; + + mln_ch_value(I,P) cp_ima; void init(const I&); bool inqueue_p_wrt_input_p(const V& input_p); @@ -56,7 +58,11 @@ bool inqueue_p_wrt_input_n(const V& input_n); void process(const P&, const P&); - mln_ch_value(I,P) cp_ima; + void init_(const I& input) { initialize(cp_ima, input); } + bool inqueue_p_wrt_input_p_(const V& input_p) { return input_p == true; } + void init_p_(unsigned p) { cp_ima.element(p) = cp_ima.point_at_index(p); } + bool inqueue_p_wrt_input_n_(const V& input_n) { return input_n == false; } + void process_(unsigned p, unsigned n) { cp_ima.element(n) = cp_ima.element(p); } }; Index: mln/transform/internal/distance_functor.hh --- mln/transform/internal/distance_functor.hh (revision 3304) +++ mln/transform/internal/distance_functor.hh (working copy) @@ -31,6 +31,8 @@ /// \file mln/transform/internal/distance_functor.hh /// /// Distance functor. +/// +/// \todo Move all functors into their corresponding file. # include <mln/core/macros.hh> @@ -56,6 +58,13 @@ void init_p(const P&); bool inqueue_p_wrt_input_n(const V& input_n); void process(const P&, const P&); + + + void init_(const I&) {} + bool inqueue_p_wrt_input_p_(const V& input_p) { return input_p == true; } + void init_p_(unsigned) {} + bool inqueue_p_wrt_input_n_(const V& input_n) { return input_n == false; } + void process_(unsigned, unsigned) {} }; Index: mln/canvas/distance_geodesic.hh --- mln/canvas/distance_geodesic.hh (revision 3304) +++ mln/canvas/distance_geodesic.hh (working copy) @@ -1,5 +1,4 @@ -// Copyright (C) 2008, 2009 EPITA Research and Development Laboratory -// (LRDE) +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) // // This file is part of the Olena Library. This library is free // software; you can redistribute it and/or modify it under the terms @@ -32,14 +31,16 @@ /// \file mln/canvas/distance_geodesic.hh /// /// Discrete geodesic distance canvas. -/// -/// \todo add fast version. Use extension(input) = true and extension(map) = 0. # include <mln/core/concept/image.hh> # include <mln/core/concept/neighborhood.hh> -# include <mln/core/site_set/p_queue_fast.hh> # include <mln/core/routine/duplicate.hh> + +# include <mln/core/site_set/p_queue_fast.hh> +# include <queue> + # include <mln/data/fill.hh> +# include <mln/extension/adjust_fill.hh> namespace mln @@ -56,20 +57,59 @@ F& functor); + # ifndef MLN_INCLUDE_ONLY + // Tests. + + namespace internal + { + + template <typename I, typename N, typename D, + typename F> + void + distance_geodesic_tests(const Image<I>& input_, const Neighborhood<N>& nbh_, D max, + F& functor) + { + const I& input = exact(input_); + const N& nbh = exact(nbh_); + + mln_precondition(input.is_valid()); + mln_precondition(nbh.is_valid()); + + (void) input; + (void) nbh; + (void) max; + (void) functor; + } + + + } // of namespace mln::canvas::internal + + + + // Implementations. + + namespace impl + { + + namespace generic + { + template <typename I, typename N, typename D, typename F> mln_ch_value(I, D) distance_geodesic(const Image<I>& input_, const Neighborhood<N>& nbh_, D max, F& functor) { - trace::entering("canvas::distance_geodesic"); + trace::entering("canvas::impl::generic::distance_geodesic"); const I& input = exact(input_); const N& nbh = exact(nbh_); + internal::distance_geodesic_tests(input, nbh, max, functor); + mln_precondition(input.is_valid()); mln_precondition(nbh.is_valid()); @@ -81,14 +121,17 @@ // Initialization. { + trace::entering("initialization"); + functor.init(input); // <-- init data::fill(dmap, max); + mln_piter(I) p(input.domain()); mln_niter(N) n(nbh, p); for_all(p) if (functor.inqueue_p_wrt_input_p(input(p))) // <-- inqueue_p_wrt_input_p { - functor.init_p(p); + functor.init_p(p); // <-- init_p dmap(p) = 0; for_all(n) if (input.domain().has(n) && @@ -98,10 +141,12 @@ break; } } + trace::exiting("initialization"); } // Propagation. { + trace::entering("propagation"); P p; mln_niter(N) n(nbh, p); while (! q.is_empty()) @@ -121,12 +166,179 @@ q.push(n); } } + trace::exiting("propagation"); } - trace::exiting("canvas::distance_geodesic"); + trace::exiting("canvas::impl::generic::distance_geodesic"); + return dmap; + } + + } // of namespace mln::canvas::impl::generic + + + + // Fastest version. + + template <typename I, typename N, typename D, + typename F> + mln_ch_value(I, D) + distance_geodesic_fastest(const Image<I>& input_, const Neighborhood<N>& nbh_, D max, + F& functor) + { + trace::entering("canvas::impl::distance_geodesic_fastest"); + + const I& input = exact(input_); + const N& nbh = exact(nbh_); + + internal::distance_geodesic_tests(input, nbh, max, functor); + + extension::adjust(input, nbh); + + mln_ch_value(I, D) dmap; // Distance map is aux data. + initialize(dmap, input); + + std::queue<unsigned> q; + + // Initialization. + { + trace::entering("initialization"); + + functor.init_(input); // <-- init + data::fill(dmap, max); + // For the extension to be ignored: + extension::fill(input, true); + extension::fill(dmap, D(0)); + + mln_pixter(const I) p(input); + mln_nixter(const I, N) n(p, nbh); + for_all(p) + if (functor.inqueue_p_wrt_input_p_(p.val())) // <-- inqueue_p_wrt_input_p + { + functor.init_p_(p); // <-- init_p + dmap.element(p.offset()) = 0; + for_all(n) + if (functor.inqueue_p_wrt_input_n_(n.val())) // <-- inqueue_p_wrt_input_n + { + q.push(p.offset()); + break; + } + } + + trace::exiting("initialization"); + } + + // Propagation. + { + trace::entering("propagation"); + unsigned p; + + util::array<int> dp = offsets_wrt(input, nbh); + const unsigned n_nbhs = dp.nelements(); + + while (! q.empty()) + { + p = q.front(); + q.pop(); + if (dmap.element(p) == max) + // Saturation so stop. + break; + for (unsigned i = 0; i < n_nbhs; ++i) + { + unsigned n = p + dp[i]; + if (dmap.element(n) == max) + { + dmap.element(n) = dmap.element(p) + 1; + functor.process_(p, n); // <- process + q.push(n); + } + } + } + trace::exiting("propagation"); + } + + trace::exiting("canvas::impl::distance_geodesic_fastest"); return dmap; } + + } // of namespace mln::canvas::impl + + + + // Dispatch. + + namespace internal + { + + template <typename I, typename N, typename D, + typename F> + inline + mln_ch_value(I, D) + distance_geodesic_dispatch(metal::false_, + const Image<I>& input, const Neighborhood<N>& nbh, D max, + F& functor) + { + return impl::generic::distance_geodesic(input, nbh, max, + functor); + } + + template <typename I, typename N, typename D, + typename F> + inline + mln_ch_value(I, D) + distance_geodesic_dispatch(metal::true_, + const Image<I>& input, const Neighborhood<N>& nbh, D max, + F& functor) + { + return impl::distance_geodesic_fastest(input, nbh, max, functor); +// return impl::generic::distance_geodesic(input, nbh, max, +// functor); + } + + template <typename I, typename N, typename D, + typename F> + inline + mln_ch_value(I, D) + distance_geodesic_dispatch(const Image<I>& input, const Neighborhood<N>& nbh, D max, + F& functor) + { + enum { + test = mlc_equal(mln_trait_image_speed(I), + trait::image::speed::fastest)::value + && + mln_is_simple_neighborhood(N)::value + }; + return distance_geodesic_dispatch(metal::bool_<test>(), + input, nbh, max, + functor); + } + + + } // of namespace mln::canvas::internal + + + + // Facade. + + template <typename I, typename N, typename D, + typename F> + inline + mln_ch_value(I, D) + distance_geodesic(const Image<I>& input, const Neighborhood<N>& nbh, D max, + F& functor) + { + trace::entering("canvas::distance_geodesic"); + + internal::distance_geodesic_tests(input, nbh, max, functor); + + mln_ch_value(I,D) output; + output = internal::distance_geodesic_dispatch(input, nbh, max, functor); + + trace::exiting("canvas::distance_geodesic"); + return output; + } + + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::canvas Index: tests/transform/Makefile.am --- tests/transform/Makefile.am (revision 3304) +++ tests/transform/Makefile.am (working copy) @@ -3,11 +3,15 @@ include $(top_srcdir)/milena/tests/tests.mk check_PROGRAMS = \ + bench_closest_point_geodesic \ + closest_point_geodesic \ distance_front \ distance_geodesic \ influence_zone_front \ influence_zone_geodesic +bench_closest_point_geodesic_SOURCES = bench_closest_point_geodesic.cc +closest_point_geodesic_SOURCES = closest_point_geodesic.cc distance_front_SOURCES = distance_front.cc distance_geodesic_SOURCES = distance_geodesic.cc influence_zone_front_SOURCES = influence_zone_front.cc Index: tests/transform/bench_closest_point_geodesic.cc --- tests/transform/bench_closest_point_geodesic.cc (revision 0) +++ tests/transform/bench_closest_point_geodesic.cc (revision 0) @@ -0,0 +1,63 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +/// \file tests/transform/bench_closest_point_geodesic.cc +/// +/// Test on mln::transform::closest_point_geodesic. + +#include <cstdlib> + +#include <mln/core/image/image3d.hh> +#include <mln/core/alias/neighb3d.hh> +#include <mln/data/fill.hh> +#include <mln/opt/at.hh> +#include <mln/transform/closest_point_geodesic.hh> + + +int main() +{ + using namespace mln; + using value::int_u8; + + const unsigned + nslis = 100, + nrows = 250, + ncols = 250; + image3d<bool> input(nslis, nrows, ncols); + data::fill(input, false); + for (unsigned i = 0; i < 100; ++i) + opt::at(input, + std::rand() % nslis, + std::rand() % nrows, + std::rand() % ncols) = true; + + trace::quiet = false; + + image3d<point3d> output = transform::closest_point_geodesic(input, + c6(), + mln_max(unsigned)); +} Index: tests/transform/closest_point_geodesic.cc --- tests/transform/closest_point_geodesic.cc (revision 0) +++ tests/transform/closest_point_geodesic.cc (revision 0) @@ -0,0 +1,54 @@ +// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// +// This file is part of the Olena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +/// \file tests/transform/closest_point_geodesic.cc +/// +/// Test on mln::transform::closest_point_geodesic. + +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/data/fill.hh> +#include <mln/debug/println.hh> +#include <mln/opt/at.hh> + +#include <mln/transform/closest_point_geodesic.hh> + + +int main() +{ + using namespace mln; + using value::int_u8; + + image2d<bool> input(9, 9); + data::fill(input, false); + opt::at(input, 2, 4) = true; + opt::at(input, 4, 2) = true; + + image2d<point2d> output = transform::closest_point_geodesic(input, c4(), int_u8(4)); + debug::println(output); +}