Fast, Exact Euclidean Distance Transform, Linear Complexity
Posted: 2012-07-31T17:02:58-07:00
Use https://github.com/ImageMagick/ImageMagick/discussions instead.
https://imagemagick.com/discourse-server/
https://imagemagick.com/discourse-server/viewtopic.php?t=21557
Code: Select all
struct EuclideanMetric
{
static inline int f (int x_i, int gi) noexcept
{
return (x_i*x_i)+gi*gi;
}
static inline int sep (int i, int u, int gi, int gu, int) noexcept
{
return (u*u - i*i + gu*gu - gi*gi) / (2*(u-i));
}
};
struct ManhattanMetric
{
static inline int f (int x_i, int gi) noexcept
{
return abs (x_i) + gi;
}
static inline int sep (int i, int u, int gi, int gu, int inf) noexcept
{
if (gu >= gi + u - i)
return inf;
else if (gi > gu + u - i)
return -inf;
else
return (gu - gi + u + i) / 2;
}
};
struct ChessMetric
{
static inline int f (int x_i, int gi) noexcept
{
return jmax (abs (x_i), gi);
}
static inline int sep (int i, int u, int gi, int gu, int) noexcept
{
if (gi < gu)
return jmax (i+gu, (i+u)/2);
else
return jmin (u-gi, (i+u)/2);
}
};
template <class Functor, class BoolImage, class Metric>
static void calculate (Functor f, BoolImage test, int const m, int n, Metric metric)
{
std::vector <int> g (m * n);
int const inf = m + n;
// phase 1
{
for (int x = 0; x < m; ++x)
{
g [x] = test (x, 0) ? 0 : inf;
// scan 1
for (int y = 1; y < n; ++y)
{
int const ym = y*m;
g [x+ym] = test (x, y) ? 0 : 1 + g [x+ym-m];
}
// scan 2
for (int y = n-2; y >=0; --y)
{
int const ym = y*m;
if (g [x+ym+m] < g [x+ym])
g [x+ym] = 1 + g[x+ym+m];
}
}
}
// phase 2
{
std::vector <int> s (jmax (m, n));
std::vector <int> t (jmax (m, n));
for (int y = 0; y < n; ++y)
{
int q = 0;
s [0] = 0;
t [0] = 0;
int const ym = y*m;
// scan 3
for (int u = 1; u < m; ++u)
{
while (q >= 0 && metric.f (t[q]-s[q], g[s[q]+ym]) > metric.f (t[q]-u, g[u+ym]))
q--;
if (q < 0)
{
q = 0;
s [0] = u;
}
else
{
int const w = 1 + metric.sep (s[q], u, g[s[q]+ym], g[u+ym], inf);
if (w < m)
{
++q;
s[q] = u;
t[q] = w;
}
}
}
// scan 4
for (int u = m-1; u >= 0; --u)
{
int const d = metric.f (u-s[q], g[s[q]+ym]);
f (u, y, d);
if (u == t[q])
--q;
}
}
}
}