More refactoring to medial axis and gap fill, more robust

This commit is contained in:
Alessandro Ranellucci 2016-05-20 06:24:05 +02:00
parent b068616366
commit 9e8022f6f6
13 changed files with 133 additions and 58 deletions

View File

@ -1,4 +1,4 @@
use Test::More tests => 21;
use Test::More tests => 23;
use strict;
use warnings;
@ -8,7 +8,7 @@ BEGIN {
}
use Slic3r;
use List::Util qw(first sum);
use List::Util qw(first sum none);
use Slic3r::Geometry qw(epsilon scale unscale scaled_epsilon Y);
use Slic3r::Test;
@ -91,7 +91,20 @@ if (0) {
is scalar(@$res), 1, 'medial axis of a narrow rectangle with an extra vertex is still a single line';
ok unscale($res->[0]->length) >= (200-100 - (120-100)) - epsilon, 'medial axis has still a reasonable length';
ok !(grep { abs($_ - scale 150) < scaled_epsilon } map $_->[Y], map @$_, @$res2), "extra vertices don't influence medial axis";
}
{
my $expolygon = Slic3r::ExPolygon->new(
Slic3r::Polygon->new([1185881,829367],[1421988,1578184],[1722442,2303558],[2084981,2999998],[2506843,3662186],[2984809,4285086],[3515250,4863959],[4094122,5394400],[4717018,5872368],[5379210,6294226],[6075653,6656769],[6801033,6957229],[7549842,7193328],[8316383,7363266],[9094809,7465751],[9879211,7500000],[10663611,7465750],[11442038,7363265],[12208580,7193327],[12957389,6957228],[13682769,6656768],[14379209,6294227],[15041405,5872366],[15664297,5394401],[16243171,4863960],[16758641,4301424],[17251579,3662185],[17673439,3000000],[18035980,2303556],[18336441,1578177],[18572539,829368],[18750748,0],[19758422,0],[19727293,236479],[19538467,1088188],[19276136,1920196],[18942292,2726179],[18539460,3499999],[18070731,4235755],[17539650,4927877],[16950279,5571067],[16307090,6160437],[15614974,6691519],[14879209,7160248],[14105392,7563079],[13299407,7896927],[12467399,8159255],[11615691,8348082],[10750769,8461952],[9879211,8500000],[9007652,8461952],[8142729,8348082],[7291022,8159255],[6459015,7896927],[5653029,7563079],[4879210,7160247],[4143447,6691519],[3451331,6160437],[2808141,5571066],[2218773,4927878],[1687689,4235755],[1218962,3499999],[827499,2748020],[482284,1920196],[219954,1088186],[31126,236479],[0,0],[1005754,0]),
);
my $res = $expolygon->medial_axis(scale 1.324888, scale 0.25);
is scalar(@$res), 1, 'medial axis of a semicircumference is a single line';
# check whether turns are all CCW or all CW
my @lines = @{$res->[0]->lines};
my @angles = map { $lines[$_-1]->ccw($lines[$_]->b) } 1..$#lines;
ok !!(none { $_ < 0 } @angles) || (none { $_ > 0 } @angles),
'all medial axis segments of a semicircumference have the same orientation';
}
{
@ -136,7 +149,7 @@ if (0) {
{
my $expolygon = Slic3r::ExPolygon->new(Slic3r::Polygon->new_scale(
[50, 100],
[300, 102],
[1000, 102],
[50, 104],
));
my $res = $expolygon->medial_axis(scale 4, scale 0.5);

View File

@ -174,17 +174,10 @@ void
ExPolygon::medial_axis(double max_width, double min_width, ThickPolylines* polylines) const
{
// init helper object
Slic3r::Geometry::MedialAxis ma(max_width, min_width);
ma.expolygon = this;
Slic3r::Geometry::MedialAxis ma(max_width, min_width, this);
ma.lines = this->lines();
// populate list of segments for the Voronoi diagram
ma.lines = this->contour.lines();
for (Polygons::const_iterator hole = this->holes.begin(); hole != this->holes.end(); ++hole) {
Lines lines = hole->lines();
ma.lines.insert(ma.lines.end(), lines.begin(), lines.end());
}
// compute the Voronoi diagram
// compute the Voronoi diagram and extract medial axis polylines
ThickPolylines pp;
ma.build(&pp);
@ -195,15 +188,14 @@ ExPolygon::medial_axis(double max_width, double min_width, ThickPolylines* polyl
svg.Close();
*/
// find the maximum width returned
/* Find the maximum width returned; we're going to use this for validating and
filtering the output segments. */
double max_w = 0;
if (!pp.empty()) {
std::vector<double> widths;
for (ThickPolylines::const_iterator it = pp.begin(); it != pp.end(); ++it)
widths.insert(widths.end(), it->width.begin(), it->width.end());
max_w = *std::max_element(widths.begin(), widths.end());
}
for (ThickPolylines::const_iterator it = pp.begin(); it != pp.end(); ++it)
max_w = fmaxf(max_w, *std::max_element(it->width.begin(), it->width.end()));
/* Loop through all returned polylines in order to extend their endpoints to the
expolygon boundaries */
bool removed = false;
for (size_t i = 0; i < pp.size(); ++i) {
ThickPolyline& polyline = pp[i];
@ -502,4 +494,15 @@ ExPolygon::lines() const
return lines;
}
std::string
ExPolygon::dump_perl() const
{
std::ostringstream ret;
ret << "[" << this->contour.dump_perl();
for (Polygons::const_iterator h = this->holes.begin(); h != this->holes.end(); ++h)
ret << "," << h->dump_perl();
ret << "]";
return ret.str();
}
}

View File

@ -42,6 +42,7 @@ class ExPolygon
void triangulate_pp(Polygons* polygons) const;
void triangulate_p2t(Polygons* polygons) const;
Lines lines() const;
std::string dump_perl() const;
};
}

View File

@ -324,7 +324,7 @@ MedialAxis::build(ThickPolylines* polylines)
if (edge->is_secondary() || edge->is_infinite()) continue;
// don't re-validate twins
if (seen_edges.find(&*edge) != seen_edges.end()) continue;
if (seen_edges.find(&*edge) != seen_edges.end()) continue; // TODO: is this needed?
seen_edges.insert(&*edge);
seen_edges.insert(edge->twin());
@ -445,49 +445,66 @@ MedialAxis::validate_edge(const VD::edge_type* edge)
}
// 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);
const VD::cell_type* cell_l = edge->cell();
const VD::cell_type* cell_r = edge->twin()->cell();
const Line &segment_l = this->retrieve_segment(cell_l);
const Line &segment_r = this->retrieve_segment(cell_r);
/* Calculate thickness of the section at both the endpoints of this edge.
/*
SVG svg("edge.svg");
svg.draw(*this->expolygon);
svg.draw(line);
svg.draw(segment_l, "red");
svg.draw(segment_r, "blue");
svg.Close();
*/
/* Calculate thickness of the cross-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
located on the left side. (segment_l).
This edge's twin goes around segment_r. Thus, segment_r is
oriented in the same direction as our main edge, and segment_l 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,
We used to only consider the (half-)distances to segment_r, and that works
whenever segment_l and segment_r 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
(our very short edge represents such visibility).
Both w0 and w1 can be calculated either towards cell_l or cell_r with equal
results by Voronoi definition.
When cell_l or cell_r 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 w0 = cell_r->contains_segment()
? line.a.distance_to(segment_r)*2
: line.a.distance_to(this->retrieve_endpoint(cell_r))*2;
coordf_t w1 = cell1->contains_segment()
? line.b.perp_distance_to(segment1)*2
: line.b.distance_to(this->retrieve_endpoint(cell1))*2;
coordf_t w1 = cell_l->contains_segment()
? line.b.distance_to(segment_l)*2
: line.b.distance_to(this->retrieve_endpoint(cell_l))*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());
if (cell_l->contains_segment() && cell_r->contains_segment()) {
// calculate the relative angle between the two boundary segments
double angle = fabs(segment_r.orientation() - segment_l.orientation());
if (angle > PI) angle = 2*PI - angle;
assert(angle >= 0 && angle <= PI);
// 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;
// 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 (PI - angle > PI/8) {
// angle is not narrow enough
// only apply this filter to segments that are not too short otherwise their
// angle could possibly be not meaningful
if (w0 < SCALED_EPSILON || w1 < SCALED_EPSILON || line.length() >= this->min_width)
return false;
}
} else {
if (w0 < SCALED_EPSILON || w1 < SCALED_EPSILON)
return false;
}
if (w0 < this->min_width && w1 < this->min_width)

View File

@ -46,8 +46,8 @@ class MedialAxis {
const ExPolygon* expolygon;
double max_width;
double min_width;
MedialAxis(double _max_width, double _min_width)
: max_width(_max_width), min_width(_min_width), expolygon(NULL) {};
MedialAxis(double _max_width, double _min_width, const ExPolygon* _expolygon = NULL)
: max_width(_max_width), min_width(_min_width), expolygon(_expolygon) {};
void build(ThickPolylines* polylines);
void build(Polylines* polylines);

View File

@ -212,6 +212,12 @@ Line::intersection(const Line& line, Point* intersection) const
return false; // not intersecting
}
double
Line::ccw(const Point& point) const
{
return point.ccw(*this);
}
Pointf3
Linef3::intersect_plane(double z) const
{

View File

@ -44,6 +44,7 @@ class Line
void extend_end(double distance);
void extend_start(double distance);
bool intersection(const Line& line, Point* intersection) const;
double ccw(const Point& point) const;
};
class ThickLine : public Line

View File

@ -128,6 +128,19 @@ MultiPoint::intersection(const Line& line, Point* intersection) const
return false;
}
std::string
MultiPoint::dump_perl() const
{
std::ostringstream ret;
ret << "[";
for (Points::const_iterator p = this->points.begin(); p != this->points.end(); ++p) {
ret << p->dump_perl();
if (p != this->points.end()-1) ret << ",";
}
ret << "]";
return ret.str();
}
Points
MultiPoint::_douglas_peucker(const Points &points, const double tolerance)
{

View File

@ -37,6 +37,7 @@ class MultiPoint
void append(const Points &points);
void append(const Points::const_iterator &begin, const Points::const_iterator &end);
bool intersection(const Line& line, Point* intersection) const;
std::string dump_perl() const;
static Points _douglas_peucker(const Points &points, const double tolerance);
};

View File

@ -92,7 +92,7 @@ PerimeterGenerator::process()
// the following offset2 ensures almost nothing in @thin_walls is narrower than $min_width
// (actually, something larger than that still may exist due to mitering or other causes)
coord_t min_width = scale_(this->ext_perimeter_flow.nozzle_diameter / 2);
coord_t min_width = scale_(this->ext_perimeter_flow.nozzle_diameter / 3);
ExPolygons expp = offset2_ex(diffpp, -min_width/2, +min_width/2);
// the maximum thickness of our thin wall area is equal to the minimum thickness of a single loop

View File

@ -26,6 +26,14 @@ Point::wkt() const
return ss.str();
}
std::string
Point::dump_perl() const
{
std::ostringstream ss;
ss << "[" << this->x << "," << this->y << "]";
return ss.str();
}
void
Point::scale(double factor)
{
@ -313,6 +321,14 @@ Pointf::wkt() const
return ss.str();
}
std::string
Pointf::dump_perl() const
{
std::ostringstream ss;
ss << "[" << this->x << "," << this->y << "]";
return ss.str();
}
void
Pointf::scale(double factor)
{

View File

@ -38,6 +38,7 @@ class Point
};
bool operator==(const Point& rhs) const;
std::string wkt() const;
std::string dump_perl() const;
void scale(double factor);
void translate(double x, double y);
void translate(const Vector &vector);
@ -87,6 +88,7 @@ class Pointf
return Pointf(unscale(p.x), unscale(p.y));
};
std::string wkt() const;
std::string dump_perl() const;
void scale(double factor);
void translate(double x, double y);
void translate(const Vectorf &vector);

View File

@ -41,6 +41,8 @@
%code{% RETVAL = Polyline(*THIS); %};
Clone<Point> normal();
Clone<Point> vector();
double ccw(Point* point)
%code{% RETVAL = THIS->ccw(*point); %};
%{
Line*