Use a function instead of an accumulator for the attribute computation and add...

Use a function instead of an accumulator for the attribute computation and add height and dynamic attribute

doc/source/morpho/watershed_hierarchy.rst

@@ -6,14 +6,14 @@ Hierarchical Watershed
 .. cpp:namespace:: mln::morpho
 
 
-.. cpp:function:: auto watershed_hierarchy(Image f, Accumulator acc, Neighborhood nbh, F dist);
+.. cpp:function:: auto watershed_hierarchy(Image f, AttributeFunction attribute_func, Neighborhood nbh, F dist);
 
     Compute the watershed hierarchy and returns a pair `(tree, node_map)`.
     See :doc:`component_tree` for more information about the representation of tree.
     The implementation is based on [Naj13]_.
 
     :param input: The input image
-    :param acc: The accumulator that define the attribute computation
+    :param attribute_func: The function that define the attribute computation
     :param nbh: The neighborhood relation
     :param dist: The function weighting the edges between two pixels.
     :return: A pair `(tree, node_map)` where *tree* is of type ``component_tree<std::invoke_result_t<F, image_value_t<Image>, image_value_t<Image>>>`` and *node_map* is a mapping between the image pixels and the node of the tree.
@@ -75,8 +75,10 @@ This example is used to generate the grayscale lena watershed hierarchy by area
         mln::image2d<uint8_t> input = ...;
 
         // Compute the watershed hierarchy by area
-        auto area_acc = mln::accu::features::count<>();
-        auto [tree, node_map] = mln::morpho::watershed_hierarchy(input, area_acc, mln::c4);
+        auto area_attribute_func = [](auto tree, auto node_map) -> std::vector<unsigned long> {
+            return tree.compute_attribute_on_points(node_map, mln::accu::features::count<>());
+        };
+        auto [tree, node_map] = mln::morpho::watershed_hierarchy(input, area_attribute_func, mln::c4);
 
         // Compute an attribute (for example the average pixels value at each node, as below)
         auto mean = tree.compute_attribute_on_values(node_map, input, mln::accu::accumulators::mean<uint8_t>());

doc/source/snippets/watershed_hierarchy_example.cpp

@@ -13,7 +13,10 @@ template <typename V>
 void process_example(const mln::image2d<V>& img, const std::string& output_filename, const double threshold)
 {
   // 2. Build the watershed hierarchy
-  auto [t, nm] = mln::morpho::watershed_hierarchy(img, mln::accu::features::count<>(), mln::c4);
+  auto area_attribute_func = [](auto tree, auto nm) -> std::vector<unsigned long> {
+    return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
+  };
+  auto [t, nm] = mln::morpho::watershed_hierarchy(img, area_attribute_func, mln::c4);
 
   // 3. Compute the mean attribute
   auto mean = t.compute_attribute_on_values(nm, img, mln::accu::accumulators::mean<V>());

pylene/include/mln/morpho/alphatree.hpp

@@ -316,7 +316,7 @@ namespace mln::morpho
       return {std::move(t), std::move(node_map)};
     }
 
-    template <class I, class N, class F, bool HQ = true,
+    template <bool HQ = true, class I, class N, class F,
               class M = edge_t<image_point_t<I>, std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>>
     std::pair<component_tree<std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>, image_ch_value_t<I, int>> //
     __alphatree(I input, N nbh, F distance, bool canonize_tree = true, bool compute_flatzones = true,

pylene/include/mln/morpho/watershed_hierarchy.hpp

@@ -5,17 +5,18 @@
 
 namespace mln::morpho
 {
-  // FIXME Pass a std::function instead of an accumulator ?
-
   /// Compute the watershed hierarchy of an image
   ///
   /// \param input The input image
-  /// \param acc The accumulator that define the attribute computation
+  /// \param attribute_func The function that define the attribute computation
   /// \param neighborhood The neighborhood relation
   /// \param distance Distance function
-  template <class I, class Accu, class N, class F = mln::functional::l2dist_t<>>
-  std::pair<component_tree<typename accu::result_of<Accu, image_point_t<I>>::type>, image_ch_value_t<I, int>> //
-  watershed_hierarchy(I input, Accu acc, N nbh, F distance = F{});
+  template <class I, class A, class N, class F = mln::functional::l2dist_t<>>
+  std::pair<component_tree<typename std::invoke_result_t<
+                A, component_tree<std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>,
+                image_ch_value_t<I, int>>::value_type>,
+            image_ch_value_t<I, int>> //
+  watershed_hierarchy(I input, A attribute_func, N nbh, F distance = F{});
 
 
   /******************************************/
@@ -77,14 +78,121 @@ namespace mln::morpho
     }
   } // namespace internal
 
-  template <class I, class Accu, class N, class F>
-  std::pair<component_tree<typename accu::result_of<Accu, image_point_t<I>>::type>, image_ch_value_t<I, int>> //
-  watershed_hierarchy(I input, Accu acc, N nbh, F distance)
+  template <typename W, class I>
+  std::vector<W> height_attribute(component_tree<W> tree, I node_map)
+  {
+    int n         = static_cast<int>(tree.parent.size());
+    int nb_leaves = static_cast<int>(node_map.domain().size());
+
+    std::vector<W> deepest_altitude(n, std::numeric_limits<W>::max());
+    for (int i = n - 1; i >= 0; --i)
+    {
+      int parent = tree.parent[i];
+
+      if (i > (n - nb_leaves - 1))
+        deepest_altitude[i] = tree.values[parent];
+
+      deepest_altitude[parent] = std::min(deepest_altitude[parent], deepest_altitude[i]);
+    }
+
+    std::vector<W> height(n);
+    for (int i = n - 1; i >= 0; --i)
+      height[i] = tree.values[tree.parent[i]] - deepest_altitude[i];
+
+    return height;
+  }
+
+  template <typename W, class I>
+  std::vector<bool> extrema_attribute(component_tree<W> tree, I node_map)
+  {
+    int n         = static_cast<int>(tree.parent.size());
+    int nb_leaves = static_cast<int>(node_map.domain().size());
+
+    std::vector<bool> extrema(n, true);
+    std::fill_n(extrema.end() - nb_leaves, nb_leaves, false);
+
+    for (int i = (n - nb_leaves - 1); i > 0; --i)
+    {
+      int parent = tree.parent[i];
+
+      bool same_weight = tree.values[i] == tree.values[parent];
+      extrema[parent]  = extrema[parent] && same_weight && extrema[i];
+      extrema[i]       = extrema[i] && !same_weight;
+    }
+
+    return extrema;
+  }
+
+  template <typename W, class I>
+  std::vector<W> dynamic_attribute(component_tree<W> tree, I node_map)
+  {
+    int n         = static_cast<int>(tree.parent.size());
+    int nb_leaves = static_cast<int>(node_map.domain().size());
+
+    std::vector<W>   deepest_altitude(n, std::numeric_limits<W>::max());
+    std::vector<int> path_to_minima(n, -1);
+
+    // Compute deepest altitude and path to deepest minima
+    for (int i = (n - nb_leaves - 1); i >= 0; --i)
+    {
+      int parent = tree.parent[i];
+
+      // Deepest non leaf node
+      if (deepest_altitude[i] == std::numeric_limits<W>::max())
+        deepest_altitude[i] = tree.values[i];
+
+      if (deepest_altitude[i] < deepest_altitude[parent])
+      {
+        deepest_altitude[parent] = deepest_altitude[i];
+        path_to_minima[parent]   = i;
+      }
+    }
+
+    auto             extrema = extrema_attribute(tree, node_map);
+    std::vector<int> nearest_minima(n, -1);
+
+    std::vector<W> dynamic(n);
+    dynamic[0] = tree.values[0] - deepest_altitude[0];
+
+    for (int i = 1; i < n; ++i)
+    {
+      int parent = tree.parent[i];
+
+      if (i > (n - nb_leaves - 1))
+      {
+        if (nearest_minima[parent] != -1)
+          dynamic[i] = dynamic[nearest_minima[parent]];
+        else
+          dynamic[i] = 0;
+
+        continue;
+      }
+
+      if (extrema[i])
+        nearest_minima[i] = i;
+      else
+        nearest_minima[i] = nearest_minima[parent];
+
+      if (i == path_to_minima[parent])
+        dynamic[i] = dynamic[parent];
+      else
+        dynamic[i] = tree.values[parent] - deepest_altitude[i];
+    }
+
+    return dynamic;
+  }
+
+  template <class I, class A, class N, class F>
+  std::pair<component_tree<typename std::invoke_result_t<
+                A, component_tree<std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>,
+                image_ch_value_t<I, int>>::value_type>,
+            image_ch_value_t<I, int>> //
+  watershed_hierarchy(I input, A attribute_func, N nbh, F distance)
   {
     std::vector<internal::edge_t<image_point_t<I>, std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>> mst;
-    auto [tree, nm] = internal::__alphatree<I, N, F, false>(input, nbh, distance, false, false, &mst);
+    auto [tree, nm] = internal::__alphatree<false>(input, nbh, distance, false, false, &mst);
 
-    auto attribute = tree.compute_attribute_on_points(nm, acc);
+    auto attribute = attribute_func(tree, nm);
 
     auto node_count = tree.parent.size();
     mln::for_each(nm, [node_count](int& id) { id = static_cast<int>(node_count) - id - 1; });

tests/morpho/watershed_hierarchy.cpp

@@ -8,6 +8,75 @@
 
 #include <gtest/gtest.h>
 
+TEST(Morpho, HeightAttribute)
+{
+  mln::image2d<uint8_t> input = {
+      {163, 112, 42, 121, 112},  //
+      {42, 121, 1, 42, 255},     //
+      {1, 112, 121, 112, 121},   //
+      {112, 255, 42, 1, 42},     //
+      {121, 112, 121, 112, 163}, //
+  };
+
+  std::vector<float> expected_attribute = {
+      134.f, 134.f, 125.f, 70.f, 70.f, 70.f, 70.f, 70.f, 61.f, 38.f, 0.f, 29.f, 38.f, 0.f, 42.f, 0.f, 0.f,
+      0.f,   42.f,  0.f,   0.f,  0.f,  70.f, 0.f,  0.f,  0.f,  0.f,  0.f, 0.f,  0.f,  0.f, 0.f,  0.f, 0.f,
+      0.f,   0.f,   0.f,   0.f,  0.f,  0.f,  0.f,  0.f,  0.f,  0.f,  0.f, 0.f,  0.f,  0.f, 0.f};
+
+  // Build BPT
+  auto [tree, nm] = mln::morpho::internal::__alphatree<false>(
+      input, mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); }, false,
+      false);
+
+  auto attribute = mln::morpho::height_attribute(tree, nm);
+
+  ASSERT_EQ(expected_attribute.size(), attribute.size());
+
+  for (std::size_t i = 0; i < expected_attribute.size(); ++i)
+    EXPECT_EQ(expected_attribute[i], attribute[i]);
+}
+
+TEST(Morpho, DynamicAttribute)
+{
+  mln::image2d<uint8_t> input = {
+      {163, 112, 42, 121, 112},  //
+      {42, 121, 1, 42, 255},     //
+      {1, 112, 121, 112, 121},   //
+      {112, 255, 42, 1, 42},     //
+      {121, 112, 121, 112, 163}, //
+  };
+
+  // Build BPT
+  auto [tree, nm] = mln::morpho::internal::__alphatree<false>(
+      input, mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); }, false,
+      false);
+
+  std::vector<bool> expected_extrema_attribute = {false, false, false, false, false, false, false, false, false, true,
+                                                  false, true,  true,  false, true,  false, false, false, true,  false,
+                                                  false, false, true,  false, false, false, false, false, false, false,
+                                                  false, false, false, false, false, false, false, false, false, false,
+                                                  false, false, false, false, false, false, false, false, false};
+
+  auto extrema_attribute = mln::morpho::extrema_attribute(tree, nm);
+
+  ASSERT_EQ(expected_extrema_attribute.size(), extrema_attribute.size());
+
+  for (std::size_t i = 0; i < expected_extrema_attribute.size(); ++i)
+    EXPECT_EQ(expected_extrema_attribute[i], extrema_attribute[i]);
+
+  std::vector<float> expected_dynamic_attribute = {
+      134.f, 134.f, 134.f, 70.f, 70.f, 70.f, 70.f, 134.f, 70.f,  38.f,  38.f,  29.f,  38.f, 29.f, 134.f, 134.f, 134.f,
+      134.f, 70.f,  70.f,  70.f, 70.f, 70.f, 70.f, 0,     134.f, 134.f, 134.f, 134.f, 38.f, 38.f, 38.f,  0,     134.f,
+      70.f,  70.f,  70.f,  70.f, 38.f, 0,    29.f, 29.f,  70.f,  38.f,  70.f,  70.f,  29.f, 70.f, 0};
+
+  auto dynamic_attribute = mln::morpho::dynamic_attribute(tree, nm);
+
+  ASSERT_EQ(expected_dynamic_attribute.size(), dynamic_attribute.size());
+
+  for (std::size_t i = 0; i < expected_dynamic_attribute.size(); ++i)
+    EXPECT_EQ(expected_dynamic_attribute[i], dynamic_attribute[i]);
+}
+
 TEST(Morpho, AreaWatershedHierarchyGray)
 {
   mln::image2d<uint8_t> input = {
@@ -25,8 +94,11 @@ TEST(Morpho, AreaWatershedHierarchyGray)
                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
-      input, mln::accu::features::count<>(), mln::c4,
-      [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); });
+      input,
+      [](auto tree, auto nm) -> std::vector<unsigned long> {
+        return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
+      },
+      mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); });
 
   ASSERT_EQ(expected_parent.size(), tree.parent.size());
   ASSERT_EQ(expected_values.size(), tree.values.size());
@@ -56,8 +128,11 @@ TEST(Morpho, AreaWatershedHierarchyGrayHQ)
                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
-      input, mln::accu::features::count<>(), mln::c4,
-      [](const auto& a, const auto& b) -> std::uint8_t { return mln::functional::l2dist_t<>()(a, b); });
+      input,
+      [](auto tree, auto nm) -> std::vector<unsigned long> {
+        return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
+      },
+      mln::c4, [](const auto& a, const auto& b) -> std::uint8_t { return mln::functional::l2dist_t<>()(a, b); });
 
   ASSERT_EQ(expected_parent.size(), tree.parent.size());
   ASSERT_EQ(expected_values.size(), tree.values.size());
@@ -87,8 +162,11 @@ TEST(Morpho, AreaWatershedHierarchyRGB)
                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
-      input, mln::accu::features::count<>(), mln::c4,
-      [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); });
+      input,
+      [](auto tree, auto nm) -> std::vector<unsigned long> {
+        return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
+      },
+      mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); });
 
   ASSERT_EQ(expected_parent.size(), tree.parent.size());
   ASSERT_EQ(expected_values.size(), tree.values.size());
@@ -118,8 +196,11 @@ TEST(Morpho, AreaWatershedHierarchyGrayC8)
                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
-      input, mln::accu::features::count<>(), mln::c8,
-      [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); });
+      input,
+      [](auto tree, auto nm) -> std::vector<unsigned long> {
+        return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
+      },
+      mln::c8, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); });
 
   ASSERT_EQ(expected_parent.size(), tree.parent.size());
   ASSERT_EQ(expected_values.size(), tree.values.size());
@@ -147,8 +228,11 @@ TEST(Morpho, AreaWatershedHierarchy3DImage)
                                                 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 
   auto [tree, _] = mln::morpho::watershed_hierarchy(
-      input, mln::accu::features::count<>(), mln::c26,
-      [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); });
+      input,
+      [](auto tree, auto nm) -> std::vector<unsigned long> {
+        return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
+      },
+      mln::c26, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); });
 
   ASSERT_EQ(expected_parent.size(), tree.parent.size());
   ASSERT_EQ(expected_values.size(), tree.values.size());