Merge branch 'development/mtos' into 'next'

Implementation of the Multivariate Tree of Shapes

See merge request !123

doc/source/core/pad.rst

@@ -119,7 +119,7 @@ This file provides functions to pad a rectangular buffer or an image either:
                   :param input: The input ndimage (or ndimensional buffer).
                   :param output: The input ndimage (or ndimensional buffer).
                   :param mode: The padding policy
-                  :param value: The value used to fill the array (used if mode is `PAD_CONSTANT`
+                  :param value: The value used to fill the array (used if mode is `PAD_CONSTANT`)
 
                   **Second version**
 
@@ -127,7 +127,7 @@ This file provides functions to pad a rectangular buffer or an image either:
                   :param out: The output ndimensional buffer.
                   :param roi: The roi of the output buffer
                   :param mode: The padding policy
-                  :param value: The value used to fill the array (used if mode is `PAD_CONSTANT`
+                  :param value: The value used to fill the array (used if mode is `PAD_CONSTANT`)
 
 
 
@@ -143,7 +143,7 @@ This file provides functions to pad a rectangular buffer or an image either:
                      mln::image2d<uint8_t> input(iroi);
                      mln::image2d<uint8_t> out(oroi);
 
-                     mln::copy_and_pad(input, out, mln::PAD_CONSTANT, 69);
+                     mln::copy_pad(input, out, mln::PAD_CONSTANT, 69);
 
 
 

doc/source/morpho.rst

@@ -47,6 +47,7 @@ Component Trees & Hierarchical Representations
     morpho/component_tree
     morpho/maxtree
     morpho/tos
+    morpho/multivariate_component_tree
     morpho/alphatree
     morpho/watershed_hierarchy
 

doc/source/morpho/component_tree.rst

@@ -295,7 +295,7 @@ to make an horizontal cut of this tree.
     :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(image2d<int> node_map, ranges::span<double> values) const

doc/source/morpho/multivariate_component_tree.rst

+Multivariate Component Tree
+===========================
+
+The component trees have been extended for multivariate images [Car19]_. Right
+now, only the multivariate tree of shapes [Car15]_ (MToS) has been implemented in
+Pylene.
+
+* Include :file:`<mln/morpho/mtos.hpp>`
+
+.. cpp:namespace:: mln::morpho
+
+.. cpp:function:: auto mtos(image2d<rgb8> ima, point2d pstart);
+
+    Compute the multivariate tree of shapes and returns a pair `(tree, node_map)`.
+    See :doc:`component_tree` for more information about the representation of tree.
+
+    :param ima: The input image in RGB format.
+    :param pstart: The rooting point
+    :return: A tree of type ``component_tree<void>`` (no values are related to the nodes of the tree since they do not have a natural value) and a map from image point to node tree.
+
+Notes
+-----
+
+* Before computing the MToS, the user should add a border to the image, with the values of this border set to the median value of the border of the original image.
+* The resulting nodemap has a domain size of ``4d - 3`` with ``d`` the input image domain.
+* The MToS not having values related to the nodes, the user has to compute a value for each node, such as the mean of each node (as shown below in example with the ``mean_node_accu`` accumulator).
+
+Example
+-------
+
+This example computes a grain filter, which removes all the node having an area inferior to 100.
+
+    ::
+
+        #include <mln/accu/accumulators/count.hpp>
+        #include <mln/core/extension/padding.hpp>
+        #include <mln/core/image/ndimage.hpp>
+        #include <mln/morpho/mtos.hpp>
+        #include <mln/io/imread.hpp>
+
+        // Function to reduce the nodemap to the original image domain
+        mln::image2d<int> reduce_nodemap(mln::image2d<int> n);
+
+        // Function to get the median of the border values
+        mln::rgb8 get_median_on_border(mln::image2d<mln::rgb8> ima);
+
+        // Accumulator to compute the mean of the pixel values of each node, without taking into account the values of the holes
+        struct mean_node_accu : mln::Accumulator<mean_node_accu>
+        {
+          using result_type = decltype(std::declval<mln::rgb8>() + std::declval<mln::rgb8>());
+        public:
+          void take(const mln::rgb8& v)
+          {
+            m_sum += v;
+            m_count++;
+          }
+
+          void take(const mean_node_accu&) {}
+
+          result_type to_result() const { return m_count > 1 ? static_cast<result_type>(m_sum / m_count) : m_sum; }
+
+        private:
+          result_type m_sum = {0, 0, 0};
+          int m_count = 0;
+        };
+
+        int main(void)
+        {
+            mln::image2d<mln::rgb8> ima;
+            mln::io::imread("lena.ppm", ima);
+
+            // Adding a border
+            const auto median = get_median_on_border(ima);
+            ima.inflate_domain(1);
+            constexpr int borders[2][2] = {{1, 1}, {1, 1}};
+            mln::pad(ima, mln::PAD_CONSTANT, borders, median);
+
+            // Compute the MToS
+            auto [t, nm] = mln::morpho::mtos(to_process, {-1, -1}); // The rooting point is in the added border
+
+            // Reduce the nodemap
+            nm = reduce_nodemap(nm);
+
+            // Remove the border
+            ima.inflate_domain(-1);
+            nm.inflate_domain(-1);
+
+            // Compute the area of each node of the tree
+            auto area = t.compute_attribute_on_points(nm, mln::accu::accumulators::count<int>());
+
+            // Perform a grain filter on the tree
+            t.filter(mln::morpho::CT_FILTER_DIRECT, nm, [&area](int n) { return area[n] >= 100; });
+
+            // Compute the mean of the connected component of each nodes
+            auto mean = t.compute_attribute_on_values(nm, ima, mean_node_accu());
+
+            // Reconstruct the tree
+            auto rec = t.reconstruct_from(nm, ranges::make_span(mean.data(), mean.size()));
+
+            return 0;
+        }
+
+.. list-table::
+
+    * - .. image:: /images/depth_map.png
+           :width: 100%
+
+      - .. image:: /images/mtos_rec.png
+           :width: 100%
+
+    * - The depth map resulting of the fusion of the trees (see [Car15]_ for more details)
+      - The reconstructed image from the filtered tree
+
+References
+----------
+
+.. [Car19] Edwin Carlinet and Thierry Géraud (2019). Introducing Multivariate Connected Openings and Closings. *International Symposium on Mathematical Morphology and Its Applications to Signal and Image Processing.* Springer, Cham. 215-227
+.. [Car15] Edwin Carlinet and Thierry Géraud (2015). MToS: A tree of shapes for multivariate images. *IEEE Transactions on Image Processing 24.12* 5330-5342
\ No newline at end of file

doc/source/morpho/watershed_hierarchy.rst

 Hierarchical Watershed
-==========
+======================
 
 * Include :file:`<mln/morpho/watershed_hierarchy.hpp>`
 

doc/source/snippets/CMakeLists.txt

@@ -64,6 +64,8 @@ add_image("watershed_hierarchy_example" "${PYLENE_IMAGE_DIR}/lena.pgm" watershed
 
 add_image("saliency_example" "${PYLENE_IMAGE_DIR}/lena.pgm" saliency_watershed.png)
 
+add_image("mtos_example" "${PYLENE_IMAGE_DIR}/lena.ppm" depth_map.png mtos_rec.png)
+
 add_custom_target(build-images
                   DEPENDS "${DOCUMENTATION_IMAGES}")
 
@@ -88,5 +90,6 @@ add_executable(area_filter area_filter.cpp)
 add_executable(cdt cdt.cpp)
 add_executable(first_start_1 first_start_1.cpp)
 add_executable(intro-1 intro-1.cpp)
+add_executable(mtos_example mtos_example.cpp)
 
 target_compile_definitions(erosion-cli PRIVATE BOOST_ALL_NO_LIB)
\ No newline at end of file

doc/source/snippets/mtos_example.cpp

+#include <mln/accu/accumulators/count.hpp>
+#include <mln/accu/accumulators/max.hpp>
+#include <mln/core/algorithm/accumulate.hpp>
+#include <mln/core/colors.hpp>
+#include <mln/core/extension/padding.hpp>
+#include <mln/core/image/ndimage.hpp>
+#include <mln/core/image/view/channel.hpp>
+#include <mln/core/image/view/transform.hpp>
+#include <mln/io/imread.hpp>
+#include <mln/io/imsave.hpp>
+#include <mln/morpho/tos.hpp>
+
+#include <mln/morpho/private/satmaxtree.hpp>
+#include <mln/morpho/private/trees_fusion.hpp>
+
+#include <iostream>
+#include <vector>
+
+#include "lut.hpp"
+
+namespace
+{
+  struct mean_node_accu : mln::Accumulator<mean_node_accu>
+  {
+    using result_type = decltype(std::declval<mln::rgb8>() + std::declval<mln::rgb8>());
+
+  public:
+    void take(const mln::rgb8& v)
+    {
+      m_sum += v;
+      m_count++;
+    }
+
+    void take(const mean_node_accu&) {}
+
+    result_type to_result() const { return m_count > 1 ? static_cast<result_type>(m_sum / m_count) : m_sum; }
+
+  private:
+    result_type m_sum   = {0, 0, 0};
+    int         m_count = 0;
+  };
+
+  mln::rgb8 get_median_on_border(mln::image2d<mln::rgb8> ima)
+  {
+    std::vector<mln::rgb8> border;
+    border.reserve(2 * ima.width() + 2 * ima.height() - 4);
+
+    for (int i = 0; i < ima.width(); i++)
+    {
+      border.push_back(ima(mln::point2d{i, 0}));
+      border.push_back(ima(mln::point2d{i, ima.height() - 1}));
+    }
+
+    for (int i = 1; i < ima.height() - 1; i++)
+    {
+      border.push_back(ima(mln::point2d{0, i}));
+      border.push_back(ima(mln::point2d{ima.width() - 1, i}));
+    }
+
+    std::sort(border.begin(), border.end(), [](const mln::rgb8& a, const mln::rgb8& b) {
+      int i = 0;
+      while (i < 3 && a[i] == b[i])
+        i++;
+      return i < 3 && a[i] < b[i];
+    });
+    return border[border.size() / 2];
+  }
+
+  /// \brief Reduce the size of the nodemap resulting from the MToS computation
+  ///        to fit the original image one
+  mln::image2d<int> reduce_nodemap(mln::image2d<int> n)
+  {
+    const mln::point2d pmin(n.domain().tl());
+    const mln::point2d pmax(n.domain().br());
+    mln::image2d<int>  res(mln::box2d{{(pmin[0] - 3) / 4, (pmin[1] - 3) / 4}, {(pmax[0] + 3) / 4, (pmax[1] + 3) / 4}});
+
+    for (int y = 0; y < res.height(); y++)
+    {
+      int* in  = n.buffer() + y * 4 * n.stride(1);
+      int* out = res.buffer() + y * res.stride();
+      for (int x = 0; x < res.width(); x++)
+      {
+        *out = *in;
+        in += 4;
+        out++;
+      }
+    }
+
+    return res;
+  }
+} // namespace
+
+int main(int argc, char* argv[])
+{
+  if (argc < 4)
+  {
+    std::cerr << "Invalid number of arguments\nUsage: " << argv[0]
+              << " input_filename depth_map_filename rec_filename\n";
+    return 1;
+  }
+
+  mln::image2d<mln::rgb8> ima;
+  mln::io::imread(argv[1], ima);
+
+  const auto med = get_median_on_border(ima);
+  ima.inflate_domain(1);
+  constexpr int borders[2][2] = {{1, 1}, {1, 1}};
+  mln::pad(ima, mln::PAD_CONSTANT, borders, med);
+
+  mln::morpho::component_tree<> trees[3];
+  mln::image2d<int>             nodemaps[3];
+  std::vector<int>              depths[3];
+
+  for (int c = 0; c < 3; c++)
+  {
+    std::tie(trees[c], nodemaps[c]) = mln::morpho::tos(mln::view::channel(ima, c), {-1, -1});
+    depths[c]                       = trees[c].compute_depth();
+  }
+
+  const auto [gos, tree_to_graph] = mln::morpho::details::compute_inclusion_graph(trees, nodemaps, depths, 3);
+  auto depth_map                  = mln::morpho::details::compute_depth_map(gos, tree_to_graph, nodemaps);
+  {
+    std::uint16_t max_depth = mln::accumulate(depth_map, mln::accu::accumulators::max<std::uint16_t>());
+    auto          normalized_depth =
+        mln::view::transform(depth_map, [&max_depth](std::uint16_t a) -> float { return (float)a / (float)max_depth; });
+    auto heat_depth = mln::view::transform(normalized_depth, heat_lut);
+    mln::io::imsave(heat_depth, argv[2]);
+  }
+
+  auto [t, nm] = mln::morpho::details::satmaxtree(depth_map, {-1, -1});
+  nm           = reduce_nodemap(nm);
+  ima.inflate_domain(-1);
+  nm.inflate_domain(-1);
+
+  auto area = t.compute_attribute_on_points(nm, mln::accu::accumulators::count<int>());
+  t.filter(mln::morpho::CT_FILTER_DIRECT, nm, [&area](int n) { return area[n] >= 100; });
+  auto mean = t.compute_attribute_on_values(nm, ima, mean_node_accu());
+  auto rec  = t.reconstruct_from(nm, ranges::make_span(mean.data(), mean.size()));
+  mln::io::imsave(mln::view::cast<mln::rgb8>(rec), argv[3]);
+
+  return 0;
+}
\ No newline at end of file

pylene/CMakeLists.txt

@@ -81,12 +81,15 @@ target_sources(Pylene-core PRIVATE
                src/io/imprint.cpp
                src/morpho/block_running_max.cpp
                src/morpho/component_tree.cpp
+               src/morpho/filters2d.cpp
                src/morpho/hvector.cpp
                src/morpho/hvector_unbounded.cpp
                src/morpho/immersion.cpp
                src/morpho/maxtree.cpp
+               src/morpho/mtos.cpp
+               src/morpho/satmaxtree.cpp
+               src/morpho/trees_fusion.cpp
                src/morpho/unionfind.cpp
-               src/morpho/filters2d.cpp
                )
 
 add_library(Pylene-io-freeimage)

pylene/include/mln/morpho/mtos.hpp

+#pragma once
+
+#include <mln/core/colors.hpp>
+#include <mln/morpho/component_tree.hpp>
+
+#include <tuple>
+
+namespace mln::morpho
+{
+  /// \brief Compute the multivariate tree of shapes
+  /// \param ima The input image
+  /// \param pstart The rooting point
+  /// \return A pair (tree, nodemap)
+  std::pair<component_tree<>, image2d<int>> mtos(image2d<rgb8> ima, point2d pstart);
+}
\ No newline at end of file

pylene/include/mln/morpho/private/satmaxtree.hpp

+#pragma once
+
+#include <mln/morpho/component_tree.hpp>
+
+#include <tuple>
+
+namespace mln::morpho::details
+{
+  /// \brief Compute the hole-filled maxtree of a depth map
+  /// \param depth_map The input depth map
+  /// \param pstart The rooting point of the satmaxtree
+  /// \return A pair (tree, nodemap) where tree is the resulting satmaxtree (values are not kept)
+  ///         and a map image point -> tree node
+  std::pair<component_tree<>, image2d<int>> satmaxtree(image2d<std::uint16_t> depth_map, point2d pstart);
+}
\ No newline at end of file

pylene/include/mln/morpho/private/trees_fusion.hpp

+#pragma once
+
+#include <mln/morpho/component_tree.hpp>
+
+#include <set>
+#include <tuple>
+#include <vector>
+
+namespace mln::morpho::details
+{
+  /// \brief Compute the inclusion graph of the trees given in input
+  /// \param trees The input trees
+  /// \param nodemaps The nodemaps related to each tree
+  /// \param depths The depths attribute of each tree
+  /// \param The number of tree given in argument
+  /// \return A pair (graph, tree_to_graph) where graph is the inclusion graph and tree_to_graph
+  ///         is a map (tree_index, tree_node) -> graph_node
+  std::pair<std::vector<std::set<int>>, std::vector<std::vector<int>>>
+  compute_inclusion_graph(component_tree<>* trees, image2d<int>* nodemaps, std::vector<int>* depths, int ntrees);
+
+  /// \brief Compute the depth map of an inclusion graph
+  /// \param graph The inclusion graph
+  /// \param tree_to_graph A map (tree_index, tree_node) -> graph_node
+  /// \param nodemaps The nodemaps of each tree which built the inclusion graph
+  /// \return The depth map of the inclusion graph (The depth should not exceed 2^16)
+  /// FIXME: For future optimization, the arguments should be change to
+  /// (const std::vector<std::set<int>>&    graph, const std::vector<int>& tree_to_graph, image2d<int>* nodemaps, const
+  /// std::vector<int>& tree_shapes)
+  image2d<std::uint16_t> compute_depth_map(const std::vector<std::set<int>>&    graph,
+                                           const std::vector<std::vector<int>>& tree_to_graph, image2d<int>* nodemaps);
+} // namespace mln::morpho::details
\ No newline at end of file

pylene/src/morpho/mtos.cpp

+#include <mln/core/image/view/channel.hpp>
+#include <mln/morpho/mtos.hpp>
+#include <mln/morpho/private/satmaxtree.hpp>
+#include <mln/morpho/private/trees_fusion.hpp>
+#include <mln/morpho/tos.hpp>
+
+namespace mln::morpho
+{
+  std::pair<component_tree<>, image2d<int>> mtos(image2d<rgb8> ima, point2d pstart)
+  {
+    mln::morpho::component_tree<> trees[3];
+    mln::image2d<int>             nodemaps[3];
+    std::vector<int>              depths[3];
+
+    for (int c = 0; c < 3; c++)
+    {
+      std::tie(trees[c], nodemaps[c]) = mln::morpho::tos(mln::view::channel(ima, c), pstart);
+      depths[c]                       = trees[c].compute_depth();
+    }
+
+    const auto [gos, tree_to_graph] = mln::morpho::details::compute_inclusion_graph(trees, nodemaps, depths, 3);
+    auto depth_map                  = mln::morpho::details::compute_depth_map(gos, tree_to_graph, nodemaps);
+    auto [t, nm]                    = mln::morpho::details::satmaxtree(depth_map, pstart);
+
+    return {std::move(t), std::move(nm)};
+  }
+} // namespace mln::morpho
\ No newline at end of file

pylene/src/morpho/satmaxtree.cpp

+#include <mln/morpho/private/satmaxtree.hpp>
+
+#include <mln/morpho/tos.hpp>
+
+namespace mln::morpho::details
+{
+  std::pair<component_tree<>, image2d<int>> satmaxtree(image2d<std::uint16_t> depth_map, point2d pstart)
+  {
+    auto [t, nm] = mln::morpho::tos(depth_map, pstart);
+
+    const int                  num_node = static_cast<int>(t.parent.size());
+    std::vector<std::uint16_t> delta(num_node, 0);
+    for (int i = 1; i < num_node; i++)
+    {
+      delta[i] = delta[t.parent[i]];
+      if (t.values[i] <= t.values[t.parent[i]])
+        delta[i] += t.values[t.parent[i]] - t.values[i];
+    }
+    for (int i = 1; i < num_node; i++)
+      t.values[i] += delta[i];
+
+    std::vector<bool> pred(t.parent.size());
+    for (int n = 1; n < (int)t.parent.size(); n++)
+      pred[n] = t.values[n] > t.values[t.parent[n]] || t.parent[n] == 0;
+    t.filter(CT_FILTER_DIRECT, nm, [&pred](int n) { return pred[n]; });
+
+    return {std::move(t), std::move(nm)};
+  }
+} // namespace mln::morpho::details
\ No newline at end of file

pylene/src/morpho/trees_fusion.cpp

+#include <mln/morpho/private/trees_fusion.hpp>
+
+#include <numeric>
+#include <stack>
+
+namespace mln::morpho::details
+{
+  namespace
+  {
+    std::vector<int> smallest_enclosing_shape(const component_tree<>& ti, image2d<int> ni, const component_tree<>& tj,
+                                              image2d<int> nj, const std::vector<int>& depth)
+    {
+      const auto lca = [&depth, &tj](int a, int b) {
+        if (a < 0)
+          return b;
+        if (b < 0)
+          return a;
+        while (depth[a] < depth[b])
+          b = tj.parent[b];
+        while (depth[b] < depth[a])
+          a = tj.parent[a];
+        while (a != b)
+        {
+          a = tj.parent[a];
+          b = tj.parent[b];
+        }
+        return a;
+      };
+
+      std::vector<int> res(ti.parent.size(), -1);
+      mln_foreach (auto p, ni.domain())
+        res[ni(p)] = lca(res[ni(p)], nj(p));
+
+      for (int n = (int)ti.parent.size() - 1; n > 0; n--)
+        res[ti.parent[n]] = lca(res[ti.parent[n]], res[n]);
+
+      return res;
+    }
+
+    std::vector<int> topo_sort(const std::vector<std::set<int>>& g)
+    {
+      static constexpr int UNMARKED  = 0;
+      static constexpr int TEMPORARY = 1;
+      static constexpr int PERMANENT = 2;
+
+      std::vector<int> res(g.size());
+      int              cur = static_cast<int>(g.size()) - 1;
+
+      std::vector<int> visited(g.size(), UNMARKED);
+      std::stack<int>  st;
+
+      for (int i = 0; i < (int)g.size(); i++)
+      {
+        st.push(i);
+
+        while (!st.empty())
+        {
+          const auto v = st.top();
+          st.pop();
+          if (visited[v] == PERMANENT)
+            continue;
+          else if (visited[v] == TEMPORARY)
+          {
+            res[cur--] = v;
+            visited[v] = PERMANENT;
+          }
+          else
+          {
+            visited[v] = TEMPORARY;
+            st.push(v);
+            for (int e : g[v])
+              st.push(e);
+          }
+        }
+      }
+
+      return res;
+    }
+  } // namespace
+
+  std::pair<std::vector<std::set<int>>, std::vector<std::vector<int>>>
+  compute_inclusion_graph(component_tree<>* trees, image2d<int>* nodemaps, std::vector<int>* depths, int ntrees)
+  {
+    std::vector<std::vector<int>> tree_to_graph(ntrees); // Link tree node -> graph node
+    std::vector<std::set<int>>    graph(1); // The graph (the container of out vertices is a set since there can only be