October 2009 - Olena-patches

[PATCH 7/7] Use mln::world::inter_pixel instead of cplx2d.hh.
by Roland Levillain 19 Oct '09

19 Oct '09

* apps/graph-morpho/morpho.hh: Here. * apps/graph-morpho/make_complex2d.hh (unmake_complex2d): Adjust as well. * apps/graph-morpho/cplx2d.hh: Remove symlink. --- milena/ChangeLog | 9 +++++++ milena/apps/graph-morpho/cplx2d.hh | 1 - milena/apps/graph-morpho/make_complex2d.hh | 5 ++- milena/apps/graph-morpho/morpho.hh | 34 +++++++++++++++------------ 4 files changed, 31 insertions(+), 18 deletions(-) delete mode 120000 milena/apps/graph-morpho/cplx2d.hh diff --git a/milena/ChangeLog b/milena/ChangeLog index c3cf6b5..6c649f8 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,14 @@ 2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + Use mln::world::inter_pixel instead of cplx2d.hh. + + * apps/graph-morpho/morpho.hh: Here. + * apps/graph-morpho/make_complex2d.hh (unmake_complex2d): + Adjust as well. + * apps/graph-morpho/cplx2d.hh: Remove symlink. + +2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + Rename 2D inter-pixel neighborhoods. * mln/world/inter_pixel/neighb2d.hh diff --git a/milena/apps/graph-morpho/cplx2d.hh b/milena/apps/graph-morpho/cplx2d.hh deleted file mode 120000 index c08d52c..0000000 --- a/milena/apps/graph-morpho/cplx2d.hh +++ /dev/null @@ -1 +0,0 @@ -../../sandbox/theo/esiee/laurent/ismm09/cplx2d.hh \ No newline at end of file diff --git a/milena/apps/graph-morpho/make_complex2d.hh b/milena/apps/graph-morpho/make_complex2d.hh index a79c798..e1c0c20 100644 --- a/milena/apps/graph-morpho/make_complex2d.hh +++ b/milena/apps/graph-morpho/make_complex2d.hh @@ -30,6 +30,7 @@ /// \brief Cubical 2-complex creation from a 2D image of pixels. # include <mln/core/image/image2d.hh> +# include <mln/world/inter_pixel/dim2/is_pixel.hh> /** \brief Create an binary 2D image representing a cubical 2-complex by doubling the resolution of the image \a input images to insert @@ -109,8 +110,8 @@ unmake_complex2d(const mln::Image& input_) mln_precondition(input.ncols() % 2 == 1); // Create a (morpher) image of the pixels of INPUT. - typedef image_if< const I, cplx2d::predicate_t > J; - J input_pixels = input | cplx2d::is_pixel; + typedef image_if< const I, world::inter_pixel::dim2::is_pixel > J; + J input_pixels = input | world::inter_pixel::dim2::is_pixel(); /* FIXME: The construction of OUTPUT is obvioulsy not generic, since it expects I's domain to provide the interface of an mln::box2d. diff --git a/milena/apps/graph-morpho/morpho.hh b/milena/apps/graph-morpho/morpho.hh index d0c1a1d..9a786da 100644 --- a/milena/apps/graph-morpho/morpho.hh +++ b/milena/apps/graph-morpho/morpho.hh @@ -27,7 +27,7 @@ # define APPS_GRAPH_MORPHO_MORPHO_HH /** \file apps/graph-morpho/morpho.hh - \brief Morphological operators on graphs (1-complexes). + \brief Morphological operators on graphs. Reference: @@ -39,20 +39,22 @@ # include <mln/core/alias/complex_image.hh> # include <mln/core/image/image2d.hh> +# include <mln/core/image/dmorph/image_if.hh> + # include <mln/core/routine/duplicate.hh> # include <mln/core/site_set/p_n_faces_piter.hh> # include <mln/core/image/complex_neighborhoods.hh> # include <mln/core/image/complex_neighborhood_piter.hh> +# include <mln/world/inter_pixel/dim2/is_pixel.hh> +# include <mln/world/inter_pixel/dim2/is_edge.hh> +# include <mln/world/inter_pixel/neighb2d.hh> + # include <mln/data/paste.hh> # include "apps/graph-morpho/image_if_large.hh" -// FIXME: Careful, cplx2d.hh is a symlink to a file in Théo's sandbox, -// and might be changed. -# include "apps/graph-morpho/cplx2d.hh" - // FIXME: Instead of providing several implementation, move specific // parts (neighborhoods, etc.) to a graph_traits class, and write // generic version of combine, dilation_e2v, erosion_v2e, etc. @@ -107,8 +109,10 @@ namespace impl mln::image2d<T> output; mln::initialize(output, vertices); mln::data::fill(output, false); - mln::data::paste(vertices | mln::cplx2d::is_pixel, output); - mln::data::paste(edges | mln::cplx2d::is_edge, output); + mln::data::paste(vertices | mln::world::inter_pixel::dim2::is_pixel(), + output); + mln::data::paste(edges | mln::world::inter_pixel::dim2::is_edge(), + output); return output; } } @@ -351,8 +355,8 @@ namespace impl { mln::image2d<T> output(input.domain()); mln::data::fill(output, false); - mln::data::paste(dilation(input || mln::cplx2d::is_pixel, - mln::cplx2d::p2e()), + mln::data::paste(dilation(input || mln::world::inter_pixel::dim2::is_pixel(), + mln::world::inter_pixel::v2e()), output); return output; } @@ -366,8 +370,8 @@ namespace impl { mln::image2d<T> output(input.domain()); mln::data::fill(output, false); - mln::data::paste(erosion(input || mln::cplx2d::is_edge, - mln::cplx2d::e2p()), + mln::data::paste(erosion(input || mln::world::inter_pixel::dim2::is_edge(), + mln::world::inter_pixel::e2v()), output); return output; } @@ -381,8 +385,8 @@ namespace impl { mln::image2d<T> output(input.domain()); mln::data::fill(output, false); - mln::data::paste(erosion(input || mln::cplx2d::is_pixel, - mln::cplx2d::p2e()), + mln::data::paste(erosion(input || mln::world::inter_pixel::dim2::is_pixel(), + mln::world::inter_pixel::v2e()), output); return output; } @@ -396,8 +400,8 @@ namespace impl { mln::image2d<T> output(input.domain()); mln::data::fill(output, false); - mln::data::paste(dilation(input || mln::cplx2d::is_edge, - mln::cplx2d::e2p()), + mln::data::paste(dilation(input || mln::world::inter_pixel::dim2::is_edge(), + mln::world::inter_pixel::e2v()), output); return output; } -- 1.6.4.4

1 0

[PATCH 2/3] Improve the integration of the topological watershed transform.
by Roland Levillain 19 Oct '09

19 Oct '09

* 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 <roland(a)lrde.epita.fr> + 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 <roland(a)lrde.epita.fr> + 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 &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, - 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, - 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 &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, + 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, + 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

1 0

[PATCH 1/3] Import Alexandre's implementation of the topological watershed.
by Roland Levillain 19 Oct '09

19 Oct '09

* mln/morpho/watershed/topo_wst.hh: New file, imported from sandbox/abraham/mln/morpho/topo_wst.hh. * tests/morpho/watershed/test_watershed_topo.cc: New file, imported from sandbox/abraham/tests/morpho/test_watershed_topo.cc. --- milena/ChangeLog | 9 + milena/mln/morpho/watershed/topo_wst.hh | 768 ++++++++++++++++++++ .../tests/morpho/watershed/test_watershed_topo.cc | 42 ++ 3 files changed, 819 insertions(+), 0 deletions(-) create mode 100644 milena/mln/morpho/watershed/topo_wst.hh create mode 100644 milena/tests/morpho/watershed/test_watershed_topo.cc diff --git a/milena/ChangeLog b/milena/ChangeLog index 19c3edd..9789d98 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,3 +1,12 @@ +2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + + Import Alexandre's implementation of the topological watershed. + + * mln/morpho/watershed/topo_wst.hh: New file, imported from + sandbox/abraham/mln/morpho/topo_wst.hh. + * tests/morpho/watershed/test_watershed_topo.cc: New file, + imported from sandbox/abraham/tests/morpho/test_watershed_topo.cc. + 2009-09-22 Guillaume Lazzara <lazzara(a)lrde.epita.fr> Add new upscaling algorithms. diff --git a/milena/mln/morpho/watershed/topo_wst.hh b/milena/mln/morpho/watershed/topo_wst.hh new file mode 100644 index 0000000..7e9fb64 --- /dev/null +++ b/milena/mln/morpho/watershed/topo_wst.hh @@ -0,0 +1,768 @@ +// 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 &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, + 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, + 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/tests/morpho/watershed/test_watershed_topo.cc b/milena/tests/morpho/watershed/test_watershed_topo.cc new file mode 100644 index 0000000..db6d0f9 --- /dev/null +++ b/milena/tests/morpho/watershed/test_watershed_topo.cc @@ -0,0 +1,42 @@ +#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; +} -- 1.6.4.4

1 0

[PATCH 3/6] configure.ac: Configure milena/apps/constrained-connectivity/Makefile.
by Roland Levillain 19 Oct '09

19 Oct '09

--- ChangeLog | 5 +++++ configure.ac | 1 + 2 files changed, 6 insertions(+), 0 deletions(-) diff --git a/ChangeLog b/ChangeLog index 10d4808..3e4343e 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,8 @@ +2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + + * configure.ac: Configure + milena/apps/constrained-connectivity/Makefile. + 2009-09-22 Guillaume Lazzara <lazzara(a)lrde.epita.fr> * configure.ac: Configure milena/tests/upscaling. diff --git a/configure.ac b/configure.ac index 4aa31d4..daea625 100644 --- a/configure.ac +++ b/configure.ac @@ -406,6 +406,7 @@ AC_CONFIG_FILES([ milena/apps/Makefile milena/apps/mesh-segm-skel/Makefile milena/apps/graph-morpho/Makefile + milena/apps/constrained-connectivity/Makefile ]) # Configure tests. -- 1.6.4.4

1 0

[PATCH 5/6] configure.ac: Configure milena/apps/constrained-connectivity/test-constrained-connectivity.
by Roland Levillain 19 Oct '09

19 Oct '09

--- ChangeLog | 5 +++++ configure.ac | 4 +++- 2 files changed, 8 insertions(+), 1 deletions(-) diff --git a/ChangeLog b/ChangeLog index 3e4343e..cfca396 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,6 +1,11 @@ 2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> * configure.ac: Configure + milena/apps/constrained-connectivity/test-constrained-connectivity. + +2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + + * configure.ac: Configure milena/apps/constrained-connectivity/Makefile. 2009-09-22 Guillaume Lazzara <lazzara(a)lrde.epita.fr> diff --git a/configure.ac b/configure.ac index daea625..e325477 100644 --- a/configure.ac +++ b/configure.ac @@ -422,10 +422,12 @@ AC_CONFIG_FILES([milena/apps/mesh-segm-skel/test-mesh-complex-segm], [chmod +x milena/apps/mesh-segm-skel/test-mesh-complex-segm]) AC_CONFIG_FILES([milena/apps/mesh-segm-skel/test-mesh-complex-max-curv-segm], [chmod +x milena/apps/mesh-segm-skel/test-mesh-complex-max-curv-segm]) - AC_CONFIG_FILES([milena/apps/mesh-segm-skel/test-mesh-complex-skel], [chmod +x milena/apps/mesh-segm-skel/test-mesh-complex-skel]) +AC_CONFIG_FILES([milena/apps/constrained-connectivity/test-constrained-connectivity], + [chmod +x milena/apps/constrained-connectivity/test-constrained-connectivity]) + # Flags for apps. AC_ARG_VAR([APPS_CXXFLAGS]) # We want fast binaries for apps. -- 1.6.4.4

1 0

[PATCH 1/6] Use Milena's topological WST instead of the one from the sandboxes.
by Roland Levillain 19 Oct '09

19 Oct '09

* roland/constrained-connectivity.cc: Here. * roland/alexandre: Remove symlink. --- milena/sandbox/ChangeLog | 7 +++++++ milena/sandbox/roland/alexandre | 1 - milena/sandbox/roland/constrained-connectivity.cc | 7 +++---- 3 files changed, 10 insertions(+), 5 deletions(-) delete mode 120000 milena/sandbox/roland/alexandre diff --git a/milena/sandbox/ChangeLog b/milena/sandbox/ChangeLog index d087c07..7f9554a 100644 --- a/milena/sandbox/ChangeLog +++ b/milena/sandbox/ChangeLog @@ -1,5 +1,12 @@ 2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + Use Milena's topological WST instead of the one from the sandboxes. + + * roland/constrained-connectivity.cc: Here. + * roland/alexandre: Remove symlink. + +2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + Catch up with new 2D inter-pixel neighborhoods names. * roland/constrained-connectivity.cc (main): s/e2c/e2v/. diff --git a/milena/sandbox/roland/alexandre b/milena/sandbox/roland/alexandre deleted file mode 120000 index 5edb140..0000000 --- a/milena/sandbox/roland/alexandre +++ /dev/null @@ -1 +0,0 @@ -../abraham/mln/morpho/ \ No newline at end of file diff --git a/milena/sandbox/roland/constrained-connectivity.cc b/milena/sandbox/roland/constrained-connectivity.cc index 7dee1d6..59009b7 100644 --- a/milena/sandbox/roland/constrained-connectivity.cc +++ b/milena/sandbox/roland/constrained-connectivity.cc @@ -76,8 +76,7 @@ #include <mln/world/inter_pixel/compute.hh> #include <mln/world/inter_pixel/neighb2d.hh> -// From Alexandre Abraham's sandbox. -#include <alexandre/topo_wst.hh> +#include <mln/morpho/watershed/topological.hh> #include <mln/morpho/tree/compute_attribute_image.hh> #include <mln/accu/stat/min.hh> @@ -129,10 +128,10 @@ int main(int argc, char* argv[]) debug::println("g:", g); // Compute a topological watershed transform on this gradient. - typedef morpho::topo_wst<g_t, world::inter_pixel::dbl_neighb2d> tree_t; + typedef morpho::watershed::topo_wst<g_t, world::inter_pixel::dbl_neighb2d> tree_t; tree_t tree(g, world::inter_pixel::e2e()); tree.go(); - mln_VAR(w, morpho::topo_watershed(tree)); + mln_VAR(w, morpho::watershed::topological(tree)); debug::println("w:", w); // Computing the set of values of W. -- 1.6.4.4

1 0

[PATCH 2/6] Move the experiment on constrained connectivity out of my sandbox.
by Roland Levillain 19 Oct '09

19 Oct '09

* apps/constrained-connectivity/: New directory. * sandbox/roland/constrained-connectivity.cc, * sandbox/roland/soille.pgm: Move files... * apps/constrained-connectivity/constrained-connectivity.cc, * apps/constrained-connectivity/soille.pgm: ...here. * sandbox/roland/Makefile: Adjust. * apps/constrained-connectivity/Makefile.am: New. * apps/Makefile.am (SUBDIRS): Add constrained-connectivity. --- milena/ChangeLog | 15 +++++++++++++++ milena/apps/Makefile.am | 2 +- .../{ => constrained-connectivity}/Makefile.am | 15 +++++++++++---- .../constrained-connectivity.cc | 0 .../constrained-connectivity}/soille.pgm | 0 milena/sandbox/roland/Makefile | 10 ++++------ 6 files changed, 31 insertions(+), 11 deletions(-) copy milena/apps/{ => constrained-connectivity}/Makefile.am (58%) rename milena/{sandbox/roland => apps/constrained-connectivity}/constrained-connectivity.cc (100%) rename milena/{sandbox/roland => apps/constrained-connectivity}/soille.pgm (100%) diff --git a/milena/ChangeLog b/milena/ChangeLog index 6c649f8..17cd2c7 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,20 @@ 2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + Move the experiment on constrained connectivity out of my sandbox. + + * apps/constrained-connectivity/: New directory. + * sandbox/roland/constrained-connectivity.cc, + * sandbox/roland/soille.pgm: + Move files... + * apps/constrained-connectivity/constrained-connectivity.cc, + * apps/constrained-connectivity/soille.pgm: + ...here. + * sandbox/roland/Makefile: Adjust. + * apps/constrained-connectivity/Makefile.am: New. + * apps/Makefile.am (SUBDIRS): Add constrained-connectivity. + +2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + Use mln::world::inter_pixel instead of cplx2d.hh. * apps/graph-morpho/morpho.hh: Here. diff --git a/milena/apps/Makefile.am b/milena/apps/Makefile.am index 2150657..dc39cff 100644 --- a/milena/apps/Makefile.am +++ b/milena/apps/Makefile.am @@ -17,4 +17,4 @@ ## Process this file through Automake to produce Makefile.in. -SUBDIRS = mesh-segm-skel graph-morpho +SUBDIRS = mesh-segm-skel graph-morpho constrained-connectivity diff --git a/milena/apps/Makefile.am b/milena/apps/constrained-connectivity/Makefile.am similarity index 58% copy from milena/apps/Makefile.am copy to milena/apps/constrained-connectivity/Makefile.am index 2150657..e27a18e 100644 --- a/milena/apps/Makefile.am +++ b/milena/apps/constrained-connectivity/Makefile.am @@ -1,4 +1,4 @@ -# Copyright (C) 2008, 2009 EPITA Research and Development Laboratory (LRDE). +# Copyright (C) 2009 EPITA Research and Development Laboratory (LRDE). # # This file is part of Olena. # @@ -13,8 +13,15 @@ # # You should have received a copy of the GNU General Public License # along with Olena. If not, see <http://www.gnu.org/licenses/>. -# -## Process this file through Automake to produce Makefile.in. +# Find Milena headers. +AM_CPPFLAGS = -I$(top_srcdir)/milena -I$(top_builddir)/milena +# Produce fast code. +APPS_CXXFLAGS = @APPS_CXXFLAGS@ +AM_CXXFLAGS = $(APPS_CXXFLAGS) + +noinst_PROGRAMS = constrained-connectivity +constrained_connectivity_SOURCES = constrained-connectivity.cc -SUBDIRS = mesh-segm-skel graph-morpho +# The sample image form Pierre Soille's PAMI 2008 article. +EXTRA_DIST = soille.pgm diff --git a/milena/sandbox/roland/constrained-connectivity.cc b/milena/apps/constrained-connectivity/constrained-connectivity.cc similarity index 100% rename from milena/sandbox/roland/constrained-connectivity.cc rename to milena/apps/constrained-connectivity/constrained-connectivity.cc diff --git a/milena/sandbox/roland/soille.pgm b/milena/apps/constrained-connectivity/soille.pgm similarity index 100% rename from milena/sandbox/roland/soille.pgm rename to milena/apps/constrained-connectivity/soille.pgm diff --git a/milena/sandbox/roland/Makefile b/milena/sandbox/roland/Makefile index d95c7a5..7ef48bf 100644 --- a/milena/sandbox/roland/Makefile +++ b/milena/sandbox/roland/Makefile @@ -6,22 +6,20 @@ img_dir = $(milena_dir)/img CPPFLAGS = -I. -I$(milena_dir) CXXFLAGS = -ggdb -Wall -Werror -PROGRAMS = double min-max constrained-connectivity +PROGRAMS = double min-max all: $(PROGRAMS) -check: check-double check-min-max check-constrained-connectivity +check: check-double check-min-max check-double: double ./$< check-min-max: min-max ./$< -check-constrained-connectivity: constrained-connectivity soille.pgm - ./$< soille.pgm -CLEANFILES = double min-max constrained-connectivity +CLEANFILES = double min-max clean: rm -f $(CLEANFILES) .PHONY: all check clean -.PHONY: check-double check-min-max check-constrained-connectivity +.PHONY: check-double check-min-max -- 1.6.4.4

1 0

[PATCH 2/7] Fix documentation in mln/world/inter_pixel/neighb2d.hh.
by Roland Levillain 19 Oct '09

19 Oct '09

* mln/world/inter_pixel/neighb2d.hh (mln::world::inter_pixel::e2c, mln::world::inter_pixel::e2e): Here. Wrap long lines. --- milena/ChangeLog | 9 +++++++++ milena/mln/world/inter_pixel/neighb2d.hh | 10 ++++++---- 2 files changed, 15 insertions(+), 4 deletions(-) diff --git a/milena/ChangeLog b/milena/ChangeLog index eb1d50a..db4d604 100644 --- a/milena/ChangeLog +++ b/milena/ChangeLog @@ -1,5 +1,14 @@ 2009-09-23 Roland Levillain <roland(a)lrde.epita.fr> + Fix documentation in mln/world/inter_pixel/neighb2d.hh. + + * mln/world/inter_pixel/neighb2d.hh + (mln::world::inter_pixel::e2c, mln::world::inter_pixel::e2e): + Here. + Wrap long lines. + +2009-09-23 Roland Levillain <roland(a)lrde.epita.fr> + Include headers missing from mln/world/inter_pixel/compute.hh. * mln/world/inter_pixel/compute.hh: Include diff --git a/milena/mln/world/inter_pixel/neighb2d.hh b/milena/mln/world/inter_pixel/neighb2d.hh index c39bd2f..21233ef 100644 --- a/milena/mln/world/inter_pixel/neighb2d.hh +++ b/milena/mln/world/inter_pixel/neighb2d.hh @@ -47,10 +47,10 @@ namespace mln /// Double neighborhood used for inter-pixel images. typedef neighb< win::multiple<window2d, dim2::is_row_odd> > dbl_neighb2d; - /// C4 neighborhood on pixels centered on an edge. + // Edge to neighborhood pixels. const dbl_neighb2d& e2c(); - /// C8 neighborhood on edges centered on an edge. + // Edge to neighboring edges. const dbl_neighb2d& e2e(); @@ -66,7 +66,8 @@ namespace mln 1, 0, 1, 0, 0, 0 }; - static dbl_neighb2d nbh = make::double_neighb2d(dim2::is_row_odd(), e2c_h, e2c_v); + static dbl_neighb2d nbh = + make::double_neighb2d(dim2::is_row_odd(), e2c_h, e2c_v); return nbh; } @@ -86,7 +87,8 @@ namespace mln 0, 1, 0, 1, 0, 0, 0, 0, 0, 0 }; - static dbl_neighb2d nbh = make::double_neighb2d(dim2::is_row_odd(), e2e_h, e2e_v); + static dbl_neighb2d nbh = + make::double_neighb2d(dim2::is_row_odd(), e2e_h, e2e_v); return nbh; } -- 1.6.4.4

1 0

[PATCH 3/3] Annotate Alexandre's sandbox.
by Roland Levillain 19 Oct '09

19 Oct '09

* abraham/README: New. --- milena/sandbox/ChangeLog | 6 ++++++ milena/sandbox/abraham/README | 11 +++++++++++ 2 files changed, 17 insertions(+), 0 deletions(-) create mode 100644 milena/sandbox/abraham/README diff --git a/milena/sandbox/ChangeLog b/milena/sandbox/ChangeLog index cb4de1e..a069dc9 100644 --- a/milena/sandbox/ChangeLog +++ b/milena/sandbox/ChangeLog @@ -1,3 +1,9 @@ +2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + + Annotate Alexandre's sandbox. + + * abraham/README: New. + 2009-09-22 Guillaume Lazzara <lazzara(a)lrde.epita.fr> Fix invalid output domain in hq*x algorithms. diff --git a/milena/sandbox/abraham/README b/milena/sandbox/abraham/README new file mode 100644 index 0000000..a7372d7 --- /dev/null +++ b/milena/sandbox/abraham/README @@ -0,0 +1,11 @@ +These files have been integrated into Milena: + + mln/morpho/topo_wst.hh + -> ../../mln/morpho/watershed/topological.hh + + tests/morpho/test_watershed_topo.cc + -> ../../tests/morpho/watershed/topological.cc + +The versions contained in this sandbox should be considered obsolete +and should not be used nor modified anymore. Actually, removing them +might be a good option. -- 1.6.4.4

1 0

[PATCH 4/7] Stop using cplx2d.hh in the constrained connectivity prototype.
by Roland Levillain 19 Oct '09

19 Oct '09

* roland/constrained-connectivity.cc: Here. Use Milena's mln::world::inter_pixel instead. Adjust. * roland/theo: Remove symlink. --- milena/sandbox/ChangeLog | 9 +++++ milena/sandbox/roland/constrained-connectivity.cc | 40 +++++++++++++-------- milena/sandbox/roland/theo | 1 - 3 files changed, 34 insertions(+), 16 deletions(-) delete mode 120000 milena/sandbox/roland/theo diff --git a/milena/sandbox/ChangeLog b/milena/sandbox/ChangeLog index a069dc9..52e6334 100644 --- a/milena/sandbox/ChangeLog +++ b/milena/sandbox/ChangeLog @@ -1,5 +1,14 @@ 2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + Stop using cplx2d.hh in the constrained connectivity prototype. + + * roland/constrained-connectivity.cc: Here. + Use Milena's mln::world::inter_pixel instead. + Adjust. + * roland/theo: Remove symlink. + +2009-09-24 Roland Levillain <roland(a)lrde.epita.fr> + Annotate Alexandre's sandbox. * abraham/README: New. diff --git a/milena/sandbox/roland/constrained-connectivity.cc b/milena/sandbox/roland/constrained-connectivity.cc index af5e49a..fb14d2c 100644 --- a/milena/sandbox/roland/constrained-connectivity.cc +++ b/milena/sandbox/roland/constrained-connectivity.cc @@ -45,7 +45,7 @@ \li Assign values of a gradient to the edges (e.g., |v1 - v2|). \li Compute a topological watershed transform on this gradient. \li Thresholding the values of this watershed gives - alpha-connected component (alpha-cc) described by Pierre Soille. + alpha-connected components (alpha-cc) described by Pierre Soille. Part 2. \li For each edge of the watershed (line graph-based) image, @@ -58,7 +58,7 @@ values from the image of max values. \li Thresholding the watershed image (using a parameter alpha) and the height image (using a parameter omega) gives - the (alpha, omega)-connected component (alpha-omega cc) + the (alpha, omega)-connected components (alpha-omega cc) described by Pierre Soille. */ #include <cstdio> @@ -71,8 +71,10 @@ #include <mln/pw/all.hh> -// From Théo's sandbox. -#include <theo/cplx2d.hh> +#include <mln/fun/vv2v/diff_abs.hh> +#include <mln/world/inter_pixel/immerse.hh> +#include <mln/world/inter_pixel/compute.hh> +#include <mln/world/inter_pixel/neighb2d.hh> // From Alexandre Abraham's sandbox. #include <alexandre/topo_wst.hh> @@ -84,6 +86,8 @@ #include <mln/io/pgm/load.hh> #include <mln/debug/println.hh> +#include <mln/core/var.hh> + int main(int argc, char* argv[]) { @@ -101,8 +105,8 @@ int main(int argc, char* argv[]) io::pgm::load(input, argv[1]); // Double its resolution. - image2d<int_u8> f(input.nrows() * 2, input.ncols() * 2); - mln_piter_(image2d<int_u8>) p_ima(f.domain()); + image2d<int_u8> f_(input.nrows() * 2, input.ncols() * 2); + mln_piter_(image2d<int_u8>) p_ima(f_.domain()); for_all(p_ima) { /* This conversion from a ``piter'' type to point2d is required, @@ -111,17 +115,22 @@ int main(int argc, char* argv[]) etc.). */ point2d p_ima_ = p_ima; point2d p_f(p_ima_.row() / 2, p_ima_.col() / 2); - f(p_ima) = input(p_f); + f_(p_ima) = input(p_f); } - debug::println(f); + debug::println("f_:", f_); + + image_if<image2d<int_u8>, world::inter_pixel::is_pixel> f = + world::inter_pixel::immerse(f_); + debug::println("f:", f); // Compute the associated line graph gradient. - mln_VAR(g, cplx2d::f_to_g(f) ); + mln_VAR(g, world::inter_pixel::compute (f, fun::vv2v::diff_abs<int_u8>())); + debug::println("g:", g); // Compute a topological watershed transform on this gradient. - typedef morpho::topo_wst<g_t, cplx2d::dbl_neighb2d> tree_t; - tree_t tree(g, cplx2d::e2e()); + typedef morpho::topo_wst<g_t, world::inter_pixel::dbl_neighb2d> tree_t; + tree_t tree(g, world::inter_pixel::e2e()); tree.go(); mln_VAR(w, morpho::topo_watershed(tree)); debug::println("w:", w); @@ -158,7 +167,7 @@ int main(int argc, char* argv[]) morpho::tree::data. */ typedef p_array<tree_t::site> sites_t; sites_t sites = data::sort_psites_decreasing(w); - morpho::tree::data<w_t, sites_t> t(w, sites, cplx2d::e2e()); + morpho::tree::data<w_t, sites_t> t(w, sites, world::inter_pixel::e2e()); // Create initial images for min and max values on sites (not components). mln_ch_value_(w_t, accu::stat::min<int_u8>) init_min_val; @@ -174,15 +183,16 @@ int main(int argc, char* argv[]) vertices/pixels). We have to convert the coordinates of V to its equivalent in F's domain to get the values on vertices. */ mln_piter_(w_t) e(w.domain()); - mln_niter_(cplx2d::dbl_neighb2d) v_g(cplx2d::e2p(), e); + mln_niter_(world::inter_pixel::dbl_neighb2d) + v_g(world::inter_pixel::e2c(), e); for_all(e) for_all(v_g) { // Same remark as above avour piter to point2d conversions. point2d v_g_ = v_g; point2d v_f(v_g_.row() / 2, v_g_.col() / 2); - init_min_val(e).take(f(v_f)); - init_max_val(e).take(f(v_f)); + init_min_val(e).take(f_(v_f)); + init_max_val(e).take(f_(v_f)); } // Attribute images of min and max values on components. accu::stat::min<int_u8> min_accu; diff --git a/milena/sandbox/roland/theo b/milena/sandbox/roland/theo deleted file mode 120000 index 7612684..0000000 --- a/milena/sandbox/roland/theo +++ /dev/null @@ -1 +0,0 @@ -../theo/esiee/laurent/ismm09 \ No newline at end of file -- 1.6.4.4

1 0