Commit 3afcc66a by Edwin Carlinet

Added immersion ToS algorithm.

parent 3ce2c29e
 ... ... @@ -62,6 +62,7 @@ target_sources(Pylene PRIVATE src/io/imread.cpp src/io/imprint.cpp src/morpho/maxtree.cpp src/morpho/immersion.cpp src/morpho/component_tree.cpp src/morpho/hvector.cpp ) ... ...
 #pragma once #include #include #include #include #include namespace mln::morpho::experimental::details { /// \brief Immerse a 2d or 3d image into twice-as-big 2d image /// where the in-between pixel are intervals. /// /// Given 4 2D-pixels: /// a b /// c d /// The resulting image is the intervaled valued image: /// /// [a] [ (a∧b) (a∨b)] [b] /// [(a∧c) (a∨c)] [(a∧b∧c∧d) (a∨c∨d∨d)] [(d∧b) (d∨b)] /// [c] [ (c∧d) (c∨d)] [d] /// /// \param[in] ima The input 2D image /// \return An intervaled valued 2D image template std::pair, image_concrete_t> // immersion(I ima); /******************************************/ /**** Implementation ****/ /******************************************/ namespace impl { // Generic implementation template [[gnu::noinline]] // void line_immersion_T(I& f, P p, int xmin, int xmax, image_concrete_t& inf, image_concrete_t& sup) { if (xmin >= xmax) return; p.x() = xmin; auto q = 2 * p; auto a = f(p); inf(q) = a; sup(q) = a; for (int x = xmin + 1; x < xmax; ++x) { p.x() = x; auto prev = a; a = f(p); auto [m, M] = std::minmax(a, prev); q.x() = 2 * x - 1; inf(q) = m; sup(q) = M; q.x() = 2 * x; inf(q) = a; sup(q) = a; } } template [[gnu::noinline]] // void line_interpolation2d_T(int y, int xmin, int xmax, I& inf, I& sup) { for (int x = xmin; x < xmax; ++x) { inf({x, y}) = std::min(inf({x, y - 1}), inf({x, y + 1})); sup({x, y}) = std::max(sup({x, y - 1}), sup({x, y + 1})); } } // 2D generic version of immersion template void immersion_T(I& f, const mln::experimental::box2d& roi, image_concrete_t& inf, image_concrete_t& sup) { if (roi.empty()) return; int xmin = roi.tl().x(); int ymin = roi.tl().y(); int xmax = roi.br().x(); int ymax = roi.br().y(); // First line line_immersion_T(f, mln::experimental::point2d{0, ymin}, xmin, xmax, inf, sup); for (int y = ymin + 1; y < ymax; ++y) { line_immersion_T(f, mln::experimental::point2d{0, y}, xmin, xmax, inf, sup); line_interpolation2d_T(2 * y - 1, xmin * 2, xmax * 2 - 1, inf, sup); } } // Dispatch // \{ // For 2d-buffer images template void immersion(mln::experimental::image2d& f, mln::experimental::box2d, mln::experimental::image2d& inf, mln::experimental::image2d& sup) { mln_entering("mln::morpho::experimental::details::immersion (2d-buffer)") immersion_impl_table_t impl; immersion_ndimage(f, inf, sup, &impl); } // For 3d-buffer images (To be implemented) template void immersion(mln::experimental::image3d& f, mln::experimental::box3d, mln::experimental::image3d& inf, mln::experimental::image3d& sup) { mln_entering("mln::morpho::experimental::details::immersion (3d-buffer)") immersion_impl_table_t impl; immersion_ndimage(f, inf, sup, &impl); } // Fallback for 2d-images template void immersion(mln::experimental::Image& f, mln::experimental::box2d roi, image_concrete_t& inf, image_concrete_t& sup) { mln_entering("mln::morpho::experimental::details::immersion (generic)") immersion_T(static_cast(f), roi, inf, sup); } // Fallback for 3d-images (To be implemented) template void immersion(mln::experimental::Image& f, mln::experimental::box3d roi, image_concrete_t& inf, image_concrete_t& sup); // \} } // namespace impl template std::pair, image_concrete_t> // immersion(I ima) { static_assert(mln::is_a()); static_assert(std::is_same_v, mln::experimental::box2d> || std::is_same_v, mln::experimental::box3d>, "Input domain must be a box2d or a box3d"); int dim = image_point_t::ndim; // Compute the extended domain auto dom = ima.domain(); for (int d = 0; d < dim; ++d) { dom.tl()[d] = dom.tl()[d] * 2; dom.br()[d] = dom.br()[d] * 2 - 1; } image_concrete_t inf(dom); image_concrete_t sup(dom); impl::immersion(ima, ima.domain(), inf, sup); return { std::move(inf), std::move(sup) }; } }
 #pragma once #include #include #include #include namespace mln::morpho::experimental::details { /// Immersion algorithm for raw-contiguous buffer /// /// Given an input buffer of size \c n /// [a, b, c, d, e], /// it computes: /// inf = [a; min(a,b); b; min(b,c); c; min(c,d); d; min(d,e); e] /// sup = [a; max(a,b); b; max(b,c); c; max(c,d); d; max(d,e); e] template void buffer_immersion(const T* __restrict input, std::size_t n, T* __restrict inf, T* __restrict sup); /// Interpolation algorithm for raw-contiguous buffer. /// Compute the PW min/max between two buffer /// /// Given two input buffers /// inf = [a; min(a,b); b; min(b,c); c; min(c,d); d; min(d,e); e] /// sup = [a; max(a,b); b; max(b,c); c; max(c,d); d; max(d,e); e] template void buffer_interpolation_max(const T* __restrict A, const T* __restrict B, std::size_t n, T* __restrict out); template void buffer_interpolation_min(const T* __restrict A, const T* __restrict B, std::size_t n, T* __restrict out); struct immersion_impl_table_base_t { virtual void immersion(void* input, std::size_t n, void* inf, void* sup) = 0; virtual void interpolation_max(void* A, void* B, std::size_t n, void* out) = 0; virtual void interpolation_min(void* A, void* B, std::size_t n, void* out) = 0; }; /// Type-erased immersion algorithm for any buffer-encoded image void immersion_ndimage(ndbuffer_image& input, ndbuffer_image& inf, ndbuffer_image& sup, immersion_impl_table_base_t* impl); template struct immersion_impl_table_t : immersion_impl_table_base_t { void immersion(void* input, std::size_t n, void* inf, void* sup) final { buffer_immersion(static_cast(input), n, static_cast(inf), static_cast(sup)); } void interpolation_max(void* A, void* B, std::size_t n, void* out) final { buffer_interpolation_max(static_cast(A), static_cast(B), n, static_cast(out)); } void interpolation_min(void* A, void* B, std::size_t n, void* out) final { buffer_interpolation_min(static_cast(A), static_cast(B), n, static_cast(out)); } }; /******************************************/ /**** Implementation ****/ /******************************************/ template void buffer_immersion(const T* __restrict f, std::size_t n, T* __restrict inf, T* __restrict sup) { if (n == 0) return; inf[0] = f[0]; sup[0] = f[0]; for (std::size_t x = 1; x < n; ++x) { auto [m, M] = std::minmax(f[x-1], f[x]); inf[2 * x - 1] = m; sup[2 * x - 1] = M; inf[2 * x] = f[x]; sup[2 * x] = f[x]; } } template void buffer_interpolation_max(const T* __restrict A, const T* __restrict B, std::size_t n, T* __restrict out) { for (std::size_t i = 0; i < n; ++i) out[i] = std::max(A[i], B[i]); } template void buffer_interpolation_min(const T* __restrict A, const T* __restrict B, std::size_t n, T* __restrict out) { for (std::size_t i = 0; i < n; ++i) out[i] = std::min(A[i], B[i]); } }
 #include namespace mln::morpho::experimental::details { namespace { void immersion_image2d(std::byte* i_buffer, // std::byte* inf_buffer, // std::byte* sup_buffer, // std::size_t width, // std::size_t height, // std::ptrdiff_t i_stride, // std::ptrdiff_t inf_stride, // std::ptrdiff_t sup_stride, // immersion_impl_table_base_t* impl) { impl->immersion(i_buffer, width, inf_buffer, sup_buffer); for (std::size_t y = 1; y < height; ++y) { impl->immersion(i_buffer + y * i_stride, width, inf_buffer + 2 * y * inf_stride, sup_buffer + 2 * y * sup_stride); impl->interpolation_min(inf_buffer + (2 * (y - 1)) * inf_stride, inf_buffer + (2 * y) * inf_stride, 2 * width - 1, inf_buffer + (2 * y - 1) * inf_stride); impl->interpolation_max(sup_buffer + (2 * (y - 1)) * sup_stride, sup_buffer + (2 * y) * sup_stride, 2 * width - 1, sup_buffer + (2 * y - 1) * sup_stride); } } void interpolate_image2d_min(std::byte* buffer, std::size_t width, std::size_t height, std::ptrdiff_t line_stride, std::ptrdiff_t slice_stride, immersion_impl_table_base_t* impl) { std::byte* prev = buffer - slice_stride; std::byte* next = buffer + slice_stride; for (std::size_t y = 0; y < height; ++y) impl->interpolation_min(prev + y * line_stride, next + y * line_stride, width, buffer + y * line_stride); } void interpolate_image2d_max(std::byte* buffer, std::size_t width, std::size_t height, std::ptrdiff_t line_stride, std::ptrdiff_t slice_stride, immersion_impl_table_base_t* impl) { std::byte* prev = buffer - slice_stride; std::byte* next = buffer + slice_stride; for (std::size_t y = 0; y < height; ++y) impl->interpolation_max(prev + y * line_stride, next + y * line_stride, width, buffer + y * line_stride); } } // namespace void immersion_ndimage(ndbuffer_image& input, ndbuffer_image& inf, ndbuffer_image& sup, immersion_impl_table_base_t* impl) { [[maybe_unused]] int pdim = input.pdim(); assert(pdim == 2 || pdim == 3); assert(input.pdim() == pdim); assert(inf.pdim() == pdim); assert(sup.pdim() == pdim); const int depth = input.depth(); const int height = input.height(); const int width = input.width(); std::byte* i_buffer = input.buffer(); std::byte* inf_buffer = inf.buffer(); std::byte* sup_buffer = sup.buffer(); const std::ptrdiff_t i_stride = input.byte_stride(); const std::ptrdiff_t inf_stride = inf.byte_stride(); const std::ptrdiff_t sup_stride = sup.byte_stride(); immersion_image2d(i_buffer, inf_buffer, sup_buffer, width, height, i_stride, inf_stride, sup_stride, impl); for (int z = 1; z < depth; ++z) { immersion_image2d(i_buffer + z * input.byte_stride(2), // inf_buffer + (2 * z) * inf.byte_stride(2), // sup_buffer + (2 * z) * sup.byte_stride(2), // width, height, i_stride, inf_stride, sup_stride, impl); interpolate_image2d_min(inf_buffer + (2 * z - 1) * inf.byte_stride(2), 2 * width -1, 2 * height -1, inf_stride, inf.byte_stride(2), impl); interpolate_image2d_max(sup_buffer + (2 * z - 1) * sup.byte_stride(2), 2 * width -1, 2 * height -1, sup_stride, sup.byte_stride(2), impl); } } }
 #include #include #include #include #include #include #include #include #include #include #include #include "tos_tests_helper.hpp" #include #include #include using mln::uint8; template void test_propagation(const mln::Image& f, const mln::Image& ref) TEST(ToSImmersion, twodimensional) { using P = typename I::size_type; std::vector

indexes; int max_depth; auto g = mln::morpho::ToS::impl::immersion(f); mln_point(I) pmin = 0; auto ord = mln::morpho::ToS::impl::propagation(g, pmin, max_depth, &indexes); auto ima_ref = exact(ref); int expected_max_depth = mln::accumulate(ima_ref, mln::accu::features::max<>()); const mln::experimental::image2d f = {{0, 1, 2}, // {3, 0, 1}}; ASSERT_IMAGES_EQ(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]; })); } const mln::experimental::image2d ref_inf = {{0, 0, 1, 1, 2}, // {0, 0, 0, 0, 1}, // {3, 0, 0, 0, 1}}; template void test_construction(const mln::Image& f, const mln::Image& ref_parent, const mln::Image& ref_roots) { auto tree = mln::morpho::tos(f); const mln::experimental::image2d ref_sup = {{0, 1, 1, 2, 2}, // {3, 3, 1, 2, 2}, // {3, 3, 0, 1, 1}}; mln_ch_value(J, mln_point(J)) roots; auto par = tree_as_parent_image(tree, roots); ASSERT_IMAGES_EQ(par, ref_parent); ASSERT_IMAGES_EQ(roots, ref_roots); auto [inf, sup] = mln::morpho::experimental::details::immersion(f); ASSERT_IMAGES_EQ_EXP(ref_inf, inf); ASSERT_IMAGES_EQ_EXP(ref_sup, sup); } TEST(ToSImmersion, twodimensional) TEST(ToSImmersion, twodimensional_generic) { const mln::image2d f = {{0, 1, 2}, {3, 0, 1}}; using namespace mln::view::ops; const mln::experimental::image2d f = {{0, 1, 2}, // {3, 0, 1}}; auto fprime = f * uint8_t(1); const mln::experimental::image2d ref_inf = {{0, 0, 1, 1, 2}, // {0, 0, 0, 0, 1}, // {3, 0, 0, 0, 1}}; using R = mln::morpho::ToS::impl::irange; const mln::experimental::image2d ref_sup = {{0, 1, 1, 2, 2}, // {3, 3, 1, 2, 2}, // {3, 3, 0, 1, 1}}; const mln::image2d ref = {{{0, 0}, {0, 1}, {1, 1}, {1, 2}, {2, 2}}, {{0, 3}, {0, 3}, {0, 1}, {0, 2}, {1, 2}}, {{3, 3}, {0, 3}, {0, 0}, {0, 1}, {1, 1}}}; auto g = mln::morpho::ToS::impl::immersion(f); ASSERT_IMAGES_EQ(ref, g); auto [inf, sup] = mln::morpho::experimental::details::immersion(fprime); ASSERT_IMAGES_EQ_EXP(ref_inf, inf); ASSERT_IMAGES_EQ_EXP(ref_sup, sup); } TEST(ToSImmersion, threedimensional) { const mln::image3d f = { {{0, 1, 2}, {3, 0, 1}}, const mln::experimental::image3d f = { {{0, 1, 2}, {3, 0, 1}}, // {{2, 0, 4}, {0, 0, 2}}, }; using R = mln::morpho::ToS::impl::irange; const mln::image3d ref = { {{{0, 0}, {0, 1}, {1, 1}, {1, 2}, {2, 2}}, {{0, 3}, {0, 3}, {0, 1}, {0, 2}, {1, 2}}, {{3, 3}, {0, 3}, {0, 0}, {0, 1}, {1, 1}}}, {{{0, 2}, {0, 2}, {0, 1}, {0, 4}, {2, 4}}, {{0, 3}, {0, 3}, {0, 1}, {0, 4}, {1, 4}}, {{0, 3}, {0, 3}, {0, 0}, {0, 2}, {1, 2}}}, {{{2, 2}, {0, 2}, {0, 0}, {0, 4}, {4, 4}}, {{0, 2}, {0, 2}, {0, 0}, {0, 4}, {2, 4}}, {{0, 0}, {0, 0}, {0, 0}, {0, 2}, {2, 2}}}, const mln::experimental::image3d ref_inf = { {{0, 0, 1, 1, 2}, // {0, 0, 0, 0, 1}, // {3, 0, 0, 0, 1}}, // {{0, 0, 0, 0, 2}, // {0, 0, 0, 0, 1}, // {0, 0, 0, 0, 1}}, // {{2, 0, 0, 0, 4}, // {0, 0, 0, 0, 2}, // {0, 0, 0, 0, 2}}, // }; auto g = mln::morpho::ToS::impl::immersion(f); ASSERT_IMAGES_EQ(ref, g); } template class ToSPropagation : public ::testing::Test { }; typedef ::testing::Types TestValueTypes; TYPED_TEST_CASE(ToSPropagation, TestValueTypes); TYPED_TEST(ToSPropagation, saddle_point) { const mln::image2d f = {{0, 0, 0, 0}, // {0, 2, 0, 0}, // {0, 0, 2, 0}, // {0, 0, 0, 0}}; const mln::image2d 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::image2d f = {{0, 0, 0, 0}, // {0, 2, 3, 0}, // {0, 0, 2, 0}, // {1, 2, 0, 0}}; const mln::image2d 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::image3d f = {{{0, 0, 0}, // {0, 2, 3}, // {0, 0, 2}}, // {{1, 0, 1}, // {0, 2, 1}, // {1, 0, 2}}}; const mln::image3d ref = {{{0, 0, 0, 0, 0}, // {0, 0, 0, 0, 0}, // {0, 0, 2, 2, 3}, // {0, 0, 0, 0, 2}, //