Implement maxtree tests.

fixtures/ImageCompare/src/image_compare.cpp

@@ -37,19 +37,13 @@ namespace fixtures::ImageCompare::experimental::impl
 
 
     mln::sample_type_id id = a.sample_type();
-
-    // If the comparison function is not provided, try memcmp
+    // If the comparison function is not provided, we use memcmp (type are supposed to be trivially comparable
     if (!linecmp_fn)
-    {
-      if (id == mln::sample_type_id::OTHER)
-        return ::testing::AssertionFailure() << "A and B value type are not trivially comparable."; // Cannot be memcmp
       linecmp_fn = std::memcmp;
-    }
-
-    std::size_t n = a.width() * mln::get_sample_type_id_traits(id).size();
 
-    auto a_buf = A.buffer();
-    auto b_buf = B.buffer();
+    std::size_t n     = a.width() * mln::get_sample_type_id_traits(id).size();
+    auto        a_buf = A.buffer();
+    auto        b_buf = B.buffer();
 
     for (int w = 0; w < A.size(3); ++w)
       for (int z = 0; z < A.size(2); ++z)

pylene/CMakeLists.txt

@@ -61,6 +61,8 @@ target_sources(Pylene PRIVATE
                src/io/io.cpp
                src/io/imread.cpp
                src/io/imprint.cpp
+               src/morpho/maxtree.cpp
+               src/morpho/component_tree.cpp
                )
 
 # Compiler configurations

pylene/include/mln/morpho/experimental/canvas/depthfirst.hpp

@@ -6,7 +6,7 @@
 
 
 #include <range/v3/functional.hpp>
-#include <bitset>
+
 
 namespace  mln::morpho::experimental::canvas
 {
@@ -29,7 +29,7 @@ namespace  mln::morpho::experimental::canvas
   //
   // viz.on_done(Level l, Point p): called when a point has been processed
   //
-  // viz.has_node_at_level(Level l) : true if there is a node at this level in branch
+  // viz.has_level_flooding_started(Level l) : true if this level is already being flooded
 
   template <class I, class N, class DFVisitor, class Proj>
   void depthfirst(I& f, N nbh, DFVisitor& viz, image_point_t<I> start)
@@ -57,22 +57,19 @@ namespace  mln::morpho::experimental::canvas
     queue.push(f(start), start);
     status(start) = INQUEUE;
 
-    constexpr int nvalues = 1 << 8;
-    std::bitset<nvalues> has_level;
 
     // Flooding function turned non-recursive: flood(p, current_level = f(p))
     {
       P    p             = start;
       auto current_level = f(start);
 
-    flood_new_level:
-      has_level.set(current_level);
+    flood:
       viz.on_flood_start(current_level, p);
 
-    flood_flat_zone:
+    keep_flooding:
       for (auto n : nbh(p))
       {
-        if (status(n) != NONE)
+        if (status.at(n) != NONE)
           continue;
 
         // Insert n INQUEUE
@@ -87,7 +84,7 @@ namespace  mln::morpho::experimental::canvas
         // Otherwise, process it, (do not remove p from stack)
         current_level = nval;
         p             = n;
-        goto flood_new_level;
+        goto flood;
       }
 
       // All the neighbors have been seen, p is DONE
@@ -96,23 +93,24 @@ namespace  mln::morpho::experimental::canvas
       queue.pop();
 
       // If the queue gets empty, we have processed the whole image
-      if (!queue.empty())
+      if (queue.empty())
+        goto end_flooding;
+
       {
         auto old_level = current_level;
         std::tie(current_level, p) = queue.top();
         if (current_level == old_level)
-          goto flood_flat_zone;
-
-        viz.on_flood_end(old_level);
-        has_level.reset(old_level);
-
-        if (has_level.test(current_level))
-          goto flood_flat_zone;
-
-        goto flood_new_level;
+          goto keep_flooding;
+        viz.on_flood_end(old_level, current_level);
       }
 
-    // End: there is no more point to process
+      if (viz.has_level_flooding_started(current_level))
+        goto keep_flooding;
+      else
+        goto flood;
+
+      // End: there is no more point to process
+    end_flooding:
       viz.on_flood_end(current_level);
     }
   }

pylene/include/mln/morpho/experimental/component_tree.hpp

+#pragma once
+
+
+#include <mln/core/image/view/zip.hpp>
+#include <mln/core/rangev3/foreach.hpp>
+#include <mln/accu/accumulator.hpp>
+
+#include <vector>
+
+#include <range/v3/view/span.hpp>
+
+namespace mln::morpho::experimental
+{
+
+  // namespace details
+  // {
+  //   void filter_treeT(const int* parent, std::size_t n, e_filtering strategy, F predicate);
+  //   void filter_tree(const int* parent, std::size_t n, e_filtering strategy, std::function<bool(int)> predicate);
+  // }
+
+  enum ct_filtering
+  {
+    CT_FILTER_DIRECT,
+    CT_FILTER_MIN,
+    CT_FILTER_MAX,
+    CT_FILTER_SUBTRACTIVE
+  };
+
+
+  template <class T = void>
+  class component_tree;
+
+
+
+  template <>
+  class component_tree<void>
+  {
+  public:
+    using node_id_type = int;
+    using node_map_t = int;
+
+
+
+
+    template <class F>
+    void filter(ct_filtering strategy, F pred);
+
+    template <class I, class F>
+    void filter(ct_filtering strategy, I node_map, F pred);
+
+
+
+    /// \brief Compute the depth attribute over a tree
+    std::vector<int> compute_depth_map() const;
+
+
+    /// \brief Compute attribute on values
+    template <class I, class J, class Accu>
+    void compute_attribute_on_values(I node_map, J input, Accu acc);
+
+    /// \brief Compute attribute on values
+    template <class I, class J, class Accu>
+    void compute_attribute_on_pixels(I node_map, J input, Accu acc);
+
+
+    /// \brief Compute attribute on points
+    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);
+
+
+    /// \brief Reconstruct an image from an attribute map
+    template <class I, class V>
+    image_ch_value_t<I, V> reconstruct_from(I node_map, ::ranges::span<V> values) const;
+
+    template <class I, class F>
+    void update_node_map(I node_map, F pred) const;
+
+    /// \brief Filter an image
+
+    // template <class I, class J, class F>
+    // void filter(e_filtering strategy, I node_map, J input, ::ranges::span<V> values, F pred) const;
+
+
+
+    using node_t = int;
+    std::vector<node_t> parent;
+
+  private:
+    void filter_direct(const std::vector<bool>& pred);
+
+    template <class F>
+    void filter_direct_T(F pred);
+  };
+
+
+  template <class T>
+  class component_tree : public component_tree<void>
+  {
+  public:
+
+    template <class I>
+    image_ch_value_t<I, T> reconstruct(I node_map)
+    {
+      return this->reconstruct_from(std::move(node_map), ::ranges::make_span(values.data(), values.size()));
+    }
+
+
+    std::vector<T> values;
+  };
+
+
+  /******************************************/
+  /****          Implementation          ****/
+  /******************************************/
+
+  template <class F>
+  void component_tree<void>::filter_direct_T(F pred)
+  {
+    std::size_t n = parent.size();
+
+    for (std::size_t i = 1; i < n; ++i)
+      if (!pred(parent[i]))
+        parent[i] = parent[parent[i]];
+  }
+
+
+  template <class I, class F>
+  void component_tree<void>::update_node_map(I node_map, F pred) const
+  {
+    mln_foreach_new(auto& id, node_map.new_values())
+    {
+      if (!pred(id))
+        id = this->parent[id];
+    }
+  }
+
+
+
+  template <class F>
+  void component_tree<void>::filter(ct_filtering strategy, F pred)
+  {
+    std::size_t n = parent.size();
+
+    // Node pass status
+    std::vector<bool> pass;
+
+    if (strategy == CT_FILTER_DIRECT)
+    {
+      this->filter_direct_T(pred);
+    }
+    else if (strategy == CT_FILTER_MIN)
+    {
+      // A node is removed if one ancestror is remove
+      pass.resize(n);
+      pass[0] = true;
+
+
+      // Propagate downward
+      for (std::size_t i = 1; i < n; ++i)
+      {
+        pass[i] = pred(i) && pass[parent[i]];
+        if (!pass[parent[i]])
+          parent[i] = parent[parent[i]];
+      }
+    }
+    else if (strategy == CT_FILTER_MAX)
+    {
+      // A node is removed if all its descandants are removed
+      pass.assign(n, false);
+
+      // Propagate upward
+      for (std::size_t i = n - 1; i > 0; --i)
+      {
+        pass[i] = pass[i] || pred(i);
+        pass[parent[i]] = pass[parent[i]] || pass[i];
+      }
+
+      this->filter_direct(pass);
+    }
+  }
+
+  template <class I, class F>
+  void component_tree<void>::filter(ct_filtering strategy, I node_map, F pred)
+  {
+    std::size_t n = parent.size();
+
+    // Node pass status
+    std::vector<bool> pass;
+
+    if (strategy == CT_FILTER_DIRECT)
+    {
+      this->filter_direct_T(pred);
+      this->update_node_map(node_map, pred);
+      return;
+    }
+    else if (strategy == CT_FILTER_MIN)
+    {
+      // A node is removed if one ancestror is removed
+      pass.resize(n);
+      pass[0] = true;
+
+
+      // Propagate downward
+      for (std::size_t i = 1; i < n; ++i)
+      {
+        pass[i] = pred(i) && pass[parent[i]];
+        if (!pass[parent[i]])
+          parent[i] = parent[parent[i]];
+      }
+    }
+    else if (strategy == CT_FILTER_MAX)
+    {
+      // A node is removed if all its descandants are removed
+      pass.assign(n, false);
+
+      // Propagate upward
+      for (std::size_t i = n - 1; i > 0; --i)
+      {
+        pass[i] = pass[i] || pred(i);
+        pass[parent[i]] = pass[parent[i]] || pass[i];
+      }
+      this->filter_direct(pass);
+    }
+
+    this->update_node_map(node_map, [&pass](int x) { return pass[x]; });
+  }
+
+
+
+
+
+  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)
+  {
+    static_assert(mln::is_a<I, mln::experimental::Image>());
+    using R = typename accu::result_of<Accu, image_point_t<I>>::type;
+
+    auto a = accu::make_accumulator(acc, image_point_t<I>());
+
+    std::vector<decltype(a)> attr(parent.size(), a);
+
+    // Accumulate for each point
+    mln_foreach_new(auto px, node_map.new_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]);
+
+
+    // Extract values
+    std::vector<R> out(n);
+    for (std::size_t i = 0; i < n; ++i)
+      out[i] = attr[i].to_result();
+
+    return out;
+  }
+
+
+
+  template <class I, class V>
+  image_ch_value_t<I, V> component_tree<void>::reconstruct_from(I node_map, ::ranges::span<V> values) const
+  {
+    image_ch_value_t<I, V> out = imchvalue<V>(node_map);
+
+    auto zz = mln::view::zip(node_map, out);
+    mln_foreach_new ((auto&& [node_id, val]), zz.new_values())
+      val = values[node_id];
+
+    return out;
+  }
+
+
+} // namespace mln::morpho::experimental

pylene/include/mln/morpho/experimental/maxtree.hpp

+#pragma once
+
+#include <mln/core/concept/new/neighborhoods.hpp>
+#include <mln/core/concept/new/images.hpp>
+
+#include <mln/morpho/experimental/canvas/depthfirst.hpp>
+#include <mln/morpho/experimental/component_tree.hpp>
+#include <mln/core/dontcare.hpp>
+#include <mln/core/algorithm/for_each.hpp>
+
+#include <vector>
+
+namespace mln::morpho::experimental
+{
+
+  template <class I, class N>
+  std::pair<component_tree<image_value_t<I>>, image_ch_value_t<I, component_tree<>::node_id_type>> //
+  maxtree(I input, N nbh);
+
+
+  /******************************************/
+  /****          Implementation          ****/
+  /******************************************/
+
+
+  namespace details
+  {
+    template <class I>
+    struct maxtree_visitor
+    {
+      using point_t  = image_point_t<I>;
+      using level_t = int;
+      using node_id_t = component_tree<>::node_id_type;
+
+      constexpr void on_flood_start(level_t level, dontcare_t) noexcept
+      {
+        nodes.push_back({-1});
+        lvlnodes.push_back(level);
+        roots.push_back({static_cast<node_id_t>(nodes.size() - 1), level});
+      }
+
+      void on_flood_end(level_t old_level, level_t new_level) noexcept
+      {
+        // Attach to parent
+        auto current = roots.back();
+        roots.pop_back();
+
+        assert(old_level == current.level);
+        assert(old_level > new_level);
+
+
+        auto par = roots.back();
+        assert(par.level <= new_level);
+
+        if (par.level != new_level)
+        {
+          nodes.push_back({-1});
+          lvlnodes.push_back(new_level);
+          par = {static_cast<node_id_t>(nodes.size() - 1), new_level};
+          roots.push_back(par);
+        }
+
+        nodes[current.node_root_id] = par.node_root_id;
+      }
+
+      void on_flood_end(dontcare_t) noexcept
+      {
+        // Attach to parent
+        node_id_t root = roots.back().node_root_id;
+        roots.pop_back();
+
+        assert(roots.empty());
+        nodes[root] = -1;
+      }
+
+      void on_done(level_t, point_t p) { node_map(p) = roots.back().node_root_id; }
+
+
+      constexpr bool has_level_flooding_started(dontcare_t) const noexcept { return true; }
+
+      struct root_info
+      {
+        node_id_t node_root_id;
+        level_t   level;
+      };
+
+
+
+
+      std::vector<root_info>                        roots;
+      std::vector<int>                              nodes;
+      std::vector<image_value_t<I>>                 lvlnodes;
+      image_ch_value_t<I, node_id_t>                node_map;
+    };
+  } // namespace details
+
+  namespace details
+  {
+
+    /// Permute the array \p tab considering the permutation array \p permut
+    ///
+    /// \param count Number of elements
+    /// \param size Size of the elements (in bytes)
+    void permute_array(int* permut, void* tab, std::size_t count, std::size_t size);
+
+
+    /// Sort the array \c par so that par[i] < i
+    /// where `par` encodes a DAG relation
+    /// The parent of the root node is supposed to be -1;
+    void permute_parent(int* par, int* permut, std::size_t n);
+
+
+    template <class I, class V>
+    [[gnu::noinline]] //
+    void
+    permute_parent_and_node_map(I&& node_map, int* parent, V* levels, std::size_t n)
+    {
+      auto permutation_arr = std::make_unique<int[]>(n);
+
+      permute_parent(parent, permutation_arr.get(), n);
+      permute_array(permutation_arr.get(), levels, n, sizeof(V));
+      mln::for_each(node_map, [perm = permutation_arr.get()](int& i) { i = perm[i]; });
+    }
+  } // namespace details
+
+  template <class I, class N>
+  std::pair<component_tree<image_value_t<I>>, image_ch_value_t<I, component_tree<>::node_id_type>> //
+  maxtree(I input, N nbh)
+  {
+    using V = image_value_t<I>;
+    constexpr std::size_t kStackReserve = 256;
+
+    details::maxtree_visitor<I> viz;
+    viz.roots.reserve(kStackReserve);
+    mln::resize(viz.node_map, input);
+
+    // FIXME: implement ::ranges::first() that return the first element of the domain
+    image_point_t<I> pstart = input.domain().tl();
+    canvas::depthfirst(input, nbh, viz, pstart);
+
+    // Reoder the parent array
+    details::permute_parent_and_node_map(viz.node_map, viz.nodes.data(), viz.lvlnodes.data(), viz.nodes.size());
+
+    component_tree<V> ct;
+    ct.parent = std::move(viz.nodes);
+    ct.values = std::move(viz.lvlnodes);
+
+    return { std::move(ct), std::move(viz.node_map) };
+  }
+
+
+
+} // namespace mln::morpho::experimental

pylene/include/mln/morpho/experimental/private/hpqueue.hpp

@@ -48,7 +48,7 @@ namespace  mln::morpho::experimental::detail
         m_current_level = std::max(m_current_level, level);
     }
 
-
+    bool              has_key(int level) const noexcept { return !base::empty(level); }
     bool              empty() const noexcept { return base::empty(m_current_level); }
 
     std::pair<int, P> top() const noexcept { return {m_current_level, base::front(m_current_level)}; }

pylene/include/mln/morpho/experimental/private/pqueue.hpp

@@ -24,6 +24,7 @@ namespace mln::morpho::experimental::details
     void                            pop();
     std::pair<key_type, value_type> top() const;
     bool                            empty() const;
+    bool                            has_key(const key_type& priority) const;
 
   private:
     static_assert(!std::is_signed_v<key_type>, "Must not be signed");
@@ -78,5 +79,12 @@ namespace mln::morpho::experimental::details
     return m_delegate.top();
   }
 
+  template <class I, bool reverse>
+  inline bool pqueue_fifo<I, reverse>::has_key(const key_type& k) const
+  {
+    return m_delegate.has_key(k);
+  }
+
+
 
 } // namespace mln::morpho::experimental::details

pylene/src/morpho/component_tree.cpp

+#include <mln/morpho/experimental/component_tree.hpp>
+
+
+namespace mln::morpho::experimental
+{