First impl of the canvas.

pylene/CMakeLists.txt

@@ -4,7 +4,7 @@ find_package(Boost 1.58 REQUIRED)
 find_package(FreeImage REQUIRED)
 find_package(TBB)
 
-find_package(range-v3 0.9.1 REQUIRED)
+find_package(range-v3 0.9.1 REQUIRED CONFIG)
 find_package(fmt 6.0 REQUIRED)
 
 set(PYLENE_USE_TBB YES CACHE BOOL "Set to NO to disable use of TBB and parallelization")

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

+#pragma once
+
+#include <mln/core/concept/new/images.hpp>
+#include <mln/core/extension/border_management.hpp>
+#include <mln/morpho/experimental/private/pqueue.hpp>
+
+
+#include <range/v3/functional.hpp>
+#include <bitset>
+
+namespace  mln::morpho::experimental::canvas
+{
+
+
+
+  template <class I, class N, class DFVisitor, class Proj = ::ranges::cpp20::identity>
+  void depthfirst(I& f, N nbh, DFVisitor& viz, image_point_t<I> start);
+
+  /******************************************/
+  /****          Implementation          ****/
+  /******************************************/
+
+  // The visitor must implement
+  //
+  //
+  // viz.on_flood_start(Level l, Point p): called when we start flooding a preak component at level l of p
+  //
+  // viz.on_flood_end(Level l): called when we end flooding the the peak component
+  //
+  // 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
+
+  template <class I, class N, class DFVisitor, class Proj>
+  void depthfirst(I& f, N nbh, DFVisitor& viz, image_point_t<I> start)
+  {
+    enum st { NONE = 0, INQUEUE = 1, DONE = 2 };
+    using P = image_point_t<I>;
+    //using V = image_value_t<I>;
+
+    // FIXME
+    image_build_error_code err = IMAGE_BUILD_OK;
+
+    auto status = imchvalue<uint8_t>(f) //
+                      .adjust(nbh)
+                      .set_init_value(NONE)
+                      .get_status(&err)
+                      .build();
+
+    if (err != IMAGE_BUILD_OK)
+      throw std::runtime_error("Unable to have an extension.");
+
+    status.extension().fill(DONE);
+
+    mln::morpho::experimental::details::pqueue_fifo<I> queue(f);
+
+    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);
+      viz.on_flood_start(current_level, p);
+
+    flood_flat_zone:
+      for (auto n : nbh(p))
+      {
+        if (status(n) != NONE)
+          continue;
+
+        // Insert n INQUEUE
+        auto nval = f(n);
+        status(n) = INQUEUE;
+        queue.push(nval, n);
+
+        // If the neighbor is lower, postpone the neighbor
+        if (nval <= current_level)
+          continue;
+
+        // Otherwise, process it, (do not remove p from stack)
+        current_level = nval;
+        p             = n;
+        goto flood_new_level;
+      }
+
+      // All the neighbors have been seen, p is DONE
+      status(p) = DONE;
+      viz.on_done(current_level, p);
+      queue.pop();
+
+      // If the queue gets empty, we have processed the whole image
+      if (!queue.empty())
+      {
+        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;
+      }
+
+    // End: there is no more point to process
+      viz.on_flood_end(current_level);
+    }
+  }
+}

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

+#pragma once
+
+#include <mln/morpho/experimental/private/hvector.hpp>
+
+namespace  mln::morpho::experimental::detail
+{
+
+  /// \class hpqueue
+  /// \brief Hierarchical Priority Queue
+  /// Set of points ordered by increasing priority (highest first by default, lowest first if reverse=true) into a list of points
+  ///
+  /// \tparam N                 Number of levels
+  /// \tparam P                 Type of point
+  /// \tparam LinkImage         Type of the link image
+  template <int N, class P, class LinkImage, bool reverse>
+  class hpqueue : private hlinked_lists<N, P, LinkImage>
+  {
+    using base = hlinked_lists<N, P, LinkImage>;
+
+  public:
+    template <class J>
+    hpqueue(J&& f)
+      : hlinked_lists<N, P, LinkImage>(std::forward<J>(f))
+    {
+    }
+
+    // Push in last position at given level
+    void push_last(int level, P p) noexcept
+    {
+      mln_precondition(0 <= level && level < N);
+
+      base::push_back(level, p);
+      if constexpr (reverse)
+        m_current_level = std::min(m_current_level, level);
+      else
+        m_current_level = std::max(m_current_level, level);
+    }
+
+    // Push in first position at givel level
+    void push_first(int level, P p) noexcept
+    {
+      mln_precondition(0 <= level && level < N);
+
+      base::push_front(level, p);
+      if constexpr (reverse)
+        m_current_level = std::min(m_current_level, level);
+      else
+        m_current_level = std::max(m_current_level, 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)}; }
+
+    void pop() noexcept;
+
+  private:
+    int m_current_level = reverse ? (N - 1) : 0;
+  };
+
+
+
+  /******************************************/
+  /****          Implementation          ****/
+  /******************************************/
+
+  template <int N, class P, class LinkImage, bool reverse>
+  inline void hpqueue<N, P, LinkImage, reverse>::pop() noexcept
+  {
+    base::pop_front(m_current_level);
+    if (!base::empty(m_current_level))
+      return;
+
+    // Try to go down
+    if constexpr (reverse)
+    {
+      while (base::empty(m_current_level) && m_current_level < (N - 1))
+        m_current_level++;
+    }
+    else
+    {
+      while (base::empty(m_current_level) && m_current_level > 0)
+        m_current_level--;
+    }
+  }
+
+
+} // namespace mln::morpho::experimental::detail

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

+#pragma once
+
+#include <mln/core/concept/new/images.hpp>
+
+#include <mln/core/assert.hpp>
+
+namespace  mln::morpho::experimental::detail
+{
+
+
+  /// \brief Set of points ordered hierarchically into a list of points
+  ///
+  /// \tparam N                 Number of levels
+  /// \tparam P                 Type of point
+  /// \tparam LinkImage         Type of the link image
+  template <int N, class P, class LinkImage>
+  class hlinked_lists
+  {
+    static_assert(N <= (1 << 24), "Too many numbers of level");
+    static_assert(std::is_convertible_v<P, image_point_t<LinkImage>>);
+
+    template <class J>
+    hlinked_lists(J&& f);
+
+
+    void                push_front(int level, P p) noexcept;
+    void                push_back(int level, P p) noexcept;
+    P                   pop_front(int level) noexcept;
+
+    P                   back(int level) const noexcept;
+    P                   front(int level) const noexcept;
+    bool                empty(int level) const noexcept;
+
+
+  private:
+    struct node_t
+    {
+      P    head;
+      P    tail;
+      int  size = 0;
+    };
+
+    std::array<node_t, N> m_lists;
+    LinkImage             m_next;
+  };
+
+  /******************************************/
+  /****          Implementation          ****/
+  /******************************************/
+
+  template <int N, class P, class LinkImage>
+  template <class J>
+  inline hlinked_lists<N, P, LinkImage>::hlinked_lists(J&& f)
+    : m_next(std::forward<J>(f), image_build_params{})
+  {
+  }
+
+  template <int N, class P, class LinkImage>
+  inline bool hlinked_lists<N, P, LinkImage>::empty(int level) const noexcept
+  {
+    mln_precondition(level < N);
+    return m_lists[level].size == 0;
+  }
+
+
+  template <int N, class P, class LinkImage>
+  inline void hlinked_lists<N, P, LinkImage>::push_front(int level, P p) noexcept
+  {
+    mln_precondition(level < N);
+    if (m_lists[level].size == 0)
+    {
+      m_lists[level].head  = p;
+      m_lists[level].tail  = p;
+    }
+    else
+    {
+      m_next(p)            = m_lists[level].head;
+      m_lists[level].head  = p;
+    }
+
+    m_lists[level].size++;
+  }
+
+
+  template <int N, class P, class LinkImage>
+  inline void hlinked_lists<N, P, LinkImage>::push_back(int level, P p) noexcept
+  {
+    mln_precondition(level < N);
+    if (m_lists[level].size == 0)
+    {
+      m_lists[level].head  = p;
+      m_lists[level].tail  = p;
+    }
+    else
+    {
+      m_next(m_lists[level].tail) = p;
+      m_lists[level].tail         = p;
+    }
+    m_lists[level].size++;
+  }
+
+  template <int N, class P, class LinkImage>
+  inline P hlinked_lists<N, P, LinkImage>::pop_front(int level) noexcept
+  {
+    mln_precondition(level < N);
+    mln_precondition(m_lists[level].size > 0 && "Empty list");
+    P head = m_lists[level].head;
+    m_lists[level].head = m_next(head);
+    m_lists[level].size--;
+
+    return head;
+  }
+
+
+  template <int N, class P, class LinkImage>
+  inline P hlinked_lists<N, P, LinkImage>::front(int level) const noexcept
+  {
+    mln_precondition(level < N);
+    mln_precondition(m_lists[level].size > 0 && "Empty list");
+    return m_lists[level].head;
+  }
+
+  template <int N, class P, class LinkImage>
+  inline P hlinked_lists<N, P, LinkImage>::back(int level) const noexcept
+  {
+    mln_precondition(level < N);
+    mln_precondition(m_lists[level].size > 0 && "Empty list");
+    return m_lists[level].tail;
+  }
+
+
+} // namespace mln::morpho::experimental::detail

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

@@ -2,14 +2,14 @@
 
 #include <mln/core/concept/new/images.hpp>
 
-#include <mln/morpho/experimental/private/pqueue_hqueue_fifo.hpp>
 
+#include <mln/morpho/experimental/private/hpqueue.hpp>
 
 namespace mln::morpho::experimental::details
 {
 
   /// Provides a priority queue of points with a the fifo property
-  template <class I>
+  template <class I, bool reversed = false>
   class pqueue_fifo
   {
     using key_type   = image_value_t<I>;
@@ -20,50 +20,63 @@ namespace mln::morpho::experimental::details
     pqueue_fifo(J&& f);
 
     void                            push(const key_type& priority, const value_type& element);
+    void                            push_first(const key_type& priority, const value_type& element);
     void                            pop();
     std::pair<key_type, value_type> top() const;
     bool                            empty() const;
 
   private:
+    static_assert(!std::is_signed_v<key_type>, "Must not be signed");
     static_assert(value_traits<key_type>::quant <= 16, "Only low quantized type supported.");
 
-    pqueue_hqueue_fifo<image_concrete_t<I>> m_delegate;
+    static constexpr int nvalues = std::numeric_limits<key_type>::max() + 1;
+
+    detail::hpqueue<nvalues, value_type, image_ch_value_t<I, value_type>, reversed> m_delegate;
   };
 
+
   /******************************************/
   /****          Implementation          ****/
   /******************************************/
 
-  template <class I>
+  template <class I, bool reverse>
   template <class J>
-  inline pqueue_fifo<I>::pqueue_fifo(J&& f)
+  inline pqueue_fifo<I, reverse>::pqueue_fifo(J&& f)
     : m_delegate{std::forward<J>(f)}
   {
   }
 
-  template <class I>
-  inline void pqueue_fifo<I>::push(const key_type& k, const value_type& v)
+  template <class I, bool reverse>
+  inline void pqueue_fifo<I, reverse>::push(const key_type& k, const value_type& v)
   {
-    m_delegate.push(k, v);
+    m_delegate.push_last(k, v);
   }
 
-  template <class I>
-  inline void pqueue_fifo<I>::pop()
+  template <class I, bool reverse>
+  inline void pqueue_fifo<I, reverse>::push_first(const key_type& k, const value_type& v)
+  {
+    m_delegate.push_first(k, v);
+  }
+
+
+  template <class I, bool reverse>
+  inline void pqueue_fifo<I, reverse>::pop()
   {
     m_delegate.pop();
   }
 
-  template <class I>
-  inline bool pqueue_fifo<I>::empty() const
+  template <class I, bool reverse>
+  inline bool pqueue_fifo<I, reverse>::empty() const
   {
     return m_delegate.empty();
   }
 
 
-  template <class I>
-  inline auto pqueue_fifo<I>::top() const -> std::pair<key_type, value_type>
+  template <class I, bool reverse>
+  inline auto pqueue_fifo<I, reverse>::top() const -> std::pair<key_type, value_type>
   {
     return m_delegate.top();
   }
 
+
 } // namespace mln::morpho::experimental::details

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

@@ -44,7 +44,7 @@ namespace mln::morpho::experimental
       // Pixels in the border gets the status 0 (deja vu)
       // Pixels in the queue get -1
       // Pixels not in the queue get -2
-      mln::morpho::experimental::details::pqueue_fifo<I> pqueue(input);
+      mln::morpho::experimental::details::pqueue_fifo<I, /* reversed = */ true> pqueue(input);
       {
         output.extension().fill(kWaterline);
 

tests/morpho/CMakeLists.txt

@@ -20,3 +20,4 @@ add_core_test(${test_prefix}watershed        watershed.cpp)
 add_core_test(${test_prefix}area_filter      area_filter.cpp)
 add_core_test(${test_prefix}dynamic_filter   dynamic_filter.cpp)
 add_core_test(${test_prefix}ToS              tos.cpp tos_tests_helper.cpp)
+add_core_test(${test_prefix}depth_first      depthfirst.cpp)

tests/morpho/depthfirst.cpp

+#include <mln/io/experimental/imprint.hpp>
+#include <mln/morpho/experimental/canvas/depthfirst.hpp>
+
+#include <mln/core/image/experimental/ndimage.hpp>
+#include <mln/core/neighborhood/c4.hpp>
+#include <mln/core/image/view/operators.hpp>
+#include <mln/core/algorithm/all_of.hpp>
+
+#include <fixtures/ImageCompare/image_compare.hpp>
+
+#include <gtest/gtest.h>
+#include <mln/io/experimental/imprint.hpp>
+#include <functional>
+#include <fmt/core.h>
+
+struct ordervisitor
+{
+  void on_flood_start(int, mln::experimental::point2d) noexcept { ++nb_component_entering; }
+  void on_flood_end(int) noexcept { ++nb_component_exiting; }
+  void on_done(int, mln::experimental::point2d p) { out(p) = count++; cnt(p)++; }
+
+  int                                 nb_component_entering = 0;
+  int                                 nb_component_exiting = 0;
+  int                                 count = 0;
+  mln::experimental::image2d<uint8_t> out;
+  mln::experimental::image2d<int>     cnt;
+
+};
+
+
+TEST(Morpho, depthfirst_max)
+{
+  using namespace mln::view::ops;
+
+  const mln::experimental::image2d<uint8_t> input = {{10, 11, 11, 15, 16, 11, +2}, //
+                                                     {+2, 10, 10, 10, 10, 10, 10}, //
+                                                     {18, +2, 18, 19, 18, 14, +6}, //
+                                                     {16, +2, 16, 10, 10, 10, 10}, //
+                                                     {18, 16, 18, +2, +2, +2, +2}};
+
+
+  ordervisitor viz;
+  viz.out.resize(input.domain());
+  viz.cnt.resize(input.domain());
+
+  mln::morpho::experimental::canvas::depthfirst(input, mln::experimental::c4, viz, {0,0});
+
+  //mln::io::experimental::imprint(input >= 2); fmt::print("\n");
+  //mln::io::experimental::imprint(input >= 6); fmt::print("\n");
+  //mln::io::experimental::imprint(input >= 10); fmt::print("\n");
+  //mln::io::experimental::imprint(input >= 11); fmt::print("\n");
+  //mln::io::experimental::imprint(input >= 14); fmt::print("\n");
+  //mln::io::experimental::imprint(input >= 15); fmt::print("\n");
+  //mln::io::experimental::imprint(input >= 16); fmt::print("\n");
+  //mln::io::experimental::imprint(input >= 18); fmt::print("\n");
+  //mln::io::experimental::imprint(input >= 19); fmt::print("\n");
+
+
+  // Counting = anti-leveling
+  mln_foreach_new(auto p, input.domain())
+    for (auto q : mln::experimental::c4(p))
+      if (input.domain().has(q))
+      {
+        if (input(p) > input(q)) {
+          ASSERT_LT(viz.out(p), viz.out(q));
+        }
+      }
+
+  EXPECT_EQ(viz.nb_component_entering, viz.nb_component_exiting);
+  EXPECT_EQ(viz.nb_component_exiting, 13);
+  EXPECT_TRUE(mln::all_of(viz.cnt == 1));
+  //ASSERT_IMAGES_EQ_EXP(ref, res);
+}