Commit a946bee1 authored by Baptiste Esteban's avatar Baptiste Esteban
Browse files

Merge branch 'development/alphatree' into 'next'

Improve alphatree in the library

See merge request !116
parents e1be240d 038e22a5
Pipeline #28327 passed with stages
in 26 minutes and 24 seconds
This diff is collapsed.
\documentclass[svgnames,tikz]{standalone}
\usepackage{pgfplots}
\begin{document}
\begin{tikzpicture}
\draw[y=0.7cm, step=1] (0, 0) grid (4, 4);
\draw[dashed, red] (0, 4) -- (4, 4);
\draw[dashed, red] (0, 6) -- (4, 6);
\draw[dashed, red] (0, 8) -- (4, 8);
\node[red] () at (-0.3, 4) {0};
\node[red] () at (-0.3, 6) {1};
\node[red] () at (-0.3, 8) {2};
\node[draw, circle, fill=white] (0) at (2, 8) {0};
\node[draw, circle, fill=white] (1) at (0.5, 6) {1};
\node[draw, circle, fill=white] (2) at (2, 6) {2};
\node[draw, circle, fill=white] (3) at (3.5, 4) {3};
\node[draw, circle, fill=white] (4) at (0.5, 4) {4};
\node[draw, circle, fill=white] (5) at (2, 4) {5};
\draw[->, >=latex] (1) -- (0);
\draw[->, >=latex] (2) -- (0);
\draw[->, >=latex] (3) -- (0);
\draw[->, >=latex] (4) -- (1);
\draw[->, >=latex] (5) -- (2);
\draw[blue, dashed] (2, 3.5) -- (2.5, 2.5);
\draw[blue, dashed] (2, 3.5) -- (2.5, 1.5);
\draw[blue, dashed] (2, 3.5) -- (1.5, 2.5);
\node () at (2, -0.3) {Node map};
\node () at (2, 8.6) {Component Tree};
\end{tikzpicture}
\end{document}
\ No newline at end of file
Alpha Tree
==========
Include :file:`<mln/morpho/alphatree.hpp>`
* Include :file:`<mln/morpho/alphatree.hpp>`
.. cpp:namespace:: mln::morpho
.. cpp:function:: auto alphatree(Image f, Neighborhood nbh);
.. cpp:function:: auto alphatree(Image f, Neighborhood nbh, F dist);
Compute the alphatree (also known as quasi-flat zone hierarchy) and returns a pair
Compute the alpha tree (also known as quasi-flat zone hierarchy) and returns a pair
`(tree, node_map)`. See :doc:`component_tree` for more information about the
representation of tree.
representation of tree. The implementation is based on the Kruskal algorithm [Naj13]_.
:param input: The input image
:param nbh: The neighborhood
:return: A pair `(tree, node_map)` where *tree* is of type
`component_tree<image_value_t<I>>` and `node_map` is a mapping from each point to the
node *id* in the tree.
:precondition: ``f.domain().has(pinf)``
:param nbh: The neighborhood
:param dist: The function weighting the edges between two pixels. For grayscale images, this is generally the absolute difference between two values.
: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
* ``image_value_t<I>`` is :cpp:concept:`std::totally_ordered`
* ``std::invoke_result_t<F, image_value_t<Image>, image_value_t<Image>>`` is :cpp:concept:`std::totally_ordered`
Definition
----------
The Alpha tree is the tree of :math:`\alpha`-connected components. An :math:`\alpha`-connected component in an image :math:`f`, for a pixel :math:`p`, is defined as
.. math::
\alpha-CC(p) = \{p\}\ \cup\ \{q\ |\ \text{there exists a path}\ \Pi = \{\pi_0, ..., \pi_n\}\ \text{between}\ p\ \text{and}\ q\ \text{such that}\ d(f(\pi_i), f(\pi_{i+1})) \leq \alpha\}
where :math:`d` is a dissimilarity function between two pixels of the image :math:`f`.
The :math:`\alpha`-connected components form an ordered sequence when :math:`\alpha` is growing, such that for :math:`\alpha_i < \alpha_j`,
:math:`\alpha_i-CC(p) \subseteq \alpha_j-CC(p)`. Thus, the alphatree is the hierarchy where the parenthood relationship represents the inclusion of the
:math:`\alpha`-connected components.
Representation
--------------
.. image:: /figures/morpho/alphatree_repr.svg
:align: center
:width: 30%
The ``alphatree`` function returns a tree and a node map. The tree has two attributes:
* ``parent``: The parenthood relationship of the node :math:`i`, representing the inclusion relationship of the :math:`\alpha`-connected components.
* ``values``: The value of :math:`\alpha` assigned to a node :math:`i`.
Then, the node map is the relation between a pixel of the image and its related node in the tree, a leaf for the case of the alphatree.
.. rubric:: Example
The image above illustrates the representation of the alpha tree in Pylene, the parenthood relationship being illustrated in arrows, the values of alpha, assigned to each node, being in red, and the relation
between a node of the tree and a pixel of the image being represented by blue dashed lines.
Example
-------
This example is used to generate the grayscale lena cut with a threshold of 3 below.
::
#include <mln/accu/accumulators/mean.hpp>
#include <mln/morpho/alphatree.hpp>
#include <mln/core/image/ndimage.hpp>
#include <mln/core/neighborhood/c4.hpp>
mln::image2d<uint8_t> input = ...;
// Compute the alpha tree
auto [tree, node_map] = mln::morpho::alphatree(input, mln::c4);
// Compute an attribute (for example the average pixels value at each node, as below)
auto mean = t.compute_attribute_on_values(node_map, input, mln::accu::accumulators::mean<uint8_t>());
// Making an horizontal cut of the tree
const auto threshold = 3; // Threshold of the horizontal cut, that means the lowest alpha in the cut
auto nodemap_cut = t.horizontal_cut(threshold, node_map); // Return a new nodemap associated to the cut
// Labelizing the cut with the mean values of each node
auto out = t.reconstruct_from(nodemap_cut, ranges::make_span(mean)); // Using range-v3 span
.. list-table::
* - .. figure:: /images/alphatree_cut_color.png
Cut of the alpha tree with a threshold of 10
- .. figure:: /images/alphatree_cut_gray.png
Cut of the alpha tree with a threshold of 3
Notes
-----
......@@ -48,3 +98,5 @@ 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
......@@ -148,9 +148,9 @@ A :cpp:class:`component_tree` `t` has the following methods.
.. cpp:namespace-push:: template <class V> component_tree
.. cpp:function:: auto accumulate_on_points(Image node_map, Accumulator accu) const
auto accumulate_on_values(Image node_map, Image values, Accumulator accu) const
auto accumulate_on_pixels(Image node_map, Image values, Accumulator accu) const
.. cpp:function:: auto compute_attribute_on_points(Image node_map, Accumulator accu) const
auto compute_attribute_on_values(Image node_map, Image values, Accumulator accu) const
auto compute_attribute_on_pixels(Image node_map, Image values, Accumulator accu) const
Accumulate the points of each component.
......@@ -190,7 +190,7 @@ Example: computing the size (number of points) of the components.
#include <mln/accu/accumulators/count.hpp>
auto sizes = tree.accumulate_on_points(node_map, mln::accu::features::count<>());
auto sizes = tree.compute_attribute_on_points(node_map, mln::accu::features::count<>());
.. cpp:function:: auto compute_depth() const
......@@ -255,6 +255,7 @@ A :cpp:class:`component_tree` `t` has the following methods.
auto pred = [&](int id) { nchildren > 1; };
t.filter(CT_FILTER_DIRECT, pred);
Image reconstruction
--------------------
......@@ -279,7 +280,19 @@ the image from the filtered tree.
t.filter(CT_FILTER_DIRECT, [&area](int id) { return area[id] > 20; });
auto out = t.reconstruct(node_map);
Horizontal cut
--------------
When the tree is a hierarchy of partition, such as the :doc:`alphatree`, it is possible
to make an horizontal cut of this tree.
.. cpp:function:: I horizontal_cut(const T threshold, I nodemap) const
I horizontal_cut_from_levels(const T threshold, I nodemap, ::ranges::span<V> levels) const
Make an horizontal cut at threshold ``threshold`` of the tree and return the nodemap associated to the cut.
:param threshold: The threshold of the horizontal cut
: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).
A complete example
------------------
......
......@@ -57,6 +57,8 @@ add_image("area_filter"
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_custom_target(build-images
DEPENDS "${DOCUMENTATION_IMAGES}")
......@@ -70,6 +72,7 @@ target_link_libraries(doc-lib Pylene::Pylene)
link_libraries(Pylene::Pylene)
link_libraries(doc-lib)
add_executable(alphatree_example alphatree_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)
......
#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/alphatree.hpp>
#include <iostream>
template <typename V>
void process_example(const mln::image2d<V>& img, const std::string& cut_filename, const double threshold)
{
// 2. Build the alphatree
auto [t, nm] = mln::morpho::alphatree(img, 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 alphatree
auto cut_nm = t.horizontal_cut(threshold, nm);
// 5. Labelize the cut
auto cut = t.reconstruct_from(cut_nm, ::ranges::make_span(mean));
// 5. Save the output cut
mln::io::imsave(mln::view::cast<V>(cut), cut_filename);
}
int main(int argc, char* argv[])
{
if (argc < 3)
{
std::cerr << "Invalid number of argument\nUsage: " << argv[0] << " input_filename cut_filename\n";
return 1;
}
auto in_filename = std::string(argv[1]);
auto cut_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, cut_filename, 3);
}
else if (in.sample_type() == mln::sample_type_id::RGB8)
{
const auto* img = in.cast_to<mln::rgb8, 2>();
process_example(*img, cut_filename, 10);
}
else
{
std::cerr << "Unhandled sample type format\n";
return 1;
}
return 0;
}
\ No newline at end of file
......@@ -51,7 +51,7 @@ namespace mln
};
public:
extend_by_projection_view_base<I, Proj>(I ima, Proj proj)
extend_by_projection_view_base(I ima, Proj proj)
: base_t{std::move(ima)}
, m_proj{std::move(proj)}
{
......
......@@ -4,13 +4,15 @@
#include <mln/core/concepts/image.hpp>
#include <mln/core/concepts/neighborhood.hpp>
#include <mln/core/functional_ops.hpp>
#include <mln/core/algorithm/for_each.hpp>
#include <mln/morpho/component_tree.hpp>
#include <mln/core/functional_ops.hpp>
#include <mln/morpho/canvas/unionfind.hpp>
#include <mln/morpho/component_tree.hpp>
#include <mln/morpho/private/directional_hqueue.hpp>
#include <type_traits>
#include <algorithm>
#include <type_traits>
#include <range/v3/functional/concepts.hpp>
......@@ -33,26 +35,71 @@ namespace mln::morpho
/******************************************/
namespace internal
{
/// \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>
class alphatree_edges;
template <class P, class W>
struct alphatree_edge_t
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>
{
P a;
P b;
W w;
public:
void push(int dir, W w, P p) { m_cont.insert(dir, w, p); }
std::tuple<P, P, W> pop() { return m_cont.pop(); }
W top() const { return m_cont.current_level(); }
bool empty() const { return m_cont.empty(); }
void on_finish_insert() {}
private:
details::directional_hqueue<P, N, W> m_cont;
};
template <typename P, typename W>
struct edge_t
{
P p;
P q;
W w;
};
template <typename P, typename N, typename W>
class alphatree_edges
{
public:
void push(int dir, W w, P p) { m_cont.push_back({p, p + cn.after_offsets()[dir], w}); }
std::tuple<P, P, W> pop()
{
assert(m_current < m_cont.size());
const auto e = m_cont[m_current++];
return {std::move(e.p), std::move(e.q), std::move(e.w)};
}
W top() const
{
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; });
}
private:
static constexpr auto cn = N();
std::vector<edge_t<P, W>> m_cont;
std::size_t m_current = 0;
};
// Compute a node_id for each flat zone
template <class I, class J>
[[gnu::noinline]]
std::size_t alphatree_create_nodemap(I node_map, J zpar)
[[gnu::noinline]] std::size_t alphatree_create_nodemap(I node_map, J zpar)
{
std::size_t node_count = 0;
mln_foreach(auto px, node_map.pixels())
mln_foreach (auto px, node_map.pixels())
{
auto p = px.point();
auto rp = canvas::impl::zfindroot(zpar, p);
......@@ -66,21 +113,18 @@ namespace mln::morpho
// Compute flat zone of the image using union-find structures
// It returns the number of edges processed
template <class E, class J>
std::size_t alphatree_compute_flatzones(const E* edges, std::size_t n, J zpar)
void alphatree_compute_flatzones(E& edges, J zpar)
{
canvas::impl::union_find_init_par(zpar);
std::size_t i;
for (i = 0; i < n && edges[i].w == 0; ++i)
while (!edges.empty() && edges.top() == 0)
{
auto p = edges[i].a;
auto q = edges[i].b;
auto rp = canvas::impl::zfindroot(zpar, p);
auto rq = canvas::impl::zfindroot(zpar, q);
const auto [p, q, w] = edges.pop();
auto rp = canvas::impl::zfindroot(zpar, p);
auto rq = canvas::impl::zfindroot(zpar, q);
if (rp != rq) // Merge set
zpar(rq) = rp;
}
return i;
}
// Inverse the node array, updating the parent relation
......@@ -104,39 +148,45 @@ namespace mln::morpho
}
template <class I, class N, class F, class E>
void alphatree_compute_edges(I input, N nbh, F distance, std::vector<E>& edges)
template <class I, class N, class F, class C>
void alphatree_compute_edges(I input, N nbh, F distance, C& edges)
{
static_assert(is_a_v<I, mln::details::Image>);
auto dom = input.domain();
mln_foreach (auto p, dom)
{
int i = 0;
for (auto n : nbh.after(p))
{
if (dom.has(n))
edges.push_back({p, n, distance(input(p), input(n))});
edges.push(i, distance(input(p), input(n)), p);
++i;
}
}
edges.on_finish_insert();
}
//
//
template <class E, class I, class W>
std::size_t alphatree_compute_hierarchy(const E* edges, std::size_t n, I node_map, //
std::size_t node_count, //
std::vector<int>& par, //
std::vector<W>& levels)
template <class E, class I, class W, class M>
std::size_t alphatree_compute_hierarchy(E& edges, I node_map, //
std::size_t node_count, //
std::vector<int>& par, //
std::vector<W>& levels, //
std::vector<M>* mst)
{
static_assert(mln::is_a<I, mln::details::Image>());
std::vector<int> zpar(node_count); // Parent of ufind structure
std::vector<int> links(node_count); // Node index in the alpha tree
std::vector<int> zpar(node_count); // Parent of ufind structure
std::vector<int> links(node_count); // Node index in the alpha tree
std::iota(std::begin(zpar), std::end(zpar), 0);
std::iota(std::begin(links), std::end(links), 0);
for (std::size_t i = 0; i < n; ++i)
while (!edges.empty())
{
auto [p, q, w] = edges[i];
auto [p, q, w] = edges.pop();
auto rp = canvas::impl::zfindroot(zpar.data(), node_map(p));
auto rq = canvas::impl::zfindroot(zpar.data(), node_map(q));
......@@ -150,14 +200,18 @@ namespace mln::morpho
// Do we need to create a new node
int rp_root = links[rp];
int rq_root = links[rq];
int max_root = std::max(rp_root, rq_root);
int min_root = std::min(rp_root, rq_root);
int new_root_id;
if (levels[rq_root] == w)
if (levels[max_root] == w)
{
new_root_id = rq_root;
new_root_id = max_root;
}
else if (levels[rp_root] == w)
else if (levels[min_root] == w)
{
new_root_id = rp_root;
new_root_id = min_root;
}
else
{
......@@ -166,74 +220,120 @@ namespace mln::morpho
par.push_back(new_root_id);
levels.push_back(w);
}
par[rp_root] = new_root_id;
par[rq_root] = new_root_id;
links[rp] = new_root_id;
if (mst != nullptr)
mst->push_back({p, q, w});
}
}
return node_count;
}
} // namespace internal
template <class I, class N, class F>
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)
{
static_assert(mln::is_a<I, mln::details::Image>());
static_assert(mln::is_a<N, mln::details::Neighborhood>());
static_assert(::ranges::cpp20::invocable<F, image_value_t<I>, image_value_t<I>>);
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::vector<int> canonized_par;
std::vector<W> canonized_levels;
using V = image_value_t<I>;
using P = image_point_t<I>;
using W = std::invoke_result_t<F, V, V>;
// Root initialization
canonized_par.push_back(0);
canonized_levels.push_back(levels[0]);
static_assert(::concepts::totally_ordered<W>);
std::vector<int> translation_map(node_count);
using edge_t = internal::alphatree_edge_t<P, W>;
translation_map[0] = 0;
int count = 1;
// Build canonized component tree
for (std::size_t i = 1; i < node_count; ++i)
{
if (levels[i] != levels[par[i]]) // Keep the node: Update tree
{
translation_map[i] = count++;
canonized_par.push_back(translation_map[par[i]]);
canonized_levels.push_back(levels[i]);
}
else // Deleted node: retrieve the parent translation
translation_map[i] = translation_map[par[i]];
}
return {canonized_par, canonized_levels};
}
// 1. Get the list of edges
std::vector<edge_t> edges;
internal::alphatree_compute_edges(std::move(input), std::move(nbh), std::move(distance), edges);
template <class I, class N, class F,
class M = edge_t<image_point_t<I>, std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>>
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, std::vector<M>* mst = nullptr, bool canonize_tree = true)
{
static_assert(mln::is_a<I, mln::details::Image>());
static_assert(mln::is_a<N, mln::details::Neighborhood>());
static_assert(::ranges::cpp20::invocable<F, image_value_t<I>, image_value_t<I>>);
// Sort the edges by increasing weights
std::sort(std::begin(edges), std::end(edges), [](const auto& a, const auto& b) { return a.w < b.w; });
using V = image_value_t<I>;
using P = image_point_t<I>;
using W = std::invoke_result_t<F, V, V>;
std::size_t edges_null_count;
std::size_t flatzones_count;
image_ch_value_t<I, int> node_map = imchvalue<int>(input).set_init_value(-1);
{
image_ch_value_t<I, P> zpar = imchvalue<P>(input);
// 2. Compute flat zone of the image
edges_null_count = internal::alphatree_compute_flatzones(edges.data(), edges.size(), zpar);
static_assert(::concepts::totally_ordered<W>);