Commit 775acc37 authored by Edwin Carlinet's avatar Edwin Carlinet
Browse files

Merge branch 'development/watershed' into 'development/ranges'

migrate watershed implementation

See merge request !76
parents f52d1a7c 7374639f
Pipeline #13799 failed with stages
in 91 minutes and 51 seconds
......@@ -15,6 +15,7 @@
#include <mln/morpho/experimental/median_filter.hpp>
#include <mln/morpho/experimental/opening.hpp>
#include <mln/morpho/experimental/reconstruction.hpp>
#include <mln/morpho/experimental/watershed.hpp>
#include <mln/labeling/experimental/local_extrema.hpp>
......@@ -216,4 +217,16 @@ BENCHMARK_F(BMMorpho, minima)(benchmark::State& st)
}
BENCHMARK_F(BMMorpho, watershed)(benchmark::State& st)
{
auto f = [](const image_t& input, image_t&) {
int nlabel;
mln::morpho::experimental::watershed<int8_t>(input, mln::experimental::c4, nlabel);
return nlabel;
};
this->run(st, f);
}
BENCHMARK_MAIN();
......@@ -15,6 +15,7 @@
#include <mln/morpho/structural/dilate.hpp>
#include <mln/morpho/structural/erode.hpp>
#include <mln/morpho/structural/opening.hpp>
#include <mln/morpho/watershed.hpp>
#include <mln/labeling/local_extrema.hpp>
......@@ -162,5 +163,15 @@ BENCHMARK_F(BMMorpho, minima)(benchmark::State& st)
this->run(st, f);
}
BENCHMARK_F(BMMorpho, watershed)(benchmark::State& st)
{
auto f = [](const image_t& input, image_t&) {
int nlabel;
mln::morpho::watershed<int8_t>(input, mln::c4, nlabel);
return nlabel;
};
this->run(st, f);
}
BENCHMARK_MAIN();
......@@ -614,9 +614,9 @@ namespace mln
for (int k = 0; k < N; ++k)
{
coords[k] = this->__axes(k).domain_begin;
shp[k] = static_cast<short>(this->__axes(k).domain_end - this->__axes(k).domain_begin);
strides[k] = this->__axes(k).byte_stride / sizeof(T);
coords[k] = this->__axes(k).domain_begin;
shp[N - k - 1] = static_cast<short>(this->__axes(k).domain_end - this->__axes(k).domain_begin);
strides[N - k - 1] = this->__axes(k).byte_stride;
}
auto ptr = Impl::get_pointer(this->__info(), coords);
......
......@@ -52,7 +52,7 @@ namespace mln::io::experimental
template <class I>
struct impl_t : impl_base_t
{
impl_t(I& x) { this->m_ima = &x; }
impl_t(I& x) { this->m_ima = (void*)(&x); }
~impl_t() final = default;
std::size_t formatted_size(P p) const final { return fmt::formatted_size("{}", at(p)); }
void print(P p, int width) const final { fmt::print("{:>{}d}", at(p), width); }
......
#pragma once
#include <mln/core/concept/new/images.hpp>
#include <mln/morpho/experimental/private/pqueue_hqueue_fifo.hpp>
namespace mln::morpho::experimental::details
{
/// Provides a priority queue of points with a the fifo property
template <class I>
class pqueue_fifo
{
using key_type = image_value_t<I>;
using value_type = image_point_t<I>;
public:
template <class J>
pqueue_fifo(J&& f);
void push(const key_type& priority, const value_type& element);
void pop();
std::pair<key_type, value_type> top() const;
bool empty() const;
private:
static_assert(value_traits<key_type>::quant <= 16, "Only low quantized type supported.");
pqueue_hqueue_fifo<image_concrete_t<I>> m_delegate;
};
/******************************************/
/**** Implementation ****/
/******************************************/
template <class I>
template <class J>
inline pqueue_fifo<I>::pqueue_fifo(J&& f)
: m_delegate{std::forward<J>(f)}
{
}
template <class I>
inline void pqueue_fifo<I>::push(const key_type& k, const value_type& v)
{
m_delegate.push(k, v);
}
template <class I>
inline void pqueue_fifo<I>::pop()
{
m_delegate.pop();
}
template <class I>
inline bool pqueue_fifo<I>::empty() const
{
return m_delegate.empty();
}
template <class I>
inline auto pqueue_fifo<I>::top() const -> std::pair<key_type, value_type>
{
return m_delegate.top();
}
} // namespace mln::morpho::experimental::details
#pragma once
#include <mln/core/concept/new/images.hpp>
#include <array>
#include <limits>
#include <type_traits>
namespace mln::morpho::experimental::details
{
template <class I>
class pqueue_hqueue_fifo
{
using key_type = image_value_t<I>;
using value_type = image_point_t<I>;
static_assert(std::is_integral_v<key_type>);
static_assert(!std::is_signed<key_type>::value, "Key type must be unsigned.");
static_assert(sizeof(key_type) <= 2, "Only low-quantized value supported.");
static_assert(sizeof(key_type) < sizeof(int), "Unsupported architecure.");
public:
template <class J>
pqueue_hqueue_fifo(J&& f);
void push(const key_type& k, const value_type& v);
void pop();
std::pair<key_type, value_type> top() const;
bool empty() const;
private:
static constexpr int nvalues = std::numeric_limits<key_type>::max() + 1;
struct node_t
{
int size = 0;
value_type head;
value_type tail;
};
// Pointer to head and tail element in the queue
std::array<node_t, nvalues> m_queues;
// Next pointer image
image_ch_value_t<I, value_type> m_next;
int m_current_level;
};
template <class I>
template <class J>
inline pqueue_hqueue_fifo<I>::pqueue_hqueue_fifo(J&& f)
: m_next(std::forward<J>(f), image_build_params{})
, m_current_level{std::numeric_limits<key_type>::max()}
{
}
template <class I>
inline void pqueue_hqueue_fifo<I>::pop()
{
mln_precondition(m_queues[m_current_level].size > 0);
auto& head = m_queues[m_current_level].head;
head = m_next(head);
m_queues[m_current_level].size--;
if (m_queues[m_current_level].size > 0)
return;
// Try to go up
{
do {
m_current_level++;
} while (m_current_level < std::numeric_limits<key_type>::max() && m_queues[m_current_level].size == 0);
}
}
template <class I>
inline void pqueue_hqueue_fifo<I>::push(const key_type& k, const value_type& v)
{
if (m_queues[k].size == 0)
{
m_queues[k].head = v;
m_queues[k].tail = v;
}
else
{
m_next(m_queues[k].tail) = v;
m_queues[k].tail = v;
}
m_queues[k].size++;
if (k < m_current_level)
m_current_level = k;
}
template <class I>
inline auto pqueue_hqueue_fifo<I>::top() const -> std::pair<key_type, value_type>
{
mln_precondition(m_queues[m_current_level].size > 0);
return std::make_pair(static_cast<key_type>(m_current_level), m_queues[m_current_level].head);
}
template <class I>
inline bool pqueue_hqueue_fifo<I>::empty() const
{
return m_queues[m_current_level].size == 0;
}
} // namespace mln::morpho::experimental::details
#pragma once
#include <mln/core/algorithm/for_each.hpp>
#include <mln/core/extension/fill.hpp>
#include <mln/core/image/image.hpp>
#include <mln/core/image/view/value_extended.hpp>
#include <mln/core/neighborhood/neighborhood.hpp>
#include <mln/core/trace.hpp>
#include <mln/labeling/experimental/local_extrema.hpp>
#include <mln/morpho/experimental/private/pqueue.hpp>
namespace mln::morpho::experimental
{
template <class Label_t, class InputImage, class Neighborhood>
image_ch_value_t<std::remove_reference_t<InputImage>, Label_t> //
watershed(InputImage&& ima, Neighborhood&& nbh, int& nlabel);
/******************************************/
/**** Implementation ****/
/******************************************/
namespace impl
{
template <class I, class N, class O>
int watershed(I input, N nbh, O output)
{
using Label_t = image_value_t<O>;
// 1. Labelize minima (note that output is initialized to -1)
const int nlabel = mln::labeling::experimental::impl::local_minima(input, nbh, output, std::less<Label_t>());
constexpr int kUnlabeled = -2;
constexpr int kInqueue = -1;
constexpr int kWaterline = 0;
// 2. inset neighbors inqueue
// Pixels in the border gets the status 0 (deja vu)
// Pixels in the queue get -1
// Pixels not in the queue get -2
mln::morpho::experimental::details::pqueue_fifo<I> pqueue(input);
{
output.extension().fill(kWaterline);
mln_foreach_new(auto px, output.new_pixels())
{
// Not a local minimum => early exit
if (px.val() != 0)
continue;
bool is_local_min_neighbor = false;
for (auto nx : nbh(px))
if (nx.val() > 0) // p is neighbhor to a local minimum
{
is_local_min_neighbor = true;
break;
}
if (is_local_min_neighbor)
{
px.val() = kInqueue;
pqueue.push(input(px.point()), px.point());
}
else
{
px.val() = kUnlabeled;
}
}
}
// 3. flood from minima
{
while (!pqueue.empty())
{
auto [level, p] = pqueue.top();
auto pxOut = output.new_pixel(p);
mln_assertion(pxOut.val() == kInqueue);
pqueue.pop();
// Check if there is a single marked neighbor
Label_t common_label = kWaterline;
bool has_single_adjacent_marker = false;
for (auto nxOut : nbh(pxOut))
{
int nlbl = nxOut.val();
if (nlbl <= 0)
continue;
else if (common_label == kWaterline)
{
common_label = nlbl;
has_single_adjacent_marker = true;
}
else if (nlbl != common_label)
{
has_single_adjacent_marker = false;
break;
}
}
if (!has_single_adjacent_marker)
{
// If there are multiple labels => waterline
pxOut.val() = kWaterline;
}
else
{
// If a single label, it gets labeled
// Add neighbors in the queue
pxOut.val() = common_label;
for (auto q : nbh(p))
{
auto nlbl = output.at(q);
if (nlbl == kUnlabeled)
{
pqueue.push(input(q), q);
output(q) = kInqueue;
}
}
}
}
}
// 3. Label all unlabeled pixels
{
mln::for_each(output, [](auto& v) {
if (v < 0)
v = kWaterline;
});
}
return nlabel;
}
}
template <class Label_t, class InputImage, class Neighborhood>
image_ch_value_t<std::remove_reference_t<InputImage>, Label_t> //
watershed(InputImage&& ima, Neighborhood&& nbh, int& nlabel)
{
using I = std::remove_reference_t<InputImage>;
using N = std::remove_reference_t<Neighborhood>;
static_assert(mln::is_a<I, mln::experimental::Image>());
static_assert(mln::is_a<N, mln::experimental::Neighborhood>());
static_assert(std::is_integral<Label_t>::value, "The label type must integral.");
static_assert(std::is_signed<Label_t>::value, "The label type must be signed.");
mln_entering("mln::morpho::watershed");
constexpr Label_t kUninitialized = -1;
image_build_error_code err = IMAGE_BUILD_OK;
auto output = imchvalue<Label_t>(ima) //
.adjust(nbh)
.set_init_value(kUninitialized)
.get_status(&err)
.build();
if (err == IMAGE_BUILD_OK)
{
nlabel = impl::watershed(ima, nbh, output);
}
else
{
mln::trace::warn("[Performance] The extension is not wide enough");
auto out = view::value_extended(output, kUninitialized);
nlabel = impl::watershed(ima, nbh, out);
}
return output;
}
} // namespace mln::morpho::experimental
#include <mln/core/image/image2d.hpp>
#include <mln/core/neighb2d.hpp>
#include <mln/morpho/watershed.hpp>
#include <mln/morpho/experimental/watershed.hpp>
#include <mln/core/image/experimental/ndimage.hpp>
#include <mln/core/neighborhood/c4.hpp>
#include <fixtures/ImageCompare/image_compare.hpp>
......@@ -11,38 +12,38 @@ using namespace mln;
TEST(Morpho, watershed_gray)
{
const mln::image2d<mln::uint8> input = {{2, 2, 2, 2, 2}, //
{40, 30, 30, 30, 40}, //
{40, 20, 20, 20, 40}, //
{40, 40, 20, 40, 40}, //
{1, 5, 20, 5, 1}};
const mln::image2d<mln::int16> ref = {{1, 1, 1, 1, 1}, //
{1, 1, 1, 1, 1}, //
{0, 1, 1, 1, 0}, //
{2, 0, 1, 0, 3}, //
{2, 2, 0, 3, 3}};
const mln::experimental::image2d<uint8_t> input = {{2, 2, 2, 2, 2}, //
{40, 30, 30, 30, 40}, //
{40, 20, 20, 20, 40}, //
{40, 40, 20, 40, 40}, //
{1, 5, 20, 5, 1}};
const mln::experimental::image2d<int16_t> ref = {{1, 1, 1, 1, 1}, //
{1, 1, 1, 1, 1}, //
{0, 1, 1, 1, 0}, //
{2, 0, 1, 0, 3}, //
{2, 2, 0, 3, 3}};
int nlabel;
auto res = mln::morpho::watershed<mln::int8>(input, mln::c4, nlabel);
ASSERT_IMAGES_EQ(ref, res);
auto res = mln::morpho::experimental::watershed<int16_t>(input, mln::experimental::c4, nlabel);
ASSERT_IMAGES_EQ_EXP(ref, res);
}
TEST(Morpho, watershed_thick)
{
const mln::image2d<mln::uint8> input = {{0, 10, 0, 10, 0}, //
{0, 10, 10, 10, 10}, //
{10, 10, 10, 10, 10}, //
{0, 10, 10, 10, 10}, //
{0, 10, 0, 10, 0}};
const mln::image2d<mln::int16> ref = {{1, 0, 2, 0, 3}, //
{1, 1, 0, 3, 3}, //
{0, 0, 0, 0, 0}, //
{4, 4, 0, 6, 6}, //
{4, 0, 5, 0, 6}};
const mln::experimental::image2d<uint8_t> input = {{0, 10, 0, 10, 0}, //
{0, 10, 10, 10, 10}, //
{10, 10, 10, 10, 10}, //
{0, 10, 10, 10, 10}, //
{0, 10, 0, 10, 0}};
const mln::experimental::image2d<int16_t> ref = {{1, 0, 2, 0, 3}, //
{1, 1, 0, 3, 3}, //
{0, 0, 0, 0, 0}, //
{4, 4, 0, 6, 6}, //
{4, 0, 5, 0, 6}};
int nlabel;
auto res = mln::morpho::watershed<mln::int16>(input, mln::c4, nlabel);
ASSERT_IMAGES_EQ(ref, res);
auto res = mln::morpho::experimental::watershed<int16_t>(input, mln::experimental::c4, nlabel);
ASSERT_IMAGES_EQ_EXP(ref, res);
}
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