Clean up g2 computation and modularisation.

      *  apps/g2/compute_ctos.cpp,
      *  apps/g2/compute_ctos.hpp: New.
      *  apps/g2/g2-maxdepth.cpp: Support for other pinf.
      *  apps/g2/satmaxtree-cli.cpp: Externalize 2-Face removal.
      *  apps/g2/satmaxtree.cpp,
      *  apps/g2/satmaxtree.hpp: Clean facades and add 2F removal method.
      *  mln/morpho/tos/ctos.hpp: Add some overloads.

apps/g2/compute_ctos.cpp

+#include <mln/core/neighb2d.hpp>
+#include <mln/core/dontcare.hpp>
+#include <mln/morpho/tos/ctos.hpp>
+#include <mln/morpho/component_tree/reconstruction.hpp>
+
+
+#include <apps/tos/addborder.hpp>
+#include <apps/tos/Kinterpolate.hpp>
+#include <apps/tos/topology.hpp>
+#include "compute_g2.hpp"
+#include "satmaxtree.hpp"
+
+#include <boost/graph/dag_shortest_paths.hpp>
+#include <boost/graph/transpose_graph.hpp>
+#include <boost/property_map/function_property_map.hpp>
+#include <tbb/parallel_for.h>
+
+namespace mln
+{
+
+  template <class MyGraph>
+  boost::vector_property_map<unsigned>
+  compute_graph_depth(const MyGraph& g)
+  {
+    mln_entering("Compute graph depth");
+
+    boost::vector_property_map<unsigned> depth(boost::num_vertices(g));
+
+    auto one = [](mln::dontcare_t) -> int{ return 1; };
+    auto w = boost::make_function_property_map<typename MyGraph::edge_descriptor, int, decltype(one)>(one);
+
+    typename MyGraph::vertex_descriptor root = boost::vertex(0, g);
+    depth[root] = 0;
+
+    MyGraph gT;
+    boost::transpose_graph(g, gT);
+    boost::dag_shortest_paths(gT, root, boost::weight_map(w)
+                              .distance_map(depth)
+                              .distance_compare(std::greater<int> ())
+                              .distance_inf(-1)
+                              .distance_zero(0)
+                            );
+
+
+    mln_exiting();
+    return depth;
+  }
+
+  template <class MyGraph, class ValueMap, class BinaryFunction, class ValueType>
+  void
+  write_vmap_to_image(const MyGraph& g, const tree_t* t, const tlink_t* tlink,
+                      const ValueMap& vmap, BinaryFunction op, ValueType init,
+                      image2d<ValueType>& out)
+  {
+    (void) g;
+
+    mln_pixter(px, out);
+    mln_forall(px)
+    {
+      ValueType w = init;
+      for (int i = 0; i < NTREE; ++i)
+        {
+          tree_t::node_type tnode = t[i].get_node_at(px->index());
+          typename MyGraph::vertex_descriptor gnode = tlink[i][tnode];
+          w = op(w, vmap[gnode]);
+        }
+      px->val() = w;
+    }
+
+  }
+
+
+  tree_t compute_ctos(const image2d<rgb8>& input,
+                      image2d<uint16>* imdepth)
+  {
+    typedef rgb8 value_t;
+    enum { NTREE = value_t::ndim };
+
+
+    image2d<value_t> f = addborder_marginal(input);
+    image2d<value_t> F = immerse_k1(f);
+
+    point2d pmin = {0,0};
+
+    /// Compute the marginal tree
+    tree_t t[NTREE];
+    tbb::parallel_for(0, (int)NTREE, [&t,&f,pmin](int i){
+        t[i] = morpho::cToS_pinf(imtransform(f, [i](value_t x) {
+              return x[i]; }), c4, pmin);
+      });
+
+    /// Compute the graph
+    Graph<NTREE> g2;
+    std::array<property_map<tree_t, typename MyGraph::vertex_descriptor>, NTREE> tlink;
+    std::tie(g2, tlink) = compute_g2<NTREE>(t);
+
+
+    /// Compute depth
+    boost::vector_property_map<unsigned> gdepth;
+    gdepth = compute_graph_depth(g2);
+
+
+    /// Compute the pw depth image from graph
+    image2d<uint16> maxdepth;
+    resize(maxdepth, F);
+    write_vmap_to_image(g2, t, &tlink[0], gdepth,
+                        functional::max_t<unsigned> (), uint16(0), maxdepth);
+
+    /// Compute the saturated maxtree
+    tree_t tw;
+    property_map<tree_t, uint16> depth2;
+    std::tie(tw, depth2) = satmaxtree(maxdepth, pmin);
+
+    /// Remove the 2F from the tree
+    tw = tree_keep_2F(tw);
+
+    /// Outputs the depth image
+    if (imdepth) {
+      resize(*imdepth, maxdepth);
+      morpho::reconstruction(tw, depth2, *imdepth);
+    }
+    return tw;
+  }
+
+}

apps/g2/compute_ctos.hpp

+#ifndef APPS_G2_COMPUTE_CTOS_HPP
+# define APPS_G2_COMPUTE_CTOS_HPP
+
+# include <mln/core/image/image2d.hpp>
+# include <mln/core/colors.hpp>
+# include <mln/core/grays.hpp>
+# include <mln/morpho/component_tree/component_tree.hpp>
+
+namespace mln
+{
+
+
+  /// \brief A convenient wrapper around:
+  /// * The Graph of shapes computation + depth attribute valuation
+  /// * The Saturated Maxtree algorithm
+  ///
+  /// \param input The input image
+  /// \param[out] depth if not NULL, store the depth image inside.
+  morpho::component_tree<unsigned, image2d<unsigned> >
+  compute_ctos(const image2d<rgb8>& input,
+               image2d<uint16>* depth = nullptr);
+
+
+}
+
+#endif // ! APPS_G2_COMPUTE_CTOS_HPP

apps/g2/g2-maxdepth.cpp

@@ -341,10 +341,17 @@ int main(int argc, char** argv)
   typedef value_t::value_type V;
 
   /// Compute the marginal ToS
+  bool ASK = false;
+  point2d pmin{0,0};
+  if (ASK) {
+    std::cout << "Point a l'infini: " << std::endl;
+    std::cin >> pmin[0] >> pmin[1];
+  }
+
   tree_t t[NTREE];
 
-  tbb::parallel_for(0, (int)NTREE, [&t,&f](int i){
-      t[i] = morpho::cToS(imtransform(f, [i](value_t x) { return x[i]; }), c4);
+  tbb::parallel_for(0, (int)NTREE, [&t,&f,pmin](int i){
+      t[i] = morpho::cToS_pinf(imtransform(f, [i](value_t x) { return x[i]; }), c4, pmin);
     });
 
   /// Compute the graph

apps/g2/satmaxtree-cli.cpp

@@ -10,53 +10,6 @@
 #include <apps/tos/topology.hpp>
 #include "satmaxtree.hpp"
 
-namespace mln
-{
-
-  template <class P>
-  morpho::component_tree<P, image2d<P> >
-  tree_keep_2F(const morpho::component_tree<P, image2d<P> >& tree)
-  {
-    morpho::component_tree<P, image2d<P> > out;
-
-    auto newdata = out._get_data();
-    auto olddata = tree._get_data();
-
-    // 1. Copy the point2node map
-    box2d olddom = olddata->m_pmap.domain();
-    box2d dom;
-    dom.pmin = olddom.pmin / 2;
-    dom.pmax = (olddom.pmax + 1) / 2;
-    newdata->m_pmap.resize(dom);
-    copy(olddata->m_pmap | sbox2d{olddom.pmin, olddom.pmax, {2,2}},
-         newdata->m_pmap);
-
-    // 2. Copy the node
-    newdata->m_nodes = olddata->m_nodes;
-
-    // 3. Copy the point set and update node first point index/
-    newdata->m_S.resize(dom.size());
-    unsigned j = 0;
-    for (unsigned i = 0; i < olddata->m_S.size(); ++i)
-      {
-        P p = olddata->m_S[i];
-        point2d pt = olddata->m_pmap.point_at_index(p);
-        if (K1::is_face_2(pt))
-          {
-            newdata->m_S[j] = newdata->m_pmap.index_of_point(pt/2);
-            auto node = tree.get_node_at(p);
-            if (node.get_first_point_id() == i)
-              newdata->m_nodes[node.id()].m_point_index = j;
-            ++j;
-          }
-      }
-    // 4. Do not forget the sentinel
-    newdata->m_nodes[out.npos()].m_point_index = j;
-
-    return out.get_subtree(tree.get_root_id());
-  }
-
-}
 
 void usage(char** argv)
 {

apps/g2/satmaxtree.cpp

@@ -3,6 +3,7 @@
 #include <mln/morpho/component_tree/accumulate.hpp>
 #include <mln/morpho/component_tree/filtering.hpp>
 #include <apps/tos/Kinterpolate.hpp>
+#include <apps/tos/topology.hpp>
 #include "satmaxtree.hpp"
 
 namespace mln
@@ -12,13 +13,21 @@ namespace mln
     morpho::component_tree<unsigned, image2d<unsigned> >,
     property_map<morpho::component_tree<unsigned, image2d<unsigned> >, uint16>
     >
-  satmaxtree(const image2d<uint16>& f)
+  satmaxtree(const image2d<uint16>& f, point2d pmin)
   {
     mln_entering("mln::satmaxtree");
 
     typedef morpho::component_tree<unsigned, image2d<unsigned> > T;
 
-    T tree = morpho::cToS(f, c4);
+    // bool ASK = true;
+    // point2d pmin{0,0};
+    // if (ASK) {
+    //   std::cout << "Point a l'infini: " << std::endl;
+    //   std::cin >> pmin[0] >> pmin[1];
+    //   pmin *= 2;
+    // }
+    T tree = morpho::cToS_pinf(f, c4, pmin);
+
     image2d<uint16> F = immerse_k1(f, 69);
     property_map<T, uint16> vmap = morpho::make_attribute_map_from_image(tree, F);
     vmap[tree.npos()] = 0;
@@ -60,4 +69,51 @@ namespace mln
     return {std::move(tree), std::move(vmap)};
   }
 
+
+
+
+  morpho::component_tree<unsigned, image2d<unsigned> >
+  tree_keep_2F(const morpho::component_tree<unsigned, image2d<unsigned> >& tree)
+  {
+    typedef unsigned P;
+    morpho::component_tree<P, image2d<P> > out;
+
+    auto newdata = out._get_data();
+    auto olddata = tree._get_data();
+
+    // 1. Copy the point2node map
+    box2d olddom = olddata->m_pmap.domain();
+    box2d dom;
+    dom.pmin = olddom.pmin / 2;
+    dom.pmax = (olddom.pmax + 1) / 2;
+    newdata->m_pmap.resize(dom);
+    copy(olddata->m_pmap | sbox2d{olddom.pmin, olddom.pmax, {2,2}},
+         newdata->m_pmap);
+
+    // 2. Copy the node
+    newdata->m_nodes = olddata->m_nodes;
+
+    // 3. Copy the point set and update node first point index/
+    newdata->m_S.resize(dom.size());
+    unsigned j = 0;
+    for (unsigned i = 0; i < olddata->m_S.size(); ++i)
+      {
+        P p = olddata->m_S[i];
+        point2d pt = olddata->m_pmap.point_at_index(p);
+        if (K1::is_face_2(pt))
+          {
+            newdata->m_S[j] = newdata->m_pmap.index_of_point(pt/2);
+            auto node = tree.get_node_at(p);
+            if (node.get_first_point_id() == i)
+              newdata->m_nodes[node.id()].m_point_index = j;
+            ++j;
+          }
+      }
+    // 4. Do not forget the sentinel
+    newdata->m_nodes[out.npos()].m_point_index = j;
+
+    return out.get_subtree(tree.get_root_id());
+  }
+
+
 }

apps/g2/satmaxtree.hpp

@@ -7,11 +7,19 @@
 namespace mln
 {
 
+  /// \brief Compute the saturated maxtree of an image
   std::pair<
     morpho::component_tree<unsigned, image2d<unsigned> >,
     property_map<morpho::component_tree<unsigned, image2d<unsigned> >, uint16>
     >
-  satmaxtree(const image2d<uint16>& f);
+  satmaxtree(const image2d<uint16>& f, point2d pinf = {0,0});
+
+
+
+  /// \brief Remove the 2F in the tree.
+  morpho::component_tree<unsigned, image2d<unsigned> >
+  tree_keep_2F(const morpho::component_tree<unsigned, image2d<unsigned>>& tree);
+
 }
 
 

mln/morpho/tos/ctos.hpp

@@ -56,6 +56,10 @@ namespace mln
     morpho::component_tree<typename I::size_type, mln_ch_value(I, unsigned)>
     cToS(const Image<I>& ima, const Neighborhood& nbh, const Compare& cmp = Compare () );
 
+    template <typename I, typename Neighborhood, typename Compare = std::less<mln_value(I)> >
+    morpho::component_tree<typename I::size_type, mln_ch_value(I, unsigned)>
+    cToS_pinf(const Image<I>& ima, const Neighborhood& nbh, mln_point(I) pmin,
+              const Compare& cmp = Compare ());
 
     /***********************************************/
     /* Same as before but using priority proagation */
@@ -167,7 +171,7 @@ namespace mln
         mln_entering("mln::morpho::tos - propagation");
 
         pset_t W(K, cmp);
-        size_type p = f.index_of_point(pmin);
+        size_type p = f.index_of_point(pmin * 2);
         W.insert(p);
         pmap[p] = PROCESSED;
         K[p] = f[p].lower;
@@ -329,6 +333,14 @@ namespace mln
       return cToS(ima, nbh, pmin, cmp, internal::equiv<Compare> (cmp));
     }
 
+    template <typename I, typename Neighborhood, typename Compare>
+    morpho::component_tree<typename I::size_type, mln_ch_value(I, unsigned)>
+    cToS_pinf(const Image<I>& ima, const Neighborhood& nbh, mln_point(I) pmin,
+              const Compare& cmp)
+    {
+      return cToS(ima, nbh, pmin, cmp, internal::equiv<Compare> (cmp));
+    }
+
 
     template <typename I,
               typename Neighborhood,