Skip to content
GitLab
Commit ce609a55 authored by Quentin Kaci's avatar Quentin Kaci
Browse files

Provide default attribute computation for hierarchical watershed

parent a9c060a7
Hide whitespace changes
Inline Side-by-side
pylene/include/mln/morpho/watershed_hierarchy.hpp
View file @ ce609a55
......@@ -12,11 +12,22 @@ namespace mln::morpho
/// \param neighborhood The neighborhood relation
/// \param distance Distance function
template <class I, class A, class N, class F = mln::functional::l2dist_t<>>
std::pair<component_tree<typename std::invoke_result_t<
A, component_tree<std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>,
image_ch_value_t<I, int>>::value_type>,
image_ch_value_t<I, int>> //
watershed_hierarchy(I input, A attribute_func, N nbh, F distance = F{});
auto watershed_hierarchy(I input, A attribute_func, N nbh, F distance = F{});
enum WatershedAttribute
{
HEIGHT = 0,
DYNAMIC
};
/// Compute the watershed hierarchy of an image
///
/// \param input The input image
/// \param watershed_attribute Enumeration for default attribute computation
/// \param neighborhood The neighborhood relation
/// \param distance Distance function
template <class I, class N, class F = mln::functional::l2dist_t<>>
auto watershed_hierarchy(I input, WatershedAttribute watershed_attribute, N nbh, F distance = F{});
/******************************************/
......@@ -26,6 +37,110 @@ namespace mln::morpho
namespace internal
{
template <typename W, class I>
std::vector<W> height_attribute(component_tree<W> tree, I node_map)
{
int n = static_cast<int>(tree.parent.size());
int nb_leaves = static_cast<int>(node_map.domain().size());
std::vector<W> deepest_altitude(n, std::numeric_limits<W>::max());
for (int i = n - 1; i >= 0; --i)
{
int parent = tree.parent[i];
if (i > (n - nb_leaves - 1))
deepest_altitude[i] = tree.values[parent];
deepest_altitude[parent] = std::min(deepest_altitude[parent], deepest_altitude[i]);
}
std::vector<W> height(n);
for (int i = n - 1; i >= 0; --i)
height[i] = tree.values[tree.parent[i]] - deepest_altitude[i];
return height;
}
template <typename W, class I>
std::vector<bool> extrema_attribute(component_tree<W> tree, I node_map)
{
int n = static_cast<int>(tree.parent.size());
int nb_leaves = static_cast<int>(node_map.domain().size());
std::vector<bool> extrema(n, true);
std::fill_n(extrema.end() - nb_leaves, nb_leaves, false);
for (int i = (n - nb_leaves - 1); i > 0; --i)
{
int parent = tree.parent[i];
bool same_weight = tree.values[i] == tree.values[parent];
extrema[parent] = extrema[parent] && same_weight && extrema[i];
extrema[i] = extrema[i] && !same_weight;
}
return extrema;
}
template <typename W, class I>
std::vector<W> dynamic_attribute(component_tree<W> tree, I node_map)
{
int n = static_cast<int>(tree.parent.size());
int nb_leaves = static_cast<int>(node_map.domain().size());
std::vector<W> deepest_altitude(n, std::numeric_limits<W>::max());
std::vector<int> path_to_minima(n, -1);
// Compute deepest altitude and path to deepest minima
for (int i = (n - nb_leaves - 1); i >= 0; --i)
{
int parent = tree.parent[i];
// Deepest non leaf node
if (deepest_altitude[i] == std::numeric_limits<W>::max())
deepest_altitude[i] = tree.values[i];
if (deepest_altitude[i] < deepest_altitude[parent])
{
deepest_altitude[parent] = deepest_altitude[i];
path_to_minima[parent] = i;
}
}
auto extrema = internal::extrema_attribute(tree, node_map);
std::vector<int> nearest_minima(n, -1);
std::vector<W> dynamic(n);
dynamic[0] = tree.values[0] - deepest_altitude[0];
for (int i = 1; i < n; ++i)
{
int parent = tree.parent[i];
if (i > (n - nb_leaves - 1))
{
if (nearest_minima[parent] != -1)
dynamic[i] = dynamic[nearest_minima[parent]];
else
dynamic[i] = 0;
continue;
}
if (extrema[i])
nearest_minima[i] = i;
else
nearest_minima[i] = nearest_minima[parent];
if (i == path_to_minima[parent])
dynamic[i] = dynamic[parent];
else
dynamic[i] = tree.values[parent] - deepest_altitude[i];
}
return dynamic;
}
template <typename W, typename A>
std::vector<A> get_computed_attribute(const component_tree<W>& tree, const std::vector<A>& attribute, int nb_leaves)
{
......@@ -78,125 +193,37 @@ namespace mln::morpho
}
} // namespace internal
template <typename W, class I>
std::vector<W> height_attribute(component_tree<W> tree, I node_map)
{
int n = static_cast<int>(tree.parent.size());
int nb_leaves = static_cast<int>(node_map.domain().size());
std::vector<W> deepest_altitude(n, std::numeric_limits<W>::max());
for (int i = n - 1; i >= 0; --i)
{
int parent = tree.parent[i];
if (i > (n - nb_leaves - 1))
deepest_altitude[i] = tree.values[parent];
deepest_altitude[parent] = std::min(deepest_altitude[parent], deepest_altitude[i]);
}
std::vector<W> height(n);
for (int i = n - 1; i >= 0; --i)
height[i] = tree.values[tree.parent[i]] - deepest_altitude[i];
return height;
}
template <typename W, class I>
std::vector<bool> extrema_attribute(component_tree<W> tree, I node_map)
template <class I, class A, class N, class F>
auto watershed_hierarchy(I input, A attribute_func, N nbh, F distance)
{
int n = static_cast<int>(tree.parent.size());
int nb_leaves = static_cast<int>(node_map.domain().size());
std::vector<bool> extrema(n, true);
std::fill_n(extrema.end() - nb_leaves, nb_leaves, false);
std::vector<edge_t<image_point_t<I>, std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>> mst;
auto [tree, nm] = internal::__alphatree<false>(input, nbh, distance, false, false, &mst);
for (int i = (n - nb_leaves - 1); i > 0; --i)
{
int parent = tree.parent[i];
auto attribute = attribute_func(tree, nm);
bool same_weight = tree.values[i] == tree.values[parent];
extrema[parent] = extrema[parent] && same_weight && extrema[i];
extrema[i] = extrema[i] && !same_weight;
}
auto node_count = tree.parent.size();
mln::for_each(nm, [node_count](int& id) { id = static_cast<int>(node_count) - id - 1; });
return extrema;
return internal::watershed<I, N>(tree, nm, mst, attribute, input.domain().size());
}
template <typename W, class I>
std::vector<W> dynamic_attribute(component_tree<W> tree, I node_map)
template <class I, class N, class F>
auto watershed_hierarchy(I input, WatershedAttribute watershed_attribute, N nbh, F distance)
{
int n = static_cast<int>(tree.parent.size());
int nb_leaves = static_cast<int>(node_map.domain().size());
std::vector<W> deepest_altitude(n, std::numeric_limits<W>::max());
std::vector<int> path_to_minima(n, -1);
// Compute deepest altitude and path to deepest minima
for (int i = (n - nb_leaves - 1); i >= 0; --i)
{
int parent = tree.parent[i];
// Deepest non leaf node
if (deepest_altitude[i] == std::numeric_limits<W>::max())
deepest_altitude[i] = tree.values[i];
using W = std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>;
std::function<std::vector<W>(component_tree<W>, image_ch_value_t<I, int>)> attribute_func;
if (deepest_altitude[i] < deepest_altitude[parent])
{
deepest_altitude[parent] = deepest_altitude[i];
path_to_minima[parent] = i;
}
}
auto extrema = extrema_attribute(tree, node_map);
std::vector<int> nearest_minima(n, -1);
std::vector<W> dynamic(n);
dynamic[0] = tree.values[0] - deepest_altitude[0];
for (int i = 1; i < n; ++i)
switch (watershed_attribute)
{
int parent = tree.parent[i];
if (i > (n - nb_leaves - 1))
{
if (nearest_minima[parent] != -1)
dynamic[i] = dynamic[nearest_minima[parent]];
else
dynamic[i] = 0;
continue;
}
if (extrema[i])
nearest_minima[i] = i;
else
nearest_minima[i] = nearest_minima[parent];
if (i == path_to_minima[parent])
dynamic[i] = dynamic[parent];
else
dynamic[i] = tree.values[parent] - deepest_altitude[i];
case HEIGHT:
attribute_func = [](auto tree, auto nm) -> std::vector<W> { return internal::height_attribute(tree, nm); };
break;
case DYNAMIC:
attribute_func = [](auto tree, auto nm) -> std::vector<W> { return internal::dynamic_attribute(tree, nm); };
break;
}
return dynamic;
return watershed_hierarchy(input, attribute_func, nbh, distance);
}
template <class I, class A, class N, class F>
std::pair<component_tree<typename std::invoke_result_t<
A, component_tree<std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>,
image_ch_value_t<I, int>>::value_type>,
image_ch_value_t<I, int>> //
watershed_hierarchy(I input, A attribute_func, N nbh, F distance)
{
std::vector<edge_t<image_point_t<I>, std::invoke_result_t<F, image_value_t<I>, image_value_t<I>>>> mst;
auto [tree, nm] = internal::__alphatree<false>(input, nbh, distance, false, false, &mst);
auto attribute = attribute_func(tree, nm);
auto node_count = tree.parent.size();
mln::for_each(nm, [node_count](int& id) { id = static_cast<int>(node_count) - id - 1; });
return internal::watershed<I, N>(tree, nm, mst, attribute, input.domain().size());
}
} // namespace mln::morpho
\ No newline at end of file
tests/morpho/alphatree.cpp
View file @ ce609a55
......@@ -54,7 +54,7 @@ TEST(Morpho, AlphaTree)
{107, 73, 125, 157, 117}, //
};
auto [t, node_map] = mln::morpho::alphatree(
ima, mln::c4, [](const auto& a, const auto& b) -> std::uint8_t { return mln::functional::l2dist_t<>()(a, b); });
ima, mln::c4, [](const auto& a, const auto& b) -> std::uint8_t { return mln::functional::l2dist(a, b); });
mln::image2d<int> ref_0 = {
{0, 1, 2, 3, 4}, //
......@@ -161,7 +161,7 @@ TEST(Morpho, AlphaTreeMST)
std::vector<E> mst;
auto [t, _] = mln::morpho::internal::__alphatree(
ima, mln::c4, [](const auto& a, const auto& b) -> W { return mln::functional::l2dist_t<>()(a, b); }, true, true,
ima, mln::c4, [](const auto& a, const auto& b) -> W { return mln::functional::l2dist(a, b); }, true, true,
&mst);
for (std::size_t i = 0; i < expected_mst.size(); ++i)
......@@ -182,7 +182,7 @@ TEST(Morpho, AlphaTreeRGB8Uint16Distance)
};
auto [t, nm] = mln::morpho::alphatree(
ima, mln::c4, [](const auto& a, const auto& b) -> std::uint16_t { return mln::functional::l2dist_t<>()(a, b); });
ima, mln::c4, [](const auto& a, const auto& b) -> std::uint16_t { return mln::functional::l2dist(a, b); });
const mln::image2d<int> ref_0 = {
{0, 0, 1}, //
......@@ -216,7 +216,7 @@ TEST(Morpho, AlphaTreeRGB8FloatDistance)
};
auto [t, nm] = mln::morpho::alphatree(
ima, mln::c4, [](const auto& a, const auto& b) -> std::float_t { return mln::functional::l2dist_t<>()(a, b); });
ima, mln::c4, [](const auto& a, const auto& b) -> std::float_t { return mln::functional::l2dist(a, b); });
const mln::image2d<int> ref_0 = {
{0, 0, 1}, //
......@@ -252,7 +252,7 @@ TEST(Morpho, AlphaTree8C)
};
auto [t, nm] = mln::morpho::alphatree(
ima, mln::c8, [](const auto& a, const auto& b) -> std::uint32_t { return mln::functional::l2dist_t<>()(a, b); });
ima, mln::c8, [](const auto& a, const auto& b) -> std::uint32_t { return mln::functional::l2dist(a, b); });
mln::image2d<int> ref_0 = {
{0, 1, 2, 3, 4}, //
......@@ -285,7 +285,7 @@ TEST(Morpho, AlphaTree8CHQUEUE)
};
auto [t, nm] = mln::morpho::alphatree(
ima, mln::c8, [](const auto& a, const auto& b) -> std::uint8_t { return mln::functional::l2dist_t<>()(a, b); });
ima, mln::c8, [](const auto& a, const auto& b) -> std::uint8_t { return mln::functional::l2dist(a, b); });
mln::image2d<int> ref_0 = {
{0, 1, 2, 3, 4}, //
......@@ -316,7 +316,7 @@ TEST(Morpho, AlphaTree3DImage)
};
auto [t, nm] = mln::morpho::alphatree(
ima, mln::c26, [](const auto& a, const auto& b) -> std::uint32_t { return mln::functional::l2dist_t<>()(a, b); });
ima, mln::c26, [](const auto& a, const auto& b) -> std::uint32_t { return mln::functional::l2dist(a, b); });
const mln::image3d<int> ref_0 = {
{{0, 1, 1}, {0, 1, 1}, {2, 3, 4}}, //
......@@ -343,7 +343,7 @@ TEST(Morpho, AlphaTree3DImageHQUEUE)
};
auto [t, nm] = mln::morpho::alphatree(
ima, mln::c26, [](const auto& a, const auto& b) -> std::uint8_t { return mln::functional::l2dist_t<>()(a, b); });
ima, mln::c26, [](const auto& a, const auto& b) -> std::uint8_t { return mln::functional::l2dist(a, b); });
const mln::image3d<int> ref_0 = {
{{0, 1, 1}, {0, 1, 1}, {2, 3, 4}}, //
......@@ -388,7 +388,7 @@ TEST(Morpho, AlphaTreeCutMeanLabelization)
auto val = t.compute_attribute_on_values(nm, ima, mln::accu::accumulators::mean<std::uint8_t, std::uint8_t>());
// Renaming t and nm because clang does not allow to capture structured bindings
auto make_cut = [&lt = t, &lnm = nm, &val](const typename decltype(t.values)::value_type threshold) {
auto make_cut = [&lt = t, &lnm = nm, &val](auto threshold) {
auto lbl = lt.horizontal_cut(threshold, lnm);
return lt.reconstruct_from(lbl, ::ranges::make_span(val));
};
......
tests/morpho/watershed_hierarchy.cpp
View file @ ce609a55
......@@ -25,10 +25,10 @@ TEST(Morpho, HeightAttribute)
// Build BPT
auto [tree, nm] = mln::morpho::internal::__alphatree<false>(
input, mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); }, false,
input, mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); }, false,
false);
auto attribute = mln::morpho::height_attribute(tree, nm);
auto attribute = mln::morpho::internal::height_attribute(tree, nm);
ASSERT_EQ(expected_attribute.size(), attribute.size());
......@@ -48,7 +48,7 @@ TEST(Morpho, DynamicAttribute)
// Build BPT
auto [tree, nm] = mln::morpho::internal::__alphatree<false>(
input, mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); }, false,
input, mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); }, false,
false);
std::vector<bool> expected_extrema_attribute = {false, false, false, false, false, false, false, false, false, true,
......@@ -57,7 +57,7 @@ TEST(Morpho, DynamicAttribute)
false, false, false, false, false, false, false, false, false, false,
false, false, false, false, false, false, false, false, false};
auto extrema_attribute = mln::morpho::extrema_attribute(tree, nm);
auto extrema_attribute = mln::morpho::internal::extrema_attribute(tree, nm);
ASSERT_EQ(expected_extrema_attribute.size(), extrema_attribute.size());
......@@ -69,7 +69,7 @@ TEST(Morpho, DynamicAttribute)
134.f, 70.f, 70.f, 70.f, 70.f, 70.f, 70.f, 0, 134.f, 134.f, 134.f, 134.f, 38.f, 38.f, 38.f, 0, 134.f,
70.f, 70.f, 70.f, 70.f, 38.f, 0, 29.f, 29.f, 70.f, 38.f, 70.f, 70.f, 29.f, 70.f, 0};
auto dynamic_attribute = mln::morpho::dynamic_attribute(tree, nm);
auto dynamic_attribute = mln::morpho::internal::dynamic_attribute(tree, nm);
ASSERT_EQ(expected_dynamic_attribute.size(), dynamic_attribute.size());
......@@ -98,7 +98,70 @@ TEST(Morpho, AreaWatershedHierarchyGray)
[](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) -> float { return mln::functional::l2dist_t<>()(a, b); });
mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); });
ASSERT_EQ(expected_parent.size(), tree.parent.size());
ASSERT_EQ(expected_values.size(), tree.values.size());
ASSERT_EQ(tree.parent.size(), tree.values.size());
for (std::size_t i = 0; i < expected_parent.size(); ++i)
{
EXPECT_EQ(expected_parent[i], tree.parent[i]);
EXPECT_EQ(expected_values[i], tree.values[i]);
}
}
TEST(Morpho, DynamicWatershedHierarchyGray)
{
mln::image2d<uint8_t> input = {
{163, 112, 42, 121, 112}, //
{42, 121, 1, 42, 255}, //
{1, 112, 121, 112, 121}, //
{112, 255, 42, 1, 42}, //
{121, 112, 121, 112, 163}, //
};
std::vector<int> expected_parent = {0, 0, 1, 0, 3, 4, 5, 6, 1, 8, 9, 2, 7, 11, 13,
14, 15, 1, 13, 2, 19, 16, 17, 0, 23, 20, 21, 13};
std::vector<unsigned long> expected_values = {70, 38, 29, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
auto [tree, _] = mln::morpho::watershed_hierarchy(
input,
mln::morpho::WatershedAttribute::DYNAMIC,
mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); });
ASSERT_EQ(expected_parent.size(), tree.parent.size());
ASSERT_EQ(expected_values.size(), tree.values.size());
ASSERT_EQ(tree.parent.size(), tree.values.size());
for (std::size_t i = 0; i < expected_parent.size(); ++i)
{
EXPECT_EQ(expected_parent[i], tree.parent[i]);
EXPECT_EQ(expected_values[i], tree.values[i]);
}
}
TEST(Morpho, HeightWatershedHierarchyGray)
{
mln::image2d<uint8_t> input = {
{163, 112, 42, 121, 112}, //
{42, 121, 1, 42, 255}, //
{1, 112, 121, 112, 121}, //
{112, 255, 42, 1, 42}, //
{121, 112, 121, 112, 163}, //
};
std::vector<int> expected_parent = {0, 0, 1, 0, 3, 4, 5, 6, 1, 8, 9, 2, 7, 11, 13, 14, 15, 1, 13, 2, 19, 16, 17, 0, 23, 20, 21, 13};
std::vector<unsigned long> expected_values = {70, 38, 29, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
auto [tree, _] = mln::morpho::watershed_hierarchy(
input,
mln::morpho::WatershedAttribute::HEIGHT,
mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); });
ASSERT_EQ(expected_parent.size(), tree.parent.size());
ASSERT_EQ(expected_values.size(), tree.values.size());
......@@ -132,7 +195,7 @@ TEST(Morpho, AreaWatershedHierarchyGrayHQ)
[](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::uint8_t { return mln::functional::l2dist_t<>()(a, b); });
mln::c4, [](const auto& a, const auto& b) -> std::uint8_t { return mln::functional::l2dist(a, b); });
ASSERT_EQ(expected_parent.size(), tree.parent.size());
ASSERT_EQ(expected_values.size(), tree.values.size());
......@@ -166,7 +229,7 @@ TEST(Morpho, AreaWatershedHierarchyRGB)
[](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) -> float { return mln::functional::l2dist_t<>()(a, b); });
mln::c4, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); });
ASSERT_EQ(expected_parent.size(), tree.parent.size());
ASSERT_EQ(expected_values.size(), tree.values.size());
......@@ -200,7 +263,7 @@ TEST(Morpho, AreaWatershedHierarchyGrayC8)
[](auto tree, auto nm) -> std::vector<size_t> {
return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
},
mln::c8, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); });
mln::c8, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); });
ASSERT_EQ(expected_parent.size(), tree.parent.size());
ASSERT_EQ(expected_values.size(), tree.values.size());
......@@ -232,7 +295,7 @@ TEST(Morpho, AreaWatershedHierarchy3DImage)
[](auto tree, auto nm) -> std::vector<size_t> {
return tree.compute_attribute_on_points(nm, mln::accu::features::count<>());
},
mln::c26, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist_t<>()(a, b); });
mln::c26, [](const auto& a, const auto& b) -> float { return mln::functional::l2dist(a, b); });
ASSERT_EQ(expected_parent.size(), tree.parent.size());
ASSERT_EQ(expected_values.size(), tree.values.size());
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment