Commit 6c7f1e58 authored by Edwin Carlinet's avatar Edwin Carlinet
Browse files

Implement disc decomposition by periodic lines.

parent 7dc04336
#ifndef MLN_CORE_SE_BALL_HPP
#define MLN_CORE_SE_BALL_HPP
#include <mln/core/neighborhood/dyn_neighborhood.hpp>
#include <vector>
namespace mln
{
namespace se
{
/// Ball structuring element.
///
/// A ball of radius \$r\$ stands is centered in (0,0) and has all the
/// points whose (euclidean) distance is lower than \$r\$
struct ball2d;
/// \brief Helper function to make a Ball
ball2d make_ball2d(float r);
/******************************************/
/**** Implementation ****/
/******************************************/
struct ball2d : dyn_neighborhood_base<dynamic_neighborhood_tag, ball2d>
{
using is_incremental = std::true_type;
using dec_type = dyn_neighborhood<std::vector<point2d>, dynamic_neighborhood_tag>;
using inc_type = dyn_neighborhood<std::vector<point2d>, dynamic_neighborhood_tag>;
ball2d() = default;
private:
ball2d(const std::vector<point2d>& dpts, const std::vector<point2d>& vdec, const std::vector<point2d>& vinc)
: dpoints(dpts), m_dec(vdec), m_inc(vinc)
{
}
public:
static ball2d make(float r_)
{
typedef point2d::value_type P;
mln_precondition(r_ >= 0);
int k = (int)r_;
int d = 2 * k + 1;
float r2 = r_ * r_;
std::vector<point2d> vinc, vdec, dpoints;
dpoints.reserve(d * d);
vinc.reserve(d);
vdec.reserve(d);
for (int i = -k; i <= k; ++i)
for (int j = -k; j <= k; ++j)
if (sqr(i) + sqr(j) <= r2)
{
point2d p = {(P)i, (P)j};
dpoints.push_back(p);
}
{
int n = (int)dpoints.size();
for (int i = 0; i < n; i++)
{
point2d p = dpoints[i];
vdec.push_back(p + point2d{0, -1}); // before begin of the line
while (i + 1 < n and dpoints[i][0] == dpoints[i + 1][0])
++i;
vinc.push_back(dpoints[i]); // last element of the line
}
}
return ball2d(dpoints, vdec, vinc);
}
inc_type inc() const { return m_inc; }
dec_type dec() const { return m_dec; }
const auto& offsets() const { return dpoints; }
const std::vector<point2d> dpoints;
private:
dec_type m_dec;
inc_type m_inc;
};
inline ball2d make_ball2d(float r) { return ball2d::make(r); }
} // end of namespace mln::se
} // end of namespace mln
#endif //! MLN_CORE_SE_BALL_HPP
#ifndef MLN_CORE_SE_DISC_HPP
# define MLN_CORE_SE_DISC_HPP
# include <mln/core/neighborhood/neighborhood_base.hpp>
# include <mln/core/se/periodic_line2d.hpp>
# include <vector>
/// \file
namespace mln
{
namespace se
{
/// Create a disc of a given radius r with or without approximation.
///
///
/// The extent of the structuring will be 2*⌊r⌋+1. If an n-approximation is
/// given, the radial decomposition of the disc is used when possible. The
/// approximiation is a 2-n side polygon. If 0, the exact euclidean disc is
/// used, that are all points \f$ \{p \mid |p| \le r\} \f$
struct disc
#ifndef MLN_DOXYGEN
: dyn_neighborhood_base<dynamic_neighborhood_tag, disc>
#endif
{
using is_incremental = std::true_type;
using is_decomposable = std::true_type;
using dec_type = dyn_neighborhood<std::vector<point2d>, dynamic_neighborhood_tag>;
using inc_type = dyn_neighborhood<std::vector<point2d>, dynamic_neighborhood_tag>;
/// Constructor
///
/// \param radius The radius r of the disc.
/// \param approximation Must be 0 or 8 (default: 8)
disc(float radius, int approximation = 8);
/// \brief Return a range of SE for decomposition
std::vector<periodic_line2d> decompose() const;
/// \brief Return a range of SE offsets
std::vector<point2d> offsets() const;
/// \brief True if the SE is decomposable (i.e. constructed with no-approximation)
bool decomposable() const;
/// \brief A WNeighborhood to be added when used incrementally
inc_type inc() const;
/// \brief A WNeighborhood to be substracted when used incrementally
dec_type dec() const;
/// \brief Return the radius of the disc.
float radius() const { return m_radius; }
// protected to allow testing
static std::array<int, 3> _compute_decomposition_coeff(int radius, int nlines);
private:
float m_radius;
int m_nlines; // number of periodic lines for decomposition (0 for the euclidean disc)
};
/******************************************/
/**** Implementation ****/
/******************************************/
inline
disc::disc(float radius, int approximation)
: m_radius(radius),
m_nlines(approximation)
{
assert(m_radius >= 0);
assert(m_nlines == 0 || m_nlines == 2 || m_nlines == 4 || m_nlines == 8);
}
inline
bool disc::decomposable() const
{
return m_nlines > 0;
}
inline
std::array<int, 3> disc::_compute_decomposition_coeff(int radius, int nlines)
{
(void) nlines;
// We decompose a disc of radius r
// D = k₁ * Lθ₁ + k₂ * Lθ₂ + ... + kₙ * Lθₙ
// with k₁ = k₂ = kₙ
// Those are coefficient for periodic lines:
// k0: coeff for 0°, 90° lines
// k1: coeff for 45°, 135° lines
// k2: coeff for 30°, 60°, 120°, 150° lines
int k0 = 0;
int k1 = 0;
int k2 = 0;
if (radius > 0)
{
// Empirically, we have found the best k0, k1, k2 values
// that minimizes the error with the real euclidean disc.
// The linear approximation of k0(r) and k2(r) have a +- 1 error,
// they are close enough to use them. So:
// k0(r) = a0 * t + b0
// k2(r) = a2 * t + b2
const float a0 = 0.22498089194617754;
const float b0 = -0.74997778133276127;
const float a2 = 0.092868889904065763;
const float b2 = 0.03471904979442568;
k0 = (radius < 11) ? (1 + (radius+1) % 2) : static_cast<int>(std::round(a0 * radius + b0));
k2 = (radius < 7) ? 0 : static_cast<int>(std::round(a2 * radius + b2));
// Compute the missing value k1
// k0 += 1 => radial extent += 1
// k1 += 1 => radial extent += 2
// k2 += 1 => radial extent += 6
int r = radius - k0 - 6 * k2;
if (r > 0)
{
k1 = r / 2;
k0 += r % 2;
}
}
return {{k0, k1, k2}};
}
inline
std::vector<periodic_line2d> disc::decompose() const
{
mln_precondition(m_nlines > 0);
std::vector<periodic_line2d> lines;
lines.reserve(m_nlines);
std::array<int, 3> k = _compute_decomposition_coeff(static_cast<int>(m_radius), m_nlines);
const point2d se[] = {{0,1}, {1,0}, {1, 1}, {1, -1}, {1, 2}, {2, 1}, {2, -1}, {1, -2}};
lines.push_back(periodic_line2d(se[0], k[0]));
lines.push_back(periodic_line2d(se[1], k[0]));
if (k[1] > 0)
{
lines.push_back(periodic_line2d(se[2], k[1]));
lines.push_back(periodic_line2d(se[3], k[1]));
}
if (k[2] > 0)
{
lines.push_back(periodic_line2d(se[4], k[2]));
lines.push_back(periodic_line2d(se[5], k[2]));
lines.push_back(periodic_line2d(se[6], k[2]));
lines.push_back(periodic_line2d(se[7], k[2]));
}
return lines;
}
inline
std::vector<point2d> disc::offsets() const
{
typedef point2d::value_type P;
int r = static_cast<int>(m_radius);
int extent = 2 * r + 1;
float radius_sqr = m_radius * m_radius;
std::vector<point2d> dpoints;
dpoints.reserve(extent * extent);
for (int i = -r; i <= r; ++i)
for (int j = -r; j <= r; ++j)
if (i*i + j*j <= radius_sqr)
{
point2d p = {(P)i, (P)j};
dpoints.push_back(p);
}
return dpoints;
}
inline
disc::dec_type disc::dec() const
{
typedef point2d::value_type P;
const int r = static_cast<int>(m_radius);
const int extent = 2 * r + 1;
const float radius_sqr = m_radius * m_radius;
std::vector<point2d> vdec;
vdec.reserve(extent);
for (int y = -r; y <= r; ++y)
{
for (int x = -r; x <= 0; ++x)
if (y*y + x*x <= radius_sqr)
{
point2d p = {(P)y, (P)(x - 1)}; // before begin of the line
vdec.push_back(p);
break;
}
}
return dec_type(std::move(vdec));
}
inline
disc::inc_type disc::inc() const
{
typedef point2d::value_type P;
const int r = static_cast<int>(m_radius);
const int extent = 2 * r + 1;
const float radius_sqr = m_radius * m_radius;
std::vector<point2d> vinc;
vinc.reserve(extent);
for (int y = -r; y <= r; ++y)
{
for (int x = r; x >= 0; --x)
if (y*y + x*x <= radius_sqr)
{
point2d p = {(P)y, (P)x}; // last point of the line
vinc.push_back(p);
break;
}
}
return inc_type(std::move(vinc));
}
/*
//This the matlab decomposition, the approximation is worst than ours
inline
std::vector<periodic_line2d> disc::decompose() const
{
std::vector<periodic_line2d> lines;
lines.reserve(m_nlines);
const point2d se[] = {{0,1}, {1,0}, {1, 1}, {1, -1}, {1, 2}, {2, 1}, {2, -1}, {1, -2}};
// We decompose a disc of radius r
// D = k₁ * Lθ₁ + k₂ * Lθ₂ + ... + kₙ * Lθₙ
// with k₁ = k₂ = kₙ
// D is a 2n polygon with mininimal radius
// rmin = k.tan⁻¹(π/2n)
// rmax = k.sin⁻¹(π/2n)
float angle = M_PI / (2 * m_nlines);
float se_objective_radius = 2 * m_radius / (1.f / std::sin(angle) + 1.f / std::tan(angle));
int current_radius = 0;
for (auto dp : se)
{
// Compute k, the number of repetition of the se.
float se_elementary_length = std::hypot(dp[0], dp[1]);
int k = static_cast<int>(se_objective_radius / se_elementary_length);
if (k > 0)
{
int se_size = 2 * k + 1;
lines.push_back(periodic_line2d(-dp * k, dp, se_size));
current_radius += std::abs(dp[0]) * k;
}
}
// The current radius is too small
// We pad it with pure horizontal and vertical lines
{
int k = static_cast<int>(m_radius) - current_radius;
int se_size = 2 * k + 1;
lines.push_back(periodic_line2d(point2d{-1,0} * k, point2d{1,0}, se_size));
lines.push_back(periodic_line2d(point2d{0,-1} * k, point2d{0,1}, se_size));
}
return lines;
}
*/
} // end of namespace mln::se
} // end of namespace mln
#endif //!MLN_CORE_SE_DISC_HPP
......@@ -48,3 +48,6 @@ add_core_test(${test_prefix}paste algorithm/paste.cpp)
add_core_test(${test_prefix}equal algorithm/equal.cpp)
add_core_test(${test_prefix}accumulate algorithm/accumulate.cpp)
add_core_test(${test_prefix}sort_indexes algorithm/sort_indexes.cpp)
# test SE
add_core_test(${test_prefix}disc se/disc.cpp)
#include <mln/core/image/image2d.hpp>
#include <mln/core/se/disc.hpp>
#include <boost/format.hpp>
#include <gtest/gtest.h>
void naive_dilate(mln::image2d<bool>& f, const mln::se::periodic_line2d& se)
{
mln::image2d<bool> g;
mln::resize(g, f).init(false);
mln_pixter(px, f);
mln_pixter(pxout, g);
mln_iter(nxout, se(pxout));
mln_forall(px, pxout)
if (px->val())
{
mln_forall(nxout)
if (f.domain().has(nxout->point()))
nxout->val() = true;
}
f = std::move(g);
}
mln::image2d<bool> draw_ball_by_decomposition(float radius, int extent, int& computed_extent)
{
mln::image2d<bool> f(mln::box2d{{(short)(-extent), (short)(-extent)}, {(short)(extent+1), (short)(extent+1)}}, 3, false);
mln::se::disc ball(radius, 8);
auto ses = ball.decompose();
f.at(0,0) = true;
for (auto line : ses)
naive_dilate(f, line);
int k = 0;
for (int i = -extent; i <= extent; ++i)
k += f.at(0, i);
computed_extent = k;
return f;
}
float compute_disc_error(const mln::image2d<bool>& f, float radius)
{
int nerror = 0;
mln_pixter(px, f);
mln_foreach(auto px, f.pixels())
{
bool ref = l2norm(px.point()) <= radius;
nerror += (ref != px.val());
}
return nerror / radius;
}
// Best coef tuples (computed by brute force) from radius r=1 to 199
int best_decomposition_coeff[199][3] = {
{1, 0, 0}, {2, 0, 0}, {1, 1, 0}, {2, 1, 0}, {1, 2, 0}, {2, 2, 0}, {1, 0, 1}, {2, 3, 0},
{1, 1, 1}, {2, 1, 1}, {1, 2, 1}, {2, 2, 1}, {3, 2, 1}, {2, 0, 2}, {3, 0, 2}, {2, 1, 2},
{3, 1, 2}, {4, 4, 1}, {3, 2, 2}, {4, 2, 2}, {3, 0, 3}, {4, 3, 2}, {5, 3, 2}, {4, 1, 3},
{5, 1, 3}, {6, 4, 2}, {5, 2, 3}, {6, 5, 2}, {5, 3, 3}, {6, 3, 3}, {5, 1, 4}, {6, 4, 3},
{7, 4, 3}, {6, 2, 4}, {7, 2, 4}, {8, 5, 3}, {7, 3, 4}, {8, 6, 3}, {9, 6, 3}, {8, 4, 4},
{9, 7, 3}, {8, 5, 4}, {9, 5, 4}, {10, 5, 4}, {9, 3, 5}, {10, 6, 4}, {9, 4, 5}, {10, 7, 4},
{11, 7, 4}, {10, 5, 5}, {11, 5, 5}, {12, 8, 4}, {11, 6, 5}, {12, 6, 5}, {11, 4, 6}, {12, 7, 5},
{13, 7, 5}, {12, 5, 6}, {13, 8, 5}, {12, 6, 6}, {13, 6, 6}, {14, 9, 5}, {13, 7, 6}, {14, 7, 6},
{15, 10, 5}, {14, 8, 6}, {15, 8, 6}, {14, 6, 7}, {15, 9, 6}, {16, 9, 6}, {15, 7, 7}, {16, 10, 6},
{15, 8, 7}, {16, 8, 7}, {17, 8, 7}, {16, 9, 7}, {17, 9, 7}, {16, 7, 8}, {17, 10, 7}, {18, 10, 7},
{17, 8, 8}, {18, 11, 7}, {17, 9, 8}, {18, 9, 8}, {19, 9, 8}, {18, 10, 8}, {19, 10, 8}, {20, 10, 8},
{19, 11, 8}, {20, 11, 8}, {19, 9, 9}, {20, 9, 9}, {21, 12, 8}, {20, 10, 9}, {21, 10, 9}, {20, 11, 9},
{21, 11, 9}, {22, 11, 9}, {21, 12, 9}, {22, 12, 9}, {21, 10, 10}, {22, 10, 10}, {23, 13, 9}, {22, 11, 10},
{23, 11, 10}, {24, 14, 9}, {23, 12, 10}, {24, 12, 10}, {23, 10, 11}, {24, 13, 10}, {25, 13, 10}, {24, 11, 11},
{25, 14, 10}, {24, 12, 11}, {25, 12, 11}, {26, 15, 10}, {25, 13, 11}, {26, 13, 11}, {25, 11, 12}, {26, 14, 11},
{27, 14, 11}, {26, 12, 12}, {27, 15, 11}, {28, 15, 11}, {27, 13, 12}, {28, 13, 12}, {27, 14, 12}, {28, 14, 12},
{29, 14, 12}, {28, 15, 12}, {29, 15, 12}, {28, 13, 13}, {29, 16, 12}, {30, 16, 12}, {29, 14, 13}, {30, 14, 13},
{31, 17, 12}, {30, 15, 13}, {31, 15, 13}, {30, 16, 13}, {31, 16, 13}, {30, 14, 14}, {31, 14, 14}, {32, 17, 13},
{31, 15, 14}, {32, 15, 14}, {33, 18, 13}, {32, 16, 14}, {33, 16, 14}, {34, 19, 13}, {33, 17, 14}, {34, 17, 14},
{33, 15, 15}, {34, 18, 14}, {35, 18, 14}, {34, 16, 15}, {35, 19, 14}, {34, 17, 15}, {35, 17, 15}, {36, 20, 14},
{35, 18, 15}, {36, 18, 15}, {35, 16, 16}, {36, 19, 15}, {37, 19, 15}, {36, 17, 16}, {37, 20, 15}, {38, 20, 15},
{37, 18, 16}, {38, 18, 16}, {37, 19, 16}, {38, 19, 16}, {37, 17, 17}, {38, 20, 16}, {39, 20, 16}, {38, 18, 17},
{39, 18, 17}, {40, 21, 16}, {39, 19, 17}, {40, 19, 17}, {41, 22, 16}, {40, 20, 17}, {41, 20, 17}, {40, 21, 17},
{41, 21, 17}, {42, 21, 17}, {41, 19, 18}, {42, 22, 17}, {41, 20, 18}, {42, 20, 18}, {43, 23, 17}, {42, 21, 18},
{43, 21, 18}, {42, 19, 19}, {43, 22, 18}, {44, 22, 18}, {43, 20, 19}, {44, 23, 18}, {45, 23, 18}};
TEST(Disc, decomposition_8_check_extent)
{
for (int r = 1; r < 200; ++r)
{
auto c = mln::se::disc::_compute_decomposition_coeff(r, 8);
int extent = c[0] + 2 * c[1] + 6 * c[2];
ASSERT_EQ(r, extent);
}
}
TEST(Disc, decomposition_8_check_close_to_best_coef)
{
for (int r = 1; r < 200; ++r)
{
auto k = mln::se::disc::_compute_decomposition_coeff(r, 8);
ASSERT_NEAR(best_decomposition_coeff[r-1][0], k[0], 2);
ASSERT_NEAR(best_decomposition_coeff[r-1][1], k[1], 3);
ASSERT_NEAR(best_decomposition_coeff[r-1][2], k[2], 1);
}
}
TEST(Disc, decomposition_8_errors_does_not_degenerate)
{
int max_radius = MLN_HAS_DEBUG ? 40 : 100; // Allow processing further in release
for (int radius = 5; radius <= max_radius; ++radius)
{
int extent;
auto f = draw_ball_by_decomposition(radius, max_radius, extent);
ASSERT_EQ(2 * radius + 1, extent) << " for radius " << radius;
ASSERT_LE(compute_disc_error(f, radius), 6.5) << " for radius " << radius;
}
}
TEST(Disc, euclidean_disc_is_not_decomposable)
{
mln::se::disc d(5, 0);
ASSERT_FALSE(d.decomposable());
}
......@@ -6,32 +6,68 @@
#include <mln/io/imread.hpp>
#include <mln/morpho/structural/dilate.hpp>
#include <mln/core/se/disc.hpp>
#include <tests/helpers.hpp>
#include <gtest/gtest.h>
using namespace mln;
TEST(Morpho, dilation_dilate_0)
TEST(Dilation, Disc_euclidean)
{
mln::box2d domain{ {0,0}, {21, 21}};
mln::image2d<mln::uint8> input(domain, 0);
const mln::image2d<mln::uint8> ref = {
{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, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0},
{0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0},
{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0},
{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0},
{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0},
{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0},
{0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0},
{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0},
{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0},
{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0},
{0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0},
{0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0},
{0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0},
{0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},
};
input.at(10, 10) = 1;
auto output = mln::morpho::structural::dilate(input, mln::se::disc(9, 0));
ASSERT_IMAGES_EQ(ref, output);
}
TEST(Morpho, Generic_wide_enough_extension)
{
image2d<uint8> ima;
io::imread(MLN_IMG_PATH "small.pgm", ima);
{ // Fast: border wide enough
auto win = make_rectangle2d(7, 7);
auto out = morpho::structural::dilate(ima, win);
mln::se::disc se(3, 0);
auto out = morpho::structural::dilate(ima, se);
ASSERT_TRUE(all(out >= ima)); // extensive
}
}
// Border is not wide enough => use morpher for bound checking
TEST(Morpho, dilation_dilate_1)
TEST(Dilation, Generic_with_too_small_extension)
{
image2d<uint8> ima;
io::imread(MLN_IMG_PATH "small.pgm", ima);