38bb7d2950
the code more readable and to highlight where PrintObject::region_volumes are actually set and consumed. 2) Replaced Slic3r::clamp() with std::clamp(). They differ in the order of their parameters, thus hopefully no new bugs were introduced. 3) Some refactoring of MultiMaterialSegmentation for efficiency.
1636 lines
81 KiB
C++
1636 lines
81 KiB
C++
// Polygon offsetting using Voronoi diagram prodiced by boost::polygon.
|
|
|
|
#include "VoronoiOffset.hpp"
|
|
#include "libslic3r.h"
|
|
|
|
#include <cmath>
|
|
|
|
// #define VORONOI_DEBUG_OUT
|
|
|
|
#include <boost/polygon/detail/voronoi_ctypes.hpp>
|
|
|
|
#ifdef VORONOI_DEBUG_OUT
|
|
#include <libslic3r/VoronoiVisualUtils.hpp>
|
|
#endif
|
|
|
|
namespace Slic3r {
|
|
namespace Voronoi {
|
|
|
|
namespace detail {
|
|
// Intersect a circle with a ray, return the two parameters.
|
|
// Currently used for unbounded Voronoi edges only.
|
|
double first_circle_segment_intersection_parameter(
|
|
const Vec2d ¢er, const double r, const Vec2d &pt, const Vec2d &v)
|
|
{
|
|
const Vec2d d = pt - center;
|
|
#ifndef NDEBUG
|
|
// Start point should be inside, end point should be outside the circle.
|
|
double d0 = (pt - center).norm();
|
|
double d1 = (pt + v - center).norm();
|
|
assert(d0 < r + SCALED_EPSILON);
|
|
assert(d1 > r - SCALED_EPSILON);
|
|
#endif /* NDEBUG */
|
|
const double a = v.squaredNorm();
|
|
const double b = 2. * d.dot(v);
|
|
const double c = d.squaredNorm() - r * r;
|
|
std::pair<int, std::array<double, 2>> out;
|
|
double u = b * b - 4. * a * c;
|
|
assert(u > - EPSILON);
|
|
double t;
|
|
if (u <= 0) {
|
|
// Degenerate to a single closest point.
|
|
t = - b / (2. * a);
|
|
assert(t >= - EPSILON && t <= 1. + EPSILON);
|
|
return std::clamp(t, 0., 1.);
|
|
} else {
|
|
u = sqrt(u);
|
|
out.first = 2;
|
|
double t0 = (- b - u) / (2. * a);
|
|
double t1 = (- b + u) / (2. * a);
|
|
// One of the intersections shall be found inside the segment.
|
|
assert((t0 >= - EPSILON && t0 <= 1. + EPSILON) || (t1 >= - EPSILON && t1 <= 1. + EPSILON));
|
|
if (t1 < 0.)
|
|
return 0.;
|
|
if (t0 > 1.)
|
|
return 1.;
|
|
return (t0 > 0.) ? t0 : t1;
|
|
}
|
|
}
|
|
|
|
struct Intersections
|
|
{
|
|
int count;
|
|
Vec2d pts[2];
|
|
};
|
|
|
|
// Return maximum two points, that are at distance "d" from both points
|
|
Intersections point_point_equal_distance_points(const Point &pt1, const Point &pt2, const double d)
|
|
{
|
|
// Calculate the two intersection points.
|
|
// With the help of Python package sympy:
|
|
// res = solve([(x - cx)**2 + (y - cy)**2 - d**2, x**2 + y**2 - d**2], [x, y])
|
|
// ccode(cse((res[0][0], res[0][1], res[1][0], res[1][1])))
|
|
// where cx, cy is the center of pt1 relative to pt2,
|
|
// d is distance from the line and the point (0, 0).
|
|
// The result is then shifted to pt2.
|
|
auto cx = double(pt1.x() - pt2.x());
|
|
auto cy = double(pt1.y() - pt2.y());
|
|
double cl = cx * cx + cy * cy;
|
|
double discr = 4. * d * d - cl;
|
|
if (discr < 0.) {
|
|
// No intersection point found, the two circles are too far away.
|
|
return Intersections { 0, { Vec2d(), Vec2d() } };
|
|
}
|
|
// Avoid division by zero if a gets too small.
|
|
bool xy_swapped = std::abs(cx) < std::abs(cy);
|
|
if (xy_swapped)
|
|
std::swap(cx, cy);
|
|
double u;
|
|
int cnt;
|
|
if (discr == 0.) {
|
|
cnt = 1;
|
|
u = 0;
|
|
} else {
|
|
cnt = 2;
|
|
u = 0.5 * cx * sqrt(cl * discr) / cl;
|
|
}
|
|
double v = 0.5 * cy - u;
|
|
double w = 2. * cy;
|
|
double e = 0.5 / cx;
|
|
double f = 0.5 * cy + u;
|
|
Intersections out { cnt, { Vec2d(-e * (v * w - cl), v),
|
|
Vec2d(-e * (w * f - cl), f) } };
|
|
if (xy_swapped) {
|
|
std::swap(out.pts[0].x(), out.pts[0].y());
|
|
std::swap(out.pts[1].x(), out.pts[1].y());
|
|
}
|
|
out.pts[0] += pt2.cast<double>();
|
|
out.pts[1] += pt2.cast<double>();
|
|
|
|
assert(std::abs((out.pts[0] - pt1.cast<double>()).norm() - d) < SCALED_EPSILON);
|
|
assert(std::abs((out.pts[1] - pt1.cast<double>()).norm() - d) < SCALED_EPSILON);
|
|
assert(std::abs((out.pts[0] - pt2.cast<double>()).norm() - d) < SCALED_EPSILON);
|
|
assert(std::abs((out.pts[1] - pt2.cast<double>()).norm() - d) < SCALED_EPSILON);
|
|
return out;
|
|
}
|
|
|
|
// Return maximum two points, that are at distance "d" from both the line and point.
|
|
Intersections line_point_equal_distance_points(const Line &line, const Point &ipt, const double d)
|
|
{
|
|
assert(line.a != ipt && line.b != ipt);
|
|
// Calculating two points of distance "d" to a ray and a point.
|
|
// Point.
|
|
Vec2d pt = ipt.cast<double>();
|
|
Vec2d lv = (line.b - line.a).cast<double>();
|
|
double l2 = lv.squaredNorm();
|
|
Vec2d lpv = (line.a - ipt).cast<double>();
|
|
double c = cross2(lpv, lv);
|
|
if (c < 0) {
|
|
lv = - lv;
|
|
c = - c;
|
|
}
|
|
|
|
// Line equation (ax + by + c - d * sqrt(l2)).
|
|
auto a = - lv.y();
|
|
auto b = lv.x();
|
|
// Line point shifted by -ipt is on the line.
|
|
assert(std::abs(lpv.x() * a + lpv.y() * b + c) < SCALED_EPSILON);
|
|
// Line vector (a, b) points towards ipt.
|
|
assert(a * lpv.x() + b * lpv.y() < - SCALED_EPSILON);
|
|
|
|
#ifndef NDEBUG
|
|
{
|
|
// Foot point of ipt on line.
|
|
Vec2d ft = Geometry::foot_pt(line, ipt);
|
|
// Center point between ipt and line, its distance to both line and ipt is equal.
|
|
Vec2d centerpt = 0.5 * (ft + pt) - pt;
|
|
double dcenter = 0.5 * (ft - pt).norm();
|
|
// Verify that the center point
|
|
assert(std::abs(centerpt.x() * a + centerpt.y() * b + c - dcenter * sqrt(l2)) < SCALED_EPSILON * sqrt(l2));
|
|
}
|
|
#endif // NDEBUG
|
|
|
|
// Calculate the two intersection points.
|
|
// With the help of Python package sympy:
|
|
// res = solve([a * x + b * y + c - d * sqrt(a**2 + b**2), x**2 + y**2 - d**2], [x, y])
|
|
// ccode(cse((res[0][0], res[0][1], res[1][0], res[1][1])))
|
|
// where (a, b, c, d) is the line equation, not normalized (vector a,b is not normalized),
|
|
// d is distance from the line and the point (0, 0).
|
|
// The result is then shifted to ipt.
|
|
|
|
double dscaled = d * sqrt(l2);
|
|
double s = c * (2. * dscaled - c);
|
|
if (s < 0.)
|
|
// Distance of pt from line is bigger than 2 * d.
|
|
return Intersections { 0 };
|
|
double u;
|
|
int cnt;
|
|
// Avoid division by zero if a gets too small.
|
|
bool xy_swapped = std::abs(a) < std::abs(b);
|
|
if (xy_swapped)
|
|
std::swap(a, b);
|
|
if (s == 0.) {
|
|
// Distance of pt from line is 2 * d.
|
|
cnt = 1;
|
|
u = 0.;
|
|
} else {
|
|
// Distance of pt from line is smaller than 2 * d.
|
|
cnt = 2;
|
|
u = a * sqrt(s) / l2;
|
|
}
|
|
double e = dscaled - c;
|
|
double f = b * e / l2;
|
|
double g = f - u;
|
|
double h = f + u;
|
|
Intersections out { cnt, { Vec2d((- b * g + e) / a, g),
|
|
Vec2d((- b * h + e) / a, h) } };
|
|
if (xy_swapped) {
|
|
std::swap(out.pts[0].x(), out.pts[0].y());
|
|
std::swap(out.pts[1].x(), out.pts[1].y());
|
|
}
|
|
out.pts[0] += pt;
|
|
out.pts[1] += pt;
|
|
|
|
assert(std::abs(Geometry::ray_point_distance<Vec2d>(line.a.cast<double>(), (line.b - line.a).cast<double>(), out.pts[0]) - d) < SCALED_EPSILON);
|
|
assert(std::abs(Geometry::ray_point_distance<Vec2d>(line.a.cast<double>(), (line.b - line.a).cast<double>(), out.pts[1]) - d) < SCALED_EPSILON);
|
|
assert(std::abs((out.pts[0] - ipt.cast<double>()).norm() - d) < SCALED_EPSILON);
|
|
assert(std::abs((out.pts[1] - ipt.cast<double>()).norm() - d) < SCALED_EPSILON);
|
|
return out;
|
|
}
|
|
|
|
// Double vertex equal to a coord_t point after conversion to double.
|
|
template<typename VertexType>
|
|
inline bool vertex_equal_to_point(const VertexType &vertex, const Point &ipt)
|
|
{
|
|
// Convert ipt to doubles, force the 80bit FPU temporary to 64bit and then compare.
|
|
// This should work with any settings of math compiler switches and the C++ compiler
|
|
// shall understand the memcpies as type punning and it shall optimize them out.
|
|
#if 1
|
|
using ulp_cmp_type = boost::polygon::detail::ulp_comparison<double>;
|
|
ulp_cmp_type ulp_cmp;
|
|
static constexpr int ULPS = boost::polygon::voronoi_diagram_traits<double>::vertex_equality_predicate_type::ULPS;
|
|
return ulp_cmp(vertex.x(), double(ipt.x()), ULPS) == ulp_cmp_type::EQUAL &&
|
|
ulp_cmp(vertex.y(), double(ipt.y()), ULPS) == ulp_cmp_type::EQUAL;
|
|
#else
|
|
volatile double u = static_cast<double>(ipt.x());
|
|
volatile double v = vertex.x();
|
|
if (u != v)
|
|
return false;
|
|
u = static_cast<double>(ipt.y());
|
|
v = vertex.y();
|
|
return u == v;
|
|
#endif
|
|
};
|
|
bool vertex_equal_to_point(const VD::vertex_type *vertex, const Point &ipt)
|
|
{ return vertex_equal_to_point(*vertex, ipt); }
|
|
|
|
double dist_to_site(const Lines &lines, const VD::cell_type &cell, const Vec2d &point)
|
|
{
|
|
const Line &line = lines[cell.source_index()];
|
|
return cell.contains_point() ?
|
|
(((cell.source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line.a : line.b).cast<double>() - point).norm() :
|
|
(Geometry::foot_pt<Vec2d>(line.a.cast<double>(), (line.b - line.a).cast<double>(), point) - point).norm();
|
|
};
|
|
|
|
bool on_site(const Lines &lines, const VD::cell_type &cell, const Vec2d &pt)
|
|
{
|
|
const Line &line = lines[cell.source_index()];
|
|
auto on_contour = [&pt](const Point &ipt) { return detail::vertex_equal_to_point(pt, ipt); };
|
|
if (cell.contains_point()) {
|
|
return on_contour(contour_point(cell, line));
|
|
} else {
|
|
assert(! (on_contour(line.a) && on_contour(line.b)));
|
|
return on_contour(line.a) || on_contour(line.b);
|
|
}
|
|
};
|
|
|
|
// For a Voronoi segment starting with voronoi_point1 and ending with voronoi_point2,
|
|
// defined by a bisector of Voronoi sites pt1_site and pt2_site (line)
|
|
// find two points on the Voronoi bisector, that delimit regions with dr/dl measure
|
|
// lower / higher than threshold_dr_dl.
|
|
//
|
|
// Linear segment from voronoi_point1 to return.first and
|
|
// linear segment from return.second to voronoi_point2 have dr/dl measure
|
|
// higher than threshold_dr_dl.
|
|
// If such respective segment does not exist, then return.first resp. return.second is nan.
|
|
std::pair<Vec2d, Vec2d> point_point_dr_dl_thresholds(
|
|
// Two Voronoi sites
|
|
const Point &pt1_site, const Point &pt2_site,
|
|
// End points of a Voronoi segment
|
|
const Vec2d &voronoi_point1, const Vec2d &voronoi_point2,
|
|
// Threshold of the skeleton function, where alpha is an angle of a sharp convex corner with the same dr/dl.
|
|
const double threshold_tan_alpha_half)
|
|
{
|
|
// sympy code to calculate +-x
|
|
// of a linear bisector of pt1_site, pt2_site parametrized with pt + x * v, |v| = 1
|
|
// where dr/dl = threshold_dr_dl
|
|
// equals d|pt1_site - pt + x * v| / dx = threshold_dr_dl
|
|
//
|
|
// y = sqrt(x^2 + d^2)
|
|
// dy = diff(y, x)
|
|
// solve(dy - c, x)
|
|
//
|
|
// Project voronoi_point1/2 to line_site.
|
|
Vec2d dir_y = (pt2_site - pt1_site).cast<double>();
|
|
Vec2d dir_x = Vec2d(- dir_y.y(), dir_y.x()).normalized();
|
|
Vec2d cntr = 0.5 * (pt1_site.cast<double>() + pt2_site.cast<double>());
|
|
double t1 = (voronoi_point1 - cntr).dot(dir_x);
|
|
double t2 = (voronoi_point2 - cntr).dot(dir_x);
|
|
if (t1 > t2) {
|
|
t1 = -t1;
|
|
t2 = -t2;
|
|
dir_x = - dir_x;
|
|
}
|
|
auto x = 0.5 * dir_y.norm() * threshold_tan_alpha_half;
|
|
static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
|
|
auto out = std::make_pair(Vec2d(nan, nan), Vec2d(nan, nan));
|
|
if (t2 > -x && t1 < x) {
|
|
// Intervals overlap.
|
|
dir_x *= x;
|
|
out.first = (t1 < -x) ? cntr - dir_x : voronoi_point1;
|
|
out.second = (t2 > +x) ? cntr + dir_x : voronoi_point2;
|
|
}
|
|
return out;
|
|
}
|
|
|
|
// For a Voronoi segment starting with voronoi_point1 and ending with voronoi_point2,
|
|
// defined by a bisector of Voronoi sites pt_site and line site (parabolic arc)
|
|
// find two points on the Voronoi parabolic arc, that delimit regions with dr/dl measure
|
|
// lower / higher than threshold_dr_dl.
|
|
//
|
|
// Parabolic arc from voronoi_point1 to return.first and
|
|
// parabolic arc from return.second to voronoi_point2 have dr/dl measure
|
|
// higher than threshold_dr_dl.
|
|
// If such respective segment does not exist, then return.first resp. return.second is nan.
|
|
std::pair<Vec2d, Vec2d> point_segment_dr_dl_thresholds(
|
|
// Two Voronoi sites
|
|
const Point &pt_site, const Line &line_site,
|
|
// End points of a Voronoi segment
|
|
const Vec2d &voronoi_point1, const Vec2d &voronoi_point2,
|
|
// Threshold of the skeleton function, where alpha is an angle of a sharp convex corner with the same dr/dl.
|
|
const double threshold_tan_alpha_half)
|
|
{
|
|
// sympy code to calculate +-x
|
|
// of a parabola y = ax^2 + b
|
|
// where dr/dl = threshold_dr_dl
|
|
//
|
|
// a = 1 / (4 * b)
|
|
// y = a*x**2 + b
|
|
// dy = diff(y, x)
|
|
// solve(dy / sqrt(1 + dy**2) - c, x)
|
|
//
|
|
// Foot point of the point site on the line site.
|
|
Vec2d ft = Geometry::foot_pt(line_site, pt_site);
|
|
// Minimum distance of the bisector (parabolic arc) from the two sites, squared.
|
|
Vec2d dir_pt_ft = pt_site.cast<double>() - ft;
|
|
double b = 0.5 * dir_pt_ft.norm();
|
|
static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
|
|
auto out = std::make_pair(Vec2d(nan, nan), Vec2d(nan, nan));
|
|
{
|
|
// +x, -x are the two parameters along the line_site, where threshold_tan_alpha_half is met.
|
|
double x = 2. * b * threshold_tan_alpha_half;
|
|
// Project voronoi_point1/2 to line_site.
|
|
Vec2d dir_x = (line_site.b - line_site.a).cast<double>().normalized();
|
|
double t1 = (voronoi_point1 - ft).dot(dir_x);
|
|
double t2 = (voronoi_point2 - ft).dot(dir_x);
|
|
if (t1 > t2) {
|
|
t1 = -t1;
|
|
t2 = -t2;
|
|
dir_x = - dir_x;
|
|
}
|
|
if (t2 > -x && t1 < x) {
|
|
// Intervals overlap.
|
|
bool t1_valid = t1 < -x;
|
|
bool t2_valid = t2 > +x;
|
|
// Direction of the Y axis of the parabola.
|
|
Vec2d dir_y(- dir_x.y(), dir_x.x());
|
|
// Orient the Y axis towards the point site.
|
|
if (dir_y.dot(dir_pt_ft) < 0.)
|
|
dir_y = - dir_y;
|
|
// Equation of the parabola: y = b + a * x^2
|
|
double a = 0.25 / b;
|
|
dir_x *= x;
|
|
dir_y *= b + a * x * x;
|
|
out.first = t1_valid ? ft - dir_x + dir_y : voronoi_point1;
|
|
out.second = t2_valid ? ft + dir_x + dir_y : voronoi_point2;
|
|
}
|
|
}
|
|
return out;
|
|
}
|
|
|
|
std::pair<Vec2d, Vec2d> point_point_skeleton_thresholds(
|
|
// Two Voronoi sites
|
|
const Point &pt1_site, const Point &pt2_site,
|
|
// End points of a Voronoi segment
|
|
const Vec2d &voronoi_point1, const Vec2d &voronoi_point2,
|
|
// Threshold of the skeleton function.
|
|
const double tan_alpha_half)
|
|
{
|
|
// Project voronoi_point1/2 to line_site.
|
|
Vec2d dir_y = (pt2_site - pt1_site).cast<double>();
|
|
Vec2d dir_x = Vec2d(- dir_y.y(), dir_y.x()).normalized();
|
|
Vec2d cntr = 0.5 * (pt1_site.cast<double>() + pt2_site.cast<double>());
|
|
double t1 = (voronoi_point1 - cntr).dot(dir_x);
|
|
double t2 = (voronoi_point2 - cntr).dot(dir_x);
|
|
if (t1 > t2) {
|
|
t1 = -t1;
|
|
t2 = -t2;
|
|
dir_x = - dir_x;
|
|
}
|
|
auto x = 0.5 * dir_y.norm() * tan_alpha_half;
|
|
static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
|
|
auto out = std::make_pair(Vec2d(nan, nan), Vec2d(nan, nan));
|
|
if (t2 > -x && t1 < x) {
|
|
// Intervals overlap.
|
|
dir_x *= x;
|
|
out.first = (t1 < -x) ? cntr - dir_x : voronoi_point1;
|
|
out.second = (t2 > +x) ? cntr + dir_x : voronoi_point2;
|
|
}
|
|
return out;
|
|
}
|
|
|
|
std::pair<Vec2d, Vec2d> point_segment_skeleton_thresholds(
|
|
// Two Voronoi sites
|
|
const Point &pt_site, const Line &line_site,
|
|
// End points of a Voronoi segment
|
|
const Vec2d &voronoi_point1, const Vec2d &voronoi_point2,
|
|
// Threshold of the skeleton function.
|
|
const double threshold_cos_alpha)
|
|
{
|
|
// Foot point of the point site on the line site.
|
|
Vec2d ft = Geometry::foot_pt(line_site, pt_site);
|
|
// Minimum distance of the bisector (parabolic arc) from the two sites, squared.
|
|
Vec2d dir_pt_ft = pt_site.cast<double>() - ft;
|
|
// Distance of Voronoi point site from the Voronoi line site.
|
|
double l = dir_pt_ft.norm();
|
|
static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
|
|
auto out = std::make_pair(Vec2d(nan, nan), Vec2d(nan, nan));
|
|
// +x, -x are the two parameters along the line_site, where threshold is met.
|
|
double r = l / (1. + threshold_cos_alpha);
|
|
double x2 = r * r - Slic3r::sqr(l - r);
|
|
double x = sqrt(x2);
|
|
// Project voronoi_point1/2 to line_site.
|
|
Vec2d dir_x = (line_site.b - line_site.a).cast<double>().normalized();
|
|
double t1 = (voronoi_point1 - ft).dot(dir_x);
|
|
double t2 = (voronoi_point2 - ft).dot(dir_x);
|
|
if (t1 > t2) {
|
|
t1 = -t1;
|
|
t2 = -t2;
|
|
dir_x = - dir_x;
|
|
}
|
|
if (t2 > -x && t1 < x) {
|
|
// Intervals overlap.
|
|
bool t1_valid = t1 < -x;
|
|
bool t2_valid = t2 > +x;
|
|
// Direction of the Y axis of the parabola.
|
|
Vec2d dir_y(- dir_x.y(), dir_x.x());
|
|
// Orient the Y axis towards the point site.
|
|
if (dir_y.dot(dir_pt_ft) < 0.)
|
|
dir_y = - dir_y;
|
|
// Equation of the parabola: y = b + a * x^2
|
|
double a = 0.5 / l;
|
|
dir_x *= x;
|
|
dir_y *= 0.5 * l + a * x2;
|
|
out.first = t1_valid ? ft - dir_x + dir_y : voronoi_point1;
|
|
out.second = t2_valid ? ft + dir_x + dir_y : voronoi_point2;
|
|
}
|
|
return out;
|
|
}
|
|
|
|
} // namespace detail
|
|
|
|
#ifndef NDEBUG
|
|
namespace debug
|
|
{
|
|
// Verify that twin halfedges are stored next to the other in vd.
|
|
bool verify_twin_halfedges_successive(const VD &vd, const Lines &lines)
|
|
{
|
|
for (size_t i = 0; i < vd.num_edges(); i += 2) {
|
|
const VD::edge_type &e = vd.edges()[i];
|
|
const VD::edge_type &e2 = vd.edges()[i + 1];
|
|
assert(e.twin() == &e2);
|
|
assert(e2.twin() == &e);
|
|
assert(e.is_secondary() == e2.is_secondary());
|
|
if (e.is_secondary()) {
|
|
assert(e.cell()->contains_point() != e2.cell()->contains_point());
|
|
const VD::edge_type &ex = (e.cell()->contains_point() ? e : e2);
|
|
// Verify that the Point defining the cell left of ex is an end point of a segment
|
|
// defining the cell right of ex.
|
|
const Line &line0 = lines[ex.cell()->source_index()];
|
|
const Line &line1 = lines[ex.twin()->cell()->source_index()];
|
|
const Point &pt = (ex.cell()->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b;
|
|
assert(pt == line1.a || pt == line1.b);
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
bool verify_inside_outside_annotations(const VD &vd)
|
|
{
|
|
// Verify that "Colors" are set at all Voronoi entities.
|
|
for (const VD::vertex_type &v : vd.vertices()) {
|
|
assert(! v.is_degenerate());
|
|
assert(vertex_category(v) != VertexCategory::Unknown);
|
|
}
|
|
for (const VD::edge_type &e : vd.edges())
|
|
assert(edge_category(e) != EdgeCategory::Unknown);
|
|
for (const VD::cell_type &c : vd.cells()) {
|
|
// Unfortunately denegerate cells could be created, which reference a null edge.
|
|
// https://github.com/boostorg/polygon/issues/47
|
|
assert(c.is_degenerate() || cell_category(c) != CellCategory::Unknown);
|
|
}
|
|
|
|
// Verify consistency between markings of Voronoi cells, edges and verticies.
|
|
for (const VD::cell_type &cell : vd.cells()) {
|
|
if (cell.is_degenerate()) {
|
|
// Unfortunately denegerate cells could be created, which reference a null edge.
|
|
// https://github.com/boostorg/polygon/issues/47
|
|
continue;
|
|
}
|
|
const VD::edge_type *first_edge = cell.incident_edge();
|
|
const VD::edge_type *edge = first_edge;
|
|
CellCategory cc = cell_category(cell);
|
|
size_t num_vertices_on_contour = 0;
|
|
size_t num_vertices_inside = 0;
|
|
size_t num_vertices_outside = 0;
|
|
size_t num_edges_point_to_contour = 0;
|
|
size_t num_edges_point_inside = 0;
|
|
size_t num_edges_point_outside = 0;
|
|
do {
|
|
{
|
|
EdgeCategory ec = edge_category(edge);
|
|
switch (ec) {
|
|
case EdgeCategory::PointsInside:
|
|
assert(edge->vertex0() != nullptr && edge->vertex1() != nullptr);
|
|
++ num_edges_point_inside; break;
|
|
case EdgeCategory::PointsOutside:
|
|
// assert(edge->vertex0() != nullptr);
|
|
++ num_edges_point_outside; break;
|
|
case EdgeCategory::PointsToContour:
|
|
assert(edge->vertex1() != nullptr);
|
|
++ num_edges_point_to_contour; break;
|
|
default:
|
|
assert(false);
|
|
}
|
|
}
|
|
{
|
|
VertexCategory vc = (edge->vertex1() == nullptr) ? VertexCategory::Outside : vertex_category(edge->vertex1());
|
|
switch (vc) {
|
|
case VertexCategory::Inside:
|
|
++ num_vertices_inside; break;
|
|
case VertexCategory::Outside:
|
|
++ num_vertices_outside; break;
|
|
case VertexCategory::OnContour:
|
|
++ num_vertices_on_contour; break;
|
|
default:
|
|
assert(false);
|
|
}
|
|
}
|
|
{
|
|
const VD::cell_type *cell_other = edge->twin()->cell();
|
|
const CellCategory cc_other = cell_category(cell_other);
|
|
assert(cc_other != CellCategory::Unknown);
|
|
switch (cc) {
|
|
case CellCategory::Boundary:
|
|
assert(cc_other != CellCategory::Boundary || cell_other->contains_segment());
|
|
break;
|
|
case CellCategory::Inside:
|
|
assert(cc_other == CellCategory::Inside || cc_other ==CellCategory::Boundary);
|
|
break;
|
|
case CellCategory::Outside:
|
|
assert(cc_other == CellCategory::Outside || cc_other == CellCategory::Boundary);
|
|
break;
|
|
default:
|
|
assert(false);
|
|
break;
|
|
}
|
|
}
|
|
edge = edge->next();
|
|
} while (edge != first_edge);
|
|
|
|
switch (cc) {
|
|
case CellCategory::Boundary:
|
|
assert(cell.contains_segment());
|
|
assert(num_edges_point_to_contour == 2);
|
|
assert(num_vertices_on_contour == 2);
|
|
assert(num_vertices_inside > 0);
|
|
assert(num_vertices_outside > 0);
|
|
assert(num_edges_point_inside > 0);
|
|
assert(num_edges_point_outside > 0);
|
|
break;
|
|
case CellCategory::Inside:
|
|
assert(num_vertices_on_contour <= 1);
|
|
assert(num_edges_point_to_contour <= 1);
|
|
assert(num_vertices_inside > 0);
|
|
assert(num_vertices_outside == 0);
|
|
assert(num_edges_point_inside > 0);
|
|
assert(num_edges_point_outside == 0);
|
|
break;
|
|
case CellCategory::Outside:
|
|
assert(num_vertices_on_contour <= 1);
|
|
assert(num_edges_point_to_contour <= 1);
|
|
assert(num_vertices_inside == 0);
|
|
assert(num_vertices_outside > 0);
|
|
assert(num_edges_point_inside == 0);
|
|
assert(num_edges_point_outside > 0);
|
|
break;
|
|
default:
|
|
assert(false);
|
|
break;
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
bool verify_vertices_on_contour(const VD &vd, const Lines &lines)
|
|
{
|
|
for (const VD::edge_type &edge : vd.edges()) {
|
|
const VD::vertex_type *v = edge.vertex0();
|
|
if (v != nullptr) {
|
|
bool on_contour = vertex_category(v) == VertexCategory::OnContour;
|
|
assert(detail::on_site(lines, *edge.cell(), vertex_point(v)) == on_contour);
|
|
assert(detail::on_site(lines, *edge.twin()->cell(), vertex_point(v)) == on_contour);
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
bool verify_signed_distances(const VD &vd, const Lines &lines, const std::vector<double> &signed_distances)
|
|
{
|
|
for (const VD::edge_type &edge : vd.edges()) {
|
|
const VD::vertex_type *v = edge.vertex0();
|
|
double d = (v == nullptr) ? std::numeric_limits<double>::max() : signed_distances[v - &vd.vertices().front()];
|
|
if (v == nullptr || vertex_category(v) == VertexCategory::Outside)
|
|
assert(d > 0.);
|
|
else if (vertex_category(v) == VertexCategory::OnContour)
|
|
assert(d == 0.);
|
|
else
|
|
assert(d < 0.);
|
|
if (v != nullptr) {
|
|
double err = std::abs(detail::dist_to_site(lines, *edge.cell(), vertex_point(v)) - std::abs(d));
|
|
double err2 = std::abs(detail::dist_to_site(lines, *edge.twin()->cell(), vertex_point(v)) - std::abs(d));
|
|
assert(err < SCALED_EPSILON);
|
|
assert(err2 < SCALED_EPSILON);
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
bool verify_offset_intersection_points(const VD &vd, const Lines &lines, const double offset_distance, const std::vector<Vec2d> &offset_intersection_points)
|
|
{
|
|
const VD::edge_type *front_edge = &vd.edges().front();
|
|
const double d = std::abs(offset_distance);
|
|
for (const VD::edge_type &edge : vd.edges()) {
|
|
const Vec2d &p = offset_intersection_points[&edge - front_edge];
|
|
if (edge_offset_has_intersection(p)) {
|
|
double err = std::abs(detail::dist_to_site(lines, *edge.cell(), p) - d);
|
|
double err2 = std::abs(detail::dist_to_site(lines, *edge.twin()->cell(), p) - d);
|
|
assert(err < SCALED_EPSILON);
|
|
assert(err2 < SCALED_EPSILON);
|
|
}
|
|
}
|
|
return true;
|
|
}
|
|
|
|
}
|
|
#endif // NDEBUG
|
|
|
|
void reset_inside_outside_annotations(VD &vd)
|
|
{
|
|
for (const VD::vertex_type &v : vd.vertices())
|
|
set_vertex_category(const_cast<VD::vertex_type&>(v), VertexCategory::Unknown);
|
|
for (const VD::edge_type &e : vd.edges())
|
|
set_edge_category(const_cast<VD::edge_type&>(e), EdgeCategory::Unknown);
|
|
for (const VD::cell_type &c : vd.cells())
|
|
set_cell_category(const_cast<VD::cell_type&>(c), CellCategory::Unknown);
|
|
}
|
|
|
|
void annotate_inside_outside(VD &vd, const Lines &lines)
|
|
{
|
|
assert(debug::verify_twin_halfedges_successive(vd, lines));
|
|
|
|
reset_inside_outside_annotations(vd);
|
|
|
|
#ifdef VORONOI_DEBUG_OUT
|
|
BoundingBox bbox;
|
|
{
|
|
bbox.merge(get_extents(lines));
|
|
bbox.min -= (0.01 * bbox.size().cast<double>()).cast<coord_t>();
|
|
bbox.max += (0.01 * bbox.size().cast<double>()).cast<coord_t>();
|
|
}
|
|
static int irun = 0;
|
|
++ irun;
|
|
dump_voronoi_to_svg(debug_out_path("voronoi-offset-initial-%d.svg", irun).c_str(), vd, Points(), lines);
|
|
#endif // VORONOI_DEBUG_OUT
|
|
|
|
// Set a VertexCategory, verify validity of the operation.
|
|
auto annotate_vertex = [](const VD::vertex_type *vertex, VertexCategory new_vertex_category) {
|
|
#ifndef NDEBUG
|
|
VertexCategory vc = vertex_category(vertex);
|
|
assert(vc == VertexCategory::Unknown || vc == new_vertex_category);
|
|
assert(new_vertex_category == VertexCategory::Inside ||
|
|
new_vertex_category == VertexCategory::Outside ||
|
|
new_vertex_category == VertexCategory::OnContour);
|
|
#endif // NDEBUG
|
|
set_vertex_category(const_cast<VD::vertex_type*>(vertex), new_vertex_category);
|
|
};
|
|
|
|
// Set an EdgeCategory, verify validity of the operation.
|
|
auto annotate_edge = [](const VD::edge_type *edge, EdgeCategory new_edge_category) {
|
|
#ifndef NDEBUG
|
|
EdgeCategory ec = edge_category(edge);
|
|
assert(ec == EdgeCategory::Unknown || ec == new_edge_category);
|
|
switch (new_edge_category) {
|
|
case EdgeCategory::PointsInside:
|
|
assert(edge->vertex0() != nullptr);
|
|
assert(edge->vertex1() != nullptr);
|
|
break;
|
|
case EdgeCategory::PointsOutside:
|
|
// assert(edge->vertex0() != nullptr);
|
|
break;
|
|
case EdgeCategory::PointsToContour:
|
|
assert(edge->vertex1() != nullptr);
|
|
break;
|
|
default:
|
|
assert(false);
|
|
}
|
|
#endif // NDEBUG
|
|
set_edge_category(const_cast<VD::edge_type*>(edge), new_edge_category);
|
|
};
|
|
|
|
// Set a CellCategory, verify validity of the operation.
|
|
// Handle marking of boundary cells (first time the cell is marked as outside, the other time as inside).
|
|
// Returns true if the current cell category was modified.
|
|
auto annotate_cell = [](const VD::cell_type *cell, CellCategory new_cell_category) -> bool {
|
|
CellCategory cc = cell_category(cell);
|
|
assert(cc == CellCategory::Inside || cc == CellCategory::Outside || cc == CellCategory::Boundary || cc == CellCategory::Unknown);
|
|
assert(new_cell_category == CellCategory::Inside || new_cell_category == CellCategory::Outside || new_cell_category == CellCategory::Boundary);
|
|
switch (cc) {
|
|
case CellCategory::Unknown:
|
|
// Old category unknown, just write the new category.
|
|
break;
|
|
case CellCategory::Outside:
|
|
if (new_cell_category == CellCategory::Inside)
|
|
new_cell_category = CellCategory::Boundary;
|
|
break;
|
|
case CellCategory::Inside:
|
|
if (new_cell_category == CellCategory::Outside)
|
|
new_cell_category = CellCategory::Boundary;
|
|
break;
|
|
case CellCategory::Boundary:
|
|
return false;
|
|
}
|
|
if (cc != new_cell_category) {
|
|
set_cell_category(const_cast<VD::cell_type*>(cell), new_cell_category);
|
|
return true;
|
|
}
|
|
return false;
|
|
};
|
|
|
|
// The next loop is supposed to annotate the "on input contour" vertices, but due to
|
|
// merging very close Voronoi vertices together by boost::polygon Voronoi generator
|
|
// the condition may not always be met. It should be safe to mark all Voronoi very close
|
|
// to the input contour as on contour.
|
|
for (const VD::edge_type &edge : vd.edges()) {
|
|
const VD::vertex_type *v = edge.vertex0();
|
|
if (v != nullptr) {
|
|
bool on_contour = detail::on_site(lines, *edge.cell(), vertex_point(v));
|
|
#ifndef NDEBUG
|
|
bool on_contour2 = detail::on_site(lines, *edge.twin()->cell(), vertex_point(v));
|
|
assert(on_contour == on_contour2);
|
|
#endif // NDEBUG
|
|
if (on_contour)
|
|
annotate_vertex(v, VertexCategory::OnContour);
|
|
}
|
|
}
|
|
|
|
// One side of a secondary edge is always on the source contour. Mark these vertices as OnContour.
|
|
// See the comment at the loop before, the condition may not always be met.
|
|
for (const VD::edge_type &edge : vd.edges()) {
|
|
if (edge.is_secondary() && edge.vertex0() != nullptr) {
|
|
assert(edge.is_linear());
|
|
assert(edge.cell()->contains_point() != edge.twin()->cell()->contains_point());
|
|
// The point associated with the point site shall be equal with one vertex of this Voronoi edge.
|
|
const Point &pt_on_contour = edge.cell()->contains_point() ? contour_point(*edge.cell(), lines) : contour_point(*edge.twin()->cell(), lines);
|
|
auto on_contour = [&pt_on_contour](const VD::vertex_type *v) { return detail::vertex_equal_to_point(v, pt_on_contour); };
|
|
if (edge.vertex1() == nullptr) {
|
|
assert(on_contour(edge.vertex0()));
|
|
annotate_vertex(edge.vertex0(), VertexCategory::OnContour);
|
|
} else {
|
|
// Only one of the two vertices may lie on input contour.
|
|
const VD::vertex_type *v0 = edge.vertex0();
|
|
const VD::vertex_type *v1 = edge.vertex1();
|
|
#ifndef NDEBUG
|
|
VertexCategory v0_category = vertex_category(v0);
|
|
VertexCategory v1_category = vertex_category(v1);
|
|
assert(v0_category != VertexCategory::OnContour || v1_category != VertexCategory::OnContour);
|
|
assert(! (on_contour(v0) && on_contour(v1)));
|
|
#endif // NDEBUG
|
|
if (on_contour(v0))
|
|
annotate_vertex(v0, VertexCategory::OnContour);
|
|
else {
|
|
assert(on_contour(v1));
|
|
annotate_vertex(v1, VertexCategory::OnContour);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
assert(debug::verify_vertices_on_contour(vd, lines));
|
|
|
|
for (const VD::edge_type &edge : vd.edges())
|
|
if (edge.vertex1() == nullptr) {
|
|
// Infinite Voronoi edge separating two Point sites or a Point site and a Segment site.
|
|
// Infinite edge is always outside and it references at least one valid vertex.
|
|
assert(edge.is_infinite());
|
|
assert(edge.is_linear());
|
|
assert(edge.vertex0() != nullptr);
|
|
const VD::cell_type *cell = edge.cell();
|
|
const VD::cell_type *cell2 = edge.twin()->cell();
|
|
// A Point-Segment secondary Voronoi edge touches the input contour, a Point-Point Voronoi
|
|
// edge does not.
|
|
assert(edge.is_secondary() ? (cell->contains_segment() != cell2->contains_segment()) :
|
|
(cell->contains_point() == cell2->contains_point()));
|
|
annotate_edge(&edge, EdgeCategory::PointsOutside);
|
|
// Opposite edge of an infinite edge is certainly not active.
|
|
annotate_edge(edge.twin(), edge.is_secondary() ? EdgeCategory::PointsToContour : EdgeCategory::PointsOutside);
|
|
annotate_vertex(edge.vertex0(), edge.is_secondary() ? VertexCategory::OnContour : VertexCategory::Outside);
|
|
// edge.vertex1() is null, it is implicitely outside.
|
|
if (cell->contains_segment())
|
|
std::swap(cell, cell2);
|
|
// State of a cell containing a boundary point is certainly outside.
|
|
assert(cell->contains_point());
|
|
annotate_cell(cell, CellCategory::Outside);
|
|
assert(edge.is_secondary() == cell2->contains_segment());
|
|
annotate_cell(cell2, cell2->contains_point() ? CellCategory::Outside : CellCategory::Boundary);
|
|
} else if (edge.vertex0() != nullptr) {
|
|
assert(edge.is_finite());
|
|
const VD::cell_type *cell = edge.cell();
|
|
const Line *line = cell->contains_segment() ? &lines[cell->source_index()] : nullptr;
|
|
if (line == nullptr) {
|
|
cell = edge.twin()->cell();
|
|
line = cell->contains_segment() ? &lines[cell->source_index()] : nullptr;
|
|
}
|
|
// Only one of the two vertices may lie on input contour.
|
|
assert(! edge.is_linear() || vertex_category(edge.vertex0()) != VertexCategory::OnContour || vertex_category(edge.vertex1()) != VertexCategory::OnContour);
|
|
// Now classify the Voronoi vertices and edges as inside outside, if at least one Voronoi
|
|
// site is a Segment site.
|
|
// Inside / outside classification of Point - Point Voronoi edges will be done later
|
|
// by a propagation (seed fill).
|
|
if (line) {
|
|
const VD::vertex_type *v1 = edge.vertex1();
|
|
const VD::cell_type *cell2 = (cell == edge.cell()) ? edge.twin()->cell() : edge.cell();
|
|
assert(v1 != nullptr);
|
|
VertexCategory v0_category = vertex_category(edge.vertex0());
|
|
VertexCategory v1_category = vertex_category(edge.vertex1());
|
|
bool on_contour = v0_category == VertexCategory::OnContour || v1_category == VertexCategory::OnContour;
|
|
#ifndef NDEBUG
|
|
if (! on_contour && cell == edge.cell() && edge.twin()->cell()->contains_segment()) {
|
|
// Constrained bisector of two segments. Vojtech is not quite sure whether the Voronoi generator is robust enough
|
|
// to always connect at least one secondary edge to an input contour point. Catch it here.
|
|
assert(edge.is_linear());
|
|
// OnContour state of this edge is not known yet.
|
|
const Point *pt_on_contour = nullptr;
|
|
// If the two segments share a point, then one end of the current Voronoi edge shares this point as well.
|
|
// A bisector may not necessarily connect to the source contour. Find pt_on_contour if it exists.
|
|
const Line &line2 = lines[cell2->source_index()];
|
|
if (line->a == line2.b)
|
|
pt_on_contour = &line->a;
|
|
else if (line->b == line2.a)
|
|
pt_on_contour = &line->b;
|
|
if (pt_on_contour) {
|
|
const VD::vertex_type *v0 = edge.vertex0();
|
|
auto on_contour = [&pt_on_contour](const VD::vertex_type *v) {
|
|
return std::abs(v->x() - pt_on_contour->x()) < 0.5001 &&
|
|
std::abs(v->y() - pt_on_contour->y()) < 0.5001;
|
|
};
|
|
assert(! on_contour(v0) && ! on_contour(v1));
|
|
}
|
|
}
|
|
#endif // NDEBUG
|
|
if (on_contour && v1_category == VertexCategory::OnContour) {
|
|
// Skip secondary edge pointing to a contour point.
|
|
annotate_edge(&edge, EdgeCategory::PointsToContour);
|
|
} else {
|
|
// v0 is certainly not on the input polygons.
|
|
// Is v1 inside or outside the input polygons?
|
|
// The Voronoi vertex coordinate is in doubles, calculate orientation in doubles.
|
|
Vec2d l0(line->a.cast<double>());
|
|
Vec2d lv((line->b - line->a).cast<double>());
|
|
double side = cross2(Vec2d(v1->x(), v1->y()) - l0, lv);
|
|
// No Voronoi edge could connect two vertices of input polygons.
|
|
assert(side != 0.);
|
|
auto vc = side > 0. ? VertexCategory::Outside : VertexCategory::Inside;
|
|
annotate_vertex(v1, vc);
|
|
auto ec = vc == VertexCategory::Outside ? EdgeCategory::PointsOutside : EdgeCategory::PointsInside;
|
|
annotate_edge(&edge, ec);
|
|
// Annotate the twin edge and its vertex. As the Voronoi edge may never cross the input
|
|
// contour, the twin edge and its vertex will share the property of edge.
|
|
annotate_vertex(edge.vertex0(), on_contour ? VertexCategory::OnContour : vc);
|
|
annotate_edge(edge.twin(), on_contour ? EdgeCategory::PointsToContour : ec);
|
|
assert(cell->contains_segment());
|
|
annotate_cell(cell, on_contour ? CellCategory::Boundary :
|
|
(vc == VertexCategory::Outside ? CellCategory::Outside : CellCategory::Inside));
|
|
annotate_cell(cell2, (on_contour && cell2->contains_segment()) ? CellCategory::Boundary :
|
|
(vc == VertexCategory::Outside ? CellCategory::Outside : CellCategory::Inside));
|
|
}
|
|
}
|
|
}
|
|
|
|
assert(debug::verify_vertices_on_contour(vd, lines));
|
|
|
|
// Now most Voronoi vertices, edges and cells are annotated, with the exception of some
|
|
// edges separating two Point sites, their cells and vertices.
|
|
// Perform one round of expansion marking Voronoi edges and cells next to boundary cells.
|
|
std::vector<const VD::cell_type*> cell_queue;
|
|
for (const VD::edge_type &edge : vd.edges()) {
|
|
assert((edge_category(edge) == EdgeCategory::Unknown) == (edge_category(edge.twin()) == EdgeCategory::Unknown));
|
|
if (edge_category(edge) == EdgeCategory::Unknown) {
|
|
assert(edge.is_finite());
|
|
const VD::cell_type &cell = *edge.cell();
|
|
const VD::cell_type &cell2 = *edge.twin()->cell();
|
|
assert(cell.contains_point() && cell2.contains_point());
|
|
CellCategory cc = cell_category(cell);
|
|
CellCategory cc2 = cell_category(cell2);
|
|
assert(cc != CellCategory::Boundary && cc2 != CellCategory::Boundary);
|
|
CellCategory cc_new = cc;
|
|
if (cc_new == CellCategory::Unknown)
|
|
cc_new = cc2;
|
|
else
|
|
assert(cc2 == CellCategory::Unknown || cc == cc2);
|
|
if (cc_new == CellCategory::Unknown) {
|
|
VertexCategory vc = vertex_category(edge.vertex0());
|
|
assert(vc != VertexCategory::OnContour);
|
|
if (vc != VertexCategory::Unknown)
|
|
cc_new = (vc == VertexCategory::Outside) ? CellCategory::Outside : CellCategory::Inside;
|
|
}
|
|
if (cc_new != CellCategory::Unknown) {
|
|
VertexCategory vc = (cc_new == CellCategory::Outside) ? VertexCategory::Outside : VertexCategory::Inside;
|
|
annotate_vertex(edge.vertex0(), vc);
|
|
annotate_vertex(edge.vertex1(), vc);
|
|
auto ec_new = (cc_new == CellCategory::Outside) ? EdgeCategory::PointsOutside : EdgeCategory::PointsInside;
|
|
annotate_edge(&edge, ec_new);
|
|
annotate_edge(edge.twin(), ec_new);
|
|
if (cc != cc_new) {
|
|
annotate_cell(&cell, cc_new);
|
|
cell_queue.emplace_back(&cell);
|
|
}
|
|
if (cc2 != cc_new) {
|
|
annotate_cell(&cell2, cc_new);
|
|
cell_queue.emplace_back(&cell2);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
assert(debug::verify_vertices_on_contour(vd, lines));
|
|
|
|
// Do a final seed fill over Voronoi cells and unmarked Voronoi edges.
|
|
while (! cell_queue.empty()) {
|
|
const VD::cell_type *cell = cell_queue.back();
|
|
const CellCategory cc = cell_category(cell);
|
|
assert(cc == CellCategory::Outside || cc == CellCategory::Inside);
|
|
cell_queue.pop_back();
|
|
const VD::edge_type *first_edge = cell->incident_edge();
|
|
const VD::edge_type *edge = first_edge;
|
|
const auto ec_new = (cc == CellCategory::Outside) ? EdgeCategory::PointsOutside : EdgeCategory::PointsInside;
|
|
do {
|
|
EdgeCategory ec = edge_category(edge);
|
|
if (ec == EdgeCategory::Unknown) {
|
|
assert(edge->cell()->contains_point() && edge->twin()->cell()->contains_point());
|
|
annotate_edge(edge, ec_new);
|
|
annotate_edge(edge->twin(), ec_new);
|
|
const VD::cell_type *cell2 = edge->twin()->cell();
|
|
CellCategory cc2 = cell_category(cell2);
|
|
assert(cc2 == CellCategory::Unknown || cc2 == cc);
|
|
if (cc2 != cc) {
|
|
annotate_cell(cell2, cc);
|
|
cell_queue.emplace_back(cell2);
|
|
}
|
|
} else {
|
|
assert(edge->vertex0() == nullptr || vertex_category(edge->vertex0()) != VertexCategory::Unknown);
|
|
assert(edge->vertex1() == nullptr || vertex_category(edge->vertex1()) != VertexCategory::Unknown);
|
|
assert(edge_category(edge->twin()) != EdgeCategory::Unknown);
|
|
assert(cell_category(edge->cell()) != CellCategory::Unknown);
|
|
assert(cell_category(edge->twin()->cell()) != CellCategory::Unknown);
|
|
}
|
|
edge = edge->next();
|
|
} while (edge != first_edge);
|
|
}
|
|
|
|
assert(debug::verify_vertices_on_contour(vd, lines));
|
|
assert(debug::verify_inside_outside_annotations(vd));
|
|
}
|
|
|
|
std::vector<double> signed_vertex_distances(const VD &vd, const Lines &lines)
|
|
{
|
|
// vd shall be annotated.
|
|
assert(debug::verify_inside_outside_annotations(vd));
|
|
|
|
std::vector<double> out(vd.vertices().size(), 0.);
|
|
const VD::vertex_type *first_vertex = &vd.vertices().front();
|
|
for (const VD::vertex_type &vertex : vd.vertices()) {
|
|
const VertexCategory vc = vertex_category(vertex);
|
|
double dist;
|
|
if (vc == VertexCategory::OnContour) {
|
|
dist = 0.;
|
|
} else {
|
|
const VD::edge_type *first_edge = vertex.incident_edge();
|
|
const VD::edge_type *edge = first_edge;
|
|
const VD::cell_type *point_cell = nullptr;
|
|
do {
|
|
if (edge->cell()->contains_point()) {
|
|
point_cell = edge->cell();
|
|
break;
|
|
}
|
|
edge = edge->rot_next();
|
|
} while (edge != first_edge);
|
|
if (point_cell == 0) {
|
|
// Project vertex onto a contour segment.
|
|
const Line &line = lines[edge->cell()->source_index()];
|
|
dist = Geometry::ray_point_distance<Vec2d>(
|
|
line.a.cast<double>(), (line.b - line.a).cast<double>(), vertex_point(vertex));
|
|
} else {
|
|
// Distance to a contour point.
|
|
dist = (contour_point(*point_cell, lines).cast<double>() - vertex_point(vertex)).norm();
|
|
}
|
|
if (vc == VertexCategory::Inside)
|
|
dist = - dist;
|
|
}
|
|
out[&vertex - first_vertex] = dist;
|
|
}
|
|
|
|
assert(debug::verify_signed_distances(vd, lines, out));
|
|
|
|
return out;
|
|
}
|
|
|
|
std::vector<Vec2d> edge_offset_contour_intersections(
|
|
const VD &vd,
|
|
const Lines &lines,
|
|
const std::vector<double> &vertex_distances,
|
|
double offset_distance)
|
|
{
|
|
// vd shall be annotated.
|
|
assert(debug::verify_inside_outside_annotations(vd));
|
|
|
|
bool outside = offset_distance > 0;
|
|
if (! outside)
|
|
offset_distance = - offset_distance;
|
|
assert(offset_distance > 0.);
|
|
|
|
const VD::vertex_type *first_vertex = &vd.vertices().front();
|
|
const VD::edge_type *first_edge = &vd.edges().front();
|
|
static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
|
|
// By default none edge has an intersection with the offset curve.
|
|
std::vector<Vec2d> out(vd.num_edges(), Vec2d(nan, 0.));
|
|
|
|
for (const VD::edge_type &edge : vd.edges()) {
|
|
size_t edge_idx = &edge - first_edge;
|
|
if (edge_offset_has_intersection(out[edge_idx]) || out[edge_idx].y() != 0.)
|
|
// This edge was already classified.
|
|
continue;
|
|
|
|
const VD::vertex_type *v0 = edge.vertex0();
|
|
const VD::vertex_type *v1 = edge.vertex1();
|
|
if (v0 == nullptr) {
|
|
assert(vertex_category(v1) == VertexCategory::OnContour || vertex_category(v1) == VertexCategory::Outside);
|
|
continue;
|
|
}
|
|
|
|
double d0 = (v0 == nullptr) ? std::numeric_limits<double>::max() : vertex_distances[v0 - first_vertex];
|
|
double d1 = (v1 == nullptr) ? std::numeric_limits<double>::max() : vertex_distances[v1 - first_vertex];
|
|
assert(d0 * d1 >= 0.);
|
|
if (! outside) {
|
|
d0 = - d0;
|
|
d1 = - d1;
|
|
}
|
|
#ifndef NDEBUG
|
|
{
|
|
double err = std::abs(detail::dist_to_site(lines, *edge.cell(), vertex_point(v0)) - std::abs(d0));
|
|
double err2 = std::abs(detail::dist_to_site(lines, *edge.twin()->cell(), vertex_point(v0)) - std::abs(d0));
|
|
assert(err < SCALED_EPSILON);
|
|
assert(err2 < SCALED_EPSILON);
|
|
if (v1 != nullptr) {
|
|
double err3 = std::abs(detail::dist_to_site(lines, *edge.cell(), vertex_point(v1)) - std::abs(d1));
|
|
double err4 = std::abs(detail::dist_to_site(lines, *edge.twin()->cell(), vertex_point(v1)) - std::abs(d1));
|
|
assert(err3 < SCALED_EPSILON);
|
|
assert(err4 < SCALED_EPSILON);
|
|
}
|
|
}
|
|
#endif // NDEBUG
|
|
double dmin, dmax;
|
|
if (d0 < d1)
|
|
dmin = d0, dmax = d1;
|
|
else
|
|
dmax = d0, dmin = d1;
|
|
// Offset distance may be lower than dmin, but never higher than dmax.
|
|
// Don't intersect an edge at dmax
|
|
// 1) To avoid zero edge length, zero area offset contours.
|
|
// 2) To ensure that the offset contours that cross a Voronoi vertex are traced consistently
|
|
// at one side of the offset curve only.
|
|
if (offset_distance >= dmax)
|
|
continue;
|
|
|
|
// Edge candidate, intersection points were not calculated yet.
|
|
assert(v0 != nullptr);
|
|
const VD::cell_type *cell = edge.cell();
|
|
const VD::cell_type *cell2 = edge.twin()->cell();
|
|
const Line &line0 = lines[cell->source_index()];
|
|
const Line &line1 = lines[cell2->source_index()];
|
|
size_t edge_idx2 = edge.twin() - first_edge;
|
|
if (v1 == nullptr) {
|
|
assert(edge.is_infinite());
|
|
assert(edge.is_linear());
|
|
// Unconstrained edges have always montonous distance.
|
|
assert(d0 != d1);
|
|
if (offset_distance > dmin) {
|
|
// There is certainly an intersection with the offset curve.
|
|
if (cell->contains_point() && cell2->contains_point()) {
|
|
assert(! edge.is_secondary());
|
|
const Point &pt0 = contour_point(*cell, line0);
|
|
const Point &pt1 = contour_point(*cell2, line1);
|
|
// pt is inside the circle (pt0, offset_distance), (pt + dir) is certainly outside the same circle.
|
|
Vec2d dir = Vec2d(double(pt0.y() - pt1.y()), double(pt1.x() - pt0.x())) * (2. * offset_distance);
|
|
Vec2d pt(v0->x(), v0->y());
|
|
double t = detail::first_circle_segment_intersection_parameter(Vec2d(pt0.x(), pt0.y()), offset_distance, pt, dir);
|
|
assert(t > 0.);
|
|
out[edge_idx] = pt + t * dir;
|
|
} else {
|
|
// Infinite edges could not be created by two segment sites.
|
|
assert(cell->contains_point() != cell2->contains_point());
|
|
// Linear edge goes through the endpoint of a segment.
|
|
assert(edge.is_secondary());
|
|
const Point &ipt = cell->contains_segment() ? contour_point(*cell2, line1) : contour_point(*cell, line0);
|
|
#ifndef NDEBUG
|
|
if (cell->contains_segment()) {
|
|
const Point &pt1 = contour_point(*cell2, line1);
|
|
assert(pt1 == line0.a || pt1 == line0.b);
|
|
} else {
|
|
const Point &pt0 = contour_point(*cell, line0);
|
|
assert(pt0 == line1.a || pt0 == line1.b);
|
|
}
|
|
assert((vertex_point(v0) - ipt.cast<double>()).norm() < SCALED_EPSILON);
|
|
#endif /* NDEBUG */
|
|
// Infinite edge starts at an input contour, therefore there is always an intersection with an offset curve.
|
|
const Line &line = cell->contains_segment() ? line0 : line1;
|
|
assert(line.a == ipt || line.b == ipt);
|
|
out[edge_idx] = ipt.cast<double>() + offset_distance * Vec2d(line.b.y() - line.a.y(), line.a.x() - line.b.x()).normalized();
|
|
}
|
|
} else if (offset_distance == dmin)
|
|
out[edge_idx] = vertex_point(v0);
|
|
// The other edge of an unconstrained edge starting with null vertex shall never be intersected. Mark it as visited.
|
|
out[edge_idx2].y() = 1.;
|
|
} else {
|
|
assert(edge.is_finite());
|
|
bool done = false;
|
|
// Bisector of two line segments, distance along the bisector is linear.
|
|
bool bisector = cell->contains_segment() && cell2->contains_segment();
|
|
assert(edge.is_finite());
|
|
// or a secondary line, again the distance along the secondary line is linear and starts at the contour (zero distance).
|
|
if (bisector || edge.is_secondary()) {
|
|
assert(edge.is_linear());
|
|
#ifndef NDEBUG
|
|
if (edge.is_secondary()) {
|
|
assert(cell->contains_point() != cell2->contains_point());
|
|
// One of the vertices is on the input contour.
|
|
assert((vertex_category(edge.vertex0()) == VertexCategory::OnContour) !=
|
|
(vertex_category(edge.vertex1()) == VertexCategory::OnContour));
|
|
assert(dmin == 0.);
|
|
}
|
|
#endif // NDEBUG
|
|
if (! bisector || (dmin != dmax && offset_distance >= dmin)) {
|
|
double t = (offset_distance - dmin) / (dmax - dmin);
|
|
t = std::clamp(t, 0., 1.);
|
|
if (d1 < d0) {
|
|
out[edge_idx2] = Slic3r::lerp(vertex_point(v1), vertex_point(v0), t);
|
|
// mark visited
|
|
out[edge_idx].y() = 1.;
|
|
} else {
|
|
out[edge_idx] = Slic3r::lerp(vertex_point(v0), vertex_point(v1), t);
|
|
// mark visited
|
|
out[edge_idx2].y() = 1.;
|
|
}
|
|
done = true;
|
|
}
|
|
} else {
|
|
// Point - Segment or Point - Point edge, distance along this Voronoi edge may not be monotonous:
|
|
// The distance function along the Voronoi edge may either have one maximum at one vertex and one minimum at the other vertex,
|
|
// or it may have two maxima at the vertices and a minimum somewhere along the Voronoi edge, and this Voronoi edge
|
|
// may be intersected twice by an offset curve.
|
|
//
|
|
// Tracing an offset curve accross Voronoi regions with linear edges of montonously increasing or decrasing distance
|
|
// to a Voronoi region is stable in a sense, that if the distance of Voronoi vertices is calculated correctly, there will
|
|
// be maximum one intersection of an offset curve found at each Voronoi edge and tracing these intersections shall
|
|
// produce a set of closed curves.
|
|
//
|
|
// Having a non-monotonous distance function between the Voronoi edge end points may lead to splitting of offset curves
|
|
// at these Voronoi edges. If a Voronoi edge is classified as having no intersection at all while it has some,
|
|
// the extracted offset curve will contain self intersections at this Voronoi edge.
|
|
//
|
|
// If on the other side the two intersection points are found by a numerical error even though none should be found, then
|
|
// it may happen that it would not be possible to connect these two points into a closed loop, which is likely worse
|
|
// than the issue above.
|
|
//
|
|
// While it is likely not possible to avoid all the numerical issues, one shall strive for the highest numerical robustness.
|
|
assert(cell->contains_point() || cell2->contains_point());
|
|
size_t num_intersections = 0;
|
|
bool point_vs_segment = cell->contains_point() != cell2->contains_point();
|
|
const Point &pt0 = cell->contains_point() ? contour_point(*cell, line0) : contour_point(*cell2, line1);
|
|
// Project p0 to line segment <v0, v1>.
|
|
Vec2d p0(v0->x(), v0->y());
|
|
Vec2d p1(v1->x(), v1->y());
|
|
Vec2d px(pt0.x(), pt0.y());
|
|
const Point *pt1 = nullptr;
|
|
Vec2d dir;
|
|
if (point_vs_segment) {
|
|
const Line &line = cell->contains_segment() ? line0 : line1;
|
|
dir = (line.b - line.a).cast<double>();
|
|
} else {
|
|
pt1 = &contour_point(*cell2, line1);
|
|
// Perpendicular to the (pt1 - pt0) direction.
|
|
dir = Vec2d(double(pt0.y() - pt1->y()), double(pt1->x() - pt0.x()));
|
|
}
|
|
double s0 = (p0 - px).dot(dir);
|
|
double s1 = (p1 - px).dot(dir);
|
|
if (offset_distance >= dmin) {
|
|
// This Voronoi edge is intersected by the offset curve just once.
|
|
// There may be numerical issues if dmin is close to the minimum of the non-monotonous distance function.
|
|
num_intersections = 1;
|
|
} else {
|
|
// This Voronoi edge may not be intersected by the offset curve, or it may split the offset curve
|
|
// into two loops. First classify this edge robustly with regard to the Point-Segment bisector or Point-Point bisector.
|
|
double dmin_new;
|
|
bool found = false;
|
|
if (point_vs_segment) {
|
|
if (s0 * s1 <= 0.) {
|
|
// One end of the Voronoi edge is on one side of the Point-Segment bisector, the other end of the Voronoi
|
|
// edge is on the other side of the bisector, therefore with a high probability we should find a minimum
|
|
// of the distance to a nearest site somewhere inside this Voronoi edge (at the intersection of the bisector
|
|
// and the Voronoi edge.
|
|
const Line &line = cell->contains_segment() ? line0 : line1;
|
|
dmin_new = 0.5 * (Geometry::foot_pt<Vec2d>(line.a.cast<double>(), dir, px) - px).norm();
|
|
found = true;
|
|
}
|
|
} else {
|
|
// Point-Point Voronoi sites.
|
|
if (s0 * s1 <= 0.) {
|
|
// This Voronoi edge intersects connection line of the two Point sites.
|
|
dmin_new = 0.5 * (pt1->cast<double>() - px).norm();
|
|
found = true;
|
|
}
|
|
}
|
|
if (found) {
|
|
assert(dmin_new < dmax + SCALED_EPSILON);
|
|
assert(dmin_new < dmin + SCALED_EPSILON);
|
|
if (dmin_new <= offset_distance) {
|
|
// 1) offset_distance > dmin_new -> two new distinct intersection points are found.
|
|
// 2) offset_distance == dmin_new -> one duplicate point is found.
|
|
// If 2) is ignored, then two tangentially touching offset curves are created.
|
|
// If not ignored, then the two offset curves merge at this double point.
|
|
// We should merge the contours while pushing the the two copies of the tangent point away a bit.
|
|
dmin = dmin_new;
|
|
num_intersections = (offset_distance > dmin) + 1;
|
|
}
|
|
}
|
|
}
|
|
if (num_intersections > 0) {
|
|
detail::Intersections intersections;
|
|
if (point_vs_segment) {
|
|
assert(cell->contains_point() || cell2->contains_point());
|
|
intersections = detail::line_point_equal_distance_points(cell->contains_segment() ? line0 : line1, pt0, offset_distance);
|
|
} else {
|
|
intersections = detail::point_point_equal_distance_points(pt0, *pt1, offset_distance);
|
|
}
|
|
// The functions above perform complex calculations in doubles, therefore the results may not be quite accurate and
|
|
// the number of intersections found may not be in accord to the number of intersections expected from evaluating simpler expressions.
|
|
// Adjust the result to the number of intersection points expected.
|
|
if (num_intersections == 2) {
|
|
switch (intersections.count) {
|
|
case 0:
|
|
// No intersection found even though one or two were expected to be found.
|
|
// Not trying to find the intersection means that we may produce offset curves, that intersect at this Voronoi edge.
|
|
//FIXME We are fine with that for now, but we may try to create artificial split points in further revisions.
|
|
break;
|
|
case 1:
|
|
// Tangential point found.
|
|
//FIXME We are fine with that for now, but we may try to create artificial split points in further revisions.
|
|
break;
|
|
default:
|
|
{
|
|
// Two intersection points found. Sort them.
|
|
assert(intersections.count == 2);
|
|
double q0 = (intersections.pts[0] - px).dot(dir);
|
|
double q1 = (intersections.pts[1] - px).dot(dir);
|
|
// Both Voronoi edge end points and offset contour intersection points should be separated by the bisector.
|
|
assert(q0 * q1 <= 0.);
|
|
assert(s0 * s1 <= 0.);
|
|
// Sort the intersection points by dir.
|
|
if ((q0 < q1) != (s0 < s1))
|
|
std::swap(intersections.pts[0], intersections.pts[1]);
|
|
}
|
|
}
|
|
} else {
|
|
assert(num_intersections == 1);
|
|
switch (intersections.count) {
|
|
case 0:
|
|
// No intersection found. This should not happen.
|
|
// Create one artificial intersection point by repeating the dmin point, which is supposed to be
|
|
// close to the minimum.
|
|
intersections.pts[0] = (dmin == d0) ? p0 : p1;
|
|
intersections.count = 1;
|
|
break;
|
|
case 1:
|
|
// One intersection found. This is a tangential point. Use it.
|
|
break;
|
|
default:
|
|
// Two intersections found.
|
|
// Now decide which of the point fall on this Voronoi edge.
|
|
assert(intersections.count == 2);
|
|
double q0 = (intersections.pts[0] - px).dot(dir);
|
|
double q1 = (intersections.pts[1] - px).dot(dir);
|
|
// Offset contour intersection points should be separated by the bisector.
|
|
assert(q0 * q1 <= 0);
|
|
double s = (dmax == d0) ? s0 : s1;
|
|
bool take_2nd = (s > 0.) ? q1 > q0 : q1 < q0;
|
|
if (take_2nd)
|
|
intersections.pts[0] = intersections.pts[1];
|
|
-- intersections.count;
|
|
}
|
|
}
|
|
assert(intersections.count > 0);
|
|
if (intersections.count == 2) {
|
|
out[edge_idx] = intersections.pts[1];
|
|
out[edge_idx2] = intersections.pts[0];
|
|
done = true;
|
|
} else if (intersections.count == 1) {
|
|
if (d1 < d0)
|
|
std::swap(edge_idx, edge_idx2);
|
|
out[edge_idx] = intersections.pts[0];
|
|
out[edge_idx2].y() = 1.;
|
|
done = true;
|
|
}
|
|
}
|
|
}
|
|
if (! done)
|
|
out[edge_idx].y() = out[edge_idx2].y() = 1.;
|
|
}
|
|
}
|
|
|
|
assert(debug::verify_offset_intersection_points(vd, lines, offset_distance, out));
|
|
|
|
return out;
|
|
}
|
|
|
|
Polygons offset(
|
|
const Geometry::VoronoiDiagram &vd,
|
|
const Lines &lines,
|
|
const std::vector<double> &signed_vertex_distances,
|
|
double offset_distance,
|
|
double discretization_error)
|
|
{
|
|
#ifdef VORONOI_DEBUG_OUT
|
|
BoundingBox bbox;
|
|
{
|
|
bbox.merge(get_extents(lines));
|
|
bbox.min -= (0.01 * bbox.size().cast<double>()).cast<coord_t>();
|
|
bbox.max += (0.01 * bbox.size().cast<double>()).cast<coord_t>();
|
|
}
|
|
static int irun = 0;
|
|
++ irun;
|
|
{
|
|
Lines helper_lines;
|
|
for (const VD::edge_type &edge : vd.edges())
|
|
if (edge_category(edge) == (offset_distance > 0 ? EdgeCategory::PointsOutside : EdgeCategory::PointsInside) &&
|
|
edge.vertex0() != nullptr) {
|
|
const VD::vertex_type *v0 = edge.vertex0();
|
|
const VD::vertex_type *v1 = edge.vertex1();
|
|
Vec2d pt1(v0->x(), v0->y());
|
|
Vec2d pt2;
|
|
if (v1 == nullptr) {
|
|
// Unconstrained edge. Calculate a trimmed position.
|
|
assert(edge.is_linear());
|
|
const VD::cell_type *cell = edge.cell();
|
|
const VD::cell_type *cell2 = edge.twin()->cell();
|
|
const Line &line0 = lines[cell->source_index()];
|
|
const Line &line1 = lines[cell2->source_index()];
|
|
if (cell->contains_point() && cell2->contains_point()) {
|
|
const Point &pt0 = (cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b;
|
|
const Point &pt1 = (cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b;
|
|
// Direction vector of this unconstrained Voronoi edge.
|
|
Vec2d dir(double(pt0.y() - pt1.y()), double(pt1.x() - pt0.x()));
|
|
pt2 = Vec2d(v0->x(), v0->y()) + dir.normalized() * scale_(10.);
|
|
} else {
|
|
// Infinite edges could not be created by two segment sites.
|
|
assert(cell->contains_point() != cell2->contains_point());
|
|
// Linear edge goes through the endpoint of a segment.
|
|
assert(edge.is_secondary());
|
|
const Point &ipt = cell->contains_segment() ?
|
|
((cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b) :
|
|
((cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b);
|
|
// Infinite edge starts at an input contour, therefore there is always an intersection with an offset curve.
|
|
const Line &line = cell->contains_segment() ? line0 : line1;
|
|
assert(line.a == ipt || line.b == ipt);
|
|
// dir is perpendicular to line.
|
|
Vec2d dir(line.a.y() - line.b.y(), line.b.x() - line.a.x());
|
|
assert(dir.norm() > 0.);
|
|
if (((line.a == ipt) == cell->contains_point()) == (v0 == nullptr))
|
|
dir = - dir;
|
|
pt2 = ipt.cast<double>() + dir.normalized() * scale_(10.);
|
|
}
|
|
} else {
|
|
pt2 = Vec2d(v1->x(), v1->y());
|
|
// Clip the line by the bounding box, so that the coloring of the line will be visible.
|
|
Geometry::liang_barsky_line_clipping(pt1, pt2, BoundingBoxf(bbox.min.cast<double>(), bbox.max.cast<double>()));
|
|
}
|
|
helper_lines.emplace_back(Line(Point(pt1.cast<coord_t>()), Point(((pt1 + pt2) * 0.5).cast<coord_t>())));
|
|
}
|
|
dump_voronoi_to_svg(debug_out_path("voronoi-offset-candidates1-%d.svg", irun).c_str(), vd, Points(), lines, Polygons(), helper_lines);
|
|
}
|
|
#endif // VORONOI_DEBUG_OUT
|
|
|
|
std::vector<Vec2d> edge_points = edge_offset_contour_intersections(vd, lines, signed_vertex_distances, offset_distance);
|
|
const VD::edge_type *front_edge = &vd.edges().front();
|
|
|
|
#ifdef VORONOI_DEBUG_OUT
|
|
Lines helper_lines;
|
|
{
|
|
for (const VD::edge_type &edge : vd.edges())
|
|
if (edge_offset_has_intersection(edge_points[&edge - front_edge]))
|
|
helper_lines.emplace_back(Line(Point(edge.vertex0()->x(), edge.vertex0()->y()), Point(edge_points[&edge - front_edge].cast<coord_t>())));
|
|
dump_voronoi_to_svg(debug_out_path("voronoi-offset-candidates2-%d.svg", irun).c_str(), vd, Points(), lines, Polygons(), helper_lines);
|
|
}
|
|
#endif // VORONOI_DEBUG_OUT
|
|
|
|
auto next_offset_edge = [&edge_points, front_edge](const VD::edge_type *start_edge) -> const VD::edge_type* {
|
|
for (const VD::edge_type *edge = start_edge->next(); edge != start_edge; edge = edge->next())
|
|
if (edge_offset_has_intersection(edge_points[edge->twin() - front_edge]))
|
|
return edge->twin();
|
|
// assert(false);
|
|
return nullptr;
|
|
};
|
|
|
|
const bool inside_offset = offset_distance < 0.;
|
|
if (inside_offset)
|
|
offset_distance = - offset_distance;
|
|
|
|
// Track the offset curves.
|
|
Polygons out;
|
|
double angle_step = 2. * acos((offset_distance - discretization_error) / offset_distance);
|
|
double cos_threshold = cos(angle_step);
|
|
static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
|
|
for (size_t seed_edge_idx = 0; seed_edge_idx < vd.num_edges(); ++ seed_edge_idx) {
|
|
Vec2d last_pt = edge_points[seed_edge_idx];
|
|
if (edge_offset_has_intersection(last_pt)) {
|
|
const VD::edge_type *start_edge = &vd.edges()[seed_edge_idx];
|
|
const VD::edge_type *edge = start_edge;
|
|
Polygon poly;
|
|
do {
|
|
// find the next edge
|
|
const VD::edge_type *next_edge = next_offset_edge(edge);
|
|
#ifdef VORONOI_DEBUG_OUT
|
|
if (next_edge == nullptr) {
|
|
Lines hl = helper_lines;
|
|
append(hl, to_lines(Polyline(poly.points)));
|
|
dump_voronoi_to_svg(debug_out_path("voronoi-offset-open-loop-%d.svg", irun).c_str(), vd, Points(), lines, Polygons(), hl);
|
|
}
|
|
#endif // VORONOI_DEBUG_OUT
|
|
assert(next_edge);
|
|
//std::cout << "offset-output: "; print_edge(edge); std::cout << " to "; print_edge(next_edge); std::cout << "\n";
|
|
// Interpolate a circular segment or insert a linear segment between edge and next_edge.
|
|
const VD::cell_type *cell = edge->cell();
|
|
// Mark the edge / offset curve intersection point as consumed.
|
|
Vec2d p1 = last_pt;
|
|
Vec2d p2 = edge_points[next_edge - front_edge];
|
|
edge_points[next_edge - front_edge].x() = nan;
|
|
#ifndef NDEBUG
|
|
{
|
|
double err = detail::dist_to_site(lines, *cell, p1) - offset_distance;
|
|
double err2 = detail::dist_to_site(lines, *cell, p2) - offset_distance;
|
|
#ifdef VORONOI_DEBUG_OUT
|
|
if (std::max(err, err2) >= SCALED_EPSILON) {
|
|
Lines helper_lines;
|
|
dump_voronoi_to_svg(debug_out_path("voronoi-offset-incorrect_pt-%d.svg", irun).c_str(), vd, Points(), lines, Polygons(), to_lines(poly));
|
|
}
|
|
#endif // VORONOI_DEBUG_OUT
|
|
assert(std::abs(err) < SCALED_EPSILON);
|
|
assert(std::abs(err2) < SCALED_EPSILON);
|
|
}
|
|
#endif /* NDEBUG */
|
|
if (cell->contains_point()) {
|
|
// Discretize an arc from p1 to p2 with radius = offset_distance and discretization_error.
|
|
// The extracted contour is CCW oriented, extracted holes are CW oriented.
|
|
// The extracted arc will have the same orientation. As the Voronoi regions are convex, the angle covered by the arc will be convex as well.
|
|
const Line &line0 = lines[cell->source_index()];
|
|
const Vec2d ¢er = ((cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b).cast<double>();
|
|
const Vec2d v1 = p1 - center;
|
|
const Vec2d v2 = p2 - center;
|
|
bool ccw = cross2(v1, v2) > 0;
|
|
double cos_a = v1.dot(v2);
|
|
double norm = v1.norm() * v2.norm();
|
|
assert(norm > 0.);
|
|
if (cos_a < cos_threshold * norm) {
|
|
// Angle is bigger than the threshold, therefore the arc will be discretized.
|
|
cos_a /= norm;
|
|
assert(cos_a > -1. - EPSILON && cos_a < 1. + EPSILON);
|
|
double angle = acos(std::max(-1., std::min(1., cos_a)));
|
|
size_t n_steps = size_t(ceil(angle / angle_step));
|
|
double astep = angle / n_steps;
|
|
if (! ccw)
|
|
astep *= -1.;
|
|
double a = astep;
|
|
for (size_t i = 1; i < n_steps; ++ i, a += astep) {
|
|
double c = cos(a);
|
|
double s = sin(a);
|
|
Vec2d p = center + Vec2d(c * v1.x() - s * v1.y(), s * v1.x() + c * v1.y());
|
|
poly.points.emplace_back(Point(coord_t(p.x()), coord_t(p.y())));
|
|
}
|
|
}
|
|
}
|
|
{
|
|
Point pt_last(coord_t(p2.x()), coord_t(p2.y()));
|
|
if (poly.empty() || poly.points.back() != pt_last)
|
|
poly.points.emplace_back(pt_last);
|
|
}
|
|
edge = next_edge;
|
|
last_pt = p2;
|
|
} while (edge != start_edge);
|
|
|
|
while (! poly.empty() && poly.points.front() == poly.points.back())
|
|
poly.points.pop_back();
|
|
if (poly.size() >= 3)
|
|
out.emplace_back(std::move(poly));
|
|
}
|
|
}
|
|
|
|
return out;
|
|
}
|
|
|
|
Polygons offset(
|
|
const VD &vd,
|
|
const Lines &lines,
|
|
double offset_distance,
|
|
double discretization_error)
|
|
{
|
|
annotate_inside_outside(const_cast<VD&>(vd), lines);
|
|
std::vector<double> dist = signed_vertex_distances(vd, lines);
|
|
return offset(vd, lines, dist, offset_distance, discretization_error);
|
|
}
|
|
|
|
// Produce a list of start positions of a skeleton segment at a halfedge.
|
|
// If the whole Voronoi edge is part of the skeleton, then zero start positions are assigned
|
|
// to both ends of the edge. Position "1" shall never be assigned to a halfedge.
|
|
//
|
|
// Skeleton edges must be inside a closed polygon, therefore these edges are finite.
|
|
// A Voronoi Edge-Edge bisector is either completely part of a skeleton, or not at all.
|
|
// An infinite Voronoi Edge-Point (parabola) or Point-Point (line) bisector is split into
|
|
// a center part close to the Voronoi sites (not skeleton) and the ends (skeleton),
|
|
// though either part could be clipped by the Voronoi segment.
|
|
//
|
|
// Further filtering of the skeleton may be necessary.
|
|
std::vector<Vec2d> skeleton_edges_rough(
|
|
const VD &vd,
|
|
const Lines &lines,
|
|
// Angle threshold at a sharp convex corner, which is marked for a gap fill.
|
|
const double threshold_alpha)
|
|
{
|
|
// vd shall be annotated.
|
|
assert(debug::verify_inside_outside_annotations(vd));
|
|
|
|
const VD::edge_type *first_edge = &vd.edges().front();
|
|
static constexpr double nan = std::numeric_limits<double>::quiet_NaN();
|
|
// By default no edge is annotated as being part of the skeleton.
|
|
std::vector<Vec2d> out(vd.num_edges(), Vec2d(nan, nan));
|
|
// Threshold at a sharp corner, derived from a dot product of the sharp corner edges.
|
|
const double threshold_cos_alpha = cos(threshold_alpha);
|
|
// For sharp corners, dr/dl = sin(alpha/2). Substituting the dr/dl threshold with tan(alpha/2) threshold
|
|
// in Voronoi point - point and Voronoi point - line site functions.
|
|
const double threshold_tan_alpha_half = tan(0.5 * threshold_alpha);
|
|
|
|
for (const VD::edge_type &edge : vd.edges()) {
|
|
size_t edge_idx = &edge - first_edge;
|
|
if (
|
|
// Ignore secondary and unbounded edges, they shall never be part of the skeleton.
|
|
edge.is_secondary() || edge.is_infinite() ||
|
|
// Skip the twin edge of an edge, that has already been processed.
|
|
&edge > edge.twin() ||
|
|
// Ignore outer edges.
|
|
(edge_category(edge) != EdgeCategory::PointsInside && edge_category(edge.twin()) != EdgeCategory::PointsInside))
|
|
continue;
|
|
const VD::vertex_type *v0 = edge.vertex0();
|
|
const VD::vertex_type *v1 = edge.vertex1();
|
|
const VD::cell_type *cell = edge.cell();
|
|
const VD::cell_type *cell2 = edge.twin()->cell();
|
|
const Line &line0 = lines[cell->source_index()];
|
|
const Line &line1 = lines[cell2->source_index()];
|
|
size_t edge_idx2 = edge.twin() - first_edge;
|
|
if (cell->contains_segment() && cell2->contains_segment()) {
|
|
// Bisector of two line segments, distance along the bisector is linear,
|
|
// dr/dl is constant.
|
|
// using trigonometric identity sin^2(a) = (1-cos(2a))/2
|
|
Vec2d lv0 = (line0.b - line0.a).cast<double>();
|
|
Vec2d lv1 = (line1.b - line1.a).cast<double>();
|
|
double d = lv0.dot(lv1);
|
|
if (d < 0.) {
|
|
double cos_alpha = - d / (lv0.norm() * lv1.norm());
|
|
if (cos_alpha > threshold_cos_alpha) {
|
|
// The whole bisector is a skeleton segment.
|
|
out[edge_idx] = vertex_point(v0);
|
|
out[edge_idx2] = vertex_point(v1);
|
|
}
|
|
}
|
|
} else {
|
|
// An infinite Voronoi Edge-Point (parabola) or Point-Point (line) bisector, clipped to a finite Voronoi segment.
|
|
// The infinite bisector has a distance (skeleton radius) minimum, which is also a minimum
|
|
// of the skeleton function dr / dt.
|
|
assert(cell->contains_point() || cell2->contains_point());
|
|
if (cell->contains_point() != cell2->contains_point()) {
|
|
// Point - Segment
|
|
const Point &pt0 = cell->contains_point() ? contour_point(*cell, line0) : contour_point(*cell2, line1);
|
|
const Line &line = cell->contains_segment() ? line0 : line1;
|
|
std::tie(out[edge_idx], out[edge_idx2]) = detail::point_segment_dr_dl_thresholds(
|
|
pt0, line, vertex_point(v0), vertex_point(v1), threshold_tan_alpha_half);
|
|
} else {
|
|
// Point - Point
|
|
const Point &pt0 = contour_point(*cell, line0);
|
|
const Point &pt1 = contour_point(*cell2, line1);
|
|
std::tie(out[edge_idx], out[edge_idx2]) = detail::point_point_dr_dl_thresholds(
|
|
pt0, pt1, vertex_point(v0), vertex_point(v1), threshold_tan_alpha_half);
|
|
}
|
|
}
|
|
}
|
|
|
|
#ifdef VORONOI_DEBUG_OUT
|
|
{
|
|
static int irun = 0;
|
|
++ irun;
|
|
Lines helper_lines;
|
|
for (const VD::edge_type &edge : vd.edges())
|
|
if (&edge < edge.twin() && edge.is_finite()) {
|
|
const Vec2d &skeleton_pt = out[&edge - first_edge];
|
|
const Vec2d &skeleton_pt2 = out[edge.twin() - first_edge];
|
|
bool has_skeleton_pt = ! std::isnan(skeleton_pt.x());
|
|
bool has_skeleton_pt2 = ! std::isnan(skeleton_pt2.x());
|
|
const Vec2d &vertex_pt = vertex_point(edge.vertex0());
|
|
const Vec2d &vertex_pt2 = vertex_point(edge.vertex1());
|
|
if (has_skeleton_pt && has_skeleton_pt2) {
|
|
// Complete edge is part of the skeleton.
|
|
helper_lines.emplace_back(Line(Point(vertex_pt.x(), vertex_pt.y()), Point(vertex_pt2.x(), vertex_pt2.y())));
|
|
} else {
|
|
if (has_skeleton_pt)
|
|
helper_lines.emplace_back(Line(Point(vertex_pt2.x(), vertex_pt2.y()), Point(skeleton_pt.x(), skeleton_pt.y())));
|
|
if (has_skeleton_pt2)
|
|
helper_lines.emplace_back(Line(Point(vertex_pt.x(), vertex_pt.y()), Point(skeleton_pt2.x(), skeleton_pt2.y())));
|
|
}
|
|
}
|
|
dump_voronoi_to_svg(debug_out_path("voronoi-skeleton-edges-%d.svg", irun).c_str(), vd, Points(), lines, Polygons(), helper_lines);
|
|
}
|
|
#endif // VORONOI_DEBUG_OUT
|
|
|
|
return out;
|
|
}
|
|
|
|
} // namespace Voronoi
|
|
} // namespace Slic3r
|