[Olena-patches] milena r4554: Add fast version of dual input max tree computation and fix

URL: https://svn.lrde.epita.fr/svn/oln/trunk/milena ChangeLog: 2009-09-25 Edwin Carlinet &lt;carlinet(a)lrde.epita.fr&gt; Add fast version of dual input max tree computation and fix. * morpho/tree/all.hh: Update. * morpho/tree/compute_parent_dual_input.hh: Remove. * morpho/tree/data.hh: Fix constructors. * morpho/tree/dual_input_tree.hh: Add dispatches for generic / low quantized algorithms. * morpho/tree/impl/dual_hqueue.hh: New. * morpho/tree/impl/dual_union_find.hh: New. * morpho/tree/impl: New. --- all.hh | 1 data.hh | 105 ++++-------- dual_input_tree.hh | 75 +++++++- impl/dual_hqueue.hh | 408 ++++++++++++++++++++++++++++++++++++++++++++++++ impl/dual_union_find.hh | 361 ++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 874 insertions(+), 76 deletions(-) Index: trunk/milena/mln/morpho/tree/compute_parent_dual_input.hh (deleted) =================================================================== Index: trunk/milena/mln/morpho/tree/dual_input_tree.hh =================================================================== --- trunk/milena/mln/morpho/tree/dual_input_tree.hh (revision 4553) +++ trunk/milena/mln/morpho/tree/dual_input_tree.hh (revision 4554) @@ -25,14 +25,16 @@ /// \file /// -/// Compute a canonized (parenthood) component tree from a dual input. - +/// Compute a canonized component tree from a dual input. #ifndef MLN_MORPHO_TREE_DUAL_INPUT_TREE_HH # define MLN_MORPHO_TREE_DUAL_INPUT_TREE_HH # include <mln/data/sort_psites.hh> + # include <mln/morpho/tree/data.hh> +# include <mln/morpho/tree/impl/dual_union_find.hh> +# include <mln/morpho/tree/impl/dual_hqueue.hh> namespace mln { @@ -43,44 +45,95 @@ namespace tree { - template <typename I, typename M, typename N> + /// Compute the dual input max tree using mask-based connectivity. + /// + /// \param[in] f The original image. + /// \param[in] m The connectivity mask. + /// \param[in] nbh The neighborhood of the mask. + /// + /// \return The computed tree. + /// + template <typename I, typename N> inline - morpho::tree::data< I, p_array<mln_psite(I)> > + data< I, p_array<mln_psite(I)> > dual_input_max_tree(const Image<I>& f, - const Image<M>& m, + const Image<I>& m, const Neighborhood<N>& nbh); # ifndef MLN_INCLUDE_ONLY - template <typename I, typename M, typename N> + namespace internal + { + + template <typename I, typename N> + inline + data< I, p_array<mln_psite(I)> > + dual_input_max_tree_dispatch(trait::image::quant::any, + const I& f, + const I& m, + const N& nbh) + { + typedef p_array<mln_psite(I)> S; + typedef data<I,S> tree_t; + + S s_f = mln::data::sort_psites_increasing(f); + S s_m = mln::data::sort_psites_increasing(m); + + tree_t tree = impl::generic::dual_union_find(f, m, s_f, s_m, nbh); + return tree; + } + + template <typename I, typename N> + inline + data< I, p_array<mln_psite(I)> > + dual_input_max_tree_dispatch(trait::image::quant::low, + const I& f, + const I& m, + const N& nbh) + { + typedef p_array<mln_psite(I)> S; + typedef data<I,S> tree_t; + + tree_t tree = impl::dual_hqueue(f, m, nbh); + return tree; + } + + } // end of namespace mln::morpho::tree::internal + + + // Facades. + template <typename I, typename N> inline morpho::tree::data< I, p_array<mln_psite(I)> > dual_input_max_tree(const Image<I>& f_, - const Image<M>& m_, + const Image<I>& m_, const Neighborhood<N>& nbh_) { trace::entering("morpho::tree::dual_input_max_tree"); const I& f = exact(f_); - const M& m = exact(m_); + const I& m = exact(m_); const N& nbh = exact(nbh_); mln_precondition(f.is_valid()); mln_precondition(m.is_valid()); mln_precondition(nbh.is_valid()); + mln_precondition(f.domain() == m.domain()); typedef p_array<mln_psite(I)> S; typedef data<I,S> tree_t; - S s_f = mln::data::sort_psites_increasing(f); - S s_m = mln::data::sort_psites_increasing(m); - tree_t tree(f, m, s_f, s_m, nbh); + tree_t tree = internal::dual_input_max_tree_dispatch(mln_trait_image_quant(I)(), f, m, nbh); trace::exiting("morpho::tree::dual_input_max_tree"); return tree; } + + + + # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::morpho::tree Index: trunk/milena/mln/morpho/tree/impl/dual_hqueue.hh =================================================================== --- trunk/milena/mln/morpho/tree/impl/dual_hqueue.hh (revision 0) +++ trunk/milena/mln/morpho/tree/impl/dual_hqueue.hh (revision 4554) @@ -0,0 +1,408 @@ +// Copyright (C) 2008, 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_TREE_IMPL_COMPUTE_PARENT_DUAL_HQUEUE_HH +# define MLN_MORPHO_TREE_IMPL_COMPUTE_PARENT_DUAL_HQUEUE_HH + +/// +/// \brief The dual input tree algorithm specialized for low quantized image. +/// +/// This implementation is based on P. Salembier algorithm using +/// hierachical queues. This implies a low-quantized input image so +/// that the number of queues is limited. +/// +/// TODO: Think about how to extend f domain in a more +/// generic way. The actual implementation doubles the size of the +/// first dimension. It implies a boxed domain. +/// +/// TODO: Use the less functor. The actual implementation is for max-tree. +/// +/// TODO: During the canonization pass, we build the tree site set +/// from the sorted site set of f, so that we compute twice f +/// histogram (can be avoided). + +# include <mln/data/sort_psites.hh> +# include <mln/data/fill.hh> + +# include <mln/geom/nsites.hh> +# include <mln/geom/translate.hh> + +# include <mln/morpho/tree/data.hh> + +# include <mln/util/hqueues.hh> +# include <mln/util/ord.hh> + +# include <mln/value/value_array.hh> +# include <mln/value/set.hh> + +# include <mln/util/timer.hh> + +namespace mln +{ + + namespace morpho + { + + namespace tree + { + + namespace impl + { + + /// Compute a tree using hqueues. + /// + /// \param f The original image. + /// \param m The connectivity mask. + /// \param nbh The neighborhood of the mask. + /// + /// \return The tree structure. + /// + template <typename I, typename N> + inline + data<I, p_array<mln_psite(I)> > + dual_hqueue(const Image<I>& f, + const Image<I>& m, + const Neighborhood<N>& nbh); + + } // end of namespace mln::morpho::tree::impl + + +# ifndef MLN_INCLUDE_ONLY + + namespace internal + { + + template <typename I, typename N, class E> + struct shared_flood_args + { + typedef mln_psite(I) P; + typedef mln_value(I) V; + typedef p_array<P> S; + + const I& f; + const I& m; + const N& nbh; + mln_ch_value(I, P)& parent; + + // Aux data + util::hqueues<P, V>& hqueues; + const E& extend; // site -> site functor. + value::value_array<V, bool> is_node_at_level; + value::value_array<V, P> node_at_level; + mln_ch_value(I, bool) deja_vu; + const value::set<V>& vset; + + shared_flood_args(const I& f_, + const I& m_, + const N& neigb_, + mln_ch_value(I, P)& parent_, + util::hqueues<mln_psite(I), V>& hqueues_, + const E& extend_) + : f (f_), + m (m_), + nbh (neigb_), + parent (parent_), + hqueues (hqueues_), + extend (extend_), + is_node_at_level (false), + vset (hqueues.vset()) + { + initialize(deja_vu, f); + mln::data::fill(deja_vu, false); + } + }; + + template <typename I> + inline + histo::array<mln_value(I)> + compute_histo(const I& f, const I& m, mln_value(I)& hmin, mln_psite(I)& pmin) + { + histo::array<mln_value(I)> hm = histo::compute(m); + const histo::array<mln_value(I)> hf = histo::compute(f); + + { // Retrieve hmin. + unsigned i = 0; + while (hm[i] == 0) + ++i; + hmin = hm.vset()[i]; + } + + // Merge histograms. + for (unsigned i = 0; i < hm.nvalues(); ++i) + hm[i] += hf[i]; + + // Retrieve pmin. + mln_piter(I) p(m.domain()); + for (p.start(); m(p) != hmin; p.next()) + ; + mln_assertion(p.is_valid()); + pmin = p; + + return hm; + } + + // Site -> site functor: give for all p in Domain(f), its + // equivalence in the extended domain. + // TODO: make it generic. It works only on boxed domain. + template <typename I> + struct extend + { + extend(const mln_psite(I)::delta& dp) + : dp_ (dp) + { + } + + mln_psite(I) operator() (const mln_psite(I)& p) const + { + return p + dp_; + } + + private: + const mln_psite(I)::delta dp_; + }; + + } // end of namespace mln::morpho::tree::internal + + namespace impl + { + + template <typename I, typename N, typename E> + unsigned + flood(internal::shared_flood_args<I, N, E>& args, const unsigned h_idx) + { + mln_assertion(args.is_node_at_level[h_idx]); + + while (!args.hqueues[h_idx].empty()) + { + mln_psite(I) p = args.hqueues[h_idx].pop_front(); + unsigned p_idx = args.vset.index_of(args.f(p)); + + if (p_idx != h_idx) + { // Intensity mismatch: irregular case. + mln_psite(I) pext = args.extend(p); + args.parent(pext) = args.node_at_level[h_idx]; + + if (p_idx > h_idx) // Singleton with parent at h. + args.parent(p) = args.node_at_level[h_idx]; + else + { + if (!args.is_node_at_level[p_idx]) + { + args.is_node_at_level[p_idx] = true; + args.node_at_level[p_idx] = p; + } + } + } + + if (p_idx <= h_idx) + { + if (!args.f.domain().has(args.node_at_level[p_idx]) || + util::ord_strict(p, args.node_at_level[p_idx])) + { // Regular case but a representative provided by the extension. + args.parent(args.node_at_level[p_idx]) = p; + args.node_at_level[p_idx] = p; + //args.parent(p) = p; + } + args.parent(p) = args.node_at_level[p_idx]; + } + + + // Process the neighbors + mln_niter(N) n(args.nbh, p); + for_all(n) + if (args.f.domain().has(n) && !args.deja_vu(n)) + { + unsigned mn = args.vset.index_of(args.m(n)); + unsigned fn = args.vset.index_of(args.f(n)); + args.hqueues[mn].push(n); + args.deja_vu(n) = true; + + mln_psite(I) ext = args.extend(n); + // Create a node at c. + { + mln_psite(I) node = (fn == mn) ? n : ext; + if (!args.is_node_at_level[mn]) + { + args.is_node_at_level[mn] = true; + args.node_at_level[mn] = node; + } + } + + while (mn > h_idx) + mn = flood(args, mn); + } + } + + // Retrieve dad. + args.is_node_at_level[h_idx] = false; + unsigned c = h_idx; + while (c > 0 && !args.is_node_at_level[c]) + --c; + + mln_psite(I) x = args.node_at_level[h_idx]; + if (c > 0) + args.parent(x) = args.node_at_level[c]; + else + args.parent(x) = x; + + return c; + } + + template <typename I, typename N> + inline + data< I, p_array<mln_psite(I)> > + dual_hqueue(const Image<I>& f_, + const Image<I>& m_, + const Neighborhood<N>& neibh_) + { + trace::entering("mln::morpho::tree::impl::dual_hqueue"); + + const I& f = exact(f_); + const I& m = exact(m_); + const N& nbh = exact(neibh_); + + typedef mln_psite(I) P; + typedef p_array<mln_psite(I)> S; + + util::timer tm; + tm.start(); + + // Histo. + mln_psite(I) pmin; + mln_value(I) hmin; + const histo::array<mln_value(I)> histo = internal::compute_histo(f, m, hmin, pmin); + util::hqueues<P, mln_value(I)> hqueues(histo); + + mln_psite(I)::delta dp(literal::zero); + mln_domain(I) d_ext = f.domain(); + mln_domain(I) d = f.domain(); + + // Extend the domain. + dp[0] = d.pmax()[0] - d.pmin()[0] + 1; + d.pmax() += dp; + d_ext.pmin() += dp; + d_ext.pmax() += dp; + + // Data. + mln_concrete(I) fext; + mln_ch_value(I, P) parent; + p_array<mln_psite(I)> s; + + // Initialization. + fext = geom::translate(m, dp.to_vec(), f, d); + initialize(parent, fext); + s.reserve(geom::nsites(fext)); + + // Process. + internal::extend<I> extend(dp); + internal::shared_flood_args< I, N, internal::extend<I> > + args(f, m, nbh, parent, hqueues, extend); + + unsigned r = args.vset.index_of(hmin); + args.deja_vu(pmin) = true; + args.hqueues[r].push(pmin); + args.node_at_level[r] = (f(pmin) == hmin) ? pmin : extend(pmin); + args.is_node_at_level[r] = true; + flood(args, r); + + // Attach the nodes under hmin. + unsigned i = r; + do + { + if (args.is_node_at_level[i]) + { + parent(args.node_at_level[r]) = args.node_at_level[i]; + r = i; + } + } + while (i-- > 0); + parent(args.node_at_level[r]) = args.node_at_level[r]; //root + + // Canonization and make tree site set. + { + mln_ch_value(I, bool) deja_vu(d_ext); + mln::data::fill(deja_vu, false); + + p_array<mln_psite(I)> s_f = mln::data::sort_psites_increasing(f); + mln_fwd_piter(S) p(s_f); // Forward. + for_all(p) + { + P x = p; + P q = parent(p); + + // Case: l: m <---- m <---- f + // Or + // Case l1: m <----- f impossible. + // | + // l2: m + mln_assertion(!(d_ext.has(q) && fext(p) == fext(q) && d_ext.has(parent(q)) && q != parent(q))); + + while (d_ext.has(q) && !deja_vu(q) && (fext(q) != fext(parent(q)) || q == parent(q))) + { + s.append(q); + deja_vu(q) = true; + x = q; + q = parent(q); + } + + if (d_ext.has(q) && fext(q) == fext(parent(q)) && q != parent(q)) + { + q = (parent(x) = parent(q)); + mln_assertion(f.domain().has(q)); + } + + if (fext(q) == fext(parent(q))) + parent(x) = parent(q); + + s.append(p); + + mln_assertion((q = parent(p)) == parent(q) || fext(q) != fext(parent(q))); + } + + } + + std::cout << "Construction de l'arbre en " << tm << " s." << std::endl; + + data<I, S> tree(fext, parent, s); + + trace::exiting("mln::morpho::tree::impl::dual_hqueue"); + + return tree; + } + + } // end of namespace mln::morpho::tree::impl + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::morpho::tree + + } // end of namespace mln::morpho + +} // end of namespace mln + +#endif // !MLN_MORPHO_TREE_COMPUTE_PARENT_DUAL_HQUEUE_HH + + + Index: trunk/milena/mln/morpho/tree/impl/dual_union_find.hh =================================================================== --- trunk/milena/mln/morpho/tree/impl/dual_union_find.hh (revision 0) +++ trunk/milena/mln/morpho/tree/impl/dual_union_find.hh (revision 4554) @@ -0,0 +1,361 @@ +// Copyright (C) 2008, 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_TREE_IMLP_DUAL_UNION_FIND_HH +# define MLN_MORPHO_TREE_IMLP_DUAL_UNION_FIND_HH + +/// +/// \brief The generic dual input tree algorithm for high quantized image. +/// +/// This implementation is based on tarjan's union method, so that +/// image quantization does not impact on the computation time. +/// +/// TODO: Think about how to extend f domain in a more +/// generic way. The actual implementation doubles the size of the +/// first dimension. It implies a boxed domain. +/// +/// TODO: Use the less functor. The actual implementation is for max-tree. + +# include <mln/data/fill.hh> + +# include <mln/geom/nsites.hh> +# include <mln/geom/translate.hh> + +# include <mln/morpho/tree/data.hh> + +# include <mln/util/timer.hh> + +namespace mln +{ + + namespace morpho + { + + namespace tree + { + + namespace impl + { + + namespace generic + { + + /// Compute a tree using union-find method. + /// + /// \param f The original image. + /// \param m The connectivity mask. + /// \param s_f The sorted site set of \p f + /// \param s_m The sorted site set of \p m. + /// \param nbh The neighborhood of the mask. + /// + /// \return The tree structure. + /// + template <typename I, typename S, typename N> + data< I, p_array<mln_psite(I)> > + dual_union_find(const Image<I>& f, + const Image<I>& m, + const Site_Set<S>& s_f, + const Site_Set<S>& s_m, + const Neighborhood<N>& nbh); + + } // end of namespace mln::morpho::tree::impl::generic + + } // end of namespace mln::morpho::tree::impl + +# ifndef MLN_INCLUDE_ONLY + + namespace internal + { + // For benchmark purpose. More than 50% of the time is spent + // in find_root function. + static util::timer t_prop; + + + template <typename I> + inline + mln_psite(I) find_root(I& zpar, + const mln_psite(I)& p) + { + if (zpar(p) == p) + return p; + + t_prop.resume(); + mln_psite(I) x = zpar(p); + while (zpar(x) != x) + x = zpar(x); + + mln_psite(I) tmp; + for (mln_psite(I) y = p; y != x; y = tmp) + { + tmp = zpar(y); + zpar(y) = x; + } + t_prop.stop(); + + return x; + } + + template <typename I> + inline + void + update_m_parent(const I& f, + mln_ch_value(I, mln_psite(I))& parent, + mln_ch_value(I, bool)& deja_vu, + p_array< mln_psite(I) >& sset, + const mln_domain(I)& d_ext, + const mln_psite(I)& p) + { + typedef mln_psite(I) P; + + P q = parent(p); + P x = parent(q); + + mln_assertion(d_ext.has(q)); + + while (d_ext.has(x) && f(q) == f(x) && q != x) + { + q = x; + x = parent(q); + } + + if (!d_ext.has(x)) + { + if (f(x) == f(parent(x))) + x = (parent(q) = parent(x)); + if (f(q) != f(x)) + { + x = q; + if (!deja_vu(q)) + { + deja_vu(q) = true; + sset.append(q); + } + } + + } + else + { + if (x != q) + { + update_m_parent(f, parent, deja_vu, sset, d_ext, q); + x = q; + } + if (!deja_vu(q)) + { + deja_vu(q) = true; + sset.append(q); + } + } + + for (P i = p, tmp = parent(i); i != q; i = tmp, tmp = parent(i)) + parent(i) = x; + } + + } // end of namespace mln::morpho::tree::internal + + namespace impl + { + + namespace generic + { + + + template <typename I, typename S, typename N> + data< I, p_array<mln_psite(I)> > + dual_union_find(const Image<I>& f_, + const Image<I>& m_, + const Site_Set<S>& s_f_, + const Site_Set<S>& s_m_, + const Neighborhood<N>& nbh_) + { + trace::entering("morpho::tree::impl::generic::dual_union_find"); + + util::timer tm; + tm.start(); + internal::t_prop.reset(); + + typedef mln_psite(I) P; + typedef unsigned L; + + const I& f = exact(f_); + const I& m = exact(m_); + const S& s_f = exact(s_f_); + const S& s_m = exact(s_m_); + const N& nbh = exact(nbh_); + + // Aux data. + mln_psite(I)::delta dp(literal::zero); + mln_domain(I) d_f = f.domain(); + mln_domain(I) d_ext = f.domain(); // translate dp + mln_domain(I) d = f.domain(); + + // Extend the domain. + dp[0] = d.pmax()[0] - d.pmin()[0] + 1; + d.pmax() += dp; + d_ext.pmin() += dp; + d_ext.pmax() += dp; + + // Data. + mln_concrete(I) fext; + mln_ch_value(I, P) parent; + p_array<mln_psite(I)> s; + + // Initialization. + fext = geom::translate(m, dp.to_vec(), f, d); + initialize(parent, fext); + s.reserve(geom::nsites(fext)); + + //Process + { + // Aux data. + mln_ch_value(I, bool) deja_vu; + mln_ch_value(I, P) zpar; + + initialize(deja_vu, fext); + initialize(zpar, fext); + mln::data::fill(deja_vu, false); + + mln_bkd_piter(S) p_f(s_f); // Backward. + mln_bkd_piter(S) p_m(s_m); // Backward. + p_f.start(); + p_m.start(); + + // Main loop. + while (p_m.is_valid() || p_f.is_valid()) + { + mln_bkd_piter(S)& it = (!p_f.is_valid() || (p_m.is_valid() && f(p_f) <= m(p_m))) ? p_m : p_f; + + P p = it; + P ext = p + dp; + + // std::cout << "-------------------" << std::endl; + // std::cout << "Take " << p << " of value " << (&it == &p_m ? m(p) : f(p)) + // << " from " << (&it == &p_m ? "mask" : "f") << std::endl; + // debug::println("Parent: ", parent); + // debug::println("Zpar: ", zpar); + + mln_assertion(!(deja_vu(p) && deja_vu(ext))); + if (deja_vu(ext)) // Seen by mask before f. + { + mln_assertion(m(p) >= f(p)); + // Make set. + parent(p) = p; + zpar(p) = p; + + P r = internal::find_root(zpar, ext); + zpar(r) = p; + parent(r) = p; + + deja_vu(p) = true; + } + else if (deja_vu(p)) // Seen by f before mask. + { + mln_assertion(f(p) > m(p)); + parent(p) = ext; + zpar(p) = ext; + parent(ext) = ext; + zpar(ext) = ext; + + mln_niter(N) n(nbh, ext); + for_all(n) + if (d_ext.has(n) && deja_vu(n)) + { + P r = internal::find_root(zpar, n); + //std::cout << "Root: " << r << std::endl; + if (r != ext) + { + parent(r) = ext; + zpar(r) = ext; + } + } + deja_vu(ext) = true; + } + else if (f(p) <= m(p)) // First time p encountered. + { + zpar(ext) = ext; + mln_niter(N) n(nbh, ext); + for_all(n) + if (d_ext.has(n) && deja_vu(n)) + { + P r = internal::find_root(zpar, n); + if (r != ext) + { + zpar(r) = ext; + parent(r) = ext; + } + } + deja_vu(ext) = true; + } + else + { + deja_vu(p) = true; + } + it.next(); + } + } + std::cout << ">> MAJ zpar: " << internal::t_prop << " s" << std::endl; + std::cout << "Parent construction: " << tm << " s" << std::endl; + tm.restart(); + + // Canonization + { + mln_ch_value(I, bool) deja_vu(d_ext); + mln::data::fill(deja_vu, false); + mln_fwd_piter(S) p(s_f); // Forward. + for_all(p) + { + P q = parent(p); + if (!f.domain().has(q)) + internal::update_m_parent(fext, parent, deja_vu, s, d_ext, p); + else if (fext(parent(q)) == f(q)) + parent(p) = parent(q); + s.append(p); + + mln_assertion((q = parent(p)) == parent(q) || fext(q) != fext(parent(q))); + } + } + std::cout << "Canonization: " << tm << " s" << std::endl; + + //mln_postcondition(internal::compute_parent_postconditions(fext, s, parent)); + + tree::data<I, p_array<mln_psite(I)> > tree(fext, parent, s); + trace::exiting("morpho::tree::impl::generic::dual_union_find"); + + return tree; + } + + } // end of namespace mln::morpho::tree::impl::generic + + } // end of namespace mln::morpho::tree::impl + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::morpho::tree + + } // end of namespace mln::morpho + +} // end of namespace mln + +#endif // !MLN_MORPHO_TREE_IMLP_DUAL_UNION_FIND_HH Index: trunk/milena/mln/morpho/tree/all.hh =================================================================== --- trunk/milena/mln/morpho/tree/all.hh (revision 4553) +++ trunk/milena/mln/morpho/tree/all.hh (revision 4554) @@ -48,6 +48,7 @@ # include <mln/morpho/tree/compute_attribute_image.hh> # include <mln/morpho/tree/compute_parent.hh> +# include <mln/morpho/tree/dual_input_tree.hh> # include <mln/morpho/tree/data.hh> # include <mln/morpho/tree/max.hh> # include <mln/morpho/tree/utils.hh> Index: trunk/milena/mln/morpho/tree/data.hh =================================================================== --- trunk/milena/mln/morpho/tree/data.hh (revision 4553) +++ trunk/milena/mln/morpho/tree/data.hh (revision 4554) @@ -32,7 +32,6 @@ /// of S container iterator) to go faster. # include <mln/morpho/tree/compute_parent.hh> -# include <mln/morpho/tree/compute_parent_dual_input.hh> # include <mln/core/site_set/p_array.hh> # include <mln/core/internal/site_set_iterator_base.hh> # include <mln/core/internal/piter_identity.hh> @@ -138,18 +137,18 @@ typedef mln::morpho::tree::depth1st_piter<self_> depth1st_piter; - /// Constructor. + /// Standard constructor. template <typename N> - data(const Image<I>& f, const Site_Set<S>& s, const Neighborhood<N>& nbh); - - /// Constructor. - template <typename N, typename M> data(const Image<I>& f, - const Image<M>& m, - const Site_Set<S>& s_f, - const Site_Set<S>& s_m, + const Site_Set<S>& s, const Neighborhood<N>& nbh); + /// Special constructor where the parent computation has + /// delegated to an external function. (To handle special case + /// of connectivity for example). + data(const Image<I>& f, + const parent_t& parent, + const Site_Set<S>& s); /// \{ Base function-related materials @@ -205,13 +204,18 @@ unsigned nroots() const; protected: - const function& f_; - const sites_t s_; + void compute_children_(); + + protected: mln_ch_value(I, mln_psite(I)) parent_; // Parent image. mln_ch_value(I, nodes_t) children_; // Children image. - nodes_t nodes_; - leaves_t leaves_; - unsigned nroots_; + + function f_; // f image containing values of the tree nodes. + sites_t s_; // Sorted site set of the tree sites. (domain(f_) includes s_). + + nodes_t nodes_; // Sorted node set. + leaves_t leaves_; // Sorted leaf set. + unsigned nroots_; // For non-contigous domain image purpose. }; @@ -403,43 +407,15 @@ { template <typename I, typename S> - template <typename N, typename M> inline data<I, S>::data(const Image<I>& f, - const Image<M>& m, - const Site_Set<S>& s_f, - const Site_Set<S>& s_m, - const Neighborhood<N>& nbh) - : f_(exact(f)), - s_(exact(s_f)) - { - const N& nbh_ = exact(nbh); - - // Compute parent image. - parent_ = morpho::tree::compute_parent_dual_input(f_, m, nbh_, s_, s_m); - initialize(children_, f); - - // Store tree nodes. - nroots_ = 0; - mln_bkd_piter(S) p(s_); - for_all(p) - { - if (f_(parent_(p)) != f_(p)) - { - nodes_.insert(p); - children_(parent_(p)).insert(p); - if (is_a_leaf(p)) - leaves_.insert(p); - } - else if (parent_(p) == p) //it's a root. + const mln_ch_value(I, mln_psite(I))& parent, + const Site_Set<S>& s) + : parent_ (parent), + f_ (exact(f)), + s_ (exact(s)) { - nodes_.insert(p); - if (is_a_leaf(p)) // One pixel image... - leaves_.insert(p); - ++nroots_; - } - } - mln_assertion(leaves_.nsites() > 0); + compute_children_(); } @@ -447,15 +423,25 @@ template <typename I, typename S> template <typename N> inline - data<I,S>::data(const Image<I>& f, const Site_Set<S>& s, const Neighborhood<N>& nbh) + data<I,S>::data(const Image<I>& f, + const Site_Set<S>& s, + const Neighborhood<N>& nbh) : f_(exact(f)), s_(exact(s)) { - const N& nbh_ = exact(nbh); - // Compute parent image. + const N& nbh_ = exact(nbh); parent_ = morpho::tree::compute_parent(f_, nbh_, s_); - initialize(children_, f); + + compute_children_(); + } + + template <typename I, typename S> + inline + void + data<I, S>::compute_children_() + { + initialize(children_, f_); // Store tree nodes. nroots_ = 0; @@ -480,24 +466,13 @@ mln_assertion(leaves_.nsites() > 0); } + template <typename I, typename S> inline mln_rvalue_(mln_ch_value(I, mln_psite(I))) data<I,S>::parent(const mln_psite(I)& p) const { - // HERE - -// mln_precondition(is_valid()); -// mln_rvalue(parent_t) tmp; -// if (! parent_.domain().has(p)) -// { -// std::cout << "KO next" << std::endl; -// std::cout << p << std::endl; -// std::cout << parent_.domain() << std::endl; -// } -// else -// std::cout << "ok next" << std::endl; -// mln_precondition(parent_.domain().has(p)); + mln_precondition(parent_.domain().has(p)); return parent_(p); }