migrate watershed implementation

bench/BMMorphoBase.cpp

@@ -15,6 +15,7 @@
 #include <mln/morpho/experimental/median_filter.hpp>
 #include <mln/morpho/experimental/opening.hpp>
 #include <mln/morpho/experimental/reconstruction.hpp>
+#include <mln/morpho/experimental/watershed.hpp>
 
 #include <mln/labeling/experimental/local_extrema.hpp>
 
@@ -216,4 +217,16 @@ BENCHMARK_F(BMMorpho, minima)(benchmark::State& st)
 }
 
 
+BENCHMARK_F(BMMorpho, watershed)(benchmark::State& st)
+{
+  auto f = [](const image_t& input, image_t&) {
+    int nlabel;
+    mln::morpho::experimental::watershed<int8_t>(input, mln::experimental::c4, nlabel);
+    return nlabel;
+  };
+  this->run(st, f);
+}
+
+
+
 BENCHMARK_MAIN();

bench/BMMorphoBaseRef.cpp

@@ -15,6 +15,7 @@
 #include <mln/morpho/structural/dilate.hpp>
 #include <mln/morpho/structural/erode.hpp>
 #include <mln/morpho/structural/opening.hpp>
+#include <mln/morpho/watershed.hpp>
 
 #include <mln/labeling/local_extrema.hpp>
 
@@ -162,5 +163,15 @@ BENCHMARK_F(BMMorpho, minima)(benchmark::State& st)
   this->run(st, f);
 }
 
+BENCHMARK_F(BMMorpho, watershed)(benchmark::State& st)
+{
+  auto f = [](const image_t& input, image_t&) {
+    int nlabel;
+    mln::morpho::watershed<int8_t>(input, mln::c4, nlabel);
+    return nlabel;
+  };
+  this->run(st, f);
+}
+
 
 BENCHMARK_MAIN();

pylene/include/mln/core/image/experimental/private/ndbuffer_image.hpp

@@ -615,8 +615,8 @@ namespace mln
     for (int k = 0; k < N; ++k)
     {
       coords[k]          = this->__axes(k).domain_begin;
-      shp[k]     = static_cast<short>(this->__axes(k).domain_end - this->__axes(k).domain_begin);
-      strides[k] = this->__axes(k).byte_stride / sizeof(T);
+      shp[N - k - 1]     = static_cast<short>(this->__axes(k).domain_end - this->__axes(k).domain_begin);
+      strides[N - k - 1] = this->__axes(k).byte_stride;
     }
 
     auto ptr = Impl::get_pointer(this->__info(), coords);

pylene/include/mln/io/experimental/imprint.hpp

@@ -52,7 +52,7 @@ namespace mln::io::experimental
       template <class I>
       struct impl_t : impl_base_t
       {
-        impl_t(I& x) { this->m_ima = &x; }
+        impl_t(I& x) { this->m_ima = (void*)(&x); }
         ~impl_t() final = default;
         std::size_t formatted_size(P p) const final { return fmt::formatted_size("{}", at(p)); }
         void        print(P p, int width) const final { fmt::print("{:>{}d}", at(p), width); }

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

+#pragma once
+
+#include <mln/core/concept/new/images.hpp>
+
+#include <mln/morpho/experimental/private/pqueue_hqueue_fifo.hpp>
+
+
+namespace mln::morpho::experimental::details
+{
+
+  /// Provides a priority queue of points with a the fifo property
+  template <class I>
+  class pqueue_fifo
+  {
+    using key_type   = image_value_t<I>;
+    using value_type = image_point_t<I>;
+
+  public:
+    template <class J>
+    pqueue_fifo(J&& f);
+
+    void                            push(const key_type& priority, const value_type& element);
+    void                            pop();
+    std::pair<key_type, value_type> top() const;
+    bool                            empty() const;
+
+  private:
+    static_assert(value_traits<key_type>::quant <= 16, "Only low quantized type supported.");
+
+    pqueue_hqueue_fifo<image_concrete_t<I>> m_delegate;
+  };
+
+  /******************************************/
+  /****          Implementation          ****/
+  /******************************************/
+
+  template <class I>
+  template <class J>
+  inline pqueue_fifo<I>::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)
+  {
+    m_delegate.push(k, v);
+  }
+
+  template <class I>
+  inline void pqueue_fifo<I>::pop()
+  {
+    m_delegate.pop();
+  }
+
+  template <class I>
+  inline bool pqueue_fifo<I>::empty() const
+  {
+    return m_delegate.empty();
+  }
+
+
+  template <class I>
+  inline auto pqueue_fifo<I>::top() const -> std::pair<key_type, value_type>
+  {
+    return m_delegate.top();
+  }
+
+} // namespace mln::morpho::experimental::details

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

+#pragma once
+
+#include <mln/core/concept/new/images.hpp>
+
+#include <array>
+#include <limits>
+#include <type_traits>
+
+namespace mln::morpho::experimental::details
+{
+
+  template <class I>
+  class pqueue_hqueue_fifo
+  {
+    using key_type   = image_value_t<I>;
+    using value_type = image_point_t<I>;
+
+    static_assert(std::is_integral_v<key_type>);
+    static_assert(!std::is_signed<key_type>::value, "Key type must be unsigned.");
+    static_assert(sizeof(key_type) <= 2, "Only low-quantized value supported.");
+    static_assert(sizeof(key_type) < sizeof(int), "Unsupported architecure.");
+
+
+  public:
+    template <class J>
+    pqueue_hqueue_fifo(J&& f);
+
+    void                            push(const key_type& k, const value_type& v);
+    void                            pop();
+    std::pair<key_type, value_type> top() const;
+    bool                            empty() const;
+
+  private:
+    static constexpr int nvalues = std::numeric_limits<key_type>::max() + 1;
+    struct node_t
+    {
+      int        size = 0;
+      value_type head;
+      value_type tail;
+    };
+
+    // Pointer to head and tail element in the queue
+    std::array<node_t, nvalues> m_queues;
+
+    // Next pointer image
+    image_ch_value_t<I, value_type> m_next;
+
+    int m_current_level;
+  };
+
+  template <class I>
+  template <class J>
+  inline pqueue_hqueue_fifo<I>::pqueue_hqueue_fifo(J&& f)
+    : m_next(std::forward<J>(f), image_build_params{})
+    , m_current_level{std::numeric_limits<key_type>::max()}
+  {
+  }
+
+  template <class I>
+  inline void pqueue_hqueue_fifo<I>::pop()
+  {
+    mln_precondition(m_queues[m_current_level].size > 0);
+
+    auto& head = m_queues[m_current_level].head;
+    head       = m_next(head);
+    m_queues[m_current_level].size--;
+
+    if (m_queues[m_current_level].size > 0)
+      return;
+
+    // Try to go up
+    {
+      do {
+        m_current_level++;
+      } while (m_current_level < std::numeric_limits<key_type>::max() && m_queues[m_current_level].size == 0);
+    }
+  }
+
+  template <class I>
+  inline void pqueue_hqueue_fifo<I>::push(const key_type& k, const value_type& v)
+  {
+    if (m_queues[k].size == 0)
+    {
+      m_queues[k].head = v;
+      m_queues[k].tail = v;
+    }
+    else
+    {
+      m_next(m_queues[k].tail) = v;
+      m_queues[k].tail         = v;
+    }
+    m_queues[k].size++;
+
+    if (k < m_current_level)
+      m_current_level = k;
+  }
+
+  template <class I>
+  inline auto pqueue_hqueue_fifo<I>::top() const -> std::pair<key_type, value_type>
+  {
+    mln_precondition(m_queues[m_current_level].size > 0);
+    return std::make_pair(static_cast<key_type>(m_current_level), m_queues[m_current_level].head);
+  }
+
+  template <class I>
+  inline bool pqueue_hqueue_fifo<I>::empty() const
+  {
+    return m_queues[m_current_level].size == 0;
+  }
+
+} // namespace mln::morpho::experimental::details

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

+#pragma once
+
+#include <mln/core/algorithm/for_each.hpp>
+#include <mln/core/extension/fill.hpp>
+#include <mln/core/image/image.hpp>
+#include <mln/core/image/view/value_extended.hpp>
+#include <mln/core/neighborhood/neighborhood.hpp>
+#include <mln/core/trace.hpp>
+
+#include <mln/labeling/experimental/local_extrema.hpp>
+#include <mln/morpho/experimental/private/pqueue.hpp>
+
+
+namespace mln::morpho::experimental
+{
+
+  template <class Label_t, class InputImage, class Neighborhood>
+  image_ch_value_t<std::remove_reference_t<InputImage>, Label_t> //
+  watershed(InputImage&& ima, Neighborhood&& nbh, int& nlabel);
+
+
+  /******************************************/
+  /****          Implementation          ****/
+  /******************************************/
+
+  namespace impl
+  {
+
+
+    template <class I, class N, class O>
+    int watershed(I input, N nbh, O output)
+    {
+      using Label_t = image_value_t<O>;
+
+      // 1. Labelize minima (note that output is initialized to -1)
+      const int nlabel = mln::labeling::experimental::impl::local_minima(input, nbh, output, std::less<Label_t>());
+
+
+      constexpr int kUnlabeled = -2;
+      constexpr int kInqueue   = -1;
+      constexpr int kWaterline = 0;
+
+      // 2. inset neighbors inqueue
+      // 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);
+      {
+        output.extension().fill(kWaterline);
+
+        mln_foreach_new(auto px, output.new_pixels())
+        {
+          // Not a local minimum => early exit
+          if (px.val() != 0)
+            continue;
+
+          bool is_local_min_neighbor = false;
+          for (auto nx : nbh(px))
+            if (nx.val() > 0) // p is neighbhor to a local minimum
+            {
+              is_local_min_neighbor = true;
+              break;
+            }
+          if (is_local_min_neighbor)
+          {
+            px.val() = kInqueue;
+            pqueue.push(input(px.point()), px.point());
+          }
+          else
+          {
+            px.val() = kUnlabeled;
+          }
+        }
+      }
+
+      // 3. flood from minima
+      {
+        while (!pqueue.empty())
+        {
+          auto [level, p] = pqueue.top();
+
+          auto pxOut = output.new_pixel(p);
+          mln_assertion(pxOut.val() == kInqueue);
+          pqueue.pop();
+
+          // Check if there is a single marked neighbor
+          Label_t common_label               = kWaterline;
+          bool    has_single_adjacent_marker = false;
+          for (auto nxOut : nbh(pxOut))
+          {
+            int nlbl = nxOut.val();
+            if (nlbl <= 0)
+              continue;
+            else if (common_label == kWaterline)
+            {
+              common_label               = nlbl;
+              has_single_adjacent_marker = true;
+            }
+            else if (nlbl != common_label)
+            {
+              has_single_adjacent_marker = false;
+              break;
+            }
+          }
+
+          if (!has_single_adjacent_marker)
+          {
+            // If there are multiple labels => waterline
+            pxOut.val() = kWaterline;
+          }
+          else
+          {
+            // If a single label, it gets labeled
+            // Add neighbors in the queue
+            pxOut.val() = common_label;
+            for (auto q : nbh(p))
+            {
+              auto nlbl = output.at(q);
+              if (nlbl == kUnlabeled)
+              {
+                pqueue.push(input(q), q);
+                output(q) = kInqueue;
+              }
+            }
+          }
+        }
+      }
+
+      // 3. Label all unlabeled pixels
+      {
+        mln::for_each(output, [](auto& v) {
+          if (v < 0)
+            v = kWaterline;
+        });
+      }
+
+      return nlabel;
+    }
+  }
+
+  template <class Label_t, class InputImage, class Neighborhood>
+  image_ch_value_t<std::remove_reference_t<InputImage>, Label_t> //
+  watershed(InputImage&& ima, Neighborhood&& nbh, int& nlabel)
+  {
+    using I = std::remove_reference_t<InputImage>;
+    using N = std::remove_reference_t<Neighborhood>;
+
+    static_assert(mln::is_a<I, mln::experimental::Image>());
+    static_assert(mln::is_a<N, mln::experimental::Neighborhood>());
+
+    static_assert(std::is_integral<Label_t>::value, "The label type must integral.");
+    static_assert(std::is_signed<Label_t>::value, "The label type must be signed.");
+
+    mln_entering("mln::morpho::watershed");
+
+    constexpr Label_t kUninitialized = -1;
+
+    image_build_error_code err = IMAGE_BUILD_OK;
+
+    auto output = imchvalue<Label_t>(ima) //
+      .adjust(nbh)
+      .set_init_value(kUninitialized)
+      .get_status(&err)
+      .build();
+
+    if (err == IMAGE_BUILD_OK)
+    {
+      nlabel = impl::watershed(ima, nbh, output);
+    }
+    else
+    {
+      mln::trace::warn("[Performance] The extension is not wide enough");
+      auto out = view::value_extended(output, kUninitialized);
+      nlabel   = impl::watershed(ima, nbh, out);
+    }
+    return output;
+  }
+
+} // namespace mln::morpho::experimental

tests/morpho/watershed.cpp

-#include <mln/core/image/image2d.hpp>
-#include <mln/core/neighb2d.hpp>
-#include <mln/morpho/watershed.hpp>
+#include <mln/morpho/experimental/watershed.hpp>
+
+#include <mln/core/image/experimental/ndimage.hpp>
+#include <mln/core/neighborhood/c4.hpp>
 
 #include <fixtures/ImageCompare/image_compare.hpp>
 
@@ -11,38 +12,38 @@ using namespace mln;
 
 TEST(Morpho, watershed_gray)
 {
-  const mln::image2d<mln::uint8> input = {{2, 2, 2, 2, 2},      //
+  const mln::experimental::image2d<uint8_t> input = {{2, 2, 2, 2, 2},      //
                                                      {40, 30, 30, 30, 40}, //
                                                      {40, 20, 20, 20, 40}, //
                                                      {40, 40, 20, 40, 40}, //
                                                      {1, 5, 20, 5, 1}};
 
-  const mln::image2d<mln::int16> ref = {{1, 1, 1, 1, 1}, //
+  const mln::experimental::image2d<int16_t> ref = {{1, 1, 1, 1, 1}, //
                                                    {1, 1, 1, 1, 1}, //
                                                    {0, 1, 1, 1, 0}, //
                                                    {2, 0, 1, 0, 3}, //
                                                    {2, 2, 0, 3, 3}};
 
   int  nlabel;
-  auto res = mln::morpho::watershed<mln::int8>(input, mln::c4, nlabel);
-  ASSERT_IMAGES_EQ(ref, res);
+  auto res = mln::morpho::experimental::watershed<int16_t>(input, mln::experimental::c4, nlabel);
+  ASSERT_IMAGES_EQ_EXP(ref, res);
 }
 
 TEST(Morpho, watershed_thick)
 {
-  const mln::image2d<mln::uint8> input = {{0, 10, 0, 10, 0},    //
+  const mln::experimental::image2d<uint8_t> input = {{0, 10, 0, 10, 0},    //
                                                      {0, 10, 10, 10, 10},  //
                                                      {10, 10, 10, 10, 10}, //
                                                      {0, 10, 10, 10, 10},  //
                                                      {0, 10, 0, 10, 0}};
 
-  const mln::image2d<mln::int16> ref = {{1, 0, 2, 0, 3}, //
+  const mln::experimental::image2d<int16_t> ref = {{1, 0, 2, 0, 3}, //
                                                    {1, 1, 0, 3, 3}, //
                                                    {0, 0, 0, 0, 0}, //
                                                    {4, 4, 0, 6, 6}, //
                                                    {4, 0, 5, 0, 6}};
 
   int  nlabel;
-  auto res = mln::morpho::watershed<mln::int16>(input, mln::c4, nlabel);
-  ASSERT_IMAGES_EQ(ref, res);
+  auto res = mln::morpho::experimental::watershed<int16_t>(input, mln::experimental::c4, nlabel);
+  ASSERT_IMAGES_EQ_EXP(ref, res);
 }