3D-printing/OrcaSlicer
1
0
Fork 0
mirror of https://github.com/SoftFever/OrcaSlicer.git synced 2025-10-31 04:31:15 -06:00
Code Issues Projects Releases Packages Wiki Activity Actions 5
OrcaSlicer/src/libslic3r/VoronoiOffset.cpp
Vojtech Bubnik b101a8e266 Fixes of the offset curves from Voronoi diagram. 
The offset curve extractor is already quite usable,
though singular cases are still not covered yet
when the offset curve intersects or nearly intersects
a Voronoi vertex.

Removal of the PRINTF_ZU "%zu" Visual Studio printf compatibility macro.
Fixes of a contours self intersection test for collinear segments.
SVG exporter now exports white background, so that the GNOME Eye viewer is usable.
2020-06-16 13:15:48 +02:00

855 lines
45 KiB
C++
Raw Blame History

// Polygon offsetting using Voronoi diagram prodiced by boost::polygon.
#include "VoronoiOffset.hpp"
#include <cmath>
// #define VORONOI_DEBUG_OUT
#ifdef VORONOI_DEBUG_OUT
#include <libslic3r/VoronoiVisualUtils.hpp>
#endif
namespace Slic3r {
using VD = Geometry::VoronoiDiagram;
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 &center, const double r, const Vec2d &pt, const Vec2d &v)
{
const Vec2d d = pt - center;
#ifndef NDEBUG
double d0 = (pt - center).norm();
double d1 = (pt + v - center).norm();
assert(r < std::max(d0, d1) + 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 Slic3r::clamp(0., 1., t);
} 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;
}
} // namespace detail
Polygons voronoi_offset(
const Geometry::VoronoiDiagram &vd,
const Lines &lines,
double offset_distance,
double discretization_error)
{
#ifndef NDEBUG
// Verify that twin halfedges are stored next to the other in vd.
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);
}
}
#endif // NDEBUG
enum class EdgeState : unsigned char {
// Initial state, don't know.
Unknown,
// This edge will certainly not be intersected by the offset curve.
Inactive,
// This edge will certainly be intersected by the offset curve.
Active,
// This edge will possibly be intersected by the offset curve.
Possible
};
enum class CellState : unsigned char {
// Initial state, don't know.
Unknown,
// Inactive cell is inside for outside curves and outside for inside curves.
Inactive,
// Active cell is outside for outside curves and inside for inside curves.
Active,
// Boundary cell is intersected by the input segment, part of it is active.
Boundary
};
// Mark edges with outward vertex pointing outside the polygons, thus there is a chance
// that such an edge will have an intersection with our desired offset curve.
bool outside = offset_distance > 0.;
std::vector<EdgeState> edge_state(vd.num_edges(), EdgeState::Unknown);
std::vector<CellState> cell_state(vd.num_cells(), CellState::Unknown);
const VD::edge_type *front_edge = &vd.edges().front();
const VD::cell_type *front_cell = &vd.cells().front();
auto set_edge_state_initial = [&edge_state, front_edge](const VD::edge_type *edge, EdgeState new_edge_type) {
EdgeState &edge_type = edge_state[edge - front_edge];
assert(edge_type == EdgeState::Unknown || edge_type == new_edge_type);
assert(new_edge_type == EdgeState::Possible || new_edge_type == EdgeState::Inactive);
edge_type = new_edge_type;
};
auto set_edge_state_final = [&edge_state, front_edge](const size_t edge_id, EdgeState new_edge_type) {
EdgeState &edge_type = edge_state[edge_id];
assert(edge_type == EdgeState::Possible || edge_type == new_edge_type);
assert(new_edge_type == EdgeState::Active || new_edge_type == EdgeState::Inactive);
edge_type = new_edge_type;
};
auto set_cell_state = [&cell_state, front_cell](const VD::cell_type *cell, CellState new_cell_type) -> bool {
CellState &cell_type = cell_state[cell - front_cell];
assert(cell_type == CellState::Active || cell_type == CellState::Inactive || cell_type == CellState::Boundary || cell_type == CellState::Unknown);
assert(new_cell_type == CellState::Active || new_cell_type == CellState::Inactive || new_cell_type == CellState::Boundary);
switch (cell_type) {
case CellState::Unknown:
break;
case CellState::Active:
if (new_cell_type == CellState::Inactive)
new_cell_type = CellState::Boundary;
break;
case CellState::Inactive:
if (new_cell_type == CellState::Active)
new_cell_type = CellState::Boundary;
break;
case CellState::Boundary:
return false;
}
if (cell_type != new_cell_type) {
cell_type = new_cell_type;
return true;
}
return false;
};
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 has at least one valid vertex.
assert(edge.vertex0() != nullptr);
set_edge_state_initial(&edge, outside ? EdgeState::Possible : EdgeState::Inactive);
// Opposite edge of an infinite edge is certainly not active.
set_edge_state_initial(edge.twin(), EdgeState::Inactive);
if (edge.is_secondary()) {
// edge.vertex0() must lie on source contour.
const VD::cell_type *cell = edge.cell();
const VD::cell_type *cell2 = edge.twin()->cell();
if (cell->contains_segment())
std::swap(cell, cell2);
// State of a cell containing a boundary point is known.
assert(cell->contains_point());
set_cell_state(cell, outside ? CellState::Active : CellState::Inactive);
// State of a cell containing a boundary edge is Boundary.
assert(cell2->contains_segment());
set_cell_state(cell2, CellState::Boundary);
}
} else if (edge.vertex0() != nullptr) {
// Finite edge.
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;
}
if (line) {
const VD::vertex_type *v1 = edge.vertex1();
const VD::cell_type *cell2 = (cell == edge.cell()) ? edge.twin()->cell() : edge.cell();
assert(v1);
const Point *pt_on_contour = nullptr;
if (cell == edge.cell() && edge.twin()->cell()->contains_segment()) {
// Constrained bisector of two segments.
// If the two segments share a point, then one end of the current Voronoi edge shares this point as well.
// 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;
} else if (edge.is_secondary()) {
assert(edge.is_linear());
// One end of the current Voronoi edge shares a point of a contour.
assert(edge.cell()->contains_point() != edge.twin()->cell()->contains_point());
const Line &line2 = lines[cell2->source_index()];
pt_on_contour = &((cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line2.a : line2.b);
}
if (pt_on_contour) {
// One end of the current Voronoi edge shares a point of a contour.
// Find out which one it is.
const VD::vertex_type *v0 = edge.vertex0();
Vec2d vec0(v0->x() - pt_on_contour->x(), v0->y() - pt_on_contour->y());
Vec2d vec1(v1->x() - pt_on_contour->x(), v1->y() - pt_on_contour->y());
double d0 = vec0.squaredNorm();
double d1 = vec1.squaredNorm();
assert(std::min(d0, d1) < SCALED_EPSILON * SCALED_EPSILON);
if (d0 < d1) {
// v0 is equal to pt.
} else {
// Skip secondary edge pointing to a contour point.
set_edge_state_initial(&edge, EdgeState::Inactive);
continue;
}
}
Vec2d l0(line->a.cast<double>());
Vec2d lv((line->b - line->a).cast<double>());
double side = cross2(lv, Vec2d(v1->x(), v1->y()) - l0);
bool edge_active = outside ? (side < 0.) : (side > 0.);
set_edge_state_initial(&edge, edge_active ? EdgeState::Possible : EdgeState::Inactive);
assert(cell->contains_segment());
set_cell_state(cell,
pt_on_contour ? CellState::Boundary :
edge_active ? CellState::Active : CellState::Inactive);
set_cell_state(cell2,
(pt_on_contour && cell2->contains_segment()) ?
CellState::Boundary :
edge_active ? CellState::Active : CellState::Inactive);
}
}
{
// Perform one round of expansion marking Voronoi edges and cells next to boundary cells as active / inactive.
std::vector<const VD::cell_type*> cell_queue;
for (const VD::edge_type &edge : vd.edges())
if (edge_state[&edge - front_edge] == EdgeState::Unknown) {
assert(edge.cell()->contains_point() && edge.twin()->cell()->contains_point());
// Edge separating two point sources, not yet classified as inside / outside.
CellState cs = cell_state[edge.cell() - front_cell];
CellState cs2 = cell_state[edge.twin()->cell() - front_cell];
if (cs != CellState::Unknown || cs2 != CellState::Unknown) {
if (cs == CellState::Unknown) {
cs = cs2;
if (set_cell_state(edge.cell(), cs))
cell_queue.emplace_back(edge.cell());
} else if (set_cell_state(edge.twin()->cell(), cs))
cell_queue.emplace_back(edge.twin()->cell());
EdgeState es = (cs == CellState::Active) ? EdgeState::Possible : EdgeState::Inactive;
set_edge_state_initial(&edge, es);
set_edge_state_initial(edge.twin(), es);
} else {
const VD::edge_type *e = edge.twin()->rot_prev();
do {
EdgeState es = edge_state[e->twin() - front_edge];
if (es != EdgeState::Unknown) {
assert(es == EdgeState::Possible || es == EdgeState::Inactive);
set_edge_state_initial(&edge, es);
CellState cs = (es == EdgeState::Possible) ? CellState::Active : CellState::Inactive;
if (set_cell_state(edge.cell(), cs))
cell_queue.emplace_back(edge.cell());
if (set_cell_state(edge.twin()->cell(), cs))
cell_queue.emplace_back(edge.twin()->cell());
break;
}
e = e->rot_prev();
} while (e != edge.twin());
}
}
// 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 CellState cs = cell_state[cell - front_cell];
cell_queue.pop_back();
const VD::edge_type *first_edge = cell->incident_edge();
const VD::edge_type *edge = cell->incident_edge();
EdgeState es = (cs == CellState::Active) ? EdgeState::Possible : EdgeState::Inactive;
do {
if (set_cell_state(edge->twin()->cell(), cs)) {
set_edge_state_initial(edge, es);
set_edge_state_initial(edge->twin(), es);
cell_queue.emplace_back(edge->twin()->cell());
}
edge = edge->next();
} while (edge != first_edge);
}
}
if (! outside)
offset_distance = - offset_distance;
#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_state[&edge - front_edge] == EdgeState::Possible) {
const VD::vertex_type *v0 = edge.vertex0();
const VD::vertex_type *v1 = edge.vertex1();
assert(v0 != nullptr);
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_offset_point(vd.num_edges(), Vec2d());
const double offset_distance2 = offset_distance * offset_distance;
for (const VD::edge_type &edge : vd.edges()) {
assert(edge_state[&edge - front_edge] != EdgeState::Unknown);
size_t edge_idx = &edge - front_edge;
if (edge_state[edge_idx] == EdgeState::Possible) {
// Edge candidate, intersection points were not calculated yet.
const VD::vertex_type *v0 = edge.vertex0();
const VD::vertex_type *v1 = edge.vertex1();
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() - front_edge;
if (v1 == nullptr) {
assert(edge.is_infinite());
assert(edge.is_linear());
assert(edge_state[edge_idx2] == EdgeState::Inactive);
if (cell->contains_point() && cell2->contains_point()) {
assert(! edge.is_secondary());
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;
double dmin2 = (Vec2d(v0->x(), v0->y()) - pt0.cast<double>()).squaredNorm();
assert(dmin2 >= SCALED_EPSILON * SCALED_EPSILON);
if (dmin2 <= offset_distance2) {
// There shall be an intersection of this unconstrained edge with the offset curve.
// Direction vector of this unconstrained Voronoi edge.
Vec2d dir(double(pt0.y() - pt1.y()), double(pt1.x() - pt0.x()));
Vec2d pt(v0->x(), v0->y());
double t = detail::first_circle_segment_intersection_parameter(Vec2d(pt0.x(), pt0.y()), offset_distance, pt, dir);
edge_offset_point[edge_idx] = pt + t * dir;
set_edge_state_final(edge_idx, EdgeState::Active);
} else
set_edge_state_final(edge_idx, EdgeState::Inactive);
} 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);
#ifndef NDEBUG
if (cell->contains_segment()) {
const Point &pt1 = (cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b;
assert((pt1.x() == line0.a.x() && pt1.y() == line0.a.y()) ||
(pt1.x() == line0.b.x() && pt1.y() == line0.b.y()));
} else {
const Point &pt0 = (cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b;
assert((pt0.x() == line1.a.x() && pt0.y() == line1.a.y()) ||
(pt0.x() == line1.b.x() && pt0.y() == line1.b.y()));
}
assert((Vec2d(v0->x(), v0->y()) - 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);
edge_offset_point[edge_idx] = ipt.cast<double>() + offset_distance * Vec2d(line.b.y() - line.a.y(), line.a.x() - line.b.x()).normalized();
set_edge_state_final(edge_idx, EdgeState::Active);
}
// The other edge of an unconstrained edge starting with null vertex shall never be intersected.
set_edge_state_final(edge_idx2, EdgeState::Inactive);
} else if (edge.is_secondary()) {
assert(edge.is_linear());
assert(cell->contains_point() != cell2->contains_point());
const Line &line0 = lines[edge.cell()->source_index()];
const Line &line1 = lines[edge.twin()->cell()->source_index()];
const Point &pt = cell->contains_point() ?
((cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b) :
((cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b);
const Line &line = cell->contains_segment() ? line0 : line1;
assert(pt == line.a || pt == line.b);
assert((pt.cast<double>() - Vec2d(v0->x(), v0->y())).norm() < SCALED_EPSILON);
Vec2d dir(v1->x() - v0->x(), v1->y() - v0->y());
double l2 = dir.squaredNorm();
if (offset_distance2 <= l2) {
edge_offset_point[edge_idx] = pt.cast<double>() + (offset_distance / sqrt(l2)) * dir;
set_edge_state_final(edge_idx, EdgeState::Active);
} else {
set_edge_state_final(edge_idx, EdgeState::Inactive);
}
set_edge_state_final(edge_idx2, EdgeState::Inactive);
} else {
// Finite edge has valid points at both sides.
bool done = false;
if (cell->contains_segment() && cell2->contains_segment()) {
// This edge is a bisector of two line segments. Project v0, v1 onto one of the line segments.
Vec2d pt(line0.a.cast<double>());
Vec2d dir(line0.b.cast<double>() - pt);
Vec2d vec0 = Vec2d(v0->x(), v0->y()) - pt;
Vec2d vec1 = Vec2d(v1->x(), v1->y()) - pt;
double l2 = dir.squaredNorm();
assert(l2 > 0.);
double dmin = (dir * (vec0.dot(dir) / l2) - vec0).squaredNorm();
double dmax = (dir * (vec1.dot(dir) / l2) - vec1).squaredNorm();
bool flip = dmin > dmax;
if (flip)
std::swap(dmin, dmax);
if (offset_distance2 >= dmin && offset_distance2 <= dmax) {
// Intersect. Maximum one intersection will be found.
// This edge is a bisector of two line segments. Distance to the input polygon increases/decreases monotonically.
dmin = sqrt(dmin);
dmax = sqrt(dmax);
assert(offset_distance > dmin - EPSILON && offset_distance < dmax + EPSILON);
double ddif = dmax - dmin;
if (ddif == 0.) {
// line, line2 are exactly parallel. This is a singular case, the offset curve should miss it.
} else {
if (flip) {
std::swap(edge_idx, edge_idx2);
std::swap(v0, v1);
}
double t = clamp(0., 1., (offset_distance - dmin) / ddif);
edge_offset_point[edge_idx] = Vec2d(lerp(v0->x(), v1->x(), t), lerp(v0->y(), v1->y(), t));
set_edge_state_final(edge_idx, EdgeState::Active);
set_edge_state_final(edge_idx2, EdgeState::Inactive);
done = true;
}
}
} else {
assert(cell->contains_point() || cell2->contains_point());
bool point_vs_segment = cell->contains_point() != cell2->contains_point();
const Point &pt0 = cell->contains_point() ?
((cell->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line0.a : line0.b) :
((cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b);
// 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());
double d0 = (p0 - px).squaredNorm();
double d1 = (p1 - px).squaredNorm();
double dmin = std::min(d0, d1);
double dmax = std::max(d0, d1);
bool has_intersection = false;
bool possibly_two_points = false;
if (offset_distance2 <= dmax) {
if (offset_distance2 >= dmin) {
has_intersection = true;
} else {
double dmin_new = dmin;
if (point_vs_segment) {
// Project on the source segment.
const Line &line = cell->contains_segment() ? line0 : line1;
const Vec2d pt_line = line.a.cast<double>();
const Vec2d v_line = (line.b - line.a).cast<double>();
double t0 = (p0 - pt_line).dot(v_line);
double t1 = (p1 - pt_line).dot(v_line);
double tx = (px - pt_line).dot(v_line);
if ((tx >= t0 && tx <= t1) || (tx >= t1 && tx <= t0)) {
// Projection of the Point site falls between the projections of the Voronoi edge end points
// onto the Line site.
Vec2d ft = pt_line + (tx / v_line.squaredNorm()) * v_line;
dmin_new = (ft - px).squaredNorm() * 0.25;
}
} else {
// Point-Point Voronoi sites. Project point site onto the current Voronoi edge.
Vec2d v = p1 - p0;
auto l2 = v.squaredNorm();
assert(l2 > 0);
auto t = v.dot(px - p0);
if (t >= 0. && t <= l2) {
// Projection falls onto the Voronoi edge. Calculate foot point and distance.
Vec2d ft = p0 + (t / l2) * v;
dmin_new = (ft - px).squaredNorm();
}
}
assert(dmin_new < dmax + SCALED_EPSILON);
assert(dmin_new < dmin + SCALED_EPSILON);
if (dmin_new < dmin) {
dmin = dmin_new;
has_intersection = possibly_two_points = offset_distance2 >= dmin;
}
}
}
if (has_intersection) {
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 {
const Point &pt1 = (cell2->source_category() == boost::polygon::SOURCE_CATEGORY_SEGMENT_START_POINT) ? line1.a : line1.b;
intersections = detail::point_point_equal_distance_points(pt0, pt1, offset_distance);
}
// If the span of distances of start / end point / foot point to the point site indicate an intersection,
// we should find one.
assert(intersections.count > 0);
if (intersections.count == 2) {
// Now decide which points fall on this Voronoi edge.
// Tangential points (single intersection) are ignored.
if (possibly_two_points) {
Vec2d v = p1 - p0;
double l2 = v.squaredNorm();
double t0 = v.dot(intersections.pts[0] - p0);
double t1 = v.dot(intersections.pts[1] - p0);
if (t0 > t1) {
std::swap(t0, t1);
std::swap(intersections.pts[0], intersections.pts[1]);
}
// Remove points outside of the line range.
if (t0 < 0. || t0 > l2) {
if (t1 < 0. || t1 > l2)
intersections.count = 0;
else {
-- intersections.count;
t0 = t1;
intersections.pts[0] = intersections.pts[1];
}
} else if (t1 < 0. || t1 > l2)
-- intersections.count;
} else {
// Take the point furthest from the end points of the Voronoi edge or a Voronoi parabolic arc.
double d0 = std::max((intersections.pts[0] - p0).squaredNorm(), (intersections.pts[0] - p1).squaredNorm());
double d1 = std::max((intersections.pts[1] - p0).squaredNorm(), (intersections.pts[1] - p1).squaredNorm());
if (d0 > d1)
intersections.pts[0] = intersections.pts[1];
-- intersections.count;
}
assert(intersections.count > 0);
if (intersections.count == 2) {
set_edge_state_final(edge_idx, EdgeState::Active);
set_edge_state_final(edge_idx2, EdgeState::Active);
edge_offset_point[edge_idx] = intersections.pts[1];
edge_offset_point[edge_idx2] = intersections.pts[0];
done = true;
} else if (intersections.count == 1) {
if (d1 < d0)
std::swap(edge_idx, edge_idx2);
set_edge_state_final(edge_idx, EdgeState::Active);
set_edge_state_final(edge_idx2, EdgeState::Inactive);
edge_offset_point[edge_idx] = intersections.pts[0];
done = true;
}
}
}
}
if (! done) {
set_edge_state_final(edge_idx, EdgeState::Inactive);
set_edge_state_final(edge_idx2, EdgeState::Inactive);
}
}
}
}
#ifndef NDEBUG
for (const VD::edge_type &edge : vd.edges()) {
assert(edge_state[&edge - front_edge] == EdgeState::Inactive || edge_state[&edge - front_edge] == EdgeState::Active);
// None of a new edge candidate may start with null vertex.
assert(edge_state[&edge - front_edge] == EdgeState::Inactive || edge.vertex0() != nullptr);
assert(edge_state[edge.twin() - front_edge] == EdgeState::Inactive || edge.twin()->vertex0() != nullptr);
}
#endif // NDEBUG
#ifdef VORONOI_DEBUG_OUT
{
Lines helper_lines;
for (const VD::edge_type &edge : vd.edges())
if (edge_state[&edge - front_edge] == EdgeState::Active)
helper_lines.emplace_back(Line(Point(edge.vertex0()->x(), edge.vertex0()->y()), Point(edge_offset_point[&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_state, 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_state[edge->twin() - front_edge] == EdgeState::Active)
return edge->twin();
// assert(false);
return nullptr;
};
#ifndef NDEBUG
auto dist_to_site = [&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();
};
#endif /* NDEBUG */
// Track the offset curves.
Polygons out;
double angle_step = 2. * acos((offset_distance - discretization_error) / offset_distance);
double cos_threshold = cos(angle_step);
for (size_t seed_edge_idx = 0; seed_edge_idx < vd.num_edges(); ++ seed_edge_idx)
if (edge_state[seed_edge_idx] == EdgeState::Active) {
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 helper_lines;
dump_voronoi_to_svg(debug_out_path("voronoi-offset-open-loop-%d.svg", irun).c_str(), vd, Points(), lines, Polygons(), to_lines(poly));
}
#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();
edge_state[next_edge - front_edge] = EdgeState::Inactive;
Vec2d p1 = edge_offset_point[edge - front_edge];
Vec2d p2 = edge_offset_point[next_edge - front_edge];
#ifndef NDEBUG
{
double err = dist_to_site(*cell, p1) - offset_distance;
double err2 = dist_to_site(*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 &center = ((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;
} while (edge != start_edge);
out.emplace_back(std::move(poly));
}
return out;
}
} // namespace Slic3r
Reference in a new issue View git blame Copy permalink