Add meyer watershed algorithm (only small grayscale value for now).

mln/morpho/watershed.hpp

+#ifndef MLN_MORPHO_WATERSHED_HPP
+# define MLN_MORPHO_WATERSHED_HPP
+
+# include <mln/core/trace.hpp>
+# include <mln/core/image/image.hpp>
+# include <mln/core/neighborhood/neighborhood.hpp>
+# include <mln/core/extension/extension.hpp>
+# include <mln/core/extension/fill.hpp>
+# include <mln/morpho/private/pqueue.hpp>
+# include <mln/labeling/local_extrema.hpp>
+
+namespace mln
+{
+  namespace morpho
+  {
+
+    template <class Label_t, class I, class N>
+    mln_ch_value(I, Label_t) watershed(const Image<I>& ima, const Neighborhood<N>& nbh, int& nlabel);
+
+
+    /******************************************/
+    /****          Implementation          ****/
+    /******************************************/
+
+    namespace details
+    {
+
+
+      template <class I, class N, class O>
+      int watershed(const I& input, const N& nbh, O& output)
+      {
+        using Label_t = mln_value(O);
+        using V = mln_value(I);
+
+
+        // 1. Labelize minima (note that output in initialized to -1)
+        const int nlabel = mln::labeling::details::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
+        pqueue_fifo<I> pqueue(input);
+        {
+          mln::extension::fill(output, kWaterline);
+
+          mln_pixter(px, output);
+          mln_iter(nx, nbh(px));
+          mln_forall(px)
+          {
+            if (px->val() == 0)
+            {
+              bool is_local_min_neighbor = false;
+              mln_forall(nx)
+                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
+        {
+          mln_cpixel(I) pxIn(&input);
+          mln_pixel(O) pxOut(&output);
+          mln_iter(nxIn, nbh(pxIn));
+          mln_iter(nxOut, nbh(pxOut));
+
+          while (!pqueue.empty())
+          {
+            mln_point(I) p;
+            V level;
+            std::tie(level, p) = pqueue.top();
+
+            pxOut = output.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;
+            mln_forall(nxOut)
+            {
+              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;
+              pxIn = input.pixel(p);
+              mln_forall(nxIn, nxOut)
+              {
+                auto nlbl = nxOut->val();
+                if (nlbl == kUnlabeled)
+                {
+                  pqueue.push(nxIn->val(), nxIn->point());
+                  nxOut->val() = kInqueue;
+                }
+              }
+            }
+          }
+        }
+
+        // 3. Label all unlabeled pixels
+        {
+          mln_foreach(auto px, output.pixels())
+            if (px.val() < 0)
+              px.val() = kWaterline;
+        }
+
+        return nlabel;
+      }
+    }
+
+
+    template <class Label_t, class I, class N>
+    mln_ch_value(I, Label_t) watershed(const Image<I>& ima_, const Neighborhood<N>& nbh_, int& nlabel)
+    {
+      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");
+
+      const I& ima = exact(ima_);
+      const N& nbh = exact(nbh_);
+
+      constexpr Label_t kUninitialized = -1;
+      mln_ch_value(I, Label_t) output = imchvalue<Label_t>(ima).adjust(nbh).init(kUninitialized);
+      if (extension::need_adjust(output, nbh))
+      {
+        auto out = extension::add_value_extension(output, kUninitialized);
+        nlabel = details::watershed(ima, nbh, out);
+      }
+      else
+      {
+        mln::extension::fill(output, kUninitialized);
+        nlabel = details::watershed(ima, nbh, output);
+      }
+      return output;
+    }
+
+
+  } // end of namespace mln::morpho
+} // end of namespace mln
+
+#endif //!MLN_MORPHO_WATERSHED_HPP

tests/morpho/CMakeLists.txt

@@ -14,6 +14,7 @@ add_core_test(${test_prefix}opening          opening.cpp)
 add_core_test(${test_prefix}extinction       extinction.cpp)
 add_core_test(${test_prefix}median_filter    median_filter.cpp)
 add_core_test(${test_prefix}hit_or_miss      hit_or_miss.cpp)
+add_core_test(${test_prefix}watershed        watershed.cpp)
 add_core_test(${test_prefix}ToS              tos.cpp tos_tests_helper.cpp)
 
 target_link_libraries(${test_prefix}saturate    ${FreeImage_LIBRARIES})
@@ -22,4 +23,4 @@ target_link_libraries(${test_prefix}erode       ${FreeImage_LIBRARIES})
 target_link_libraries(${test_prefix}gradient    ${FreeImage_LIBRARIES})
 target_link_libraries(${test_prefix}opening     ${FreeImage_LIBRARIES})
 target_link_libraries(${test_prefix}median_filter ${FreeImage_LIBRARIES})
-
+target_link_libraries(${test_prefix}watershed    ${FreeImage_LIBRARIES})

tests/morpho/watershed.cpp

+#include <mln/core/image/image2d.hpp>
+#include <mln/core/neighb2d.hpp>
+#include <mln/morpho/watershed.hpp>
+
+#include <gtest/gtest.h>
+#include <tests/helpers.hpp>
+
+
+using namespace mln;
+
+TEST(Morpho, watershed_gray)
+{
+  const mln::image2d<mln::uint8> 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},
+    { 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);
+}
+
+TEST(Morpho, watershed_thick)
+{
+  const mln::image2d<mln::uint8> 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},
+    { 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);
+}
+
+