[Olena-patches] 1967: Add a new Watershed algorithm.

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Alexandre Abraham &lt;abraham(a)lrde.epita.fr&gt; Add a new Watershed algorithm. * abraham/tests/morpho/test_component_tree.cc: New. * abraham/tests/morpho/images: New (test images). * abraham/tests/morpho/images/+test_watershed.pgm: New. * abraham/tests/morpho/images/result_m_watershed.pgm: New. * abraham/tests/morpho/images/test_component_tree.pgm: New. * abraham/tests/morpho/images/test_2.pgm: New. * abraham/tests/morpho/images/test_3.pgm: New. * abraham/tests/morpho/images/test_4.pgm: New. * abraham/tests/morpho/images/result_topo_watershed.pgm: New. * abraham/tests/morpho/images/test_watershed.pgm: New. * abraham/tests/morpho/images/+irm6.pgm: New. * abraham/tests/morpho/images/test_component_mapping.pgm: New. * abraham/tests/morpho/test_watershed.cc: New. * abraham/tests/morpho/test_watershed_topo.cc: New. * abraham/tests/morpho/ref/src/lib/lwshedtopo.c: uncomment specific pieces of code to show debug. * abraham/tests/morpho/Makefile: Fix. * abraham/tests/io/tikz/Makefile: Fix. * abraham/mln/morpho/basic_najman.hh: Add the new watershed from Pink library. * abraham/mln/morpho/test_component_tree.cc: Remove. * abraham/mln/morpho/test_watershed.cc: Remove. * abraham/mln/morpho/test_watershed_topo.cc: Remove. * abraham/mln/morpho/Makefile: Remove. * abraham/mln/io/tikz/save.hh: Add the ability to tag values on images. * abraham/img/dots.pgm: New. mln/io/tikz/save.hh | 13 - mln/morpho/Makefile | 22 - mln/morpho/basic_najman.hh | 396 +++++++++++++++++++++++++++++----- mln/morpho/test_component_tree.cc | 230 ------------------- mln/morpho/test_watershed.cc | 103 -------- mln/morpho/test_watershed_topo.cc | 60 ----- tests/io/tikz/Makefile | 1 tests/morpho/Makefile | 1 tests/morpho/ref/src/lib/lwshedtopo.c | 140 ++++++------ tests/morpho/test_component_tree.cc | 2 tests/morpho/test_watershed.cc | 2 tests/morpho/test_watershed_topo.cc | 9 12 files changed, 424 insertions(+), 555 deletions(-) Index: abraham/tests/morpho/test_component_tree.cc --- abraham/tests/morpho/test_component_tree.cc (revision 1935) +++ abraham/tests/morpho/test_component_tree.cc (working copy) @@ -8,7 +8,7 @@ #include <mln/io/pgm/load.hh> #include <mln/io/pgm/save.hh> -#include "basic_najman.hh" +#include <mln/morpho/basic_najman.hh> #include <string> #include <iostream> Index: abraham/tests/morpho/images/+test_watershed.pgm Cannot display: file marked as a binary type. svn:mime-type = abraham/tests/morpho/images/+irm6.pgm Property changes on: abraham/tests/morpho/images/+test_watershed.pgm ___________________________________________________________________ Name: svn:mime-type + abraham/tests/morpho/images/+irm6.pgm Index: abraham/tests/morpho/images/result_m_watershed.pgm Cannot display: file marked as a binary type. svn:mime-type = abraham/tests/morpho/images/+irm6.pgm Property changes on: abraham/tests/morpho/images/result_m_watershed.pgm ___________________________________________________________________ Name: svn:mime-type + abraham/tests/morpho/images/+irm6.pgm Index: abraham/tests/morpho/images/test_component_tree.pgm Cannot display: file marked as a binary type. svn:mime-type = abraham/tests/morpho/images/+irm6.pgm Property changes on: abraham/tests/morpho/images/test_component_tree.pgm ___________________________________________________________________ Name: svn:mime-type + abraham/tests/morpho/images/+irm6.pgm Index: abraham/tests/morpho/images/test_2.pgm Cannot display: file marked as a binary type. svn:mime-type = abraham/tests/morpho/images/+irm6.pgm Property changes on: abraham/tests/morpho/images/test_2.pgm ___________________________________________________________________ Name: svn:mime-type + abraham/tests/morpho/images/+irm6.pgm Index: abraham/tests/morpho/images/test_3.pgm Cannot display: file marked as a binary type. svn:mime-type = abraham/tests/morpho/images/+irm6.pgm Property changes on: abraham/tests/morpho/images/test_3.pgm ___________________________________________________________________ Name: svn:mime-type + abraham/tests/morpho/images/+irm6.pgm Index: abraham/tests/morpho/images/test_4.pgm Cannot display: file marked as a binary type. svn:mime-type = abraham/tests/morpho/images/+irm6.pgm Property changes on: abraham/tests/morpho/images/test_4.pgm ___________________________________________________________________ Name: svn:mime-type + abraham/tests/morpho/images/+irm6.pgm Index: abraham/tests/morpho/images/result_topo_watershed.pgm Cannot display: file marked as a binary type. svn:mime-type = abraham/tests/morpho/images/+irm6.pgm Property changes on: abraham/tests/morpho/images/result_topo_watershed.pgm ___________________________________________________________________ Name: svn:mime-type + abraham/tests/morpho/images/+irm6.pgm Index: abraham/tests/morpho/images/test_watershed.pgm Cannot display: file marked as a binary type. svn:mime-type = abraham/tests/morpho/images/+irm6.pgm Property changes on: abraham/tests/morpho/images/test_watershed.pgm ___________________________________________________________________ Name: svn:mime-type + abraham/tests/morpho/images/+irm6.pgm Index: abraham/tests/morpho/images/+irm6.pgm Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: abraham/tests/morpho/images/+irm6.pgm ___________________________________________________________________ Name: svn:mime-type + application/octet-stream Index: abraham/tests/morpho/images/test_component_mapping.pgm Cannot display: file marked as a binary type. svn:mime-type = abraham/tests/morpho/images/+irm6.pgm Property changes on: abraham/tests/morpho/images/test_component_mapping.pgm ___________________________________________________________________ Name: svn:mime-type + abraham/tests/morpho/images/+irm6.pgm Index: abraham/tests/morpho/test_watershed.cc --- abraham/tests/morpho/test_watershed.cc (revision 1935) +++ abraham/tests/morpho/test_watershed.cc (working copy) @@ -7,7 +7,7 @@ #include <mln/io/pgm/load.hh> #include <mln/io/pgm/save.hh> -#include "basic_najman.hh" +#include <mln/morpho/basic_najman.hh> #include <string> #include <iostream> Index: abraham/tests/morpho/test_watershed_topo.cc --- abraham/tests/morpho/test_watershed_topo.cc (revision 1935) +++ abraham/tests/morpho/test_watershed_topo.cc (working copy) @@ -9,7 +9,7 @@ #include <mln/io/pgm/load.hh> #include <mln/io/pgm/save.hh> -#include "basic_najman.hh" +#include <mln/morpho/basic_najman.hh> #include <string> #include <iostream> @@ -32,7 +32,8 @@ // io::pgm::load(input, "./images/test_watershed.pgm"); // io::pgm::load(input, "./images/little_test.pgm"); - io::pgm::load(input, "./images/test_2.pgm"); + io::pgm::load(input, "./images/test_4.pgm"); + // io::pgm::load(input, "../../img/dots.pgm"); // io::pgm::load(input, "./images/+irm6.pgm"); // io::pgm::load(input, "./images/lena_light.pgm"); Index: abraham/tests/morpho/ref/src/lib/lwshedtopo.c --- abraham/tests/morpho/ref/src/lib/lwshedtopo.c (revision 1966) +++ abraham/tests/morpho/ref/src/lib/lwshedtopo.c (working copy) @@ -266,13 +266,13 @@ } } } -#ifdef _DEBUG_ + //#ifdef _DEBUG_ for (i=0; i<logn; i++) { for (j=0; j<nbRepresent; j++) printf("M[%d][%d] = %d - ", i, j, Minim[i][j]); printf("\n"); } -#endif + //#endif return Minim; } @@ -491,10 +491,10 @@ if (CT->tabnodes[c].nbsons > 1) // noeud { } -#ifdef PARANO + //#ifdef PARANO else printf("%s : ERREUR COMPOSANTE BRANCHE!!!\n", F_NAME); -#endif + //#endif // empile les c-voisins de x non marques MASSIF ni EN_FAH switch (connex) { Index: abraham/tests/morpho/Makefile --- abraham/tests/morpho/Makefile (revision 1966) +++ abraham/tests/morpho/Makefile (working copy) @@ -16,4 +16,6 @@ .PHONY: clean distclean ## Depedencies -tikz.o: ../../../mln/io/tikz/save.hh \ No newline at end of file +test_watershed_topo.cc: ../../mln/morpho/basic_najman.hh +test_watershed.cc: ../../mln/morpho/basic_najman.hh +test_componenet_tree.cc: ../../mln/morpho/basic_najman.hh \ No newline at end of file Index: abraham/tests/io/tikz/Makefile --- abraham/tests/io/tikz/Makefile (revision 1966) +++ abraham/tests/io/tikz/Makefile (working copy) @@ -21,7 +21,6 @@ rm -rf *~ rm -rf *.log rm -rf *.aux - rm -rf distclean: clean rm -rf $(EXEC) Index: abraham/mln/morpho/basic_najman.hh --- abraham/mln/morpho/basic_najman.hh (revision 1966) +++ abraham/mln/morpho/basic_najman.hh (working copy) @@ -7,6 +7,7 @@ #include <mln/math/sqr.hh> #include <queue> +#include <vector> #include <set> namespace mln @@ -158,14 +159,15 @@ - template <typename I, typename N> + template <typename I, typename N, typename B> class my_lesser_psite { public: typedef mln_psite(I) psite; my_lesser_psite(const Image<I>& ima, - const Neighborhood<N>& nbh); + const Neighborhood<N>& nbh, + const Image<B>& enqueued); /// Is \a x lesser than \a y? bool operator()(const psite& x, const psite& y); @@ -173,30 +175,33 @@ private: const I& ima_; const N& nbh_; + const B& enqueued_; }; - template <typename I, typename N> - my_lesser_psite<I, N> + template <typename I, typename N, typename B> + my_lesser_psite<I, N, B> make_my_lesser_psite(const Image<I>& ima, - const Neighborhood<N>& nbh); + const Neighborhood<N>& nbh, + const Image<B>& enqueued); # ifndef MLN_INCLUDE_ONLY - template <typename I, typename N> - my_lesser_psite<I, N>::my_lesser_psite(const Image<I>& ima, - const Neighborhood<N>& nbh) - : ima_ (exact(ima)), nbh_(exact(nbh)) + template <typename I, typename N, typename B> + my_lesser_psite<I, N, B>::my_lesser_psite(const Image<I>& ima, + const Neighborhood<N>& nbh, + const Image<B>& enqueued) + : ima_ (exact(ima)), nbh_(exact(nbh)), enqueued_(exact(enqueued)) { } - template <typename I, typename N> + template <typename I, typename N, typename B> bool - my_lesser_psite<I, N>::operator()(const mln_psite(I)& x, const mln_psite(I)& y) + my_lesser_psite<I, N, B>::operator()(const mln_psite(I)& x, const mln_psite(I)& y) { -#ifdef V1 +#if 0 if (ima_(x) == ima_(y)) { @@ -212,27 +217,47 @@ if (ima_.has(b) and ima_(b) > ima_(y)) ++my; - return mx > my; + return mx < my; } -#else +#endif +#if 0 if (ima_(x) == ima_(y)) { mln_value(I) mx = ima_(x), my = ima_(y); mln_niter(N) a(nbh_, x); for_all(a) - if (ima_.has(a) and ima_(a) > ima_(mx)) - mx = a; + if (ima_.has(a) and ima_(a) > mx) + mx = ima_(a); + + mln_niter(N) b(nbh_, y); + for_all(b) + if (ima_.has(b) and ima_(b) > my) + my = ima_(b); + + return mx < my; + } +#endif + +#if 1 + if (ima_(x) == ima_(y)) + { + int mx = 0, my = 0; + + mln_niter(N) a(nbh_, x); + for_all(a) + if (ima_.has(a) and enqueued_(a)) + ++mx; mln_niter(N) b(nbh_, y); for_all(b) - if (ima_.has(b) and ima_(b) > ima_(my)) - my = b; + if (ima_.has(b) and enqueued_(b)) + ++my; - return ima_(x) < ima(y); + return mx > my; } @@ -243,16 +268,154 @@ } - template <typename I, typename N> - my_lesser_psite<I, N> + template <typename I, typename N, typename B> + my_lesser_psite<I, N, B> make_my_lesser_psite(const Image<I>& ima, - const Neighborhood<N>& nbh) + const Neighborhood<N>& nbh, + const Image<B>& enqueued) { - return my_lesser_psite<I, N>(ima); + return my_lesser_psite<I, N, B>(ima, nbh, enqueued); } # endif // ! MLN_INCLUDE_ONLY + + + + + + + // We need a custom queue in fact... + + // #define NPRIO 65536 + + template <typename P, typename V> + class fah_queue { + + public : + + // Push + void push(P point, + V level) + { + if (!util_ || (level < level_)) + level_ = level; + + ++util_; + if (util_ > max_util_) + max_util_ = util_; + + if (is_empty(level)) + queues[level] = new pqueue(); + + queues[level]->push(point); + + } // void push(P point, V level) + + + // Top + P top () + { + if (!util_) { + std::cerr << "Fah vide!" << std::endl; + exit (1); // FIXME! + } + return top(level_); + } // P top () + + + // Pop + void pop () + { + if (!util_) { + std::cerr << "Fah vide!" << std::endl; + return; + } + pop(level_); + } // void pop () + + + // Is Empty + bool is_empty () { + return !util_; + } // bool is_empty () + + + // Ctor + fah_queue() : + level_(0), + util_(0), + max_util_(0) + { } // fah_queue() + + + // Dtor + ~fah_queue() + { + typename qmap::iterator i; + + for (i = queues.begin(); i != queues.end(); ++i) + delete (i->second); + } + + + // Print + + + // Flush + + + private : + + V level_; /* niveau a partir duquel des listes existent */ + int util_; /* nombre de points courant dans la fah */ + int max_util_; /* nombre de points utilises max (au cours du temps) */ + + typedef std::queue <P> pqueue; + typedef std::map <V, pqueue *> qmap; + qmap queues; + + + // Is Empty + bool is_empty (V level) { + return ((queues.find(level) == queues.end()) || queues[level]->empty()); + } // bool is_empty (V level) + + + // Top + P top (V level) + { + if (is_empty(level)) + { + std::cerr << "erreur Fah vide au niveau " << level << std::endl; + exit (1); // FIXME! + } + return queues[level]->front(); + } + + + // Pop + void pop (V level) + { + if (is_empty(level)) + { + std::cerr << "erreur Fah vide au niveau " << level << std::endl; + return; + } + + util_--; + queues[level]->pop(); + + typename qmap::iterator i = queues.find(level); + for (; i != queues.end() && i->second->empty(); ++i) + ; + + level_ = i->first; + } + + }; // class fah_queue + + } // end of namespace mln::util @@ -713,8 +876,6 @@ mln_niter(N) q(nbh, p); p_set<psite> v; - psite pmax = p; - for_all(q) if (pima.has(q) && pima(q) > pima(p)) v.insert(Par_node(q)); @@ -726,16 +887,17 @@ return v[0]; psite - c = max(v); + c = max(v), + cmax = c; typename p_set<psite>::fwd_piter it(v); for_all(it) { // Can't remove the point from the set - if (it.to_point() == c) + if (it.to_point() == cmax) continue; - psite c1 = lca(c, it.to_point()); + psite c1 = lca_opt(c, it.to_point()); if (c1 != it.to_point()) c = c1; } @@ -919,7 +1081,9 @@ { init(); + std::cout << "Build component tree..." << std::endl; BuildComponentTree(); + std::cout << " done" << std::endl; mln_piter(I) p (pima.domain()); for_all(p) @@ -927,7 +1091,17 @@ revert_tree(Root); + std::cout << "Build euler tour..." << std::endl;; + build_euler_tour(); + std::cout << " done" << std::endl; + + std::cout << "Build minim table..." << std::endl;; + build_minim(); + std::cout << " done" << std::endl; + + std::cout << "Watersheding..." << std::endl;; topo_watershed(); + std::cout << " done" << std::endl; for_all(p) pima(p) = 255 - pima(p); @@ -945,15 +1119,8 @@ // Not used // level::fill(isproc, false); - - typedef - std::priority_queue< psite, std::vector<psite>, util::lesser_psite<I> > - ordered_queue_type; - - ordered_queue_type l(util::make_lesser_psite(pima)); - - - + util::fah_queue < psite, mln_value(I) > l; + mln_value(I) max = mln_max(mln_value(I)); // Flag C-maxima level::fill(cmax, false); @@ -962,16 +1129,6 @@ if (nodes(Par_node(it.to_point())).children.npoints() == 0) cmax(it.to_point()) = true; -#ifdef SLOW - - // Enqueue all - for_all(it) - { - enqueued(it.to_point()) = true; - l.push(it.to_point()); - } -#else // !SLOW - // Optimisation : enqueue minima's neighbours level::fill(enqueued, false); for_all(it) @@ -981,16 +1138,16 @@ if (cmax.has(q) && cmax(q)) { enqueued(it.to_point()) = true; - l.push(it.to_point()); + + l.push(it.to_point(), max - pima(it.to_point())); + break; } } -#endif // SLOW - // Main loop - while(!l.empty()) + while(!l.is_empty()) { psite x = l.top(); l.pop(); @@ -998,11 +1155,6 @@ psite c = w_constructible(x); - if (x == psite(3, 3)) - { - std::cout << "w-c : " << c << std::endl; - io::pgm::save(pima, "int.pgm"); - } if (c != psite(-1, -1)) { @@ -1017,14 +1169,15 @@ { } else - std::cerr << "ERREUR COMPOSANTE BRANCHE" << std::endl; + std::cerr << "ERREUR COMPOSANTE BRANCHE " << nodes(c).children.npoints() << std::endl; mln_niter(N) q(nbh, x); for_all(q) if (pima.has(q) && !cmax(q) && !enqueued(q)) { enqueued(q) = true; - l.push(q); + + l.push(q, max - pima(q)); } } // if (c != psite(-1, -1)) } // while(!l.empty()) @@ -1032,6 +1185,133 @@ for_all(it) pima(it.to_point()) = nodes(Par_node(it.to_point())).level; } + + + // Optimized LCA Algorithm + + int *euler; + int *depth; + int ctree_size; + std::map<psite, int> pos; + psite *repr; + + int *minim; + int **Minim; + + void compute_ctree_size (psite p) + { + ctree_size += 1; + node& n = nodes(p); + + typename p_array<mln_psite(I)>::fwd_piter it(n.children); + for_all(it) + compute_ctree_size(it.to_point()); + } + + void build_euler_tour_rec(psite 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_psite(I)>::fwd_piter it(n.children); + for_all(it) + { + build_euler_tour_rec(it.to_point(), position, d+1); + depth[position] = d; // Not optimized + euler[position] = pos[p]; + position++; + } + + } + + void 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 psite[size]; + + int position = 0; + build_euler_tour_rec(Root, position, 0); + } + + void 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 () + + psite lca_opt (psite a, psite 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)]]]; + } + } + }; // struct basic_najman }; // namespace morpho Index: abraham/mln/io/tikz/save.hh --- abraham/mln/io/tikz/save.hh (revision 1966) +++ abraham/mln/io/tikz/save.hh (working copy) @@ -165,9 +165,7 @@ { file << "\\path (" << p.col() << "," << max_row - p.row() << ") node[fill="; write_value(file, ima(p)); - file << "] (p_" << p.row() << "_" << p.col() << ") { \\color{"; - file << "black"; - file << "}" << "};" << std::endl; + file << "] (p_" << p.row() << "_" << p.col() << ") { }" << "};" << std::endl; } } @@ -184,6 +182,10 @@ min_col = geom::min_col(ima), max_col = geom::max_col(ima); + typedef mln_value(I) ival; + ival mid = (mln_max(ival) - mln_min(ival)) / 2; + mid += mln_min(ival); + point2d p; for (p.row() = min_row; p.row() <= max_row; ++p.row()) for (p.col() = min_col; p.col() <= max_col; ++p.col()) @@ -191,7 +193,10 @@ file << "\\path (" << p.col() << "," << max_row - p.row() << ") node[fill="; write_value(file, ima(p)); file << "] (p_" << p.row() << "_" << p.col() << ") { \\color{"; + if (ima(p) > mid) file << "black"; + else + file << "white"; file << "}" << val(p) << "};" << std::endl; } } Index: abraham/img/dots.pgm Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: abraham/img/dots.pgm ___________________________________________________________________ Name: svn:mime-type + application/octet-stream