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

19 Jan 2009

https://svn.lrde.epita.fr/svn/oln/trunk/milena/sandbox


Index: ChangeLog
from  Thierry Geraud  <thierry.geraud@lrde.epita.fr>

	Update Laurent's ISMM code.

	* laurent/ismm2009.cc: Update.
	* laurent/ismm2009.v1.cc: New.
	* laurent/ismm2009.hh (is_not_point): New.
	(find_root): Remove.

 ismm2009.cc    |  687 +++++++++++++++++++++++++++++++++++++++++++++++++---
 ismm2009.hh    |   33 --
 ismm2009.v1.cc |  745 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 1403 insertions(+), 62 deletions(-)

Index: laurent/ismm2009.cc

--- laurent/ismm2009.cc	(revision 3162)
+++ laurent/ismm2009.cc	(working copy)
@@ -2,16 +2,216 @@
 
 #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/util/ord_pair.hh>
 #include <mln/debug/println.hh>
 #include <mln/core/routine/duplicate.hh>
 
+# include <mln/core/site_set/p_set.hh>
+# include <mln/core/site_set/p_queue.hh>
+# include <mln/core/site_set/p_priority.hh>
+
 #include <mln/labeling/compute.hh>
 #include <mln/accu/count.hh>
 
 
+/*
+
+   TO-DO list:
+   -----------
+
+   - having a heap for every lr (support merge)
+
+   - handling adjacency on the fly (instead of having an image)
+
+ */
+
+
+namespace mln
+{
+
+
+  // Region label parenthood routines.
+
+
+  template <typename L>
+  inline
+  void make_sets(std::vector<L>& lpar, L l_max)
+  {
+    for (L l = 1; l <= l_max; ++l)
+      lpar[l] = l; // Make-set.
+  }
+
+  template <typename L>
+  inline
+  L find_root_l(std::vector<L>& lpar, L l)
+  {
+    if (lpar[l] == l)
+      return l;
+    else
+      return lpar[l] = find_root_l(lpar, lpar[l]);
+  }
+
+
+
+  // Edge tree.
+
+
+  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));
+  }
+
+
+  // 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 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());
+	  l_ = lpar[l_];
+	}
+    }
+
+    L lr = find_root_l(lpar, l);
+    Q& q = qs[lr];
+
+    typedef mln_element(Q) E;
+
+    while (! q.is_empty())
+      {
+	const E& e = q.front();
+	L l1 = wst(p1_from_e(e)),
+	  l2 = wst(p2_from_e(e));
+	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.
+	  {
+	    q.pop();
+	    continue;
+	  }
+
+	found = true;
+	return e;
+      }
+
+    found = false;
+    return E(0,0); // FIXME: null edge.
+  }
+
+
+  // Compute an LCA; reference(?) code.
+
+  template <typename L, typename E, typename Epar>
+  std::vector< std::vector<E> >
+  compute_lca(const L& l_max,
+	      const std::vector<E>& edge,
+	      const Epar& epar)
+  {
+    mln_ch_value(Epar, std::vector<L>) chl_; // Children (known as array).
+    {
+      initialize(chl_, epar);
+      E e;
+      for (L l = 1; l <= l_max; ++l)
+	{
+	  e = edge[l];
+	  do
+	    {
+	      chl_(e).push_back(l);
+	      e = epar(e);
+	    }
+	  while (e != epar(e));
+	  chl_(e).push_back(l);
+	}
+      // e is the root edge.
+      mln_invariant(chl_(e).size() == l_max); // All labels are children.
+    }
+
+    /*
+    // Display children tree.
+    {
+      std::cout << "chl_ tree: " << std::endl;
+      for (L l = 1; l <= l_max; ++l)
+	{
+	  E e_ = edge[l];
+	  std::cout << "l = " << l << "  =>  ";
+	  do
+	    {
+	      std::cout << e_ << " : ";
+	      for (unsigned i = 0; i < chl_(e_).size(); ++i)
+		std::cout << chl_(e_)[i] << ' ';
+	      std::cout << " ->  ";
+	      e_ = epar(e_);
+	    }
+	  while (epar(e_) != e_);
+	  std::cout << e_ << " : ";
+	  for (unsigned i = 0; i < chl_(e_).size(); ++i)
+	    std::cout << chl_(e_)[i] << ' ';
+	  std::cout << std::endl;
+	}
+    }
+    */
+
+    mln_ch_value(Epar, std::set<L>) chl; // Children (as a set).
+    {
+      initialize(chl, epar);
+    
+      mln_piter(Epar) e(chl.domain());
+      for_all(e)
+	if (chl_(e).size() != 0)
+	  chl(e).insert(chl_(e).begin(), chl_(e).end());
+    }
+
+    unsigned size = l_max.next();
+
+    std::vector< std::vector<E> > lca(size);
+    for (L l = 1; l <= l_max; ++l)
+      {
+	lca[l].resize(size);
+
+	// Diagonal.
+	lca[l][l] = edge[l]; // Itself.
+
+	// Superior right.
+	for (L l2 = l.next(); l2 <= l_max; ++l2)
+	  {
+	    E e = edge[l];
+	    while (chl(e).find(l2) == chl(e).end())
+	      e = epar(e);
+	    lca[l][l2] = e;
+	}
+      }
+
+    return lca;
+  }
+
+
+} // mln
+
+
+
 
 void usage(char* argv[])
 {
@@ -22,11 +222,12 @@
 
 
 
+
 int main(int argc, char* argv[])
 {
   using namespace mln;
   using value::int_u8;
-  using value::label_8;
+  using value::label_16;
 
   if (argc != 2)
     usage(argv);
@@ -34,27 +235,37 @@
   image2d<int_u8> raw_f;
   io::pgm::load(raw_f, argv[1]);
 
-  image2d<int_u8> f_ = add_edges(raw_f);
+  typedef image2d<int_u8> Complex;
+
+  typedef mln_pset_(Complex) full_D;
+  Complex f_ = add_edges(raw_f);
 
   mln_VAR(f, f_ | is_pixel);
   // debug::println("f:", f);
 
   mln_VAR(g, f_ | is_edge);
 
-  typedef label_8 L; // Type of labels.
+  typedef mln_value_(g_t) T; // Type of edge values.
+  typedef mln_psite_(g_t) E; // Type of edges.
+  typedef label_16 L;        // Type of labels.
 
   L n_basins;
   mln_VAR( wst_g,
 	   f_to_wst_g(f, g, e2p(), e2e(), n_basins) );
 
-  // debug::println("g:", g);
+  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()));
-  debug::println("adjacency:", adj | (pw::value(wst_g) == pw::cst(0)));
+  debug::println("adj:", adj);
+  mln_VAR(adj_line, adj | (pw::value(wst_g) == pw::cst(0)));
+  debug::println("adjacency:", adj_line);
 
 
 
@@ -71,33 +282,89 @@
 				c4().win()),
 		wst_g_full);
   }
-  debug::println("wst(g) fully valuated:", wst_g_full);
 
 
   // 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);
+  // --- debug::println("wst(g) line:", wst_g | is_line);
 
+  mln_VAR(g_line, g | is_line);
+  debug::println("g | wst(g) line:", g_line);
 
 
-  // Get the smallest edge out of every basin.
 
-  std::vector<point2d> edge = smallest_edges(extend(wst_g | is_line, wst_g_full),
-					     n_basins, e2p(), g);
-  for (L l = 1; l <= n_basins; ++l)
-    std::cout << int(l) << ": " << edge[l] << std::endl;
 
 
+  // Priority queue -> edge.
+
+  typedef p_priority< T, p_queue<E> > Q;
+  util::array<Q> q(n_basins.next());
+
+  {
+    mln_piter_(g_t) e(g.domain());
+    for_all(e)
+      if (wst_g(e) == 0)
+	{
+	  L l1 = wst_g_full(p1_from_e(e)),
+	    l2 = wst_g_full(p2_from_e(e));
+	  if (l1 > l2)
+	    std::swap(l1, l2);
+	  mln_invariant(l1 < l2);
+	  q[l1].push(mln_max(T) - g(e), e);
+	  q[l2].push(mln_max(T) - g(e), e);
+	}
+    // --- for (L l = 1; l <= n_basins; ++l)
+    // ---   std::cout << "q["<< l << "] = " << q[l] << std::endl;
+  }
+
+
+
+
+
+  // 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.
+  std::vector<L> lpar(n_basins.next());
+  make_sets(lpar, n_basins);
+
+
+  // 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.
+
+
+//   // Test the 'get_smallest_edge' routine.
 //   {
-//     // Test the result with an alternate code.
-//     std::vector<point2d> edge_alt = smallest_edges_alt(wst_g, n_basins, e2e(), g);
+//     bool found;
 //     for (L l = 1; l <= n_basins; ++l)
-//       mln_invariant(edge_alt[l] == edge[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.
+
+
+  // Information "label l -> edge e".
+
+  std::vector<E> edge(n_basins.next());
+  for (L l = 1; l <= n_basins; ++l)
+    edge[l] = null;
+
+
+
+
 
 
 
   // Compute an attribute per region.
+  // --------------------------------
 
   typedef unsigned A;
   util::array<A> a = labeling::compute(accu::meta::count(),
@@ -109,49 +376,385 @@
   std::vector<pair_t> v = sort_attributes(a, n_basins); // Sort regions.
 
 
-  std::cout << "attributes =  ";
+  // Display attributes.
+  {
+    std::cout << "(attribute, label)  =  { ";
   for (unsigned i = 1; i <= n_basins; ++i)
-    std::cout << v[i].first   // Attribute (increasing).
-	      << ':'
+      std::cout << '(' << v[i].first   // Attribute (increasing).
+		<< ','
 	      << v[i].second  // Region label.
-	      << " - ";
-  std::cout << std::endl;
+		<< ") ";
+    std::cout << '}' << std::endl
+	      << std::endl;
+  }
 
 
-  std::vector<L> lpar(n_basins.next());
-  for (L l = 1; l <= n_basins; ++l)
-    lpar[l] = l; // Make-set.
 
 
-  util::array<A> a_merged = a;
+  // Go!
+  // --------------------------------
+
+
+  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!
+  }
+
+
+
+  // 'aa' is the result attribute image; it is defined on the complex
+  // (so it has edges, pixels, and 0-face-points).
+
+  image2d<A> aa;
+  initialize(aa, f_);
+  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.
-	lr = find_root(lpar, l);
+      L l = v[i].second; // Region label.
+      bool found;
+      E e = get_smallest_edge(q, l, wst_g_full, lpar, found);
+
+      if (! found)
+	{
+	  // 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;
+	}
 
-      point2d e = edge[l]; // FIXME: Use the heap!
+      L l1 = wst_g_full(p1_from_e(e)),
+	l2 = wst_g_full(p2_from_e(e)),
+	l1r = find_root_l(lpar, l1),
+	l2r = find_root_l(lpar, l2);
 
-      mln_invariant(wst_g_full(p1_from_e(e)) == l ||
-		    wst_g_full(p2_from_e(e)) == l);
-      L l2 = ( wst_g_full(p1_from_e(e)) == l ?
-	       wst_g_full(p2_from_e(e)) :
-	       wst_g_full(p1_from_e(e)) ),
-	l2r = find_root(lpar, l2);
-
-      if (lr == l2r)
-	continue; // Already merged.
-      if (l2r < lr)
-	std::swap(lr, l2r);
-      mln_invariant(l2r > lr);
-      lpar[lr] = l2r;
-      a_merged[l2r] += lr; // FIXME: We want accumulators here!
+      aa(e) = a[l];
+
+      // Consistency tests.
+      {
+	mln_invariant(a[l] == v[i].first);
+	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)));
     }
 
+      /*
+      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 <= n_basins; ++i)
     {
       L l = v[i].second;
-      std::cout << l << " -> " << lpar[l] << std::endl;
+	      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"
+
+
+
+  // Displaying the "edge tree".
+
+  {
+    p_set<E> leaves;
+
+
+    // Region l -> edge e.
+
+    std::cout << "region:(edge)  =  ";
+    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;
+
+    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_));
+	  e_ = epar(e_);
+	}
+      while (epar(e_) != e_);
+      std::cout << e_ << std::endl;
+    }
+    std::cout << std::endl;
+
+
+    // Edge parenthood goes with 'a' increasing:
+
+    std::cout << "edge tree new attributes: " << std::endl;
+    for_all(e)
+    {
+      E e_ = e;
+      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_);
+      std::cout << aa(e_) << std::endl;
+    }
+    std::cout << std::endl;
+
+  } // end of tree display.
+
+
+
+  // Testing both that all regions and all "smallest" edges have
+  // merged (assumption: connected domain!).
+
+  {
+    L l_1 = v[1].second,
+      l_root = find_root_l(lpar, l_1);
+    mln_invariant(edge[l_1] != null);
+    E e_root = find_root_e(z_epar, edge[l_1]);
+    for (unsigned i = 2; i <= n_basins; ++i)
+      {
+	L l = v[i].second;
+ 	mln_invariant(find_root_l(lpar, l) == l_root);
+	mln_invariant(edge[l] != null);
+  	mln_invariant(find_root_e(z_epar, edge[l]) == e_root);
+     }
+  }
+
+
+  // Reminders:
+  debug::println("wst(g) fully valuated:", wst_g_full);
+  debug::println("g_line:", g_line);
+
+
+
+  // +---------------------------------------------------------------+
+  // | We want an "image" of the "merging value of a" over all edges |
+  // | (not only the set of 'e' used to merge).                      |
+  // +---------------------------------------------------------------+
+
+
+  {
+
+    // Dealing with the watershed line.
+
+    mln_VAR(aa_line, aa | is_line);
+
+    {
+      // First, processing the 1D-faces (edges) of the watershed line.
+
+      std::vector< std::vector<E> > lca;
+      lca = compute_lca(n_basins, edge, epar); // lca is auxiliary data.
+
+      mln_piter_(g_line_t) e(g_line.domain());
+      for_all(e)
+      {
+	L l1 = adj_line(e).first(),
+	  l2 = adj_line(e).second();
+	mln_invariant(l1 != 0 && l2 > l1);
+	E e_ = lca[l1][l2];
+	mln_invariant(is_line(e_));
+	mln_invariant(aa(e_) != 0); // aa is valid at e_.
+
+	// The attribute value propagates from the lca to the current edge
+	// of the line:
+	aa(e) = aa(e_);
+      }
+
+      debug::println("aa | wst(g) line:", (aa | is_edge) | is_line);
+
+      // Second, processing the 0D-faces (points) of the watershed line.
+
+      data::paste(morpho::dilation(extend(aa_line | (pw::value(aa_line) == pw::cst(0)),
+					  aa_line),
+				   c4().win()),
+		  aa_line);
+
+      /*
+
+	// Version with "aa(x) = 0" to separate sheds (a shed =
+	// adjacency between a couple of basins).
+
+      mln_piter_(aa_line_t) pxl(aa_line.domain());
+      mln_niter_(neighb2d) ne(c4(), pxl); // FIXME: "pxl -> nbhoring edges" is explictly c4()...
+      for_all(pxl)
+      {
+	if (aa(pxl) != 0)
+	  continue;
+	unsigned count = 0;
+	A na, na2;
+	for_all(ne)
+	  if (aa_line.has(ne))
+	    {
+	      mln_invariant(aa(ne) != 0);
+	      if (count == 0)
+		na  = aa(ne); // First attribute value encountered.
+	      else
+		na2 = aa(ne); // Second (or more) attr value encountered.
+	      ++count;
+	    }
+	if (count == 2)
+	  {
+	    mln_invariant(na == na2); // Both attribute values have to be the same.
+	    aa(pxl) = na;
+	  }
+      }
+      */
+
+
+    }
+
+    debug::println("aa | line:", aa_line);
+
+
+    // Now dealing with all elements of basins:
+
+    //   -  2D-faces, i.e., original pixels;
+    //   -  1D-faces, i.e., edges within regions;
+    //   -  0D-faces, i.e., points within regions.
+
+    mln_VAR(aa_basins, aa | (pw::value(wst_g_full) != pw::cst(0)));
+
+    {
+      mln_piter_(aa_basins_t) p(aa_basins.domain());
+      for_all(p)
+	aa(p) = a[wst_g_full(p)];
+    }
+
+    debug::println("aa | basins", aa_basins);
+
+  }
+
+
+
+  debug::println("aa:", aa);
+
 
 
 } // end of main
Index: laurent/ismm2009.v1.cc
--- laurent/ismm2009.v1.cc	(revision 0)
+++ laurent/ismm2009.v1.cc	(revision 0)
@@ -0,0 +1,745 @@
+#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/util/ord_pair.hh>
+#include <mln/debug/println.hh>
+#include <mln/core/routine/duplicate.hh>
+
+# include <mln/core/site_set/p_set.hh>
+# include <mln/core/site_set/p_queue.hh>
+# include <mln/core/site_set/p_priority.hh>
+
+#include <mln/labeling/compute.hh>
+#include <mln/accu/count.hh>
+
+
+
+namespace mln
+{
+
+
+  // Region label parenthood routines.
+
+
+  template <typename L>
+  inline
+  void make_sets(std::vector<L>& lpar, L l_max)
+  {
+    for (L l = 1; l <= l_max; ++l)
+      lpar[l] = l; // Make-set.
+  }
+
+  template <typename L>
+  inline
+  L find_root_l(std::vector<L>& lpar, L l)
+  {
+    if (lpar[l] == l)
+      return l;
+    else
+      return lpar[l] = find_root_l(lpar, lpar[l]);
+  }
+
+
+
+  // Edge tree.
+
+
+  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));
+  }
+
+
+  // Get smallest edge.
+
+  template <typename Q, typename L, typename W>
+  mln_element(Q)
+  get_smallest_edge(Q& q, L l, const W& wst, std::vector<L>& lpar, bool& found)
+  {
+    L lr = find_root_l(lpar, l);
+    typedef mln_element(Q) E;
+
+    mln_fwd_piter(Q) e(q);
+    for_all(e)
+    {
+      L l1 = wst(p1_from_e(e)),
+	l2 = wst(p2_from_e(e));
+      mln_invariant(l1 != l2);
+      L l1r = find_root_l(lpar, l1),
+	l2r = find_root_l(lpar, l2);
+      
+      if (l1r == l2r)
+	// 'e' is an internal edge => forget it.
+	continue;
+      
+      if (l1r == lr || l2r == lr)
+	{
+	  found = true;
+	  return e;
+	}
+    }
+
+    found = false;
+    return E();
+  }
+
+
+  // Compute an LCA; reference(?) code.
+
+  template <typename L, typename E, typename Epar>
+  image2d<E> compute_lca(const L& l_max,
+			 const std::vector<E>& edge,
+			 const Epar& epar)
+  {
+    mln_ch_value(Epar, std::vector<L>) chl_; // Children (known as array).
+    {
+      initialize(chl_, epar);
+      E e;
+      for (L l = 1; l <= l_max; ++l)
+	{
+	  e = edge[l];
+	  do
+	    {
+	      chl_(e).push_back(l);
+	      e = epar(e);
+	    }
+	  while (e != epar(e));
+	  chl_(e).push_back(l);
+	}
+      // e is the root edge.
+      mln_invariant(chl_(e).size() == l_max); // All labels are children.
+    }
+
+    /*
+    // Display children tree.
+    {
+      std::cout << "chl_ tree: " << std::endl;
+      for (L l = 1; l <= l_max; ++l)
+	{
+	  E e_ = edge[l];
+	  std::cout << "l = " << l << "  =>  ";
+	  do
+	    {
+	      std::cout << e_ << " : ";
+	      for (unsigned i = 0; i < chl_(e_).size(); ++i)
+		std::cout << chl_(e_)[i] << ' ';
+	      std::cout << " ->  ";
+	      e_ = epar(e_);
+	    }
+	  while (epar(e_) != e_);
+	  std::cout << e_ << " : ";
+	  for (unsigned i = 0; i < chl_(e_).size(); ++i)
+	    std::cout << chl_(e_)[i] << ' ';
+	  std::cout << std::endl;
+	}
+    }
+    */
+
+    mln_ch_value(Epar, std::set<L>) chl; // Children (as a set).
+    {
+      initialize(chl, epar);
+    
+      mln_piter(Epar) e(chl.domain());
+      for_all(e)
+	if (chl_(e).size() != 0)
+	  chl(e).insert(chl_(e).begin(), chl_(e).end());
+    }
+
+
+    box2d b = make::box2d(1, 1, l_max, l_max);
+    image2d<E> lca(b);
+    data::fill(lca, E(0,0));
+    // Superior right.
+    for (int l1 = 1; l1 <= l_max; ++l1)
+      for (int l2 = l1 + 1; l2 <= l_max; ++l2)
+	{
+	  E e = edge[l1];
+	  while (chl(e).find(l2) == chl(e).end())
+	    e = epar(e);
+	  lca.at_(l1, l2) = e;
+	}
+    // Diagonal.
+    for (int l = 1; l <= l_max; ++l)
+      lca.at_(l, l) = edge[l]; // Itself.
+
+    // debug::println(lca);
+
+    return lca;
+  }
+
+
+} // mln
+
+
+
+
+void usage(char* argv[])
+{
+  std::cerr << "usage: " << argv[0] << " input.pgm" << std::endl;
+  std::cerr << "Laurent ISMM 2009 scheme." << std::endl;
+  abort();
+}
+
+
+
+
+int main(int argc, char* argv[])
+{
+  using namespace mln;
+  using value::int_u8;
+  using value::label_16;
+
+  if (argc != 2)
+    usage(argv);
+
+  image2d<int_u8> raw_f;
+  io::pgm::load(raw_f, argv[1]);
+
+  typedef image2d<int_u8> Complex;
+
+  typedef mln_pset_(Complex) full_D;
+  Complex f_ = add_edges(raw_f);
+
+  mln_VAR(f, f_ | is_pixel);
+  // debug::println("f:", f);
+
+  mln_VAR(g, f_ | is_edge);
+
+  typedef mln_value_(g_t) T; // Type of edge values.
+  typedef mln_psite_(g_t) E; // Type of edges.
+  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;
+
+  debug::println("g:", g);
+  debug::println("wst(g):", wst_g);
+
+
+
+  // Just to see things.
+  //
+  mln_VAR(adj, adjacency(wst_g, e2e()));
+  debug::println("adj:", adj);
+  mln_VAR(adj_line, adj | (pw::value(wst_g) == pw::cst(0)));
+  debug::println("adjacency:", adj_line);
+
+
+
+  image2d<L> wst_g_full = wst_g.unmorph_();
+  {
+    // edges (1D-faces) -> pixels (2D-faces)
+    mln_VAR(w_pixels, wst_g_full | is_pixel);
+    data::paste(morpho::dilation(extend(w_pixels, pw::value(wst_g_full)),
+				 c4().win()),
+		wst_g_full);
+    // edges (1D-faces) -> points (0D-faces)
+    mln_VAR(w_points, wst_g_full | is_point);
+    data::paste(morpho::erosion(extend(w_points, pw::value(wst_g_full)),
+				c4().win()),
+		wst_g_full);
+  }
+
+
+  // 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);
+
+  mln_VAR(g_line, g | is_line);
+  debug::println("g | wst(g) line:", g_line);
+
+
+
+
+
+  // Priority queue -> edge.
+
+  typedef p_priority< T, p_queue<E> > Q;
+  Q q;
+
+  {
+    mln_piter_(g_t) e(g.domain());
+    for_all(e)
+      if (wst_g(e) == 0)
+	q.push(mln_max(T) - g(e), e);
+    // --- std::cout << "q = " << q << std::endl;
+  }
+
+  // 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.
+  std::vector<L> lpar(n_basins.next());
+  make_sets(lpar, n_basins);
+
+  // 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.
+
+
+//   // 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.
+
+
+  // Information "label l -> edge e".
+
+  std::vector<E> edge(n_basins.next());
+  for (L l = 1; l <= n_basins; ++l)
+    edge[l] = null;
+
+
+
+
+
+
+
+  // Compute an attribute per region.
+  // --------------------------------
+
+  typedef unsigned A;
+  util::array<A> a = labeling::compute(accu::meta::count(),
+				       g,
+				       wst_g,
+				       n_basins);
+
+  typedef std::pair<A,L> pair_t;
+  std::vector<pair_t> v = sort_attributes(a, n_basins); // Sort regions.
+
+
+  // Display attributes.
+  {
+    std::cout << "(attribute, label)  =  { ";
+    for (unsigned i = 1; i <= n_basins; ++i)
+      std::cout << '(' << v[i].first   // Attribute (increasing).
+		<< ','
+		<< v[i].second  // Region label.
+		<< ") ";
+    std::cout << '}' << std::endl
+	      << std::endl;
+  }
+
+
+
+
+  // Go!
+  // --------------------------------
+
+
+  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!
+  }
+
+
+
+  // 'aa' is the result attribute image; it is defined on the complex
+  // (so it has edges, pixels, and 0-face-points).
+
+  image2d<A> aa;
+  initialize(aa, f_);
+  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);
+
+      if (! found)
+	{
+	  // 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 l1 = wst_g_full(p1_from_e(e)),
+	l2 = wst_g_full(p2_from_e(e)),
+	l1r = find_root_l(lpar, l1),
+	l2r = find_root_l(lpar, l2);
+
+      aa(e) = a[l];
+
+      // Consistency tests.
+      {
+	mln_invariant(a[l] == v[i].first);
+	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)));
+      }
+
+      /*
+      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;
+
+	/*
+	// Displaying lpar.
+	{
+	  std::cout << "lpar =  ";
+	  for (unsigned i = 1; i <= n_basins; ++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"
+
+
+
+  // Displaying the "edge tree".
+
+  {
+    p_set<E> leaves;
+
+
+    // Region l -> edge e.
+
+    std::cout << "region:(edge)  =  ";
+    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;
+
+    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_));
+	  e_ = epar(e_);
+	}
+      while (epar(e_) != e_);
+      std::cout << e_ << std::endl;
+    }
+    std::cout << std::endl;
+
+
+    // Edge parenthood goes with 'a' increasing:
+
+    std::cout << "edge tree new attributes: " << std::endl;
+    for_all(e)
+    {
+      E e_ = e;
+      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_);
+      std::cout << aa(e_) << std::endl;
+    }
+    std::cout << std::endl;
+
+  } // end of tree display.
+
+
+
+  // Testing both that all regions and all "smallest" edges have
+  // merged (assumption: connected domain!).
+
+  {
+    L l_1 = v[1].second,
+      l_root = find_root_l(lpar, l_1);
+    mln_invariant(edge[l_1] != null);
+    E e_root = find_root_e(z_epar, edge[l_1]);
+    for (unsigned i = 2; i <= n_basins; ++i)
+      {
+	L l = v[i].second;
+ 	mln_invariant(find_root_l(lpar, l) == l_root);
+	mln_invariant(edge[l] != null);
+  	mln_invariant(find_root_e(z_epar, edge[l]) == e_root);
+     }
+  }
+
+
+  // Reminders:
+  debug::println("wst(g) fully valuated:", wst_g_full);
+  debug::println("g_line:", g_line);
+
+
+
+  // +---------------------------------------------------------------+
+  // | We want an "image" of the "merging value of a" over all edges |
+  // | (not only the set of 'e' used to merge).                      |
+  // +---------------------------------------------------------------+
+
+
+  {
+
+    // Dealing with the watershed line.
+
+    mln_VAR(aa_line, aa | is_line);
+
+    {
+      // First, processing the 1D-faces (edges) of the watershed line.
+
+      image2d<E> lca = compute_lca(n_basins, edge, epar); // lca is auxiliary data.
+
+      mln_piter_(g_line_t) e(g_line.domain());
+      for_all(e)
+      {
+	E e_ = lca.at_(adj_line(e).first(), adj_line(e).second());
+	mln_invariant(is_line(e_));
+	mln_invariant(aa(e_) != 0); // aa is valid at e_.
+
+	// The attribute value propagates from the lca to the current edge
+	// of the line:
+	aa(e) = aa(e_);
+      }
+
+      debug::println("aa | wst(g) line:", (aa | is_edge) | is_line);
+
+      // Second, processing the 0D-faces (points) of the watershed line.
+
+      mln_piter_(aa_line_t) pxl(aa_line.domain());
+      mln_niter_(neighb2d) ne(c4(), pxl); // FIXME: "pxl -> nbhoring edges" is explictly c4()...
+      for_all(pxl)
+      {
+	if (aa(pxl) != 0)
+	  continue;
+	unsigned count = 0;
+	A na, na2;
+	for_all(ne)
+	  if (aa_line.has(ne))
+	    {
+	      mln_invariant(aa(ne) != 0);
+	      if (count == 0)
+		na  = aa(ne); // First attribute value encountered.
+	      else
+		na2 = aa(ne); // Second (or more) attr value encountered.
+	      ++count;
+	    }
+	if (count == 2)
+	  {
+	    mln_invariant(na == na2); // Both attribute values have to be the same.
+	    aa(pxl) = na;
+	  }
+      }
+
+    }
+
+    debug::println("aa | line:", aa_line);
+
+
+    // Now dealing with all elements of basins:
+
+    //   -  2D-faces, i.e., original pixels;
+    //   -  1D-faces, i.e., edges within regions;
+    //   -  0D-faces, i.e., points within regions.
+
+    mln_VAR(aa_basins, aa | (pw::value(wst_g_full) != pw::cst(0)));
+
+    {
+      mln_piter_(aa_basins_t) p(aa_basins.domain());
+      for_all(p)
+	aa(p) = a[wst_g_full(p)];
+    }
+
+    debug::println("aa | basins", aa_basins);
+
+  }
+
+
+
+  debug::println("aa:", aa);
+
+
+
+} // end of main
+
+
+
+
+/*
+
+
+for (i=0; i<V; i++) // Initialiser les heaps de fibonacci
+{
+  fh_setcmp(G->hp[i], cmp); // Chaque region a une heap de ses edges
+}
+
+forall edges z that separates two regions v and w
+{
+  fh_insert(G->hp[v], (void *)(z));   // Ajouter les edges dans les heaps
+  fh_insert(G->hp[w], (void *)(z));
+}
+
+UFinit(G->V); // Initialiser l'union-find
+
+// Parcourir les regions par attribut croissant
+for (j=0; j<V-1; j++)
+{
+  i = find(j);
+
+  do
+    { // trouver l'edge minimum qui sorte de la region
+      e = fh_extractmin(G->hp[i]);
+    }
+  while ((UFfind(e->v, e->w)) && (e !=NULL));
+
+  if (e != NULL)
+    { // Normalement, e != NULL, sinon c'est un BIG pb!!!
+      int ui, uj, uk;
+      ui = find(e->v);
+      uj = find(e->w);
+      uk = UFunion(e->v,e->w); // Merger les regions
+      if (uk == ui)
+	{ // et merger les edges
+	  G->hp[ui] = fh_union(G->hp[ui], G->hp[uj]);
+	}
+      else
+	{
+	  G->hp[uj] = fh_union(G->hp[uj], G->hp[ui]);
+	}
+      mst[k] = *e; // Garder l'arete
+      SaliencyWeight[k] = attribut[ui];// Poids dans la nouvelle hierarchie
+      OldWeight[k] = e->weight; // Poids dans l'ancienne hierarchie (inutile)
+      k++;
+    }
+
+  // Calcul de la sortie
+  Pour toutes les edges separantes z=(x,y)
+    {
+	S[z] = max {SaliencyWeight[k] | sur le chemin reliant x a y dans mst}
+    }
+}
+
+
+ */
Index: laurent/ismm2009.hh
--- laurent/ismm2009.hh	(revision 3162)
+++ laurent/ismm2009.hh	(working copy)
@@ -49,6 +49,12 @@
     return p.row() % 2 && p.col() % 2;
   }
   
+  inline
+  bool is_not_point(const point2d& p)
+  {
+    return ! is_point(p);
+  }
+  
 
   // Neighborhoods.
 
@@ -268,19 +274,6 @@
   }
 
 
-  // Find root.
-
-  template <typename L>
-  inline
-  L find_root(std::vector<L>& par, L l)
-  {
-    if (par[l] == l)
-      return l;
-    else
-      return par[l] = find_root(par, par[l]);
-  }
-
-
 
   // Display.
 
@@ -334,13 +327,13 @@
     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);
-    }
+//     // 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;
   }

    