Implement area connected filter

doc/source/morpho.rst

@@ -15,7 +15,15 @@ Structural Morphological Operation
    morpho/rank_filter
    morpho/median_filter
    morpho/gradient
+
+Geodesic transformations
+************************
+
+.. toctree::
+   :maxdepth: 1
+
    morpho/opening_by_reconstruction
+   morpho/area_filter
 
 
 Segmentation

doc/source/morpho/area_filter.rst

+Opening and Closing by Area
+===========================
+
+Include :file:`<mln/morpho/area_filter.hpp>`
+
+.. cpp:namespace:: mln::morpho
+
+.. cpp:function:: \
+    Image{I} concrete_t<I> area_opening(I f, Neighborhood nbh, int area, Compare cmp)
+    Image{I} concrete_t<I> area_closing(I f, Neighborhood nbh, int area)
+
+
+    On binary images, the area connected opening that preserves connected
+    components (blobs) of a minimum size. On grayscale images, it extracts
+    connected image objects of higher intensity values than the surrounding
+    objects. (The area closing is its complementary operator).
+
+
+   :param f: Input image 𝑓
+   :param nbh:  Elementary structuring element.
+   :param area: The minimal size of the blobs in order to be preserved
+
+   :return: An image whose type is deduced from the input image
+
+   :exception: N/A
+
+
+Notes
+-----
+
+
+Complexity
+----------
+
+
+Example: dense objects detection
+--------------------------------
+
+.. list-table::
+
+   *  -  .. figure:: /images/blobs2_binary.png
+
+            (a) Original image
+
+
+      -  .. figure:: /images/morpho_area_filter_dilated.png
+
+            (b) Dilated of the original image (a)
+
+   *  -  .. figure:: /images/morpho_area_filter_dilated.png
+
+            (c) Result of the area opening of (b)
+
+
+      -  .. figure:: /images/morpho_area_filter_out.png
+
+            (d) Input image masked by (c)
+
+
+.. literalinclude:: /snippets/area_filter.cpp
+   :start-after: M2_START
+   :end-before: M2_END
+   :language: cpp
+
+Given an original image. A dilation with a small disc allows to connect
+objects and we remove small connected components with an area opening. 
+Finally, we just have to mask the input with the mask to get the objects in dense regions.
+

doc/source/snippets/CMakeLists.txt

@@ -45,6 +45,14 @@ add_image("reconstruction"
   morpho_reconstruction_markers.png
   morpho_reconstruction_rec.png
   morpho_reconstruction_out.png)
+add_image("area_filter"
+  "${DOCUMENTATION_IMAGE_DIR}/blobs2_binary.png"
+  morpho_area_filter_dilated.png
+  morpho_area_filter_opening.png
+  morpho_area_filter_out.png)
+
+
+
 add_image("blobs_watershed" "${DOCUMENTATION_IMAGE_DIR}/blobs_binary.png" blobs_distance_transform.png blobs_segmentation.png)
 
 
@@ -59,4 +67,6 @@ 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)
+add_executable(area_filter area_filter.cpp)
 add_executable(cdt cdt.cpp)
+ 

doc/source/snippets/area_filter.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/morpho/experimental/dilation.hpp>
+#include <mln/morpho/experimental/area_filter.hpp>
+#include <mln/core/image/view/operators.hpp>
+#include <ratio>
+
+int main(int argc, char** argv)
+{
+  if (argc < 5)
+  {
+    std::cerr << "Usage: " << argv[0] << " input dilated mask output\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(2);
+  auto dil = mln::morpho::experimental::dilation(input, disc);
+
+  // Filtering
+  auto mask = mln::morpho::experimental::area_opening(dil, mln::experimental::c8, 2000);
+
+  // Mask
+  auto out = mask && input;
+   // #M2_END
+
+  // Save
+  mln::io::experimental::imsave(dil, argv[2]);
+  mln::io::experimental::imsave(mask, argv[3]);
+  mln::io::experimental::imsave(out, argv[4]);
+}

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

+#pragma once
+
+#include <mln/core/concept/new/images.hpp>
+#include <mln/core/concept/new/neighborhoods.hpp>
+
+#include <mln/morpho/experimental/canvas/unionfind.hpp>
+
+
+
+namespace mln::morpho::experimental
+{
+
+  /// \brief Compute the area algebraic opening of an image.
+  /// \param ima The input image
+  /// \param nbh The neighborhood
+  /// \param area The grain size
+  /// \param cmp A strict total ordering on values
+  template <class I, class N, class Compare = std::less<image_value_t<std::remove_reference_t<I>>>>
+  image_concrete_t<std::remove_reference_t<I>> //
+  area_opening(I&& ima, const N& nbh, int area, Compare cmp = {});
+
+  /// \brief Compute the area algebraic closing of an image.
+  /// \param ima The input image
+  /// \param nbh The neighborhood
+  /// \param area The grain size
+  template <class I, class N>
+  image_concrete_t<std::remove_reference_t<I>> //
+  area_closing(I&& ima, const N& nbh, int area);
+
+  /******************************/
+  /*** Implementation          **/
+  /******************************/
+
+  namespace impl
+  {
+    template <class I>
+    struct area_filter_ufind_visitor
+    {
+      using P = image_point_t<I>;
+
+      void on_make_set(P p) noexcept { m_count(p) = 1; }
+
+      void on_union(P p, P q) noexcept { m_count(q) += m_count(p); }
+
+      void on_finish(P p, P root) noexcept { m_ima(p) = m_ima(root); }
+
+      bool test(P p) const noexcept { return m_count(p) >= m_area; }
+
+      area_filter_ufind_visitor(I& input, int area)
+        : m_ima(input)
+        , m_area(area)
+      {
+        m_count = imchvalue<int>(input);
+      }
+
+
+    private:
+      I&                       m_ima;
+      int                      m_area;
+      image_ch_value_t<I, int> m_count;
+    };
+
+
+    template <class I, class N, class Compare>
+    void area_opening_inplace(I& ima, const N& nbh, int area, Compare cmp)
+    {
+      area_filter_ufind_visitor<I> viz(ima, area);
+      mln::morpho::experimental::canvas::union_find(ima, nbh, viz, cmp);
+    }
+
+  }
+
+
+  template <class InputImage, class N, class Compare>
+  image_concrete_t<std::remove_reference_t<InputImage>> //
+  area_opening(InputImage&& ima, const N& nbh, int area, Compare cmp)
+  {
+    using I = std::remove_reference_t<InputImage>;
+
+    static_assert(mln::is_a<I, mln::experimental::Image>());
+    static_assert(mln::is_a<N, mln::experimental::Neighborhood>());
+
+    mln_entering("mln::morpho::area_opening");
+    image_concrete_t<I> out = clone(ima);
+    impl::area_opening_inplace(out, nbh, area, std::move(cmp));
+
+    return out;
+    }
+
+    template <class InputImage, class N>
+    image_concrete_t<std::remove_reference_t<InputImage>> //
+    area_closing(InputImage&& ima, const N& nbh, int area)
+    {
+      using I = std::remove_reference_t<InputImage>;
+
+      static_assert(mln::is_a<I, mln::experimental::Image>());
+      static_assert(mln::is_a<N, mln::experimental::Neighborhood>());
+
+      mln_entering("mln::morpho::area_closing");
+      image_concrete_t<I> out = clone(ima);
+      impl::area_opening_inplace(out, nbh, area, std::greater<image_value_t<I>>());
+
+      return out;
+    }
+
+} // namespace mln::moprho::experimental

tests/morpho/CMakeLists.txt

@@ -17,4 +17,5 @@ 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}area_filter      area_filter.cpp)
 add_core_test(${test_prefix}ToS              tos.cpp tos_tests_helper.cpp)

tests/morpho/area_filter.cpp

+#include <mln/morpho/experimental/area_filter.hpp>
+
+#include <mln/core/image/experimental/ndimage.hpp>
+#include <mln/core/neighborhood/c4.hpp>
+#include <mln/core/neighborhood/c8.hpp>
+
+#include <fixtures/ImageCompare/image_compare.hpp>
+
+#include <mln/io/experimental/imprint.hpp>
+#include <gtest/gtest.h>
+
+
+using namespace mln;
+
+TEST(Morpho, area_opening_grayscale)
+{
+  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::experimental::image2d<uint8_t> ref = {{+2, +2, +2, +2, +2}, //
+                                                     {30, 30, 30, 30, 30}, //
+                                                     {30, 20, 20, 20, 30}, //
+                                                     {30, 30, 20, 30, 30}, //
+                                                     {+1, +5, 20, +5, +1}};
+
+    auto res = mln::morpho::experimental::area_opening(input, mln::experimental::c4, 5);
+    ASSERT_IMAGES_EQ_EXP(ref, res);
+  }
+
+  {
+    const mln::experimental::image2d<uint8_t> ref = {{+2, +2, +2, +2, +2}, //
+                                                     {20, 20, 20, 20, 20}, //
+                                                     {20, 20, 20, 20, 20}, //
+                                                     {20, 20, 20, 20, 20}, //
+                                                     {+1, +5, 20, +5, +1}};
+
+    auto res = mln::morpho::experimental::area_opening(input, mln::experimental::c4, 12);
+    ASSERT_IMAGES_EQ_EXP(ref, res);
+  }
+}
+
+TEST(Morpho, area_opening_binary)
+{
+  const mln::experimental::image2d<bool> input = {{0, 1, 1, 0, 0, 0, 1}, //
+                                                  {1, 0, 1, 0, 0, 1, 1},
+                                                  {0, 1, 0, 0, 1, 1, 1}};
+
+
+  {
+    const mln::experimental::image2d<bool> ref = {{0, 1, 1, 0, 0, 0, 1}, //
+                                                  {0, 0, 1, 0, 0, 1, 1},
+                                                  {0, 0, 0, 0, 1, 1, 1}};
+
+    auto res = mln::morpho::experimental::area_opening(input, mln::experimental::c4, 3);
+    ASSERT_IMAGES_EQ_EXP(ref, res);
+  }
+
+  {
+    const mln::experimental::image2d<bool> ref = {{0, 0, 0, 0, 0, 0, 1}, //
+                                                  {0, 0, 0, 0, 0, 1, 1},
+                                                  {0, 0, 0, 0, 1, 1, 1}};
+
+
+    auto res = mln::morpho::experimental::area_opening(input, mln::experimental::c4, 5);
+    ASSERT_IMAGES_EQ_EXP(ref, res);
+  }
+
+
+  {
+    const mln::experimental::image2d<bool> ref = {{0, 1, 1, 0, 0, 0, 1}, //
+                                                  {1, 0, 1, 0, 0, 1, 1},
+                                                  {0, 1, 0, 0, 1, 1, 1}};
+
+    auto res = mln::morpho::experimental::area_opening(input, mln::experimental::c8, 5);
+    ASSERT_IMAGES_EQ_EXP(ref, res);
+  }
+
+
+  {
+    const mln::experimental::image2d<bool> ref = {{0, 0, 0, 0, 0, 0, 1}, //
+                                                  {0, 0, 0, 0, 0, 1, 1},
+                                                  {0, 0, 0, 0, 1, 1, 1}};
+
+    auto res = mln::morpho::experimental::area_opening(input, mln::experimental::c8, 6);
+    ASSERT_IMAGES_EQ_EXP(ref, res);
+  }
+}