3257: Update Laurent's code.

3 Feb 2009

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


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

	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" );
   }

    

Thierry Geraud

tags

participants (1)