[Olena-patches] 3257: Update Laurent's code.

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Thierry Geraud &lt;thierry.geraud(a)lrde.epita.fr&gt; Update Laurent's code. * theo/esiee/laurent/ismm09/main.cc: Augment. * theo/esiee/laurent/ismm09/lca.hh: New. * theo/color/segment.cc: Use 16bit labels. color/segment.cc | 12 - esiee/laurent/ismm09/lca.hh | 173 +++++++++++++++ esiee/laurent/ismm09/main.cc | 489 ++++++++++++++++++++++++++++++++++++------- 3 files changed, 595 insertions(+), 79 deletions(-) Index: theo/esiee/laurent/ismm09/main.cc --- theo/esiee/laurent/ismm09/main.cc (revision 3256) +++ theo/esiee/laurent/ismm09/main.cc (working copy) @@ -174,10 +174,15 @@ return e + (is_row_odd(e) ? down : right); } + + // Function-Object "e -> (l1, l2)". + + struct e_to_labels_t + { template <typename W, typename L> inline void - e_to_labels(const W& w, const point2d& e, L& l1, L& l2) + operator()(const W& w, const point2d& e, L& l1, L& l2) const { mln_precondition(w(e) == 0); l1 = 0; @@ -199,6 +204,7 @@ std::swap(l1, l2); mln_postcondition(l2 >= l1); } + }; } // end of namespace mln::cplx2d @@ -207,21 +213,402 @@ template <typename A, typename L> util::array<L> - sort_by_increasing_attributes(const util::array<A>& a, L n_basins) + sort_by_increasing_attributes(const util::array<A>& a, L l_max) { typedef std::pair<A,L> pair_t; std::vector<pair_t> v; - v.reserve(n_basins); - for (L l = 1; l <= n_basins; ++l) + v.reserve(l_max.next()); + + v.push_back(pair_t(mln_min(A), 0)); // First elt, even after sorting. + for (L l = 1; l <= l_max; ++l) v.push_back(pair_t(a[l], l)); std::sort(v.begin(), v.end()); - util::array<L> l(n_basins); - for (unsigned i = 0; i < n_basins; ++i) - l[i] = v[i].second; + util::array<L> ls(l_max.next()); + for (unsigned i = 1; i <= l_max; ++i) + ls[i] = v[i].second; + + return ls; + } + + + // Find-Root for a region labeled with 'l'. + + template <typename L> + inline + L find_root_l(util::array<L>& lpar, L l) + { + if (lpar[l] == l) return l; + else + return lpar[l] = find_root_l(lpar, lpar[l]); + } + + + + template <typename I> + inline + mln_psite(I) + find_root_e(I& z_epar, mln_psite(I) e) + { + if (z_epar(e) == e) + return e; + else + return z_epar(e) = find_root_e(z_epar, z_epar(e)); + } + + + + // Test emptiness of the queue of a set of regions. + + template <typename Qs, typename L> + bool test_q_emptiness(const Qs& qs, L l, util::array<L>& lpar) + { + L l_ = l; + while (lpar[l_] != l_) + { + if (! qs[l_].is_empty()) + return false; + l_ = lpar[l_]; + } + return true; + } + + + // Get smallest edge. + + template <typename Q_, + typename W, typename L, typename F, + typename E> + void + get_smallest_edge(Q_& q_, // in-out + L l, const W& w, util::array<L>& lpar, const F& e_to_labels, // in + E& e) // out + { + typedef mln_element(Q_) Q; // q_ is an array of queues with type Q. + + // Test that, when a region has merged, its edge queue has been + // emptied. + mln_invariant(test_q_emptiness(q_, l, lpar)); + + L lr = find_root_l(lpar, l); + Q& q = q_[lr]; + + while (! q.is_empty()) + { + e = q.pop_front(); + + L l1, l2; + e_to_labels(w, e, // input + l1, l2); // output + + mln_invariant(l1 != l2); + + L l1r = find_root_l(lpar, l1), + l2r = find_root_l(lpar, l2); + + mln_invariant(l1r == lr || l2r == lr); + + if (l1r == l2r) + // 'e' is an internal edge => forget it. + continue; + + // Otherwise 'e' has been found. + return; + } + + mln_invariant(0); // We should not be here! + } + + + + + // ################################################################### + // ################################################################### + // ################################################################### + // ################################################################### + // ################################################################### + + + + template <typename W, + typename G, + typename L, typename A, + typename F> + void + compute_pseudo_tree(const W& w, + const G& g, + const util::array<L>& ls, + const A& a, + const F& e_to_labels) + { + + typedef mln_value(G) T; // <--- Type of edge values. + typedef mln_psite(G) E; // <--- Type of edges. + + const L l_max = ls.nelements() - 1; + + mln_VAR( g_line, g | (pw::value(w) == 0) ); + + + // Edges -> Priority queue. + + typedef p_priority< T, p_queue<E> > Q; + util::array<Q> q(l_max.next()); + + { + L l1, l2; + + mln_piter(g_line_t) e(g_line.domain()); + for_all(e) + { + mln_invariant(w(e) == 0); + e_to_labels(w, e, // input + l1, l2); // output + q[l1].push(mln_max(T) - g(e), e); + q[l2].push(mln_max(T) - g(e), e); + } + } + + + + // Initialization. + // ----------------------------------- + + + + // Information "label l -> edge e". + + E null = E(0,0); // Impossible value. FIXME: lack of genericity. + + + util::array<E> edge(l_max.next()); + for (L l = 0; l <= l_max; ++l) + edge[l] = null; + + + util::array<L> lpar(l_max.next()); + for (L l = 0; l <= l_max; ++l) + lpar[l] = l; // Make-Set. + + + // To know if an edge is out of a given region (label l), we + // maintain the information about region merging using an + // union-find structure named "lpar". + // + // In the following "lpar[l]" is shortly denoted by lr, meaning + // l-root. + + + // Given a region R with label l and an edge e = (l1, l2) from the + // priority queue, we get know if that edge is out of the set of + // regions containing l: we shall have l1r = lr or l2r = lr. + + // Please note that an edge (l1, l2) is internal to a set of + // regions if l1r = l2r. + + + + + mln_ch_value(g_line_t, E) + epar, // Edge forest. + z_epar; // Auxiliary data: edge forest with compression and balancing. + + { + initialize(epar, g_line); + initialize(z_epar, g_line); + mln_piter(g_line_t) e(g_line.domain()); + for_all(e) + { + // Make-Set. + epar(e) = e; + z_epar(e) = e; + } + debug::println("epar (init):", epar); // epar(e) == e so we depict the edges! + } + + + + // GO GO GO !!! + + + for (unsigned i = 1; i < l_max; ++i) + { + L l = ls[i]; // Region label. + E e; + get_smallest_edge(q, // in-out + l, w, lpar, e_to_labels, // in + e); // out + // FIXME: the call below is performed above!!! + L l1, l2; + e_to_labels(w, e, l1, l2); + + L l1r = find_root_l(lpar, l1), + l2r = find_root_l(lpar, l2); + + // Consistency tests. + { + if (i > 1) + { + L former_l = ls[i-1]; + mln_invariant(a[l] >= a[former_l]); + } + mln_invariant(epar(e) == e); // 'e' has not been processed yet. + mln_invariant(l1 != l2); // It is a valid edge, i.e., between 2 different regions... + mln_invariant(l1r != l2r); // ...that are not already merged. + + L lr = find_root_l(lpar, l); + mln_invariant(lr == l // Either l is not deja-vu + || ((lr == l1r && lr != l2r) || // or l belongs to l1r xor l2r. + (lr == l2r && lr != l1r))); + } + + + + // aa(e) = a[l]; // FIXME: Re-activate. + + + + /* + std::cout << "l = " << l + << " e = " << e + << " (l1, l2) = (" << l1 << ", " << l2 << ") => " + << " merging R" << l1r << " and R" << l2r + << std::endl; + */ + + + // Merging both regions. + { + if (l2r < l1r) + std::swap(l1r, l2r); + mln_invariant(l2r > l1r); + lpar[l1r] = l2r; + + /* + std::cout << "q[l1r] = " << q[l1r] << std::endl; + std::cout << "q[l2r] = " << q[l2r] << std::endl; + */ + + + q[l2r].insert(q[l1r]); + q[l1r].clear(); + + + /* + std::cout << "q[l2r] = " << q[l2r] << std::endl + << std::endl; + */ + + + /* + // Displaying lpar. + { + std::cout << "lpar = "; + for (unsigned i = 1; i <= l_max; ++i) + { + L l = v[i].second; + std::cout << l << "->" << lpar[l] << " "; + } + std::cout << std::endl; + } + */ + } + + E e1 = edge[l1], + e2 = edge[l2]; + + if (e1 == null && e2 == null) + { + // New edge-component (singleton) + // Put differently: new edge-node! + edge[l1] = e; + edge[l2] = e; + // after: + // e + // / \ + // l1 l2 + } + else if (e1 != null && e2 != null) + { + // Two trees shall merge. + E e1r = find_root_e(z_epar, e1), + e2r = find_root_e(z_epar, e2); + mln_invariant(e1r != e2r); // Otherwise, there's a bug! + // before: + // e1r e2r + // |... |... + // e1 e2 + // / \ / \ + // l1 . l2 . + epar(e1r) = e; z_epar(e1r) = e; + epar(e2r) = e; z_epar(e2r) = e; + // after: + // e + // / \ + // e1r e2r + // |... |... + // e1 e2 + // / \ / \ + // l1 . l2 . + } + else if (e1 != null && e2 == null) + { + E e1r = find_root_e(z_epar, e1); + // before: + // e1r + // |... + // e1 null + // / \ / + // l1 . l2 + epar(e1r) = e; z_epar(e1r) = e; + edge[l2] = e; + // after: + // e + // / \ + // e1r l2 + // |... + // e1 + // / \ + // l1 . + } + else + { + mln_invariant(e1 == null && e2 != null); + E e2r = find_root_e(z_epar, e2); + epar(e2r) = e; z_epar(e2r) = e; + edge[l1] = e; + } + + } // end of "for every region with increasing attribute" + + + { + // Display edge tree. + + mln_ch_value(g_line_t, bool) deja_vu; + initialize(deja_vu, g_line); + data::fill(deja_vu, false); + + std::cout << "edge tree: " << std::endl; + for (L l = 1; l <= l_max; ++l) + { + std::cout << l << ": "; + E e = edge[l]; + while (! deja_vu(e)) + { + std::cout << e << " -> "; + deja_vu(e) = true; + e = epar(e); + } + std::cout << e << std::endl; + } + } + + + } @@ -232,9 +619,8 @@ void usage(char* argv[]) { - std::cerr << "usage: " << argv[0] << " input.pgm echo output.pgm" << std::endl; - std::cerr << "Laurent ISMM 2009 scheme with saliency map as output." << std::endl; - std::cerr << " echo = 0 or 1." << std::endl; + std::cerr << "usage: " << argv[0] << " input.pgm" << std::endl; + std::cerr << "Laurent ISMM 2009 scheme." << std::endl; abort(); } @@ -247,7 +633,7 @@ using value::int_u8; using value::label_16; - if (argc != 4) + if (argc != 2) usage(argv); // f: regular image. @@ -256,6 +642,8 @@ io::pgm::load(f, argv[1]); + cplx2d::e_to_labels_t e_to_labels; + // g: weights on edges. mln_VAR(g, cplx2d::f_to_g(f) ); @@ -271,8 +659,8 @@ // w: watershed labeling on edges. typedef label_16 L; // <--- Type of labels. - L n_basins; - mln_VAR( w, morpho::meyer_wst(g, nbh_g, n_basins) ); + L l_max; + mln_VAR( w, morpho::meyer_wst(g, nbh_g, l_max) ); debug::println("w:", w); @@ -299,17 +687,17 @@ util::array<A> a = labeling::compute(accu::meta::count(), g, // image of values w, // image of labels - n_basins); + l_max); + + util::array<L> ls = sort_by_increasing_attributes(a, l_max); - util::array<L> l = sort_by_increasing_attributes(a, n_basins); { - /* - std::cout << "l:" << std::endl; - for (unsigned i = 0; i < l.nelements(); ++i) - std::cout << l[i] << "(" << a[l[i]] << ") "; - std::cout << std::endl; - */ + std::cout << "ls:" << std::endl; + for (unsigned i = 1; i <= l_max; ++i) + std::cout << ls[i] << "(" << a[ls[i]] << ") "; + std::cout << std::endl + << std::endl; } @@ -320,66 +708,13 @@ // for_all(e) // if (w(e) == 0) // { -// cplx2d::e_to_labels(w, e, l1, l2); + // e_to_labels(w, e, l1, l2); // std::cout << e << ':' << l1 << ',' << l2 << std::endl; // } // } - - // Edges -> Priority queue. - - typedef p_priority< T, p_queue<E> > Q; - util::array<Q> q(n_basins.next()); - - { - L l1, l2; - mln_piter_(g_t) e(g.domain()); - for_all(e) - if (w(e) == 0) - { - cplx2d::e_to_labels(w, e, // input - l1, l2); // output - q[l1].push(mln_max(T) - g(e), e); - q[l2].push(mln_max(T) - g(e), e); - } - } - - - // Information "label l -> edge e". - - E null = E(0,0); // Impossible value. - - util::array<E> edge(n_basins.next()); - for (L l = 0; l <= n_basins; ++l) - edge[l] = null; - - - // Initialization. - - util::array<L> lpar(n_basins.next()); - for (L l = 0; l <= n_basins; ++l) - lpar[l] = l; // Make-Set. - - - - mln_ch_value_(g_line_t, E) - epar, // Edge forest. - z_epar; // Auxiliary data: edge forest with compression and balancing. - - { - initialize(epar, g_line); - initialize(z_epar, g_line); - mln_piter_(g_line_t) e(g_line.domain()); - for_all(e) - { - // Make-Set. - epar(e) = e; - z_epar(e) = e; - } - debug::println("all edges:", epar); // epar(e) == e so we depict the edges! - } - - + compute_pseudo_tree(w, g, ls, a, + e_to_labels); } Index: theo/esiee/laurent/ismm09/lca.hh --- theo/esiee/laurent/ismm09/lca.hh (revision 0) +++ theo/esiee/laurent/ismm09/lca.hh (revision 0) @@ -0,0 +1,173 @@ + +namespace mln +{ + + + // From parent image to children image. + + template <typename I, typename E, typename L> + mln_ch_value(I, std::vector<mln_psite(I)>) + compute_children(const I& epar, const std::vector<E>& edge, L l_max, std::vector<E>& roots) + { + typedef std::vector<mln_psite(I)> C; // Children. + mln_ch_value(I,C) chl; + initialize(chl, epar); + + mln_ch_value(I,bool) deja_vu; + initialize(deja_vu, epar); + data::fill(deja_vu, false); + + for (L l = 1; l <= l_max; ++l) + { + E e = edge[l]; + while (deja_vu(e) == false) + { + deja_vu(e) = true; + if (epar(e) != e) + { + chl(epar(e)).push_back(e); + e = epar(e); + } + else + roots.push_back(e); + } + } + + return chl; + } + + + /// Compute an LCA + /// Reference article: The LCA Problem Revisited, M. A. Bender and M. + /// Farach-Colton, May 2000 + /// + /// See also Topological Watershed, Alexandre Abraham, LRDE CSI report + /// 0821, June 2008 + /// + template <typename L, typename I, typename E> + class lca_t + { + + public: + + /// Constructor + lca_t(const L& l_max_, + const I& edge_children, + const std::vector<E>& roots) + { + unsigned l_max = l_max_.next(); + euler_tour_edge_.resize(2 * l_max + 1); + euler_tour_depth_.resize(2 * l_max + 1); + initialize (euler_position_, edge_children); + + // BUILD EULER TOUR + std::stack<E> to_treat; + for (unsigned i = 0; i < roots.size(); ++i) + to_treat.push(roots[i]); + + unsigned euler_edge = 0; + mln_ch_value(I, bool) deja_vu; + initialize(deja_vu, edge_children); + data::fill(deja_vu, false); + + while (!to_treat.empty()) + { + E e = to_treat.top(); + to_treat.pop(); + + euler_tour_edge_[euler_edge] = e; + if (deja_vu(e)) + euler_tour_depth_[euler_edge] = euler_tour_depth_[euler_edge - 1] - 1; + else + euler_tour_depth_[euler_edge] = euler_tour_depth_[euler_edge - 1] + 1; + + if (!deja_vu(e)) + { + for (int j = edge_children(e).size() - 1; j >= 0; --j) + { + to_treat.push(e); + to_treat.push(edge_children(e)[j]); + } + deja_vu(e) = true; + euler_position_(e) = euler_edge; + } + + ++euler_edge; + } + + // BUILD MINIMAS + int size = 2 * l_max - 1; + int logn = (int)(ceil(log((double)(size))/log(2.0))); + + minim_ = new int [logn * size]; // FIXME : Think about freeing this + Minim_ = new int* [logn]; + + Minim_[0] = minim_; + + // Init + for (unsigned i = 0; i < size - 1; ++i) + if (euler_tour_depth_[i] < euler_tour_depth_[i+1]) + { + Minim_[0][i] = i; + } else { + Minim_[0][i] = i+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) + { + if (euler_tour_depth_[Minim_[j - 1][i]] <= euler_tour_depth_[Minim_[j - 1][i + k]]) + Minim_[j][i] = Minim_[j - 1][i]; + else + Minim_[j][i] = Minim_[j - 1][i + k]; + } + } + } + } + + + /// Destructor + ~lca_t() + { + delete[] Minim_; + delete[] minim_; + } + + + const E& operator() (const E& a, const E& b) + { + if (a == b) + return a; + + int m = euler_position_(a), + n = euler_position_(b), + k; + + if (m > n) + std::swap(m, n); + + k = (int)(log((double)(n - m))/log(2.)); + + if (euler_tour_depth_[Minim_[k][m]] + < euler_tour_depth_[Minim_[k][n - (1 << k)]]) + return euler_tour_edge_[Minim_[k][m]]; + else + return euler_tour_edge_[Minim_[k][n - (1 << k)]]; + } + + private: + int *minim_; + int **Minim_; + std::vector<E> euler_tour_edge_; + std::vector<unsigned> euler_tour_depth_; + mln_ch_value(I, unsigned) euler_position_; + }; + + +} // mln Index: theo/color/segment.cc --- theo/color/segment.cc (revision 3256) +++ theo/color/segment.cc (working copy) @@ -227,6 +227,13 @@ } + value::int_u8 L_to_int_u8(unsigned l) + { + return l == 0 ? + 0 : // wshed line + 1 + (l - 1) % 255; // basin + } + } // mln @@ -297,7 +304,7 @@ // Watershed transform. - typedef value::label_8 L; + typedef value::label_16 L; L nbasins; mln_ch_value_(f_t, L) w = morpho::meyer_wst(g, e2e(), nbasins); @@ -330,7 +337,8 @@ w_all); } - io::pgm::save(w_all, "temp_w_all.pgm"); + io::pgm::save( level::transform(w_all, convert::to_fun(L_to_int_u8)), + "temp_w_all.pgm" ); }