[Olena-patches] [PATCH 2/3] Improve the integration of the topological watershed transform.

* mln/morpho/watershed/topo_wst.hh: Rename as... * mln/morpho/watershed/topological.hh: ...this. Add documentation header. * mln/morpho/watershed/all.hh: Include mln/morpho/watershed/topological.hh. * tests/morpho/watershed/test_watershed_topo.cc: Adjust. Revamp. Add copyright and documentation headers. Rename as... * tests/morpho/watershed/topological.cc: ...this. * tests/morpho/watershed/Makefile.am (check_PROGRAMS): Add topological. (topological_SOURCES): New. --- milena/ChangeLog | 18 + milena/mln/morpho/watershed/all.hh | 2 + milena/mln/morpho/watershed/topo_wst.hh | 768 ------------------- milena/mln/morpho/watershed/topological.hh | 789 ++++++++++++++++++++ milena/tests/morpho/watershed/Makefile.am | 4 +- .../tests/morpho/watershed/test_watershed_topo.cc | 42 - .../morpho/watershed/topological.cc} | 49 +- 7 files changed, 832 insertions(+), 840 deletions(-) delete mode 100644 milena/mln/morpho/watershed/topo_wst.hh create mode 100644 milena/mln/morpho/watershed/topological.hh delete mode 100644 milena/tests/morpho/watershed/test_watershed_topo.cc copy milena/{mln/morpho/watershed/all.hh => tests/morpho/watershed/topological.cc} (66%) diff --git a/milena/ChangeLog b/milena/ChangeLog index 9789d98..58b58ff 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,23 @@ 2009-09-24 Roland Levillain &lt;roland(a)lrde.epita.fr&gt; + Improve the integration of the topological watershed transform. + + * mln/morpho/watershed/topo_wst.hh: Rename as... + * mln/morpho/watershed/topological.hh: ...this. + Add documentation header. + * mln/morpho/watershed/all.hh: Include + mln/morpho/watershed/topological.hh. + * tests/morpho/watershed/test_watershed_topo.cc: Adjust. + Revamp. + Add copyright and documentation headers. + Rename as... + * tests/morpho/watershed/topological.cc: ...this. + * tests/morpho/watershed/Makefile.am (check_PROGRAMS): Add + topological. + (topological_SOURCES): New. + +2009-09-24 Roland Levillain &lt;roland(a)lrde.epita.fr&gt; + Import Alexandre's implementation of the topological watershed. * mln/morpho/watershed/topo_wst.hh: New file, imported from diff --git a/milena/mln/morpho/watershed/all.hh b/milena/mln/morpho/watershed/all.hh index a07e75a..8ce71b2 100644 --- a/milena/mln/morpho/watershed/all.hh +++ b/milena/mln/morpho/watershed/all.hh @@ -58,6 +58,8 @@ namespace mln # include <mln/morpho/watershed/flooding.hh> +# include <mln/morpho/watershed/topological.hh> + # include <mln/morpho/watershed/superpose.hh> diff --git a/milena/mln/morpho/watershed/topo_wst.hh b/milena/mln/morpho/watershed/topo_wst.hh deleted file mode 100644 index 7e9fb64..0000000 --- a/milena/mln/morpho/watershed/topo_wst.hh +++ /dev/null @@ -1,768 +0,0 @@ -// Copyright (C) 2008, 2009 EPITA Research and Development Laboratory (LRDE) -// -// This file is part of the Milena Library. This library is free -// software; you can redistribute it and/or modify it under the terms -// of the GNU General Public License version 2 as published by the -// Free Software Foundation. -// -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this library; see the file COPYING. If not, write to -// the Free Software Foundation, 51 Franklin Street, Fifth Floor, -// Boston, MA 02111-1307, USA. -// -// As a special exception, you may use this file as part of a free -// software library without restriction. Specifically, if other files -// instantiate templates or use macros or inline functions from this -// file, or you compile this file and link it with other files to -// produce an executable, this file does not by itself cause the -// resulting executable to be covered by the GNU General Public -// License. This exception does not however invalidate any other -// reasons why the executable file might be covered by the GNU General -// Public License. - -#ifndef MLN_MORPHO_TOPO_WST_HH -# define MLN_MORPHO_TOPO_WST_HH - -# include <vector> -# include <map> -# include <queue> - -# include <mln/core/site_set/p_set.hh> -# include <mln/core/site_set/p_priority.hh> -# include <mln/core/site_set/p_queue_fast.hh> - -# include <mln/util/ord.hh> - -# include <mln/data/sort_psites.hh> -# include <mln/data/fill.hh> - - -namespace mln -{ - namespace morpho - { - - template <class I> - mln_site(I) min (const Image<I> &ima_, p_set<mln_site(I)>& components) - { - const I& ima = exact(ima_); - - if (components.nsites() == 0) - return mln_site(I)(-1, -1); // FIXME - - typename p_set<mln_site(I)>::fwd_piter it(components); - mln_site(I) min = components[0]; - - for_all(it) - if (ima(it) < ima(min)) - min = it; - - return min; - } - - template <class I> - mln_site(I) max (p_set<mln_site(I)>& components) - { - if (components.nsites() == 0) - return mln_site(I)(-1, -1); // FIXME - - typename p_set<mln_site(I)>::fwd_piter it(components); - mln_site(I) max = components[0]; - - for_all(it) - if (ima(it) > ima(max)) - max = it; - - return max; - } - - - // Actually, this structure is a tree, despite its confusing name. - template <class I, class N> - struct topo_wst - { - - // Typedefs - // -------- - typedef mln_site(I) site; - typedef I image_t; - typedef N neighborhood_t; - - // Component tree management - // ------------------------- - private: - struct node { - mln_value(I) level; - int area; - int highest; - typedef mln_site(I) site; - // Can't modify the sites in a p_array - // p_array<mln_site(I)> children; - std::vector<site> children; - - void addChildren(const node& n) - { - // typename p_array<mln_site(I)>::fwd_piter it(n.children); - // for (it.start(); - // it.is_valid(); - // it.next()) - // children.append(it.to_site()); - for (unsigned i=0; i < n.children.size(); ++i) - children.push_back(n.children[i]); - } - - void addChild(const site p) - { - // children.append(n); - children.push_back(p); - } - }; // struct node - - public: - - mln_ch_value(I, site) Par_node; - mln_ch_value(I, site) Par_tree; - - mln_ch_value(I, int) Rnk_tree; - mln_ch_value(I, int) Rnk_node; - mln_ch_value(I, site) subtreeRoot; - mln_ch_value(I, node) nodes; - site Root; - - private: - - void MakeSet_tree(site x); - void MakeSet_node(site x); - site Find_tree(site x); - site Find_node(site x); - void BuildComponentTree(); - site MergeNode(site& node1, site& node2); - site Link_tree(site x, site y); - site Link_node(site x, site y); - node MakeNode(int level); - - private: - mln_ch_value(I, bool) isproc; - - // Ctor - public: - topo_wst(const Image<I>& i, - const Neighborhood<N>& n); - - public: - const I &ima; - const N &nbh; - - public: - void go(); - - private: - site min (p_set<site>& components); - site max (p_set<site>& components); - bool highest_fork (p_set<site>& components); - bool highest_fork (p_set<site>& components, site &r); - - // Optimized LCA Algorithm - - public: - site lca (site a, site b); - - private: - int *euler; - int *depth; - int ctree_size; - std::map<site, int, mln::util::ord<site> > pos; - site *repr; - - int *minim; - int **Minim; - - void compute_ctree_size (site p); - void build_euler_tour_rec(site p, int &position, int d); - void build_euler_tour(); - void build_minim (); - void removeOneSonNodes(site *p, mln_ch_value(I, site) &newPar_node); - void compressTree(); - }; // struct topo_wst - - - -# ifndef MLN_INCLUDE_ONLY - - // Ctor - template <class I, class N> - topo_wst<I, N>::topo_wst(const Image<I>& i, - const Neighborhood<N>& n) - : ima(exact(i)), nbh(exact(n)) - { - initialize(Par_node, i); - initialize(Par_tree, i); - initialize(Rnk_tree, i); - initialize(Rnk_node, i); - initialize(subtreeRoot, i); - initialize(nodes, i); - initialize(isproc, i); - } - - template <class I, class N> - void topo_wst<I, N>::MakeSet_tree(site x) - { - Par_tree(x) = x; - Rnk_tree(x) = 0; - } - - template <class I, class N> - void topo_wst<I, N>::MakeSet_node(site x) - { - Par_node(x) = x; - Rnk_node(x) = 0; - } - - template <class I, class N> - typename topo_wst<I, N>::site topo_wst<I, N>::Find_tree(site x) - { - if (Par_tree(x) != x) - Par_tree(x) = Find_tree(Par_tree(x)); - return Par_tree(x); - } - - template <class I, class N> - typename topo_wst<I, N>::site topo_wst<I, N>::Find_node(site x) - { - if (Par_node(x) != x) - Par_node(x) = Find_node(Par_node(x)); - return Par_node(x); - } - - template <class I, class N> - void topo_wst<I, N>::go() - { - BuildComponentTree(); - compressTree(); - - build_euler_tour(); - build_minim(); - } - - template <class I, class N> - void topo_wst<I, N>::BuildComponentTree() - { - // Init: - - // Sort the sites in increasing order - p_array<mln_site(I)> S; - S = data::sort_psites_increasing(ima); - - // Clear the marker map - data::fill(isproc, false); - for (int ip = 0; ip < int(S.nsites()); ++ip) - { - site p = S[ip]; - MakeSet_node(p); - MakeSet_tree(p); - // if (subtreeRoot.hold(p)) - subtreeRoot(p) = p; - // if (nodes.hold(p)) - nodes(p) = MakeNode(ima(p)); - } - - - - typename p_array<mln_site(I)>::fwd_piter ip(S); - for_all(ip) - { - site p = ip.to_site(); - - site curCanonicalElt = Find_tree(p); - site curNode = Find_node(subtreeRoot(curCanonicalElt)); - - // FIXME: Should be `n' instead of `q'. - mln_niter(N) q(nbh, ip); - for_all(q) - if (ima.has(q) and isproc(q) and ima(q) <= ima(p)) - { - site adjCanonicalElt = Find_tree(q); - site adjNode = Find_node(subtreeRoot(adjCanonicalElt)); - if (curNode != adjNode) - { - if (nodes(curNode).level == nodes(adjNode).level) - curNode = MergeNode(adjNode, curNode); - else - { - nodes(curNode).addChild(adjNode); - nodes(curNode).area += nodes(adjNode).area; - nodes(curNode).highest += nodes(adjNode).highest; - } - } - - curCanonicalElt = Link_tree(adjCanonicalElt, curCanonicalElt); - subtreeRoot(curCanonicalElt) = curNode; - } - isproc(p) = true; - } - // Pour garder une map de correspondance site <-> local_root - // for (int ip = 0; ip < int(S.size()); ++ip) - // { - // site p = S[ip]; - // M(p) = Find_node(p); - // } - - mln_piter(I) r(Par_node.domain()); - for_all(r) - Par_node(r) = Find_node(r); - - // Find the ``first'' site of ima, according to the forward - // traversal order. - mln_fwd_piter(I) rp(Par_node.domain());; - rp.start(); - - Root = subtreeRoot(Find_tree(Find_node(rp))); - } - - - template <class I, class N> - typename topo_wst<I, N>::site topo_wst<I, N>::MergeNode(site& node1, - site& node2) - { - site tmpNode = Link_node(node1, node2); - site tmpNode2; - if (tmpNode == node2) - { - nodes(node2).addChildren(nodes(node1)); - tmpNode2 = node1; - } - else - { - nodes(node1).addChildren(nodes(node2)); - tmpNode2 = node2; - } - nodes(tmpNode).area += nodes(tmpNode2).area; - nodes(tmpNode).highest += nodes(tmpNode2).highest; - return tmpNode; - } - - template <class I, class N> - typename topo_wst<I, N>::site topo_wst<I, N>::Link_tree(site x, site y) - { - if (Rnk_tree(x) > Rnk_tree(y)) - std::swap(x, y); - else - if (Rnk_tree(x) == Rnk_tree(y)) - Rnk_tree(y) += 1; - Par_tree(x) = y; - return y; - } - - template <class I, class N> - typename topo_wst<I, N>::site topo_wst<I, N>::Link_node(site x, site y) - { - if (Rnk_node(x) > Rnk_node(y)) - std::swap(x, y); - else - if (Rnk_node(x) == Rnk_node(y)) - Rnk_node(y) += 1; - Par_node(x) = y; - return y; - } - - template <class I, class N> - typename topo_wst<I, N>::node topo_wst<I, N>::MakeNode(int level) - { - node n; - n.level = level; - n.area = 1; - n.highest = level; - return n; - } - - - template <class I, class N> - bool topo_wst<I, N>::highest_fork (p_set<site>& components, site &r) - { - if (components.nsites() == 0) - { - std::cerr << "highest fork : empty set" << std::endl; - return false; - } - - site - m = min(components), - hfork = m; - - typename p_set<site>::fwd_piter it(components); - for_all(it) - { - // Can't remove the site from the set - if (it.to_site() == m) - continue; - - site c = lca(hfork, it.to_site()); - if (c != it.to_site()) - hfork = c; - } - - if (nodes(m).level == nodes(hfork).level) - return false; - - r = hfork; - return true; - } - - template <class I, class N> - bool topo_wst<I, N>::highest_fork (p_set<site>& components) { - site i; - return highest_fork(components, i); - } - - template <class I, class N> - void topo_wst<I, N>::compute_ctree_size (site p) - { - ctree_size += 1; - node& n = nodes(p); - - // typename p_array<mln_site(I)>::fwd_piter it(n.children); - // for_all(it) - // compute_ctree_size(it.to_site()); - - for (unsigned i=0; i < n.children.size(); ++i) - compute_ctree_size(n.children[i]); - } - - - template <class I, class N> - void topo_wst<I, N>::build_euler_tour_rec(site p, int &position, int d) - { - if (pos.find(p) == pos.end()) - pos[p] = position; - - repr[position] = p; - depth[position] = d; - euler[position] = pos[p]; - ++position; - node& n = nodes(p); - - // typename p_array<mln_site(I)>::fwd_piter it(n.children); - // for_all(it) - // { - // build_euler_tour_rec(it.to_site(), position, d+1); - // depth[position] = d; // Not optimized - // euler[position] = pos[p]; - // repr[position] = p; // Pas necessaire? - // ++position; - // } - - for (unsigned i=0; i < n.children.size(); ++i) - { - build_euler_tour_rec(n.children[i], position, d+1); - depth[position] = d; // Not optimized - euler[position] = pos[p]; - repr[position] = p; // Pas necessaire? - ++position; - } - } - - - template <class I, class N> - void topo_wst<I, N>::build_euler_tour () - { - ctree_size = 0; - compute_ctree_size(Root); - - int size = 2 * ctree_size - 1; - - // FIXME : free this - euler = new int[size]; - depth = new int[size]; - repr = new site[size]; - - int position = 0; - build_euler_tour_rec(Root, position, 0); - } - - - template <class I, class N> - void topo_wst<I, N>::build_minim () - { - // compute_tree_size(); // Already done in build_euler_tour - int size = 2 * ctree_size - 1; - int logn = (int)(ceil(log((double)(size))/log(2.0))); - // minim.reserve(size); - minim = new int [logn * size]; // FIXME : Think about freeing this - Minim = new int* [logn]; - - Minim[0] = minim; - - // Init - for (int i = 0; i < size - 1; ++i) - if (depth[euler[i]] < depth[euler[i+1]]) { - Minim[0][i] = i; - } else { - Minim[0][i] = i+1; - } - - // Minim[0][size - 1] = size - 1; - - int k; - for (int j = 1; j < logn; j++) { - k = 1 << (j - 1); - Minim[j] = &minim[j * size]; - for (int i = 0; i < size; i++) { - if ((i + (k << 1)) >= size) { - //Minim[j][i] = size - 1; - } - else { - if (depth[euler[Minim[j - 1][i]]] - <= depth[euler[Minim[j - 1][i + k]]]) - Minim[j][i] = Minim[j - 1][i]; - else - Minim[j][i] = Minim[j - 1][i + k]; - } - } - } - - } // void build_minim () - - - template <class I, class N> - typename topo_wst<I, N>::site topo_wst<I, N>::lca (site a, site b) - { - int - m = pos[a], - n = pos[b], - k; - - if (m == n) - return repr[m]; - - if (m > n) - { - k = n; - n = m; - m = k; - } - - k = (int)(log((double)(n - m))/log(2.)); - - if (depth[euler[Minim[k][m]]] < depth[euler[Minim[k][n - (1 << k)]]]) { - return repr[euler[Minim[k][m]]]; - } else { - return repr[euler[Minim[k][n - (1 << k)]]]; - } - } - - - template <class I, class N> - void topo_wst<I, N>::removeOneSonNodes(site *p, - mln_ch_value(I, site) &newPar_node) - { - node &n = nodes(*p); - - if (n.children.size() == 1) // this node has 1 son, delete it - { - n.area = -1; - newPar_node(*p) = n.children[0]; - *p = n.children[0]; - removeOneSonNodes(p, newPar_node); - } - else // there is more than one son, recursive call - { - for (unsigned i = 0; i < n.children.size(); ++i) - removeOneSonNodes(&(n.children[i]), newPar_node); - } - } - - - template <class I, class N> - void topo_wst<I, N>::compressTree() - { - mln_ch_value(I, site) newPar_node; - initialize(newPar_node, Par_node); - - // Remove the nodes with one son - removeOneSonNodes(&Root, newPar_node); - - // Update the references on deleted nodes - mln_piter(I) p(Par_node.domain()); - for_all(p) - while (nodes(Par_node(p)).area == -1) - Par_node(p) = newPar_node(Par_node(p)); - } - - template <class T> - bool w_constructible(T &tree, typename T::site p, typename T::site &r) - { - - typedef typename T::image_t I; - typedef typename T::neighborhood_t N; - - const I &ima = exact(tree.ima); - const N &nbh = exact(tree.nbh); - - // FIXME: Should be `n' instead of `q'. - mln_niter(N) q(nbh, p); - p_set<mln_site(I)> v; - - for_all(q) - // FIXME: Shouldn't it be: `ima.has(q)' instead of - // `ima.domain().has(q)'? - if (ima.domain().has(q) && ima(q) < ima(p)) - v.insert(tree.Par_node(q)); - - if (v.nsites() == 0) - return false; - - if (v.nsites() == 1) - { - r = v[0]; - return true; - } - - mln_site(I) c = min(ima, v); - mln_site(I) cmin = c; - - typename p_set<mln_site(I)>::fwd_piter it(v); - for_all(it) - { - // Can't remove the site from the set - if (it.to_site() == cmin) - continue; - - mln_site(I) c1 = tree.lca(c, it); - - if (c1 != it) - c = c1; - } - - if (tree.nodes(c).level >= ima(p)) - return false; - - r = c; - return true; - } - - template <class T> - bool w_constructible(T &tree, typename T::site p) { - typename T::site r; - return w_constructible(tree, p, r); - } - - - - template <class T> - typename T::image_t topo_watershed(T &tree) - { - - typedef typename T::image_t I; - typedef typename T::neighborhood_t N; - - I ima = exact(tree.ima); - const N &nbh = exact(tree.nbh); - - // Maxima components - mln_ch_value(I, bool) cmax; - initialize(cmax, ima); - - // Mark enqueued sites - mln_ch_value(I, bool) enqueued; - initialize(enqueued, ima); - - p_priority< mln_value(I), p_queue_fast<mln_site(I)> > l; - // p_queue < site > m; - std::queue<mln_site(I)> m; - mln_value(I) min = mln_min(mln_value(I)); - - std::cout << "Init" << std::endl; - - // Flag C-maxima - data::fill(cmax, false); - data::fill(enqueued, false); - - mln_piter(I) it(tree.Par_node.domain()); - for_all(it) - // if (nodes(Par_node(it.to_site())).children.nsites() == 0) - if (tree.nodes(tree.Par_node(it)).children.size() == 0) - { - cmax(it) = true; - m.push(it); - } - - while (!m.empty()) - { - // FIXME: Should be `n' instead of `q'. - mln_niter(N) q(nbh, m.front()); - // FIXME: Shouldn't it be: `cmax.has(q)' instead of - // `cmax.domain().has(q)'? - for_all(q) - if (cmax.domain().has(q) && !cmax(q) && !enqueued(q)) - { - enqueued(q) = true; - l.push(mln_max(mln_value(I)) - ima(q), q); - } - m.pop(); - } - - - std::cout << "end" << std::endl; - - // Main loop - while (!l.is_empty()) - { - mln_site(I) x = l.front(); - l.pop(); - enqueued(x) = false; - - mln_site(I) c; - bool is_w = w_constructible(tree, x, c); - - if (is_w) - { - ima(x) = tree.nodes(c).level; - tree.Par_node(x) = c; - - // if (nodes(c).children.nsites() == 0) - if (tree.nodes(c).children.size() == 0) - cmax(x) = true; - else - // if (nodes(c).children.nsites() > 1) - if (tree.nodes(c).children.size() == 1) - std::cerr << "ERREUR COMPOSANTE BRANCHE " - << tree.nodes(c).children.size() << std::endl; - - - // FIXME: Should be `n' instead of `q'. - mln_niter(N) q(nbh, x); - // FIXME: Shouldn't it be: `ima.has(q)' instead of - // `ima.domain().has(q)'? - for_all(q) - if (ima.domain().has(q) && !cmax(q) && !enqueued(q)) - { - enqueued(q) = true; - l.push(mln_max(mln_value(I)) - ima(q), q); // FIXME: - // Just - // invert - // the - // priority. - } - } - } // while(!l.empty()) - - for_all(it) - ima(it) = tree.nodes(tree.Par_node(it)).level; - - return ima; - } - - -# endif // MLN_INCLUDE_ONLY - - - }; // namespace morpho - -}; // namespace mln - -#endif // ! MLN_MORPHO_TOPO_WST_HH diff --git a/milena/mln/morpho/watershed/topological.hh b/milena/mln/morpho/watershed/topological.hh new file mode 100644 index 0000000..8fe890c --- /dev/null +++ b/milena/mln/morpho/watershed/topological.hh @@ -0,0 +1,789 @@ +// Copyright (C) 2008, 2009 EPITA Research and Development Laboratory (LRDE) +// +// This file is part of the Milena Library. This library is free +// software; you can redistribute it and/or modify it under the terms +// of the GNU General Public License version 2 as published by the +// Free Software Foundation. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this library; see the file COPYING. If not, write to +// the Free Software Foundation, 51 Franklin Street, Fifth Floor, +// Boston, MA 02111-1307, USA. +// +// As a special exception, you may use this file as part of a free +// software library without restriction. Specifically, if other files +// instantiate templates or use macros or inline functions from this +// file, or you compile this file and link it with other files to +// produce an executable, this file does not by itself cause the +// resulting executable to be covered by the GNU General Public +// License. This exception does not however invalidate any other +// reasons why the executable file might be covered by the GNU General +// Public License. + +#ifndef MLN_MORPHO_WATERSHED_TOPOLOGICAL_HH +# define MLN_MORPHO_WATERSHED_TOPOLOGICAL_HH + +/// \file +/// +/// \brief Topological Watershed Transform (WST) algorithm. +/// +/// Reference: +/// Michel Couprie, Laurent Najman and Gilles Bertrand. +/// Quasi-linear algorithms for the topological watershed. In: +/// Journal of Mathematical Imaging and Vision (JMIV), vol. 22, +/// no. 2-3, pages 231--249, 2005. + +# include <vector> +# include <map> +# include <queue> + +# include <mln/core/site_set/p_set.hh> +# include <mln/core/site_set/p_priority.hh> +# include <mln/core/site_set/p_queue_fast.hh> + +# include <mln/util/ord.hh> + +# include <mln/data/sort_psites.hh> +# include <mln/data/fill.hh> + + +// FIXME: Clean up and document. + +namespace mln +{ + + namespace morpho + { + + namespace watershed + { + + template <class I> + mln_site(I) min (const Image<I> &ima_, p_set<mln_site(I)>& components) + { + const I& ima = exact(ima_); + + if (components.nsites() == 0) + return mln_site(I)(-1, -1); // FIXME + + typename p_set<mln_site(I)>::fwd_piter it(components); + mln_site(I) min = components[0]; + + for_all(it) + if (ima(it) < ima(min)) + min = it; + + return min; + } + + template <class I> + mln_site(I) max (p_set<mln_site(I)>& components) + { + if (components.nsites() == 0) + return mln_site(I)(-1, -1); // FIXME + + typename p_set<mln_site(I)>::fwd_piter it(components); + mln_site(I) max = components[0]; + + for_all(it) + if (ima(it) > ima(max)) + max = it; + + return max; + } + + + // Actually, this structure is a tree, despite its confusing name. + template <class I, class N> + struct topo_wst + { + // Typedefs + // -------- + typedef mln_site(I) site; + typedef I image_t; + typedef N neighborhood_t; + + // Component tree management + // ------------------------- + private: + struct node { + mln_value(I) level; + int area; + int highest; + typedef mln_site(I) site; + // Can't modify the sites in a p_array + // p_array<mln_site(I)> children; + std::vector<site> children; + + void addChildren(const node& n) + { + // typename p_array<mln_site(I)>::fwd_piter it(n.children); + // for (it.start(); + // it.is_valid(); + // it.next()) + // children.append(it.to_site()); + for (unsigned i=0; i < n.children.size(); ++i) + children.push_back(n.children[i]); + } + + void addChild(const site p) + { + // children.append(n); + children.push_back(p); + } + }; // struct node + + public: + + mln_ch_value(I, site) Par_node; + mln_ch_value(I, site) Par_tree; + + mln_ch_value(I, int) Rnk_tree; + mln_ch_value(I, int) Rnk_node; + mln_ch_value(I, site) subtreeRoot; + mln_ch_value(I, node) nodes; + site Root; + + private: + + void MakeSet_tree(site x); + void MakeSet_node(site x); + site Find_tree(site x); + site Find_node(site x); + void BuildComponentTree(); + site MergeNode(site& node1, site& node2); + site Link_tree(site x, site y); + site Link_node(site x, site y); + node MakeNode(int level); + + private: + mln_ch_value(I, bool) isproc; + + // Ctor + public: + topo_wst(const Image<I>& i, + const Neighborhood<N>& n); + + public: + const I &ima; + const N &nbh; + + public: + void go(); + + private: + site min (p_set<site>& components); + site max (p_set<site>& components); + bool highest_fork (p_set<site>& components); + bool highest_fork (p_set<site>& components, site &r); + + // Optimized LCA Algorithm + + public: + site lca (site a, site b); + + private: + int *euler; + int *depth; + int ctree_size; + std::map<site, int, mln::util::ord<site> > pos; + site *repr; + + int *minim; + int **Minim; + + void compute_ctree_size (site p); + void build_euler_tour_rec(site p, int &position, int d); + void build_euler_tour(); + void build_minim (); + void removeOneSonNodes(site *p, mln_ch_value(I, site) &newPar_node); + void compressTree(); + }; + + /// Compute a toological watershed transform from \a tree. + template <class T> + typename T::image_t + topological(T &tree); + + + +# ifndef MLN_INCLUDE_ONLY + + // Ctor + template <class I, class N> + topo_wst<I, N>::topo_wst(const Image<I>& i, + const Neighborhood<N>& n) + : ima(exact(i)), nbh(exact(n)) + { + initialize(Par_node, i); + initialize(Par_tree, i); + initialize(Rnk_tree, i); + initialize(Rnk_node, i); + initialize(subtreeRoot, i); + initialize(nodes, i); + initialize(isproc, i); + } + + template <class I, class N> + void topo_wst<I, N>::MakeSet_tree(site x) + { + Par_tree(x) = x; + Rnk_tree(x) = 0; + } + + template <class I, class N> + void topo_wst<I, N>::MakeSet_node(site x) + { + Par_node(x) = x; + Rnk_node(x) = 0; + } + + template <class I, class N> + typename topo_wst<I, N>::site topo_wst<I, N>::Find_tree(site x) + { + if (Par_tree(x) != x) + Par_tree(x) = Find_tree(Par_tree(x)); + return Par_tree(x); + } + + template <class I, class N> + typename topo_wst<I, N>::site topo_wst<I, N>::Find_node(site x) + { + if (Par_node(x) != x) + Par_node(x) = Find_node(Par_node(x)); + return Par_node(x); + } + + template <class I, class N> + void topo_wst<I, N>::go() + { + BuildComponentTree(); + compressTree(); + + build_euler_tour(); + build_minim(); + } + + template <class I, class N> + void topo_wst<I, N>::BuildComponentTree() + { + // Init: + + // Sort the sites in increasing order + p_array<mln_site(I)> S; + S = data::sort_psites_increasing(ima); + + // Clear the marker map + data::fill(isproc, false); + for (int ip = 0; ip < int(S.nsites()); ++ip) + { + site p = S[ip]; + MakeSet_node(p); + MakeSet_tree(p); + // if (subtreeRoot.hold(p)) + subtreeRoot(p) = p; + // if (nodes.hold(p)) + nodes(p) = MakeNode(ima(p)); + } + + + + typename p_array<mln_site(I)>::fwd_piter ip(S); + for_all(ip) + { + site p = ip.to_site(); + + site curCanonicalElt = Find_tree(p); + site curNode = Find_node(subtreeRoot(curCanonicalElt)); + + // FIXME: Should be `n' instead of `q'. + mln_niter(N) q(nbh, ip); + for_all(q) + if (ima.has(q) and isproc(q) and ima(q) <= ima(p)) + { + site adjCanonicalElt = Find_tree(q); + site adjNode = Find_node(subtreeRoot(adjCanonicalElt)); + if (curNode != adjNode) + { + if (nodes(curNode).level == nodes(adjNode).level) + curNode = MergeNode(adjNode, curNode); + else + { + nodes(curNode).addChild(adjNode); + nodes(curNode).area += nodes(adjNode).area; + nodes(curNode).highest += nodes(adjNode).highest; + } + } + + curCanonicalElt = Link_tree(adjCanonicalElt, curCanonicalElt); + subtreeRoot(curCanonicalElt) = curNode; + } + isproc(p) = true; + } + // Pour garder une map de correspondance site <-> local_root + // for (int ip = 0; ip < int(S.size()); ++ip) + // { + // site p = S[ip]; + // M(p) = Find_node(p); + // } + + mln_piter(I) r(Par_node.domain()); + for_all(r) + Par_node(r) = Find_node(r); + + // Find the ``first'' site of ima, according to the forward + // traversal order. + mln_fwd_piter(I) rp(Par_node.domain());; + rp.start(); + + Root = subtreeRoot(Find_tree(Find_node(rp))); + } + + + template <class I, class N> + typename topo_wst<I, N>::site topo_wst<I, N>::MergeNode(site& node1, + site& node2) + { + site tmpNode = Link_node(node1, node2); + site tmpNode2; + if (tmpNode == node2) + { + nodes(node2).addChildren(nodes(node1)); + tmpNode2 = node1; + } + else + { + nodes(node1).addChildren(nodes(node2)); + tmpNode2 = node2; + } + nodes(tmpNode).area += nodes(tmpNode2).area; + nodes(tmpNode).highest += nodes(tmpNode2).highest; + return tmpNode; + } + + template <class I, class N> + typename topo_wst<I, N>::site topo_wst<I, N>::Link_tree(site x, site y) + { + if (Rnk_tree(x) > Rnk_tree(y)) + std::swap(x, y); + else + if (Rnk_tree(x) == Rnk_tree(y)) + Rnk_tree(y) += 1; + Par_tree(x) = y; + return y; + } + + template <class I, class N> + typename topo_wst<I, N>::site topo_wst<I, N>::Link_node(site x, site y) + { + if (Rnk_node(x) > Rnk_node(y)) + std::swap(x, y); + else + if (Rnk_node(x) == Rnk_node(y)) + Rnk_node(y) += 1; + Par_node(x) = y; + return y; + } + + template <class I, class N> + typename topo_wst<I, N>::node topo_wst<I, N>::MakeNode(int level) + { + node n; + n.level = level; + n.area = 1; + n.highest = level; + return n; + } + + + template <class I, class N> + bool topo_wst<I, N>::highest_fork (p_set<site>& components, site &r) + { + if (components.nsites() == 0) + { + std::cerr << "highest fork : empty set" << std::endl; + return false; + } + + site + m = min(components), + hfork = m; + + typename p_set<site>::fwd_piter it(components); + for_all(it) + { + // Can't remove the site from the set + if (it.to_site() == m) + continue; + + site c = lca(hfork, it.to_site()); + if (c != it.to_site()) + hfork = c; + } + + if (nodes(m).level == nodes(hfork).level) + return false; + + r = hfork; + return true; + } + + template <class I, class N> + bool topo_wst<I, N>::highest_fork (p_set<site>& components) { + site i; + return highest_fork(components, i); + } + + template <class I, class N> + void topo_wst<I, N>::compute_ctree_size (site p) + { + ctree_size += 1; + node& n = nodes(p); + + // typename p_array<mln_site(I)>::fwd_piter it(n.children); + // for_all(it) + // compute_ctree_size(it.to_site()); + + for (unsigned i=0; i < n.children.size(); ++i) + compute_ctree_size(n.children[i]); + } + + + template <class I, class N> + void topo_wst<I, N>::build_euler_tour_rec(site p, int &position, int d) + { + if (pos.find(p) == pos.end()) + pos[p] = position; + + repr[position] = p; + depth[position] = d; + euler[position] = pos[p]; + ++position; + node& n = nodes(p); + + // typename p_array<mln_site(I)>::fwd_piter it(n.children); + // for_all(it) + // { + // build_euler_tour_rec(it.to_site(), position, d+1); + // depth[position] = d; // Not optimized + // euler[position] = pos[p]; + // repr[position] = p; // Pas necessaire? + // ++position; + // } + + for (unsigned i=0; i < n.children.size(); ++i) + { + build_euler_tour_rec(n.children[i], position, d+1); + depth[position] = d; // Not optimized + euler[position] = pos[p]; + repr[position] = p; // Pas necessaire? + ++position; + } + } + + + template <class I, class N> + void topo_wst<I, N>::build_euler_tour () + { + ctree_size = 0; + compute_ctree_size(Root); + + int size = 2 * ctree_size - 1; + + // FIXME : free this + euler = new int[size]; + depth = new int[size]; + repr = new site[size]; + + int position = 0; + build_euler_tour_rec(Root, position, 0); + } + + + template <class I, class N> + void topo_wst<I, N>::build_minim () + { + // compute_tree_size(); // Already done in build_euler_tour + int size = 2 * ctree_size - 1; + int logn = (int)(ceil(log((double)(size))/log(2.0))); + // minim.reserve(size); + minim = new int [logn * size]; // FIXME : Think about freeing this + Minim = new int* [logn]; + + Minim[0] = minim; + + // Init + for (int i = 0; i < size - 1; ++i) + if (depth[euler[i]] < depth[euler[i+1]]) { + Minim[0][i] = i; + } else { + Minim[0][i] = i+1; + } + + // Minim[0][size - 1] = size - 1; + + int k; + for (int j = 1; j < logn; j++) { + k = 1 << (j - 1); + Minim[j] = &minim[j * size]; + for (int i = 0; i < size; i++) { + if ((i + (k << 1)) >= size) { + //Minim[j][i] = size - 1; + } + else { + if (depth[euler[Minim[j - 1][i]]] + <= depth[euler[Minim[j - 1][i + k]]]) + Minim[j][i] = Minim[j - 1][i]; + else + Minim[j][i] = Minim[j - 1][i + k]; + } + } + } + + } // void build_minim () + + + template <class I, class N> + typename topo_wst<I, N>::site topo_wst<I, N>::lca (site a, site b) + { + int + m = pos[a], + n = pos[b], + k; + + if (m == n) + return repr[m]; + + if (m > n) + { + k = n; + n = m; + m = k; + } + + k = (int)(log((double)(n - m))/log(2.)); + + if (depth[euler[Minim[k][m]]] < depth[euler[Minim[k][n - (1 << k)]]]) { + return repr[euler[Minim[k][m]]]; + } else { + return repr[euler[Minim[k][n - (1 << k)]]]; + } + } + + + template <class I, class N> + void topo_wst<I, N>::removeOneSonNodes(site *p, + mln_ch_value(I, site) &newPar_node) + { + node &n = nodes(*p); + + if (n.children.size() == 1) // this node has 1 son, delete it + { + n.area = -1; + newPar_node(*p) = n.children[0]; + *p = n.children[0]; + removeOneSonNodes(p, newPar_node); + } + else // there is more than one son, recursive call + { + for (unsigned i = 0; i < n.children.size(); ++i) + removeOneSonNodes(&(n.children[i]), newPar_node); + } + } + + + template <class I, class N> + void topo_wst<I, N>::compressTree() + { + mln_ch_value(I, site) newPar_node; + initialize(newPar_node, Par_node); + + // Remove the nodes with one son + removeOneSonNodes(&Root, newPar_node); + + // Update the references on deleted nodes + mln_piter(I) p(Par_node.domain()); + for_all(p) + while (nodes(Par_node(p)).area == -1) + Par_node(p) = newPar_node(Par_node(p)); + } + + template <class T> + bool w_constructible(T &tree, typename T::site p, typename T::site &r) + { + + typedef typename T::image_t I; + typedef typename T::neighborhood_t N; + + const I &ima = exact(tree.ima); + const N &nbh = exact(tree.nbh); + + // FIXME: Should be `n' instead of `q'. + mln_niter(N) q(nbh, p); + p_set<mln_site(I)> v; + + for_all(q) + // FIXME: Shouldn't it be: `ima.has(q)' instead of + // `ima.domain().has(q)'? + if (ima.domain().has(q) && ima(q) < ima(p)) + v.insert(tree.Par_node(q)); + + if (v.nsites() == 0) + return false; + + if (v.nsites() == 1) + { + r = v[0]; + return true; + } + + mln_site(I) c = min(ima, v); + mln_site(I) cmin = c; + + typename p_set<mln_site(I)>::fwd_piter it(v); + for_all(it) + { + // Can't remove the site from the set + if (it.to_site() == cmin) + continue; + + mln_site(I) c1 = tree.lca(c, it); + + if (c1 != it) + c = c1; + } + + if (tree.nodes(c).level >= ima(p)) + return false; + + r = c; + return true; + } + + template <class T> + bool w_constructible(T &tree, typename T::site p) { + typename T::site r; + return w_constructible(tree, p, r); + } + + + + template <class T> + typename T::image_t topological(T &tree) + { + + typedef typename T::image_t I; + typedef typename T::neighborhood_t N; + + I ima = exact(tree.ima); + const N &nbh = exact(tree.nbh); + + // Maxima components + mln_ch_value(I, bool) cmax; + initialize(cmax, ima); + + // Mark enqueued sites + mln_ch_value(I, bool) enqueued; + initialize(enqueued, ima); + + p_priority< mln_value(I), p_queue_fast<mln_site(I)> > l; + // p_queue < site > m; + std::queue<mln_site(I)> m; + mln_value(I) min = mln_min(mln_value(I)); + + std::cout << "Init" << std::endl; + + // Flag C-maxima + data::fill(cmax, false); + data::fill(enqueued, false); + + mln_piter(I) it(tree.Par_node.domain()); + for_all(it) + // if (nodes(Par_node(it.to_site())).children.nsites() == 0) + if (tree.nodes(tree.Par_node(it)).children.size() == 0) + { + cmax(it) = true; + m.push(it); + } + + while (!m.empty()) + { + // FIXME: Should be `n' instead of `q'. + mln_niter(N) q(nbh, m.front()); + // FIXME: Shouldn't it be: `cmax.has(q)' instead of + // `cmax.domain().has(q)'? + for_all(q) + if (cmax.domain().has(q) && !cmax(q) && !enqueued(q)) + { + enqueued(q) = true; + l.push(mln_max(mln_value(I)) - ima(q), q); + } + m.pop(); + } + + + std::cout << "end" << std::endl; + + // Main loop + while (!l.is_empty()) + { + mln_site(I) x = l.front(); + l.pop(); + enqueued(x) = false; + + mln_site(I) c; + bool is_w = w_constructible(tree, x, c); + + if (is_w) + { + ima(x) = tree.nodes(c).level; + tree.Par_node(x) = c; + + // if (nodes(c).children.nsites() == 0) + if (tree.nodes(c).children.size() == 0) + cmax(x) = true; + else + // if (nodes(c).children.nsites() > 1) + if (tree.nodes(c).children.size() == 1) + std::cerr << "ERREUR COMPOSANTE BRANCHE " + << tree.nodes(c).children.size() << std::endl; + + + // FIXME: Should be `n' instead of `q'. + mln_niter(N) q(nbh, x); + // FIXME: Shouldn't it be: `ima.has(q)' instead of + // `ima.domain().has(q)'? + for_all(q) + if (ima.domain().has(q) && !cmax(q) && !enqueued(q)) + { + enqueued(q) = true; + l.push(mln_max(mln_value(I)) - ima(q), q); // FIXME: + // Just + // invert + // the + // priority. + } + } + } // while(!l.empty()) + + for_all(it) + ima(it) = tree.nodes(tree.Par_node(it)).level; + + return ima; + } + +# endif // MLN_INCLUDE_ONLY + + + } // end of namespace mln::morpho::watershed + + } // end of namespace mln::morpho + +} // end of namespace mln + +#endif // ! MLN_MORPHO_WATERSHED_TOPOLOGICAL_HH diff --git a/milena/tests/morpho/watershed/Makefile.am b/milena/tests/morpho/watershed/Makefile.am index be39ac7..d384f5e 100644 --- a/milena/tests/morpho/watershed/Makefile.am +++ b/milena/tests/morpho/watershed/Makefile.am @@ -21,10 +21,12 @@ include $(top_srcdir)/milena/tests/tests.mk check_PROGRAMS = \ flooding \ - superpose + superpose \ + topological flooding_SOURCES = flooding.cc superpose_SOURCES = superpose.cc +topological_SOURCES = topological.cc TESTS = $(check_PROGRAMS) diff --git a/milena/tests/morpho/watershed/test_watershed_topo.cc b/milena/tests/morpho/watershed/test_watershed_topo.cc deleted file mode 100644 index db6d0f9..0000000 --- a/milena/tests/morpho/watershed/test_watershed_topo.cc +++ /dev/null @@ -1,42 +0,0 @@ -#include <mln/core/image/image2d.hh> -#include <mln/core/alias/window2d.hh> -#include <mln/core/alias/neighb2d.hh> -#include <mln/core/site_set/p_set.hh> - -#include <mln/value/int_u8.hh> -#include <mln/data/compare.hh> - -#include <mln/io/pgm/load.hh> -#include <mln/io/pgm/save.hh> - -#include <mln/morpho/topo_wst.hh> -#include <string> -#include <iostream> - -int main (int argc, const char * argv[]) -{ - using namespace mln; - using value::int_u8; - - typedef image2d<int_u8> image2dint; - - if (argc < 2) { - std::cerr << "usage: " << argv[0] << " in.pgm [other_files.pgm]" << std::endl; - return 1; - } - - for (int i = 1; i < argc; ++i) { - image2dint ima; - io::pgm::load(ima, argv[i]); - - morpho::topo_wst< image2d<int_u8>, neighb2d> n(ima, c4()); - - n.go(); - - std::string name(argv[i]); - name.erase(name.length() - 4); - - io::pgm::save(morpho::topo_watershed(n), name.append("_wsd.pgm")); - } - return 0; -} diff --git a/milena/mln/morpho/watershed/all.hh b/milena/tests/morpho/watershed/topological.cc similarity index 66% copy from milena/mln/morpho/watershed/all.hh copy to milena/tests/morpho/watershed/topological.cc index a07e75a..728b6c8 100644 --- a/milena/mln/morpho/watershed/all.hh +++ b/milena/tests/morpho/watershed/topological.cc @@ -23,42 +23,33 @@ // exception does not however invalidate any other reasons why the // executable file might be covered by the GNU General Public License. -#ifndef MLN_MORPHO_WATERSHED_ALL_HH -# define MLN_MORPHO_WATERSHED_ALL_HH - /// \file /// -/// File that includes all morphological watershed routines. +/// Exercise the Topological Watershed Transform (WST). +#include <mln/value/int_u8.hh> -namespace mln -{ +#include <mln/core/image/image2d.hh> +#include <mln/core/alias/window2d.hh> +#include <mln/core/alias/neighb2d.hh> - namespace morpho - { +#include <mln/morpho/watershed/topological.hh> - /// Namespace of morphological watershed routines. - namespace watershed - { +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> - /// Namespace of morphological watershed routines - /// implementations. - namespace watershed - { +#include "tests/data.hh" - /// Namespace of morphological watershed routines generic - /// implementations. - namespace generic - {} - } - } - } +int main() +{ + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> ima_t; + ima_t ima; + io::pgm::load(ima, MLN_IMG_DIR "/small.pgm"); + morpho::watershed::topo_wst< ima_t, neighb2d> n(ima, c4()); + n.go(); + io::pgm::save(morpho::watershed::topological(n), "topo_wst.pgm"); } - - -# include <mln/morpho/watershed/flooding.hh> -# include <mln/morpho/watershed/superpose.hh> - - -#endif // ! MLN_MORPHO_WATERSHED_ALL_HH -- 1.6.4.4