* mln/core/image/all.hh: Update. * mln/morpho/tree/all.hh: Update. * mln/tag/skeleton.hh: Update. * mln/morpho/tree/component_tree.hh: Routines to build component trees. * mln/morpho/tree/compute_attribute_image.hh: Routines to compute attribute images. * mln/morpho/tree/dual_input_tree.hh: Routines to build dual-input trees. * mln/morpho/tree/filter/direct.hh, * mln/morpho/tree/filter/max.hh, * mln/morpho/tree/filter/min.hh, * mln/morpho/tree/filter/subtractive.hh: Filtering routines. * mln/morpho/tree/impl/hqueue.hh, * mln/morpho/tree/impl/hqueue_fast.hh, * mln/morpho/tree/impl/union_find.hh, * mln/morpho/tree/impl/union_find_fast.hh: Tree computation algorithms. * mln/morpho/tree/leaf_last.hh, * mln/morpho/tree/propagate_node.hh: Tree management routines. --- milena/mln/core/image/all.hh | 1 + milena/mln/morpho/tree/all.hh | 16 +- milena/mln/morpho/tree/component_tree.hh | 86 +++- milena/mln/morpho/tree/compute_attribute_image.hh | 281 ++++------ milena/mln/morpho/tree/dual_input_tree.hh | 79 +--- milena/mln/morpho/tree/filter/direct.hh | 35 +- milena/mln/morpho/tree/filter/max.hh | 62 ++- milena/mln/morpho/tree/filter/min.hh | 59 ++- milena/mln/morpho/tree/filter/subtractive.hh | 52 ++- milena/mln/morpho/tree/impl/dual_hqueue.hh | 604 ++++++++++++-------- milena/mln/morpho/tree/impl/hqueue.hh | 428 ++++++++++++++ milena/mln/morpho/tree/impl/hqueue_fast.hh | 391 +++++++++++++ .../edwin => }/mln/morpho/tree/impl/union_find.hh | 2 + milena/mln/morpho/tree/impl/union_find_fast.hh | 396 +++++++++++++ milena/mln/morpho/tree/leaf_last.hh | 111 ++++ milena/mln/morpho/tree/propagate_node.hh | 189 ++----- milena/mln/tag/skeleton.hh | 1 + milena/mln/util/hqueues.hh | 32 +- 18 files changed, 2050 insertions(+), 775 deletions(-) create mode 100644 milena/mln/morpho/tree/impl/hqueue.hh create mode 100644 milena/mln/morpho/tree/impl/hqueue_fast.hh copy milena/{sandbox/edwin => }/mln/morpho/tree/impl/union_find.hh (99%) create mode 100644 milena/mln/morpho/tree/impl/union_find_fast.hh create mode 100644 milena/mln/morpho/tree/leaf_last.hh diff --git a/milena/mln/core/image/all.hh b/milena/mln/core/image/all.hh index ee6a15e..1cd2773 100644 --- a/milena/mln/core/image/all.hh +++ b/milena/mln/core/image/all.hh @@ -40,6 +40,7 @@ // Files. +# include <mln/core/image/attribute_image.hh> # include <mln/core/image/ch_piter.hh> # include <mln/core/image/complex_image.hh> # include <mln/core/image/complex_neighborhood_piter.hh> diff --git a/milena/mln/morpho/tree/all.hh b/milena/mln/morpho/tree/all.hh index d088ea0..c6d93ee 100644 --- a/milena/mln/morpho/tree/all.hh +++ b/milena/mln/morpho/tree/all.hh @@ -30,7 +30,6 @@ /// /// File that includes all morphological tree-related routines. - namespace mln { @@ -46,11 +45,14 @@ namespace mln } // end of namespace mln -# 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> +# include <mln/morpho/tree/component_tree.hh> +# include <mln/morpho/tree/compute_parent.hh +# include <mln/morpho/tree/compute_attribute_image.hh> +# include <mln/morpho/tree/leaf_last.hh> +# include <mln/morpho/tree/propagate_node.hh> + +# include <mln/morpho/tree/filter/all.hh> + -#endif // ! MLN_MORPHO_TREE_ALL_HH +#endif // !MLN_MORPHO_TREE_ALL_HH diff --git a/milena/mln/morpho/tree/component_tree.hh b/milena/mln/morpho/tree/component_tree.hh index f45986d..dc9063d 100644 --- a/milena/mln/morpho/tree/component_tree.hh +++ b/milena/mln/morpho/tree/component_tree.hh @@ -32,8 +32,10 @@ #ifndef MLN_MORPHO_TREE_COMPONENT_TREE_HH # define MLN_MORPHO_TREE_COMPONENT_TREE_HH -# include <mln/data/sort_psites.hh> -# include <mln/morpho/tree/data.hh> +# include <mln/tag/tree.hh> +# include <mln/util/ctree/ctree.hh> +# include <mln/morpho/tree/impl/union_find_fast.hh> +# include <mln/morpho/tree/impl/hqueue_fast.hh> namespace mln { @@ -53,7 +55,7 @@ namespace mln /// template <typename I, typename N> inline - data< I, p_array<mln_psite(I)> > + util::ctree::ctree<I> min_tree(const Image<I>& f, const Neighborhood<N>& nbh); /// Compute a canonized max-tree. @@ -65,15 +67,73 @@ namespace mln /// template <typename I, typename N> inline - data< I, p_array<mln_psite(I)> > + util::ctree::ctree<I> max_tree(const Image<I>& f, const Neighborhood<N>& nbh); # ifndef MLN_INCLUDE_ONLY + namespace internal + { + + template <typename I, typename N> + inline + util::ctree::ctree<I> + max_tree_dispatch(trait::image::quant::any, + const I& f, + const N& nbh) + { + util::ctree::ctree<I> tree = + morpho::tree::impl::generic::union_find_fast(tag::tree::max, f, nbh); + + return tree; + } + + template <typename I, typename N> + inline + util::ctree::ctree<I> + min_tree_dispatch(trait::image::quant::any, + const I& f, + const N& nbh) + { + util::ctree::ctree<I> tree = + morpho::tree::impl::generic::union_find_fast(tag::tree::min, f, nbh); + + return tree; + } + + + template <typename I, typename N> + inline + util::ctree::ctree<I> + max_tree_dispatch(trait::image::quant::low, + const I& f, + const N& nbh) + { + util::ctree::ctree<I> tree = + morpho::tree::impl::hqueue_fast(tag::tree::max, f, nbh); + + return tree; + } + + template <typename I, typename N> + inline + util::ctree::ctree<I> + min_tree_dispatch(trait::image::quant::low, + const I& f, + const N& nbh) + { + util::ctree::ctree<I> tree = + morpho::tree::impl::hqueue_fast(tag::tree::min, f, nbh); + + return tree; + } + + } + template <typename I, typename N> inline - data< I, p_array<mln_psite(I)> > + util::ctree::ctree<I> min_tree(const Image<I>& f_, const Neighborhood<N>& nbh_) { trace::entering("morpho::tree::min_tree"); @@ -84,11 +144,8 @@ namespace mln mln_precondition(f.is_valid()); mln_precondition(nbh.is_valid()); - typedef p_array<mln_psite(I)> S; - typedef data<I,S> tree_t; - - S s = mln::data::sort_psites_decreasing(f); - tree_t tree(f, s, nbh); + util::ctree::ctree<I> tree = + internal::min_tree_dispatch(mln_trait_image_quant(I)(), f, nbh); trace::exiting("morpho::tree::min_tree"); return tree; @@ -96,7 +153,7 @@ namespace mln template <typename I, typename N> inline - data< I, p_array<mln_psite(I)> > + util::ctree::ctree<I> max_tree(const Image<I>& f_, const Neighborhood<N>& nbh_) { trace::entering("morpho::tree::max_tree"); @@ -107,11 +164,8 @@ namespace mln mln_precondition(f.is_valid()); mln_precondition(nbh.is_valid()); - typedef p_array<mln_psite(I)> S; - typedef data<I,S> tree_t; - - S s = mln::data::sort_psites_increasing(f); - tree_t tree(f, s, nbh); + util::ctree::ctree<I> tree = + internal::max_tree_dispatch(mln_trait_image_quant(I)(), f, nbh); trace::exiting("morpho::tree::max_tree"); return tree; diff --git a/milena/mln/morpho/tree/compute_attribute_image.hh b/milena/mln/morpho/tree/compute_attribute_image.hh index 9126bf1..6e5e891 100644 --- a/milena/mln/morpho/tree/compute_attribute_image.hh +++ b/milena/mln/morpho/tree/compute_attribute_image.hh @@ -23,22 +23,20 @@ // 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_COMPUTE_ATTRIBUTE_IMAGE_HH -# define MLN_MORPHO_TREE_COMPUTE_ATTRIBUTE_IMAGE_HH - -/// \file +#ifndef COMPUTE_ATTRIBUTE_IMAGE_HH +# define COMPUTE_ATTRIBUTE_IMAGE_HH +/// +/// +/// \brief Compute an attribute on a morphological tree. /// -/// Compute a canonized tree from an image. +/// \todo Specialize for fast image. /// -/// \todo Specialize for low quant (and try fastest). -# include <mln/core/routine/duplicate.hh> -# include <mln/core/concept/image.hh> -# include <mln/morpho/tree/data.hh> -# include <mln/trait/accumulators.hh> -# include <mln/util/pix.hh> -# include <mln/data/fill.hh> +# include <mln/core/concept/tree.hh> +# include <mln/core/concept/accumulator.hh> +# include <mln/core/image/attribute_image.hh> +# include <mln/data/fill.hh> namespace mln { @@ -49,213 +47,132 @@ namespace mln namespace tree { - /** - ** Compute an attribute image using tree with a parent - ** relationship between sites. In the attribute image, the - ** resulting value at a node is the 'sum' of its sub-components - ** value + the attribute value at this node. - ** - ** Warning: \p s translates the ordering related to the - ** "natural" childhood relationship. The parenthood is thus - ** inverted w.r.t. to \p s. - ** - ** It is very convenient since all processing upon the parent - ** tree are performed following \p s (in the default "forward" - ** way). - ** - ** FIXME: Put it more clearly... - ** - ** The parent result image verifies: \n - ** - p is root iff parent(p) == p \n - ** - p is a node iff either p is root or f(parent(p)) != f(p). - ** - ** \param[in] a Attribute. - ** \param[in] t Component tree. - ** \param[out] accu_image Optional argument used to store image - ** of attribute accumulator. - ** - ** \return The attribute image. - */ - template <typename A, typename T> - mln_ch_value(typename T::function, mln_result(A)) - compute_attribute_image(const Accumulator<A>& a, - const T& t, - mln_ch_value(typename T::function, A)* accu_image = 0); - - - - /** The same as compute_attribute_image but uses the values - ** stored by \p values image instead. - ** - ** \param[in] a Attribute. - ** \param[in] t Component tree. - ** \param[in] values Value image. - ** \param[out] accu_image Optional argument used to store image. - ** - ** \return - */ - template <typename A, typename T, typename V> - mln_ch_value(typename T::function, mln_result(A)) - compute_attribute_image_from(const Accumulator<A>& a, - const T& t, - const Image<V>& values, - mln_ch_value(typename T::function, A)* accu_image = 0); - + template <typename T, typename A> + attribute_image<T, mln_result(A)> + compute_attribute_image(const Tree<T>& tree, + const Accumulator<A>& accu); # ifndef MLN_INCLUDE_ONLY - // Take_as_init specialization namespace internal { - template <typename A, typename I, typename P> - void take_as_init (trait::accumulator::when_pix::use_none, A& accu, - const I& input, const P& p) - { - (void)input; - (void)p; - accu.take_as_init(); - } - template <typename A, typename I, typename P> - void take_as_init (trait::accumulator::when_pix::use_pix, A& accu, - const I& input, const P& p) + template <typename T, typename A> + inline + void + take_as_init(trait::accumulator::when_pix::use_none, + const T&, A& accu, + const mln_psite(T)&) { - accu.take_as_init(make::pix(input, p)); + accu.take(); } - template <typename A, typename I, typename P> - void take_as_init (trait::accumulator::when_pix::use_v, A& accu, - const I& input, const P& p) + template <typename T, typename A> + inline + void + take_as_init(trait::accumulator::when_pix::use_pix, + const T& tree, A& accu, + const mln_psite(T)& p) { - accu.take_as_init(input(p)); + accu.take(make::pix(p, tree.f(p))); } - template <typename A, typename I, typename P> - void take_as_init (trait::accumulator::when_pix::use_p, A& accu, - const I& input, const P& p) + + template <typename T, typename A> + inline + void + take_as_init(trait::accumulator::when_pix::use_v, + const T& tree, A& accu, + const mln_psite(T)& p) { - (void) input; - accu.take_as_init(p); + accu.take(tree.f(p)); } - - template <typename A, typename I, typename P> - void take_as_init (A& accu, const I& input, const P& p) + template <typename T, typename A> + inline + void + take_as_init(trait::accumulator::when_pix::use_p, + const T&, A& accu, + const mln_psite(T)& p) { - take_as_init (mln_trait_accumulator_when_pix(A)(), accu, input, p); + accu.take(p); } - - template <typename A, typename T, typename V> - inline - mln_ch_value(typename T::function, mln_result(A)) - compute_attribute_image(const A& a, - const T& t, - const V& values, - mln_ch_value(typename T::function, A)* accu_image = 0) + template <typename T, typename A> + void + take_as_init(const T& tree, A& a, const mln_psite(T)& p) { + take_as_init(mln_trait_accumulator_when_pix(A)(), tree, a, p); + } - typedef typename T::function I; - mln_ch_value(I, A) acc; - initialize(acc, t.f()); + } // end of namespace namespace mln::morpho::tree::internal - { - // Transmit "dynamic data" (state) of 'a' to every values of - // 'acc'. It is usually a no-op (so useless) except for a - // few accumulators, e.g., for accu::stat::rank which has the 'k' - // attribute. - mln::data::fill(acc, a); - } + namespace impl + { - { - // Initialize every attribute with the corresponding pixel. - mln_site_piter(T) p(t); - for_all(p) - take_as_init(acc(p), values, p); - } + namespace generic + { + template <typename T, typename A> + attribute_image<T, mln_result(A)> + compute_attribute_image(const T& tree, const A& accu, attribute_image<T, A>* accu_img = 0) { - mln_up_site_piter(T) p(t); - // Propagate attribute from a site to its parent. - for_all(p) - if (! t.is_root(p)) - acc(t.parent(p)).take(acc(p)); - - // Back-propagate attribute from a node to sites of its - // component. Below, p is a non-node component site and - // parent(p) is a node, that is, the site representative of - // the component p belongs to. - for_all(p) - if (! t.is_a_node(p)) - { - mln_assertion(t.is_a_node(t.parent(p))); - acc(p) = acc(t.parent(p)); - } + typedef attribute_image<T, A> I; + attribute_image<T, A> a(tree); + (void) accu_img; + + // Transmit dynamic data of accumulator (ex: rank filter). + data::fill(a, accu); + + // Initialize. Each node takes its content. + { + mln_fwd_piter(T::domain_t) p(tree.domain()); + for_all(p) + internal::take_as_init(tree, a.element(tree.node_at_(p)), p); + } + + // Transmit to parent. + for (unsigned i = tree.n_nodes() - 1; i > 0; --i) + { + unsigned q = tree.parent_at_(i); + if (q != i) // The current node is not a root. + a.element(q).take(a.element(i)); + } + + attribute_image<T, mln_result(A)> out(tree); + for (int i = 0; i < tree.n_nodes(); i++) + out.element(i) = a.element(i).to_result(); + + return out; } + } // end of namespace mln::morpho::tree::impl::generic - // Store accumulator image. - if (accu_image) - *accu_image = duplicate(acc); - - typedef typename T::function I; - mln_ch_value(I, mln_result(A)) output; - initialize(output, acc); - mln::data::fill(output, acc); - - return output; - } - } - - // Facade. - - template <typename A, typename T> - inline - mln_ch_value(typename T::function, mln_result(A)) - compute_attribute_image(const Accumulator<A>& a_, - const T& t, - mln_ch_value(typename T::function, A)* accu_image = 0) - { - trace::entering("morpho::tree::compute_attribute_image"); + } // end of namespace mln::morpho::tree::impl - mln_ch_value(typename T::function, mln_result(A)) output; - output = internal::compute_attribute_image(exact(a_), t, t.f(), - accu_image); - trace::exiting("morpho::tree::compute_attribute_image"); - return (output); - } - - template <typename A, typename T, typename V> - inline - mln_ch_value(typename T::function, mln_result(A)) - compute_attribute_image_from(const Accumulator<A>& a_, - const T& t, - const Image<V>& values, - mln_ch_value(typename T::function, A)* accu_image = 0) + // Facades + template <typename T, typename A> + attribute_image<T, mln_result(A)> + compute_attribute_image(const Tree<T>& tree_, const Accumulator<A>& acc) { - trace::entering("morpho::tree::compute_attribute_image_from"); - + trace::entering("mln::morpho::tree::compute_attribute_image"); + const T& tree = exact(tree_); + const A& accu = exact(acc); - mln_ch_value(typename T::function, mln_result(A)) output; - output = internal::compute_attribute_image(exact(a_), t, exact(values), - accu_image); + mln_precondition(tree.is_valid()); - trace::exiting("morpho::tree::compute_attribute_image_from"); - return output; + return impl::generic::compute_attribute_image(tree, accu); + trace::exiting("mln::morpho::tree::compute_attribute_image"); } - - - # endif // ! MLN_INCLUDE_ONLY - } // end of namespace mln::morpho::tree - - } // end of namespace mln::morpho + } // end of namespace mln::morpho::tree -} // end of namespace mln + } // end of namespace mln::morpho +} // end of namespace mln -#endif // ! MLN_MORPHO_TREE_COMPUTE_ATTRIBUTE_IMAGE_HH +#endif // !COMPUTE_ATTRIBUTE_IMAGE_HH diff --git a/milena/mln/morpho/tree/dual_input_tree.hh b/milena/mln/morpho/tree/dual_input_tree.hh index 08e6553..b4d4839 100644 --- a/milena/mln/morpho/tree/dual_input_tree.hh +++ b/milena/mln/morpho/tree/dual_input_tree.hh @@ -23,17 +23,12 @@ // exception does not however invalidate any other reasons why the // executable file might be covered by the GNU General Public License. -/// \file -/// -/// 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/tag/tree.hh> +# include <mln/util/ctree/ctree.hh> +# include <mln/morpho/tree/impl/union_find_fast.hh> # include <mln/morpho/tree/impl/dual_hqueue.hh> namespace mln @@ -45,67 +40,25 @@ namespace mln namespace tree { - /// Compute the dual input max tree using mask-based connectivity. + /// Compute a canonized max-tree. /// - /// \param[in] f The original image. - /// \param[in] m The connectivity mask. - /// \param[in] nbh The neighborhood of the mask. + /// \param[in] f The input image. + /// \param[in] nbh The neighborhood. /// - /// \return The computed tree. + /// \return The corresponding max-tree structure. /// template <typename I, typename N> inline - data< I, p_array<mln_psite(I)> > + util::ctree::ctree<I> dual_input_max_tree(const Image<I>& f, const Image<I>& m, const Neighborhood<N>& nbh); - # ifndef MLN_INCLUDE_ONLY - 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)> > + util::ctree::ctree<I> dual_input_max_tree(const Image<I>& f_, const Image<I>& m_, const Neighborhood<N>& nbh_) @@ -117,23 +70,16 @@ namespace mln 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; + util::ctree::ctree<I> tree = + morpho::tree::impl::dual_hqueue(tag::tree::max, f, 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"); + trace::exiting("morpho::tree::max_tree"); return tree; } - - - # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::morpho::tree @@ -142,5 +88,4 @@ namespace mln } // end of namespace mln - #endif // !MLN_MORPHO_TREE_DUAL_INPUT_TREE_HH diff --git a/milena/mln/morpho/tree/filter/direct.hh b/milena/mln/morpho/tree/filter/direct.hh index 1835e07..0fdb01c 100644 --- a/milena/mln/morpho/tree/filter/direct.hh +++ b/milena/mln/morpho/tree/filter/direct.hh @@ -27,7 +27,9 @@ # define MLN_MORPHO_TREE_FILTER_DIRECT_HH # include <mln/core/concept/function.hh> -# include <mln/morpho/tree/data.hh> +# include <mln/core/image/attribute_image.hh> +# include <mln/morpho/tree/propagate_node.hh> + /** ** \file mln/morpho/tree/filter/direct.hh @@ -51,36 +53,39 @@ namespace mln /** ** Direct non-pruning strategy. A node is removed if it does - ** not verify the predicate. The sub-components remain intact. + ** not verify the predicate. Its pixels are lowered in gray + ** level to the highest ancestor which meets the criterion, + ** and while its descendants are unaffected. ** - ** \param[in] tree Component tree. - ** \param[out] f_ Image to filter. + ** \param[in,out] attr The attribute image. ** \param[in] pred_ Filtering criterion. */ - template <typename T, typename F, typename P> + template <typename T, typename V, typename P> inline void - direct(const T& tree, Image<F>& f_, const Function_v2b<P>& pred_); + direct(attribute_image<T,V>& attr, const Function_v2b<P>& pred_); # ifndef MLN_INCLUDE_ONLY - template <typename T, typename F, typename P> + template <typename T, typename V, typename F> inline void - direct(const T& tree, Image<F>& f_, const Function_v2b<P>& pred_) + direct(attribute_image<T,V>& attr, const Function_v2b<F>& pred_) { - F& f = exact(f_); - const P& pred = exact(pred_); - trace::entering("mln::morpho::tree::filter::direct"); - mln_dn_node_piter(T) n(tree); - for_all(n) - if (!pred(n)) - f(n) = f(tree.parent(n)); + typedef attribute_image<T,V> A; + + const F& pred = exact(pred_); + const T& tree = attr.tree(); + + mln_fwd_piter(A) node(attr.domain()); + for_all(node) + if (!pred(node)) + attr(node) = attr(tree.parent(node)); trace::exiting("mln::morpho::tree::filter::direct"); } diff --git a/milena/mln/morpho/tree/filter/max.hh b/milena/mln/morpho/tree/filter/max.hh index fb0c580..96c2928 100644 --- a/milena/mln/morpho/tree/filter/max.hh +++ b/milena/mln/morpho/tree/filter/max.hh @@ -27,14 +27,16 @@ # define MLN_MORPHO_TREE_FILTER_MAX_HH # include <mln/core/concept/function.hh> -# include <mln/morpho/tree/data.hh> +# include <mln/core/image/attribute_image.hh> +# include <mln/morpho/tree/propagate_node.hh> /** ** \file mln/morpho/tree/filter/max.hh ** ** \brief Filtering with max strategy. ** -** +** See comments in min filter documentation. +** FIXME: it is possible to make it in one pass ! */ namespace mln @@ -54,51 +56,59 @@ namespace mln ** children are removed or if it does not verify the predicate ** \p pred_. ** - ** \param[in] tree Component tree. - ** \param[out] f_ Image to filter. + ** \param[in,out] attr The attribute image. ** \param[in] pred_ Filtering criterion. */ - template <typename T, typename F, typename P> + template <typename T, typename V, typename P> inline void - max(const T& tree, Image<F>& f_, const Function_v2b<P>& pred_); - + max(attribute_image<T,V>& attr, const Function_v2b<P>& pred_); # ifndef MLN_INCLUDE_ONLY - template <typename T, typename F, typename P> + template <typename T, typename V, typename P> inline void - max(const T& tree, Image<F>& f_, const Function_v2b<P>& pred_) + max(attribute_image<T,V>& attr, const Function_v2b<P>& pred_) { - F& f = exact(f_); - const P& pred = exact(pred_); - trace::entering("mln::morpho::tree::filter::max"); - mln_ch_value(F, bool) mark; - initialize(mark, f); - mln::data::fill(mark, true); + typedef attribute_image<T,V> A; + const P& pred = exact(pred_); + const T& tree = attr.tree(); + + typedef mln_ch_value(A, char) M; + M marker; + initialize(marker, attr); + mln::data::fill(marker, true); + // 1st pass up. { - mln_up_node_piter(T) n(tree); - for_all(n) - if (!mark(n)) - mark(tree.parent(n)) = false; - else if (pred(n)) + mln_bkd_piter(M) node(marker.domain()); + for_all(node) + if (!marker(node)) + marker(tree.parent(node)) = false; + else if (pred(node)) { - mark(tree.parent(n)) = false; - mark(n) = false; + marker(tree.parent(node)) = false; + marker(node) = false; } } + // 2nd pass down: propagation. { - mln_dn_node_piter(T) n(tree); - for_all(n) - if (mark(n)) - f(n) = f(tree.parent(n)); + for (int i = 0; i < tree.n_nodes(); ++i) + { + mln_invariant(tree.has_index(i)); + mln_psite(A) n(tree, i); + if (marker(n)) + { + attr(n) = attr(tree.parent(n)); + i += propagate_node_to_descendants(attr, n, attr(n)); + } + } } trace::exiting("mln::morpho::tree::filter::max"); diff --git a/milena/mln/morpho/tree/filter/min.hh b/milena/mln/morpho/tree/filter/min.hh index db46243..16d2224 100644 --- a/milena/mln/morpho/tree/filter/min.hh +++ b/milena/mln/morpho/tree/filter/min.hh @@ -27,7 +27,8 @@ # define MLN_MORPHO_TREE_FILTER_MIN_HH # include <mln/core/concept/function.hh> -# include <mln/morpho/tree/data.hh> +# include <mln/core/image/attribute_image.hh> +# include <mln/morpho/tree/propagate_node.hh> /** @@ -52,45 +53,57 @@ namespace mln /** - ** Min pruning strategy. A node is removed iif its parent is + ** Min pruning strategy. A node is removed iff its parent is ** removed or if it does not verify the predicate \p pred_. ** - ** \param[in] tree Component tree. - ** \param[out] f_ Image to filter. + ** \param[in,out] attr The attribute image. ** \param[in] pred_ Filtering criterion. */ - template <typename T, typename F, typename P> + template <typename T, typename V, typename F> inline void - min(const T& tree, Image<F>& f_, const Function_v2b<P>& pred_); + min(attribute_image<T,V>& attr, const Function_v2b<F>& pred); # ifndef MLN_INCLUDE_ONLY - template <typename T, typename F, typename P> + template <typename T, typename V, typename F> inline void - min(const T& tree, Image<F>& f_, const Function_v2b<P>& pred_) + min(attribute_image<T,V>& attr, const Function_v2b<F>& pred_) { - F& f = exact(f_); - const P& pred = exact(pred_); - trace::entering("mln::morpho::tree::filter::min"); - mln_ch_value(F, bool) mark; - initialize(mark, f); - mln::data::fill(mark, false); - - mln_dn_node_piter(T) n(tree); - for_all(n) - if (mark(tree.parent(n)) || !pred(n)) - { - f(n) = f(tree.parent(n)); - mark(n) = true; - } - + typedef attribute_image<T,V> A; + + //mlc_is(mln_psite(A), mln_argument(F))::check(); + + const F& pred = exact(pred_); + const T& tree = attr.tree(); + + // Ce qui devrait être ! + // mln_fwd_piter(A) n(attr); + // for_all(n) + // if (!pred(n)) + // { + // data::fill_with_value(attr | tree.desc(n), attr(n)); + // n.go_to_brother(); + // } + + // En attendant que le morpher subimage soit spécialisé, + // on joue avec les indexes. + for (int i = 0; i < tree.n_nodes(); ++i) + { + mln_invariant(tree.has_index(i)); + mln_psite(A) n(tree, i); + if (!pred(n)) + { + attr(n) = attr(tree.parent(n)); + i += propagate_node_to_descendants(attr, n, attr(n)); + } + } trace::exiting("mln::morpho::tree::filter::min"); } diff --git a/milena/mln/morpho/tree/filter/subtractive.hh b/milena/mln/morpho/tree/filter/subtractive.hh index 10d3043..b905f15 100644 --- a/milena/mln/morpho/tree/filter/subtractive.hh +++ b/milena/mln/morpho/tree/filter/subtractive.hh @@ -27,10 +27,8 @@ # define MLN_MORPHO_TREE_FILTER_SUBTRACTIVE_HH # include <mln/core/concept/function.hh> -# include <mln/fun/ops.hh> - -# include <mln/morpho/tree/data.hh> -# include <mln/morpho/tree/propagate_if.hh> +# include <mln/core/image/attribute_image.hh> +# include <mln/morpho/tree/propagate_node.hh> /** ** \file mln/morpho/tree/filter/subtractive.hh @@ -57,45 +55,61 @@ namespace mln ** does not verify the predicate. The sub-components values ** are set to the value of the removed component. ** - ** \param[in] tree Component tree. - ** \param[out] f_ Image to filter. + ** \param[in] ori The values as an attribute image. ** \param[in] pred_ Filtering criterion. + ** \return A copy of \p ori filtered with the subtractive strategy. */ - template <typename T, typename F, typename P> + template <typename T, typename V, typename P> inline - void - subtractive(const T& tree, Image<F>& f_, const Function_v2b<P>& pred_); - - + attribute_image<T, V> + subtractive(const attribute_image<T,V>& ori, const Function_v2b<P>& pred_); # ifndef MLN_INCLUDE_ONLY - template <typename T, typename F, typename P> + template <typename T, typename V, typename F> inline - void - subtractive(const T& tree, Image<F>& f_, const Function_v2b<P>& pred_) + attribute_image<T, V> + subtractive(const attribute_image<T,V>& ori, const Function_v2b<F>& pred_) { - F& f = exact(f_); - const P& pred = exact(pred_); - trace::entering("mln::morpho::tree::filter::subtractive"); - morpho::tree::propagate_if(tree, f, morpho::tree::desc_propagation (), !pred); + typedef attribute_image<T,V> A; + + const F& pred = exact(pred_); + const T& tree = ori.tree(); + + mln_precondition(ori.is_valid()); - mln_up_node_piter(T) n(tree); + A f; + initialize(f, ori); + + // Fixme: introduce mln_root_piter... + // and traverse all roots, not just the first one. + mln_psite(A) root = util::ctree::node<T>(tree, 0); + f(root) = ori(root); + std::cout << "Root value: " << ori(root) << std::endl; + + mln_piter(A) n(ori.domain()); for_all(n) if (!pred(n)) f(n) = f(tree.parent(n)); + else + f(n) = f(tree.parent(n)) + ori(n) - ori(tree.parent(n)); trace::exiting("mln::morpho::tree::filter::subtractive"); + + return f; } # endif // ! MLN_INCLUDE_ONLY } // end of namespace mln::morpho::tree::filter + } // end of namespace mln::morpho::tree + } // end of namespace mln::morpho + } // end of namespace mln #endif // ! MLN_MORPHO_TREE_FILTER_SUBTRACTIVE_HH diff --git a/milena/mln/morpho/tree/impl/dual_hqueue.hh b/milena/mln/morpho/tree/impl/dual_hqueue.hh index 41334f0..f00642f 100644 --- a/milena/mln/morpho/tree/impl/dual_hqueue.hh +++ b/milena/mln/morpho/tree/impl/dual_hqueue.hh @@ -32,31 +32,19 @@ /// 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/util/ctree/ctree.hh> # include <mln/data/fill.hh> +# include <mln/extension/fill.hh> +# include <mln/extension/adjust.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 @@ -79,320 +67,430 @@ namespace mln /// /// \return The tree structure. /// - template <typename I, typename N> - inline - data<I, p_array<mln_psite(I)> > - dual_hqueue(const Image<I>& f, + template <typename T, typename I, typename N> + util::ctree::ctree<I> + dual_hqueue(const tag::tree::tree_t<T>&, + 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 + namespace impl { - template <typename I, typename N, class E> - struct shared_flood_args + + template <typename I, typename N, typename E> + struct flooder_base_ : Object<E> { - typedef mln_psite(I) P; + typedef unsigned P; typedef mln_value(I) V; - typedef p_array<P> S; + // Constructor. + flooder_base_(const I& f_, const I& m_, const N& nbh_); + + // Returned values. + mln_ch_value(I, unsigned) map; + std::vector<P> location; + std::vector<V> values; + std::vector<P> parent; + std::vector<unsigned> area; + unsigned nnode; + + + protected: + enum { + nvalues = mln_card(mln_value(I)) + }; + + + // Recursive function that floods. Canevas function. + int flood(int h_idx); + + // HLevel comparaison. To be specialized for min/max tree. + bool hlevel_less(int h_idx1, int h_idx2); + + // Parent attachment function. To be specialized for min/max tree. + int attach_parent(int h_idx); + + // Data. + // \{ 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. + // Auxiliary data. + util::hqueues<unsigned, V> hqueues; + value::set<V> vset; + util::array<int> dp; + unsigned n_nbhs; + + mln_ch_value(I, bool) deja_vu; 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); - } + unsigned loc; + // \} }; - 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) + template <typename T, typename I, typename N> + struct flooder; + + template <typename I, typename N> + struct flooder<tag::tree::max_t, I, N> : + flooder_base_< I, N, flooder<tag::tree::max_t, I, N> > { - 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]; - } + typedef mln_value(I) V; + typedef flooder_base_< I, N, flooder<tag::tree::max_t, I, N> > super_; - // Merge histograms. - for (unsigned i = 0; i < hm.nvalues(); ++i) - hm[i] += hf[i]; + // Constructor. + flooder(const I& f, const I& m, const N& nbh); - // Retrieve pmin. - mln_piter(I) p(m.domain()); - for (p.start(); m(p) != hmin; p.next()) - ; - mln_assertion(p.is_valid()); - pmin = p; + // HLevel comparaison. + bool hlevel_less_(int h_idx1, int h_idx2); - return hm; - } + // Parent attachment function. + int attach_parent_(int h_idx); + + public: + using super_::parent; + using super_::area; + using super_::location; - // 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 + protected: + using super_::is_node_at_level; + using super_::node_at_level; + using super_::vset; + using super_::loc; + using super_::hqueues; + }; + + template <typename I, typename N> + struct flooder<tag::tree::min_t, I, N> : + flooder_base_< I, N, flooder<tag::tree::min_t, I, N> > { - extend(const mln_psite(I)::delta& dp) - : dp_ (dp) - { - } + typedef mln_value(I) V; + typedef flooder_base_< I, N, flooder<tag::tree::min_t, I, N> > super_; - mln_psite(I) operator() (const mln_psite(I)& p) const - { - return p + dp_; - } + // Constructor. + flooder(const I& f, const I& m, const N& nbh); - private: - const mln_psite(I)::delta dp_; + // HLevel comparaison. + bool hlevel_less_(int h_idx1, int h_idx2); + + // Parent attachment function. + unsigned attach_parent_(int h_idx); + + public: + using super_::parent; + using super_::area; + using super_::location; + + protected: + using super_::is_node_at_level; + using super_::node_at_level; + using super_::vset; + using super_::loc; + using super_::hqueues; }; - } // end of namespace mln::morpho::tree::internal - namespace impl - { + template <typename I, typename N, typename E> + inline + bool + flooder_base_<I, N, E>::hlevel_less(int h_idx1, int h_idx2) + { + E& obj = exact(*this); + return obj.hlevel_less_(h_idx1, h_idx2); + } template <typename I, typename N, typename E> - unsigned - flood(internal::shared_flood_args<I, N, E>& args, const unsigned h_idx) + inline + int + flooder_base_<I, N, E>::attach_parent(int h_idx) { - mln_assertion(args.is_node_at_level[h_idx]); + E& obj = exact(*this); + return obj.attach_parent_(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; - } - } - } + template <typename I, typename N> + inline + flooder<tag::tree::max_t, I, N>::flooder(const I& f_, const I& m_, const N& nbh_) + : flooder_base_<I, N, flooder<tag::tree::max_t, I, N> >(f_, m_, nbh_) + { + } - 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]; - } + template <typename I, typename N> + inline + flooder<tag::tree::min_t, I, N>::flooder(const I& f_, const I& m_, const N& nbh_) + : flooder_base_<I, N, flooder<tag::tree::min_t, I, N> >(f_, m_, nbh_) + { + } + template <typename I, typename N> + inline + bool + flooder<tag::tree::max_t, I, N>::hlevel_less_(int h_idx1, int h_idx2) + { + return h_idx1 < h_idx2; + } - // 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); - } - } + template <typename I, typename N> + inline + bool + flooder<tag::tree::min_t, I, N>::hlevel_less_(int h_idx1, int h_idx2) + { + return h_idx2 < h_idx1; + } - // Retrieve dad. - args.is_node_at_level[h_idx] = false; - unsigned c = h_idx; - while (c > 0 && !args.is_node_at_level[c]) + template <typename I, typename N> + inline + int + flooder<tag::tree::max_t, I, N>::attach_parent_(int h_idx) + { + mln_assertion(h_idx >= 0); + mln_assertion(is_node_at_level[h_idx]); + + + unsigned x = node_at_level[h_idx]; + int c = h_idx - 1; + + is_node_at_level[h_idx] = false; + while (c >= 0 && !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; + if (c >= 0) + { + unsigned p = node_at_level[c]; + parent[x] = p; + area[p] += area[x] + 1; + } + else // root case. + parent[x] = x; + location[x] = loc++; 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"); + // template <typename I, typename N> + // inline + // unsigned + // flooder<tag::tree::min_t, I, N>::attach_parent_(unsigned h_idx) + // { + // unsigned nv = vset.nvalues() - 1; + // unsigned c = h_idx; - const I& f = exact(f_); - const I& m = exact(m_); - const N& nbh = exact(neibh_); + // while (c < nv && !is_node_at_level[c]) + // ++c; - typedef mln_psite(I) P; - typedef p_array<mln_psite(I)> S; + // unsigned x = node_at_level[h_idx]; + // if (is_node_at_level[c]) + // { + // unsigned p = node_at_level[c]; + // parent[x] = p; + // area[p] += area[x] + 1; + // } + // else + // parent[x] = x; - util::timer tm; - tm.start(); + // location[x] = loc++; + // return c; + // } - // 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(); + template <typename I, typename N, typename E> + flooder_base_<I, N, E>::flooder_base_(const I& f_, const I& m_, const N& nbh_) + : f (f_), m(m_), nbh (nbh_), is_node_at_level (false) + { + typedef mln_value(I) V; - // Extend the domain. - dp[0] = d.pmax()[0] - d.pmin()[0] + 1; - d.pmax() += dp; - d_ext.pmin() += dp; - d_ext.pmax() += dp; + unsigned pmin_offset; + V vmin; - // 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 + // Allocate hierarchical queues and retrieve the starting point { - if (args.is_node_at_level[i]) - { - parent(args.node_at_level[r]) = args.node_at_level[i]; - r = i; - } + util::array< unsigned > h(nvalues, 0); + + pmin_offset = m.index_of_point(m.domain().pmin()); + vmin = m.element(pmin_offset); + + mln_fwd_pixter(const I) pxl(m); + for_all(pxl) + { + ++h[pxl.val()]; + if (hlevel_less(pxl.val(), vmin)) + { + vmin = pxl.val(); + pmin_offset = pxl.offset(); + } + } + + for (unsigned i = 0; i < nvalues; ++i) + hqueues[i].reserve(h[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); + // Initialize aux data + { + dp = offsets_wrt(f, nbh); + n_nbhs = dp.nelements(); + vset = value::set<V>::the(); + + extension::adjust(f, nbh); + extension::adjust(m, nbh); + + initialize(deja_vu, f); + initialize(map, f); + data::fill(deja_vu, false); + extension::fill(deja_vu, true); + + unsigned nsites = geom::nsites(f); + parent.resize(2 * nsites); + location.resize(2 * nsites); + area.resize(2 * nsites, 0); + values.resize(2 * nsites); + + nnode = 0; + loc = 0; + } + + // Start flooding + unsigned hmin = vset.index_of(vmin); + unsigned fmin = vset.index_of(f.element(pmin_offset)); + std::cout << vmin << std::endl; + hqueues[hmin].push(pmin_offset); + deja_vu.element(pmin_offset) = true; + values[nnode] = vmin; + node_at_level[fmin] = nnode++; + is_node_at_level[fmin] = true; + if (hmin != fmin) + { + values[nnode] = f.element(pmin_offset); + node_at_level[hmin] = nnode++; + is_node_at_level[hmin] = true; + } + int c = flood(hmin); + while (c != -1) // There still are nodes under minimal mask level. + c = attach_parent(c); + mln_assertion(c == -1); + } - p_array<mln_psite(I)> s_f = mln::data::sort_psites_increasing(f); - mln_fwd_piter(S) p(s_f); // Forward. - for_all(p) + template <typename I, typename N, typename E> + int + flooder_base_<I, N, E>::flood(int h_idx) + { + mln_assertion(h_idx >= 0); + mln_assertion(is_node_at_level[h_idx]); + + while (!hqueues[h_idx].empty()) { - P x = p; - P q = parent(p); + unsigned p = hqueues[h_idx].pop_front(); + int p_idx = vset.index_of(f.element(p)); + mln_assertion(vset.index_of(m.element(p)) == (unsigned)h_idx); + //mln_assertion(is_node_at_level[p_idx]); - // 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))) + if (h_idx < p_idx) // FIXME: max tree only { - s.append(q); - deja_vu(q) = true; - x = q; - q = parent(q); + // Singleton + int par_i = node_at_level[h_idx]; + int self_i = nnode++; // Create the node. + map.element(p) = self_i; + values[self_i] = m.element(p); + area[self_i] = 1; + parent[self_i] = par_i; + location[self_i] = loc++; + area[par_i]++; } + else + map.element(p) = node_at_level[p_idx]; - if (d_ext.has(q) && fext(q) == fext(parent(q)) && q != parent(q)) + // Process the neighbors + for (unsigned j = 0; j < n_nbhs; ++j) { - q = (parent(x) = parent(q)); - mln_assertion(f.domain().has(q)); + unsigned n = p + dp[j]; + + if (!deja_vu.element(n)) + { + unsigned n_hidx = vset.index_of(m.element(n)); + unsigned n_fidx = vset.index_of(f.element(n)); + hqueues[n_hidx].push(n); + deja_vu.element(n) = true; + + if (!is_node_at_level[n_hidx]) + { + values[nnode] = (m.element(n)); + node_at_level[n_hidx] = nnode++; + is_node_at_level[n_hidx] = true; + } + if (!is_node_at_level[n_fidx] && n_fidx < n_hidx) + { + values[nnode] = (f.element(n)); + node_at_level[n_fidx] = nnode++; + is_node_at_level[n_fidx] = true; + } + + while (hlevel_less(h_idx, n_hidx)) + n_hidx = flood(n_hidx); + } } + } - if (fext(q) == fext(parent(q))) - parent(x) = parent(q); + return attach_parent(h_idx); + } + + template <typename T, typename I, typename N> + inline + util::ctree::ctree<I> + dual_hqueue(const tag::tree::tree_t<T>&, + const Image<I>& f_, + const Image<I>& m_, + const Neighborhood<N>& nbh_) + { + trace::entering("mln::morpho::tree::impl::hqueue"); + const I& f = exact(f_); + const I& m = exact(m_); + const N& nbh = exact(nbh_); - s.append(p); + mln_precondition(f.is_valid()); + mln_precondition(m.is_valid()); + mln_precondition(nbh.is_valid()); - mln_assertion((q = parent(p)) == parent(q) || fext(q) != fext(parent(q))); - } + // Process + impl::flooder<T, I, N> args(f, m, nbh); - } + // Reserve + util::ctree::ctree<I> tree; + tree.reserve(f, args.nnode); - std::cout << "Construction de l'arbre en " << tm << " s." << std::endl; + // Set values + { + mln_piter(I) p(f.domain()); + for_all(p) + { + unsigned idx = args.map(p); + unsigned loc = args.nnode - args.location[idx] - 1; + tree.node_at_(p) = loc; + } - data<I, S> tree(fext, parent, s); + for (int idx = 0; idx < args.nnode; ++idx) + { + unsigned loc = args.nnode - args.location[idx] - 1; - trace::exiting("mln::morpho::tree::impl::dual_hqueue"); + tree.parent_at_(loc) = args.nnode - args.location[args.parent[idx]] - 1; + tree.f_at_(loc) = args.values[idx]; + tree.length_at_(loc) = args.area[idx]; + } + } + trace::exiting("mln::morpho::tree::impl::hqueue"); return tree; } - } // end of namespace mln::morpho::tree::impl + } // end of namespace mln::morpho::tree::impl + # endif // ! MLN_INCLUDE_ONLY diff --git a/milena/mln/morpho/tree/impl/hqueue.hh b/milena/mln/morpho/tree/impl/hqueue.hh new file mode 100644 index 0000000..cfe0f4b --- /dev/null +++ b/milena/mln/morpho/tree/impl/hqueue.hh @@ -0,0 +1,428 @@ +// 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_HQUEUE_HH +# define MLN_MORPHO_TREE_IMPL_HQUEUE_HH + +# include <mln/util/ctree/ctree.hh> +# include <mln/util/hqueues.hh> +# include <mln/value/value_array.hh> +# include <mln/geom/nsites.hh> + +# include <mln/histo/array.hh> +# include <mln/histo/compute.hh> +# include <mln/data/fill.hh> + +namespace mln +{ + + namespace morpho + { + + namespace tree + { + + namespace impl + { + + template <typename T, typename I, typename N> + inline + util::ctree::ctree<I> + hqueue(const tag::tree::tree_t<T>&, const Image<I>& f, const Neighborhood<N>& nbh); + + } // end of namespace mln::morpho::tree::impl + +# ifndef MLN_INCLUDE_ONLY + + namespace internal + { + + template <typename T, typename I> + struct op_psite_less; + + template <typename I> + struct op_psite_less< tag::tree::max_t, I> + { + op_psite_less(const I& f) : f_(f) {} + + bool operator() (const mln_psite(I)& lhs, const mln_psite(I)& rhs) const + { + return util::ord_strict(f_(lhs), f_(rhs)); + } + + private: + const I& f_; + }; + + template <typename I> + struct op_psite_less< tag::tree::min_t, I> + { + op_psite_less(const I& f) : f_(f) {} + + bool operator() (const mln_psite(I)& lhs, const mln_psite(I)& rhs) const + { + return util::ord_strict(f_(rhs), f_(lhs)); + } + + private: + const I& f_; + }; + + } + + + namespace impl + { + + + template <typename I, typename N, typename E> + struct flooder_base_ : Object<E> + { + typedef unsigned P; + typedef mln_value(I) V; + + // Constructor. + flooder_base_(const I& f_, const N& nbh_); + + // Returned values. + mln_ch_value(I, unsigned) map; + std::vector<P> location; + std::vector<P> parent; + std::vector<unsigned> area; + unsigned n_node; + + + protected: + // Recursive function that floods. Canevas function. + unsigned flood(const unsigned h_idx); + + // HLevel comparaison. To be specialized for min/max tree. + bool hlevel_less(unsigned h_idx1, unsigned h_idx2); + + // Parent attachment function. To be specialized for min/max tree. + unsigned attach_parent(unsigned h_idx); + + // Data. + // \{ + const I& f; + const N& nbh; + + // Auxiliary data. + util::hqueues<mln_psite(I), V> hqueues; + value::set<V> vset; + + mln_ch_value(I, bool) deja_vu; + value::value_array<V, bool> is_node_at_level; + value::value_array<V, P> node_at_level; + unsigned loc; + // \} + }; + + template <typename T, typename I, typename N> + struct flooder; + + template <typename I, typename N> + struct flooder<tag::tree::max_t, I, N> : + flooder_base_< I, N, flooder<tag::tree::max_t, I, N> > + { + typedef mln_value(I) V; + typedef flooder_base_< I, N, flooder<tag::tree::max_t, I, N> > super_; + + // Constructor. + flooder(const I& f_, const N& nbh_); + + // HLevel comparaison. + bool hlevel_less_(unsigned h_idx1, unsigned h_idx2); + + // Parent attachment function. + unsigned attach_parent_(unsigned h_idx); + + public: + using super_::parent; + using super_::area; + using super_::location; + + protected: + using super_::is_node_at_level; + using super_::node_at_level; + using super_::vset; + using super_::loc; + }; + + template <typename I, typename N> + struct flooder<tag::tree::min_t, I, N> : + flooder_base_< I, N, flooder<tag::tree::min_t, I, N> > + { + typedef mln_value(I) V; + typedef flooder_base_< I, N, flooder<tag::tree::min_t, I, N> > super_; + + // Constructor. + flooder(const I& f_, const N& nbh_); + + // HLevel comparaison. + bool hlevel_less_(unsigned h_idx1, unsigned h_idx2); + + // Parent attachment function. + unsigned attach_parent_(unsigned h_idx); + + public: + using super_::parent; + using super_::area; + using super_::location; + + protected: + using super_::is_node_at_level; + using super_::node_at_level; + using super_::vset; + using super_::loc; + }; + + + template <typename I, typename N, typename E> + inline + bool + flooder_base_<I, N, E>::hlevel_less(unsigned h_idx1, unsigned h_idx2) + { + E& obj = exact(*this); + return obj.hlevel_less_(h_idx1, h_idx2); + } + + template <typename I, typename N, typename E> + inline + unsigned + flooder_base_<I, N, E>::attach_parent(unsigned h_idx) + { + E& obj = exact(*this); + return obj.attach_parent_(h_idx); + } + + template <typename I, typename N> + inline + flooder<tag::tree::max_t, I, N>::flooder(const I& f_, const N& nbh_, + const histo::array< mln_value(I) >& histo_) + : flooder_base_<I, N, flooder<tag::tree::max_t, I, N> >(f_, nbh_, histo_) + { + } + + template <typename I, typename N> + inline + flooder<tag::tree::min_t, I, N>::flooder(const I& f_, const N& nbh_, + const histo::array<mln_value(I)>& histo_) + : flooder_base_<I, N, flooder<tag::tree::min_t, I, N> >(f_, nbh_, histo_) + { + } + + template <typename I, typename N> + inline + bool + flooder<tag::tree::max_t, I, N>::hlevel_less_(unsigned h_idx1, unsigned h_idx2) + { + return h_idx1 < h_idx2; + } + + template <typename I, typename N> + inline + bool + flooder<tag::tree::min_t, I, N>::hlevel_less_(unsigned h_idx1, unsigned h_idx2) + { + return h_idx2 < h_idx1; + } + + template <typename I, typename N> + inline + unsigned + flooder<tag::tree::max_t, I, N>::attach_parent_(unsigned h_idx) + { + unsigned c = h_idx; + while (c > 0 && !is_node_at_level[c]) + --c; + + unsigned x = node_at_level[h_idx]; + if (is_node_at_level[c]) // || c > 0 + { + unsigned p = node_at_level[c]; + parent[x] = p; + area[p] += area[x] + 1; + } + else + parent[x] = x; + + location[x] = loc++; + return c; + } + + template <typename I, typename N> + inline + unsigned + flooder<tag::tree::min_t, I, N>::attach_parent_(unsigned h_idx) + { + unsigned nv = vset.nvalues() - 1; + unsigned c = h_idx; + + while (c < nv && !is_node_at_level[c]) + ++c; + + unsigned x = node_at_level[h_idx]; + if (is_node_at_level[c]) // || c < nv + { + unsigned p = node_at_level[c]; + parent[x] = p; + area[p] += area[x] + 1; + } + else + parent[x] = x; + + location[x] = loc++; + return c; + } + + + template <typename I, typename N, typename E> + flooder_base_<I, N, E>::flooder_base_(const I& f_, const N& nbh_, + const histo::array< mln_value(I) >& histo_) + : f (f_), nbh (nbh_), is_node_at_level (false) + { + // Initialization + { + unsigned nsites = geom::nsites(f); + initialize(map, f); + initialize(deja_vu, f); + + parent.resize(nsites); + location.resize(nsites); + area.resize(nsites, 0); + data::fill(deja_vu, false); + + vset = value::set< mln_value(I) >::the(); + + n_node = 0; + loc = 0; + } + + // Get start value + mln_psite(I) pstart = f.domain().pmin(); + mln_value(I) hstart = f(pstart); + { + mln_piter(I) p(f.domain()); + for_all(p) + { + mln_value(I) v = f(p); + ++h[v]; + if (hlevel_less(v, hstart)) + { + hstart = v; + pstart = p; + } + } + + for (unsigned i = 0; i < nvalues; ++i) + hqueues[i].reserve(h[i]); + } + + // Start flooding + unsigned r = vset.index_of(hstart); + deja_vu(pstart) = true; + hqueues[r].push(pstart); + node_at_level[r] = n_node++; + is_node_at_level[r] = true; + + flood(r); + } + + template <typename I, typename N, typename E> + unsigned + flooder_base_<I, N, E>::flood(const unsigned h_idx) + { + mln_assertion(is_node_at_level[h_idx]); + + while (!hqueues[h_idx].empty()) + { + mln_psite(I) p = hqueues[h_idx].pop_front(); + mln_assertion(vset.index_of(f(p)) == h_idx); + + map(p) = node_at_level[h_idx]; + + // Process the neighbors + mln_niter(N) n(nbh, p); + for_all(n) + if (f.domain().has(n) && !deja_vu(n)) + { + unsigned n_idx = vset.index_of(f(n)); + hqueues[n_idx].push(n); + deja_vu(n) = true; + + if (!is_node_at_level[n_idx]) + { + node_at_level[n_idx] = n_node++; + is_node_at_level[n_idx] = true; + } + + while (hlevel_less(h_idx, n_idx)) + n_idx = flood(n_idx); + } + } + + is_node_at_level[h_idx] = false; + return attach_parent(h_idx); + } + + template <typename T, typename I, typename N> + inline + util::ctree::ctree<I> + hqueue(const tag::tree::tree_t<T>&, const Image<I>& f_, const Neighborhood<N>& nbh_) + { + trace::entering("mln::morpho::tree::impl::hqueue"); + const I& f = exact(f_); + const N& nbh = exact(nbh_); + + mln_precondition(f.is_valid()); + mln_precondition(nbh.is_valid()); + + // Process + impl::flooder<T, I, N> args(f, nbh); + // \todo handles here discontinous domains: so multi-roots trees. + + args.parent.resize(args.n_node); + args.area.resize(args.n_node); + args.location.resize(args.n_node); + + util::ctree::ctree<I> tree(f, args.map, args.location, args.parent, + args.area); + + trace::exiting("mln::morpho::tree::impl::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_IMPL_HQUEUE_HH diff --git a/milena/mln/morpho/tree/impl/hqueue_fast.hh b/milena/mln/morpho/tree/impl/hqueue_fast.hh new file mode 100644 index 0000000..645bbc4 --- /dev/null +++ b/milena/mln/morpho/tree/impl/hqueue_fast.hh @@ -0,0 +1,391 @@ +#ifndef MLN_MORPHO_TREE_IMPL_HQUEUE_FASTHH +# define MLN_MORPHO_TREE_IMPL_HQUEUE_FASTHH + +# include <mln/extension/adjust.hh> +# include <mln/extension/fill.hh> +# include <mln/util/ctree/ctree.hh> +# include <mln/util/hqueues.hh> +# include <mln/value/value_array.hh> +# include <mln/geom/nsites.hh> +# include <mln/data/fill.hh> + +namespace mln +{ + + namespace morpho + { + + namespace tree + { + + namespace impl + { + + template <typename T, typename I, typename N> + util::ctree::ctree<I> + hqueue_fast(const tag::tree::tree_t<T>&, const Image<I>& f, const Neighborhood<N>& nbh); + + } // end of namespace mln::morpho::tree::impl + +# ifndef MLN_INCLUDE_ONLY + + namespace impl + { + + template <typename I, typename N, typename E> + struct flooder_base_ : Object<E> + { + typedef unsigned P; + typedef mln_value(I) V; + + // Constructor. + flooder_base_(const I& f_, const N& nbh_); + + mln_ch_value(I, unsigned) map; + std::vector<P> location; + std::vector<P> parent; + std::vector<unsigned> area; + unsigned n_node; + + protected: + enum { + nvalues = mln_card(mln_value(I)) + }; + + // Recursive function that floods. Canevas function. + unsigned flood(const unsigned h_idx); + + // HLevel comparaison. To be specialized for min/max tree. + bool hlevel_less(unsigned h_idx1, unsigned h_idx2); + + // Parent attachment function. To be specialized for min/max tree. + unsigned attach_parent(unsigned h_idx); + + + // Data. + // \{ + const I& f; + const N& nbh; + + // Auxiliary data. + util::hqueues<unsigned, V> hqueues; + value::set<V> vset; + util::array<int> dp; + unsigned n_nbhs; + + mln_ch_value(I, bool) deja_vu; + value::value_array<V, bool> is_node_at_level; + value::value_array<V, P> node_at_level; + unsigned loc; + // \} + }; + + template <typename T, typename I, typename N> + struct flooder; + + template <typename I, typename N> + struct flooder<tag::tree::max_t, I, N> : + flooder_base_< I, N, flooder<tag::tree::max_t, I, N> > + { + typedef mln_value(I) V; + typedef flooder_base_< I, N, flooder<tag::tree::max_t, I, N> > super_; + + // Constructor. + flooder(const I& f_, const N& nbh_); + + // HLevel comparaison. + bool hlevel_less_(unsigned h_idx1, unsigned h_idx2); + + // Parent attachment function. + unsigned attach_parent_(unsigned h_idx); + + public: + using super_::parent; + using super_::area; + using super_::location; + + protected: + using super_::is_node_at_level; + using super_::node_at_level; + using super_::vset; + using super_::loc; + }; + + template <typename I, typename N> + struct flooder<tag::tree::min_t, I, N> : + flooder_base_< I, N, flooder<tag::tree::min_t, I, N> > + { + typedef mln_value(I) V; + typedef flooder_base_< I, N, flooder<tag::tree::min_t, I, N> > super_; + + // Constructor. + flooder(const I& f_, const N& nbh_); + + // HLevel comparaison. + bool hlevel_less_(unsigned h_idx1, unsigned h_idx2); + + // Parent attachment function. + unsigned attach_parent_(unsigned h_idx); + + public: + using super_::parent; + using super_::area; + using super_::location; + + protected: + using super_::is_node_at_level; + using super_::node_at_level; + using super_::vset; + using super_::loc; + }; + + template <typename I, typename N, typename E> + inline + bool + flooder_base_<I, N, E>::hlevel_less(unsigned h_idx1, unsigned h_idx2) + { + E& obj = exact(*this); + return obj.hlevel_less_(h_idx1, h_idx2); + } + + template <typename I, typename N, typename E> + inline + unsigned + flooder_base_<I, N, E>::attach_parent(unsigned h_idx) + { + E& obj = exact(*this); + return obj.attach_parent_(h_idx); + } + + template <typename I, typename N> + inline + flooder<tag::tree::max_t, I, N>::flooder(const I& f_, const N& nbh_) + : flooder_base_<I, N, flooder<tag::tree::max_t, I, N> >(f_, nbh_) + { + } + + template <typename I, typename N> + inline + flooder<tag::tree::min_t, I, N>::flooder(const I& f_, const N& nbh_) + : flooder_base_<I, N, flooder<tag::tree::min_t, I, N> >(f_, nbh_) + { + } + + template <typename I, typename N> + inline + bool + flooder<tag::tree::max_t, I, N>::hlevel_less_(unsigned h_idx1, unsigned h_idx2) + { + return h_idx1 < h_idx2; + } + + template <typename I, typename N> + inline + bool + flooder<tag::tree::min_t, I, N>::hlevel_less_(unsigned h_idx1, unsigned h_idx2) + { + return h_idx2 < h_idx1; + } + + template <typename I, typename N> + inline + unsigned + flooder<tag::tree::max_t, I, N>::attach_parent_(unsigned h_idx) + { + unsigned c = h_idx; + while (c > 0 && !is_node_at_level[c]) + --c; + + unsigned x = node_at_level[h_idx]; + if (is_node_at_level[c]) // || c > 0 + { + unsigned p = node_at_level[c]; + parent[x] = p; + area[p] += area[x] + 1; + } + else + parent[x] = x; + + location[x] = loc++; + return c; + } + + template <typename I, typename N> + inline + unsigned + flooder<tag::tree::min_t, I, N>::attach_parent_(unsigned h_idx) + { + unsigned nv = vset.nvalues() - 1; + unsigned c = h_idx; + + while (c < nv && !is_node_at_level[c]) + ++c; + + unsigned x = node_at_level[h_idx]; + if (is_node_at_level[c]) // || c < nv + { + unsigned p = node_at_level[c]; + parent[x] = p; + area[p] += area[x] + 1; + } + else + parent[x] = x; + + location[x] = loc++; + return c; + } + + + template <typename I, typename N, typename E> + flooder_base_<I, N, E>::flooder_base_(const I& f_, const N& nbh_) + : f (f_), nbh (nbh_), is_node_at_level(false) + { + // Prepare hierarchical queues and retrieve start point. + unsigned pmin; + mln_value(I) vmin; + { + util::array< unsigned > h(nvalues, 0); + + pmin = f.index_of_point(f.domain().pmin()); + vmin = f.element(pmin); + + mln_fwd_pixter(const I) pxl(f); + for_all(pxl) + { + ++h[pxl.val()]; + if (hlevel_less(pxl.val(), vmin)) + { + vmin = pxl.val(); + pmin = pxl.offset(); + } + } + + for (unsigned i = 0; i < nvalues; ++i) + hqueues[i].reserve(h[i]); + } + + // Initialize aux data + { + dp = offsets_wrt(f, nbh); + n_nbhs = dp.nelements(); + vset = value::set< mln_value(I) >::the(); + + extension::adjust(f, nbh); + + initialize(deja_vu, f); + initialize(map, f); + data::fill(deja_vu, false); + extension::fill(deja_vu, true); + + unsigned nsites = geom::nsites(f); + parent.resize(nsites); + location.resize(nsites); + area.resize(nsites, 0); + + n_node = 0; + loc = 0; + } + + // Start flooding + unsigned hmin = vset.index_of(vmin); + deja_vu.element(pmin) = true; + hqueues[hmin].push(pmin); + node_at_level[hmin] = n_node++; + is_node_at_level[hmin] = true; + + flood(hmin); + } + + template <typename I, typename N, typename E> + unsigned + flooder_base_<I, N, E>::flood(const unsigned h_idx) + { + mln_assertion(is_node_at_level[h_idx]); + + while (!hqueues[h_idx].empty()) + { + unsigned p = hqueues[h_idx].pop_front(); + mln_assertion(vset.index_of(f.element(p)) == h_idx); + + map.element(p) = node_at_level[h_idx]; + + // Process the neighbors + for (unsigned j = 0; j < n_nbhs; ++j) + { + unsigned n = p + dp[j]; + + if (!deja_vu.element(n)) + { + unsigned n_idx = vset.index_of(f.element(n)); + hqueues[n_idx].push(n); + deja_vu.element(n) = true; + + if (!is_node_at_level[n_idx]) + { + node_at_level[n_idx] = n_node++; + is_node_at_level[n_idx] = true; + } + + while (hlevel_less(h_idx, n_idx)) + n_idx = flood(n_idx); + } + } + } + + is_node_at_level[h_idx] = false; + return attach_parent(h_idx); + } + + template <typename T, typename I, typename N> + inline + util::ctree::ctree<I> + hqueue_fast(const tag::tree::tree_t<T>&, const Image<I>& f_, const Neighborhood<N>& nbh_) + { + trace::entering("mln::morpho::tree::impl::hqueue_fast"); + const I& f = exact(f_); + const N& nbh = exact(nbh_); + + mln_precondition(f.is_valid()); + mln_precondition(nbh.is_valid()); + + // Process + impl::flooder<T, I, N> args(f, nbh); + // \todo handles here discontinous domains: so multi-roots trees. + + // Create tree and fill it. + util::ctree::ctree<I> tree; + { + tree.reserve(f, args.n_node); + + mln_pixter(const I) pix(f); + for_all(pix) + { + int idx = args.map.element(pix.offset()); + int loc = args.n_node - args.location[idx] - 1; + tree.data_->map_.element(pix.offset()) = loc; + tree.data_->values_[loc] = pix.val(); + } + + for (unsigned idx = 0; idx < args.n_node; ++idx) + { + int loc = args.n_node - args.location[idx] - 1; + tree.data_->parent_[loc] = args.n_node - args.location[args.parent[idx]] - 1; + tree.data_->length_[loc] = args.area[idx]; + } + } + + trace::exiting("mln::morpho::tree::impl::hqueue_fast"); + 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_IMPL_HQUEUE_FAST_HH diff --git a/milena/sandbox/edwin/mln/morpho/tree/impl/union_find.hh b/milena/mln/morpho/tree/impl/union_find.hh similarity index 99% copy from milena/sandbox/edwin/mln/morpho/tree/impl/union_find.hh copy to milena/mln/morpho/tree/impl/union_find.hh index f2aee82..8485727 100644 --- a/milena/sandbox/edwin/mln/morpho/tree/impl/union_find.hh +++ b/milena/mln/morpho/tree/impl/union_find.hh @@ -163,7 +163,9 @@ namespace mln } } + # endif // ! MLN_INCLUDE_ONLY + } } diff --git a/milena/mln/morpho/tree/impl/union_find_fast.hh b/milena/mln/morpho/tree/impl/union_find_fast.hh new file mode 100644 index 0000000..80a521d --- /dev/null +++ b/milena/mln/morpho/tree/impl/union_find_fast.hh @@ -0,0 +1,396 @@ +// 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_UNION_FIND_FAST_HH +# define MLN_MORPHO_TREE_IMPL_UNION_FIND_FAST_HH + +/// +/// + + +namespace mln +{ + + namespace morpho + { + + namespace tree + { + + namespace impl + { + + namespace generic + { + + /// The fast version of union-find algorithm. + /// \pre The image \p f must be fastest + /// + /// \param tree_tag The type of tree to build (min or max-tree). + /// \param f Source image used to build the component tree. + /// \param nbh Neighborhood used to build the component tree. + /// + /// \return A component tree object. + /// + template <typename T, typename I, typename N> + util::ctree::ctree<I> + union_find_fast(const tag::tree::tree_t<T>& tree_tag, + const Image<I>& f, + const Neighborhood<N>& nbh); + } + + } + +# ifndef MLN_INCLUDE_ONLY + + namespace internal + { + + template <typename I> + inline + void + sort_sites_set_radix(const tag::tree::max_t&, const I& f, util::array<unsigned>& s, util::array<unsigned>& levels) + { + typedef mln_value(I) V; + static const value::set<V>& vset = value::set<V>::the(); + static const int n = mln_card(V); + + util::array<unsigned> h(n, 0); + s.resize(f.domain().nsites()); + levels.resize(n + 1); + + mln_fwd_pixter(const I) pxl(f); + for_all(pxl) + ++h[vset.index_of(pxl.val())]; + + for (int i = 1; i <= n; ++i) + { + levels[i] = levels[i - 1] + h[i - 1]; + h[i - 1] = levels[i - 1]; + } + + // computing output data + for_all(pxl) + { + int i = vset.index_of(pxl.val()); + s[h[i]++] = pxl.offset(); + } + } + + template <typename I> + static inline + void + sort_sites_set_radix(const tag::tree::min_t&, const I& f, util::array<unsigned>& s, util::array<unsigned>& levels) + { + typedef mln_value(I) V; + static const value::set<V>& vset = value::set<V>::the(); + static const int n = mln_card(V); + + util::array<unsigned> h(n, 0); + s.resize(f.domain().nsites()); + levels.resize(n + 1); + + mln_fwd_pixter(const I) pxl(f); + for_all(pxl) + ++h[n - 1 - vset.index_of(pxl.val())]; + + for (int i = 1; i <= n; ++i) + { + levels[i] = levels[i - 1] + h[i - 1]; + h[i - 1] = levels[i - 1]; + } + + // computing output data + for_all(pxl) + { + int i = vset.index_of(pxl.val()); + s[h[n - 1 - i]++] = pxl.offset(); + } + } + + template <typename T, typename I> + struct sort_op_less_; + + template <typename I> + struct sort_op_less_<tag::tree::max_t, I> + { + sort_op_less_ (const I& f_) : f(f_) {} + + bool operator() (unsigned e1, unsigned e2) + { + return f.element(e1) < f.element(e2) || (f.element(e1) == f.element(e2) && e1 < e2); + } + + private: + const I& f; + }; + + template <typename I> + struct sort_op_less_<tag::tree::min_t, I> + { + sort_op_less_ (const I& f_) : f(f_) {} + + bool operator() (unsigned e1, unsigned e2) + { + return f.element(e2) < f.element(e1) || (f.element(e1) == f.element(e2) && e2 < e1); + } + + private: + const I& f; + }; + + template <typename T, typename I> + inline + void + sort_sites_set_std(const T&, const I& f, util::array<unsigned>& s) + { + s.reserve(f.domain().nsites()); + + mln_pixter(const I) pix(f); + for_all(pix) + s.append(pix.offset()); + + std::vector<unsigned>& v = s.hook_std_vector_(); + std::sort(v.begin(), v.end(), sort_op_less_<T, I>(f)); + } + + template <typename T, typename I> + static inline + void + sort_sites_set_dispatch(const mln::trait::value::quant::low&, const T& tree_, const I& f, util::array<unsigned>& s, util::array<unsigned>& levels) + { + trace::entering("morpho::tree::internal::sort_sites_set_radix"); + sort_sites_set_radix(tree_, f, s, levels); + trace::exiting("morpho::tree::internal::sort_sites_set_radix"); + } + + template <typename T, typename I> + static inline + void + sort_sites_set_dispatch(const mln::trait::value::quant::any&, const T& tree_, const I& f, util::array<unsigned>& s, util::array<unsigned>& levels) + { + (void) levels; + trace::entering("morpho::tree::internal::sort_sites_set_std"); + sort_sites_set_std(tree_, f, s); + trace::exiting("morpho::tree::internal::sort_sites_set_std"); + } + + + template <typename T, typename I> + static inline + void + sort_sites_set(const T& tree_, const I& f, util::array<unsigned>& s, util::array<unsigned>& levels) + { + sort_sites_set_dispatch(mln_trait_value_quant(mln_value(I)) (), tree_, f, s, levels); + } + + template <typename T> + static inline + unsigned + zfind_root(T& zpar, unsigned x) + { + mlc_equal(mln_value(T), unsigned)::check(); + if (zpar.element(x) == x) + return x; + else + return zpar.element(x) = zfind_root(zpar, zpar.element(x)); + } + + + } // end of namespace mln::morpho::tree::internal + + namespace impl + { + + namespace generic + { + + template <typename I, typename F> + static inline + void + union_find_fast_process2_(const I& f, const util::array<int>& dp, + F& parent, F& area, F& zpar, + unsigned& nnodes, unsigned p) + { + // Make set. + const unsigned n_nbhs = dp.nelements(); + + parent.element(p) = p; + zpar.element(p) = p; + + for (unsigned j = 0; j < n_nbhs; ++j) + { + unsigned n = p + dp[j]; + + if (zpar.element(n) == 0) // not deja-vu + continue; + + unsigned r = internal::zfind_root(zpar, n); + if (r != p) + { + parent.element(r) = p; + zpar.element(r) = p; + area.element(p) += area.element(r); + if (f.element(r) != f.element(p)) + ++area.element(p); + else + --nnodes; + } + } + } + + + template <typename I, typename F> + static inline + unsigned + union_find_fast_process(const mln::trait::value::quant::low&, const I& f, const util::array<int>& dp, + F& parent, F& area, F& zpar, + const util::array<unsigned>& s, const util::array<unsigned>& levels) + { + trace::entering("morpho::tree::impl::union_find_fast_process_low"); + static const unsigned nvalues = mln_card(mln_value(I)); + + unsigned nnodes = f.domain().nsites(); + + for (int v = nvalues - 1; v >= 0; --v) // Reverse + { + for (int i = levels[v + 1] - 1; i >= (int)levels[v]; --i) // Reverse + union_find_fast_process2_(f, dp, parent, area, zpar, nnodes, s[i]); + + // Fast root compression + for (unsigned i = levels[v]; i < levels[v + 1]; ++i) // Direct + { + unsigned p = s[i]; + zpar.element(p) = zpar.element(zpar.element(p)); + } + } + + trace::exiting("morpho::tree::impl::union_find_fast_process_low"); + return nnodes; + } + + template <typename I, typename F> + static inline + unsigned + union_find_fast_process(const trait::value::quant::any&, const I& f, const util::array<int>& dp, + F& parent, F& area, F& zpar, + const util::array<unsigned>& s, const util::array<unsigned>&) + { + trace::entering("morpho::tree::impl::union_find_fast_process_any"); + + unsigned nnodes = f.domain().nsites(); + for (int i = s.size() - 1; i >= 0; --i) // Reverse + union_find_fast_process2_(f, dp, parent, area, zpar, nnodes, s[i]); + + trace::exiting("morpho::tree::impl::union_find_fast_process_any"); + return nnodes; + } + + + template <typename T, typename I, typename N> + util::ctree::ctree<I> + union_find_fast(const tag::tree::tree_t<T>&, + const Image<I>& f_, + const Neighborhood<N>& nbh_) + { + trace::entering("morpho::tree::impl::generic::union_find_fast"); + + typedef p_array<mln_psite(I)> S; + typedef mln_value(I) V; + + const I& f = exact(f_); + const N& nbh = exact(nbh_); + util::array< unsigned > s; + util::array< unsigned > levels; + internal::sort_sites_set(T (), f, s, levels); + + util::array<int> dp = offsets_wrt(f, nbh); + + // Auxiliary data. + mln_ch_value(I, unsigned) parent; + mln_ch_value(I, unsigned) area; + mln_ch_value(I, unsigned) zpar; + + initialize(parent, f); + initialize(area, f); + initialize(zpar, f); + + // Initialization. + data::fill(zpar, 0); + data::fill(area, 0); + + unsigned nnodes = + union_find_fast_process(mln_trait_value_quant(V) (), f, dp, parent, area, zpar, s, levels); + + util::ctree::ctree<I> tree; + tree.reserve(f, nnodes); + + unsigned root_offset = 0; + for (unsigned i = 0; i < s.size(); ++i) // Forward + { + unsigned p = s[i]; + unsigned q = parent.element(p); + + if (f.element(parent.element(q)) == f.element(q)) // Canonization + q = (parent.element(p) = parent.element(q)); + + if (f.element(p) == f.element(q) && p != q) // Not a node. + { + mln_assertion(q == parent.element(p)); + tree.data_->map_.element(p) = + tree.data_->map_.element(q); + } + else + { + unsigned& offset = (p == q) ? root_offset : area.element(q); + + // Insert Node. + mln_assertion(offset < nnodes); + tree.data_->map_.element(p) = offset; + tree.data_->parent_[offset] = tree.data_->map_.element(q); + tree.data_->values_[offset] = f.element(p); + tree.data_->length_[offset] = area.element(p); + area.element(p) = offset + 1; + offset += tree.data_->length_[offset] + 1; + } + } + + trace::exiting("morpho::tree::impl::generic::union_find_fast"); + 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_IMPL_UNION_FIND_FAST_HH diff --git a/milena/mln/morpho/tree/leaf_last.hh b/milena/mln/morpho/tree/leaf_last.hh new file mode 100644 index 0000000..2f283d5 --- /dev/null +++ b/milena/mln/morpho/tree/leaf_last.hh @@ -0,0 +1,111 @@ +// 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_LEAF_LAST_HH +# define MLN_MORPHO_TREE_LEAF_LAST_HH + +# include <mln/accu/stat/max.hh> +# include <mln/util/ctree/ctree.hh> +# include <mln/morpho/tree/propagate_node.hh> +# include <mln/convert/to_p_array.hh> + +namespace mln +{ + + namespace morpho + { + + namespace tree + { + + template <typename T, typename V> + util::array<typename T::node_t> + leaf_last(const attribute_image<T, V>& ima); + +# ifndef MLN_INCLUDE_ONLY + + namespace internal + { + + template <typename T, typename V> + struct less_ + { + typedef typename T::node_t node_t; + + less_(const attribute_image<T, V>& img) : + f_ (img) + { + } + + bool operator() (const node_t& lhs, const node_t& rhs) + { + return f_(lhs) > f_(rhs); + } + + private: + const attribute_image<T, V>& f_; + }; + } + + template <typename T, typename V> + util::array<typename T::node_t> + leaf_last(const attribute_image<T, V>& ima) + { + typedef attribute_image<T, V> A; + typedef typename T::node_t node_t; + + p_array<node_t> sorted_sites = convert::to_p_array(ima); + std::vector<node_t>& hook = sorted_sites.hook_std_vector_(); + std::sort(hook.begin(), hook.end(), internal::less_<T, V> (ima)); + + mln_ch_value(A, char) deja_vu; + initialize(deja_vu, ima); + data::fill(deja_vu, 0); + + util::array<node_t> result; + + mln_piter(p_array<node_t>) n(sorted_sites); + for_all(n) + { + if (deja_vu(n)) + continue; + result.append(n); + deja_vu(n) = true; // useless + propagate_node_to_descendants<T, char>(deja_vu, n, 1); + propagate_node_to_ancestors<T, char>(deja_vu, n, 1); + } + + return result; + } + +# endif // ! MLN_INCLUDE_ONLY + + } // end of namespace mln::morpho::tree + + } // end of namespace mln::morpho + +} // end of namespace mln + +#endif // !MLN_MORPHO_TREE_LEAF_LAST_HH diff --git a/milena/mln/morpho/tree/propagate_node.hh b/milena/mln/morpho/tree/propagate_node.hh index 7778905..76f2410 100644 --- a/milena/mln/morpho/tree/propagate_node.hh +++ b/milena/mln/morpho/tree/propagate_node.hh @@ -26,14 +26,27 @@ #ifndef MLN_MORPHO_TREE_PROPAGATE_NODE_HH # define MLN_MORPHO_TREE_PROPAGATE_NODE_HH -# include <mln/core/concept/image.hh> -# include <mln/core/macros.hh> -# include <mln/morpho/tree/data.hh> - /// \file /// /// Functions to propagate node in the tree. +/// +/// \todo: The new implementation permits: + to propagate toward +/// subcomponents very fast (data::fill)\ + to propagate toward +/// supernodes as fast as the previous component\ + the attribute +/// image is still valid without the use of propagate_representative +/// that does not make sense any more. +/// +/// We are able in both cases to do something like: +/// + data::fill(attr_img | asc(n), v) +/// + data::fill(attr_img | desc(n), v) +/// +/// However the subimage morpher is not specialized for fast image and +/// does not take benefits in the fact that: +/// + desc(n) is a p_run (so can be memseté) +/// + attr_img can accessed by indexed (element() method for fast-image) +# include <mln/core/image/attribute_image.hh> +# include <algorithm> namespace mln { @@ -44,166 +57,52 @@ namespace mln namespace tree { - /*! - ** Propagate a value \p v from a node \p n to its descendants. - ** - ** \param[in] n Node to propagate. - ** \param[in] t Component tree used for propagation. - ** \param[in] a_ Attribute image where values are propagated. - ** \param[in] v Value to propagate. - ** \param[out] nb_leaves Optional. Store the number of leaves in - ** the component. - */ - template <typename T, typename A> - void - propagate_node_to_descendants(mln_psite(A) n, - const T& t, - Image<A>& a_, - const mln_value(A)& v, - unsigned* nb_leaves = 0); - - /*! - ** Propagate the node's value to its descendants. - ** - ** \param[in] n Node to propagate. - ** \param[in] t Component tree used for propagation. - ** \param[in] a_ Attribute image where values are propagated. - ** \param[out] nb_leaves Optional. Store the number of leaves in - ** the component. - */ - template <typename T, typename A> - inline - void - propagate_node_to_descendants(mln_psite(A)& n, - const T& t, - Image<A>& a_, - unsigned* nb_leaves = 0); - - - /*! - ** Propagate a value \p v from a node \p n to its ancestors. - ** - ** \param[in] n Node to propagate. - ** \param[in] t Component tree used for propagation. - ** \param[in] a_ Attribute image where values are propagated. - ** \param[in] v Value to propagate. - */ - template <typename T, typename A> - void - propagate_node_to_ancestors(mln_psite(A) n, - const T& t, - Image<A>& a_, - const mln_value(A)& v); - - /*! - ** Propagate the node's value to its ancestors. - ** - ** \param[in] n Node to propagate. - ** \param[in] t Component tree used for propagation. - ** \param[in,out] a_ Attribute image where values are propagated. - */ - template <typename T, typename A> - inline - void - propagate_node_to_ancestors(mln_psite(A) n, - const T& t, - Image<A>& a_); - - -# ifndef MLN_INCLUDE_ONLY - - /* Descendants propagation */ - - template <typename T, typename A> + template <typename T, typename V> inline - void - propagate_node_to_descendants(mln_psite(A) n, - const T& t, - Image<A>& a_, - const mln_value(A)& v, - unsigned* nb_leaves = 0) + unsigned + propagate_node_to_descendants(attribute_image<T, V>& a, + const typename T::node_t& n, + const V& value) { - A& a = exact(a_); + typedef typename T::nodes_t desc_t; mln_precondition(a.is_valid()); - mln_precondition(a.domain() == t.f().domain()); - mln_precondition(a.domain().has(n)); - - - if (!t.is_a_node(n)) // Get the representant. - n = t.parent(n); - mln_assertion(t.is_a_node(n)); - - typename T::depth1st_piter pp(t, n); - - pp.start(); // We don't set n to v. + mln_precondition(a.has(n)); - if (nb_leaves) - *nb_leaves += t.is_a_leaf(pp); + // (n € desc(n)) but we don't set its value. + unsigned i = a.tree().desc(n).nsites(); - for (pp.next(); pp.is_valid(); pp.next()) + if (i > 1) { - a(pp) = v; - if (nb_leaves && t.is_a_leaf(pp)) - ++(*nb_leaves); + mln_assertion(a.tree().has_index(n.index() + i)); + V* v = a.buffer() + n.index(); + std::fill(v + 1, v + i, value); } + return i - 1; } - - template <typename T, typename A> + template <typename T, typename V> inline - void - propagate_node_to_descendants(mln_psite(A) n, - const T& t, - Image<A>& a_, - unsigned* nb_leaves = 0) - - { - A& a = exact(a_); - propagate_node_to_descendants(n, t, a, a(n), nb_leaves); - } - - - /* Ancestors propagation */ - - template <typename T, typename A> - void - propagate_node_to_ancestors(mln_psite(A) n, - const T& t, - Image<A>& a_, - const mln_value(A)& v) + unsigned + propagate_node_to_ancestors(attribute_image<T, V>& a, + typename T::node_t n, + const V& v) { - A& a = exact(a_); mln_precondition(a.is_valid()); - mln_precondition(a.domain() == t.f().domain()); - mln_precondition(a.domain().has(n)); + mln_precondition(a.has(n)); - if (!t.is_a_node(n)) // Get the representant. - n = t.parent(n); - mln_assertion(t.is_a_node(n)); - - if (t.is_root(n)) - return; + if (a.tree().is_a_root(n)) + return 0; + unsigned i = 0; do { - n = t.parent(n); + n = a.tree().parent(n); a(n) = v; - } while (!t.is_root(n)); - - } + ++i; + } while (!a.tree().is_a_root(n)); - template <typename T, typename A> - inline - void - propagate_node_to_ancestors(mln_psite(A) n, - const T& t, - Image<A>& a_) - { - A& a = exact(a_); - propagate_node_to_ancestors(n, t, a, a(n)); + return i; } -# endif // ! MLN_INCLUDE_ONLY - } // end of namespace mln::morpho::tree } // end of namespace mln::morpho diff --git a/milena/mln/tag/skeleton.hh b/milena/mln/tag/skeleton.hh index f446317..bca9b4f 100644 --- a/milena/mln/tag/skeleton.hh +++ b/milena/mln/tag/skeleton.hh @@ -44,6 +44,7 @@ namespace mln template <typename F> struct function_ { typedef F param; }; template <typename G> struct graph_ { typedef G param; }; template <typename I> struct image_ { typedef I param; }; + template <typename T> struct tree_ { typedef T param; }; template <typename N> struct neighb_ { typedef N param; }; template <typename P> struct psite_ { typedef P param; }; template <typename S> struct domain_ { typedef S param; }; diff --git a/milena/mln/util/hqueues.hh b/milena/mln/util/hqueues.hh index 151c71c..35c7a50 100644 --- a/milena/mln/util/hqueues.hh +++ b/milena/mln/util/hqueues.hh @@ -55,6 +55,7 @@ namespace mln nvalues = mln_card(T) }; + hqueues(); hqueues(const histo::array<T>& h); const p_queue_fast<P>& operator[](unsigned i) const; @@ -66,38 +67,29 @@ namespace mln const mln::value::set<T>& vset() const; protected: - void pre_allocate_(unsigned i); - - const histo::array<T>& h_; const mln::value::set<T>& s_; - std::vector<bool> allocated_; std::vector< p_queue_fast<P> >queues_; }; # ifndef MLN_INCLUDE_ONLY - template <typename P, typename T> inline - hqueues<P,T>::hqueues(const histo::array<T>& h) - : h_ (h), - s_ (mln::value::set<T>::the()), - allocated_ (nvalues, false), - queues_ (nvalues) + hqueues<P,T>::hqueues() + : s_ (mln::value::set<T>::the()) { + queues_.resize(nvalues); } + template <typename P, typename T> inline - void - hqueues<P,T>::pre_allocate_(unsigned i) + hqueues<P,T>::hqueues(const histo::array<T>& h) + : s_ (mln::value::set<T>::the()) { - mln_precondition(i < nvalues); - if (!allocated_[i]) - { - queues_[i].reserve(h_[i]); - allocated_[i] = true; - } + queues_.resize(nvalues); + for (unsigned i = 0; i < nvalues; ++i) + queues_[i].reserve(h[i]); } @@ -107,7 +99,6 @@ namespace mln hqueues<P,T>::operator[](unsigned i) const { mln_precondition(i < nvalues); - pre_allocate_(i); return queues_[i]; } @@ -117,7 +108,6 @@ namespace mln hqueues<P,T>::operator[](unsigned i) { mln_precondition(i < nvalues); - pre_allocate_(i); return queues_[i]; } @@ -127,7 +117,6 @@ namespace mln hqueues<P,T>::operator()(const T& v) const { unsigned i = s_.index_of(v); - pre_allocate_(i); return queues_[i]; } @@ -137,7 +126,6 @@ namespace mln hqueues<P,T>::operator()(const T& v) { unsigned i = s_.index_of(v); - pre_allocate_(i); return queues_[i]; } -- 1.5.6.5