[Olena-patches] 1871: Add functional M-Watershed algorithm and tests.

https://svn.lrde.epita.fr/svn/oln/trunk/milena Index: ChangeLog from Alexandre Abraham &lt;abraham(a)lrde.epita.fr&gt; Add functional M-Watershed algorithm and tests. * sandbox/abraham/morpho/basic_najman.hh (mln::morpho::basic_najman::lca) (mln::morpho::basic_najman::highest_fork) (mln::morpho::basic_najman::m-destructible) (mln::morpho::basic_najman::w-destructible) (mln::morpho::basic_najman::m-watershed) (mln::morpho::basic_najman::w-watershed): New methods. * sandbox/abraham/morpho/test.cc: Add debug tests. * sandbox/abraham/morpho/test_component_tree.cc: New, Add tests for component tree algorithms. * sandbox/abraham/morpho/images/result_m_watershed.pgm: New, Reference for m-watershed algorithm. * sandbox/abraham/morpho/images/result_topo_watershed.pgm: New, Reference for w-watershed algorithm. * sandbox/abraham/morpho/images/test_watershed.pgm: New, Test image for watershed algorithms. * sandbox/abraham/morpho/images/test_component_mapping.pgm: New, Reference image for component mapping. * sandbox/abraham/morpho/test_watershed.cc: New, Add tests for watershed algorithms. * sandbox/abraham/morpho/Makefile: Add rules to compile new tests. Makefile | 16 +- basic_najman.hh | 295 +++++++++++++++++++++++++++++++------- images/test_component_mapping.pgm | 5 test.cc | 15 - test_component_tree.cc | 224 ++++++++++++++++++++++++++++ test_watershed.cc | 55 +++++++ 6 files changed, 542 insertions(+), 68 deletions(-) Index: sandbox/abraham/morpho/basic_najman.hh --- sandbox/abraham/morpho/basic_najman.hh (revision 1870) +++ sandbox/abraham/morpho/basic_najman.hh (working copy) @@ -1,6 +1,10 @@ #include <mln/level/sort_psites.hh> #include <mln/level/fill.hh> #include <mln/core/image2d.hh> +#include <mln/core/p_set.hh> +//#include <mln/util/greater_psite.hh> +#include <queue> +#include <set> namespace mln { @@ -28,19 +32,19 @@ it.next()) children.append(it.to_psite()); } + void addChild(const psite n) { children.append(n); } }; // struct node - + I pima; const Image<I>& ima; const Neighborhood<N>& nbh; mln_ch_value(I, psite) Par_node; - mln_ch_value(I, psite) Par_tree; - // Par_node is handled by par (from ap_maxtree) + mln_ch_value(I, int) Rnk_tree; mln_ch_value(I, int) Rnk_node; mln_ch_value(I, psite) subtreeRoot; @@ -52,7 +56,8 @@ basic_najman(const Image<I>& i, const Neighborhood<N>& n) - : ima(i), + : pima(exact(i)), + ima(i), nbh(n), Par_node(exact(i).domain(), exact(i).border()), Par_tree(exact(i).domain(), exact(i).border()), @@ -73,7 +78,7 @@ void MakeSet_node(psite x) { - Par_node(x) = x; // was Par_node(x) = x; + Par_node(x) = x; Rnk_node(x) = 0; } @@ -91,24 +96,6 @@ y = memo; } -// bool is_root(const psite& x) const -// { -// return Par_node(x) == x; -// } - -// bool is_level_root(const psite& x) const -// { -// return is_root(x) || exact(ima)(Par_node(x)) < exact(ima)(x); -// } - -// psite find_level_root(const psite& x) -// { -// if (is_level_root(x)) -// return x; -// else -// return Par_node(x) = find_level_root(Par_node(x)); -// } - psite Find_node(psite x) { if (Par_node(x) != x) @@ -148,7 +135,6 @@ for (int ip = 0; ip < int(S.npoints()); ++ip) { psite p = S[ip]; - // isproc(p) = true; psite curTree = Find_tree(p); psite curNode = Find_node(subtreeRoot(curTree)); @@ -162,36 +148,22 @@ psite adjNode = Find_node(subtreeRoot(adjTree)); if (curNode != adjNode) { - std::cout << "curNode != adjNode; " << nodes(curNode).level << " vs. " << nodes(adjNode).level << std::endl; if (nodes(curNode).level == nodes(adjNode).level) { - // curNode = MergeNode(adjNode, curNode); - psite tmpNode = MergeNode(adjNode, curNode); - /*if (tmpNode == curNode) - nodes(curNode).addChildren(nodes(adjNode)); - else - nodes(adjNode).addChildren(nodes(curNode)); - nodes(adjNode) = nodes(curNode);*/ - curNode = tmpNode; + curNode = MergeNode(adjNode, curNode); } else { - // we have nodes[curNode].level < nodes[adjNode].level nodes(curNode).addChild(adjNode); - //apparemment NON :Par_node[adjNode] = curNode; // car adjNode devient fils de curNode nodes(curNode).area += nodes(adjNode).area; nodes(curNode).highest += nodes(adjNode).highest; } - // curTree = Link_tree(adjTree, curTree); - // subtreeRoot(curTree) = curNode; - } curTree = Link_tree(adjTree, curTree); subtreeRoot(curTree) = curNode; } isproc(p) = true; } - Root = subtreeRoot(Find_tree(Find_node(psite(0, 0)))); // Pour garder une map de correspondance point <-> local_root // for (int ip = 0; ip < int(S.size()); ++ip) // { @@ -203,6 +175,7 @@ for_all(r) Par_node(r) = Find_node(r); + Root = subtreeRoot(Find_tree(Find_node(psite(0, 0)))); } @@ -243,7 +216,7 @@ else if (Rnk_node(x) == Rnk_node(y)) Rnk_node(y) += 1; - Par_node(x) = y; // was Par_node(x) = y; + Par_node(x) = y; return y; } @@ -275,40 +248,254 @@ void print_tree(psite p) { node& n = nodes(p); - std::cout << "psite " << p << "=" << (p.row() * exact(subtreeRoot).domain().len(1) + p.col()) << " : "; + std::cout << "psite " << p << "(" << n.level << ")=" << (p.row() * exact(subtreeRoot).domain().len(1) + p.col()) << " : "; typename p_array<mln_psite(I)>::fwd_piter it(n.children); - for (it.start(); - it.is_valid(); - it.next()) + for_all(it) { psite q = it.to_psite(); std::cout << q << "=" << (q.row() * subtreeRoot.domain().len(1) + q.col()) << " "; } std::cout << std::endl; - for (it.start(); - it.is_valid(); - it.next()) + for_all(it) { print_tree(it.to_psite()); } } - void print_sup_tree() + psite lca (psite a, psite b) { - psite max = psite(0, 0); + return lca_rec(Root, Par_node(a), Par_node(b)); + } - for(unsigned int i = 0; i < nodes.domain().len(0); ++i) - for(unsigned int j = 0; j < nodes.domain().len(1); ++j) + psite lca_rec (psite cur, psite a, psite b) + { + // FIXME : naive implementation, make something better + + if (cur == a) + return a; + + if (cur == b) + return b; + + if (nodes(cur).children.npoints() == 0) + return psite (-1, -1); + + psite tmp, acc = psite(-1, -1); + int n = 0; + typename p_array<mln_psite(I)>::fwd_piter it(nodes(cur).children); + for_all(it) + { + tmp = lca_rec(it.to_psite(), a, b); + if (tmp != psite(-1, -1)) + { + if (tmp == a) { - if (nodes(psite(i, j)).area > nodes(max).area) - max = psite(i, j); + acc = tmp; + n++; + continue; } + if (tmp == b) + { + acc = tmp; + n++; + continue; + } + return tmp; + } + } + if (n == 2) + return cur; + + return acc; + } + + psite min (p_set<psite>& components) + { + if (components.npoints() == 0) + return psite(-1, -1); + + typename p_set<psite>::fwd_piter it(components); + psite min = components[0]; + + for_all(it) + if (nodes(it.to_point()).level > nodes(min).level) + min = it.to_point(); + + return min; + } + + psite highest_fork (p_set<psite>& components) + { + if (components.npoints() == 0) + return psite(-1, -1); + + psite hfork = components[0], min = hfork; + typename p_set<psite>::fwd_piter it(components); - print_tree(max); + for_all(it) + hfork = lca(hfork, it.to_point()); + + if (nodes(min).level == nodes(hfork).level) + return psite(-1, -1); + + return hfork; } + psite w_destructible(psite p) + { + mln_niter(N) q(nbh, p); + p_set<psite> v; + + for_all(q) + if (exact(ima)(q) < exact(ima)(p)) + v.append(Par_node(q)); + + if (v.npoints() == 0) + return psite(-1, -1); + + psite hf = highest_fork(v); + + if (hf == psite(-1, -1)) + return min(v); + + if (nodes(hf).level <= exact(ima)(p)) + return hf; + + return psite(-1, -1); + } + + psite m_destructible(psite p) + { + mln_niter(N) q(nbh, p); + p_set<psite> v; + + for_all(q) + if (exact(ima)(q) < exact(ima)(p)) + v.insert(Par_node(q)); + + if (v.npoints() == 0) + return psite(-1, -1); + + if (nodes(min(v)).children.npoints() != 0) + return (psite(-1, -1)); + + psite hf = highest_fork(v); + + if (hf == psite(-1, -1)) + return min(v); + + return psite(-1, -1); + } + + void m_watershed () + { + + // FIXME : make a better priority function + // typedef + // std::priority_queue< psite, std::vector<psite>, util::greater_psite<I> > + // ordered_queue_type; + + + std::priority_queue<psite> l; + + // Clear the marker map + level::fill(isproc, false); + mln_piter(I) it(exact(ima).domain()); + + for_all(it) + if (m_destructible(it.to_point()) != psite(-1, -1)) + { + l.push(it.to_point()); + isproc(it.to_point()) = true; + } + + psite p, i; + while (!l.empty()) + { + p = l.top(); + l.pop(); + isproc(p) = false; + + i = m_destructible(p); + + if (i != psite(-1, -1)) + { + pima(p) = nodes(i).level ; + Par_node(p) = i; + mln_niter(N) q(nbh, p); + + for_all(q) + if (!isproc(q)) + if (m_destructible(q) != psite(-1, -1)) + { + l.push(q); + isproc(q) = true; + } + } + } + } + + void w_watershed() + { + p_array< p_set<psite> > L; + // TODO : replace int with the type of value + mln_ch_value(I, int) K; + mln_ch_value(I, psite) H; + + mln_piter(I) it(exact(ima).domain()); + + psite i; + for_all(it) + { + psite p = it.to_point(); + + i = w_destructible(p); + if (i != psite(-1, -1)) + { + L[nodes(i).level].insert(p); + K(p) = nodes(i).level; + H(p) = i; + + typename p_array< p_set<psite> >::fwd_piter n(L); + for_all(n) + { + while (!n.empty()) + { + psite p = *n.begin(); + n.erase(p); + // TODO : replace int with the type of value + int k = nodes(p).level + 1; + if (K(p) == k) + { + pima(p) = k; + Par_node(p) = H(p); + + mln_niter(N) q(nbh, p); + for_all(q) + if (k < pima(q)) + { + psite j = w_destructible(q); + if (j != psite(-1, -1)) + K(q) = -1; + else + if (K(q) != nodes(j).level) + { + n.insert(q); + K(q) = i-1; + H(q) = j; + } + } + } + } + } + } + } + + } + + }; // struct basic_najman }; // namespace morpho Index: sandbox/abraham/morpho/test.cc --- sandbox/abraham/morpho/test.cc (revision 1870) +++ sandbox/abraham/morpho/test.cc (working copy) @@ -17,21 +17,13 @@ using namespace mln; using value::int_u8; - // component_tree< image2d<int_u8> > c; - using namespace mln; using value::int_u8; image2d<int_u8> input; io::pgm::load(input, "./images/test_component_tree.pgm"); - - for(unsigned int i = 0; i < input.domain().len(0); ++i) - for(unsigned int j = 0; j < input.domain().len(1); ++j) - input.at(i, j) = 255 - input.at(i, j); - - - morpho::basic_najman< image2d<int_u8>, neighb2d> n(input, c8()); + morpho::basic_najman< image2d<int_u8>, neighb2d> n(input, c4()); n.image_output(input); n.go(); n.image_output(n.Par_tree); @@ -42,7 +34,10 @@ std::cout << " --------------------------- " << std::endl; n.image_output(n.Rnk_node); std::cout << " --------------------------- " << std::endl; - n.print_sup_tree(); + n.print_tree(n.Root); + std::cout << " --------------------------- " << std::endl; + + n.m_watershed(); // image2d<int_u8> output = morpho::topo_wst(input, c4()); // io::pgm::save(output, "out.pgm"); Index: sandbox/abraham/morpho/test_component_tree.cc --- sandbox/abraham/morpho/test_component_tree.cc (revision 0) +++ sandbox/abraham/morpho/test_component_tree.cc (revision 0) @@ -0,0 +1,224 @@ +#include <mln/core/image2d.hh> +#include <mln/core/window2d.hh> +#include <mln/core/neighb2d.hh> +#include <mln/core/p_set.hh> + +#include <mln/value/int_u8.hh> + +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +#include "basic_najman.hh" +#include <string> +#include <iostream> + +int print_and_exit (std::string s) +{ + std::cerr << s << std::endl; + return 1; +} + +template <typename T> +void swap (T& x, T& y) +{ + T memo = x; + x = y; + y = memo; +} + +template <typename T> +class test_tree; + +template <typename T> +class test_tree +{ +public: + + test_tree(T pdata) : data(pdata), children(NULL), brother(NULL) { } + + T data; + test_tree<T> *children; + test_tree<T> *brother; + + test_tree<T>* insert_child(test_tree<T> *n) + { + test_tree<T> **p = &children; + while (*p && (**p).data < n->data) + p = &((**p).brother); + n->brother = *p; + *p = n; + return n; + } + + test_tree<T>* insert_child(T pdata) + { + test_tree<T> *n = new test_tree<T>(pdata); + return insert_child(n); + } +}; + +template <typename T> +bool compare_tree(test_tree<T>* t1, test_tree<T>* t2) +{ + if (!t1 && !t2) + return true; + + if ((!t1 || !t2) || (t1->data != t2->data)) + return false; + + return (compare_tree(t1->brother, t2->brother) + && compare_tree(t1->children, t2->children)); +} + +template <typename T> +void show_tree (test_tree<T> *t) +{ + if (!t) + return ; + + if (t->brother) + { + std::cout << t->data << " - "; + show_tree(t->brother); + } + else + std::cout << t->data << std::endl; + + if (t->children) + { + std::cout << "sons of " << t->data << " : "; + show_tree(t->children); + } +} + + +using namespace mln; +using value::int_u8; + + +test_tree<image2d<int_u8>::psite>* test_convert (morpho::basic_najman<image2d<int_u8>, neighb2d>& ima, image2d<int_u8>::psite root); + +test_tree<image2d<int_u8>::psite>* test_convert (morpho::basic_najman<image2d<int_u8>, neighb2d>& ima, image2d<int_u8>::psite root) +{ + morpho::basic_najman<image2d<int_u8>, neighb2d>::node& n = ima.nodes(ima.Par_node(root)); + + test_tree<image2d<int_u8>::psite> *res = new test_tree<image2d<int_u8>::psite>(ima.Par_node(root)); + + p_array<image2d<int_u8>::psite>::fwd_piter it (n.children); + for_all(it) + res->insert_child(test_convert(ima, it.to_psite())); + + return res; +} + +int main () +{ + typedef image2d<int_u8> image2dint; + + image2dint input, verif; + io::pgm::load(input, "./images/test_component_tree.pgm"); + io::pgm::load(verif, "./images/test_component_mapping.pgm"); + + morpho::basic_najman< image2d<int_u8>, neighb2d> n(input, c4()); + n.go(); + + /* + Component mapping : + + 3 3 3 3 9 2 2 2 2 + 1 1 1 9 7 9 2 2 2 + 1 1 9 7 7 7 9 2 2 + 1 9 6 9 7 9 5 9 2 + 1 1 9 1 7 2 7 2 2 + 1 1 1 1 7 2 2 2 2 + + + Component tree + + 9 + / \ + 7 6 + /|\ + 3 2 5 + | + 1 + */ + + // Component mapping test + + image2dint::fwd_piter it(input.domain()); + image2dint::psite c[10]; + c[1] = n.Par_node(image2dint::psite(1, 0)); + c[2] = n.Par_node(image2dint::psite(0, 5)); + c[3] = n.Par_node(image2dint::psite(0, 0)); + c[5] = n.Par_node(image2dint::psite(3, 6)); + c[6] = n.Par_node(image2dint::psite(3, 2)); + c[7] = n.Par_node(image2dint::psite(1, 4)); + c[9] = n.Par_node(image2dint::psite(0, 4)); + + bool id = true; + + for_all(it) + id = (c[verif(it.to_point())] == n.Par_node(it.to_point())); + + if (!id) + { + std::cerr << "Component mapping is fucked up!" << std::endl; + return 1; + } + + // Component tree test + + test_tree<image2dint::psite> *gen, ref(c[9]), *l1, *l2; + + // Create the ref + ref.insert_child(c[6]); + l1 = ref.insert_child(c[7]); + l1->insert_child(c[2]); + l1->insert_child(c[5]); + l2 = l1->insert_child(c[3]); + l2->insert_child(c[1]); + + // Create the generated tree + gen = test_convert(n, n.Par_node(n.Root)); + + // Compare the trees + if (!compare_tree(gen, &ref)) + { + std::cout << "Your tree :" << std::endl; + show_tree(gen); + + std::cout << "Reference tree :" << std::endl; + show_tree(&ref); + + std::cerr << "Component tree is fucked up!" << std::endl; + return 1; + } + + // Test the LCA alorithm + image2dint::psite p(0, 0), q(0, 6), r(2, 4); + if (n.lca(p, q) != image2dint::psite(2,4)) + return 1; + + if (n.lca(r, q) != image2dint::psite(2,4)) + return 1; + + // Test the Highest Fork + p_set<image2dint::psite> comp; + comp.insert(image2dint::psite(0,0)); + comp.insert(image2dint::psite(2,4)); + comp.insert(image2dint::psite(0,6)); + comp.insert(image2dint::psite(1,0)); + + if (n.highest_fork(comp) != image2dint::psite(2, 4)) + return 1; + + comp.insert(image2dint::psite(3,2)); + + if (n.highest_fork(comp) != image2dint::psite(1, 3)) + return 1; + + std::cout << "Test is successfull" << std::endl; + + return 0; +} Index: sandbox/abraham/morpho/images/result_m_watershed.pgm Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: sandbox/abraham/morpho/images/result_m_watershed.pgm ___________________________________________________________________ Name: svn:mime-type + application/octet-stream Index: sandbox/abraham/morpho/images/result_topo_watershed.pgm Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: sandbox/abraham/morpho/images/result_topo_watershed.pgm ___________________________________________________________________ Name: svn:mime-type + application/octet-stream Index: sandbox/abraham/morpho/images/test_watershed.pgm Cannot display: file marked as a binary type. svn:mime-type = application/octet-stream Property changes on: sandbox/abraham/morpho/images/test_watershed.pgm ___________________________________________________________________ Name: svn:mime-type + application/octet-stream Index: sandbox/abraham/morpho/images/test_component_mapping.pgm --- sandbox/abraham/morpho/images/test_component_mapping.pgm (revision 0) +++ sandbox/abraham/morpho/images/test_component_mapping.pgm (revision 0) @@ -0,0 +1,5 @@ +P5 +# CREATOR: GIMP PNM Filter Version 1.1 +9 6 +255 + 


 

 
   






 \ No newline at end of file Index: sandbox/abraham/morpho/test_watershed.cc --- sandbox/abraham/morpho/test_watershed.cc (revision 0) +++ sandbox/abraham/morpho/test_watershed.cc (revision 0) @@ -0,0 +1,55 @@ +#include <mln/core/image2d.hh> +#include <mln/core/window2d.hh> +#include <mln/core/neighb2d.hh> +#include <mln/core/p_set.hh> + +#include <mln/value/int_u8.hh> +#include <mln/level/compare.hh> + +#include <mln/io/pgm/load.hh> +#include <mln/io/pgm/save.hh> + +#include "basic_najman.hh" +#include <string> +#include <iostream> + +int print_and_exit (std::string s) +{ + std::cerr << s << std::endl; + return 1; +} + +int main () +{ + using namespace mln; + using value::int_u8; + + typedef image2d<int_u8> image2dint; + + image2dint input, mverif; + io::pgm::load(input, "./images/test_watershed.pgm"); + // io::pgm::load(input, "./images/lena.pgm"); + io::pgm::load(mverif, "./images/result_m_watershed.pgm"); + + morpho::basic_najman< image2d<int_u8>, neighb2d> n(input, c4()); + n.go(); + n.m_watershed(); + + bool equal = true; + + equal = n.pima == mverif; + + if (!equal) + { + std::cout << "M-Watershed broken" << std::endl; + return 1; + } + else + std::cout << "M-Watershed good as chocolate ice cream" << std::endl; + + // n.w_watershed(); + + io::pgm::save(n.pima, "out.pgm"); + + return 0; +} Index: sandbox/abraham/morpho/Makefile --- sandbox/abraham/morpho/Makefile (revision 1870) +++ sandbox/abraham/morpho/Makefile (working copy) @@ -1,8 +1,16 @@ -SOURCE=test.cc -OBJECTS=test.o CXXFLAGS=-Wall -W -I ../../.. -g -all:${OBJECTS} - g++ ${FLAGS} ${OBJECTS} -o test +test: test.o + g++ ${CXXFLAGS} test.o -o test +test_component_tree: test_component_tree.o + g++ ${CXXFLAGS} test_component_tree.o -o test_component_tree + +test_watershed: test_watershed.o + g++ ${CXXFLAGS} test_watershed.o -o test_watershed + +all: tree + +test_watershed.o: basic_najman.hh +test_component_tree.o: basic_najman.hh test.o: basic_najman.hh \ No newline at end of file