juli/OrcaSlicer
1
1
Fork 0
mirror of https://github.com/OrcaSlicer/OrcaSlicer.git synced 2026-05-31 00:52:51 +00:00
Code Issues Packages Projects Releases Wiki Activity
Files
4156b51c18ffc533994b083a4857cf813724a12f
OrcaSlicer/xs/src/libslic3r/Geometry.cpp

1004 lines
38 KiB
C++
Raw Blame History

#include "Geometry.hpp"
#include "ClipperUtils.hpp"
#include "ExPolygon.hpp"
#include "Line.hpp"
#include "PolylineCollection.hpp"
#include "clipper.hpp"
#include <algorithm>
#include <cassert>
#include <cmath>
#include <list>
#include <map>
#include <set>
#include <utility>
#include <stack>
#include <vector>
#ifdef SLIC3R_DEBUG
#include "SVG.hpp"
#endif
#ifdef SLIC3R_DEBUG
namespace boost { namespace polygon {
// The following code for the visualization of the boost Voronoi diagram is based on:
//
// Boost.Polygon library voronoi_graphic_utils.hpp header file
// Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
template <typename CT>
class voronoi_visual_utils {
public:
// Discretize parabolic Voronoi edge.
// Parabolic Voronoi edges are always formed by one point and one segment
// from the initial input set.
//
// Args:
// point: input point.
// segment: input segment.
// max_dist: maximum discretization distance.
// discretization: point discretization of the given Voronoi edge.
//
// Template arguments:
// InCT: coordinate type of the input geometries (usually integer).
// Point: point type, should model point concept.
// Segment: segment type, should model segment concept.
//
// Important:
// discretization should contain both edge endpoints initially.
template <class InCT1, class InCT2,
template<class> class Point,
template<class> class Segment>
static
typename enable_if<
typename gtl_and<
typename gtl_if<
typename is_point_concept<
typename geometry_concept< Point<InCT1> >::type
>::type
>::type,
typename gtl_if<
typename is_segment_concept<
typename geometry_concept< Segment<InCT2> >::type
>::type
>::type
>::type,
void
>::type discretize(
const Point<InCT1>& point,
const Segment<InCT2>& segment,
const CT max_dist,
std::vector< Point<CT> >* discretization) {
// Apply the linear transformation to move start point of the segment to
// the point with coordinates (0, 0) and the direction of the segment to
// coincide the positive direction of the x-axis.
CT segm_vec_x = cast(x(high(segment))) - cast(x(low(segment)));
CT segm_vec_y = cast(y(high(segment))) - cast(y(low(segment)));
CT sqr_segment_length = segm_vec_x * segm_vec_x + segm_vec_y * segm_vec_y;
// Compute x-coordinates of the endpoints of the edge
// in the transformed space.
CT projection_start = sqr_segment_length *
get_point_projection((*discretization)[0], segment);
CT projection_end = sqr_segment_length *
get_point_projection((*discretization)[1], segment);
// Compute parabola parameters in the transformed space.
// Parabola has next representation:
// f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y).
CT point_vec_x = cast(x(point)) - cast(x(low(segment)));
CT point_vec_y = cast(y(point)) - cast(y(low(segment)));
CT rot_x = segm_vec_x * point_vec_x + segm_vec_y * point_vec_y;
CT rot_y = segm_vec_x * point_vec_y - segm_vec_y * point_vec_x;
// Save the last point.
Point<CT> last_point = (*discretization)[1];
discretization->pop_back();
// Use stack to avoid recursion.
std::stack<CT> point_stack;
point_stack.push(projection_end);
CT cur_x = projection_start;
CT cur_y = parabola_y(cur_x, rot_x, rot_y);
// Adjust max_dist parameter in the transformed space.
const CT max_dist_transformed = max_dist * max_dist * sqr_segment_length;
while (!point_stack.empty()) {
CT new_x = point_stack.top();
CT new_y = parabola_y(new_x, rot_x, rot_y);
// Compute coordinates of the point of the parabola that is
// furthest from the current line segment.
CT mid_x = (new_y - cur_y) / (new_x - cur_x) * rot_y + rot_x;
CT mid_y = parabola_y(mid_x, rot_x, rot_y);
// Compute maximum distance between the given parabolic arc
// and line segment that discretize it.
CT dist = (new_y - cur_y) * (mid_x - cur_x) -
(new_x - cur_x) * (mid_y - cur_y);
dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) +
(new_x - cur_x) * (new_x - cur_x));
if (dist <= max_dist_transformed) {
// Distance between parabola and line segment is less than max_dist.
point_stack.pop();
CT inter_x = (segm_vec_x * new_x - segm_vec_y * new_y) /
sqr_segment_length + cast(x(low(segment)));
CT inter_y = (segm_vec_x * new_y + segm_vec_y * new_x) /
sqr_segment_length + cast(y(low(segment)));
discretization->push_back(Point<CT>(inter_x, inter_y));
cur_x = new_x;
cur_y = new_y;
} else {
point_stack.push(mid_x);
}
}
// Update last point.
discretization->back() = last_point;
}
private:
// Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b).
static CT parabola_y(CT x, CT a, CT b) {
return ((x - a) * (x - a) + b * b) / (b + b);
}
// Get normalized length of the distance between:
// 1) point projection onto the segment
// 2) start point of the segment
// Return this length divided by the segment length. This is made to avoid
// sqrt computation during transformation from the initial space to the
// transformed one and vice versa. The assumption is made that projection of
// the point lies between the start-point and endpoint of the segment.
template <class InCT,
template<class> class Point,
template<class> class Segment>
static
typename enable_if<
typename gtl_and<
typename gtl_if<
typename is_point_concept<
typename geometry_concept< Point<int> >::type
>::type
>::type,
typename gtl_if<
typename is_segment_concept<
typename geometry_concept< Segment<long> >::type
>::type
>::type
>::type,
CT
>::type get_point_projection(
const Point<CT>& point, const Segment<InCT>& segment) {
CT segment_vec_x = cast(x(high(segment))) - cast(x(low(segment)));
CT segment_vec_y = cast(y(high(segment))) - cast(y(low(segment)));
CT point_vec_x = x(point) - cast(x(low(segment)));
CT point_vec_y = y(point) - cast(y(low(segment)));
CT sqr_segment_length =
segment_vec_x * segment_vec_x + segment_vec_y * segment_vec_y;
CT vec_dot = segment_vec_x * point_vec_x + segment_vec_y * point_vec_y;
return vec_dot / sqr_segment_length;
}
template <typename InCT>
static CT cast(const InCT& value) {
return static_cast<CT>(value);
}
};
} } // namespace boost::polygon
#endif
using namespace boost::polygon; // provides also high() and low()
namespace Slic3r { namespace Geometry {
static bool
sort_points (Point a, Point b)
{
return (a.x < b.x) || (a.x == b.x && a.y < b.y);
}
/* This implementation is based on Andrew's monotone chain 2D convex hull algorithm */
Polygon
convex_hull(Points points)
{
assert(points.size() >= 3);
// sort input points
std::sort(points.begin(), points.end(), sort_points);
int n = points.size(), k = 0;
Polygon hull;
hull.points.resize(2*n);
// Build lower hull
for (int i = 0; i < n; i++) {
while (k >= 2 && points[i].ccw(hull.points[k-2], hull.points[k-1]) <= 0) k--;
hull.points[k++] = points[i];
}
// Build upper hull
for (int i = n-2, t = k+1; i >= 0; i--) {
while (k >= t && points[i].ccw(hull.points[k-2], hull.points[k-1]) <= 0) k--;
hull.points[k++] = points[i];
}
hull.points.resize(k);
assert( hull.points.front().coincides_with(hull.points.back()) );
hull.points.pop_back();
return hull;
}
Polygon
convex_hull(const Polygons &polygons)
{
Points pp;
for (Polygons::const_iterator p = polygons.begin(); p != polygons.end(); ++p) {
pp.insert(pp.end(), p->points.begin(), p->points.end());
}
return convex_hull(pp);
}
/* accepts an arrayref of points and returns a list of indices
according to a nearest-neighbor walk */
void
chained_path(const Points &points, std::vector<Points::size_type> &retval, Point start_near)
{
PointConstPtrs my_points;
std::map<const Point*,Points::size_type> indices;
my_points.reserve(points.size());
for (Points::const_iterator it = points.begin(); it != points.end(); ++it) {
my_points.push_back(&*it);
indices[&*it] = it - points.begin();
}
retval.reserve(points.size());
while (!my_points.empty()) {
Points::size_type idx = start_near.nearest_point_index(my_points);
start_near = *my_points[idx];
retval.push_back(indices[ my_points[idx] ]);
my_points.erase(my_points.begin() + idx);
}
}
void
chained_path(const Points &points, std::vector<Points::size_type> &retval)
{
if (points.empty()) return; // can't call front() on empty vector
chained_path(points, retval, points.front());
}
/* retval and items must be different containers */
template<class T>
void
chained_path_items(Points &points, T &items, T &retval)
{
std::vector<Points::size_type> indices;
chained_path(points, indices);
for (std::vector<Points::size_type>::const_iterator it = indices.begin(); it != indices.end(); ++it)
retval.push_back(items[*it]);
}
template void chained_path_items(Points &points, ClipperLib::PolyNodes &items, ClipperLib::PolyNodes &retval);
bool
directions_parallel(double angle1, double angle2, double max_diff)
{
double diff = fabs(angle1 - angle2);
max_diff += EPSILON;
return diff < max_diff || fabs(diff - PI) < max_diff;
}
template<class T>
bool
contains(const std::vector<T> &vector, const Point &point)
{
for (typename std::vector<T>::const_iterator it = vector.begin(); it != vector.end(); ++it) {
if (it->contains(point)) return true;
}
return false;
}
template bool contains(const ExPolygons &vector, const Point &point);
double
rad2deg(double angle)
{
return angle / PI * 180.0;
}
double
rad2deg_dir(double angle)
{
angle = (angle < PI) ? (-angle + PI/2.0) : (angle + PI/2.0);
if (angle < 0) angle += PI;
return rad2deg(angle);
}
double
deg2rad(double angle)
{
return PI * angle / 180.0;
}
void
simplify_polygons(const Polygons &polygons, double tolerance, Polygons* retval)
{
Polygons pp;
for (Polygons::const_iterator it = polygons.begin(); it != polygons.end(); ++it) {
Polygon p = *it;
p.points.push_back(p.points.front());
p.points = MultiPoint::_douglas_peucker(p.points, tolerance);
p.points.pop_back();
pp.push_back(p);
}
Slic3r::simplify_polygons(pp, retval);
}
double
linint(double value, double oldmin, double oldmax, double newmin, double newmax)
{
return (value - oldmin) * (newmax - newmin) / (oldmax - oldmin) + newmin;
}
Pointfs
arrange(size_t total_parts, Pointf part, coordf_t dist, const BoundingBoxf* bb)
{
// use actual part size (the largest) plus separation distance (half on each side) in spacing algorithm
part.x += dist;
part.y += dist;
Pointf area;
if (bb != NULL && bb->defined) {
area = bb->size();
} else {
// bogus area size, large enough not to trigger the error below
area.x = part.x * total_parts;
area.y = part.y * total_parts;
}
// this is how many cells we have available into which to put parts
size_t cellw = floor((area.x + dist) / part.x);
size_t cellh = floor((area.y + dist) / part.y);
if (total_parts > (cellw * cellh))
CONFESS("%zu parts won't fit in your print area!\n", total_parts);
// total space used by cells
Pointf cells(cellw * part.x, cellh * part.y);
// bounding box of total space used by cells
BoundingBoxf cells_bb;
cells_bb.merge(Pointf(0,0)); // min
cells_bb.merge(cells); // max
// center bounding box to area
cells_bb.translate(
(area.x - cells.x) / 2,
(area.y - cells.y) / 2
);
// list of cells, sorted by distance from center
std::vector<ArrangeItemIndex> cellsorder;
// work out distance for all cells, sort into list
for (size_t i = 0; i <= cellw-1; ++i) {
for (size_t j = 0; j <= cellh-1; ++j) {
coordf_t cx = linint(i + 0.5, 0, cellw, cells_bb.min.x, cells_bb.max.x);
coordf_t cy = linint(j + 0.5, 0, cellh, cells_bb.min.y, cells_bb.max.y);
coordf_t xd = fabs((area.x / 2) - cx);
coordf_t yd = fabs((area.y / 2) - cy);
ArrangeItem c;
c.pos.x = cx;
c.pos.y = cy;
c.index_x = i;
c.index_y = j;
c.dist = xd * xd + yd * yd - fabs((cellw / 2) - (i + 0.5));
// binary insertion sort
{
coordf_t index = c.dist;
size_t low = 0;
size_t high = cellsorder.size();
while (low < high) {
size_t mid = (low + ((high - low) / 2)) | 0;
coordf_t midval = cellsorder[mid].index;
if (midval < index) {
low = mid + 1;
} else if (midval > index) {
high = mid;
} else {
cellsorder.insert(cellsorder.begin() + mid, ArrangeItemIndex(index, c));
goto ENDSORT;
}
}
cellsorder.insert(cellsorder.begin() + low, ArrangeItemIndex(index, c));
}
ENDSORT: true;
}
}
// the extents of cells actually used by objects
coordf_t lx = 0;
coordf_t ty = 0;
coordf_t rx = 0;
coordf_t by = 0;
// now find cells actually used by objects, map out the extents so we can position correctly
for (size_t i = 1; i <= total_parts; ++i) {
ArrangeItemIndex c = cellsorder[i - 1];
coordf_t cx = c.item.index_x;
coordf_t cy = c.item.index_y;
if (i == 1) {
lx = rx = cx;
ty = by = cy;
} else {
if (cx > rx) rx = cx;
if (cx < lx) lx = cx;
if (cy > by) by = cy;
if (cy < ty) ty = cy;
}
}
// now we actually place objects into cells, positioned such that the left and bottom borders are at 0
Pointfs positions;
for (size_t i = 1; i <= total_parts; ++i) {
ArrangeItemIndex c = cellsorder.front();
cellsorder.erase(cellsorder.begin());
coordf_t cx = c.item.index_x - lx;
coordf_t cy = c.item.index_y - ty;
positions.push_back(Pointf(cx * part.x, cy * part.y));
}
if (bb != NULL && bb->defined) {
for (Pointfs::iterator p = positions.begin(); p != positions.end(); ++p) {
p->x += bb->min.x;
p->y += bb->min.y;
}
}
return positions;
}
#ifdef SLIC3R_DEBUG
// The following code for the visualization of the boost Voronoi diagram is based on:
//
// Boost.Polygon library voronoi_visualizer.cpp file
// Copyright Andrii Sydorchuk 2010-2012.
// Distributed under the Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
namespace Voronoi { namespace Internal {
typedef double coordinate_type;
typedef boost::polygon::point_data<coordinate_type> point_type;
typedef boost::polygon::segment_data<coordinate_type> segment_type;
typedef boost::polygon::rectangle_data<coordinate_type> rect_type;
// typedef voronoi_builder<int> VB;
typedef boost::polygon::voronoi_diagram<coordinate_type> VD;
typedef VD::cell_type cell_type;
typedef VD::cell_type::source_index_type source_index_type;
typedef VD::cell_type::source_category_type source_category_type;
typedef VD::edge_type edge_type;
typedef VD::cell_container_type cell_container_type;
typedef VD::cell_container_type vertex_container_type;
typedef VD::edge_container_type edge_container_type;
typedef VD::const_cell_iterator const_cell_iterator;
typedef VD::const_vertex_iterator const_vertex_iterator;
typedef VD::const_edge_iterator const_edge_iterator;
static const std::size_t EXTERNAL_COLOR = 1;
inline void color_exterior(const VD::edge_type* edge)
{
if (edge->color() == EXTERNAL_COLOR)
return;
edge->color(EXTERNAL_COLOR);
edge->twin()->color(EXTERNAL_COLOR);
const VD::vertex_type* v = edge->vertex1();
if (v == NULL || !edge->is_primary())
return;
v->color(EXTERNAL_COLOR);
const VD::edge_type* e = v->incident_edge();
do {
color_exterior(e);
e = e->rot_next();
} while (e != v->incident_edge());
}
inline point_type retrieve_point(const std::vector<segment_type> &segments, const cell_type& cell)
{
assert(cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT || cell.source_category() == SOURCE_CATEGORY_SEGMENT_END_POINT);
return (cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT) ? low(segments[cell.source_index()]) : high(segments[cell.source_index()]);
}
inline void clip_infinite_edge(const std::vector<segment_type> &segments, const edge_type& edge, coordinate_type bbox_max_size, std::vector<point_type>* clipped_edge)
{
const cell_type& cell1 = *edge.cell();
const cell_type& cell2 = *edge.twin()->cell();
point_type origin, direction;
// Infinite edges could not be created by two segment sites.
if (cell1.contains_point() && cell2.contains_point()) {
point_type p1 = retrieve_point(segments, cell1);
point_type p2 = retrieve_point(segments, cell2);
origin.x((p1.x() + p2.x()) * 0.5);
origin.y((p1.y() + p2.y()) * 0.5);
direction.x(p1.y() - p2.y());
direction.y(p2.x() - p1.x());
} else {
origin = cell1.contains_segment() ? retrieve_point(segments, cell2) : retrieve_point(segments, cell1);
segment_type segment = cell1.contains_segment() ? segments[cell1.source_index()] : segments[cell2.source_index()];
coordinate_type dx = high(segment).x() - low(segment).x();
coordinate_type dy = high(segment).y() - low(segment).y();
if ((low(segment) == origin) ^ cell1.contains_point()) {
direction.x(dy);
direction.y(-dx);
} else {
direction.x(-dy);
direction.y(dx);
}
}
coordinate_type koef = bbox_max_size / (std::max)(fabs(direction.x()), fabs(direction.y()));
if (edge.vertex0() == NULL) {
clipped_edge->push_back(point_type(
origin.x() - direction.x() * koef,
origin.y() - direction.y() * koef));
} else {
clipped_edge->push_back(
point_type(edge.vertex0()->x(), edge.vertex0()->y()));
}
if (edge.vertex1() == NULL) {
clipped_edge->push_back(point_type(
origin.x() + direction.x() * koef,
origin.y() + direction.y() * koef));
} else {
clipped_edge->push_back(
point_type(edge.vertex1()->x(), edge.vertex1()->y()));
}
}
inline void sample_curved_edge(const std::vector<segment_type> &segments, const edge_type& edge, std::vector<point_type> &sampled_edge, coordinate_type max_dist)
{
point_type point = edge.cell()->contains_point() ?
retrieve_point(segments, *edge.cell()) :
retrieve_point(segments, *edge.twin()->cell());
segment_type segment = edge.cell()->contains_point() ?
segments[edge.twin()->cell()->source_index()] :
segments[edge.cell()->source_index()];
::boost::polygon::voronoi_visual_utils<coordinate_type>::discretize(point, segment, max_dist, &sampled_edge);
}
} /* namespace Internal */ } // namespace Voronoi
static inline void dump_voronoi_to_svg(const Lines &lines, /* const */ voronoi_diagram<double> &vd, const ThickPolylines *polylines, const char *path)
{
const double scale = 0.2;
const std::string inputSegmentPointColor = "lightseagreen";
const coord_t inputSegmentPointRadius = coord_t(0.09 * scale / SCALING_FACTOR);
const std::string inputSegmentColor = "lightseagreen";
const coord_t inputSegmentLineWidth = coord_t(0.03 * scale / SCALING_FACTOR);
const std::string voronoiPointColor = "black";
const coord_t voronoiPointRadius = coord_t(0.06 * scale / SCALING_FACTOR);
const std::string voronoiLineColorPrimary = "black";
const std::string voronoiLineColorSecondary = "green";
const std::string voronoiArcColor = "red";
const coord_t voronoiLineWidth = coord_t(0.02 * scale / SCALING_FACTOR);
const bool internalEdgesOnly = false;
const bool primaryEdgesOnly = false;
BoundingBox bbox = BoundingBox(lines);
bbox.min.x -= coord_t(1. / SCALING_FACTOR);
bbox.min.y -= coord_t(1. / SCALING_FACTOR);
bbox.max.x += coord_t(1. / SCALING_FACTOR);
bbox.max.y += coord_t(1. / SCALING_FACTOR);
::Slic3r::SVG svg(path, bbox);
if (polylines != NULL)
svg.draw(*polylines, "lime", "lime", voronoiLineWidth);
// bbox.scale(1.2);
// For clipping of half-lines to some reasonable value.
// The line will then be clipped by the SVG viewer anyway.
const double bbox_dim_max = double(bbox.max.x - bbox.min.x) + double(bbox.max.y - bbox.min.y);
// For the discretization of the Voronoi parabolic segments.
const double discretization_step = 0.0005 * bbox_dim_max;
// Make a copy of the input segments with the double type.
std::vector<Voronoi::Internal::segment_type> segments;
for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++ it)
segments.push_back(Voronoi::Internal::segment_type(
Voronoi::Internal::point_type(double(it->a.x), double(it->a.y)),
Voronoi::Internal::point_type(double(it->b.x), double(it->b.y))));
// Color exterior edges.
for (voronoi_diagram<double>::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it)
if (!it->is_finite())
Voronoi::Internal::color_exterior(&(*it));
// Draw the end points of the input polygon.
for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it) {
svg.draw(it->a, inputSegmentPointColor, inputSegmentPointRadius);
svg.draw(it->b, inputSegmentPointColor, inputSegmentPointRadius);
}
// Draw the input polygon.
for (Lines::const_iterator it = lines.begin(); it != lines.end(); ++it)
svg.draw(Line(Point(coord_t(it->a.x), coord_t(it->a.y)), Point(coord_t(it->b.x), coord_t(it->b.y))), inputSegmentColor, inputSegmentLineWidth);
#if 1
// Draw voronoi vertices.
for (voronoi_diagram<double>::const_vertex_iterator it = vd.vertices().begin(); it != vd.vertices().end(); ++it)
if (! internalEdgesOnly || it->color() != Voronoi::Internal::EXTERNAL_COLOR)
svg.draw(Point(coord_t(it->x()), coord_t(it->y())), voronoiPointColor, voronoiPointRadius);
for (voronoi_diagram<double>::const_edge_iterator it = vd.edges().begin(); it != vd.edges().end(); ++it) {
if (primaryEdgesOnly && !it->is_primary())
continue;
if (internalEdgesOnly && (it->color() == Voronoi::Internal::EXTERNAL_COLOR))
continue;
std::vector<Voronoi::Internal::point_type> samples;
std::string color = voronoiLineColorPrimary;
if (!it->is_finite()) {
Voronoi::Internal::clip_infinite_edge(segments, *it, bbox_dim_max, &samples);
if (! it->is_primary())
color = voronoiLineColorSecondary;
} else {
// Store both points of the segment into samples. sample_curved_edge will split the initial line
// until the discretization_step is reached.
samples.push_back(Voronoi::Internal::point_type(it->vertex0()->x(), it->vertex0()->y()));
samples.push_back(Voronoi::Internal::point_type(it->vertex1()->x(), it->vertex1()->y()));
if (it->is_curved()) {
Voronoi::Internal::sample_curved_edge(segments, *it, samples, discretization_step);
color = voronoiArcColor;
} else if (! it->is_primary())
color = voronoiLineColorSecondary;
}
for (std::size_t i = 0; i + 1 < samples.size(); ++i)
svg.draw(Line(Point(coord_t(samples[i].x()), coord_t(samples[i].y())), Point(coord_t(samples[i+1].x()), coord_t(samples[i+1].y()))), color, voronoiLineWidth);
}
#endif
if (polylines != NULL)
svg.draw(*polylines, "blue", voronoiLineWidth);
svg.Close();
}
#endif /* SLIC3R_DEBUG */
// Euclidian distance of two boost::polygon points.
template<typename T>
T dist(const boost::polygon::point_data<T> &p1,const boost::polygon::point_data<T> &p2)
{
T dx = p2.x() - p1.x();
T dy = p2.y() - p1.y();
return sqrt(dx*dx+dy*dy);
}
// Find a foot point of "px" on a segment "seg".
template<typename segment_type, typename point_type>
inline point_type project_point_to_segment(segment_type &seg, point_type &px)
{
typedef typename point_type::coordinate_type T;
const point_type &p0 = low(seg);
const point_type &p1 = high(seg);
const point_type dir(p1.x()-p0.x(), p1.y()-p0.y());
const point_type dproj(px.x()-p0.x(), px.y()-p0.y());
const T t = (dir.x()*dproj.x() + dir.y()*dproj.y()) / (dir.x()*dir.x() + dir.y()*dir.y());
assert(t >= T(-1e-6) && t <= T(1. + 1e-6));
return point_type(p0.x() + t*dir.x(), p0.y() + t*dir.y());
}
template<typename VD, typename SEGMENTS>
inline const typename VD::point_type retrieve_cell_point(const typename VD::cell_type& cell, const SEGMENTS &segments)
{
assert(cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT || cell.source_category() == SOURCE_CATEGORY_SEGMENT_END_POINT);
return (cell.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT) ? low(segments[cell.source_index()]) : high(segments[cell.source_index()]);
}
template<typename VD, typename SEGMENTS>
inline std::pair<typename VD::coord_type, typename VD::coord_type>
measure_edge_thickness(const VD &vd, const typename VD::edge_type& edge, const SEGMENTS &segments)
{
typedef typename VD::coord_type T;
const typename VD::point_type pa(edge.vertex0()->x(), edge.vertex0()->y());
const typename VD::point_type pb(edge.vertex1()->x(), edge.vertex1()->y());
const typename VD::cell_type &cell1 = *edge.cell();
const typename VD::cell_type &cell2 = *edge.twin()->cell();
if (cell1.contains_segment()) {
if (cell2.contains_segment()) {
// Both cells contain a linear segment, the left / right cells are symmetric.
// Project pa, pb to the left segment.
const typename VD::segment_type segment1 = segments[cell1.source_index()];
const typename VD::point_type p1a = project_point_to_segment(segment1, pa);
const typename VD::point_type p1b = project_point_to_segment(segment1, pb);
return std::pair<T, T>(T(2.)*dist(pa, p1a), T(2.)*dist(pb, p1b));
} else {
// 1st cell contains a linear segment, 2nd cell contains a point.
// The medial axis between the cells is a parabolic arc.
// Project pa, pb to the left segment.
const typename VD::point_type p2 = retrieve_cell_point<VD>(cell2, segments);
return std::pair<T, T>(T(2.)*dist(pa, p2), T(2.)*dist(pb, p2));
}
} else if (cell2.contains_segment()) {
// 1st cell contains a point, 2nd cell contains a linear segment.
// The medial axis between the cells is a parabolic arc.
const typename VD::point_type p1 = retrieve_cell_point<VD>(cell1, segments);
return std::pair<T, T>(T(2.)*dist(pa, p1), T(2.)*dist(pb, p1));
} else {
// Both cells contain a point. The left / right regions are triangular and symmetric.
const typename VD::point_type p1 = retrieve_cell_point<VD>(cell1, segments);
return std::pair<T, T>(T(2.)*dist(pa, p1), T(2.)*dist(pb, p1));
}
}
// Converts the Line instances of Lines vector to VD::segment_type.
template<typename VD>
class Lines2VDSegments
{
public:
Lines2VDSegments(const Lines &alines) : lines(alines) {}
typename VD::segment_type operator[](size_t idx) const {
return typename VD::segment_type(
typename VD::point_type(typename VD::coord_type(lines[idx].a.x), typename VD::coord_type(lines[idx].a.y)),
typename VD::point_type(typename VD::coord_type(lines[idx].b.x), typename VD::coord_type(lines[idx].b.y)));
}
private:
const Lines &lines;
};
void
MedialAxis::build(ThickPolylines* polylines)
{
construct_voronoi(this->lines.begin(), this->lines.end(), &this->vd);
/*
// DEBUG: dump all Voronoi edges
{
for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) {
if (edge->is_infinite()) continue;
ThickPolyline polyline;
polyline.points.push_back(Point( edge->vertex0()->x(), edge->vertex0()->y() ));
polyline.points.push_back(Point( edge->vertex1()->x(), edge->vertex1()->y() ));
polylines->push_back(polyline);
}
return;
}
*/
typedef const VD::vertex_type vert_t;
typedef const VD::edge_type edge_t;
// collect valid edges (i.e. prune those not belonging to MAT)
// note: this keeps twins, so it inserts twice the number of the valid edges
this->valid_edges.clear();
{
std::set<const VD::edge_type*> seen_edges;
for (VD::const_edge_iterator edge = this->vd.edges().begin(); edge != this->vd.edges().end(); ++edge) {
// if we only process segments representing closed loops, none if the
// infinite edges (if any) would be part of our MAT anyway
if (edge->is_secondary() || edge->is_infinite()) continue;
// don't re-validate twins
if (seen_edges.find(&*edge) != seen_edges.end()) continue;
seen_edges.insert(&*edge);
seen_edges.insert(edge->twin());
if (!this->validate_edge(&*edge)) continue;
this->valid_edges.insert(&*edge);
this->valid_edges.insert(edge->twin());
}
}
this->edges = this->valid_edges;
// iterate through the valid edges to build polylines
while (!this->edges.empty()) {
const edge_t* edge = *this->edges.begin();
// start a polyline
ThickPolyline polyline;
polyline.points.push_back(Point( edge->vertex0()->x(), edge->vertex0()->y() ));
polyline.points.push_back(Point( edge->vertex1()->x(), edge->vertex1()->y() ));
polyline.width.push_back(this->thickness[edge].first);
polyline.width.push_back(this->thickness[edge].second);
// remove this edge and its twin from the available edges
(void)this->edges.erase(edge);
(void)this->edges.erase(edge->twin());
// get next points
this->process_edge_neighbors(edge, &polyline);
// get previous points
{
ThickPolyline rpolyline;
this->process_edge_neighbors(edge->twin(), &rpolyline);
polyline.points.insert(polyline.points.begin(), rpolyline.points.rbegin(), rpolyline.points.rend());
polyline.width.insert(polyline.width.begin(), rpolyline.width.rbegin(), rpolyline.width.rend());
polyline.endpoints.first = rpolyline.endpoints.second;
}
assert(polyline.width.size() == polyline.points.size()*2 - 2);
// prevent loop endpoints from being extended
if (polyline.first_point().coincides_with(polyline.last_point())) {
polyline.endpoints.first = false;
polyline.endpoints.second = false;
}
// append polyline to result
polylines->push_back(polyline);
}
#ifdef SLIC3R_DEBUG
{
char path[2048];
static int iRun = 0;
sprintf(path, "out/MedialAxis-%d.svg", iRun ++);
dump_voronoi_to_svg(this->lines, this->vd, polylines, path);
printf("Thick lines: ");
for (ThickPolylines::const_iterator it = polylines->begin(); it != polylines->end(); ++ it) {
ThickLines lines = it->thicklines();
for (ThickLines::const_iterator it2 = lines.begin(); it2 != lines.end(); ++ it2) {
printf("%f,%f ", it2->a_width, it2->b_width);
}
}
printf("\n");
}
#endif /* SLIC3R_DEBUG */
}
void
MedialAxis::build(Polylines* polylines)
{
ThickPolylines tp;
this->build(&tp);
polylines->insert(polylines->end(), tp.begin(), tp.end());
}
void
MedialAxis::process_edge_neighbors(const VD::edge_type* edge, ThickPolyline* polyline)
{
while (true) {
// Since rot_next() works on the edge starting point but we want
// to find neighbors on the ending point, we just swap edge with
// its twin.
const VD::edge_type* twin = edge->twin();
// count neighbors for this edge
std::vector<const VD::edge_type*> neighbors;
for (const VD::edge_type* neighbor = twin->rot_next(); neighbor != twin;
neighbor = neighbor->rot_next()) {
if (this->valid_edges.count(neighbor) > 0) neighbors.push_back(neighbor);
}
// if we have a single neighbor then we can continue recursively
if (neighbors.size() == 1) {
const VD::edge_type* neighbor = neighbors.front();
// break if this is a closed loop
if (this->edges.count(neighbor) == 0) return;
Point new_point(neighbor->vertex1()->x(), neighbor->vertex1()->y());
polyline->points.push_back(new_point);
polyline->width.push_back(this->thickness[neighbor].first);
polyline->width.push_back(this->thickness[neighbor].second);
(void)this->edges.erase(neighbor);
(void)this->edges.erase(neighbor->twin());
edge = neighbor;
} else if (neighbors.size() == 0) {
polyline->endpoints.second = true;
return;
} else {
// T-shaped or star-shaped joint
return;
}
}
}
bool
MedialAxis::validate_edge(const VD::edge_type* edge)
{
// construct the line representing this edge of the Voronoi diagram
const Line line(
Point( edge->vertex0()->x(), edge->vertex0()->y() ),
Point( edge->vertex1()->x(), edge->vertex1()->y() )
);
// discard edge if it lies outside the supplied shape
// this could maybe be optimized (checking inclusion of the endpoints
// might give false positives as they might belong to the contour itself)
if (this->expolygon != NULL) {
if (line.a.coincides_with(line.b)) {
// in this case, contains(line) returns a false positive
if (!this->expolygon->contains(line.a)) return false;
} else {
if (!this->expolygon->contains(line)) return false;
}
}
// retrieve the original line segments which generated the edge we're checking
const VD::cell_type* cell1 = edge->cell();
const VD::cell_type* cell2 = edge->twin()->cell();
const Line &segment1 = this->retrieve_segment(cell1);
const Line &segment2 = this->retrieve_segment(cell2);
/* Calculate thickness of the section at both the endpoints of this edge.
Our Voronoi edge is part of a CCW sequence going around its Voronoi cell
(segment1). This edge's twin goes around segment2. Thus, segment2 is
oriented in the same direction as our main edge, and segment1 is oriented
in the same direction as our twin edge.
We used to only consider the (half-)distances to segment2, and that works
whenever segment1 and segment2 are almost specular and facing. However,
at curves they are staggered and they only face for a very little length
(such visibility actually coincides with our very short edge). This is why
we calculate w0 and w1 this way.
When cell1 or cell2 don't refer to the segment but only to an endpoint, we
calculate the distance to that endpoint instead. */
coordf_t w0 = cell2->contains_segment()
? line.a.perp_distance_to(segment2)*2
: line.a.distance_to(this->retrieve_endpoint(cell2))*2;
coordf_t w1 = cell1->contains_segment()
? line.b.perp_distance_to(segment1)*2
: line.b.distance_to(this->retrieve_endpoint(cell1))*2;
// if this edge is the centerline for a very thin area, we might want to skip it
// in case the area is too thin
if (w0 < SCALED_EPSILON || w1 < SCALED_EPSILON) {
if (cell1->contains_segment() && cell2->contains_segment()) {
// calculate the relative angle between the two boundary segments
double angle = fabs(segment2.orientation() - segment1.orientation());
// fabs(angle) ranges from 0 (collinear, same direction) to PI (collinear, opposite direction)
// we're interested only in segments close to the second case (facing segments)
// so we allow some tolerance.
// this filter ensures that we're dealing with a narrow/oriented area (longer than thick)
// we don't run it on edges not generated by two segments (thus generated by one segment
// and the endpoint of another segment), since their orientation would not be meaningful
if (fabs(angle - PI) > PI/5) return false;
} else {
return false;
}
}
if (w0 < this->min_width && w1 < this->min_width)
return false;
if (w0 > this->max_width && w1 > this->max_width)
return false;
this->thickness[edge] = std::make_pair(w0, w1);
this->thickness[edge->twin()] = std::make_pair(w1, w0);
return true;
}
const Line&
MedialAxis::retrieve_segment(const VD::cell_type* cell) const
{
return this->lines[cell->source_index()];
}
const Point&
MedialAxis::retrieve_endpoint(const VD::cell_type* cell) const
{
const Line& line = this->retrieve_segment(cell);
if (cell->source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT) {
return line.a;
} else {
return line.b;
}
}
} }
Reference in New Issue View Git Blame Copy Permalink