February 2009 - Olena-patches

r3309: Implement fastest regional_minima
by Fabien Freling 06 Feb '09

06 Feb '09

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-02-06 Fabien Freling <freling(a)lrde.epita.fr> Implement fastest regional_minima. * fabien/level.hh: New file to implement fastest level. * fabien/regional_maxima.hh: Fix. * fabien/regional_minima.cc: New file to test regional_minima. * fabien/regional_minima.hh: Fastest implementation. --- level.hh | 188 +++++++++++++++++++++++++++++++++++++++++++++++++++++ regional_minima.cc | 57 ++++++++++++++++ regional_minima.hh | 172 ++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 417 insertions(+) Index: trunk/milena/sandbox/fabien/level.hh =================================================================== --- trunk/milena/sandbox/fabien/level.hh (revision 0) +++ trunk/milena/sandbox/fabien/level.hh (revision 3309) @@ -0,0 +1,188 @@ +// Copyright (C) 2007, 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_LABELING_LEVEL_HH +# define MLN_LABELING_LEVEL_HH + +/// \file mln/labeling/level.hh +/// +/// Connected component labeling of the image objects at a given +/// level. + +# include <mln/core/concept/image.hh> +# include <mln/core/concept/neighborhood.hh> +# include <mln/canvas/labeling.hh> +# include <mln/data/fill.hh> + + + +namespace mln +{ + + namespace labeling + { + + /*! Connected component labeling of the image objects at a given + * level. + * + * \param[in] input The input image. + * \param[in] val The level to consider for the labeling. + * \param[in] nbh The connexity of the level components. + * \param[out] nlabels The number of labels. + * \return The label image. + */ + template <typename I, typename N, typename L> + mln_ch_value(I, L) + level(const Image& input, const mln_value(I)& val, + const Neighborhood<N>& nbh, L& nlabels); + + +# ifndef MLN_INCLUDE_ONLY + + + // Tests. + + namespace internal + { + + template <typename I, typename N, typename L> + void + level_tests(const Image& input, const mln_value(I)& val, const Neighborhood<N>& nbh, + L& nlabels) + { + mln_precondition(exact(input).is_valid()); + // mln_precondition(exact(nbh).is_valid()); + + (void) input; + (void) val; + (void) nbh; + (void) nlabels; + } + + } // end of namespace mln::labeling::internal + + + + // Generic implementation. + + namespace impl + { + + + struct labeling_functor_base + { + void init() {} + + template <typename P> + bool handles(const P&) const { return true; } + + template <typename L, typename R> + bool equiv(const L&, const R&) const { return false; } + + template <typename P> + bool labels(const P&) const { return true; } + + template <typename L, typename R> + void do_no_union(const L&, const R&) {} + + template <typename P> + void init_attr(const P&) {} + + template <typename L, typename R> + void merge_attr(const L&, const R&) {} + }; + + + // Generic functor. + + template <typename I> + struct level_functor : labeling_functor_base + { + typedef mln_psite(I) P; + + const I& input; + const mln_value(I)& val; + + // Requirements from mln::canvas::labeling. + + typedef mln_pset(I) S; + + // Generic implementation + + void init() {} + bool handles(const P& p) const { return input(p) == val; } + bool equiv(const P& n, const P&) const { return input(n) == val; } + bool labels(const P&) const { return true; } + + // Fastest implementation + + void init_() {} + bool handles_(const P& p) const { return input.element(p) == val; } + bool equiv_(const P& n, const P&) const { return input.element(n) == val; } + bool labels_(const P&) const { return true; } + + + // end of Requirements. + + level_functor(const Image& input_, const mln_value(I)& val) + : input(exact(input_)), + val(val) + { + } + }; + + + + + // Facade. + + template <typename I, typename N, typename L> + mln_ch_value(I, L) + level(const Image& input, const mln_value(I)& val, const Neighborhood<N>& nbh, + L& nlabels) + { + trace::entering("labeling::level"); + + internal::level_tests(input, val, nbh, nlabels); + + mln_ch_value(I, L) output; + level_functor f(input, val); + output = canvas::labeling_video(input, nbh, nlabels, f); + + trace::exiting("labeling::level"); + return output; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::labeling + +} // end of namespace mln + + +#endif // ! MLN_LABELING_LEVEL_HH Index: trunk/milena/sandbox/fabien/regional_minima.hh =================================================================== --- trunk/milena/sandbox/fabien/regional_minima.hh (revision 0) +++ trunk/milena/sandbox/fabien/regional_minima.hh (revision 3309) @@ -0,0 +1,172 @@ +// Copyright (C) 2007, 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_LABELING_REGIONAL_MINIMA_HH +# define MLN_LABELING_REGIONAL_MINIMA_HH + +/// \file mln/labeling/regional_minima.hh +/// +/// Connected component labeling of the regional minima of an image. + +# include <mln/core/concept/image.hh> +# include <mln/core/concept/neighborhood.hh> + +# include "labeling.hh" + +# include <mln/data/fill.hh> +# include <mln/level/sort_psites.hh> + + +namespace mln +{ + + namespace labeling + { + + /*! Connected component labeling of the regional minima of an + * image. + * + * \param[in] input The input image. + * \param[in] nbh The connexity of the regional minima. + * \param[out] nlabels The number of labeled regions. + * \return The label image. + * + */ + template <typename I, typename N, typename L> + mln_ch_value(I, L) + regional_minima(const Image& input, const Neighborhood<N>& nbh, + L& nlabels); + + +# ifndef MLN_INCLUDE_ONLY + + namespace impl + { + + // Generic functor. + + template <typename I> + struct regional_minima_functor + { + typedef mln_psite(I) P; + + // requirements from mln::canvas::labeling: + + const I& input; + + // Generic implementation + + void init() { data::fill(attr, true); } + bool handles(const P&) const { return true; } + bool labels(const P& p) const { return attr(p); } + bool equiv(const P& n, const P& p) const { return input(n) == + input(p); } + void do_no_union(const P& n, const P& p) + { + // Avoid a warning about an undefined variable when NDEBUG + // is not defined. + (void)n; + + mln_invariant(input(n) < input(p)); + attr(p) = false; + } + + void init_attr(const P&) {} + void merge_attr(const P& r, const P& p) { attr(p) = attr(p) && + attr(r); } + + // Fastest implementation + + void init_() { data::fill(attr, true); } + bool handles_(unsigned p) const { return true; } + bool labels_(unsigned p) const { return attr.element(p); } + bool equiv_(unsigned n, unsigned p) const { return input.element(n) == + input.element(p); } + void do_no_union_(unsigned n, unsigned p) + { + // Avoid a warning about an undefined variable when NDEBUG + // is not defined. + (void)n; + + mln_invariant(input.element(n) < input.element(p)); + attr.element(p) = false; + } + + void init_attr_(unsigned) {} + void merge_attr_(unsigned r, unsigned p) { attr.element(p) = attr.element(p) && + attr.element(r); } + + // end of requirements + + mln_ch_value(I, bool) attr; + + regional_minima_functor(const I& input) + : input(input) + { + initialize(attr, input); + } + }; + + + } // end of namespace mln::labeling::impl + + + + // Facade. + + template <typename I, typename N, typename L> + mln_ch_value(I, L) + regional_minima(const Image& input_, const Neighborhood<N>& nbh_, + L& nlabels) + { + trace::entering("labeling::regional_minima"); + + const I& input = exact(input_); + const N& nbh = exact(nbh_); + mln_precondition(input.is_valid()); + + // FIXME: abort if L is not wide enough to encode the set of + // minima. + + typedef impl::regional_minima_functor F; + F f(exact(input)); + mln_ch_value(I, L) output = canvas::labeling_sorted(input, nbh, nlabels, + f, true); + + trace::exiting("labeling::regional_minima"); + return output; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::labeling + +} // end of namespace mln + + +#endif // ! MLN_LABELING_REGIONAL_MINIMA_HH Index: trunk/milena/sandbox/fabien/regional_maxima.hh =================================================================== Index: trunk/milena/sandbox/fabien/regional_minima.cc =================================================================== --- trunk/milena/sandbox/fabien/regional_minima.cc (revision 0) +++ trunk/milena/sandbox/fabien/regional_minima.cc (revision 3309) @@ -0,0 +1,57 @@ +// Copyright (C) 2007, 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/labeling/regional_minima.cc + * + * \brief Test on mln::labeling::regional_minima. + */ + +#include <mln/core/image/image2d.hh> +#include <mln/value/int_u8.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/io/pgm/load.hh> +#include <mln/pw/all.hh> + +#include "regional_minima.hh" + +#include <tests/data.hh> +#include <mln/debug/println.hh> + + +int main() +{ + using namespace mln; + using value::int_u8; + + trace::quiet = false; + + image2d<int_u8> lena = io::pgm::load<int_u8>(MLN_IMG_DIR "/tiny.pgm"); + + unsigned n; + labeling::regional_minima((pw::cst(255) - pw::value(lena)) | lena.domain(), + c4(), n); +}

1 0

3308: Prepare to fast the morpho wst by flooding.
by Thierry Geraud 06 Feb '09

06 Feb '09

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Thierry Geraud <thierry.geraud(a)lrde.epita.fr> Prepare to fast the morpho wst by flooding. * mln/morpho/watershed: New directory. * mln/morpho/meyer_wst.hh: Copy to... * mln/morpho/watershed/flooding.hh: ...this new file. (watershed): New namespace. (nbasins): Rename as... (n_basins): ...this. (meyer_wst): Remove overload without n_basins. (meyer_wst): Rename this only remaining routine as... (flooding): ...this. (flooding_dispatch): New. * mln/morpho/watershed/all.hh: New. * mln/morpho/all.hh: Update. * tests/morpho/watershed: New directory. * tests/morpho/watershed/flooding.cc: New. * tests/morpho/watershed/Makefile.am: New. * tests/morpho/Makefile.am: Update. mln/morpho/all.hh | 1 mln/morpho/watershed/all.hh | 65 +++++++++ mln/morpho/watershed/flooding.hh | 266 +++++++++++++++++++++++++++++-------- tests/morpho/Makefile.am | 3 tests/morpho/watershed/Makefile.am | 10 + tests/morpho/watershed/flooding.cc | 61 ++++++++ 6 files changed, 351 insertions(+), 55 deletions(-) Index: mln/morpho/watershed/flooding.hh --- mln/morpho/watershed/flooding.hh (revision 0) +++ mln/morpho/watershed/flooding.hh (working copy) @@ -25,19 +25,20 @@ // reasons why the executable file might be covered by the GNU General // Public License. -#ifndef MLN_MORPHO_MEYER_WST_HH -# define MLN_MORPHO_MEYER_WST_HH +#ifndef MLN_MORPHO_WATERSHED_FLOODING_HH +# define MLN_MORPHO_WATERSHED_FLOODING_HH -/** \file mln/morpho/meyer_wst.hh - \brief Meyer's Watershed Transform (WST) algorithm. - - The Watershed Transform algorithm from Meyer using a hierarchical - queue. - - Reference: - Fernand Meyer. Un algorithme optimal de ligne de partage des - eaux. In: Actes du 8�me Congr�s AFCET, Lyon-Villeurbanne, France - (1991), pages 847--859. */ +/// \file mln/morpho/watershed/flooding.hh +/// +/// Meyer's Watershed Transform (WST) algorithm. +/// +/// The Watershed Transform algorithm from Meyer using a hierarchical +/// queue. +/// +/// Reference: +/// Fernand Meyer. Un algorithme optimal de ligne de partage des +/// eaux. In: Actes du 8�me Congr�s AFCET, Lyon-Villeurbanne, France +/// (1991), pages 847--859. # include <mln/trait/ch_value.hh> @@ -56,59 +57,155 @@ namespace morpho { - /* FIXME: Provide also a version of the algorithm taking an image - of minima as input. */ - /* FIXME: See also the interface of the Shortest-Path Watershed + namespace watershed + { + + /* + FIXME: + Provide also a version of the algorithm taking an image + of minima as input. + + FIXME: + See also the interface of the Shortest-Path Watershed Transform, which proposes to lower-complete the image before processing it. Then, add a reference to - mln/morpho/lower_completion.hh. */ + mln/morpho/lower_completion.hh. + */ - /** \brief Meyer's Watershed Transform (WST) algorithm. + /// Meyer's Watershed Transform (WST) algorithm. + /// + /// \param[in] input The input image. + /// \param[in] nbh The connexity of markers. + /// \param[out] n_basins The number of basins. + /// + /// \li \p L is the type of labels, used to number the watershed + /// itself (with the minimal value), and the basins. + /// \li \p I is the exact type of the input image. + /// \li \p N is the exact type of the neighborhood used to express + /// \a input's connexity. - \param[in] input The input image. - \param[in] nbh The connexity of markers. - \param[out] nbasins The number of basins. - - \li \p L is the type of labels, used to number the watershed - itself (with the minimal value), and the basins. - \li \p I is the exact type of the input image. - \li \p N is the exact type of the neighborhood used to express - \a input's connexity. */ template <typename I, typename N, typename L> mln_ch_value(I, L) - meyer_wst(const Image& input, const Neighborhood<N>& nbh, - L& nbasins); + flooding(const Image& input, const Neighborhood<N>& nbh, + L& n_basins); + - /** \brief Meyer's Watershed Transform (WST) algorithm, with no - count of basins. - \param[in] input The input image. - \param[in] nbh The connexity of markers. +# ifndef MLN_INCLUDE_ONLY + - \li \p L is the type of labels, used to number the watershed - itself (with the minimal value), and the basins. - \li \p I is the exact type of the input image. - \li \p N is the exact type of the neighborhood used to express - \a input's connexity. + // Implementations. - Note that the first parameter, \p L, is not automatically - valued from the type of the actual argument during implicit - instantiation: you have to explicitly pass this parameter at - call sites. */ - template <typename L, typename I, typename N> + namespace impl + { + + namespace generic + { + + template <typename I, typename N, typename L> mln_ch_value(I, L) - meyer_wst(const Image& input, const Neighborhood<N>& nbh); + flooding(const Image& input_, const Neighborhood<N>& nbh_, + L& n_basins) + { + trace::entering("morpho::watershed::impl::generic::flooding"); + /* FIXME: Ensure the input image has scalar values. */ + + const I input = exact(input_); + const N nbh = exact(nbh_); + typedef L marker; + const marker unmarked = literal::zero; -# ifndef MLN_INCLUDE_ONLY + typedef mln_value(I) V; + const V max = mln_max(V); + + // Initialize the output with the markers (minima components). + mln_ch_value(I, marker) output = + labeling::regional_minima (input, nbh, n_basins); + + typedef mln_psite(I) psite; + + // Ordered queue. + typedef p_queue_fast<psite> Q; + p_priority<V, Q> queue; + + // In_queue structure to avoid processing sites several times. + mln_ch_value(I, bool) in_queue; + initialize(in_queue, input); + data::fill(in_queue, false); + + // Insert every neighbor P of every marked area in a + // hierarchical queue, with a priority level corresponding to + // the grey level input(P). + mln_piter(I) p(output.domain()); + mln_niter(N) n(nbh, p); + for_all (p) + if (output(p) == unmarked) + for_all(n) + if (output.domain().has(n) && output(n) != unmarked) + { + queue.push(max - input(p), p); + in_queue(p) = true; + break; + } + + /* Until the queue is empty, extract a psite P from the + hierarchical queue, at the highest priority level, that is, + the lowest level. */ + while (! queue.is_empty()) + { + psite p = queue.front(); + queue.pop(); + // Last seen marker adjacent to P. + marker adjacent_marker = unmarked; + // Has P a single adjacent marker? + bool single_adjacent_marker_p = true; + mln_niter(N) n(nbh, p); + for_all(n) + if (output.domain().has(n) && output(n) != unmarked) + { + if (adjacent_marker == unmarked) + { + adjacent_marker = output(n); + single_adjacent_marker_p = true; + } + else + if (adjacent_marker != output(n)) + { + single_adjacent_marker_p = false; + break; + } + } + /* If the neighborhood of P contains only psites with the + same label, then P is marked with this label, and its + neighbors that are not yet marked are put into the + hierarchical queue. */ + if (single_adjacent_marker_p) + { + output(p) = adjacent_marker; + for_all(n) + if (output.domain().has(n) && output(n) == unmarked + && ! in_queue(n)) + { + queue.push(max - input(n), n); + in_queue(n) = true; + } + } + } + trace::exiting("morpho::watershed::impl::generic::flooding"); + return output; + } + + + // Fastest version. template <typename I, typename N, typename L> mln_ch_value(I, L) - meyer_wst(const Image& input_, const Neighborhood<N>& nbh_, - L& nbasins) + flooding_fastest(const Image& input_, const Neighborhood<N>& nbh_, + L& n_basins) { - trace::entering("morpho::meyer_wst"); + trace::entering("morpho::watershed::impl::flooding_fastest"); /* FIXME: Ensure the input image has scalar values. */ const I input = exact(input_); @@ -122,7 +219,7 @@ // Initialize the output with the markers (minima components). mln_ch_value(I, marker) output = - labeling::regional_minima (input, nbh, nbasins); + labeling::regional_minima (input, nbh, n_basins); typedef mln_psite(I) psite; @@ -193,23 +290,84 @@ } } } - trace::exiting("morpho::meyer_wst"); + trace::exiting("morpho::watershed::impl::flooding_fastest"); return output; } - template <typename L, typename I, typename N> + + } // end of namespace mln::morpho::watershed::impl::generic + + } // end of namespace mln::morpho::watershed::impl + + + + // Dispatch. + + namespace internal + { + + template <typename I, typename N, typename L> + inline mln_ch_value(I, L) - meyer_wst(const Image& input, const Neighborhood<N>& nbh) + flooding_dispatch(metal::false_, + const Image& input, const Neighborhood<N>& nbh, L& n_basins) { - L nbasins; - return meyer_wst<L>(input, nbh, nbasins); + return impl::generic::flooding(input, nbh, n_basins); + } + + + template <typename I, typename N, typename L> + inline + mln_ch_value(I, L) + flooding_dispatch(metal::true_, + const Image& input, const Neighborhood<N>& nbh, L& n_basins) + { + return impl::generic::flooding(input, nbh, n_basins); +// return impl::generic::flooding_fastest(input, nbh, n_basins); + } + + template <typename I, typename N, typename L> + inline + mln_ch_value(I, L) + flooding_dispatch(const Image& input, const Neighborhood<N>& nbh, L& n_basins) + { + enum { + test = mlc_equal(mln_trait_image_speed(I), + trait::image::speed::fastest)::value + && + mln_is_simple_neighborhood(N)::value + }; + return flooding_dispatch(metal::bool_<test>(), + input, nbh, n_basins); + } + + } // end of namespace mln::morpho::watershed::impl + + + // Facade. + + template <typename I, typename N, typename L> + inline + mln_ch_value(I, L) + flooding(const Image& input, const Neighborhood<N>& nbh, L& n_basins) + { + trace::entering("morpho::watershed::flooding"); + + // FIXME: internal::flooding_tests(input, nbh, n_basins); + + mln_ch_value(I, L) output = internal::flooding_dispatch(input, nbh, n_basins); + + trace::exiting("morpho::watershed::flooding"); + return output; } # endif // ! MLN_INCLUDE_ONLY + } // end of namespace mln::morpho::watershed + } // end of namespace mln::morpho } // end of namespace mln -#endif // ! MLN_MORPHO_MEYER_WST_HH +#endif // ! MLN_MORPHO_WATERSHED_FLOODING_HH Property changes on: mln/morpho/watershed/flooding.hh ___________________________________________________________________ Added: svn:mergeinfo Index: mln/morpho/watershed/all.hh --- mln/morpho/watershed/all.hh (revision 0) +++ mln/morpho/watershed/all.hh (revision 0) @@ -0,0 +1,65 @@ +// 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_MORPHO_WATERSHED_ALL_HH +# define MLN_MORPHO_WATERSHED_ALL_HH + +/// \file mln/morpho/watershed/all.hh +/// +/// File that includes all morphological watershed routines. + + +namespace mln +{ + + namespace morpho + { + + /// Namespace of morphological watershed routines. + namespace watershed + { + + /// Namespace of morphological watershed routines + /// implementations. + namespace watershed + { + + /// Namespace of morphological watershed routines generic + /// implementations. + namespace generic + {} + + } + } + } +} + + +# include <mln/morpho/watershed/flooding.hh> + + +#endif // ! MLN_MORPHO_WATERSHED_ALL_HH Index: mln/morpho/all.hh --- mln/morpho/all.hh (revision 3307) +++ mln/morpho/all.hh (working copy) @@ -86,6 +86,7 @@ # include <mln/morpho/elementary/all.hh> # include <mln/morpho/tree/all.hh> +# include <mln/morpho/watershed/all.hh> Index: tests/morpho/watershed/flooding.cc --- tests/morpho/watershed/flooding.cc (revision 0) +++ tests/morpho/watershed/flooding.cc (revision 0) @@ -0,0 +1,61 @@ +// 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/morpho/watershed/flooding.cc +/// +/// Test on mln::morpho::watershed::flooding. + +#include <iostream> + +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/window2d.hh> +#include <mln/core/alias/neighb2d.hh> + +#include <mln/value/int_u8.hh> + +#include <mln/morpho/watershed/flooding.hh> + +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +#include "tests/data.hh" + + +int main() +{ + using namespace mln; + using value::int_u8; + + image2d<int_u8> input; + io::pgm::load(input, MLN_IMG_DIR "/squares.pgm"); + + typedef int_u8 L; + L nbasins; + image2d<L> output = morpho::watershed::flooding(input, c4(), nbasins); + + io::pgm::save(output, "out.pgm"); +} Index: tests/morpho/watershed/Makefile.am --- tests/morpho/watershed/Makefile.am (revision 0) +++ tests/morpho/watershed/Makefile.am (revision 0) @@ -0,0 +1,10 @@ +## Process this file through Automake to create Makefile.in -*- Makefile -*- + +include $(top_srcdir)/milena/tests/tests.mk + +check_PROGRAMS = \ + flooding + +flooding_SOURCES = flooding.cc + +TESTS = $(check_PROGRAMS) Index: tests/morpho/Makefile.am --- tests/morpho/Makefile.am (revision 3307) +++ tests/morpho/Makefile.am (working copy) @@ -4,7 +4,8 @@ SUBDIRS = \ elementary \ - tree + tree \ + watershed check_PROGRAMS = \ artificial_line_graph_image_wst \

1 0

3307: Update icp.
by Thierry Geraud 06 Feb '09

06 Feb '09

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Thierry Geraud <thierry.geraud(a)lrde.epita.fr> Update icp. * mln/registration/icp2.hh (remove_too_far_sites): Factor code and change threshold. icp2.hh | 35 ++++++++++++++++++++++++++++------- 1 file changed, 28 insertions(+), 7 deletions(-) Index: mln/registration/icp2.hh --- mln/registration/icp2.hh (revision 3306) +++ mln/registration/icp2.hh (working copy) @@ -67,6 +67,7 @@ #include <mln/core/image/extension_fun.hh> #include <mln/accu/histo.hh> +#include <mln/accu/sum.hh> #include <mln/debug/histo.hh> @@ -278,14 +279,37 @@ template <typename P, typename F> p_array - remove_too_far_sites(image3d<value::rgb8>& out, const p_array& P_bak, + remove_too_far_sites(image3d<value::rgb8>& out, const p_array& P_, const F& closest_point, const std::pair<algebra::quat,mln_vec(P)>& pair, const p_array& X, p_array& removed_set, unsigned r) { - float sd = compute_standard_deviation(P_bak, pair, closest_point); + mln_piter(p_array) p(P_); + accu::histo<value::int_u8> h; +// float sd = compute_standard_deviation(P_, pair, closest_point); + + float sd; + int d_min, d_max; + { + accu::sum<float> s, s2; + for_all(p) + { + vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; + unsigned d_i = closest_point.dmap_X_(Pk_i); + h.take(d_i); + s.take(d_i); + s2.take(d_i * d_i); + } + float mean = s / P_.nsites(); + sd = std::sqrt(s2 / P_.nsites() - mean * mean); + d_min = int(mean - sd); + d_max = int(mean + sd); + } + + std::cout << "Standard deviation = " << sd << std::endl; + std::cout << "d thresholds = " << d_min << ' ' << d_max << std::endl; p_array tmp; unsigned removed = 0; @@ -293,16 +317,13 @@ data::fill(out, literal::white); data::fill((out | X).rw(), literal::black); - accu::histo<value::int_u8> h; - mln_piter(p_array) p(P_bak); for_all(p) { vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; vec3d_f Yk_i = closest_point(Pk_i).to_vec(); - float norm = norm::l2_distance(Yk_i, Pk_i); - h.take(closest_point.dmap_X_(Pk_i)); - if (norm < sd && norm > sd / 2) + int d_i = closest_point.dmap_X_(Pk_i); + if (d_i >= d_min && d_i <= d_max) { tmp.append(p); out(Pk_i) = literal::green;

1 0

3306: Fix registration std dev computation.
by Thierry Geraud 06 Feb '09

06 Feb '09

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Thierry Geraud <thierry.geraud(a)lrde.epita.fr> Fix registration std dev computation. * mln/registration/icp2.hh (compute_standard_deviation): Fix. icp2.hh | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) Index: mln/registration/icp2.hh --- mln/registration/icp2.hh (revision 3305) +++ mln/registration/icp2.hh (working copy) @@ -32,6 +32,7 @@ /// /// Register an image over an another using the ICP algorithm. +# include <cmath> # include <algorithm> # include <mln/core/alias/vec3d.hh> @@ -114,8 +115,11 @@ const p_array& X, const F& closest_point); + + # ifndef MLN_INCLUDE_ONLY + /// Closest point functor based on map distance. template <typename P> class closest_point_with_map @@ -267,7 +271,7 @@ } float d = e_k_accu.to_result(); - sd = math::sqrt((e_k_accu.hook_value_() - 2.5 * math::sqr(d)) / P_.nsites()); + sd = std::sqrt(e_k_accu.hook_value_() / P_.nsites() - d * d); return sd; } @@ -399,7 +403,7 @@ Qk(2,3) = Spx(1,2) + Spx(2,1); Qk(3,2) = Spx(1,2) + Spx(2,1); - return mln::math::jacobi(Qk); + return math::jacobi(Qk); }

1 0

3305: Make distance-geodesic-based routines have fastest versions.
by Thierry Geraud 05 Feb '09

05 Feb '09

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Thierry Geraud <thierry.geraud(a)lrde.epita.fr> 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& 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& 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 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& 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& 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& 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& 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& 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& 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& 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); +}

1 0

r3304: Fix bugs for algebraic filter and cleaning
by Edwin Carlinet 05 Feb '09

05 Feb '09

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox ChangeLog: 2009-02-05 Edwin Carlinet <carlinet(a)lrde.epita.fr> Fix bugs for algebraic filter and cleaning. * edwin/accu.cc: Separate and make clean. * edwin/accu_trait.hh: Test use_p and use_whatever traits. * edwin/algebraic.hh: Algebraic filter. * edwin/card.hh: Card accumulator. * edwin/fredwin.cc: Remove. * edwin/leveling.hh: Leveling filter. --- accu.cc | 161 +++++++--------------------------------------------------- accu_trait.hh | 3 - algebraic.hh | 67 ++++++++++++++++++++++++ card.hh | 33 +++++++++++ leveling.hh | 99 +++++++++++++++++++++++++++++++++++ 5 files changed, 221 insertions(+), 142 deletions(-) Index: trunk/milena/sandbox/edwin/fredwin.cc (deleted) =================================================================== Index: trunk/milena/sandbox/edwin/accu.cc =================================================================== --- trunk/milena/sandbox/edwin/accu.cc (revision 3303) +++ trunk/milena/sandbox/edwin/accu.cc (revision 3304) @@ -1,155 +1,27 @@ -# include <mln/metal/equal.hh> -# include <mln/metal/if.hh> -# include <mln/metal/is_const.hh> - -# include <mln/core/concept/image.hh> -# include <mln/accu/all.hh> -# include <mln/util/pix.hh> -# include <mln/make/pix.hh> -# include "accu_trait.hh" - -using namespace mln; - -namespace mln -{ - - namespace morpho - { - - namespace accu - { - template <typename T> - struct card : public mln::accu::internal::base< unsigned, card<T> > - { - typedef T argument; - void init () { c_ = 0; }; - void take (const card<T>& elt) { ++c_; }; - void take (const T& elt) { ++c_; }; - void take (const mln_value(T)& v) { ++c_; }; - unsigned to_result() const { return c_; }; - operator unsigned () const { return c_; }; - bool is_valid () const { return true; }; - card () { init(); }; - - private: - unsigned c_; - }; - } // accu - } // morpho - - - namespace impl - { - - template <typename I, typename A> - void - algebraic(const Image& input, - Accumulator<A>& acc) - { - const I& ima = exact(input); - A& accu = exact (acc); - - mln_piter(I) p(ima.domain()); - for_all(p) - accu.take(p); - } - - // fast implementation - template <typename I, typename A> - void - leveling_fast (const Image& input, - Accumulator<A>& acc) - { - const I& ima = exact(input); - A& accu = exact (acc); - - mln_pixter(const I) px(ima); - for_all(px) - accu.take (px.val ()); - } - - // generic implementation - template <typename I, typename A> - void - leveling (const Image& input, - Accumulator<A>& acc) - { - const I& ima = exact(input); - A& accu = exact (acc); - - mln_piter(I) p(ima.domain()); - for_all(p) - accu.take (mln::make::pix(ima, p)); - } - - } // impl - - namespace internal - { - template <typename I, typename A> - void - leveling_dispatch(metal::false_, - const Image& input, - Accumulator<A>& acc) - { - impl::leveling(input, acc); - } - - template <typename I, typename A> - inline - void - leveling_dispatch(metal::true_, - const Image& input, - Accumulator<A>& acc) - { - impl::leveling_fast(input, acc); - } - - template <typename I, typename A> - inline - void - leveling_dispatch(const Image& input, - Accumulator<A>& acc) - { - enum { - test = (mlc_equal(mln_trait_image_speed(I), - trait::image::speed::fastest)::value && - mlc_equal(mln_trait_accu_when_pix(A), - trait::accu::when_pix::use_v)::value) - }; - leveling_dispatch(metal::bool_<test>(), input, acc); - } - - } // internal - -} //mln - -// Facade. - -template <typename I, typename A> -void -leveling(const Image& input, - Accumulator<A>& acc) -{ - mln::internal::leveling_dispatch(input, acc); -} - - -# include <mln/accu/all.hh> # include <mln/core/image/image2d.hh> - # include <mln/debug/iota.hh> # include <mln/debug/println.hh> # include <mln/core/var.hh> # include <mln/util/timer.hh> + +# include "accu_trait.hh" +# include "accu_trait.hh" +# include "card.hh" +# include "card.hh" +# include "algebraic.hh" +# include "algebraic.hh" +# include "leveling.hh" +# include "leveling.hh" + + int main() { using namespace mln; typedef image2d<int> I; I ima(1000, 1000); - mln::morpho::accu::card<util::pix > acc; + mln::morpho::accu::card<point2d> acc; float elapsed; mln::util::timer chrono; @@ -160,6 +32,14 @@ acc.init(); chrono.start(); for (int i = 0; i < 50; i++) + algebraic(ima, acc); + elapsed = chrono.stop(); + + std::cout << "(alge) " << elapsed << "s : " << acc.to_result() << std::endl; + /* + acc.init(); + chrono.start(); + for (int i = 0; i < 50; i++) leveling(ima, acc); elapsed = chrono.stop(); @@ -180,5 +60,6 @@ elapsed = chrono.stop(); std::cout << "(fast) " << elapsed << "s : " << acc.to_result() << std::endl; + */ } Index: trunk/milena/sandbox/edwin/leveling.hh =================================================================== --- trunk/milena/sandbox/edwin/leveling.hh (revision 0) +++ trunk/milena/sandbox/edwin/leveling.hh (revision 3304) @@ -0,0 +1,99 @@ +#ifndef LEVELING_HH_ +# define LEVELING_HH_ + + +# include "accu_trait.hh" +# include <mln/metal/equal.hh> +# include <mln/accu/all.hh> +# include <mln/core/concept/image.hh> +# include <mln/util/pix.hh> +# include <mln/make/pix.hh> + +using namespace mln; + +namespace mln +{ + namespace impl + { + // fast implementation + template <typename I, typename A> + void + leveling_fast (const Image& input, + Accumulator<A>& acc) + { + const I& ima = exact(input); + A& accu = exact (acc); + + mln_pixter(const I) px(ima); + for_all(px) + accu.take (px.val ()); + } + + // generic implementation + template <typename I, typename A> + void + leveling (const Image& input, + Accumulator<A>& acc) + { + const I& ima = exact(input); + A& accu = exact (acc); + + mln_piter(I) p(ima.domain()); + for_all(p) + accu.take (mln::make::pix(ima, p)); + } + + } // mln::impl + + namespace internal + { + template <typename I, typename A> + void + leveling_dispatch(metal::false_, + const Image& input, + Accumulator<A>& acc) + { + impl::leveling(input, acc); + } + + template <typename I, typename A> + inline + void + leveling_dispatch(metal::true_, + const Image& input, + Accumulator<A>& acc) + { + impl::leveling_fast(input, acc); + } + + template <typename I, typename A> + inline + void + leveling_dispatch(const Image& input, + Accumulator<A>& acc) + { + enum { + test = (mlc_equal(mln_trait_image_speed(I), + trait::image::speed::fastest)::value && + mlc_equal(mln_trait_accu_when_pix(A), + trait::accu::when_pix::use_v)::value) + }; + leveling_dispatch(metal::bool_<test>(), input, acc); + } + + } // mln::internal + +} //mln + + +// Facade. +template <typename I, typename A> +void +leveling(const Image& input, + Accumulator<A>& acc) +{ + mln::internal::leveling_dispatch(input, acc); +} + + +#endif /* !LEVELING_HH_ */ Index: trunk/milena/sandbox/edwin/algebraic.hh =================================================================== --- trunk/milena/sandbox/edwin/algebraic.hh (revision 0) +++ trunk/milena/sandbox/edwin/algebraic.hh (revision 3304) @@ -0,0 +1,67 @@ +#ifndef ALGEBRAIC_HH_ +# define ALGEBRAIC_HH_ + +# include "accu_trait.hh" +# include <mln/metal/equal.hh> +# include <mln/accu/all.hh> +# include <mln/core/concept/image.hh> + +namespace mln +{ + namespace impl + { + template <typename I, typename A> + void + algebraic(const Image& input, + Accumulator<A>& acc) + { + const I& ima = exact(input); + A& accu = exact (acc); + + mln_piter(I) p(ima.domain()); + for_all(p) + accu.take(p); + } + } // impl + + namespace internal + { + template <typename I, typename A> + inline + void + algebraic_dispatch(metal::true_, + const Image& input, + Accumulator<A>& acc) + { + impl::algebraic(input, acc); + } + + template <typename I, typename A> + inline + void + algebraic_dispatch(const Image& input, + Accumulator<A>& acc) + { + enum { + test = (mlc_equal(mln_trait_accu_when_pix(A), + trait::accu::when_pix::use_p)::value || + mlc_equal(mln_trait_accu_when_pix(A), + trait::accu::when_pix::use_whatever)::value) + }; + algebraic_dispatch(metal::bool_<test>(), input, acc); + } + + } // mln::internal + +} //mln + + +template <typename I, typename A> +void +algebraic(const mln::Image& input, + mln::Accumulator<A>& acc) +{ + mln::internal::algebraic_dispatch(input, acc); +} + +#endif /* !ALGEBRAIC_HH_ */ Index: trunk/milena/sandbox/edwin/card.hh =================================================================== --- trunk/milena/sandbox/edwin/card.hh (revision 0) +++ trunk/milena/sandbox/edwin/card.hh (revision 3304) @@ -0,0 +1,33 @@ +#ifndef CARD_HH_ +# define CARD_HH_ + +# include <mln/accu/all.hh> + +namespace mln +{ + namespace morpho + { + namespace accu + { + template <typename T> + struct card : public mln::accu::internal::base< unsigned, card<T> > + { + typedef T argument; + + card () { init(); }; + void init () { c_ = 0; }; + void take (const card<T>& accu) { c_ += accu.c_; }; + void take (const T& elt) { ++c_; }; + unsigned to_result() const { return c_; }; + bool is_valid () const { return true; }; + + private: + unsigned c_; + }; + } // mln::morpho::accu + } // mln::morpho +} // mln + +#endif /* !CARD_HH_ */ + + Index: trunk/milena/sandbox/edwin/accu_trait.hh =================================================================== --- trunk/milena/sandbox/edwin/accu_trait.hh (revision 3303) +++ trunk/milena/sandbox/edwin/accu_trait.hh (revision 3304) @@ -3,7 +3,6 @@ # include <mln/trait/undef.hh> - # define mln_trait_accu_when_pix(A) typename trait::accu::accu_traits<A>::when_pix @@ -41,7 +40,7 @@ template <typename T> struct accu_traits< mln::morpho::accu::card <T> > { - typedef when_pix::use_v when_pix; + typedef when_pix::use_p when_pix; }; } //accu

1 0

3303: Fix sorted labeling fastest canvas.
by Thierry Geraud 05 Feb '09

05 Feb '09

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Thierry Geraud <thierry.geraud(a)lrde.epita.fr> Fix sorted labeling fastest canvas. * fabien/labeling.hh: Fix. (labeling_video_fastest): Handle extension. (labeling_sorted_fastest): Likewise. labeling.hh | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) Index: fabien/labeling.hh --- fabien/labeling.hh (revision 3302) +++ fabien/labeling.hh (working copy) @@ -39,6 +39,7 @@ # include <mln/data/fill.hh> # include <mln/literal/zero.hh> # include <mln/convert/to_upper_window.hh> +# include <mln/extension/adjust_fill.hh> # include <mln/level/sort_psites.hh> # include <mln/level/sort_offsets.hh> @@ -233,13 +234,15 @@ labeling_video_fastest(const Image& input_, const Neighborhood<N>& nbh_, F& f, L& nlabels) { - trace::entering("canvas::impl::labeling"); + trace::entering("canvas::impl::labeling_video_fastest"); // FIXME: Test?! const I& input = exact(input_); const N& nbh = exact(nbh_); + extension::adjust(input, nbh); + // Local type. typedef mln_psite(I) P; @@ -255,6 +258,7 @@ { initialize(deja_vu, input); mln::data::fill(deja_vu, false); + extension::fill(deja_vu, false); // So that the extension is ignored. initialize(parent, input); @@ -276,7 +280,7 @@ f.init_attr(p); for_all(n) - if (input.has(n) && deja_vu(n)) + if (input.domain().has(n) && deja_vu(n)) { if (f.equiv_(n, p)) { @@ -290,8 +294,8 @@ } else f.do_no_union_(n, p); - (p) = true; } + deja_vu(p) = true; } } @@ -318,7 +322,7 @@ status = true; } - trace::exiting("canvas::impl::labeling"); + trace::exiting("canvas::impl::labeling_video_fastest"); return output; } @@ -332,13 +336,15 @@ labeling_sorted_fastest(const Image& input_, const Neighborhood<N>& nbh_, L& nlabels, const S& s, F& f) { - trace::entering("canvas::impl::labeling"); + trace::entering("canvas::impl::labeling_sorted_fastest"); // FIXME: Test?! const I& input = exact(input_); const N& nbh = exact(nbh_); + extension::adjust(input, nbh); + // Local type. typedef mln_psite(I) P; @@ -354,6 +360,7 @@ { initialize(deja_vu, input); mln::data::fill(deja_vu, false); + extension::fill(deja_vu, false); // So that the extension is ignored. initialize(parent, input); @@ -392,7 +399,7 @@ { // Do-Union. unsigned r = find_root_fastest(parent, n); - if (input.element(r) != input.element(p)) + if (r != p) { parent.element(r) = p; f.merge_attr_(r, p); @@ -432,7 +439,7 @@ status = true; } - trace::exiting("canvas::impl::labeling"); + trace::exiting("canvas::impl::labeling_sorted_fastest"); return output; }

1 0

3302: Fix missing is_valid() in complex and graph windows/neighborhoods.
by Guillaume Lazzara 05 Feb '09

05 Feb '09

* mln/accu/transform_diagonal.hh, * mln/accu/transform_directional.hh, * mln/accu/transform_snake.hh: fix wrong precondition. * mln/core/concept/neighborhood.hh, * mln/core/concept/window.hh: add static check to be sure the method is present. * mln/core/internal/complex_neighborhood_base.hh, * mln/core/internal/complex_window_p_base.hh, * mln/core/internal/graph_window_base.hh, * mln/core/internal/weighted_window_base.hh: add missing is_valid. --- milena/ChangeLog | 17 +++++++++++++++++ milena/mln/accu/transform_diagonal.hh | 2 +- milena/mln/accu/transform_directional.hh | 2 +- milena/mln/accu/transform_snake.hh | 2 +- milena/mln/core/concept/neighborhood.hh | 3 +++ milena/mln/core/concept/window.hh | 2 ++ .../mln/core/internal/complex_neighborhood_base.hh | 11 +++++++++++ milena/mln/core/internal/complex_window_p_base.hh | 12 +++++++++++- milena/mln/core/internal/graph_window_base.hh | 12 ++++++++++++ milena/mln/core/internal/weighted_window_base.hh | 13 ++++++++++++- 10 files changed, 71 insertions(+), 5 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index e6b8418..00edd1e 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,22 @@ 2009-02-05 Guillaume Lazzara <z(a)lrde.epita.fr> + Fix missing is_valid() in complex and graph windows/neighborhoods. + + * mln/accu/transform_diagonal.hh, + * mln/accu/transform_directional.hh, + * mln/accu/transform_snake.hh: fix wrong precondition. + + * mln/core/concept/neighborhood.hh, + * mln/core/concept/window.hh: add static check to be sure the method + is present. + + * mln/core/internal/complex_neighborhood_base.hh, + * mln/core/internal/complex_window_p_base.hh, + * mln/core/internal/graph_window_base.hh, + * mln/core/internal/weighted_window_base.hh: add missing is_valid. + +2009-02-05 Guillaume Lazzara <z(a)lrde.epita.fr> + Add new ICP variants. * mln/registration/icp2.hh: Add two new variants. diff --git a/milena/mln/accu/transform_diagonal.hh b/milena/mln/accu/transform_diagonal.hh index c8788e6..db29cc7 100644 --- a/milena/mln/accu/transform_diagonal.hh +++ b/milena/mln/accu/transform_diagonal.hh @@ -90,8 +90,8 @@ namespace mln const W& win = exact(win_); mln_precondition(input.is_valid()); + mln_precondition(win.is_valid()); mln_precondition(! win.is_empty()); - mln_precondition(! win.is_valid()); (void) input; (void) win; diff --git a/milena/mln/accu/transform_directional.hh b/milena/mln/accu/transform_directional.hh index 6a78fcd..e8a4dd1 100644 --- a/milena/mln/accu/transform_directional.hh +++ b/milena/mln/accu/transform_directional.hh @@ -88,8 +88,8 @@ namespace mln const W& win = exact(win_); mln_precondition(input.is_valid()); + mln_precondition(win.is_valid()); mln_precondition(! win.is_empty()); - mln_precondition(! win.is_valid()); (void) input; (void) win; diff --git a/milena/mln/accu/transform_snake.hh b/milena/mln/accu/transform_snake.hh index 98bcf93..88b6361 100644 --- a/milena/mln/accu/transform_snake.hh +++ b/milena/mln/accu/transform_snake.hh @@ -83,8 +83,8 @@ namespace mln const W& win = exact(win_); mln_precondition(input.is_valid()); + mln_precondition(win.is_valid()); mln_precondition(! win.is_empty()); - mln_precondition(! win.is_valid()); (void) input; (void) win; diff --git a/milena/mln/core/concept/neighborhood.hh b/milena/mln/core/concept/neighborhood.hh index 5547865..5a69b12 100644 --- a/milena/mln/core/concept/neighborhood.hh +++ b/milena/mln/core/concept/neighborhood.hh @@ -117,6 +117,9 @@ namespace mln typedef mln_window(E) window; bool m = (& E::win) == (& E::win); m = 0; + bool (E::*m2)() const = &E::is_valid; + m2 = 0; + # if 0 /* FIXME: Disabled, as win() can either return a const reference or a copy of the window (see documentation above). Hence the diff --git a/milena/mln/core/concept/window.hh b/milena/mln/core/concept/window.hh index 64b05fe..8b1a55c 100644 --- a/milena/mln/core/concept/window.hh +++ b/milena/mln/core/concept/window.hh @@ -185,6 +185,8 @@ namespace mln m3 = 0; unsigned (E::*m4)() const = &E::delta; m4 = 0; + bool (E::*m5)() const = &E::is_valid; + m5 = 0; } static void run(mln::trait::window::definition::unique) diff --git a/milena/mln/core/internal/complex_neighborhood_base.hh b/milena/mln/core/internal/complex_neighborhood_base.hh index 231170b..c1bde4e 100644 --- a/milena/mln/core/internal/complex_neighborhood_base.hh +++ b/milena/mln/core/internal/complex_neighborhood_base.hh @@ -105,6 +105,9 @@ namespace mln /// Create a window corresponding to this neighborhood. const window& win() const; /// \} + + /// Return true by default. + bool is_valid() const; }; @@ -122,6 +125,14 @@ namespace mln return exact(*this); } + template <unsigned D, typename G, typename F, typename B, typename E> + inline + bool + complex_neighborhood_base<D, G, F, B, E>::is_valid() const + { + return true; + } + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::internal diff --git a/milena/mln/core/internal/complex_window_p_base.hh b/milena/mln/core/internal/complex_window_p_base.hh index efc6ea9..21540ba 100644 --- a/milena/mln/core/internal/complex_window_p_base.hh +++ b/milena/mln/core/internal/complex_window_p_base.hh @@ -137,6 +137,9 @@ namespace mln bool is_empty() const; /// Is this window centered? (Always returns \c true). bool is_centered() const; + + /// Return true by default. + bool is_valid() const; /// \} }; @@ -158,8 +161,15 @@ namespace mln return true; } + template <unsigned D, typename G, typename F, typename B, typename E> + bool + complex_window_p_base<D, G, F, B, E>::is_valid() const + { + return true; + } + # endif // ! MLN_INCLUDE_ONLY - + } // end of namespace mln::internal } // end of namespace mln diff --git a/milena/mln/core/internal/graph_window_base.hh b/milena/mln/core/internal/graph_window_base.hh index 8af7d65..3883c8a 100644 --- a/milena/mln/core/internal/graph_window_base.hh +++ b/milena/mln/core/internal/graph_window_base.hh @@ -82,6 +82,9 @@ namespace mln bool is_neighbable_() const; /// \} + /// Return true by default. + bool is_valid() const; + protected: graph_window_base(); }; @@ -146,6 +149,15 @@ namespace mln return true; } + template <typename P, typename E> + inline + bool + graph_window_base<P,E>::is_valid() const + { + return true; + } + + # endif // !MLN_INCLUDE_ONLY diff --git a/milena/mln/core/internal/weighted_window_base.hh b/milena/mln/core/internal/weighted_window_base.hh index f633a53..a0380a5 100644 --- a/milena/mln/core/internal/weighted_window_base.hh +++ b/milena/mln/core/internal/weighted_window_base.hh @@ -88,6 +88,9 @@ namespace mln /// the definition is unique. bool has(const mln_dpsite(W)& dp) const; + /// return true by default. + bool is_valid() const; + protected: weighted_window_base(); }; @@ -155,7 +158,7 @@ namespace mln mln_precondition(i < this->size()); return exact(this)->win().dp(i); } - + template <typename W, typename E> inline bool @@ -168,6 +171,14 @@ namespace mln return exact(this)->win().has(dp); } + template <typename W, typename E> + inline + bool + weighted_window_base<W,E>::is_valid() const + { + return true; + } + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::internal -- 1.5.6.5

1 0

3301: Update IGR's code.
by Guillaume Lazzara 05 Feb '09

05 Feb '09

* lazzara/igr.cc: update code. Add more debug and use a new variant of ICP. --- milena/sandbox/ChangeLog | 7 ++++ milena/sandbox/lazzara/igr.cc | 67 +++++++++++++++++++++++++++++++---------- 2 files changed, 58 insertions(+), 16 deletions(-) diff --git a/milena/sandbox/ChangeLog b/milena/sandbox/ChangeLog index 7c9b399..dc34f80 100644 --- a/milena/sandbox/ChangeLog +++ b/milena/sandbox/ChangeLog @@ -1,3 +1,10 @@ +2009-02-05 Guillaume Lazzara <z(a)lrde.epita.fr> + + Update IGR's code. + + * lazzara/igr.cc: update code. Add more debug and use a new variant of + ICP. + 2009-02-05 Fabien Freling <freling(a)lrde.epita.fr> Enable the fastest version of labeling. diff --git a/milena/sandbox/lazzara/igr.cc b/milena/sandbox/lazzara/igr.cc index 1a0d498..05d8471 100644 --- a/milena/sandbox/lazzara/igr.cc +++ b/milena/sandbox/lazzara/igr.cc @@ -11,8 +11,11 @@ #include <mln/essential/3d.hh> #include <mln/core/image/slice_image.hh> +#include <mln/core/image/tr_image.hh> +#include <mln/core/image/interpolated.hh> #include <mln/io/cloud/load.hh> +#include <mln/util/timer.hh> struct threshold : mln::Function_p2b<threshold> { @@ -219,11 +222,11 @@ keep_largest_component(const mln::Image& ima) using namespace mln; image2d<bool> in_bw_cleaned = fill_holes(ima); - io::pbm::save(in_bw_cleaned, "0x_in_bw_cleaned.pbm"); + io::pbm::save(in_bw_cleaned, "in_bw_cleaned.pbm"); logical::not_inplace(in_bw_cleaned); image2d<bool> in_bw_cleaned_full = fill_holes(in_bw_cleaned); - io::pbm::save(in_bw_cleaned_full, "0x_in_bw_cleaned_full.pbm"); + io::pbm::save(in_bw_cleaned_full, "in_bw_cleaned_full.pbm"); logical::not_inplace(in_bw_cleaned_full); return in_bw_cleaned_full; @@ -244,16 +247,50 @@ get_main_object_shape(const mln::Image& in) // J ima = keep_largest_component(in_bw); J ima = keep_largest_component(in); // io::pbm::save(in_bw, "02_ima.pbm"); - io::pbm::save(in, "02_ima.pbm"); + io::pbm::save(in, "ima.pbm"); std::cout << "Compute gradient" << std::endl; J ima_grad = morpho::gradient(ima, win_c4p()); - io::pbm::save(ima_grad, "03_ima_grad.pbm"); + io::pbm::save(ima_grad, "ima_grad.pbm"); return ima_grad; } +namespace mln +{ + + namespace debug + { + template <typename I, typename T> + void + compare_registration(Image& P_, Image& X_, const T& transf) + { + I& P = exact(P_); + I& X = exact(X_); + + mln_pset(I) box = geom::bbox(X); + box.enlarge(40); + + typedef mln_ch_value(I,value::rgb8) result_t; + result_t result(box); + extension_fun<result_t,pw::cst_<mln_value(result_t)> > ext_result(result, pw::cst(value::rgb8(0,0,0))); + extension_fun<I,pw::cst_<mln_value(I)> > ext_X(X, pw::cst(false)); + + data::fill(ext_result, literal::black); + data::fill((ext_result | (pw::value(ext_X) == true)).rw(), literal::white); + + mln_VAR(ig, (P | pw::value(P) == true)); + mln_piter(ig_t) p(ig.domain()); + for_all(p) + ext_result(transf(p.to_vec())) = literal::green; + + io::ppm::save(slice(ext_result,0), "registered-1.ppm"); + } + } +} + + int main(int, char* argv[]) { using namespace mln; @@ -288,21 +325,19 @@ int main(int, char* argv[]) K in_3d = convert::to<K>(in_3d_); K ref_3d = convert::to<K>(ref_3d_); - registration::basic_closest_point<point3d> closest_point(ref_3d_); + registration::closest_point_with_map<point3d> closest_point(ref_3d_); +// registration::closest_point_basic<point3d> closest_point(ref_3d_); + + util::timer t; + t.start(); + typedef rotation<3u,float> rot_t; typedef translation<3u,float> trans_t; - composed<rot_t,trans_t> qk = registration::icp(in_3d_, ref_3d_, closest_point); + composed<trans_t,rot_t> qk = registration::icp_clean2(in_3d_, ref_3d_, closest_point); - std::cout << "* Build result image" << std::endl; - box3d box = geom::bbox(ref_3d); - box.enlarge(40); - K res(box); - data::fill(res, false); + std::cout << "igr.cc - Registration - " << t << std::endl; - mln_VAR(ig, (in_3d | pw::value(in_3d) == true)); - mln_piter_(ig_t) p(ig.domain()); - for_all(p) - res(qk(p.to_vec())) = true; + std::cout << "* Build result image" << std::endl; + debug::compare_registration(in_3d, ref_3d, qk); - io::pbm::save(slice(res,0), "04_registered.pbm"); } -- 1.5.6.5

1 0

3300: Add new ICP variants.
by Guillaume Lazzara 05 Feb '09

05 Feb '09

* mln/registration/icp2.hh: Add two new variants. --- milena/ChangeLog | 6 + milena/mln/registration/icp2.hh | 565 +++++++++++++++++++++++++++++++-------- 2 files changed, 457 insertions(+), 114 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index 8960f74..e6b8418 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,4 +1,10 @@ 2009-02-05 Guillaume Lazzara <z(a)lrde.epita.fr> + + Add new ICP variants. + + * mln/registration/icp2.hh: Add two new variants. + +2009-02-05 Guillaume Lazzara <z(a)lrde.epita.fr> Small changes in transformation related classes. diff --git a/milena/mln/registration/icp2.hh b/milena/mln/registration/icp2.hh index 2993caa..070b4d9 100644 --- a/milena/mln/registration/icp2.hh +++ b/milena/mln/registration/icp2.hh @@ -32,6 +32,8 @@ /// /// Register an image over an another using the ICP algorithm. +# include <algorithm> + # include <mln/core/alias/vec3d.hh> # include <mln/math/jacobi.hh> # include <mln/fun/x2x/all.hh> @@ -49,6 +51,24 @@ # include <mln/canvas/distance_geodesic.hh> # include <mln/pw/all.hh> +# include <mln/io/ppm/save.hh> +# include <mln/io/pbm/save.hh> +# include <mln/debug/colorize.hh> + +# include <mln/literal/black.hh> +# include <mln/literal/white.hh> +# include <mln/literal/colors.hh> + +# include <mln/core/image/slice_image.hh> +# include <mln/util/timer.hh> + +#include <mln/core/image/tr_image.hh> +#include <mln/core/image/extension_fun.hh> + +#include <mln/accu/histo.hh> +#include <mln/debug/histo.hh> + + namespace mln { @@ -61,8 +81,8 @@ namespace mln /*! Register point in \p c using a function of closest points * \p closest_point. * - * \param[in] data_P The cloud of points. - * \param[in] model_X the reference surface. + * \param[in] P_ The cloud of points. + * \param[in] X the reference surface. * \param[in] closest_point The function of closest points. * \param[out] qk The rigid transformation obtained. * @@ -73,54 +93,84 @@ namespace mln * vector as arguments. Otherwise the resulting transformation may be * wrong due to the truncation of the vector coordinate values. * - * \pre \p data_p and \p model_X must not be empty. + * \pre \p P_ and \p X must not be empty. * * Reference article: "A Method for Registration of 3-D Shapes", Paul J. * Besl and Neil D. McKay, IEEE, 2, February 1992. * */ template <typename P, typename F> - composed< rotation<P::dim,float>,translation<P::dim,float> > - icp(const p_array& data_P, - const p_array& model_X, +// composed< translation<P::dim,float>,rotation<P::dim,float> > + std::pair<algebra::quat,mln_vec(P)> + icp(const p_array& P_, + const p_array& X, const F& closest_point, const algebra::quat& initial_rot, const mln_vec(P)& initial_translation); template <typename P, typename F> - composed< rotation<P::dim,float>,translation<P::dim,float> > - icp(const p_array& data_P, - const p_array& model_X, + composed< translation<P::dim,float>,rotation<P::dim,float> > + icp(const p_array& P_, + const p_array& X, const F& closest_point); # ifndef MLN_INCLUDE_ONLY /// Closest point functor based on map distance. template <typename P> - class cp_with_map_t + class closest_point_with_map { typedef mln_image_from_grid(mln_grid(P), P) I; typedef mln_ch_value(I, P) cp_ima_t; + typedef mln_ch_value(I,value::int_u8) dmap_t; public: - cp_with_map_t(const p_array& model_X) + closest_point_with_map(const p_array& X) { - box3d box = geom::bbox(model_X); - box.enlarge(50); + box3d box = geom::bbox(X); + box.enlarge(1, box.nrows() / 2); + box.enlarge(2, box.ncols() / 2); + std::cout << "Map image defined on " << box << std::endl; typedef mln_ch_value(I, bool) model_t; model_t model(box); data::fill(model, false); - data::fill((model | model_X).rw(), true); + data::fill((model | X).rw(), true); - transform::internal::closest_point_functor f; - typedef mln_ch_value(I,value::int_u16) dmap_t; - dmap_t dmap_X = canvas::distance_geodesic(model, c6(), - mln_max(value::int_u16), - f); + transform::internal::closest_point_functor<model_t> f; + util::timer t; + t.start(); + dmap_X_ = canvas::distance_geodesic(model, c6(), + mln_max(value::int_u8), + f); + std::cout << "canvas::distance_geodesic - " << t << "s" << std::endl; cp_ima_ = f.cp_ima; + +#ifndef NDEBUG + mln_ch_value(I, bool) debug2(box); + data::fill(debug2, false); + mln_ch_value(I, value::rgb8) debug(box); + mln_piter(p_array) p(X); + for_all(p) + { + debug(p) = debug::internal::random_color(value::rgb8()); + debug2(p) = true; + } + io::pbm::save(slice(debug2,0), "debug2-a.ppm"); + + mln_piter(I) pi(cp_ima_.domain()); + for_all(pi) + { + debug(pi) = debug(cp_ima_(pi)); + debug2(pi) = debug2(cp_ima_(pi)); + } + + io::pbm::save(slice(debug2,0), "debug2-b.ppm"); + io::ppm::save(slice(debug,0), "debug.ppm"); + std::cout << "map saved" << std::endl; +#endif } mln_site(I) operator()(const mln_site(I)& p) const @@ -128,6 +178,7 @@ namespace mln return cp_ima_(p); } + dmap_t dmap_X_; private: cp_ima_t cp_ima_; }; @@ -135,40 +186,27 @@ namespace mln /// Closest point functor based on map distance. template <typename P> - class basic_closest_point + class closest_point_basic { typedef mln_image_from_grid(mln_grid(P), P) I; typedef p_array X_t; public: - basic_closest_point(const p_array& X) + closest_point_basic(const p_array& X) : X_(X) { } - float - l2_distance(const vec3d_f& vec1, const vec3d_f& vec2) const - { - typedef float D; - D d = 0; - for (unsigned i = 0; i < 3; ++i) - { - D sqr_v1_v2 = static_cast<D>(mln::math::sqr(vec1[i] - vec2[i])); - d = static_cast<D>(d + sqr_v1_v2); - } - return d; - } - mln_site(I) operator()(const vec3d_f& v) const { vec3d_f best_x = convert::to<vec3d_f>(X_[0].to_vec()); - float best_d = l2_distance(v, best_x); + float best_d = norm::l2_distance(v, best_x); mln_piter(X_t) X_i(X_); for_all(X_i) { - float d = l2_distance(v, convert::to<vec3d_f>(X_i)); + float d = norm::l2_distance(v, convert::to<vec3d_f>(X_i)); if (d < best_d) { best_d = d; @@ -183,61 +221,146 @@ namespace mln }; + template <typename P> + void + draw_last_run(const box3d& box, const p_array& kept, + const p_array& removed, const p_array& X, + const algebra::quat& qR, const vec3d_f qT) + { + typedef image3d<value::rgb8> result_t; + result_t result(box); + typedef extension_fun<result_t,pw::cst_<mln_value(result_t)> > ext_result_t; + ext_result_t ext_result(result, pw::cst(value::rgb8(0,0,0))); + + data::fill(ext_result, literal::black); + data::fill((ext_result | X).rw(), literal::white); + + mln_piter(p_array) p(kept); + for_all(p) + ext_result(qR.rotate(p.to_vec()) + qT) = literal::green; + + mln_piter(p_array) p2(removed); + for_all(p2) + ext_result(qR.rotate(p2.to_vec()) + qT) = literal::red; + + io::ppm::save(slice(ext_result,0), "registered-2.ppm"); + } + + + template <typename P, typename F> - inline - float - compute_d_k(const p_array& data_P, - const F& closest_point, - const algebra::quat& qR, - const algebra::quat& qR_old, - const vec3d_f& qT, - const vec3d_f& qT_old) + float compute_standard_deviation(const p_array& P_, + const std::pair<algebra::quat,mln_vec(P)>& pair, + const F& closest_point) { - accu::rms<vec3d_f, float> accu; - mln_piter(p_array) p(data_P); + accu::rms<vec3d_f,float> e_k_accu; + + // Standard Deviation + float sd; + mln_piter(p_array) p(P_); for_all(p) { - // yk_i - pk+1_i - vec3d_f Pk_i = qR_old.rotate(p.to_vec()) + qT_old; - vec3d_f Pk_1_i = qR.rotate(p.to_vec()) + qT; - accu.take(closest_point(Pk_i).to_vec() - Pk_1_i); + vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; + vec3d_f Yk_i = closest_point(Pk_i).to_vec(); + // yk_i - pk_i + e_k_accu.take(Yk_i - Pk_i); } - return accu.to_result(); + + float d = e_k_accu.to_result(); + sd = math::sqrt((e_k_accu.hook_value_() - 2.5 * math::sqr(d)) / P_.nsites()); + return sd; } + + template <typename P, typename F> + p_array + remove_too_far_sites(image3d<value::rgb8>& out, const p_array& P_bak, + const F& closest_point, + const std::pair<algebra::quat,mln_vec(P)>& pair, + const p_array& X, p_array& removed_set, + unsigned r) + { + float sd = compute_standard_deviation(P_bak, pair, closest_point); + std::cout << "Standard deviation = " << sd << std::endl; + + p_array tmp; + unsigned removed = 0; + + data::fill(out, literal::white); + data::fill((out | X).rw(), literal::black); + + accu::histo<value::int_u8> h; + mln_piter(p_array) p(P_bak); + for_all(p) + { + vec3d_f Pk_i = pair.first.rotate(p.to_vec()) + pair.second; + vec3d_f Yk_i = closest_point(Pk_i).to_vec(); + + float norm = norm::l2_distance(Yk_i, Pk_i); + h.take(closest_point.dmap_X_(Pk_i)); + if (norm < sd && norm > sd / 2) + { + tmp.append(p); + out(Pk_i) = literal::green; + } + else + { + ++removed; + removed_set.append(p); + out(Pk_i) = literal::red; + } + } + + std::ostringstream ss1; + ss1 << "histo_" << r << ".dat"; + debug::histo_plot(h, ss1.str()); + std::cout << h << std::endl; + + std::ostringstream ss2; + ss2 << "out_0" << r << ".ppm"; + io::ppm::save(mln::slice(out,0), ss2.str()); + + std::cout << "Points removed: " << removed << std::endl; + + return tmp; + } + + template <typename P, typename F> inline float - compute_e_k(const p_array& data_P, + compute_d_k(const p_array& P_, const F& closest_point, const algebra::quat& qR, - const vec3d_f& qT) + const algebra::quat& qR_old, + const vec3d_f& qT, + const vec3d_f& qT_old) { accu::rms<vec3d_f, float> accu; - mln_piter(p_array) p(data_P); + mln_piter(p_array) p(P_); for_all(p) { - // yk_i - pk_i - vec3d_f Pk_i = qR.rotate(p.to_vec()) + qT; - accu.take(closest_point(Pk_i).to_vec() - Pk_i); + // yk_i - pk+1_i + vec3d_f Pk_i = qR_old.rotate(p.to_vec()) + qT_old; + vec3d_f Pk_1_i = qR.rotate(p.to_vec()) + qT; + accu.take(closest_point(Pk_i).to_vec() - Pk_1_i); } return accu.to_result(); } - /// FIXME: work only for 3d images. template <typename P, typename F> algebra::quat - get_rot(const p_array& data_P, + get_rot(const p_array& P_, const vec3d_f& mu_P, const vec3d_f& mu_Yk, const F& closest_point, const algebra::quat& qR, const vec3d_f& qT) { - /// Mk: cross-covariance matrix. - algebra::mat<3u,3u,float> Mk; - mln_piter(p_array) p(data_P); + /// Spx: cross-covariance matrix. + algebra::mat<3u,3u,float> Spx; + mln_piter(p_array) p(P_); // FIXME: could we use an accu? for_all(p) @@ -245,17 +368,17 @@ namespace mln vec3d_f P_i = p.to_vec(); vec3d_f Pk_i = qR.rotate(p.to_vec()) + qT; vec3d_f Yk_i = closest_point(Pk_i).to_vec(); - Mk += make::mat(P_i - mu_P) * trans(make::mat(Yk_i - mu_Yk)); + Spx += make::mat(P_i - mu_P) * trans(make::mat(Yk_i - mu_Yk)); } - Mk /= data_P.nsites(); + Spx /= P_.nsites(); vec3d_f A; - A[0] = Mk(1,2) - Mk(2,1); - A[1] = Mk(2,0) - Mk(0,2); - A[2] = Mk(0,1) - Mk(1,0); + A[0] = Spx(1,2) - Spx(2,1); + A[1] = Spx(2,0) - Spx(0,2); + A[2] = Spx(0,1) - Spx(1,0); algebra::mat<4u,4u,float> Qk; - float t = tr(Mk); + float t = tr(Spx); Qk(0,0) = t; for (int i = 1; i < 4; ++i) @@ -264,17 +387,17 @@ namespace mln Qk(0,i) = A[i - 1]; for (int j = 1; j < 4; ++j) if (i == j) - Qk(i,j) = 2 * Mk(i - 1,i - 1) - t; + Qk(i,j) = 2 * Spx(i - 1,i - 1) - t; } - Qk(1,2) = Mk(0,1) + Mk(1,0); - Qk(2,1) = Mk(0,1) + Mk(1,0); + Qk(1,2) = Spx(0,1) + Spx(1,0); + Qk(2,1) = Spx(0,1) + Spx(1,0); - Qk(1,3) = Mk(0,2) + Mk(2,0); - Qk(3,1) = Mk(0,2) + Mk(2,0); + Qk(1,3) = Spx(0,2) + Spx(2,0); + Qk(3,1) = Spx(0,2) + Spx(2,0); - Qk(2,3) = Mk(1,2) + Mk(2,1); - Qk(3,2) = Mk(1,2) + Mk(2,1); + Qk(2,3) = Spx(1,2) + Spx(2,1); + Qk(3,2) = Spx(1,2) + Spx(2,1); return mln::math::jacobi(Qk); } @@ -284,63 +407,83 @@ namespace mln template <typename P, typename F> inline vec3d_f - get_mu_yk(const p_array& data_P, + get_mu_yk(const p_array& P_, const F& closest_point, const algebra::quat& qR, - const vec3d_f& qT) + const vec3d_f& qT, + float& e_k) { - accu::center<P,vec3d_f> accu; + accu::rms<vec3d_f,float> e_k_accu; + accu::center<P,vec3d_f> mu_yk; - mln_piter(p_array) p(data_P); + mln_piter(p_array) p(P_); for_all(p) { + // yk_i - pk_i vec3d_f Pk_i = qR.rotate(p.to_vec()) + qT; - accu.take(closest_point(Pk_i).to_vec()); + vec3d_f Yk_i = closest_point(Pk_i).to_vec(); + mu_yk.take(Yk_i); + e_k_accu.take(Yk_i - Pk_i); } - return accu.to_result(); + + e_k = e_k_accu.to_result(); + return mu_yk.to_result(); } + + /// Base version of the ICP algorithm. It is called in other variants. template <typename P, typename F> inline - composed< rotation<P::dim,float>,translation<P::dim,float> > - icp(const p_array& data_P, - const p_array& model_X, +// composed< translation<P::dim,float>,rotation<P::dim,float> > + std::pair<algebra::quat,mln_vec(P)> + icp(const p_array& P_, + const p_array& X, const F& closest_point, const algebra::quat& initial_rot, const mln_vec(P)& initial_translation) { trace::entering("registration::icp"); + (void) X; mln_precondition(P::dim == 3); - mln_precondition(!data_P.is_empty()); - mln_precondition(!model_X.is_empty()); + mln_precondition(!P_.is_empty()); + mln_precondition(!X.is_empty()); mln_precondition(!initial_rot.is_null()); typedef p_array cloud_t; - vec3d_f mu_P = set::compute(accu::center<P,vec3d_f>(), data_P); + vec3d_f mu_P = set::compute(accu::center<P,vec3d_f>(), P_); vec3d_f qT_old, qT = initial_translation; algebra::quat qR_old, qR = initial_rot; - float ek, ek_old = mln_max(float); - float dk, dk_old = mln_max(float); + float e_k, e_k_old = mln_max(float); + float d_k, d_k_old = mln_max(float); unsigned k = 0; + +# ifndef NDEBUG + box3d box = geom::bbox(X); + //FIXME: too large? + box.enlarge(1, box.nrows() / 2); + box.enlarge(2, box.ncols() / 2); + image3d<value::rgb8> debug(box); + data::fill(debug, literal::black); + data::fill((debug | X).rw(), literal::white); +# endif + do { qT_old = qT; qR_old = qR; - // Compute error ek = d(Pk,Yk) - ek = compute_e_k(data_P, closest_point, qR, qT); - /// Compute transformation /// // mu_Yk - Pk's mass center. - vec3d_f mu_Yk = get_mu_yk(data_P, closest_point, qR_old, qT_old); + // + Compute error ek = d(Pk,Yk) + vec3d_f mu_Yk = get_mu_yk(P_, closest_point, qR_old, qT_old, e_k); // quaternion qR - rotation - qR = get_rot(data_P, mu_P, mu_Yk, closest_point, qR_old, qT_old); + qR = get_rot(P_, mu_P, mu_Yk, closest_point, qR_old, qT_old); vec3d_f tmp = qR.v(); // vector qT - translation @@ -349,42 +492,236 @@ namespace mln /// End of "compute transformation" // Distance dk = d(Pk+1, Yk) - dk = compute_d_k(data_P, closest_point, qR, qR_old, qT, qT_old); - - // Check property according the related paper. - mln_assertion(0 <= dk); - mln_assertion(dk <= ek); - mln_assertion(ek <= dk_old); - mln_assertion(dk_old <= ek_old); + d_k = compute_d_k(P_, closest_point, qR, qR_old, qT, qT_old); + + +#ifndef NDEBUG + image3d<value::rgb8> tmp_ = duplicate(debug); + mln_piter(p_array) p_dbg(P_); + for_all(p_dbg) + tmp_(qR_old.rotate(p_dbg.to_vec()) + qT_old) = literal::green; + std::ostringstream ss; + ss << "reg/tmp_0"; + if (k < 10) + ss << "0"; + ss << k << ".ppm"; + io::ppm::save(mln::slice(tmp_,0), ss.str()); +#endif + + std::cout << "e_" << k << "=" << e_k << std::endl; + std::cout << "d_" << k << "=" << d_k << std::endl; + + // Check distance and error according to the related paper. + mln_assertion(0 <= d_k); + mln_assertion(d_k <= e_k); + + // Disabled because of the following 'if' +// mln_assertion(e_k <= d_k_old); +// mln_assertion(d_k_old <= e_k_old); + + if (d_k > d_k_old) + { + qR = qR_old; + qT = qT_old; + break; + } // Backing up results. - dk_old = dk; - ek_old = ek; + d_k_old = d_k; + e_k_old = e_k; ++k; } while ((k < 3) - || norm::l2((qT - qT_old)) + norm::l2((qR - qR_old).to_vec()) > 1e-5); + || norm::l2((qT - qT_old)) + norm::l2((qR - qR_old).to_vec()) > 1e-3); + + trace::exiting("registration::icp"); + return std::make_pair(qR, qT); + } + + + + /// Single call to ICP with all sites. + template <typename P, typename F> + inline + composed< translation<P::dim,float>,rotation<P::dim,float> > + icp(const p_array& P_, + const p_array& X, + const F& closest_point) + { + util::timer t; + t.start(); + + std::pair<algebra::quat,mln_vec(P)> pair = icp(P_, X, closest_point, + algebra::quat(1,0,0,0), + literal::zero); + std::cout << "icp = " << t << std::endl; typedef rotation<3u,float> rot_t; - rot_t tqR(qR); + rot_t tqR(pair.first); typedef translation<3u,float> trans_t; - trans_t tqT(qT); - composed<rot_t,trans_t> result(tqR, tqT); + trans_t tqT(pair.second); + composed<trans_t, rot_t> result(tqT, tqR); - trace::exiting("registration::icp"); return result; } + + /// Shuffle the sites in P_. + /// Use one third of P_'s sites for each run. + /// For each run, it removes sites which are too close or too far. template <typename P, typename F> inline - composed< rotation<P::dim,float>,translation<P::dim,float> > - icp(const p_array& data_P, - const p_array& model_X, - const F& closest_point) + composed< translation<P::dim,float>,rotation<P::dim,float> > + icp_clean(const p_array& P_, + const p_array& X, + const F& closest_point) { - return icp(data_P, model_X, closest_point, - algebra::quat(1,0,0,0), literal::zero); + util::timer t; + t.start(); + + // P_bak is shuffled. + p_array P_bak = P_; + std::vector<mln_element(p_array)>& v = P_bak.hook_std_vector_(); + std::random_shuffle(v.begin(), v.end()); + + // P_sub = 1/3 * P_bak; + p_array P_sub = P_bak; + P_sub.hook_std_vector_().resize(P_bak.nsites() / 3); + + unsigned r = 0; + std::pair<algebra::quat,mln_vec(P)> pair; + pair.first = algebra::quat(1,0,0,0); + pair.second = literal::zero; + box3d box = geom::bbox(X); + box.enlarge(40); + image3d<value::rgb8> out(box); + p_array removed_set; + do + { + /// Compute transformation. + pair = icp(P_sub, X, closest_point, + pair.first, + pair.second); + + pair = icp(P_sub, X, closest_point, + pair.first, + pair.second); + + P_sub = remove_too_far_sites(out, P_sub, + closest_point, pair, X, removed_set, r); + + + ++r; + + //Add more data + if (r < 3) + for (unsigned i = (P_bak.nsites() / 3) * r; + i < (P_bak.nsites() / 3) * (r + 1); ++i) + { + P_sub.append(P_bak[i]); + } + + } while (r < 4); + std::cout << "icp = " << t << std::endl; + + typedef rotation<3u,float> rot_t; + rot_t tqR(pair.first); + typedef translation<3u,float> trans_t; + trans_t tqT(pair.second); + composed<trans_t,rot_t> result(tqT, tqR); + + return result; + } + + + /// Shuffle sites in P_. + /// Do the first run with all sites. + /// For each run, remove sites which are too far or too close. + template <typename P, typename F> + inline + composed< translation<P::dim,float>,rotation<P::dim,float> > + icp_clean2(const p_array& P_, + const p_array& X, + const F& closest_point) + { + util::timer t; + t.start(); + + // P_bak is shuffled. + p_array P_bak = P_; + + unsigned r = 0; + std::pair<algebra::quat,mln_vec(P)> pair; + pair.first = algebra::quat(1,0,0,0); + pair.second = literal::zero; + box3d box = geom::bbox(X); + box.enlarge(40); + image3d<value::rgb8> out(box); + p_array removed_set; + do + { + pair = icp(P_bak, X, closest_point, + pair.first, + pair.second); + + P_bak = remove_too_far_sites(out, P_bak, + closest_point, pair, X, removed_set, r); + + ++r; + + } while (r < 4); + std::cout << "icp = " << t << std::endl; + + draw_last_run(box, P_bak, removed_set, X, pair.first, pair.second); + + typedef rotation<3u,float> rot_t; + rot_t tqR(pair.first); + typedef translation<3u,float> trans_t; + trans_t tqT(pair.second); + composed<trans_t,rot_t> result(tqT, tqR); + + return result; + } + + + + /// Run icp once with 1/10 of the sites and run it once again with the + /// resulting tranformation and all the sites. + template <typename P, typename F> + inline + composed< translation<P::dim,float>,rotation<P::dim,float> > + icp_fast(const p_array& P_, + const p_array& X, + const F& closest_point) + { + typedef std::pair<algebra::quat,mln_vec(P)> pair_t; + + p_array P_sub = P_; + std::vector<mln_element(p_array)>& v = P_sub.hook_std_vector_(); + std::random_shuffle(v.begin(), v.end()); + v.resize(P_.nsites() / 10); + + util::timer t; + t.start(); + pair_t tmp = icp(P_sub, X, closest_point, + algebra::quat(1,0,0,0), literal::zero); + + std::cout << "icp_1 - " << t << "s" << std::endl; + t.restart(); + + pair_t tmp2 = icp(P_, X, closest_point, + tmp.first, tmp.second); + + std::cout << "icp_2 - " << t << "s" << std::endl; + + typedef rotation<3u,float> rot_t; + rot_t tqR(tmp2.first); + typedef translation<3u,float> trans_t; + trans_t tqT(tmp2.second); + composed<rot_t,trans_t> result(tqR, tqT); + + return result; } # endif // ! MLN_INCLUDE_ONLY -- 1.5.6.5

1 0