Merge branch 'development/erosion' into 'development/ranges'

Basic Morphological Operators: Erosion & Dilation

See merge request !69

.gitignore

@@ -2,6 +2,7 @@
 /cmake-build-*
 /doc/source/images/*.png
 /test_package/build
+/doc/source/_build
 .vs
 .idea
 *.pyc

bench/BMDilation.cpp

-#include <mln/core/extension/fill.hpp>
-#include <mln/core/image/image2d.hpp>
-#include <mln/io/imread.hpp>
-#include <mln/io/imsave.hpp>
+
+#include <mln/core/image/experimental/ndimage.hpp>
+#include <mln/io/experimental/imread.hpp>
 
 #include <benchmark/benchmark.h>
 
 #include <mln/core/se/disc.hpp>
 #include <mln/core/se/rect2d.hpp>
-#include <mln/morpho/structural/dilate.hpp>
+#include <mln/morpho/experimental/dilation.hpp>
 
 #include <fixtures/ImagePath/image_path.hpp>
 
-using namespace mln;
+
 
 class BMDilation : public benchmark::Fixture
 {
 public:
   BMDilation()
   {
-    io::imread(fixtures::ImagePath::concat_with_filename("lena.pgm"), m_input);
-    int nr = m_input.nrows();
-    int nc = m_input.ncols();
-    resize(m_output, m_input);
-    m_bytes = nr * nc * sizeof(uint8);
+    mln::io::experimental::imread(fixtures::ImagePath::concat_with_filename("lena.pgm"), m_input);
+    int nr = m_input.width();
+    int nc = m_input.height();
+    mln::resize(m_output, m_input);
+    m_size = nr * nc;
   }
 
-  template <class SE, class Compare = mln::productorder_less<mln::uint8>>
-  void run(benchmark::State& st, const SE& se, Compare cmp = Compare())
+
+  void run(benchmark::State& st, std::function<void()> callback)
   {
     for (auto _ : st)
-      mln::morpho::structural::dilate(m_input, se, m_output, cmp);
-    st.SetBytesProcessed(st.iterations() * m_bytes);
+      callback();
+    st.SetBytesProcessed(int64_t(st.iterations()) * int64_t(m_size));
   }
 
 protected:
-  image2d<uint8> m_input;
-  image2d<uint8> m_output;
-  std::size_t    m_bytes;
+  mln::experimental::image2d<uint8_t> m_input;
+  mln::experimental::image2d<uint8_t> m_output;
+  std::size_t    m_size;
 };
 
+class slow_disc : public mln::se_facade<slow_disc>
+{
+  using Base = mln::experimental::se::disc;
+public:
+  using category     = mln::dynamic_neighborhood_tag;
+  using incremental  = std::false_type;
+  using decomposable = std::false_type;
+  using separable    = std::false_type;
+
+  explicit slow_disc(float radius, Base::approx approximation = Base::PERIODIC_LINES_8)
+    : m_disc(radius, approximation)
+  {
+  }
+
+  /// \brief Return a range of SE offsets
+  auto offsets() const { return m_disc.offsets(); }
+
+  /// \brief Return a range of SE offsets before center
+  auto before_offsets() const { return m_disc.before_offsets(); }
+
+  /// \brief Return a range of SE offsets after center
+  auto after_offsets() const { return m_disc.after_offsets(); }
+
+  /// \brief Returns the extent radius
+  int radial_extent() const { return m_disc.radial_extent(); }
+
+  /// \brief Return the input ROI for 2D box.
+  auto compute_input_region(mln::experimental::box2d roi) const { return m_disc.compute_input_region(roi); }
+
+  /// \brief Return the output ROI for 2D box.
+  auto compute_output_region(mln::experimental::box2d roi) const { return m_disc.compute_output_region(roi); }
+
+private:
+  Base m_disc;
+};
+
+
 BENCHMARK_DEFINE_F(BMDilation, EuclideanDisc_naive)(benchmark::State& st)
 {
-  int           radius = st.range(0);
-  mln::se::disc se(radius, 0);
-  this->run(st, se, std::less<mln::uint8>());
+  int       radius = st.range(0);
+  slow_disc se(radius, mln::experimental::se::disc::EXACT);
+
+  auto f = [&]() { mln::morpho::experimental::dilation(m_input, se, m_output); };
+  this->run(st, f);
 }
 
 BENCHMARK_DEFINE_F(BMDilation, EuclideanDisc_incremental)(benchmark::State& st)
 {
-  int           radius = st.range(0);
-  mln::se::disc se(radius, 0);
-  this->run(st, se);
+  int  radius = st.range(0);
+  auto se     = mln::experimental::se::disc(radius, mln::experimental::se::disc::EXACT);
+
+  auto f = [&]() { mln::morpho::experimental::dilation(m_input, se, m_output); };
+  this->run(st, f);
 }
 
 BENCHMARK_DEFINE_F(BMDilation, ApproximatedDisc)(benchmark::State& st)
 {
-  int           radius = st.range(0);
-  mln::se::disc se(radius);
-  this->run(st, se);
+  int  radius = st.range(0);
+  auto se     = mln::experimental::se::disc(radius, mln::experimental::se::disc::PERIODIC_LINES_8);
+
+  auto f = [&]() { mln::morpho::experimental::dilation(m_input, se, m_output); };
+  this->run(st, f);
 }
 
 
 BENCHMARK_DEFINE_F(BMDilation, Square)(benchmark::State& st)
 {
-  int             width = 2 * st.range(0) + 1;
-  mln::se::rect2d se(width, width);
-  this->run(st, se);
+  int  radius = st.range(0);
+  auto se     = mln::experimental::se::rect2d(2 * radius + 1, 2 * radius + 1);
+  auto f = [&]() { mln::morpho::experimental::dilation(m_input, se, m_output); };
+  this->run(st, f);
 }
 
 constexpr int max_range = 2 << 6;

doc/source/core/se/rectangle.rst

 Rectangle
 =========
 
-.. doxygenclass:: mln::experimental::se::rect2d
+.. doxygenstruct:: mln::experimental::se::rect2d
    :members:
 
 

doc/source/morpho/dilation.rst

@@ -3,22 +3,18 @@
 Dilation
 ========
 
-Include :file:`<mln/morpho/structural/dilate.hpp>`
+Include :file:`<mln/morpho/dilation.hpp>`
 
-.. cpp:namespace:: mln::morpho::structural
+.. cpp:namespace:: mln::morpho
 
 
-#. .. cpp:function:: \
-      template <class InputImage, class StructuringElement> \
-      concrete_t<InputImage> dilate(const InputImage& ima, const StructuringElement& se)
+.. cpp:function:: \
+      Image{I} concrete_t<I> dilation(I image, StructuringElement se)
+      Image{I} concrete_t<I> dilation(I image, StructuringElement se, BorderManager bm)
+      void dilation(Image image, StructuringElement se, OutputImage out)
+      void dilation(Image image, StructuringElement se, BorderManager bm, OutputImage out)
 
-#. .. cpp:function:: \
-      template <class InputImage, class StructuringElement, class Compare> \
-      concrete_t<InputImage> dilate(const InputImage& ima, const StructuringElement& se, Compare cmp)
 
-#. .. cpp:function:: \
-      template <class InputImage, class StructuringElement, class OutputImage, class Compare> \
-      void dilate(const InputImage& ima, const StructuringElement& se, OutputImage& output, Compare cmp)
 
       Dilation by a structuring element.
 
@@ -28,21 +24,20 @@ Include :file:`<mln/morpho/structural/dilate.hpp>`
       .. math::
          \delta(f)(x) = \bigvee \{ \, f(y), y \in B_x \, \}
 
-      * (2,3) If a optional \p cmp function is provided, the algorithm will internally do
-        an unqualified call to ``inf(x, y,cmp)``.The default is the product-order so
-        that it works for vectorial type as well.
+      * An optional border management may be used to manage border side-effects.
+        Only *fill* and *user* are currently supported.
 
-      * (3) If the optional ``output`` image is provided, it must be wide enough to store
-        the results (the function does not perform any resizing).
+      * If the optional ``output`` image is provided, it must be wide enough to store
+        the result (the function does not perform any resizing).
 
       :param ima: Input image 𝑓
       :param se:  Structuring element 𝐵
-      :param cmp (optional): Comparison function
+      :param bm (optional): Border manager
       :param output (optional): Output image
 
       :return:
          * (1,2) An image whose type is deduced from the input image
-         * (3\) Nothing (the output image is passed as an argument)
+         * (3,4) Nothing (the output image is passed as an argument)
 
       :exception: N/A
 
@@ -61,13 +56,13 @@ Example 1 : Dilation by a square on a gray-level image
 
 .. code-block:: cpp
 
-   #include <mln/morpho/structural/dilate.hpp>
+   #include <mln/morpho/dilation.hpp>
    #include <mln/core/wind2d.hpp>
 
    // Define a square SE of size 21x21
    auto input = ...;
-   mln::se::rect2d rect(21,21);
-   auto output = mln::morpho::structural::dilate(input, rect);
+   auto rect = mln::se::rect2d(21,21);
+   auto output = mln::morpho::dilation(input, rect);
 
 
 .. image:: /images/lena_gray.jpg

doc/source/morpho/erosion.rst

 Erosion
 =======
 
-Include :file:`<mln/morpho/structural/erode.hpp>`
+Include :file:`<mln/morpho/erosion.hpp>`
 
-#. .. cpp:function:: \
-      template <class InputImage, class StructuringElement> \
-      concrete_t<InputImage> erode(const InputImage& ima, const StructuringElement& se)
 
-#. .. cpp:function:: \
-      template <class InputImage, class StructuringElement, class Compare> \
-      concrete_t<InputImage> erode(const InputImage& ima, const StructuringElement& se, Compare cmp)
-
-#. .. cpp:function:: \
-      template <class InputImage, class StructuringElement, class OutputImage, class Compare> \
-      void erode(const InputImage& ima, const StructuringElement& se, OutputImage& output, Compare cmp)
+.. cpp:function:: \
+      Image{I} concrete_t<I> erosion(I image, StructuringElement se)
+      Image{I} concrete_t<I> erosion(I image, StructuringElement se, BorderManager bm)
+      void erosion(Image image, StructuringElement se, OutputImage out)
+      void erosion(Image image, StructuringElement se, BorderManager bm, OutputImage out)
 
       Erosion by a structuring element.
 
@@ -23,21 +18,20 @@ Include :file:`<mln/morpho/structural/erode.hpp>`
       .. math::
          \varepsilon(f)(x) = \bigwedge \{ \, f(y), y \in B_x \, \}
 
-      * (2,3) If a optional \p cmp function is provided, the algorithm will internally do
-        an unqualified call to ``inf(x, y,cmp)``.The default is the product-order so
-        that it works for vectorial type as well.
+      * An optional border management may be used to manage border side-effects.
+        Only *fill* and *user* are currently supported.
 
-      * (3) If the optional ``output`` image is provided, it must be wide enough to store
-        the results (the function does not perform any resizing).
+      * If the optional ``output`` image is provided, it must be wide enough to store
+        the result (the function does not perform any resizing).
 
       :param ima: Input image 𝑓
       :param se:  Structuring element 𝐵
-      :param cmp (optional): Comparison function
+      :param bm (optional): Border manager
       :param output (optional): Output image
 
       :return:
          * (1,2) An image whose type is deduced from the input image
-         * (3\) Nothing (the output image is passed as an argument)
+         * (3,4) Nothing (the output image is passed as an argument)
 
       :exception: N/A
 
@@ -56,13 +50,13 @@ Example 1 : Erosion by a square on a gray-level image
 
 .. code-block:: cpp
 
-   #include <mln/morpho/structural/erode.hpp>
+   #include <mln/morpho/erosion.hpp>
    #include <mln/core/wind2d.hpp>
 
    // Define a square SE of size 21x21
    auto input = ...;
-   auto rect = mln::make_rectangle2d(21,21);
-   auto output = mln::morpho::structural::erode(input, rect);
+   auto rect = mln::se::rect2d(21,21);
+   auto output = mln::morpho::erosion(input, rect);
 
 
 .. image:: /images/lena_gray.jpg

doc/source/ruminations/border_management.rst

+Taxonomie des algorithmes
+#########################
+
+* Global
+* Local (Stencil)
+* Point-wise (Map)
+
+Role du Border Manager
+======================
+
+Soit:
+
+* `WR` la ROI de travail, sur laquelle on souhaite obtenir le résultat de l'opérateur
+* `IR` la ROI d'input, dont les valeurs sont nécéssaires pour calculer les valeurs de `OR`.
+
+Dans le cas d'un opérator local, les deux valeurs sont liées par:
+
+* ``IR = se.compute_input_region(WR)``
+* ``WR = se.compute_output_region(IR)``
+
+Génerateur d'input
+******************
+
+Le border manager à alors pour rôle de faire en sorte qu'une image d'entrée soit *lisible* sur `IR`
+en laissant à l'utilisateur la possiblité de configurer la façon dont les valeurs de IR. Les configurations
+possibles de gestion sont:
+
+* user: les données proviennent de l'extension de l'image, gérées par l'utilisateur
+* fill(v):  les valeurs prendront toujours la valeur `v`
+* periodize: En 2D, ``managed(x,y) = f(x % N, x % M)`` 
+* mirror: En 2D: ``managed(x) = f(N - |(x % 2N) - N|, M - |y % 2M - M|)`` 
+
+A noter que tous les modes ne sont pas compatibles avec toutes les images.
+
+Lorsqu'on éxécute ``bm.manage(ima, se, WR)``, l'algorithme suivant est appliqué.
+
+#. Si l'image d'entrée possède une extension, et si l' extension est suffisament large (``extension.includes(IR)``), et
+   supporte le *mode* (``extension::suport_**mode**``), alors en fonction de *mode*, on éxécute:
+
+   * *rien* si *mode = user*
+   * ``extension.fill(v)`` si *mode = fill*
+   * ``extension.periodize()`` si *mode = periodize*
+   * ``extension.mirror()`` si *mode = mirror*
+
+#. Sinon si, le mode n'est pas *user* et que la gestion est *auto*, une image avec une extension artificielle est
+   ajoutée en fonction du mode:
+
+   * view::value_extended
+   * view::periodize_extended
+   * view::mirror_extended
+
+#. Sinon, c'est une erreur.
+
+Generateur d'image temporaire
+*****************************
+
+Certains algorithmes ont besoin d'images temporaires qui seront écrites. Le BM doit donc être capable de créer une
+nouvelle image avec des valeurs d'extension matérialisée en mémoire::
+
+    output = bm.create_temporary(ima, se, WR)
+
+Crée une image dont le domaine sera ``IR = se.compute_input_region(WR)`` et dont les valeurs étendues seront calculées en fonction du mode.
+
+
+Role du canvas algorithmique local
+==================================
+
+Il y a deux cas à penser:
+df
+

pylene/CMakeLists.txt

@@ -52,6 +52,7 @@ target_sources(Pylene PRIVATE
                src/core/image_format.cpp
                src/core/ndbuffer_image_data.cpp
                src/core/ndbuffer_image.cpp
+               src/core/traverse2d.cpp
                src/io/freeimage_plugin.cpp
                src/io/io.cpp
                src/io/imread.cpp

pylene/include/mln/accu/accumulators/h_infsup.hpp

@@ -65,7 +65,9 @@ namespace mln
         typedef T                                                                                   result_type;
         typedef boost::mpl::set<features::h_sup, features::h_inf, features::inf<>, features::sup<>> provides;
 
-        h_infsup_base()
+        using has_untake = std::true_type;
+
+        constexpr h_infsup_base()
           : m_inf(value_traits<T>::sup())
           , m_sup(value_traits<T>::inf())
           , m_count(0)

pylene/include/mln/accu/accumulators/infsup.hpp

@@ -209,7 +209,7 @@ namespace mln
         typedef T argument_type;
         typedef T result_type;
 
-        inf(const Compare& cmp = Compare())
+        constexpr inf(const Compare& cmp = Compare())
           : m_cmp(cmp)
           , m_inf(value_traits<T, Compare>::sup())
         {
@@ -250,7 +250,7 @@ namespace mln
         typedef T argument_type;
         typedef T result_type;
 
-        sup(const Compare& cmp = Compare())
+        constexpr sup(const Compare& cmp = Compare())
           : m_cmp(cmp)
           , m_sup(value_traits<T, Compare>::inf())
         {

pylene/include/mln/core/algorithm/copy.hpp

@@ -3,6 +3,7 @@
 #include <mln/core/image/image.hpp>
 #include <mln/core/rangev3/rows.hpp>
 #include <mln/core/rangev3/view/zip.hpp>
+#include <mln/core/trace.hpp>
 
 #include <range/v3/algorithm/copy.hpp>
 
@@ -110,6 +111,8 @@ namespace mln
     template <class InputImage, class OutputImage>
     void copy(InputImage input, OutputImage output)
     {
+      mln_entering("mln::copy");
+
       // FIXME: Add a precondition about the size of the domain ::ranges::size
       static_assert(mln::is_a<InputImage, Image>());
       static_assert(mln::is_a<OutputImage, Image>());

pylene/include/mln/core/algorithm/paste.hpp

 #pragma once
+#include <mln/core/algorithm/copy.hpp>
 #include <mln/core/concept/new/images.hpp>
+#include <mln/core/rangev3/foreach.hpp>
 #include <mln/core/rangev3/rows.hpp>
 
+#include <mln/core/trace.hpp>
+
 namespace mln
 {
 
@@ -30,23 +34,85 @@ namespace mln
   template <class InputImage, class OutputImage>
   void paste(InputImage src, OutputImage dest);
 
+  template <class InputImage, class InputRange, class OutputImage>
+  void paste(InputImage src, InputRange roi, OutputImage dest);
+
   /******************************************/
   /****          Implementation          ****/
   /******************************************/
 
+  namespace details
+  {
+    template <class I, class D = image_domain_t<I>, class = void>
+    struct is_image_clippable : std::false_type
+    {
+    };
+
+    template <class I, class D>
+    struct is_image_clippable<I, D, std::void_t<decltype(std::declval<I>().clip(std::declval<D>()))>> : std::true_type
+    {
+    };
+
+    template <class I, class D>
+    inline constexpr bool is_image_clippable_v = is_image_clippable<I, D>::value;
+  }
+
+  namespace impl
+  {
+    // Like paste but allows access to the extension
+    template <class InputImage, class InputRange, class OutputImage>
+    void paste_unsafe(InputImage src, InputRange roi, OutputImage dest)
+    {
+      mln_foreach_new(auto p, roi)
+        dest.at(p) = src.at(p);
+    }
+  }
+
+
+  template <class InputImage, class InputRange, class OutputImage>
+  void paste(InputImage src, InputRange roi, OutputImage dest)
+  {
+    mln_entering("mln::paste");
+    // FIXME: Add a precondition about the domain inclusion
+    // FIXME: check OutputImage is accessible
+    static_assert(mln::is_a<InputImage, experimental::Image>());
+    static_assert(mln::is_a<OutputImage, experimental::Image>());
+    static_assert(std::is_convertible_v<image_value_t<InputImage>, image_value_t<OutputImage>>);
+
+    if constexpr (details::is_image_clippable_v<InputImage, InputRange> && details::is_image_clippable_v<OutputImage, InputRange>)
+    {
+      mln::experimental::copy(src.clip(roi), dest.clip(roi));
+    }
+    else
+    {
+      mln_foreach_new(auto p, roi)
+        dest(p) = src(p);
+    }
+  }
+
   template <class InputImage, class OutputImage>
   void paste(InputImage src, OutputImage dest)
   {
+    mln_entering("mln::paste");
     // FIXME: Add a precondition about the domain inclusion
     // FIXME: check OutputImage is accessible
     static_assert(mln::is_a<InputImage, experimental::Image>());
     static_assert(mln::is_a<OutputImage, experimental::Image>());
     static_assert(std::is_convertible_v<image_value_t<InputImage>, image_value_t<OutputImage>>);
 
-    auto&& pixels = src.new_pixels();
-    for (auto row : ranges::rows(pixels))
-      for (auto px : row)
-        dest(px.point()) = px.val();
+    using InputRange = image_domain_t<InputImage>;
+
+    if constexpr (details::is_image_clippable_v<OutputImage, InputRange>)
+    {
+      mln::experimental::copy(src, dest.clip(src.domain()));
+    }
+    else
+    {
+      auto&& pixels = src.new_pixels();
+      for (auto row : ranges::rows(pixels))
+        for (auto px : row)
+          dest(px.point()) = px.val();
+    }
   }
 
 } // namespace mln

pylene/include/mln/core/canvas/local_accumulation.hpp

 #pragma once
 
 #include <mln/core/canvas/local_algorithm.hpp>
+#include <mln/core/box.hpp>
 
 namespace mln::canvas
 {
 
-
   template <class Accu, class SE, class I, class J,
-            bool __incremental__ = SE::incremental::value&& Accu::has_untake::value>
+            bool __incremental__ = SE::incremental::value&& Accu::has_untake::value && details::is_domain_row_contiguous<image_domain_t<I>>::value>
   class LocalAccumulation;
 
 

pylene/include/mln/core/canvas/local_algorithm.hpp

@@ -5,6 +5,7 @@
 #include <mln/core/rangev3/rows.hpp>
 #include <mln/core/rangev3/view/zip.hpp>
 #include <mln/core/trace.hpp>
+#include <mln/core/box.hpp>
 
 namespace mln::canvas
 {
@@ -48,11 +49,31 @@ namespace mln::canvas
     virtual void EvalAfterLocalLoop(image_reference_t<I> pval_i, image_reference_t<J> pval_j)  = 0;
   };
 
+  // To be incremental:
+  // * the domain contiguous on rows
+  // * the SE must incremental
+
+  // FIXME:
+  // we need to know if a domain is "row-contiguous" that is if the iteration leads to a unit move.
+  namespace details
+  {
+    template <class Domain, class = void>
+    struct is_domain_row_contiguous : std::false_type
+    {
+    };
+
+    template <class Impl>
+    struct is_domain_row_contiguous<mln::experimental::_box<Impl>> : std::true_type
+    {
+    };
+  }
+
 
   template <class SE, class I, class J>
   class IncrementalLocalAlgorithm : public LocalAlgorithm<SE, I, J>
   {
     static_assert(SE::incremental::value, "The Structuring Element must be incremental.");
+    static_assert(details::is_domain_row_contiguous<image_domain_t<I>>::value, "The domain must be quite regular.");
 
   public:
     using IncrementalLocalAlgorithm::LocalAlgorithm::LocalAlgorithm;
@@ -92,6 +113,7 @@ namespace mln::canvas
   {
     if (!this->m_se.is_incremental())
     {
+      mln::trace::warn("[Performance] The accumulator is not incremental.");
       this->LocalAlgorithm<SE, I, J>::Execute();
       return;
     }

pylene/include/mln/core/canvas/private/traverse2d.hpp

+#pragma once
+
+#include <mln/core/box.hpp>
+
+#include <functional>
+
+
+namespace mln::canvas::details
+{