Median filer implementation.

bench/BMMorphoBase.cpp

@@ -10,6 +10,7 @@
 #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 <benchmark/benchmark.h>
 
@@ -131,6 +132,7 @@ BENCHMARK_DEFINE_F(BMMorpho, Dilation_Square)(benchmark::State& st)
   this->run(st, f);
 }
 
+
 constexpr int max_range = 128;
 
 BENCHMARK_REGISTER_F(BMMorpho, Dilation_ApproximatedDisc)->RangeMultiplier(2)->Range(2, max_range);
@@ -148,6 +150,16 @@ BENCHMARK_F(BMMorpho, Opening_Disc)(benchmark::State& st)
   this->run(st, f);
 }
 
+BENCHMARK_F(BMMorpho, Median_Filter_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); };
+  this->run(st, f);
+}
+
+
+
 BENCHMARK_F(BMMorpho, Hit_or_miss_corner)(benchmark::State& st)
 {
   mln::se::experimental::mask2d se_hit = {

doc/source/morpho.rst

@@ -12,6 +12,8 @@ Structural Morphological Operation
    morpho/opening
    morpho/closing
    morpho/hit_or_miss
+   morpho/rank_filter
+   morpho/median_filter
    morpho/opening_by_reconstruction
    morpho/closing_by_reconstruction
 

doc/source/morpho/median_filter.rst

+Median filter
+=============
+
+Include :file:`<mln/morpho/median_filter.hpp>`
+
+
+.. cpp:function:: \
+    void median_filter(Image image, StructuringElement se, BorderManager bm, OutputImage out)
+    Image{I} concrete_t<I> median_filter(I image, StructuringElement se, BorderManager bm)
+
+    
+    The median filter is non-linear filter that assigns the median value in a given structuring element 𝐵.
+
+    .. math::
+        g(x) = med(\{ f(y) \in \mathcal{B}_x) \})
+
+    where `med` returns the median value of the set of pixels in the structuring element 𝑩 centered in 𝑥.
+
+    * A border management may be used to manage border side-effects.
+
+    * 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 bm: Border manager policy
+    :param output (optional): Output image
+
+    :return:
+       * (1) Nothing (the output image is passed as an argument)
+       * (2) An image whose type is deduced from the input image
+    
+    :exception: N/A
+
+
+Notes
+-----
+
+
+Complexity
+----------
+
+
+
+Example 1 : Median-filter by a square on a gray-level image
+---------------------------------------------------------
+
+.. code-block:: cpp
+
+   #include <mln/morpho/median_filter.hpp>
+   #include <mln/core/se/rect2d.hpp>
+
+   // Define a square SE of size 21x21
+   auto input = ...;
+   auto rect = mln::se::rect2d(21,21);
+
+   auto output = mln::morpho::median_filter(input, rect, mln::extension::mirror);
+
+
+.. image:: /images/lena_gray.jpg
+           :width: 49%
+
+.. image:: /images/morpho_median_1.png
+           :width: 49%

doc/source/morpho/opening_by_reconstruction.rst

@@ -22,7 +22,7 @@ Include :file:`<mln/morpho/opening_by_reconstruction.hpp>`
 
    :param f: Input image 𝑓
    :param markers:  Marker image (must have the same value type of 𝑓)
-   :param nbh:  Fundamental structiring element.
+   :param nbh:  Fundamental structuring element.
 
    :return: An image whose type is deduced from the input image
 

doc/source/morpho/rank_filter.rst

+Rank filter
+===========
+
+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)
+
+    
+    The rank filter is non-linear filter that assigns the 𝑟-th value in a given structuring element 𝐵.
+
+    .. math::
+        g(x) = r(\{ f(y) \in \mathcal{B}_x) \})
+
+    where `r` returns the 𝑟-th value of the set of pixels of the structuring element 𝑩 centered in 𝑥.
+
+    * A border management may be used to manage border side-effects.
+
+    * 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 bm: Border manager policy
+    :param output (optional): Output image
+
+    :return:
+       * (1) Nothing (the output image is passed as an argument)
+       * (2) An image whose type is deduced from the input image
+    
+    :exception: N/A
+
+
+Notes
+-----
+
+
+Complexity
+----------
+
+
+
+Example 1 : Rank-filter by a square on a gray-level image
+---------------------------------------------------------
+
+.. code-block:: cpp
+
+   #include <mln/morpho/rank_filter.hpp>
+   #include <mln/core/se/rect2d.hpp>
+
+   // Define a square SE of size 21x21
+   auto input = ...;
+   auto rect = mln::se::rect2d(21,21);
+
+   using R = std::ratio<1, 2>; // Get the median value
+   auto output = mln::morpho::rank_filter<R>(input, rect, mln::extension::mirror);
+
+
+.. image:: /images/lena_gray.jpg
+           :width: 49%
+
+.. image:: /images/morpho_median_1.png
+           :width: 49%

doc/source/snippets/CMakeLists.txt

@@ -26,6 +26,7 @@ add_image("erosion-cli;erosion;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_
 add_image("erosion-cli;dilation;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_dilation_1.png)
 add_image("erosion-cli;opening;square;21" "${PYLENE_IMAGE_DIR}/lena.pgm" morpho_opening_1.png)
 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)
 add_image("blobs_watershed" "${DOCUMENTATION_IMAGE_DIR}/blobs_binary.png" blobs_distance_transform.png blobs_segmentation.png)

doc/source/snippets/erosion-cli.cpp

-#include <mln/core/image/image2d.hpp>
+#include <mln/core/image/experimental/ndimage.hpp>
+
 #include <mln/core/se/rect2d.hpp>
-#include <mln/morpho/structural/closing.hpp>
-#include <mln/morpho/structural/dilate.hpp>
-#include <mln/morpho/structural/erode.hpp>
-#include <mln/morpho/structural/opening.hpp>
 
-#include <mln/io/imread.hpp>
-#include <mln/io/imsave.hpp>
+#include <mln/morpho/experimental/closing.hpp>
+#include <mln/morpho/experimental/dilation.hpp>
+#include <mln/morpho/experimental/erosion.hpp>
+#include <mln/morpho/experimental/opening.hpp>
+#include <mln/morpho/experimental/median_filter.hpp>
+
+#include <mln/io/experimental/imread.hpp>
+#include <mln/io/experimental/imsave.hpp>
+
 
 #include <algorithm>
 #include <boost/program_options.hpp>
@@ -50,7 +54,8 @@ enum morpho_op_type
   kErosion,
   kDilation,
   kOpening,
-  kClosing
+  kClosing,
+  kMedian,
 };
 
 std::istream& operator>>(std::istream& in, morpho_op_type& se)
@@ -66,6 +71,8 @@ std::istream& operator>>(std::istream& in, morpho_op_type& se)
     se = kOpening;
   else if (token == "closing")
     se = kClosing;
+  else if (token == "median")
+    se = kMedian;
   else
     throw po::invalid_option_value("Invalid Operator");
   return in;
@@ -116,27 +123,30 @@ int main(int argc, char** argv)
     return 1;
   }
 
-  mln::image2d<mln::uint8> input, output;
-  mln::io::imread(vm["input"].as<std::string>(), input);
+  mln::experimental::image2d<mln::uint8> input, output;
+  mln::io::experimental::imread(vm["input"].as<std::string>(), input);
 
 
-  mln::se::rect2d nbh(size, size);
+  mln::experimental::se::rect2d nbh(size, size);
 
   switch (vm["operator"].as<morpho_op_type>())
   {
   case kErosion:
-    output = mln::morpho::structural::erode(input, nbh);
+    output = mln::morpho::experimental::erosion(input, nbh);
     break;
   case kDilation:
-    output = mln::morpho::structural::dilate(input, nbh);
+    output = mln::morpho::experimental::dilation(input, nbh);
     break;
   case kOpening:
-    output = mln::morpho::structural::opening(input, nbh);
+    output = mln::morpho::experimental::opening(input, nbh);
     break;
   case kClosing:
-    output = mln::morpho::structural::closing(input, nbh);
+    output = mln::morpho::experimental::closing(input, nbh);
+    break;
+  case kMedian:
+    output = mln::morpho::experimental::median_filter(input, nbh, mln::extension::bm::mirror{});
     break;
   }
 
-  mln::io::imsave(output, vm["output"].as<std::string>());
+  mln::io::experimental::imsave(output, vm["output"].as<std::string>());
 }

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

+#pragma once
+
+#include <mln/morpho/experimental/rank_filter.hpp>
+
+#include <ratio>
+
+
+/// \file
+
+namespace mln::morpho::experimental
+{
+
+  /// \ingroup morpho
+  /// \brief Median filter
+  ///
+  /// A median filter is non-linerar filter using the value ordering to
+  /// compute the filtered value. It computes:
+  /// \f[
+  /// g(x) = med(\{ f(y) \in \mathcal{B}_x) \})
+  /// \f]
+  /// where \p med returns the median value of the set of pixels of the
+  /// structuring element 𝑩 centered in 𝑥.
+  ///
+  /// \note The algorithm internally uses the function rank_filter with
+  /// the ratio set to `std::ratio<1,2>`.
+  ///
+  /// \param ima Input image 𝑓
+  /// \param se  Structuring element
+  /// \param out (optional) Output image 𝑔 so store the result.
+  template <class InputImage, class SE, class BorderManager>
+  image_concrete_t<std::remove_reference_t<InputImage>>
+  median_filter(InputImage&& image, const mln::experimental::StructuringElement<SE>& se, BorderManager bm);
+
+  template <class InputImage, class SE, class BorderManager, class OutputImage>
+  void median_filter(InputImage&& image, const mln::experimental::StructuringElement<SE>& se, BorderManager bm,
+                     OutputImage&& out);
+
+
+  /******************************************/
+  /****          Implementation          ****/
+  /******************************************/
+
+
+  template <class InputImage, class SE, class BorderManager, class OutputImage>
+  void median_filter(InputImage&& input, const mln::experimental::StructuringElement<SE>& se, BorderManager bm,
+                     OutputImage&& out)
+  {
+    using R = std::ratio<1, 2>;
+    rank_filter<R>(std::forward<InputImage>(input), se, bm, out);
+  }
+
+
+  template <class InputImage, class SE, class BorderManager>
+  image_concrete_t<std::remove_reference_t<InputImage>>
+  median_filter(InputImage&& image, const mln::experimental::StructuringElement<SE>& se, BorderManager bm)
+  {
+    using R = std::ratio<1, 2>;
+    return rank_filter<R>(std::forward<InputImage>(image), se, bm);
+  }
+
+} // namespace mln

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

+#pragma once
+
+
+#include <mln/accu/accumulators/h_rank.hpp>
+#include <mln/core/canvas/local_accumulation.hpp>
+#include <mln/core/concept/new/images.hpp>
+#include <mln/core/concept/new/structuring_elements.hpp>
+#include <mln/core/extension/border_management.hpp>
+#include <mln/core/trace.hpp>
+
+
+namespace mln::morpho::experimental
+{
+
+  /// \brief Rank filter
+  ///
+  /// A rank filter is non-linerar filter using the value ordering to
+  /// compute the filtered value. It computes:
+  /// \f[
+  /// g(x) = r(\{ f(y) \in \mathcal{B}_x) \})
+  /// \f]
+  /// where \p r returns the 𝑟-th value of the set of pixels of the
+  /// structuring element 𝑩 centered in 𝑥.
+  ///
+  ///
+  /// \param ima Input image 𝑓
+  /// \param se  Structuring element
+  /// \param out (optional) Output image 𝑔 so store the result.
+  /// \tparam Ratio A std::ratio type that defines the rank of the filtered value.
+  template <class Ratio, class InputImage, class SE, class BorderManager>
+  image_concrete_t<std::remove_reference_t<InputImage>>
+  rank_filter(InputImage&& image, const mln::experimental::StructuringElement<SE>& se, BorderManager bm);
+
+  template <class Ratio, class InputImage, class SE, class BorderManager, class OutputImage>
+  void rank_filter(InputImage&& image, const mln::experimental::StructuringElement<SE>& se, BorderManager bm,
+                   OutputImage&& out);
+
+
+  /******************************************/
+  /****          Implementation          ****/
+  /******************************************/
+
+
+  template <class Ratio, class InputImage, class SE, class BorderManager, class OutputImage>
+  void rank_filter(InputImage&& input, const mln::experimental::StructuringElement<SE>& se, BorderManager bm,
+                   OutputImage&& out)
+  {
+    using I = std::remove_reference_t<InputImage>;
+    static_assert(mln::is_a<I, mln::experimental::Image>::value);
+
+    using V = image_value_t<I>;
+    mln_entering("mln::morpho::experimental::rank_filter");
+
+    // To enable when we can concept check that domain are comparable
+    // assert(image.domain() == out.domain());
+
+    auto [f, ses] = bm.manage(input, static_cast<const SE&>(se));
+
+    std::visit(
+        [&out](auto&& f, auto&& se) {
+          mln::accu::accumulators::h_rank<V, Ratio> accu;
+          mln::canvas::LocalAccumulation            algo(accu, se, f, out);
+          algo.Execute();
+        },
+        f, ses);
+  }
+
+
+  template <class Ratio, class InputImage, class SE, class BorderManager>
+  image_concrete_t<std::remove_reference_t<InputImage>>
+  rank_filter(InputImage&& image, const mln::experimental::StructuringElement<SE>& se, BorderManager bm)
+  {
+    using I = std::remove_reference_t<InputImage>;
+    static_assert(mln::is_a<I, mln::experimental::Image>::value);
+
+    image_concrete_t<I> out = imconcretize(image);
+    rank_filter<Ratio>(image, se, bm, out);
+    return out;
+  }
+
+} // namespace mln::morpho::experimental

tests/morpho/CMakeLists.txt

@@ -13,6 +13,7 @@ add_core_test(${test_prefix}gradient         gradient.cpp)
 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}rank_filter      rank_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)

tests/morpho/median_filter.cpp

-#include <mln/core/algorithm/fill.hpp>
-#include <mln/core/grays.hpp>
-#include <mln/core/image/image2d.hpp>
+#include <mln/morpho/experimental/median_filter.hpp>
+
+#include <mln/core/image/experimental/ndimage.hpp>
 #include <mln/core/se/rect2d.hpp>
-#include <mln/io/imread.hpp>
-#include <mln/morpho/median_filter.hpp>
+#include <mln/core/rangev3/view/zip.hpp>
+
+#include <mln/io/experimental/imread.hpp>
+#include <mln/io/experimental/imsave.hpp>
 
 #include <fixtures/ImageCompare/image_compare.hpp>
 #include <fixtures/ImagePath/image_path.hpp>
@@ -13,42 +15,38 @@
 
 #include <gtest/gtest.h>
 
-using namespace mln;
 
-image2d<uint8> naive_median(const image2d<uint8>& f, se::rect2d win, int sz)
+mln::experimental::image2d<uint8_t> naive_median(const mln::experimental::image2d<uint8_t>& f,
+                                                 mln::experimental::se::rect2d win, int sz)
 {
   mln_entering("naive_median");
 
-  image2d<uint8> g;
-  resize(g, f);
-
-  mln_pixter(px, qx, f, g);
-  mln_iter(nx, win(px));
+  mln::experimental::image2d<uint8_t> g;
+  mln::resize(g, f);
 
-  std::vector<uint8> V;
-
-  mln_forall (px, qx)
+  std::vector<uint8_t> V;
+  auto zz = mln::ranges::view::zip(f.new_pixels(), g.new_values());
+  for (auto&& r : zz.rows())
+    for (auto&& [px, vout] : r)
     {
       V.clear();
-    mln_forall (nx)
-      V.push_back(nx->val());
-    std::sort(V.begin(), V.end());
-    qx->val() = V[sz / 2];
+      for (auto qx : win(px))
+        V.push_back(qx.val());
+      std::partial_sort(V.begin(), V.begin() + sz / 2 + 1, V.end());
+      vout = V[sz / 2];
     }
-
-  mln_exiting();
   return g;
 }
 
-TEST(Morpho, median_filter_median_filter_0)
+TEST(Morpho, median_filter)
 {
-  image2d<uint8> ima;
-  io::imread(fixtures::ImagePath::concat_with_filename("lena.pgm"), ima);
+  mln::experimental::image2d<uint8_t> ima;
+  mln::io::experimental::imread(fixtures::ImagePath::concat_with_filename("lena.pgm"), ima);
 
   { // Fast: border wide enough
-    mln::se::rect2d win(7, 7);
-    auto            out  = morpho::median_filter(ima, win);
+    mln::experimental::se::rect2d win(7, 7);
+    auto                          out  = mln::morpho::experimental::median_filter(ima, win, mln::extension::bm::native::mirror());
     auto                          out2 = naive_median(ima, win, 49);
-    ASSERT_IMAGES_EQ(out2, out);
+    ASSERT_IMAGES_EQ_EXP(out2, out);
   }
 }

tests/morpho/rank_filter.cpp

+#include <mln/morpho/experimental/rank_filter.hpp>
+
+
+#include <mln/core/image/experimental/ndimage.hpp>
+#include <mln/core/se/periodic_line2d.hpp>
+#include <mln/core/algorithm/iota.hpp>
+
+
+#include <fixtures/ImageCompare/image_compare.hpp>
+
+#include <gtest/gtest.h>
+#include <ratio>
+
+TEST(Morpho, rank_filter_fill_back)
+{
+  mln::experimental::image2d<uint8_t> ima(10,5);
+  mln::iota(ima, 0);
+
+  mln::experimental::image2d<uint8_t> ref = {
+      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},           //
+      {0, 1, 2, 3, 4, 5, 6, 7, 8, 9},           //
+      {10, 11, 12, 13, 14, 15, 16, 17, 18, 19}, //
+      {10, 11, 12, 13, 14, 15, 16, 17, 18, 19}, //
+      {0, 0, 0, 0, 0, 0, 0, 0, 0, 0},           //
+  };
+
+  using R = std::ratio<1,5>;
+
+  auto vline = mln::experimental::se::periodic_line2d({0, 1}, 2);
+  auto out =
+    mln::morpho::experimental::rank_filter<R>(ima, vline, mln::extension::bm::fill(uint8_t(0)));
+
+  ASSERT_IMAGES_EQ_EXP(out, ref);
+}
+
+
+TEST(Morpho, rank_filter_fill_mirror)
+{
+  mln::experimental::image2d<uint8_t> ima(10,5);
+  mln::iota(ima, 0);
+
+  mln::experimental::image2d<uint8_t> ref = {
+      {10, 11, 12, 13, 14, 15, 16, 17, 18, 19}, //
+      {10, 11, 12, 13, 14, 15, 16, 17, 18, 19}, //
+      {20, 21, 22, 23, 24, 25, 26, 27, 28, 29}, //
+      {30, 31, 32, 33, 34, 35, 36, 37, 38, 39}, //
+      {30, 31, 32, 33, 34, 35, 36, 37, 38, 39}, //
+  };
+
+  using R = std::ratio<2,5>;
+
+  auto vline = mln::experimental::se::periodic_line2d({0, 1}, 2);
+  auto out   = mln::morpho::experimental::rank_filter<R>(ima, vline, mln::extension::bm::mirror{});
+
+  ASSERT_IMAGES_EQ_EXP(out, ref);
+}
+