Merge branch 'development/watershed-hierarchy-saliency-integration' into 'next'

Integration of watershed hierarchies and saliency map

See merge request !118

AUTHORS

@@ -10,7 +10,8 @@ Contributions
 
 Célian Gossec
 Virgile Hirtz
-
+Quentin Kaci
+Victor Simonin
 
 
 Main past authors from Milena having influenced the library

bench/BMAlphaTree.cpp

+#include <mln/io/imread.hpp>
+
+#include <benchmark/benchmark.h>
+#include <mln/core/colors.hpp>
+#include <mln/core/image/ndimage.hpp>
+#include <mln/core/neighborhood/c4.hpp>
+#include <mln/morpho/alphatree.hpp>
+
+using image_t = mln::image2d<mln::rgb8>;
+
+static const std::vector<std::string> bench_images = {"Aerial_view_of_Olbia.jpg", "Space1_20MB.jpg"};
+
+static const auto alphatree_vect_function = [](const image_t& input) {
+  mln::morpho::alphatree(input, mln::c4, [](const auto& a, const auto& b) -> double { return mln::l2dist(a, b); });
+};
+
+static const auto alphatree_hq_function = [](const image_t& input) {
+  mln::morpho::alphatree(input, mln::c4,
+                         [](const auto& a, const auto& b) -> std::uint16_t { return mln::l2dist(a, b); });
+};
+
+class BMAlphaTree : public benchmark::Fixture
+{
+public:
+  BMAlphaTree()
+  {
+    if (inputs.empty())
+    {
+      inputs = std::vector<image_t>(bench_images.size());
+      for (std::size_t i = 0; i < bench_images.size(); ++i)
+        mln::io::imread(bench_images[i], inputs[i]);
+    }
+  }
+
+  void run(benchmark::State& st, const std::function<void(const image_t& input)>& callback, std::size_t i) const
+  {
+    for (auto _ : st)
+      callback(inputs[i]);
+  }
+
+  std::vector<image_t> inputs;
+};
+
+
+BENCHMARK_F(BMAlphaTree, AlphatreeOlbiaVect)(benchmark::State& st)
+{
+  this->run(st, alphatree_vect_function, 0);
+}
+
+BENCHMARK_F(BMAlphaTree, AlphatreeOlbiaHQ)(benchmark::State& st)
+{
+  this->run(st, alphatree_hq_function, 0);
+}
+
+BENCHMARK_F(BMAlphaTree, AlphatreeSpaceVect)(benchmark::State& st)
+{
+  this->run(st, alphatree_vect_function, 1);
+}
+
+BENCHMARK_F(BMAlphaTree, AlphatreeSpaceHQ)(benchmark::State& st)
+{
+  this->run(st, alphatree_hq_function, 1);
+}
+
+BENCHMARK_MAIN();
\ No newline at end of file

bench/BMWatershedHierarchy.cpp

+#include <benchmark/benchmark.h>
+#include <mln/accu/accumulators/count.hpp>
+#include <mln/accu/accumulators/mean.hpp>
+#include <mln/core/colors.hpp>
+#include <mln/core/neighborhood/c4.hpp>
+#include <mln/io/imread.hpp>
+#include <mln/morpho/watershed_hierarchy.hpp>
+
+using image_t = mln::image2d<mln::rgb8>;
+
+static const std::vector<std::string> bench_images = {"Aerial_view_of_Olbia.jpg", "Space1_20MB.jpg"};
+
+static const auto watershed_hierarchy_by_area_function = [](const image_t& input) {
+  mln::morpho::watershed_hierarchy(
+      input,
+      [](auto tree, auto nm) -> std::vector<size_t> {
+        return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
+      },
+      mln::c4, [](const auto& a, const auto& b) -> std::float_t { return mln::l2dist(a, b); });
+};
+
+static const auto hierarchical_segmentation_function = [](const image_t& input) {
+  auto [tree, node_map] = mln::morpho::watershed_hierarchy(
+      input,
+      [](auto tree, auto nm) -> std::vector<size_t> {
+        return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
+      },
+      mln::c4, [](const auto& a, const auto& b) -> std::float_t { return mln::l2dist(a, b); });
+
+  auto mean = tree.compute_attribute_on_values(node_map, input, mln::accu::accumulators::mean<mln::rgb8>());
+
+  // Threshold cut at 10
+  auto node_map_cut = tree.horizontal_cut(10, node_map);
+  tree.reconstruct_from(node_map_cut, ranges::make_span(mean));
+};
+
+class BMWatershedHierarchy : public benchmark::Fixture
+{
+public:
+  BMWatershedHierarchy()
+  {
+    if (inputs.empty())
+    {
+      inputs = std::vector<image_t>(bench_images.size());
+      for (std::size_t i = 0; i < bench_images.size(); ++i)
+        mln::io::imread(bench_images[i], inputs[i]);
+    }
+  }
+
+  void run(benchmark::State& st, const std::function<void(const image_t& input)>& callback, std::size_t i) const
+  {
+    for (auto _ : st)
+      callback(inputs[i]);
+  }
+
+  static std::vector<image_t> inputs;
+};
+std::vector<image_t> BMWatershedHierarchy::inputs;
+
+BENCHMARK_F(BMWatershedHierarchy, WatershedHierarchyAreaOlbia)(benchmark::State& st)
+{
+  this->run(st, watershed_hierarchy_by_area_function, 0);
+}
+
+BENCHMARK_F(BMWatershedHierarchy, HierarchicalSegmentationOlbia)(benchmark::State& st)
+{
+  this->run(st, hierarchical_segmentation_function, 0);
+}
+
+BENCHMARK_F(BMWatershedHierarchy, WatershedHierarchyAreaSpace)(benchmark::State& st)
+{
+  this->run(st, watershed_hierarchy_by_area_function, 1);
+}
+
+BENCHMARK_F(BMWatershedHierarchy, HierarchicalSegmentationSpace)(benchmark::State& st)
+{
+  this->run(st, hierarchical_segmentation_function, 1);
+}
+
+
+BENCHMARK_MAIN();
\ No newline at end of file

bench/CMakeLists.txt

@@ -64,5 +64,7 @@ add_benchmark(BMMorphers                BMMorphers.cpp BMMorphers_main.cpp)
 add_benchmark(BMReference_Linear        BMReference_Linear.cpp BMReference_Linear_Reversed.cpp BMReference_Linear_main.cpp)
 add_benchmark(BMReference_Neighborhood  BMReference_Neighborhood_main.cpp)
 add_benchmark(BMBufferPrimitives        BMBufferPrimitives.cpp)
+add_benchmark(BMAlphaTree               BMAlphaTree.cpp)
+add_benchmark(BMWatershedHierarchy      BMWatershedHierarchy.cpp)
 
 ExternalData_Add_Target(fetch-external-data)

doc/source/morpho.rst

@@ -48,7 +48,7 @@ Component Trees & Hierarchical Representations
     morpho/maxtree
     morpho/tos
     morpho/alphatree
-
+    morpho/watershed_hierarchy
 
 
 

doc/source/morpho/component_tree.rst

@@ -294,6 +294,29 @@ to make an horizontal cut of this tree.
     :param nodemap: An image thats maps ``point -> node id``
     :param levels: (Optional) The altitude of each node in the tree (for example the :math:`\alpha` associated to each node for the alphatree).
 
+Saliency Computation
+--------------
+It is also possible to compute the saliency map to obtain another visualization.
+
+.. cpp:function:: auto saliency(Image node_map)
+
+    Compute and return the saliency map of the tree.
+
+    :param node_map: An image thats maps ``point -> node id``
+
+    :return: The saliency map as an image
+
+.. list-table::
+
+   * - .. image:: /images/watershed_hierarchy_area_gray.png
+          :width: 100%
+
+     - .. image:: /images/saliency_watershed.png
+          :width: 100%
+
+   * - Watershed hierarchy by area with a cut at a threshold of 25
+     - The corresponding saliency map
+
 A complete example
 ------------------
 

doc/source/morpho/watershed_hierarchy.rst

+Hierarchical Watershed
+==========
+
+* Include :file:`<mln/morpho/watershed_hierarchy.hpp>`
+
+.. cpp:namespace:: mln::morpho
+
+
+.. cpp:function:: auto watershed_hierarchy(Image input, AttributeFunction attribute_func, Neighborhood nbh, F dist);
+
+    Compute the watershed hierarchy and returns a pair `(tree, node_map)`.
+    See :doc:`component_tree` for more information about the representation of tree.
+    The implementation is based on [Naj13]_.
+
+    :param input: The input image
+    :param attribute_func: The function that define the attribute computation
+    :param nbh: The neighborhood relation
+    :param dist: The function weighting the edges between two pixels.
+    :return: A pair `(tree, node_map)` where *tree* is of type ``component_tree<std::invoke_result_t<F, image_value_t<Image>, image_value_t<Image>>>`` and *node_map* is a mapping between the image pixels and the node of the tree.
+
+    .. rubric:: Requirements
+
+    * ``std::invoke_result_t<F, image_value_t<Image>, image_value_t<Image>>`` is :cpp:concept:`std::totally_ordered`
+
+Definition
+----------
+
+Watershed hierarchies, also called attribute-based hierarchies, are constructed by filtering and re-weighting the
+Minimum Spanning Tree based on an attribute computation (area, volume, dynamic, etc.) on the Binary Partition Tree. This
+process can be imagined as the filling of branches of the Binary Partition Tree with labeled water where the attribute
+drives the speed of the filling.
+
+The tree that represent a watershed hierarchy is an Alpha tree, also known as quasi-flat zone hierarchy.
+See :doc:`alphatree` for more information.
+
+Example
+-------
+
+This example is used to generate the grayscale lena watershed hierarchy by area with a cut at a threshold of 25 below.
+
+    ::
+
+        #include <mln/accu/accumulators/count.hpp>
+        #include <mln/accu/accumulators/mean.hpp>
+        #include <mln/core/image/ndimage.hpp>
+        #include <mln/core/neighborhood/c4.hpp>
+        #include <mln/morpho/watershed_hierarchy.hpp>
+
+        mln::image2d<uint8_t> input = ...;
+
+        // Compute the watershed hierarchy by area
+        auto area_attribute_func = [](auto tree, auto node_map) -> std::vector<size_t> {
+            return tree.compute_attribute_on_points(node_map, mln::accu::features::count<>());
+        };
+        auto [tree, node_map] = mln::morpho::watershed_hierarchy(input, area_attribute_func, mln::c4);
+
+        // Compute an attribute (for example the average pixels value at each node, as below)
+        auto mean = tree.compute_attribute_on_values(node_map, input, mln::accu::accumulators::mean<uint8_t>());
+
+        // Making an horizontal cut of the tree
+        const auto threshold = 25; // Threshold of the horizontal cut, that means the lowest alpha in the cut
+        auto node_map_cut = tree.horizontal_cut(threshold, node_map); // Return a new node map associated to the cut
+
+        // Labeling the cut with the mean values of each node
+        auto out = tree.reconstruct_from(node_map_cut, ranges::make_span(mean)); // Using range-v3 span
+
+    .. list-table::
+
+        * -   .. figure:: /images/watershed_hierarchy_area_color.png
+
+                Watershed hierarchy by area with a cut at a threshold of 100
+
+          -   .. figure:: /images/watershed_hierarchy_area_gray.png
+
+                Watershed hierarchy by area with a cut at a threshold of 25
+
+Notes
+-----
+
+Complexity
+----------
+
+References
+----------
+
+.. [Naj13] Laurent Najman, Jean Cousty, and Benjamin Perret (2013). Playing with kruskal: algorithms for morphological trees in edge-weighted graphs. *International Symposium on Mathematical Morphology and Its Applications to Signal and Image Processing*. Springer, Berlin, Heidelberg. 135-146
\ No newline at end of file

doc/source/snippets/CMakeLists.txt

@@ -54,12 +54,16 @@ add_image("area_filter"
   morpho_area_filter_out.png)
 
 
-
 add_image("blobs_watershed" "${DOCUMENTATION_IMAGE_DIR}/blobs_binary.png" blobs_distance_transform.png blobs_segmentation.png)
 
 add_image("alphatree_example" "${PYLENE_IMAGE_DIR}/lena.ppm" alphatree_cut_color.png)
 add_image("alphatree_example" "${PYLENE_IMAGE_DIR}/lena.pgm" alphatree_cut_gray.png)
 
+add_image("watershed_hierarchy_example" "${PYLENE_IMAGE_DIR}/lena.ppm" watershed_hierarchy_area_color.png)
+add_image("watershed_hierarchy_example" "${PYLENE_IMAGE_DIR}/lena.pgm" watershed_hierarchy_area_gray.png)
+
+add_image("saliency_example" "${PYLENE_IMAGE_DIR}/lena.pgm" saliency_watershed.png)
+
 add_custom_target(build-images
                   DEPENDS "${DOCUMENTATION_IMAGES}")
 
@@ -73,6 +77,8 @@ link_libraries(Pylene::Core)
 link_libraries(doc-lib)
 
 add_executable(alphatree_example alphatree_example.cpp)
+add_executable(watershed_hierarchy_example watershed_hierarchy_example.cpp)
+add_executable(saliency_example saliency_example.cpp)
 add_executable(erosion-cli erosion-cli.cpp)
 add_executable(staff_lines staff_lines.cpp)
 add_executable(component_tree_1 component_tree_1.cpp)
@@ -90,4 +96,4 @@ if (TARGET Boost::program_options)
   target_link_libraries(erosion-cli PRIVATE Boost::program_options)
 elseif (TARGET Boost::Boost)
   target_link_libraries(erosion-cli PRIVATE Boost::Boost)
-endif()
+endif()
\ No newline at end of file

doc/source/snippets/saliency_example.cpp

+#include <mln/accu/accumulators/mean.hpp>
+#include <mln/core/colors.hpp>
+#include <mln/core/image/ndimage.hpp>
+#include <mln/core/image/view/cast.hpp>
+#include <mln/core/neighborhood/c4.hpp>
+#include <mln/io/imread.hpp>
+#include <mln/io/imsave.hpp>
+#include <mln/morpho/tos.hpp>
+#include <mln/morpho/watershed_hierarchy.hpp>
+
+
+#include <iostream>
+
+void process_example(const mln::image2d<uint8_t>& img, const std::string& output_filename, const double threshold)
+{
+  // 2. Build the watershed hierarchy
+  auto area_attribute_func = [](auto tree, auto nm) -> std::vector<size_t> {
+    return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
+  };
+  auto [t, nm] = mln::morpho::watershed_hierarchy(img, area_attribute_func, mln::c4);
+
+  // 3. Compute the mean attribute
+  auto mean = t.compute_attribute_on_values(nm, img, mln::accu::accumulators::mean<uint8_t>());
+
+  // 4. Compute a cut of the watershed hierarchy
+  auto cut_nm = t.horizontal_cut(threshold, nm);
+
+  // 5. Label the cut
+  auto cut = t.reconstruct_from(cut_nm, ::ranges::make_span(mean));
+
+  mln::io::imsave(mln::view::cast<uint8_t>(cut), output_filename);
+
+  mln::image2d<uint8_t> in;
+  mln::io::imread(output_filename, in);
+
+  auto [t2, node_map] = mln::morpho::tos(in, {0, 0});
+
+  auto saliency = t2.saliency(in);
+
+  // 5. Save the output
+  mln::io::imsave(mln::view::cast<uint8_t>(saliency), output_filename);
+}
+
+int main(int argc, char* argv[])
+{
+  if (argc < 3)
+  {
+    std::cerr << "Invalid number of argument\nUsage: " << argv[0] << " input_filename output_filename\n";
+    return 1;
+  }
+  auto in_filename  = std::string(argv[1]);
+  auto out_filename = std::string(argv[2]);
+
+  // 1. Read the input image
+  mln::ndbuffer_image in;
+  in = mln::io::imread(in_filename);
+
+  if (in.sample_type() == mln::sample_type_id::UINT8)
+  {
+    const auto* img = in.cast_to<std::uint8_t, 2>();
+    process_example(*img, out_filename, 25);
+  }
+  else
+  {
+    std::cerr << "Unhandled sample type format\n";
+    return 1;
+  }
+
+  return 0;
+}
\ No newline at end of file

doc/source/snippets/watershed_hierarchy_example.cpp

+#include <mln/accu/accumulators/mean.hpp>
+#include <mln/core/colors.hpp>
+#include <mln/core/image/ndimage.hpp>
+#include <mln/core/image/view/cast.hpp>
+#include <mln/core/neighborhood/c4.hpp>
+#include <mln/io/imread.hpp>
+#include <mln/io/imsave.hpp>
+#include <mln/morpho/watershed_hierarchy.hpp>
+
+#include <iostream>
+
+template <typename V>
+void process_example(const mln::image2d<V>& img, const std::string& output_filename, const double threshold)
+{
+  // 2. Build the watershed hierarchy
+  auto area_attribute_func = [](auto tree, auto nm) -> std::vector<size_t> {
+    return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
+  };
+  auto [t, nm] = mln::morpho::watershed_hierarchy(img, area_attribute_func, mln::c4);
+
+  // 3. Compute the mean attribute
+  auto mean = t.compute_attribute_on_values(nm, img, mln::accu::accumulators::mean<V>());
+
+  // 4. Compute a cut of the watershed hierarchy
+  auto cut_nm = t.horizontal_cut(threshold, nm);
+
+  // 5. Label the cut
+  auto cut = t.reconstruct_from(cut_nm, ::ranges::make_span(mean));
+
+  // 5. Save the output
+  mln::io::imsave(mln::view::cast<V>(cut), output_filename);
+}
+
+int main(int argc, char* argv[])
+{
+  if (argc < 3)
+  {
+    std::cerr << "Invalid number of argument\nUsage: " << argv[0] << " input_filename output_filename\n";
+    return 1;
+  }
+  auto in_filename  = std::string(argv[1]);
+  auto out_filename = std::string(argv[2]);
+
+  // 1. Read the input image
+  mln::ndbuffer_image in;
+  in = mln::io::imread(in_filename);
+
+  if (in.sample_type() == mln::sample_type_id::UINT8)
+  {
+    const auto* img = in.cast_to<std::uint8_t, 2>();
+    process_example(*img, out_filename, 25);
+  }
+  else if (in.sample_type() == mln::sample_type_id::RGB8)
+  {
+    const auto* img = in.cast_to<mln::rgb8, 2>();
+    process_example(*img, out_filename, 100);
+  }
+  else
+  {
+    std::cerr << "Unhandled sample type format\n";
+    return 1;
+  }
+
+  return 0;
+}
\ No newline at end of file

pylene/include/mln/morpho/alphatree.hpp

@@ -28,7 +28,7 @@ namespace mln::morpho
   template <class I, class N, class F = mln::functional::l2dist_t<>>
   std::pair<component_tree<std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>, image_ch_value_t<I, int>> //
   alphatree(I input, N nbh, F distance = F{});
-
+  
 
   /******************************************/
   /****          Implementation          ****/
@@ -39,11 +39,11 @@ namespace mln::morpho
   {
     /// \brief Canvas for the edges in the alphatree. Using different data
     /// structures related to the type of the edges.
-    template <typename P, typename N, typename W>
+    template <typename P, typename N, typename W, bool HQ>
     class alphatree_edges;
 
-    template <typename P, typename N, typename W>
-    requires(std::is_integral_v<W>&& std::is_unsigned_v<W> && sizeof(W) <= 2) class alphatree_edges<P, N, W>
+    template <typename P, typename N, typename W, bool HQ>
+    requires(std::is_integral_v<W>&& std::is_unsigned_v<W> && sizeof(W) <= 2 && HQ) class alphatree_edges<P, N, W, HQ>
     {
     public:
       void                push(int dir, W w, P p) { m_cont.insert(dir, w, p); }
@@ -64,11 +64,14 @@ namespace mln::morpho
       W w;
     };
 
-    template <typename P, typename N, typename W>
+    template <typename P, typename N, typename W, bool HQ>
     class alphatree_edges
     {
     public:
-      void                push(int dir, W w, P p) { m_cont.push_back({p, p + cn.after_offsets()[dir], w}); }
+      void push(int dir, W w, P p) { m_cont.push_back({p, p + cn.after_offsets()[dir], w}); }
+
+      void push(P p, P q, W w) { m_cont.push_back({p, q, w}); }
+
       std::tuple<P, P, W> pop()
       {
         assert(m_current < m_cont.size());
@@ -81,10 +84,13 @@ namespace mln::morpho
         assert(m_current < m_cont.size());
         return m_cont[m_current].w;
       }
+
       bool empty() const { return m_cont.size() == m_current; }
+
       void on_finish_insert()
       {
-        std::sort(m_cont.begin(), m_cont.end(), [](const edge_t<P, W>& a, const edge_t<P, W>& b) { return a.w < b.w; });
+        std::stable_sort(m_cont.begin(), m_cont.end(),
+                         [](const edge_t<P, W>& a, const edge_t<P, W>& b) { return a.w < b.w; });
       }
 
     private:
@@ -115,8 +121,6 @@ namespace mln::morpho
     template <class E, class J>
     void alphatree_compute_flatzones(E& edges, J zpar)
     {
-      canvas::impl::union_find_init_par(zpar);
-
       while (!edges.empty() && edges.top() == 0)
       {
         const auto [p, q, w] = edges.pop();
@@ -174,7 +178,8 @@ namespace mln::morpho
                                             std::size_t       node_count, //
                                             std::vector<int>& par,        //
                                             std::vector<W>&   levels,     //
-                                            std::vector<M>*   mst)
+                                            std::vector<M>*   mst,        //
+                                            bool              canonize_tree)
     {
       static_assert(mln::is_a<I, mln::details::Image>());
 
@@ -205,11 +210,11 @@ namespace mln::morpho
           int min_root = std::min(rp_root, rq_root);
 
           int new_root_id;
-          if (levels[max_root] == w)
+          if (canonize_tree && levels[max_root] == w)
           {
             new_root_id = max_root;
           }
-          else if (levels[min_root] == w)
+          else if (canonize_tree && levels[min_root] == w)
           {
             new_root_id = min_root;
           }
@@ -233,14 +238,18 @@ namespace mln::morpho
     }
 
     template <class W>
-    std::pair<std::vector<int>, std::vector<W>> canonize_component_tree(const std::vector<int>& par,    //
-                                                                        const std::vector<W>&   levels, //
-                                                                        std::size_t             node_count)
+    std::pair<std::vector<int>, std::vector<W>> canonize_component_tree(const std::vector<int>& par,        //
+                                                                        const std::vector<W>&   levels,     //
+                                                                        std::size_t             node_count, //
+                                                                        std::size_t             nb_leaves)
     {
+      assert(node_count >= nb_leaves);
+
       std::vector<int> canonized_par;
       std::vector<W>   canonized_levels;
 
       // Root initialization
+
       canonized_par.push_back(0);
       canonized_levels.push_back(levels[0]);
 
@@ -249,8 +258,11 @@ namespace mln::morpho
       translation_map[0] = 0;
       int count          = 1;
 
+      std::size_t begin_leaves = node_count - nb_leaves;
+
       // Build canonized component tree
-      for (std::size_t i = 1; i < node_count; ++i)
+
+      for (std::size_t i = 1; i < begin_leaves; ++i)
       {
         if (levels[i] != levels[par[i]]) // Keep the node: Update tree
         {
@@ -262,13 +274,53 @@ namespace mln::morpho
           translation_map[i] = translation_map[par[i]];
       }
 
+      for (std::size_t i = begin_leaves; i < node_count; ++i)
+      {
+        translation_map[i] = count++;
+        canonized_par.push_back(translation_map[par[i]]);
+        canonized_levels.push_back(levels[i]);
+      }
+
       return {canonized_par, canonized_levels};
     }
 
-    template <class I, class N, class F,