https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox Index: ChangeLog from Thierry Geraud <thierry.geraud@lrde.epita.fr> Augment code for Laurent. * theo/color/sum_pix.hh (take): Remove useless +1. * theo/color/segment.cc: Add unactivated code. * theo/color/blen_pix.hh (sum_len): New. (take): Use it. * laurent/ismm2009.cc: UPdate. * laurent/ismm2009.v0.cc: Fix utf8. * laurent/playing_with_attributes.cc: New. laurent/ismm2009.cc | 269 ++++++++++++++++++++++++++++--------- laurent/ismm2009.v0.cc | 2 laurent/playing_with_attributes.cc | 101 +++++++++++++ theo/color/blen_pix.hh | 17 ++ theo/color/segment.cc | 18 +- theo/color/sum_pix.hh | 2 6 files changed, 334 insertions(+), 75 deletions(-) Index: theo/color/sum_pix.hh --- theo/color/sum_pix.hh (revision 3168) +++ theo/color/sum_pix.hh (working copy) @@ -125,7 +125,7 @@ inline void sum_pix<P,S>::take(const argument& p) { - s_ += 1 + p.v(); + s_ += /* 1 + */ p.v(); } template <typename P, typename S> Index: theo/color/segment.cc --- theo/color/segment.cc (revision 3168) +++ theo/color/segment.cc (working copy) @@ -459,7 +459,6 @@ if (echo) debug::println("nchildren =", nchildren | t.nodes()); - // Filtering. mln_concrete(I) g; { @@ -497,14 +496,22 @@ unsigned n_regmins_g_ref; mln_ch_value(I, unsigned) regmin_g = labeling::regional_minima(g_ref, nbh, n_regmins_g_ref); - if (echo) +// if (echo) std::cout << "n_regmins(g_ref) = " << n_regmins_g_ref << std::endl << std::endl; if (g != g_ref) + { std::cerr << "OOPS: g DIFFERS FROM ref!" << std::endl << std::endl; +// debug::println("diff", (pw::value(g_ref) == pw::value(g)) | g.domain()); + +// debug::println("g_ref", g_ref); +// debug::println("g", g); +// debug::println("regmin_g", regmin_g); + } + // bool consistency = (n_regmins_g_ref + less == n_objects); // if (consistency == false) // std::cerr << "OOPS: INCONSISTENCY (BUG...)!" << std::endl @@ -570,10 +577,10 @@ // accu::count< util::pix<I> > a_; // accu::volume<I> a_; - accu::sum_pix< util::pix<I> > a_; + // accu::sum_pix< util::pix<I> > a_; -// blen_image = input_; -// accu::blen_pix<I> a_; + blen_image = input_; + accu::blen_pix<I> a_; mln_VAR(a, compute_attribute_on_nodes(a_, t)); @@ -605,7 +612,6 @@ } io::pgm::save(w_all, "temp_w_all.pgm"); - } Index: theo/color/blen_pix.hh --- theo/color/blen_pix.hh (revision 3168) +++ theo/color/blen_pix.hh (working copy) @@ -112,6 +112,19 @@ } + template <typename B> + unsigned sum_len(const Box<B>& b_) + { + const B& b = exact(b_); + typedef mln_site(B) P; + enum { n = P::dim }; + unsigned len = 0; + for (unsigned i = 0; i < n; ++i) + len += b.len(i); + return len; + } + + # ifndef MLN_INCLUDE_ONLY template <typename I> @@ -136,7 +149,7 @@ { const mln_site(I)& p = pxl.p(); accu_blen_take(b_, p); - len_ = max_len(b_.to_result()); + len_ = sum_len(b_.to_result()); } template <typename I> @@ -145,7 +158,7 @@ blen_pix<I>::take(const blen_pix<I>& other) { b_.take(other.b()); - len_ = max_len(b_.to_result()); + len_ = sum_len(b_.to_result()); } template <typename I> Index: laurent/ismm2009.cc --- laurent/ismm2009.cc (revision 3168) +++ laurent/ismm2009.cc (working copy) @@ -71,25 +71,36 @@ } - // Get smallest edge. - template <typename Qs, typename L, typename W> - mln_deduce(Qs, element, element) - get_smallest_edge(Qs& qs, L l, const W& wst, std::vector<L>& lpar, bool& found) - { - typedef mln_element(Qs) Q; // qs is an array of queues with type Q. + // Test emptiness of the queue of a set of regions. + + template <typename Qs, typename L> + bool test_q_emptiness(Qs& qs, L l, std::vector<L>& lpar) { - // Test that, when a region has merged, its edge queue has been - // emptied. L l_ = l; while (lpar[l_] != l_) { - mln_invariant(qs[l_].is_empty()); + if (! qs[l_].is_empty()) + return false; l_ = lpar[l_]; } + return true; } + + // Get smallest edge. + + template <typename Qs, typename L, typename W> + mln_deduce(Qs, element, element) + get_smallest_edge(Qs& qs, L l, const W& wst, std::vector<L>& lpar) + { + typedef mln_element(Qs) Q; // qs 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(qs, l, lpar)); + L lr = find_root_l(lpar, l); Q& q = qs[lr]; @@ -97,7 +108,8 @@ while (! q.is_empty()) { - const E& e = q.front(); + E e = q.pop_front(); + L l1 = wst(p1_from_e(e)), l2 = wst(p2_from_e(e)); mln_invariant(l1 != l2); @@ -108,16 +120,13 @@ if (l1r == l2r) // 'e' is an internal edge => forget it. - { - q.pop(); continue; - } - found = true; return e; } - found = false; + mln_invariant(0); // We should not be here! + return E(0,0); // FIXME: null edge. } @@ -208,6 +217,115 @@ } + + + // Transform attributes so that they are all different! + + + template <typename A, typename L> + void + make_unique_attribute(util::array<A>& a, + std::vector< std::pair<A,L> >& v, + L l_max, + bool echo) + { + // Display attributes. + if (echo) + { + std::cout << "(attribute, label) = { "; + for (unsigned i = 1; i <= l_max; ++i) + std::cout << '(' << v[i].first // Attribute (increasing). + << ',' + << v[i].second // Region label. + << ") "; + std::cout << '}' << std::endl + << std::endl; + } + + std::map<A,A> lut; // old value -> new value + for (unsigned l = 1; l <= l_max; ++l) + { + A old_a = v[l].first, + new_a = old_a + l - 1; + lut[old_a] = new_a; + v[l].first = new_a; // new value + } + for (unsigned l = 1; l <= l_max; ++l) + a[l] = lut[ a[l] ]; + } + + + + // Compute g from f; then transform g into g' so that all values are + // different; last return the watershed of g'. + + + template < typename F, + typename G, + typename N_e2p, + typename N_e2e, + typename L > + mln_ch_value(G, L) + f_to_wst_unique_g(F& f, // On pixels. + G& g, // On edges. + const N_e2p& e2p, + const N_e2e& e2e, + L& n_basins, + bool echo) + { + data::paste(morpho::gradient(extend(g, f), + e2p.win()), + g); + + // FIXME: Here, we want a unique value per edge! + + if (echo) + debug::println("g (before):", g); + + { + + typedef std::pair<short, short> edge_t; + typedef std::pair<mln_value(G), edge_t> ve_t; + std::vector<ve_t> tmp; + { + mln_piter(G) e(g.domain()); + for_all(e) + { + ve_t ve; + ve.first = g(e); + ve.second = edge_t(e.row(), e.col()); + tmp.push_back(ve); + } + } + std::sort(tmp.begin(), tmp.end()); + { + mln_value(G) v; + mln_site(G) e; + for (unsigned i = 0; i < tmp.size(); ++i) + { + v = tmp[i].first + i; + e.row() = tmp[i].second.first; + e.col() = tmp[i].second.second; + g(e) = v; + } + } + } + + mln_VAR( wst_g, + morpho::meyer_wst(g, e2e, n_basins) ); + +// // Test the consistency with regional minima. +// { +// L n_regmins; +// mln_VAR( regmin_g, +// labeling::regional_minima(g, e2e, n_regmins) ); +// mln_invariant(n_basins == n_regmins); +// } + + return wst_g; + } + + } // mln @@ -229,6 +347,8 @@ using value::int_u8; using value::label_16; + bool echo = true; + if (argc != 2) usage(argv); @@ -241,7 +361,8 @@ Complex f_ = add_edges(raw_f); mln_VAR(f, f_ | is_pixel); - // debug::println("f:", f); + if (echo) + debug::println("f:", f); mln_VAR(g, f_ | is_edge); @@ -250,21 +371,34 @@ typedef label_16 L; // Type of labels. L n_basins; + mln_VAR( wst_g, f_to_wst_g(f, g, e2p(), e2e(), n_basins) ); - std::cout << "n basins = " << n_basins << std::endl; + /* + UNIQUE: + + mln_VAR( wst_g, + f_to_wst_unique_g(f, g, e2p(), e2e(), n_basins, echo) ); + */ + + if (echo) + { + std::cout << "n basins = " << n_basins << std::endl; debug::println("g:", g); debug::println("wst(g):", wst_g); - + } // Just to see things. - // + mln_VAR(adj, adjacency(wst_g, e2e())); + if (echo) debug::println("adj:", adj); + mln_VAR(adj_line, adj | (pw::value(wst_g) == pw::cst(0))); + if (echo) debug::println("adjacency:", adj_line); @@ -286,9 +420,11 @@ // Depict the watershed line. mln_VAR(is_line, pw::value(wst_g_full) == pw::cst(0)); - // --- debug::println("wst(g) line:", wst_g | is_line); + if (echo) + debug::println("wst(g) line:", wst_g | is_line); mln_VAR(g_line, g | is_line); + if (echo) debug::println("g | wst(g) line:", g_line); @@ -336,17 +472,6 @@ // if l1r = l2r. -// // Test the 'get_smallest_edge' routine. -// { -// bool found; -// for (L l = 1; l <= n_basins; ++l) -// { -// E e = get_smallest_edge(q, l, wst_g_full, lpar, found); -// std::cout << l << ": e = " << e -// << " -- g(e) = " << g(e) << std::endl; -// } -// } - E null = E(0,0); // Impossible value. @@ -376,7 +501,15 @@ std::vector<pair_t> v = sort_attributes(a, n_basins); // Sort regions. + + // Maybe activate to test purpose: + /* + make_unique_attribute(a, v, n_basins, echo); + */ + + // Display attributes. + if (echo) { std::cout << "(attribute, label) = { "; for (unsigned i = 1; i <= n_basins; ++i) @@ -408,6 +541,7 @@ epar(e) = e; z_epar(e) = e; } + if (echo) debug::println("all edges:", epar); // epar(e) == e so we depict the edges! } @@ -421,19 +555,13 @@ data::fill(aa, 0); // Safety iff 0 is an invalidate attribute value! - for (unsigned i = 1; i <= n_basins; ++i) - { - L l = v[i].second; // Region label. - bool found; - E e = get_smallest_edge(q, l, wst_g_full, lpar, found); + // The last iteration (i == n_basins) is useless: all regions have + // already merged. - if (! found) + for (unsigned i = 1; i < n_basins; ++i) { - // The last iteration is useless: all regions have already - // merged. We could have written: "for i = 1 to n_basins - 1". - mln_invariant(i == n_basins); - break; - } + L l = v[i].second; // Region label. + E e = get_smallest_edge(q, l, wst_g_full, lpar); L l1 = wst_g_full(p1_from_e(e)), l2 = wst_g_full(p2_from_e(e)), @@ -568,47 +696,51 @@ - // Displaying the "edge tree". + // About the "edge tree" and attributes. { - p_set<E> leaves; - - // Region l -> edge e. - - std::cout << "region:(edge) = "; + p_set<E> leaves; for (L l = 1; l <= n_basins; ++l) { mln_invariant(edge[l] != null); leaves.insert(edge[l]); - std::cout << l << ':' << edge[l] << " "; } - std::cout << std::endl - << std::endl; + // Tests. mln_piter_(p_set<E>) e(leaves); - - - // Tree "e -> epar(e)". - - std::cout << "edge tree: " << std::endl; for_all(e) { E e_ = e; do { - std::cout << e_ << " -> "; mln_invariant(aa(e_) != 0 && aa(epar(e_)) != 0); // aa is valid - mln_invariant(aa(epar(e_)) >= aa(e_)); + mln_invariant(aa(epar(e_)) >= aa(e_)); // edge parenthood goes with 'a' increasing e_ = epar(e_); } while (epar(e_) != e_); - std::cout << e_ << std::endl; } - std::cout << std::endl; + // Display. + if (echo) + { + std::cout << "region:(edge) = "; + for (L l = 1; l <= n_basins; ++l) + std::cout << l << ':' << edge[l] << " "; + std::cout << std::endl + << std::endl; - // Edge parenthood goes with 'a' increasing: + std::cout << "edge tree: " << std::endl; + for_all(e) { + E e_ = e; + do { + std::cout << e_ << " -> "; + e_ = epar(e_); + } + while (epar(e_) != e_); + std::cout << e_ << std::endl; + } + std::cout << std::endl; std::cout << "edge tree new attributes: " << std::endl; for_all(e) @@ -617,8 +749,6 @@ do { std::cout << aa(e_) << " -> "; - // Edge parenthood goes with a increasing: - mln_invariant(aa(epar(e_)) >= aa(e_)); e_ = epar(e_); } while (epar(e_) != e_); @@ -626,7 +756,10 @@ } std::cout << std::endl; - } // end of tree display. + } // end of Display. + + + } // end of About the "edge tree" and attributes. @@ -649,8 +782,11 @@ // Reminders: + if (echo) + { debug::println("wst(g) fully valuated:", wst_g_full); debug::println("g_line:", g_line); + } @@ -687,6 +823,7 @@ aa(e) = aa(e_); } + if (echo) debug::println("aa | wst(g) line:", (aa | is_edge) | is_line); // Second, processing the 0D-faces (points) of the watershed line. @@ -730,6 +867,7 @@ } + if (echo) debug::println("aa | line:", aa_line); @@ -747,16 +885,17 @@ aa(p) = a[wst_g_full(p)]; } + if (echo) debug::println("aa | basins", aa_basins); } + if (echo) debug::println("aa:", aa); - } // end of main Index: laurent/ismm2009.v0.cc --- laurent/ismm2009.v0.cc (revision 3168) +++ laurent/ismm2009.v0.cc (working copy) @@ -387,7 +387,7 @@ debug::println("adjacency:", adj | (pw::value(wst_g) == pw::cst(0))); - /* // Délire! + /* // De'lire! { box2d b = make::box2d(1,1, n_basins, n_basins); point2d null(0, 0); Index: laurent/playing_with_attributes.cc --- laurent/playing_with_attributes.cc (revision 0) +++ laurent/playing_with_attributes.cc (revision 0) @@ -0,0 +1,101 @@ +#include "ismm2009.hh" + +#include <mln/value/int_u8.hh> +#include <mln/value/label_8.hh> +#include <mln/value/label_16.hh> + +#include <mln/io/pgm/load.hh> +#include <mln/debug/println.hh> + +#include <mln/morpho/elementary/gradient.hh> +#include <mln/morpho/gradient.hh> +#include <mln/morpho/tree/data.hh> +#include <mln/morpho/tree/compute_attribute_image.hh> +#include <mln/labeling/regional_minima.hh> + +#include <mln/accu/count.hh> + + + +void usage(char* argv[]) +{ + std::cerr << "usage: " << argv[0] << " input.pgm" << std::endl; + std::cerr << "For Laurent's ISMM 2009 scheme: min-tree, attributes, and watershed." << std::endl; + abort(); +} + + + +template <typename I, typename N, typename L> +void do_it(const I& g, const N& nbh, L& n_labels) +{ + using namespace mln; + + // regional minima + + mln_ch_value(I,L) regmin = labeling::regional_minima(g, nbh, n_labels); + debug::println("regmin(g):", regmin); + + + // watershed + + L n_basins; + mln_ch_value(I,L) w = morpho::meyer_wst(g, nbh, n_basins); + mln_invariant(n_basins == n_labels); + debug::println("w(g):", w); + + + { + typedef p_array<mln_site(I)> S; + S s = level::sort_psites_decreasing(g); + + typedef morpho::tree::data<I,S> tree_t; + tree_t t(g, s, nbh); + + accu::count< util::pix<I> > a_; // Kind of attribute. + + mln_ch_value(I,unsigned) a = morpho::tree::compute_attribute_image(a_, t); + debug::println("a:", a); + + debug::println("a | w line:", a | (pw::value(w) == pw::cst(0))); + debug::println("a | w basins:", a | (pw::value(w) != pw::cst(0))); + + debug::println("a | regmin:", a | (pw::value(regmin) != pw::cst(0))); + } + +} // end of do_it + + + +int main(int argc, char* argv[]) +{ + using namespace mln; + using value::int_u8; + + + typedef image2d<int_u8> I; + typedef value::label_8 L; + + border::thickness = 0; + + if (argc != 2) + usage(argv); + + I raw_f; + io::pgm::load(raw_f, argv[1]); + + + I f_ = add_edges(raw_f); + mln_VAR(f, f_ | is_pixel); + debug::println("f:", f); + + mln_VAR(g, f_ | is_edge); + data::paste( morpho::gradient(extend(g, f_), + e2p().win()), + g ); + debug::println("g:", g); + + L n; + do_it(g, e2e(), n); + +}