Commit 9a690ad2 authored by Edwin Carlinet's avatar Edwin Carlinet
Browse files

Add tos algorithm

parent 037e4904
Pipeline #15755 failed with stages
in 6 minutes and 56 seconds
......@@ -31,6 +31,13 @@ namespace mln::morpho::experimental::details
P front(int level) const noexcept;
bool empty(int level) const noexcept;
/// Get the level 𝑙 of the lowest non-empty queue with level ≤ 𝑙
/// such that every queue at level x with level ≤ x < 𝑙 are empty.
int lower_bound(int level) const noexcept;
/// Get the level 𝑙 of the highest non-empty queue with 𝑙 < level
/// such that every queue at level x with 𝑙 < x ≤ level are empty.
int upper_bound(int level) const noexcept;
private:
struct node_t
......@@ -128,5 +135,27 @@ namespace mln::morpho::experimental::details
return m_lists[level].tail;
}
template <int N, class P, class LinkImage>
inline int hlinked_lists<N, P, LinkImage>::lower_bound(int level) const noexcept
{
mln_precondition(level < N);
while (level < N && m_lists[level].size == 0)
++level;
return level;
}
template <int N, class P, class LinkImage>
inline int hlinked_lists<N, P, LinkImage>::upper_bound(int level) const noexcept
{
mln_precondition(level < N);
while (level >= 0 && m_lists[level].size == 0)
--level;
return level;
}
} // namespace mln::morpho::experimental::detail
#pragma once
#include <mln/core/concept/new/images.hpp>
#include <mln/core/algorithm/fill.hpp>
#include <mln/core/extension/border_management.hpp>
#include <mln/core/extension/fill.hpp>
#include <mln/core/trace.hpp>
#include <mln/core/neighborhood/c4.hpp>
#include <mln/core/neighborhood/c6.hpp>
#include <mln/morpho/experimental/private/pset.hpp>
namespace mln::morpho::experimental::details
{
/// Propagation
template <class I>
void propagation(I inf, I sup, image_ch_value_t<I, int> out, image_point_t<I> pstart, int& max_depth);
/******************************************/
/**** Implementation ****/
/******************************************/
template <class I>
void propagation(I inf, I sup, image_ch_value_t<I, int> ord, image_point_t<I> pstart, int& max_depth)
{
mln_entering("mln::morpho::details::propagation");
using P = image_point_t<I>;
using V = image_value_t<I>;
static_assert(P::ndim == 2 || P::ndim == 3, "Invalid number of dimension");
assert(inf.domain() == sup.domain());
assert(inf.domain() == ord.domain());
assert(inf.domain().has(pstart));
using connectivity_t = std::conditional_t<P::ndim == 2, mln::experimental::c4_t, mln::experimental::c6_t>;
enum
{
UNPROCESSED = -1,
PROCESSED = 0
};
// Ord will hold the depth of the pixel during
// the propagation, i.e. when it is poped from the queue
connectivity_t nbh;
if (!mln::extension::fit(ord, nbh))
throw std::runtime_error("Image extension is not wide enough");
mln::fill(ord, UNPROCESSED);
mln::extension::try_fill(ord, PROCESSED);
// if (compute_indexes)
// sorted_indexes->reserve(ord.domain().size());
pset<I> queue(inf);
auto p = pstart;
V previous_level = inf(p);
queue.insert(previous_level, p);
ord(p) = 0;
// if (compute_indexes)
// sorted_indexes->push_back(p);
int depth = 0;
while (!queue.empty())
{
V cur_level;
P p;
auto tmp = queue.try_pop(previous_level);
if (tmp.has_value())
{
cur_level = previous_level;
p = tmp.value();
}
else
{
std::tie(cur_level, p) = queue.pop(previous_level);
++depth;
}
ord(p) = depth;
// if (compute_indexes)
// sorted_indexes->push_back(p);
for (auto q : nbh(p))
{
if (ord.at(q) == UNPROCESSED)
{
auto m = inf(q);
auto M = sup(q);
if (M < cur_level)
queue.insert(M, q);
else if (cur_level < m)
queue.insert(m, q);
else
queue.insert(cur_level, q);
ord(q) = PROCESSED;
}
}
previous_level = cur_level;
}
max_depth = depth;
}
} // namespace mln::moprho::experimental::details
#pragma once
#include <mln/morpho/experimental/private/hlinked_lists.hpp>
#include <limits>
#include <optional>
namespace mln::morpho::experimental::details
{
/// \brief A point set allows to keep a ordered set of points or indexes
/// sorted by a given key (generally the pixel value)
/// For now, it only supports indexes
///
/// \warning For two entries with identical keys (levels), this data structure acts
/// like a LIFO (last in first out)
template <class I>
class pset
{
public:
using Key = image_value_t<I>;
using Point = image_point_t<I>;
template <class J>
pset(J&& ima);
/// \brief Insert a point \p p at the given \p level
void insert(int level, Point p) noexcept;
/// \brief Return and remove the closest level pair (key,point) such that:
/// * there is a pair (l₁,p₁) in the set with l₁ ≥ level then (l₁,p₁)
/// * there is a pair (l₁,p₂) in the set with l₂ < level then (l₂,p₂)
std::pair<int, Point> pop(int level) noexcept;
/// \brief Return and remove the point if the queue is not empty
std::optional<Point> try_pop(int level) noexcept;
/// \brief Return true if the set is empty
bool empty() const noexcept;
private:
static_assert(std::is_integral_v<Key>);
static_assert(!std::numeric_limits<Key>::is_signed);
static_assert(std::numeric_limits<Key>::digits <= 16);
static constexpr int nLevel = 1 << std::numeric_limits<Key>::digits;
using link_image_t = image_ch_value_t<I, Point>;
using delegate_t = hlinked_lists<nLevel, Point, link_image_t>;
delegate_t m_delegate;
std::size_t m_size;
};
/******************************************/
/**** Implementation ****/
/******************************************/
template <class I>
template <class J>
pset<I>::pset(J&& ima)
: m_delegate{std::forward<J>(ima)},
m_size{0}
{
}
template <class I>
void pset<I>::insert(int level, Point p) noexcept
{
m_delegate.push_front(level, p);
m_size++;
}
template <class I>
auto pset<I>::pop(int level) noexcept -> std::pair<int, Point>
{
assert(m_size > 0);
int l = m_delegate.lower_bound(level);
if (l == nLevel)
{
l = m_delegate.upper_bound(level);
assert(l >= 0);
}
m_size--;
return {l, m_delegate.pop_front(l)};
}
template <class I>
auto pset<I>::try_pop(int level) noexcept -> std::optional<Point>
{
if (m_delegate.empty(level))
return {};
m_size--;
return m_delegate.pop_front(level);
}
template <class I>
bool pset<I>::empty() const noexcept
{
return m_size == 0;
}
} // namespace mln::morpho::experimental::details
#pragma once
#include <mln/morpho/experimental/maxtree.hpp>
#include <mln/morpho/experimental/private/immersion.hpp>
#include <mln/morpho/experimental/private/propagation.hpp>
#include <mln/core/image/view/cast.hpp>
namespace mln::morpho::experimental
{
enum {
ToS_NodeMapOriginalSize,
ToS_NodeMapTwiceSize
};
/// \brief Compute the tree of shapes.
/// \param[in] input Input image
/// \param[in] start_point Root point
/// \return A pair (tree, node_map) encoding the tree of shapes and the mapping (pixel -> node_index)
template <class I>
[[gnu::noinline]] auto tos(I input, image_point_t<I> pstart, int processing_flags = ToS_NodeMapTwiceSize);
/******************************************/
/**** Implementation ****/
/******************************************/
template <class I>
auto tos(I input, image_point_t<I> pstart, int /* processing_flags */)
{
static_assert(mln::is_a<I, mln::experimental::Image>());
using Domain = image_domain_t<I>;
using P = image_point_t<I>;
static_assert(std::is_same_v<Domain, mln::experimental::box2d> || std::is_same_v<Domain, mln::experimental::box3d>,
"Only 2D or 3D regular domain supported");
using connectivity_t = std::conditional_t<P::ndim == 2, mln::experimental::c4_t, mln::experimental::c6_t>;
connectivity_t nbh;
int max_depth;
auto [inf, sup] = details::immersion(input);
image_ch_value_t<I, int> ord = imchvalue<int>(inf).adjust(nbh);
details::propagation(inf, sup, ord, pstart, max_depth);
if (max_depth >= (1 << 16))
throw std::runtime_error("The ToS is too deep (too many number of levels)");
auto casted = mln::view::cast<uint16_t>(ord);
return maxtree(casted, nbh);
}
} // namespace mln::morpho::experimental
#include <mln/core/image/experimental/ndimage.hpp>
#include <mln/morpho/experimental/private/immersion.hpp>
#include <mln/morpho/experimental/private/propagation.hpp>
#include <mln/morpho/experimental/tos.hpp>
#include <mln/core/image/experimental/ndimage.hpp>
#include <mln/core/image/view/operators.hpp>
#include <mln/core/algorithm/accumulate.hpp>
#include <mln/accu/accumulators/max.hpp>
#include <fixtures/ImageCompare/image_compare.hpp>
#include <gtest/gtest.h>
TEST(ToSImmersion, twodimensional)
{
const mln::experimental::image2d<uint8_t> f = {{0, 1, 2}, //
......@@ -85,3 +93,159 @@ TEST(ToSImmersion, threedimensional)
ASSERT_IMAGES_EQ_EXP(ref_inf, inf);
ASSERT_IMAGES_EQ_EXP(ref_sup, sup);
}
template <class I, class J>
void test_propagation(I& f, J& ref)
{
//std::vector<P> indexes;
int max_depth;
auto [inf, sup] = mln::morpho::experimental::details::immersion(f);
mln::image_ch_value_t<I, int> ord(inf.domain());
mln::morpho::experimental::details::propagation(inf, sup, ord, inf.domain().tl(), max_depth);
int expected_max_depth = mln::accumulate(ref, mln::accu::features::max<>());
ASSERT_IMAGES_EQ_EXP(ord, ref);
ASSERT_EQ(expected_max_depth, max_depth);
//ASSERT_TRUE(std::is_sorted(indexes.begin(), indexes.end(), [&ord](P i, P j) { return ord[i] < ord[j]; }));
}
template <class T>
class ToSPropagation : public ::testing::Test
{
};
typedef ::testing::Types<mln::uint8> TestValueTypes;
TYPED_TEST_CASE(ToSPropagation, TestValueTypes);
TYPED_TEST(ToSPropagation, saddle_point)
{
const mln::experimental::image2d<TypeParam> f = {{0, 0, 0, 0}, //
{0, 2, 0, 0}, //
{0, 0, 2, 0}, //
{0, 0, 0, 0}};
const mln::experimental::image2d<int> ref = {{0, 0, 0, 0, 0, 0, 0}, //
{0, 0, 0, 0, 0, 0, 0}, //
{0, 0, 1, 0, 0, 0, 0}, //
{0, 0, 0, 0, 0, 0, 0}, //
{0, 0, 0, 0, 1, 0, 0}, //
{0, 0, 0, 0, 0, 0, 0}, //
{0, 0, 0, 0, 0, 0, 0}};
test_propagation(f, ref);
}
TYPED_TEST(ToSPropagation, two_branches_with_parallel_flooding)
{
const mln::experimental::image2d<TypeParam> f = {{0, 0, 0, 0}, //
{0, 2, 3, 0}, //
{0, 0, 2, 0}, //
{1, 2, 0, 0}};
const mln::experimental::image2d<int> ref = {{0, 0, 0, 0, 0, 0, 0}, //
{0, 0, 0, 0, 0, 0, 0}, //
{0, 0, 2, 2, 3, 0, 0}, //
{0, 0, 0, 0, 2, 0, 0}, //
{0, 0, 0, 0, 2, 0, 0}, //
{0, 0, 0, 0, 0, 0, 0}, //
{1, 1, 2, 0, 0, 0, 0}};
test_propagation(f, ref);
}
TYPED_TEST(ToSPropagation, two_branches_with_parallel_flooding_3d)
{
const mln::experimental::image3d<TypeParam> f = {{{0, 0, 0}, //
{0, 2, 3}, //
{0, 0, 2}}, //
{{1, 0, 1}, //
{0, 2, 1}, //
{1, 0, 2}}};
const mln::experimental::image3d<int> ref = {{{0, 0, 0, 0, 0}, //
{0, 0, 0, 0, 0}, //
{0, 0, 2, 2, 3}, //
{0, 0, 0, 0, 2}, //
{0, 0, 0, 0, 2}}, //
{{0, 0, 0, 0, 0}, //
{0, 0, 0, 0, 0}, //
{0, 0, 2, 1, 1}, //
{0, 0, 0, 0, 1}, //
{0, 0, 0, 0, 2}}, //
{{1, 0, 0, 0, 1}, //
{0, 0, 0, 0, 1}, //
{0, 0, 2, 1, 1}, //
{0, 0, 0, 0, 1}, //
{1, 0, 0, 0, 2}}};
test_propagation(f, ref);
}
TYPED_TEST(ToSPropagation, maxbounds)
{
const mln::experimental::image2d<TypeParam> f = {{0, 255, 0, 255}, //
{0, 255, 0, 0}, //
{255, 255, 255, 255}, //
{0, 255, 255, 0}};
const mln::experimental::image2d<int> ref = {{0, 0, 1, 1, 2, 2, 3}, //
{0, 0, 1, 1, 2, 2, 2}, //
{0, 0, 1, 1, 2, 2, 2}, //
{0, 0, 1, 1, 1, 1, 1}, //
{1, 1, 1, 1, 1, 1, 1}, //
{1, 1, 1, 1, 1, 1, 1}, //
{2, 1, 1, 1, 1, 1, 2}};
test_propagation(f, ref);
}
TYPED_TEST(ToSPropagation, chessboard)
{
const mln::experimental::image2d<TypeParam> f = {{0, 255, 0, 255}, //
{255, 0, 255, 0}, //
{0, 255, 0, 255}, //
{255, 0, 255, 0}};
const mln::experimental::image2d<int> ref = {{0, 0, 1, 0, 0, 0, 1}, //
{0, 0, 0, 0, 0, 0, 0}, //
{1, 0, 0, 0, 1, 0, 0}, //
{0, 0, 0, 0, 0, 0, 0}, //
{0, 0, 1, 0, 0, 0, 1}, //
{0, 0, 0, 0, 0, 0, 0}, //
{1, 0, 0, 0, 1, 0, 0}};
test_propagation(f, ref);
}
TEST(ToSConstruction, saddle_point)
{
const mln::experimental::image2d<uint8_t> f = {{0, 0, 0, 0}, //
{0, 2, 0, 0}, //
{0, 0, 2, 0}, //
{0, 0, 0, 0}};
const mln::experimental::point2d a = {0, 0};
const mln::experimental::point2d b = {2, 2};
const mln::experimental::point2d c = {4, 4};
const mln::experimental::image2d<mln::experimental::point2d> ref_roots = {{a, a, a, a, a, a, a}, //
{a, a, a, a, a, a, a}, //
{a, a, b, a, a, a, a}, //
{a, a, a, a, a, a, a}, //
{a, a, a, a, c, a, a}, //
{a, a, a, a, a, a, a}, //
{a, a, a, a, a, a, a}};
const mln::experimental::image2d<mln::experimental::point2d> ref_parent = {{a, a, a, a, a, a, a}, //
{a, a, a, a, a, a, a}, //
{a, a, a, a, a, a, a}, //
{a, a, a, a, a, a, a}, //
{a, a, a, a, a, a, a}, //
{a, a, a, a, a, a, a}, //
{a, a, a, a, a, a, a}};
mln::morpho::experimental::tos(f, {0,0});
}
Markdown is supported
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