Implement Opening By Reconstruction

bench/BMMorphoBase.cpp

 #include <mln/core/algorithm/transform.hpp>
 #include <mln/core/colors.hpp>
 #include <mln/core/image/experimental/ndimage.hpp>
+
 #include <mln/core/se/disc.hpp>
 #include <mln/core/se/rect2d.hpp>
 #include <mln/core/se/mask2d.hpp>
+#include <mln/core/neighborhood/c4.hpp>
 #include <mln/io/experimental/imread.hpp>
+
 #include <mln/morpho/experimental/closing.hpp>
 #include <mln/morpho/experimental/dilation.hpp>
 #include <mln/morpho/experimental/erosion.hpp>
 #include <mln/morpho/experimental/hit_or_miss.hpp>
-#include <mln/morpho/experimental/opening.hpp>
 #include <mln/morpho/experimental/median_filter.hpp>
+#include <mln/morpho/experimental/opening.hpp>
+#include <mln/morpho/experimental/reconstruction.hpp>
 
 #include <benchmark/benchmark.h>
 
@@ -151,15 +155,31 @@ BENCHMARK_F(BMMorpho, Opening_Disc)(benchmark::State& st)
 }
 
 BENCHMARK_F(BMMorpho, Median_Filter_Disc)(benchmark::State& st)
+{
+  int  radius = 32;
+  auto se      = mln::experimental::se::rect2d(2 * radius + 1, 2 * radius + 1);
+  auto f       = [se](const image_t& input, image_t& output) {
+    mln::morpho::experimental::median_filter(input, se, mln::extension::bm::mirror{}, output);
+  };
+  this->run(st, f);
+}
+
+
+BENCHMARK_F(BMMorpho, Opening_By_Reconstruction_Disc)(benchmark::State& st)
 {
   int  radius = 32;
   auto se     = mln::experimental::se::disc(radius);
-  auto f = [se](const image_t& input, image_t& output) { mln::morpho::experimental::median_filter(input, se, mln::extension::bm::mirror{}, output); };
+  auto markers = mln::morpho::experimental::erosion(m_input, se);
+  auto f       = [m = std::move(markers)](const image_t& input, image_t& output) {
+    output = mln::morpho::experimental::opening_by_reconstruction(input, m, mln::experimental::c4);
+  };
   this->run(st, f);
 }
 
 
 
+
+
 BENCHMARK_F(BMMorpho, Hit_or_miss_corner)(benchmark::State& st)
 {
   mln::se::experimental::mask2d se_hit = {

bench/BMMorphoBaseRef.cpp

@@ -4,11 +4,16 @@
 #include <mln/core/se/disc.hpp>
 #include <mln/core/se/rect2d.hpp>
 #include <mln/core/se/mask2d.hpp>
+#include <mln/core/neighb2d.hpp>
+
 #include <mln/io/experimental/imread.hpp>
+
+#include <mln/morpho/hit_or_miss.hpp>
+#include <mln/morpho/median_filter.hpp>
+#include <mln/morpho/opening_by_reconstruction.hpp>
 #include <mln/morpho/structural/closing.hpp>
 #include <mln/morpho/structural/dilate.hpp>
 #include <mln/morpho/structural/erode.hpp>
-#include <mln/morpho/hit_or_miss.hpp>
 #include <mln/morpho/structural/opening.hpp>
 
 // [legacy]
@@ -103,12 +108,26 @@ BENCHMARK_REGISTER_F(BMMorpho, Dilation_Square)->RangeMultiplier(2)->Range(2, ma
 
 BENCHMARK_F(BMMorpho, Opening_Disc)(benchmark::State& st)
 {
-  int  radius = 32;
+  int  radius  = 32;
+  auto se      = mln::se::rect2d(2 * radius + 1, 2 * radius + 1);
+  auto f       = [se](const image_t& input, image_t& output) {
+    output = mln::morpho::structural::opening(input, se, mln::productorder_less<uint8_t>());
+  };
+  this->run(st, f);
+}
+
+BENCHMARK_F(BMMorpho, Opening_By_Reconstruction_Disc)(benchmark::State& st)
+{
+  int  radius  = 32;
   auto se     = mln::se::disc(radius);
-  auto f = [se](const image_t& input, image_t& output) { mln::morpho::structural::opening(input, se, mln::productorder_less<uint8_t>(), output); };
+  auto markers = mln::morpho::structural::erode(m_input_, se);
+  auto f       = [m = std::move(markers)](const image_t& input, image_t& output) {
+    output = mln::morpho::opening_by_reconstruction(input, m, mln::c4);
+  };
   this->run(st, f);
 }
 
+
 BENCHMARK_F(BMMorpho, Hit_or_miss_corner)(benchmark::State& st)
 {
   mln::se::mask2d se_hit = {

doc/source/morpho.rst

@@ -15,7 +15,6 @@ Structural Morphological Operation
    morpho/rank_filter
    morpho/median_filter
    morpho/opening_by_reconstruction
-   morpho/closing_by_reconstruction
 
 
 Segmentation

doc/source/morpho/closing_by_reconstruction.rst

-Closing by reconstruction
-=========================
-
-Include :file:`<mln/morpho/closing_by_reconstruction.hpp>`
-
-.. cpp:namespace:: mln::morpho
-
-
-.. cpp:function:: \
-   template <class InputImage1, class InputImage2, class Neighborhood> \
-   concrete_t<InputImage1> closing_by_reconstruction(const InputImage1& f, const InputImage2& markers, const Neighborhood& nbh)
-
-.. cpp:function:: \
-   template <class InputImage1, class InputImage2, class Neighborhood, class Compare> \
-   concrete_t<InputImage1> closing_by_reconstruction(const InputImage1& f, const InputImage2& markers, const Neighborhood& nbh, Compare cmp)
-
-   Perform the reconstruction of the image ``markers`` over the constrain
-   image ``f``. Contrary to :ref:`opening-by-reconstruction`, the *background* is
-   reconstructed instead of the *foreground*.
-
-
-   :param f: Input image 𝑓
-   :param markers:  Marker image (must have the same value type of 𝑓)
-   :param nbh:  Fundamental structuring element.
-   :param cmp (optional): Comparaison function defining a strict weak ordering of values.
-
-   :return: An image whose type is deduced from the input image
-
-   :precondition: :math:`markers > f`
-
-   :exception: N/A
-
-
-Notes
------
-
-
-Complexity
-----------
-

doc/source/morpho/hit_or_miss.rst

@@ -94,11 +94,10 @@ Hit or miss transform to detect horizontal 2px-thick line, with pattern::
            :width: 49%
 
 
-Logical or between the two previous images:
+Logical or between the two previous images::
 
-.. literalinclude:: /snippets/staff_lines.cpp
-   :lines: 50
-   :language: cpp
+  using mln::view::ops;
+  auto markers = markers1 || markers2;
 
 .. image:: /images/staff_lines.png
            :width: 49%

doc/source/morpho/median_filter.rst

@@ -43,7 +43,7 @@ Complexity
 
 
 Example 1 : Median-filter by a square on a gray-level image
----------------------------------------------------------
+-----------------------------------------------------------
 
 .. code-block:: cpp
 

doc/source/morpho/opening_by_reconstruction.rst

 .. _opening-by-reconstruction:
 
-Opening by reconstruction
-=========================
+Opening & Closing by reconstruction
+===================================
 
-Include :file:`<mln/morpho/opening_by_reconstruction.hpp>`
+Include :file:`<mln/morpho/reconstruction.hpp>`
 
 .. cpp:namespace:: mln::morpho
 
 .. cpp:function:: \
-   template <class InputImage1, class InputImage2, class Neighborhood> \
-   concrete_t<InputImage1> opening_by_reconstruction(const InputImage1& f, const InputImage2& markers, const Neighborhood& nbh)
-
-   Perform the reconstruction of the image ``markers`` under the constrain
-   image ``f``. The markers designate the parts of the image we want to
-   retain. In binary, it is equivalent to perform the *conditional*
-   dilation of the ``markers`` :math:`X`, by a structuring element 𝑩, using
-   the *reference* 𝑓 until reaching stability.
+   Image{I1}\
+   concrete_t<I1> opening_by_reconstruction(I1 f, I2 markers, Neighborhood nbh, Compare cmp)
+   Image{I1}\
+   concrete_t<I1> opening_by_reconstruction(I1 f, I2 markers, Neighborhood nbh)
+   Image{I1}\
+   concrete_t<I1> closing_by_reconstruction(I1 f, I2 markers, Neighborhood nbh)
+
+   The opening by reconstruction performs the reconstruction of the image
+   ``markers`` under the constrain image ``f``. The markers designate the parts
+   of the image we want to retain. In binary, it is equivalent to perform the
+   *conditional* dilation of the ``markers`` :math:`X`, by a structuring element
+   𝑩, using the *reference* 𝑓 until reaching stability.
 
    .. math::
       \delta^{n+1}_{f,\mathcal{B}}(X) = \delta_\mathcal{B}(\delta^{n}_{f,\mathcal{B}}(X)) \cap f
 
+   Similarly, the closing by reconstruction reconstructs the *background* instead of the *foreground*.
+
+
    :param f: Input image 𝑓
    :param markers:  Marker image (must have the same value type of 𝑓)
-   :param nbh:  Fundamental structuring element.
+   :param nbh:  Elementary structuring element.
 
    :return: An image whose type is deduced from the input image
 
-   :precondition: :math:`markers < f`
+   :precondition: :math:`markers \le f` (for (3) :math:`markers \ge f`)
 
    :exception: N/A
 
@@ -43,37 +50,88 @@ Complexity
 Example 1 : Staff lines reconstruction
 --------------------------------------
 
-.. image:: /images/staff_lines.png
-           :width: 49%
+.. list-table::
+
+   * -   .. figure:: /images/staff_lines.png
 
-.. image:: /images/staff_lines_markers.png
-           :width: 49%
+            Original image
 
-The markers have been obtained with the :cpp:func:`hit_or_miss`.
+     -   .. figure:: /images/staff_lines_markers.png
 
-* Reconstruction of the objects touching staff lines; with the foundamental SE
-  (4-connection). All objects that do not touch the staff lines are removed.
+            Markers obtained by the :doc:`hit_or_miss` transform.
+
+Given an original image and some markers obtained with the :doc:`hit_or_miss`
+transform. The geodesic reconstruction (with the 4-connection) of the original
+image by the markers give the objects touching staff lines. All objects that do
+not touch the staff lines are removed.
 
 .. literalinclude:: /snippets/staff_lines.cpp
-   :lines: 56
+   :start-after: M3_START
+   :end-before: M3_END
    :language: cpp
 
-.. image:: /images/staff_lines.png
-           :width: 49%
+.. figure:: /images/morpho_reconstruction_1.png
+   :figwidth: 49%
+   :figclass: align-center
 
-.. image:: /images/morpho_reconstruction_1.png
-           :width: 49%
+   Geodesic reconstruction from the markers.
 
-* Reconstruction of the lines only; with an horizontal SE `x-o-x`.
 
+If we want to reconstruct only the staff line only, use an horizontal SE `x-o-x`.
 
 .. literalinclude:: /snippets/staff_lines.cpp
-   :lines: 59
+   :start-after: M4_START
+   :end-before: M4_END
    :language: cpp
 
-.. image:: /images/staff_lines.png
-           :width: 49%
+.. figure:: /images/morpho_reconstruction_2.png
+   :figwidth: 49%
+   :figclass: align-center
+
+   Horizontal reconstruction from the markers.
+
+Example 2 : Dense region reconstruction
+---------------------------------------
+
+.. list-table::
+
+   *  -  .. figure:: /images/blobs2_binary.png
+
+            (a) Original image
+
+
+      -  .. figure:: /images/morpho_reconstruction_markers.png
+
+            (b) Markers from the :doc:`rank_filter`
+
+      -  .. figure:: /images/morpho_reconstruction_dilated.png
+
+            (c) Dilated of the original image (a)
+
+.. literalinclude:: /snippets/reconstruction.cpp
+   :start-after: M2_START
+   :end-before: M2_END
+   :language: cpp
+
+Given an original image. We first start with a :doc:`rank_filter` to locate
+dense region (regions with much more foreground pixels that background pixels)
+that gives us markers. Then a dilation with a small disc allows to connect
+objects. The reconstruction of the dilated image with a the markers gives a mask
+for the dense region. Finally, we just have to mask the input with the mask to
+get the objects in dense regions::
+
+   auto out = mln::clone(rec && input);
+
+.. list-table::
+
+   *  -  .. figure:: /images/morpho_reconstruction_rec.png
+   
+            (d) Reconstruction of (c) from the markers (b)
+
+      -  .. figure:: /images/morpho_reconstruction_out.png
+
+            Input (a) restricted to the mask (d)
+
+
 
-.. image:: /images/morpho_reconstruction_2.png
-           :width: 49%
 

doc/source/morpho/rank_filter.rst

@@ -6,7 +6,7 @@ Include :file:`<mln/morpho/rank_filter.hpp>`
 
 .. cpp:function:: \
     template <class Ratio> void rank_filter(Image image, StructuringElement se, BorderManager bm, OutputImage out)
-    template <class Ratio> Image{I} concrete_t<I> rank_filter(I image, StructuringElement se, BorderManager bm)
+    template <class Ratio, Image I> concrete_t<I> rank_filter(I image, StructuringElement se, BorderManager bm)
 
     
     The rank filter is non-linear filter that assigns the 𝑟-th value in a given structuring element 𝐵.

doc/source/snippets/CMakeLists.txt

@@ -28,7 +28,13 @@ add_image("erosion-cli;opening;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_
 add_image("erosion-cli;closing;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_closing_1.png)
 add_image("erosion-cli;median;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_median_1.png)
 add_image("staff_lines" "${DOCUMENTATION_IMAGE_DIR}/staff_lines.pbm"
-          morpho_hitormiss_1.png morpho_hitormiss_2.png staff_lines_markers.png morpho_reconstruction_1.png morpho_reconstruction_2.png)
+  morpho_hitormiss_1.png morpho_hitormiss_2.png staff_lines_markers.png morpho_reconstruction_1.png morpho_reconstruction_2.png)
+add_image("reconstruction"
+  "${DOCUMENTATION_IMAGE_DIR}/blobs2_binary.png"
+  morpho_reconstruction_dilated.png
+  morpho_reconstruction_markers.png
+  morpho_reconstruction_rec.png
+  morpho_reconstruction_out.png)
 add_image("blobs_watershed" "${DOCUMENTATION_IMAGE_DIR}/blobs_binary.png" blobs_distance_transform.png blobs_segmentation.png)
 
 
@@ -42,4 +48,4 @@ add_executable(erosion-cli erosion-cli.cpp)
 add_executable(staff_lines staff_lines.cpp)
 add_executable(component_tree_1 component_tree_1.cpp)
 add_executable(blobs_watershed blobs_watershed.cpp)
-
+add_executable(reconstruction reconstruction.cpp)

doc/source/snippets/reconstruction.cpp

+#include <mln/io/experimental/imread.hpp>
+#include <mln/io/experimental/imsave.hpp>
+
+
+#include <mln/core/image/experimental/ndimage.hpp>
+#include <mln/core/neighborhood/c8.hpp>
+#include <mln/core/se/disc.hpp>
+#include <mln/core/se/rect2d.hpp>
+#include <mln/morpho/experimental/dilation.hpp>
+#include <mln/morpho/experimental/erosion.hpp>
+#include <mln/morpho/experimental/rank_filter.hpp>
+#include <mln/morpho/experimental/reconstruction.hpp>
+#include <mln/core/image/view/operators.hpp>
+#include <ratio>
+
+int main(int argc, char** argv)
+{
+  if (argc < 6)
+  {
+    std::cerr << "Usage: " << argv[0] << " dilated.pbm markers.pbm reconstruction.pbm out.pbm\n";
+    return 1;
+  }
+
+  using namespace mln::view::ops;
+
+  mln::experimental::image2d<bool> input;
+
+  mln::io::experimental::imread(argv[1], input);
+
+  // #M2_START
+  // Make blobs connected
+  auto disc = mln::experimental::se::disc(4);
+  auto dil = mln::morpho::experimental::dilation(input, disc);
+
+  // Get markers for large connected components
+  auto rect = mln::experimental::se::rect2d(20, 20);
+  auto markers = mln::morpho::experimental::rank_filter<std::ratio<1,4>>(input, rect, mln::extension::bm::fill(false));
+
+  // Reconstruction of the large CC
+  auto rec = mln::morpho::experimental::opening_by_reconstruction(dil, markers, mln::experimental::c8);
+   // #M2_END
+
+  // Mask
+  auto out = mln::clone(rec && input);
+
+  // Save
+
+  mln::io::experimental::imsave(dil, argv[2]);
+  mln::io::experimental::imsave(markers, argv[3]);
+  mln::io::experimental::imsave(rec, argv[4]);
+  mln::io::experimental::imsave(out, argv[5]);
+}

doc/source/snippets/staff_lines.cpp

-#include <mln/core/image/image2d.hpp>
+#include <mln/core/image/experimental/ndimage.hpp>
 #include <mln/core/image/view/operators.hpp>
-#include <mln/core/neighb2d.hpp>
+#include <mln/core/neighborhood/c4.hpp>
+#include <mln/core/neighborhood/c2_h.hpp>
 #include <mln/core/se/mask2d.hpp>
-#include <mln/io/imread.hpp>
-#include <mln/io/imsave.hpp>
-#include <mln/morpho/hit_or_miss.hpp>
-#include <mln/morpho/opening_by_reconstruction.hpp>
+#include <mln/io/experimental/imread.hpp>
+#include <mln/io/experimental/imsave.hpp>
+#include <mln/morpho/experimental/hit_or_miss.hpp>
+#include <mln/morpho/experimental/reconstruction.hpp>
 
 
 int main(int argc, char** argv)
@@ -18,40 +19,43 @@ int main(int argc, char** argv)
     return 1;
   }
 
-  mln::image2d<bool> input;
-  mln::io::imread(argv[1], input);
+  mln::experimental::image2d<bool> input;
+  mln::io::experimental::imread(argv[1], input);
 
   // Make a dual neighborhood corresponding to
   // x x x
   // o o o
   // x x x
 
-  mln::image2d<bool> markers1, markers2;
+  mln::experimental::image2d<bool> markers1, markers2;
   { // #M1_START
-    mln::se::mask2d se_hit  = {{0, 0, 0}, {1, 1, 1}, {0, 0, 0}};
-    mln::se::mask2d se_miss = {{1, 1, 1}, {0, 0, 0}, {1, 1, 1}};
-    auto            output  = mln::morpho::hit_or_miss(input, se_hit, se_miss);
-    markers1                = output;
-    mln::io::imsave(mln::lnot(markers1), argv[2]);
+    mln::se::experimental::mask2d se_hit  = {{0, 0, 0}, {1, 1, 1}, {0, 0, 0}};
+    mln::se::experimental::mask2d se_miss = {{1, 1, 1}, {0, 0, 0}, {1, 1, 1}};
+
+    markers1 = mln::morpho::experimental::hit_or_miss(input, se_hit, se_miss);
+    mln::io::experimental::imsave(not markers1, argv[2]);
   } // #M1_END
 
   { // #M2_START
-    mln::se::mask2d se_hit  = {{0, 0, 0}, {0, 0, 0}, {1, 1, 1}, {1, 1, 1}, {0, 0, 0}};
-    mln::se::mask2d se_miss = {{0, 0, 0}, {1, 1, 1}, {0, 0, 0}, {0, 0, 0}, {1, 1, 1}};
-    auto            output  = mln::morpho::hit_or_miss(input, se_hit, se_miss);
-    markers2                = output;
-    mln::io::imsave(mln::lnot(markers2), argv[3]);
+    mln::se::experimental::mask2d se_hit  = {{0, 0, 0}, {0, 0, 0}, {1, 1, 1}, {1, 1, 1}, {0, 0, 0}};
+    mln::se::experimental::mask2d se_miss = {{0, 0, 0}, {1, 1, 1}, {0, 0, 0}, {0, 0, 0}, {1, 1, 1}};
+
+    markers2 = mln::morpho::experimental::hit_or_miss(input, se_hit, se_miss);
+    mln::io::experimental::imsave(not markers2, argv[3]);
   } // #M2_END
 
-  auto markers = mln::lor(markers1, markers2);
+  auto markers = markers1 or markers2;
 
-  mln::io::imsave(mln::lnot(markers), argv[4]);
+  mln::io::experimental::imsave(not markers, argv[4]);
 
-  // FIXME: add experimental version
-  auto markers_     = mln::lor(markers1, markers2);
-  auto all_touching = mln::morpho::opening_by_reconstruction(input, markers_, mln::c4);
-  mln::io::imsave(mln::lnot(all_touching), argv[5]);
+  // #M3_START
+  auto markers_     = markers1 or markers2;
+  auto all_touching = mln::morpho::experimental::opening_by_reconstruction(input, markers_, mln::experimental::c4);
+  mln::io::experimental::imsave(not all_touching, argv[5]);
+  // #M3_END
 
-  auto lines_only = mln::morpho::opening_by_reconstruction(input, markers_, mln::c2_h);
-  mln::io::imsave(mln::lnot(lines_only), argv[6]);
+  // #M4_START
+  auto lines_only = mln::morpho::experimental::opening_by_reconstruction(input, markers_, mln::experimental::c2_h);
+  mln::io::experimental::imsave(not lines_only, argv[6]);
+  // #M4_END
 }

pylene/include/mln/core/extension/border_management.hpp

@@ -239,6 +239,25 @@ namespace mln::extension
       mln::paste(ima, roi, output);
       return output;
     }
+
+    template <class I, class SE, class D>
+    I create_temporary_image(const SE& se, const D& roi) const
+    {
+      static_assert(mln::is_a<I, mln::experimental::Image>());
+      static_assert(mln::is_a<SE, mln::experimental::StructuringElement>());
+
+      if (m_value.type() != typeid(image_value_t<I>))
+        throw std::runtime_error(
+          "Trying to fill the border with a bad value type. Ensure that value type fits the image type.");
+
+      D input_roi = se.compute_input_region(roi);
+
+      image_build_params params;
+      params.init_value = m_value;
+
+      image_concrete_t<I> output(input_roi, params);
+      return output;
+    }
   };
 
   template <typename U>

pylene/include/mln/core/neighborhood/c2_h.hpp

+#pragma once
+
+#include <mln/core/neighborhood/private/neighborhood_facade.hpp>
+#include <mln/core/experimental/point.hpp>
+
+#include <array>
+#include <range/v3/span.hpp>
+
+namespace mln::experimental
+{
+  struct c2_h_t : neighborhood_facade<c2_h_t>
+  {
+  private:
+    using point_t = mln::experimental::ndpoint<2, std::ptrdiff_t>;
+
+  public:
+    using category     = constant_neighborhood_tag;
+    using incremental  = std::false_type;
+    using decomposable = std::false_type;
+    using separable    = std::false_type;
+
+    static constexpr ::ranges::span<const point_t, 2> offsets() { return {m_offsets.data(), 2}; }
+    static constexpr ::ranges::span<const point_t, 1> before_offsets() { return {m_offsets.data(), 1}; }
+    static constexpr ::ranges::span<const point_t, 1> after_offsets() { return {m_offsets.data() + 1, 1}; }
+
+    static constexpr int radial_extent() { return 1; }
+
+    /// \brief Return the input ROI for 2D box.
+    mln::experimental::box2d compute_input_region(mln::experimental::box2d roi) const
+    {
+      --roi.tl().x();
+      ++roi.br().x();
+      return roi;
+    }
+
+    /// \brief Return the output ROI for 2D box.
+    mln::experimental::box2d compute_output_region(mln::experimental::box2d roi) const
+    {
+      // Fixme: check emptiness
+      ++roi.tl().x();
+      --roi.br().x();
+      return roi;
+    }
+
+
+  private:
+    // clang-format off
+    static inline constexpr std::array<point_t, 2> m_offsets = {{
+        {-1, +0},          {+1, +0},
+      }};
+    // clang-format on
+  };
+
+  static constexpr inline c2_h_t c2_h = {};
+} // namespace mln::experimental

pylene/include/mln/core/neighborhood/c4.hpp

@@ -25,6 +25,21 @@ namespace mln::experimental
 
     static constexpr int radial_extent() { return 1; }
 
+    /// \brief Return the input ROI for 2D box.
+    mln::experimental::box2d compute_input_region(mln::experimental::box2d roi) const
+    {
+      roi.inflate(1);
+      return roi;
+    }
+
+    /// \brief Return the output ROI for 2D box.