[Olena-patches] 4680: Add the reconstruction on set by union-find canvas.

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Thierry Geraud &lt;thierry.geraud(a)lrde.epita.fr&gt; Add the reconstruction on set by union-find canvas. * theo/mln/morpho/reconstruction/by_dilation.hh: Layout. * theo/mln/morpho/canvas: New directory. * theo/mln/morpho/canvas/reconstruction_on_set_tiny.cc: New. * theo/mln/morpho/canvas/reconstruction_on_set.hh: New. * theo/mln/morpho/canvas/internal: New directory. * theo/mln/morpho/canvas/internal/find_root.hh: New. * theo/mln/morpho/canvas/reconstruction_on_set.cc: New. canvas/internal/find_root.hh | 89 ++++++++ canvas/reconstruction_on_set.cc | 376 +++++++++++++++++++++++++++++++++++ canvas/reconstruction_on_set.hh | 356 +++++++++++++++++++++++++++++++++ canvas/reconstruction_on_set_tiny.cc | 190 +++++++++++++++++ reconstruction/by_dilation.hh | 2 5 files changed, 1011 insertions(+), 2 deletions(-) Index: theo/mln/morpho/reconstruction/by_dilation.hh --- theo/mln/morpho/reconstruction/by_dilation.hh (revision 4679) +++ theo/mln/morpho/reconstruction/by_dilation.hh (working copy) @@ -656,7 +656,6 @@ } } - // Sequence---Forth. { data::fill(deja_vu, false); @@ -687,7 +686,6 @@ } } - // Propagation. { P p; Index: theo/mln/morpho/canvas/reconstruction_on_set_tiny.cc --- theo/mln/morpho/canvas/reconstruction_on_set_tiny.cc (revision 0) +++ theo/mln/morpho/canvas/reconstruction_on_set_tiny.cc (revision 0) @@ -0,0 +1,190 @@ +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/debug/println.hh> + +#include "reconstruction_on_set.hh" + + + +template <typename F, typename G> +struct rd_functor +{ + const F& f; + const G& g; + + rd_functor(const F& f, const G& g) + : f(f), g(g) + { + } + + mln_concrete(F) output; + typedef mln_psite(F) P; + + void set_default() + { + mln::data::fill(output, false); + } + + bool domain(const P& p) const + { + return g(p) == true; + } + + void init(const P& p) + { + output(p) = f(p); + } + + void merge(const P& p, const P& r) + { + output(p) = output(p) || output(r); + } + + void border(const P&, const P&) + { + } + + // Fastest versions. + void set_extension_() + { + mln::border::fill(g, false); + } + bool domain_(unsigned p) const + { + return g.element(p) == true; + } + void init_(unsigned p) + { + output.element(p) = f.element(p); + } + void merge_(unsigned p, unsigned r) + { + output.element(p) = output.element(p) || output.element(r); + } + bool inspect_border_() const + { + return false; + } + void border_(unsigned, unsigned) + { + } +}; + + + + + +template <typename F, typename G> +struct rd_alt_functor // /!\ alternative version +{ + const F& f; + const G& g; + + rd_alt_functor(const F& f, const G& g) + : f(f), g(g) + { + } + + mln_concrete(F) output; + typedef mln_psite(F) P; + + void set_default() + { + mln::data::fill(output, f); + } + + bool domain(const P& p) const + { + return g(p) && ! f(p); + } + + void init(const P& p) + { + output(p) = false; + } + + void merge(const P& p, const P& r) + { + output(p) = output(p) || output(r); + } + + void border(const P& p, const P& n) + { + if (f(n)) + output(p) = true; + } + + // Fastest versions. + void set_extension_() + { + mln::border::fill(f, false); + mln::border::fill(g, false); + } + bool domain_(unsigned p) const + { + return g.element(p) && ! f.element(p); + } + void init_(unsigned p) + { + output.element(p) = false; + } + void merge_(unsigned p, unsigned r) + { + output.element(p) = output.element(p) || output.element(r); + } + bool inspect_border_() const + { + return true; + } + void border_(unsigned p, unsigned n) + { + if (f.element(n)) + { + mln_invariant(f.domain().has(f.point_at_index(n))); + output.element(p) = true; + } + } +}; + + + + + +int main() +{ + using namespace mln; + + typedef image2d<bool> I; + + I f, g, o; + neighb2d nbh = c4(); + + bool gvals[9][5] = { { 0, 1, 1, 1, 0 }, + { 0, 0, 0, 1, 0 }, + { 0, 1, 1, 1, 0 }, + { 0, 1, 0, 0, 0 }, + { 0, 1, 1, 1, 0 }, + { 0, 0, 0, 1, 0 }, + { 0, 1, 1, 1, 0 }, + { 0, 1, 0, 0, 0 }, + { 0, 1, 1, 1, 0 } }; + g = make::image(gvals); + debug::println("g", g); + + bool fvals[9][5] = { { 0, 0, 0, 0, 0 }, + { 0, 0, 0, 0, 0 }, + { 0, 0, 0, 0, 0 }, + { 0, 0, 0, 0, 0 }, + { 0, 0, 1, 0, 0 }, + { 0, 0, 0, 0, 0 }, + { 0, 0, 0, 0, 0 }, + { 0, 0, 0, 0, 0 }, + { 0, 0, 0, 0, 0 } }; + f = make::image(fvals); + debug::println("f", f); + + rd_alt_functor<I,I> fun(f, g); + o = morpho::canvas::reconstruction_on_set_union_find(f, g, nbh, fun); + + debug::println("o", o); +} Index: theo/mln/morpho/canvas/reconstruction_on_set.hh --- theo/mln/morpho/canvas/reconstruction_on_set.hh (revision 0) +++ theo/mln/morpho/canvas/reconstruction_on_set.hh (revision 0) @@ -0,0 +1,356 @@ +// Copyright (C) 2009 EPITA Research and Development Laboratory (LRDE) +// +// This file is part of Olena. +// +// Olena is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free +// Software Foundation, version 2 of the License. +// +// Olena 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 Olena. If not, see <http://www.gnu.org/licenses/>. +// +// As a special exception, you may use this file as part of a free +// software project 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_CANVAS_RECONSTRUCTION_ON_SET_UNION_FIND_HH +# define MLN_MORPHO_CANVAS_RECONSTRUCTION_ON_SET_UNION_FIND_HH + +# include <mln/core/concept/image.hh> +# include <mln/core/concept/neighborhood.hh> +# include <mln/morpho/canvas/internal/find_root.hh> + +# include <mln/data/fill.hh> +# include <mln/border/equalize.hh> +# include <mln/border/fill.hh> + + + +namespace mln +{ + + namespace morpho + { + + namespace canvas + { + + + template <typename I, typename J, typename N, + typename F> + mln_concrete(I) + reconstruction_on_set_union_find(const Image<I>& f, const Image<J>& g, const Neighborhood<N>& nbh, + F& functor); + + + + +# ifndef MLN_INCLUDE_ONLY + + + // Implementations. + + namespace impl + { + + namespace generic + { + + template <typename I, typename J, typename N, + typename F> + inline + mln_concrete(I) + reconstruction_on_set_union_find(const Image<I>& f_, const Image<J>& g_, const Neighborhood<N>& nbh_, + F& functor) + { + trace::entering("morpho::canvas::impl::generic::reconstruction_on_set_union_find"); + + const I& f = exact(f_); + const J& g = exact(g_); + const N& nbh = exact(nbh_); + + typedef mln_psite(I) P; + typedef mln_value(I) V; + + // Auxiliary data. + mln_ch_value(I, bool) deja_vu; + mln_ch_value(I, P) parent; + mln_concrete(I) output; + + // Initialization. + { + initialize(output, f); + functor.output = output; + + initialize(parent, f); + initialize(deja_vu, f); + data::fill(deja_vu, false); + + functor.set_default(); // <-- set_default + } + + // First pass. + { + mln_fwd_piter(I) p(f.domain()); + mln_niter(N) n(nbh, p); + for_all(p) + if (functor.domain(p)) // <-- domain + { + // Make-Set. + parent(p) = p; + functor.init(p); // <-- init + + for_all(n) if (f.domain().has(n)) + { + // See below (*) + if (functor.domain(n)) // <-- domain + { + if (deja_vu(n)) + { + // Do-Union. + P r = internal::find_root(parent, n); + if (r != p) + { + parent(r) = p; + functor.merge(p, r); // <-- merge + } + } + } + else + { + mln_invariant(deja_vu(n) == false); + functor.border(p, n); // <-- border + } + } + deja_vu(p) = true; + } + } + + // Second pass. + { + mln_bkd_piter(I) p(f.domain()); + for_all(p) + if (functor.domain(p)) // <-- domain + if (parent(p) != p) + output(p) = output(parent(p)); + } + + // (*) Technical note. + // n in D' + // | true | false | + // deja_vu(n) +----------+--------------+ + // true | do-union | impossible | + // false | no-op | border(p, n) | + // + // The no-op is sound because /n/ will be processed as a "p in D'" later on. + + trace::exiting("morpho::canvas::impl::generic::reconstruction_on_set_union_find"); + return output; + } + + } // end of namespace morpho::canvas::impl::generic + + + + template <typename I, typename J, typename N, + typename F> + inline + mln_concrete(I) + reconstruction_on_set_union_find_fastest(const Image<I>& f_, const Image<J>& g_, const Neighborhood<N>& nbh_, + F& functor) + { + trace::entering("morpho::canvas::impl::reconstruction_on_set_union_find_fastest"); + + const I& f = exact(f_); + const J& g = exact(g_); + const N& nbh = exact(nbh_); + + typedef mln_psite(I) P; + typedef mln_value(I) V; + + // Auxiliary data. + mln_ch_value(I, bool) deja_vu; + mln_ch_value(I, unsigned) parent; + mln_concrete(I) output; + + // Initialization. + { + border::equalize(f, g, nbh.delta()); + functor.set_extension_(); // <-- set_extension + + initialize(output, f); + functor.output = output; + + initialize(parent, f); + + functor.set_default(); // <-- set_default + } + + util::array<int> dp = negative_offsets_wrt(f, nbh); + const unsigned n_nbhs = dp.nelements(); + + // First pass. + { + mln_fwd_pixter(const I) pxl(f); + for_all(pxl) + { + unsigned p = pxl.offset(); + if (functor.domain_(p)) // <-- domain + { + // Make-Set. + parent.element(p) = p; + functor.init_(p); // <-- init + + // Deja-vu part. + for (unsigned j = 0; j < n_nbhs; ++j) + { + unsigned n = p + dp[j]; + if (functor.domain_(n)) // <-- domain + { + mln_invariant(f.domain().has(f.point_at_index(n))); + // Do-Union. + unsigned r = internal::find_root_fastest(parent, n); + if (r != p) + { + // set_parent, i.e., union + parent.element(r) = p; + functor.merge_(p, r); // <-- merge + } + } + else + functor.border_(p, n); // <-- border + } + + // Non-(deja-vu) part. + if (functor.inspect_border_()) + for (unsigned j = 0; j < n_nbhs; ++j) + { + unsigned n = p - dp[j]; + if (! functor.domain_(n)) // <-- domain + functor.border_(p, n); // <-- border + } + + // (*) Technical note. + // n in D' + // | true | false | + // deja_vu(n) +----------+--------------+ + // true | do-union | border(p, n) | + // false | no-op | border(p, n) | + // + // In this 'fastest' version, deja-vu(n) does not + // mean that n is in D' so we have to handle the + // case "deja_vu(n) and n not in D'". + } + } + } + + // Second pass. + { + mln_bkd_pixter(const I) pxl(f); + for_all(pxl) + { + unsigned p = pxl.offset(); + if (functor.domain_(p)) // <-- domain + if (parent.element(p) != p) + output.element(p) = output.element(parent.element(p)); + } + } + + trace::exiting("morpho::canvas::impl::reconstruction_on_set_union_find_fastest"); + return output; + } + + + } // end of namespace mln::morpho::canvas::impl + + + + // Dispatch. + + namespace internal + { + + template <typename I, typename J, typename N, + typename F> + inline + mln_concrete(I) + reconstruction_on_set_union_find_dispatch(metal::false_, + const Image<I>& f, const Image<J>& g, const Neighborhood<N>& nbh, + F& functor) + { + return impl::generic::reconstruction_on_set_union_find(f, g, nbh, + functor); + } + + template <typename I, typename J, typename N, + typename F> + inline + mln_concrete(I) + reconstruction_on_set_union_find_dispatch(metal::true_, + const Image<I>& f, const Image<J>& g, const Neighborhood<N>& nbh, + F& functor) + { + return impl::reconstruction_on_set_union_find_fastest(f, g, nbh, + functor); + } + + template <typename I, typename J, typename N, + typename F> + inline + mln_concrete(I) + reconstruction_on_set_union_find_dispatch(const Image<I>& f, const Image<J>& g, const Neighborhood<N>& nbh, + F& functor) + { + enum { + is_fastest = (mlc_equal(mln_trait_image_speed(I), + trait::image::speed::fastest)::value && + mlc_equal(mln_trait_image_speed(J), + trait::image::speed::fastest)::value && + mln_is_simple_neighborhood(N)::value) + }; + return reconstruction_on_set_union_find_dispatch(metal::bool_<is_fastest>(), + f, g, nbh, + functor); + } + + } // end of namespace mln::morpho::canvas::internal + + + // Facade. + + template <typename I, typename J, typename N, + typename F> + inline + mln_concrete(I) + reconstruction_on_set_union_find(const Image<I>& f, const Image<J>& g, const Neighborhood<N>& nbh, + F& functor) + { + trace::entering("morpho::canvas::reconstruction_on_set_union_find"); + + mln_concrete(I) output; + output = internal::reconstruction_on_set_union_find_dispatch(f, g, nbh, functor); + + trace::exiting("morpho::canvas::reconstruction_on_set_union_find"); + return output; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::morpho::canvas + + } // end of namespace mln::morpho + +} // end of namespace mln + + +#endif // ! MLN_MORPHO_CANVAS_RECONSTRUCTION_ON_SET_UNION_FIND_HH Index: theo/mln/morpho/canvas/internal/find_root.hh --- theo/mln/morpho/canvas/internal/find_root.hh (revision 0) +++ theo/mln/morpho/canvas/internal/find_root.hh (revision 0) @@ -0,0 +1,89 @@ +// Copyright (C) 2009 EPITA Research and Development Laboratory (LRDE) +// +// This file is part of Olena. +// +// Olena is free software: you can redistribute it and/or modify it under +// the terms of the GNU General Public License as published by the Free +// Software Foundation, version 2 of the License. +// +// Olena 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 Olena. If not, see <http://www.gnu.org/licenses/>. +// +// As a special exception, you may use this file as part of a free +// software project 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_CANVAS_INTERNAL_FIND_ROOT_HH +# define MLN_MORPHO_CANVAS_INTERNAL_FIND_ROOT_HH + +# include <mln/core/concept/image.hh> + + + +namespace mln +{ + + namespace morpho + { + + namespace canvas + { + + namespace internal + { + + template <typename Par> + mln_psite(Par) + find_root(Par& parent, mln_psite(Par) x); + + template <typename Par> + unsigned + find_root_fastest(Par& parent, unsigned x); + + +# ifndef MLN_INCLUDE_ONLY + + template <typename Par> + inline + mln_psite(Par) + find_root(Par& parent, mln_psite(Par) x) + { + if (parent(x) == x) + return x; + else + return parent(x) = find_root(parent, parent(x)); + } + + template <typename Par> + inline + unsigned + find_root_fastest(Par& parent, unsigned x) + { + if (parent.element(x) == x) + return x; + else + return parent.element(x) = find_root_fastest(parent, parent.element(x)); + } + + } // end of namespace mln::morpho::canvas::internal + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::morpho::canvas + + } // end of namespace mln::morpho + +} // end of namespace mln + + +#endif // ! MLN_MORPHO_CANVAS_INTERNAL_FIND_ROOT_HH Index: theo/mln/morpho/canvas/reconstruction_on_set.cc --- theo/mln/morpho/canvas/reconstruction_on_set.cc (revision 0) +++ theo/mln/morpho/canvas/reconstruction_on_set.cc (revision 0) @@ -0,0 +1,376 @@ +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/neighb2d.hh> +#include <mln/debug/println.hh> +#include <mln/io/pbm/load.hh> +#include <mln/io/pbm/save.hh> +#include <mln/data/compare.hh> +#include <mln/util/timer.hh> + +#include "reconstruction_on_set.hh" + + +// by_dilation + +template <typename F, typename G> +struct rd_functor +{ + const F& f; + const G& g; + + rd_functor(const F& f, const G& g) + : f(f), g(g) + { + } + + mln_concrete(F) output; + typedef mln_psite(F) P; + + void set_default() + { + mln::data::fill(output, false); + } + + bool domain(const P& p) const + { + return g(p) == true; + } + + void init(const P& p) + { + output(p) = f(p); + } + + void merge(const P& p, const P& r) + { + output(p) = output(p) || output(r); + } + + void border(const P&, const P&) + { + } + + // Fastest versions. + void set_extension_() + { + mln::border::fill(g, false); + } + bool domain_(unsigned p) const + { + return g.element(p) == true; + } + void init_(unsigned p) + { + output.element(p) = f.element(p); + } + void merge_(unsigned p, unsigned r) + { + output.element(p) = output.element(p) || output.element(r); + } + bool inspect_border_() const + { + return false; + } + void border_(unsigned, unsigned) + { + } +}; + + + + + +template <typename F, typename G> +struct rd_alt_functor // /!\ alternative version +{ + const F& f; + const G& g; + + rd_alt_functor(const F& f, const G& g) + : f(f), g(g) + { + } + + mln_concrete(F) output; + typedef mln_psite(F) P; + + void set_default() + { + mln::data::fill(output, f); + } + + bool domain(const P& p) const + { + return g(p) && ! f(p); + } + + void init(const P& p) + { + output(p) = false; + } + + void merge(const P& p, const P& r) + { + output(p) = output(p) || output(r); + } + + void border(const P& p, const P& n) + { + if (f(n)) + output(p) = true; + } + + // Fastest versions. + void set_extension_() + { + mln::border::fill(f, false); + mln::border::fill(g, false); + } + bool domain_(unsigned p) const + { + return g.element(p) && ! f.element(p); + } + void init_(unsigned p) + { + output.element(p) = false; + } + void merge_(unsigned p, unsigned r) + { + output.element(p) = output.element(p) || output.element(r); + } + bool inspect_border_() const + { + return true; + } + void border_(unsigned p, unsigned n) + { + if (f.element(n)) + { + mln_invariant(f.domain().has(f.point_at_index(n))); + output.element(p) = true; + } + } +}; + + + +// by_erosion + +template <typename F, typename G> +struct re_functor +{ + const F& f; + const G& g; + + re_functor(const F& f, const G& g) + : f(f), g(g) + { + } + + mln_concrete(F) output; + typedef mln_psite(F) P; + + void set_default() + { + mln::data::fill(output, true); + } + + bool domain(const P& p) const + { + return g(p) == false; + } + + void init(const P& p) + { + output(p) = f(p); + } + + void merge(const P& p, const P& r) + { + output(p) = output(p) && output(r); + } + + void border(const P&, const P&) + { + } + + // Fastest versions. + void set_extension_() + { + mln::border::fill(g, true); + } + bool domain_(unsigned p) const + { + return g.element(p) == false; + } + void init_(unsigned p) + { + output.element(p) = f.element(p); + } + void merge_(unsigned p, unsigned r) + { + output.element(p) = output.element(p) && output.element(r); + } + bool inspect_border_() const + { + return false; + } + void border_(unsigned, unsigned) + { + } +}; + + + + +template <typename F, typename G> +struct re_alt_functor +{ + const F& f; + const G& g; + + re_alt_functor(const F& f, const G& g) + : f(f), g(g) + { + } + + mln_concrete(F) output; + typedef mln_psite(F) P; + + void set_default() + { + mln::data::fill(output, g); + } + + bool domain(const P& p) const + { + return f(p) && ! g(p); + } + + void init(const P& p) + { + output(p) = true; + } + + void merge(const P& p, const P& r) + { + output(p) = output(p) && output(r); + } + + void border(const P& p, const P& n) + { + if (! f(n)) + output(p) = false; + } + + // Fastest versions. + void set_extension_() + { + mln::border::fill(f, true); + mln::border::fill(g, true); + } + bool domain_(unsigned p) const + { + return f.element(p) && ! g.element(p); + } + void init_(unsigned p) + { + output.element(p) = true; + } + void merge_(unsigned p, unsigned r) + { + output.element(p) = output.element(p) && output.element(r); + } + bool inspect_border_() const + { + return true; + } + void border_(unsigned p, unsigned n) + { + if (! f.element(n)) + output.element(p) = false; + } +}; + + + + +int main() +{ + using namespace mln; + + typedef image2d<bool> I; + + I f, g, ref, o; + neighb2d nbh = c4(); + + io::pbm::load(f, "f_and_g.pbm"); + io::pbm::load(g, "g.pbm"); + + util::timer t; + + + // by dilation + { + rd_functor<I,I> fun(f, g); + { + t.start(); + ref = morpho::canvas::impl::generic::reconstruction_on_set_union_find(f, g, nbh, fun); + std::cout << t.stop() << std::endl; + } + { + t.start(); + o = morpho::canvas::impl::reconstruction_on_set_union_find_fastest(f, g, nbh, fun); + std::cout << t.stop() << std::endl; + std::cout << (o == ref ? "OK" : "KO") << std::endl; + } + + rd_alt_functor<I,I> fun_alt(f, g); + { + t.start(); + o = morpho::canvas::impl::generic::reconstruction_on_set_union_find(f, g, nbh, fun_alt); + std::cout << t.stop() << std::endl; + std::cout << (o == ref ? "OK" : "KO") << std::endl; + } + { + t.start(); + o = morpho::canvas::impl::reconstruction_on_set_union_find_fastest(f, g, nbh, fun_alt); + std::cout << t.stop() << std::endl; + std::cout << (o == ref ? "OK" : "KO") << std::endl; + } + + io::pbm::save(ref, "rd.pbm"); + } + + // by erosion + { + re_functor<I,I> fun(g, f); + { + t.start(); + ref = morpho::canvas::impl::generic::reconstruction_on_set_union_find(g, f, nbh, fun); + std::cout << t.stop() << std::endl; + } + { + t.start(); + o = morpho::canvas::impl::reconstruction_on_set_union_find_fastest(g, f, nbh, fun); + std::cout << t.stop() << std::endl; + std::cout << (o == ref ? "OK" : "KO") << std::endl; + } + + re_alt_functor<I,I> fun_alt(g, f); + { + t.start(); + o = morpho::canvas::impl::generic::reconstruction_on_set_union_find(g, f, nbh, fun_alt); + std::cout << t.stop() << std::endl; + std::cout << (o == ref ? "OK" : "KO") << std::endl; + } + { + t.start(); + o = morpho::canvas::impl::reconstruction_on_set_union_find_fastest(g, f, nbh, fun_alt); + std::cout << t.stop() << std::endl; + std::cout << (o == ref ? "OK" : "KO") << std::endl; + } + + io::pbm::save(ref, "re.pbm"); + } + +}