Merge branch 'development/fix-canonization' into 'next'

Fix canonization of the alphatree and watershed hierarchy

See merge request !124

doc/source/conf.py

@@ -99,13 +99,8 @@ pygments_style = 'sphinx'
 # a list of builtin themes.
 html_theme = 'sphinx_rtd_theme'
 
-
-html_context = {
-    'css_files': [
-        '_static/code.css',
-        '_static/theme_overrides.css',  # override wide tables in RTD theme
-        ],
-}
+html_static_path = ['_static']
+html_css_files = ['code.css', 'theme_overrides.css']
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the

doc/source/morpho/component_tree.rst

@@ -148,15 +148,16 @@ A :cpp:class:`component_tree` `t` has the following methods.
 
 .. cpp:namespace-push::  template <class V> component_tree
 
-.. cpp:function:: auto compute_attribute_on_points(Image node_map, Accumulator accu) const
-                  auto compute_attribute_on_values(Image node_map, Image values, Accumulator accu) const
-                  auto compute_attribute_on_pixels(Image node_map, Image values, Accumulator accu) const
+.. cpp:function:: auto compute_attribute_on_points(Image node_map, Accumulator accu, bool propagate) const
+                  auto compute_attribute_on_values(Image node_map, Image values, Accumulator accu, bool propagate) const
+                  auto compute_attribute_on_pixels(Image node_map, Image values, Accumulator accu, bool propagate) const
 
     Accumulate the points of each component.
 
     :param node_map: An image thats maps: point -> node id
     :param values: An image where values to accumlate are taken from
     :param accu: The feature to compute
+    :param propagate: Option to propagate the values to the parent (default: true)
 
     :return: A vector mapping for each node the result of accumulation.
 

doc/source/morpho/multivariate_component_tree.rst

@@ -44,26 +44,6 @@ This example computes a grain filter, which removes all the node having an area
         // Function to get the median of the border values
         mln::rgb8 get_median_on_border(mln::image2d<mln::rgb8> ima);
 
-        // Accumulator to compute the mean of the pixel values of each node, without taking into account the values of the holes
-        struct mean_node_accu : mln::Accumulator<mean_node_accu>
-        {
-          using result_type = decltype(std::declval<mln::rgb8>() + std::declval<mln::rgb8>());
-        public:
-          void take(const mln::rgb8& v)
-          {
-            m_sum += v;
-            m_count++;
-          }
-
-          void take(const mean_node_accu&) {}
-
-          result_type to_result() const { return m_count > 1 ? static_cast<result_type>(m_sum / m_count) : m_sum; }
-
-        private:
-          result_type m_sum = {0, 0, 0};
-          int m_count = 0;
-        };
-
         int main(void)
         {
             mln::image2d<mln::rgb8> ima;
@@ -92,7 +72,7 @@ This example computes a grain filter, which removes all the node having an area
             t.filter(mln::morpho::CT_FILTER_DIRECT, nm, [&area](int n) { return area[n] >= 100; });
 
             // Compute the mean of the connected component of each nodes
-            auto mean = t.compute_attribute_on_values(nm, ima, mean_node_accu());
+            auto mean = t.compute_attribute_on_values(nm, ima, mln::accu::features::mean<>(), false);
 
             // Reconstruct the tree
             auto rec = t.reconstruct_from(nm, ranges::make_span(mean.data(), mean.size()));

doc/source/snippets/mtos_example.cpp

 #include <mln/accu/accumulators/count.hpp>
 #include <mln/accu/accumulators/max.hpp>
+#include <mln/accu/accumulators/mean.hpp>
 #include <mln/core/algorithm/accumulate.hpp>
 #include <mln/core/colors.hpp>
 #include <mln/core/extension/padding.hpp>
@@ -20,26 +21,6 @@
 
 namespace
 {
-  struct mean_node_accu : mln::Accumulator<mean_node_accu>
-  {
-    using result_type = decltype(std::declval<mln::rgb8>() + std::declval<mln::rgb8>());
-
-  public:
-    void take(const mln::rgb8& v)
-    {
-      m_sum += v;
-      m_count++;
-    }
-
-    void take(const mean_node_accu&) {}
-
-    result_type to_result() const { return m_count > 1 ? static_cast<result_type>(m_sum / m_count) : m_sum; }
-
-  private:
-    result_type m_sum   = {0, 0, 0};
-    int         m_count = 0;
-  };
-
   mln::rgb8 get_median_on_border(mln::image2d<mln::rgb8> ima)
   {
     std::vector<mln::rgb8> border;
@@ -134,7 +115,7 @@ int main(int argc, char* argv[])
 
   auto area = t.compute_attribute_on_points(nm, mln::accu::accumulators::count<int>());
   t.filter(mln::morpho::CT_FILTER_DIRECT, nm, [&area](int n) { return area[n] >= 100; });
-  auto mean = t.compute_attribute_on_values(nm, ima, mean_node_accu());
+  auto mean = t.compute_attribute_on_values(nm, ima, mln::accu::features::mean<>(), false);
   auto rec  = t.reconstruct_from(nm, ranges::make_span(mean.data(), mean.size()));
   mln::io::imsave(mln::view::cast<mln::rgb8>(rec), argv[3]);
 

pylene/include/mln/morpho/alphatree.hpp

@@ -16,7 +16,6 @@
 
 #include <range/v3/functional/concepts.hpp>
 
-
 namespace mln::morpho
 {
 
@@ -237,14 +236,12 @@ namespace mln::morpho
       return node_count;
     }
 
-    template <class W>
+    template <class W, class I>
     std::pair<std::vector<int>, std::vector<W>> canonize_component_tree(const std::vector<int>& par,        //
                                                                         const std::vector<W>&   levels,     //
-                                                                        std::size_t             node_count, //
-                                                                        std::size_t             nb_leaves)
+                                                                        I node_map
+                                                                        )
     {
-      assert(node_count >= nb_leaves);
-
       std::vector<int> canonized_par;
       std::vector<W>   canonized_levels;
 
@@ -253,16 +250,13 @@ namespace mln::morpho
       canonized_par.push_back(0);
       canonized_levels.push_back(levels[0]);
 
-      std::vector<int> translation_map(node_count);
+      std::vector<int> translation_map(par.size());
 
       translation_map[0] = 0;
       int count          = 1;
 
-      std::size_t begin_leaves = node_count - nb_leaves;
-
       // Build canonized component tree
-
-      for (std::size_t i = 1; i < begin_leaves; ++i)
+      for (std::size_t i = 1; i < par.size(); ++i)
       {
         if (levels[i] != levels[par[i]]) // Keep the node: Update tree
         {
@@ -274,12 +268,7 @@ namespace mln::morpho
           translation_map[i] = translation_map[par[i]];
       }
 
-      for (std::size_t i = begin_leaves; i < node_count; ++i)
-      {
-        translation_map[i] = count++;
-        canonized_par.push_back(translation_map[par[i]]);
-        canonized_levels.push_back(levels[i]);
-      }
+      mln::for_each(node_map, [&translation_map, node_count=par.size()](int& id) { id = translation_map[static_cast<int>(node_count) - id - 1]; });
 
       return {canonized_par, canonized_levels};
     }
@@ -298,20 +287,19 @@ namespace mln::morpho
       // Parent / levels are ordered from leaves to root, we need to reverse
       internal::alphatree_reorder_nodes(par.data(), levels.data(), node_count);
 
+      component_tree<W> t;
       if (canonize_tree)
       {
-        auto [canonized_par, canonized_levels] = internal::canonize_component_tree(par, levels, node_count, nb_leaves);
-        par                                    = canonized_par;
-        levels                                 = canonized_levels;
+        auto [canonized_par, canonized_levels] = internal::canonize_component_tree(par, levels, node_map);
+        t.parent = std::move(canonized_par);
+        t.values = std::move(canonized_levels);
+      }
+      else
+      {
+        t.parent = std::move(par);
+        t.values = std::move(levels);
+        mln::for_each(node_map, [node_count](int& id) { id = static_cast<int>(node_count) - id - 1; });
       }
-
-      component_tree<W> t;
-      t.parent   = par;
-      t.values   = levels;
-      node_count = t.parent.size();
-
-      // Update the node map according to the component tree representation
-      mln::for_each(node_map, [node_count](int& id) { id = static_cast<int>(node_count) - id - 1; });
 
       return {std::move(t), std::move(node_map)};
     }

pylene/include/mln/morpho/component_tree.hpp

@@ -88,9 +88,10 @@ namespace mln::morpho
     /// \param node_map Image point -> node_id mapping
     /// \param values Image point -> value mapping
     /// \param acc Accumulator to apply on values
+    /// \param propagate Boolean to propagate the values to the parent
     template <class I, class J, class Accu>
     std::vector<typename accu::result_of<Accu, image_value_t<J>>::type> //
-    compute_attribute_on_values(I node_map, J values, Accu acc) const;
+    compute_attribute_on_values(I node_map, J values, Accu acc, bool propagate = true) const;
 
     /// \brief Compute attribute on values
     ///
@@ -99,9 +100,10 @@ namespace mln::morpho
     /// \param node_map Image point -> node_id mapping
     /// \param values Image point -> value mapping
     /// \param acc Accumulator to apply on points
+    /// \param propagate Boolean to propagate the values to the parent
     template <class I, class Accu>
     std::vector<typename accu::result_of<Accu, image_point_t<I>>::type> //
-    compute_attribute_on_points(I node_map, Accu acc) const;
+    compute_attribute_on_points(I node_map, Accu acc, bool propagate = true) const;
 
 
     /// \brief Compute attribute on pixels
@@ -111,9 +113,10 @@ namespace mln::morpho
     /// \param node_map Image point -> node_id mapping
     /// \param values Image point -> value mapping
     /// \param acc Accumulator to apply on values
+    /// \param propagate Boolean to propagate the values to the parent
     template <class I, class J, class Accu>
     std::vector<typename accu::result_of<Accu, image_pixel_t<J>>::type> //
-    compute_attribute_on_pixels(I node_map, J values, Accu acc) const;
+    compute_attribute_on_pixels(I node_map, J values, Accu acc, bool propagate = true) const;
 
     /// \brief Compute the horizontal cut of a hierarchie at level `threshold` and return a nodemap
     /// valued with the node indices of the lowest nodes satisfying levels[n] > threshold
@@ -297,7 +300,7 @@ namespace mln::morpho
 
   template <class I, class Accu>
   std::vector<typename accu::result_of<Accu, image_point_t<I>>::type>
-  component_tree<void>::compute_attribute_on_points(I node_map, Accu acc) const
+  component_tree<void>::compute_attribute_on_points(I node_map, Accu acc, bool propagate) const
   {
     mln_entering("mln::morpho::component_tree::compute_attribute_on_points");
 
@@ -313,12 +316,13 @@ namespace mln::morpho
     mln_foreach (auto px, node_map.pixels())
       attr[px.val()].take(px.point());
 
-
-    // Propgate to parent
-    std::size_t n = parent.size();
-    for (std::size_t i = n - 1; i > 0; --i)
-      attr[parent[i]].take(attr[i]);
-
+    const std::size_t n = parent.size();
+    if (propagate)
+    {
+      // Propgate to parent
+      for (std::size_t i = n - 1; i > 0; --i)
+        attr[parent[i]].take(attr[i]);
+    }
 
     // Extract values
     std::vector<R> out(n);
@@ -330,7 +334,7 @@ namespace mln::morpho
 
   template <class I, class J, class Accu>
   std::vector<typename accu::result_of<Accu, image_value_t<J>>::type> //
-  component_tree<void>::compute_attribute_on_values(I node_map, J input, Accu acc) const
+  component_tree<void>::compute_attribute_on_values(I node_map, J input, Accu acc, bool propagate) const
   {
     mln_entering("mln::morpho::component_tree::compute_attribute_on_values");
 
@@ -347,12 +351,13 @@ namespace mln::morpho
     mln_foreach ((auto [node_id, val]), zz.values())
       attr[node_id].take(val);
 
-
-    // Propgate to parent
-    std::size_t n = parent.size();
-    for (std::size_t i = n - 1; i > 0; --i)
-      attr[parent[i]].take(attr[i]);
-
+    const std::size_t n = parent.size();
+    if (propagate)
+    {
+      // Propgate to parent
+      for (std::size_t i = n - 1; i > 0; --i)
+        attr[parent[i]].take(attr[i]);
+    }
 
     // Extract values
     std::vector<R> out(n);
@@ -364,7 +369,7 @@ namespace mln::morpho
 
   template <class I, class J, class Accu>
   std::vector<typename accu::result_of<Accu, image_pixel_t<J>>::type> //
-  component_tree<void>::compute_attribute_on_pixels(I node_map, J values, Accu acc) const
+  component_tree<void>::compute_attribute_on_pixels(I node_map, J values, Accu acc, bool propagate) const
   {
     mln_entering("mln::morpho::component_tree::compute_attribute_on_pixels");
 
@@ -380,12 +385,13 @@ namespace mln::morpho
     mln_foreach (auto px, values.pixels())
       attr[node_map(px.point())].take(px);
 
-
-    // Propgate to parent
-    std::size_t n = parent.size();
-    for (std::size_t i = n - 1; i > 0; --i)
-      attr[parent[i]].take(attr[i]);
-
+    const std::size_t n = parent.size();
+    if (propagate)
+    {
+      // Propgate to parent
+      for (std::size_t i = n - 1; i > 0; --i)
+        attr[parent[i]].take(attr[i]);
+    }
 
     // Extract values
     std::vector<R> out(n);

tests/morpho/watershed_hierarchy.cpp

@@ -87,11 +87,8 @@ TEST(Morpho, AreaWatershedHierarchyGray)
       {121, 112, 121, 112, 163}, //
   };
 
-  std::vector<int> expected_parent = {0,  0,  1,  1, 0,  4, 5,  6,  7,  1, 9,  10, 2,  8, 12,
-                                      14, 15, 16, 2, 14, 3, 20, 17, 18, 3, 24, 21, 22, 14};
-
-  std::vector<unsigned long> expected_values = {6, 3, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-                                                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  std::vector<int>           expected_parent = {0, 0, 1, 1, 0, 1, 2, 2, 3, 3};
+  std::vector<unsigned long> expected_values = {6, 3, 2, 2, 0, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
       input,
@@ -121,16 +118,12 @@ TEST(Morpho, DynamicWatershedHierarchyGray)
       {121, 112, 121, 112, 163}, //
   };
 
-  std::vector<int> expected_parent = {0,  0,  1,  0, 3,  4, 5,  6,  1,  8, 9,  2, 7,  11, 13,
-                                      14, 15, 1, 13, 2, 19, 16, 17, 0, 23, 20, 21, 13};
-
-  std::vector<unsigned long> expected_values = {70, 38, 29, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-                                                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  std::vector<int>           expected_parent = {0, 0, 1, 0, 1, 2, 1, 2, 0};
+  std::vector<unsigned long> expected_values = {70, 38, 29, 0, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
-      input,
-      mln::morpho::WatershedAttribute::DYNAMIC,
-      mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); });
+      input, mln::morpho::WatershedAttribute::DYNAMIC, mln::c4,
+      [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); });
 
   ASSERT_EQ(expected_parent.size(), tree.parent.size());
   ASSERT_EQ(expected_values.size(), tree.values.size());
@@ -153,15 +146,12 @@ TEST(Morpho, HeightWatershedHierarchyGray)
       {121, 112, 121, 112, 163}, //
   };
 
-  std::vector<int> expected_parent = {0, 0, 1, 0, 3, 4, 5, 6, 1, 8, 9, 2, 7, 11, 13, 14, 15, 1, 13, 2, 19, 16, 17, 0, 23, 20, 21, 13};
-
-  std::vector<unsigned long> expected_values = {70, 38, 29, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-                                                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  std::vector<int>           expected_parent = {0, 0, 1, 0, 1, 2, 1, 2, 0};
+  std::vector<unsigned long> expected_values = {70, 38, 29, 0, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
-      input,
-      mln::morpho::WatershedAttribute::HEIGHT,
-      mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); });
+      input, mln::morpho::WatershedAttribute::HEIGHT, mln::c4,
+      [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); });
 
   ASSERT_EQ(expected_parent.size(), tree.parent.size());
   ASSERT_EQ(expected_values.size(), tree.values.size());
@@ -184,11 +174,8 @@ TEST(Morpho, AreaWatershedHierarchyGrayHQ)
       {121, 112, 121, 112, 163}, //
   };
 
-  std::vector<int> expected_parent = {0,  0,  1,  1, 0,  4, 5,  6,  7,  1, 9,  10, 2,  8, 12,
-                                      14, 15, 16, 2, 14, 3, 20, 17, 18, 3, 24, 21, 22, 14};
-
-  std::vector<unsigned long> expected_values = {6, 3, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-                                                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  std::vector<int>           expected_parent = {0, 0, 1, 1, 0, 1, 2, 2, 3, 3};
+  std::vector<unsigned long> expected_values = {6, 3, 2, 2, 0, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
       input,
@@ -218,11 +205,8 @@ TEST(Morpho, AreaWatershedHierarchyRGB)
       {{255, 75, 0}, {31, 174, 0}, {255, 75, 0}, {31, 174, 0}, {255, 146, 0}}, //
   };
 
-  std::vector<int> expected_parent = {0, 0,  0,  0,  3, 2, 5,  6, 3, 4, 3,  3,  7,  3,
-                                      9, 12, 12, 16, 3, 1, 19, 2, 5, 1, 23, 20, 21, 26};
-
-  std::vector<unsigned long> expected_values = {3, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-                                                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  std::vector<int>           expected_parent = {0, 0, 0, 0, 2, 1, 2, 1};
+  std::vector<unsigned long> expected_values = {3, 2, 2, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
       input,
@@ -252,11 +236,8 @@ TEST(Morpho, AreaWatershedHierarchyGrayC8)
       {107, 73, 125, 157, 117},  //
   };
 
-  std::vector<int> expected_parent = {0, 0,  0, 0,  1, 3, 2,  3,  1,  4,  6, 5,  7,  8,  9,
-                                      2, 11, 4, 13, 4, 5, 15, 21, 17, 19, 5, 21, 22, 17, 23};
-
-  std::vector<unsigned long> expected_values = {7, 3, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-                                                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  std::vector<int>           expected_parent = {0, 0, 0, 0, 1, 3, 2, 3, 1, 4, 2, 4, 4};
+  std::vector<unsigned long> expected_values = {7, 3, 2, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
       input,
@@ -284,11 +265,8 @@ TEST(Morpho, AreaWatershedHierarchy3DImage)
       {{11, 18, 0}, {25, 17, 11}, {9, 0, 5}}    //
   };
 
-  std::vector<int> expected_parent = {0, 0,  1,  0, 3, 3,  2, 2, 6, 1,  5,  8,  8,  7,  9, 9,
-                                      4, 10, 12, 8, 9, 18, 9, 8, 1, 17, 25, 24, 26, 28, 27};
-
-  std::vector<unsigned long> expected_values = {9, 3, 2, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
-                                                0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+  std::vector<int>           expected_parent = {0, 0, 1, 0, 3, 3, 2, 2, 1, 1};
+  std::vector<unsigned long> expected_values = {9, 3, 2, 2, 0, 0, 0, 0, 0, 0};