[Olena-patches] 3154: Make the topological Watershed a standalon function.

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Alexandre Abraham &lt;abraham(a)lrde.epita.fr&gt; Make the topological Watershed a standalon function. * abraham/tests/morpho/test_watershed_topo.cc: Updated. * abraham/tests/morpho/levelings.cc: Updated. * abraham/tests/morpho/Makefile: Add dependencies. * abraham/tests/fun/meta/red.cc: Updated. * abraham/mln/morpho/vector_median.hh: Add 2009 year. * abraham/mln/morpho/topo_wst.hh: Extract watershed method from the big structure. mln/morpho/topo_wst.hh | 419 +++++++++++++++++------------------- mln/morpho/vector_median.hh | 2 tests/fun/meta/red.cc | 3 tests/morpho/Makefile | 4 tests/morpho/levelings.cc | 41 +-- tests/morpho/test_watershed_topo.cc | 7 6 files changed, 229 insertions(+), 247 deletions(-) Index: abraham/tests/morpho/test_watershed_topo.cc --- abraham/tests/morpho/test_watershed_topo.cc (revision 3153) +++ abraham/tests/morpho/test_watershed_topo.cc (working copy) @@ -25,9 +25,7 @@ return 1; } - for (int i = 1; i < argc; ++i) - { - + for (int i = 1; i < argc; ++i) { image2dint ima; io::pgm::load(ima, argv[i]); @@ -38,8 +36,7 @@ std::string name(argv[i]); name.erase(name.length() - 4); - io::pgm::save(n.pima, name.append("_wsd.pgm")); - + io::pgm::save(morpho::topo_watershed(n), name.append("_wsd.pgm")); } return 0; } Index: abraham/tests/morpho/levelings.cc --- abraham/tests/morpho/levelings.cc (revision 3153) +++ abraham/tests/morpho/levelings.cc (working copy) @@ -3,6 +3,7 @@ #include <mln/io/ppm/load.hh> #include <mln/value/int_u8.hh> #include <mln/morpho/vector_median.hh> +#include <mln/morpho/vmt.hh> #include <mln/morpho/autarkical_leveling.hh> #include <mln/fun/meta/red.hh> #include <mln/value/rgb8.hh> @@ -10,31 +11,6 @@ #include <mln/core/image/fun_image.hh> #include <mln/core/image/violent_cast_image.hh> -namespace mln -{ - struct red_only - { - template <unsigned n> - value::rgb<n> max_() const - { - return value::rgb<n>(mln_max(typename value::rgb<n>::red_t), 0, 0); - } - - value::rgb8 max () const - { - return value::rgb8(255, 0, 0); - } - - template <unsigned n> - bool less(const value::rgb<n>& a, const value::rgb<n>& b) - { - return a.red() < b.red(); - } - }; - -} - - int main (int argc, const char * argv[]) { using namespace mln; @@ -45,8 +21,6 @@ return 1; } - typedef value::mixin<value::rgb8,red_only> Rgb; - for (int i = 1; i < argc; ++i) { border::thickness = 2; @@ -57,9 +31,18 @@ std::string name(argv[i]); name.erase(name.length() - 4); - win::rectangle2d rect(8, 8); + win::rectangle2d smallrect(13, 13); + win::rectangle2d bigrect(25, 25); + + image2d<rgb8> ima_vm = morpho::vector_median(ima, bigrect); + + io::ppm::save(ima_vm, name.append("_vector_medianed.ppm")); + + name = argv[i]; + name.erase(name.length() - 4); + - io::ppm::save(morpho::autarkical_leveling(ima, morpho::vector_median(ima, rect), rect), name.append("_autarkical_leveled.ppm")); + io::ppm::save(morpho::autarkical_leveling(ima, ima_vm, smallrect), name.append("_autarkical_leveled.ppm")); } return 0; Index: abraham/tests/morpho/Makefile --- abraham/tests/morpho/Makefile (revision 3153) +++ abraham/tests/morpho/Makefile (working copy) @@ -16,6 +16,7 @@ .PHONY: clean distclean ## Depedencies -test_watershed_topo.cc: ../../mln/morpho/basic_najman.hh +test_watershed_topo_new.cc: ../../mln/morpho/topo_wst_new.hh +test_watershed_topo.cc: ../../mln/morpho/topo_wst.hh test_watershed.cc: ../../mln/morpho/basic_najman.hh -test_componenet_tree.cc: ../../mln/morpho/basic_najman.hh \ No newline at end of file +test_component_tree.cc: ../../mln/morpho/basic_najman.hh Index: abraham/tests/fun/meta/red.cc --- abraham/tests/fun/meta/red.cc (revision 3153) +++ abraham/tests/fun/meta/red.cc (working copy) @@ -57,7 +57,10 @@ mln_piter_(mln::image2d<C>) p(i.domain()); for_all (p) + { + out(p) = 2; std::cout << out(p) << std::endl; + } } Index: abraham/mln/morpho/vector_median.hh --- abraham/mln/morpho/vector_median.hh (revision 3153) +++ abraham/mln/morpho/vector_median.hh (working copy) @@ -1,4 +1,4 @@ -// Copyright (C) 2007, 2008 EPITA Research and Development Laboratory +// Copyright (C) 2007, 2008, 2009 EPITA Research and Development Laboratory // (LRDE) // // This file is part of the Olena Library. This library is free Index: abraham/mln/morpho/topo_wst.hh --- abraham/mln/morpho/topo_wst.hh (revision 3153) +++ abraham/mln/morpho/topo_wst.hh (working copy) @@ -1,4 +1,4 @@ -// Copyright (C) 2008 EPITA Research and Development Laboratory (LRDE) +// Copyright (C) 2008, 2009 EPITA Research and Development Laboratory (LRDE) // // This file is part of the Olena Library. This library is free // software; you can redistribute it and/or modify it under the terms @@ -52,6 +52,42 @@ 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; + } + + + template <class I, class N> struct topo_wst { @@ -59,6 +95,8 @@ // Typedefs // -------- typedef mln_site(I) site; + typedef I image_t; + typedef N neighborhood_t; // Component tree management // ------------------------- @@ -88,14 +126,10 @@ // children.append(n); children.push_back(p); } - - mln_value(I) min() { return mln_min(mln_value(I)); } - mln_value(I) max() { return mln_max(mln_value(I)); } - void revert () { level = max() - level + min(); } - }; // struct node - const mln_exact(N)& nbh; + public: + mln_ch_value(I, site) Par_node; mln_ch_value(I, site) Par_tree; @@ -105,10 +139,11 @@ mln_ch_value(I, node) nodes; site Root; + private: + void MakeSet_tree(site x); void MakeSet_node(site x); site Find_tree(site x); - void swap(site& x, site& y); site Find_node(site x); void BuildComponentTree(); site MergeNode(site& node1, site& node2); @@ -116,10 +151,6 @@ site Link_node(site x, site y); node MakeNode(int level); - - // Watershed algorithm - // ------------------- - private: mln_ch_value(I, bool) isproc; @@ -128,10 +159,9 @@ topo_wst(const Image<I>& i, const Neighborhood<N>& n); - - // Result public: - mln_exact(I) pima; + const I &ima; + const N &nbh; public: void go(); @@ -141,12 +171,12 @@ site max (p_set<site>& components); bool highest_fork (p_set<site>& components); bool highest_fork (p_set<site>& components, site &r); - bool w_constructible(site p, site &r); - bool w_constructible(site p); - void topo_watershed(); - // Optimized LCA Algorithm + + public: + site lca (site a, site b); + private: int *euler; int *depth; @@ -161,7 +191,6 @@ void build_euler_tour_rec(site p, int &position, int d); void build_euler_tour(); void build_minim (); - site lca (site a, site b); void removeOneSonNodes(site *p, mln_ch_value(I, site) &newPar_node); void compressTree(); }; // struct topo_wst @@ -172,15 +201,15 @@ template <class I, class N> topo_wst<I, N>::topo_wst(const Image<I>& i, const Neighborhood<N>& n) - : nbh(exact(n)), - Par_node(exact(i).domain(), exact(i).border()), + : Par_node(exact(i).domain(), exact(i).border()), Par_tree(exact(i).domain(), exact(i).border()), Rnk_tree(exact(i).domain(), exact(i).border()), Rnk_node(exact(i).domain(), exact(i).border()), subtreeRoot(exact(i).domain(), exact(i).border()), nodes(exact(i).domain(), exact(i).border()), isproc(exact(i).domain(), exact(i).border()), - pima(exact(i)) + ima(exact(i)), + nbh(exact(n)) { } @@ -207,14 +236,6 @@ } template <class I, class N> - void topo_wst<I, N>::swap(site& x, site& y) - { - site memo = x; - x = y; - y = memo; - } - - template <class I, class N> typename topo_wst<I, N>::site topo_wst<I, N>::Find_node(site x) { if (Par_node(x) != x) @@ -227,16 +248,9 @@ { BuildComponentTree(); compressTree(); - arith::revert_inplace(pima); - // arith::revert_inplace(nodes); - mln_piter(image2d<node>) p (nodes.domain()); - for_all(p) - nodes(p).revert(); build_euler_tour(); build_minim(); - topo_watershed(); - arith::revert_inplace(pima); } template <class I, class N> @@ -246,7 +260,7 @@ // Sort the sites in increasing order p_array<mln_site(I)> S; - S = level::sort_psites_increasing(pima); + S = level::sort_psites_increasing(ima); // Clear the marker map data::fill(isproc, false); @@ -258,7 +272,7 @@ // if (subtreeRoot.hold(p)) subtreeRoot(p) = p; // if (nodes.hold(p)) - nodes(p) = MakeNode(pima(p)); + nodes(p) = MakeNode(ima(p)); } @@ -273,7 +287,7 @@ mln_niter(N) q(nbh, ip); for_all(q) - if (pima.has(q) and isproc(q) and pima(q) <= pima(p)) + if (ima.has(q) and isproc(q) and ima(q) <= ima(p)) { site adjCanonicalElt = Find_tree(q); site adjNode = Find_node(subtreeRoot(adjCanonicalElt)); @@ -333,7 +347,7 @@ typename topo_wst<I, N>::site topo_wst<I, N>::Link_tree(site x, site y) { if (Rnk_tree(x) > Rnk_tree(y)) - swap(x, y); + std::swap(x, y); else if (Rnk_tree(x) == Rnk_tree(y)) Rnk_tree(y) += 1; @@ -345,7 +359,7 @@ typename topo_wst<I, N>::site topo_wst<I, N>::Link_node(site x, site y) { if (Rnk_node(x) > Rnk_node(y)) - swap(x, y); + std::swap(x, y); else if (Rnk_node(x) == Rnk_node(y)) Rnk_node(y) += 1; @@ -363,38 +377,6 @@ return n; } - template <class I, class N> - typename topo_wst<I, N>::site topo_wst<I, N>::min (p_set<site>& components) - { - if (components.nsites() == 0) - return site(-1, -1); - - typename p_set<site>::fwd_piter it(components); - site min = components[0]; - - for_all(it) - if (pima(it.to_site()) < pima(min)) - min = it.to_site(); - - return min; - } - - template <class I, class N> - typename topo_wst<I, N>::site topo_wst<I, N>::max (p_set<site>& components) - { - if (components.nsites() == 0) - return site(-1, -1); - - typename p_set<site>::fwd_piter it(components); - site max = components[0]; - - for_all(it) - if (pima(it.to_site()) > pima(max)) - max = it.to_site(); - - return max; - } - template <class I, class N> bool topo_wst<I, N>::highest_fork (p_set<site>& components, site &r) @@ -435,140 +417,6 @@ } template <class I, class N> - bool topo_wst<I, N>::w_constructible(site p) { - site i; - return w_constructible(p, i); - } - - template <class I, class N> - bool topo_wst<I, N>::w_constructible(site p, site &r) - { - mln_niter(N) q(nbh, p); - p_set<site> v; - - for_all(q) - if (pima.has(q) && pima(q) > pima(p)) - v.insert(Par_node(q)); - - if (v.nsites() == 0) - return false; - - if (v.nsites() == 1) { - r = v[0]; - return true; - } - - site - c = max(v), - cmax = c; - - typename p_set<site>::fwd_piter it(v); - for_all(it) - { - // Can't remove the site from the set - if (it.to_site() == cmax) - continue; - - site c1 = lca(c, it.to_site()); - - if (c1 != it.to_site()) - c = c1; - } - - if (nodes(c).level <= pima(p)) - return false; - - r = c; - return true; - } - - - template <class I, class N> - void topo_wst<I, N>::topo_watershed() - { - // Maxima components - mln_ch_value(I, bool) cmax(pima.domain(), pima.border()); - - // Mark enqueued sites - mln_ch_value(I, bool) enqueued(pima.domain(), pima.border()); - - p_priority< mln_value(I), p_queue_fast<site> > l; - // p_queue < site > m; - std::queue<site> m; - mln_value(I) max = mln_max(mln_value(I)); - - std::cout << "Init" << std::endl; - - // Flag C-maxima - data::fill(cmax, false); - data::fill(enqueued, false); - - mln_piter(I) it(Par_node.domain()); - for_all(it) - // if (nodes(Par_node(it.to_site())).children.nsites() == 0) - if (nodes(Par_node(it.to_site())).children.size() == 0) - { - cmax(it.to_site()) = true; - m.push(it.to_site()); - } - - while (!m.empty()) - { - mln_niter(N) q(nbh, m.front()); - for_all(q) - if (cmax.has(q) && !cmax(q) && !enqueued(q)) - { - enqueued(q) = true; - l.push(pima(q), q); - } - m.pop(); - } - - - std::cout << "end" << std::endl; - - // Main loop - while (!l.is_empty()) - { - site x = l.front(); - l.pop(); - enqueued(x) = false; - - site c; - bool is_w = w_constructible(x, c); - - if (is_w) - { - pima(x) = nodes(c).level; - Par_node(x) = c; - - // if (nodes(c).children.nsites() == 0) - if (nodes(c).children.size() == 0) - cmax(x) = true; - else - // if (nodes(c).children.nsites() > 1) - if (nodes(c).children.size() == 1) - std::cerr << "ERREUR COMPOSANTE BRANCHE " << nodes(c).children.size() << std::endl; - - - mln_niter(N) q(nbh, x); - for_all(q) - if (pima.has(q) && !cmax(q) && !enqueued(q)) - { - enqueued(q) = true; - - l.push(pima(q), q); - } - } // if (c != site(-1, -1)) - } // while(!l.empty()) - - for_all(it) - pima(it.to_site()) = nodes(Par_node(it.to_site())).level; - } - - - - template <class I, class N> void topo_wst<I, N>::compute_ctree_size (site p) { ctree_size += 1; @@ -739,6 +587,157 @@ 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); + + mln_niter(N) q(nbh, p); + p_set<mln_site(I)> v; + + for_all(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), + 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(ima.domain(), ima.border()); + + // Mark enqueued sites + mln_ch_value(I, bool) enqueued(ima.domain(), ima.border()); + + 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()) + { + mln_niter(N) q(nbh, m.front()); + 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; + + + mln_niter(N) q(nbh, x); + 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