New functions for variable offsets of polygons / expolygons.

Test cases for the above.
Improvements of older test cases.

4 files changed, 328 insertions, 9 deletions

diff --git a/src/libslic3r/ClipperUtils.cpp b/src/libslic3r/ClipperUtils.cpp
index b863b4712..55bb4b446 100644
--- a/src/libslic3r/ClipperUtils.cpp
+++ b/src/libslic3r/ClipperUtils.cpp
@@ -107,8 +107,7 @@ void AddOuterPolyNodeToExPolygons(ClipperLib::PolyNode& polynode, ExPolygons* ex
   }
 }
  
-ExPolygons
-PolyTreeToExPolygons(ClipperLib::PolyTree& polytree)
+ExPolygons PolyTreeToExPolygons(ClipperLib::PolyTree& polytree)
 {
     ExPolygons retval;
     for (int i = 0; i < polytree.ChildCount(); ++i)
@@ -151,8 +150,7 @@ Slic3r::Polylines ClipperPaths_to_Slic3rPolylines(const ClipperLib::Paths &input
     return retval;
 }
 
-ExPolygons
-ClipperPaths_to_Slic3rExPolygons(const ClipperLib::Paths &input)
+ExPolygons ClipperPaths_to_Slic3rExPolygons(const ClipperLib::Paths &input)
 {
     // init Clipper
     ClipperLib::Clipper clipper;
@@ -167,8 +165,7 @@ ClipperPaths_to_Slic3rExPolygons(const ClipperLib::Paths &input)
     return PolyTreeToExPolygons(polytree);
 }
 
-ClipperLib::Path
-Slic3rMultiPoint_to_ClipperPath(const MultiPoint &input)
+ClipperLib::Path Slic3rMultiPoint_to_ClipperPath(const MultiPoint &input)
 {
     ClipperLib::Path retval;
     for (Points::const_iterator pit = input.points.begin(); pit != input.points.end(); ++pit)
@@ -176,8 +173,7 @@ Slic3rMultiPoint_to_ClipperPath(const MultiPoint &input)
     return retval;
 }
 
-ClipperLib::Path
-Slic3rMultiPoint_to_ClipperPath_reversed(const Slic3r::MultiPoint &input)
+ClipperLib::Path Slic3rMultiPoint_to_ClipperPath_reversed(const Slic3r::MultiPoint &input)
 {
     ClipperLib::Path output;
     output.reserve(input.points.size());
@@ -521,7 +517,7 @@ T _clipper_do(const ClipperLib::ClipType     clipType,
 
 // Fix of #117: A large fractal pyramid takes ages to slice
 // The Clipper library has difficulties processing overlapping polygons.
-// Namely, the function Clipper::JoinCommonEdges() has potentially a terrible time complexity if the output
+// Namely, the function ClipperLib::JoinCommonEdges() has potentially a terrible time complexity if the output
 // of the operation is of the PolyTree type.
 // This function implmenets a following workaround:
 // 1) Peform the Clipper operation with the output to Paths. This method handles overlaps in a reasonable time.
@@ -918,4 +914,304 @@ Polygons top_level_islands(const Slic3r::Polygons &polygons)
     return out;
 }
 
+// Outer offset shall not split the input contour into multiples. It is expected, that the solution will be non empty and it will contain just a single polygon.
+ClipperLib::Paths fix_after_outer_offset(const ClipperLib::Path &input, ClipperLib::PolyFillType filltype, bool reverse_result)
+{
+  	ClipperLib::Paths solution;
+  	if (! input.empty()) {
+		ClipperLib::Clipper clipper;
+	  	clipper.AddPath(input, ClipperLib::ptSubject, true);
+		clipper.ReverseSolution(reverse_result);
+		clipper.Execute(ClipperLib::ctUnion, solution, filltype, filltype);
+	}
+    return solution;
+}
+
+// Inner offset may split the source contour into multiple contours, but one shall not be inside the other.
+ClipperLib::Paths fix_after_inner_offset(const ClipperLib::Path &input, ClipperLib::PolyFillType filltype, bool reverse_result)
+{
+  	ClipperLib::Paths solution;
+  	if (! input.empty()) {
+		ClipperLib::Clipper clipper;
+		clipper.AddPath(input, ClipperLib::ptSubject, true);
+		ClipperLib::IntRect r = clipper.GetBounds();
+		r.left -= 10; r.top -= 10; r.right += 10; r.bottom += 10;
+		if (filltype == ClipperLib::pftPositive)
+			clipper.AddPath({ ClipperLib::IntPoint(r.left, r.bottom), ClipperLib::IntPoint(r.left, r.top), ClipperLib::IntPoint(r.right, r.top), ClipperLib::IntPoint(r.right, r.bottom) }, ClipperLib::ptSubject, true);
+		else
+			clipper.AddPath({ ClipperLib::IntPoint(r.left, r.bottom), ClipperLib::IntPoint(r.right, r.bottom), ClipperLib::IntPoint(r.right, r.top), ClipperLib::IntPoint(r.left, r.top) }, ClipperLib::ptSubject, true);
+		clipper.ReverseSolution(reverse_result);
+		clipper.Execute(ClipperLib::ctUnion, solution, filltype, filltype);
+		if (! solution.empty())
+			solution.erase(solution.begin());
+	}
+	return solution;
+}
+
+ClipperLib::Path mittered_offset_path_scaled(const Points &contour, const std::vector<float> &deltas, double miter_limit)
+{
+	assert(contour.size() == deltas.size());
+#ifndef NDEBUG
+	// Verify that the deltas are either all positive, or all negative.
+	bool positive = false;
+	bool negative = false;
+	for (float delta : deltas)
+		if (delta < 0.f)
+			negative = true;
+		else if (delta > 0.f)
+			positive = true;
+	assert(! (negative && positive));
+#endif /* NDEBUG */
+
+	ClipperLib::Path out;
+
+	if (deltas.size() > 2)
+	{
+		out.reserve(contour.size() * 2);
+
+		// Clamp miter limit to 2.
+		miter_limit = (miter_limit > 2.) ? 2. / (miter_limit * miter_limit) : 0.5;
+		
+		// perpenduclar vector
+		auto   perp = [](const Vec2d &v) -> Vec2d { return Vec2d(v.y(), - v.x()); };
+
+		// Add a new point to the output, scale by CLIPPER_OFFSET_SCALE and round to ClipperLib::cInt.
+		auto   add_offset_point = [&out](Vec2d pt) {
+			pt *= double(CLIPPER_OFFSET_SCALE);
+			pt += Vec2d(0.5 - (pt.x() < 0), 0.5 - (pt.y() < 0));
+			out.emplace_back(ClipperLib::cInt(pt.x()), ClipperLib::cInt(pt.y()));
+		};
+
+		// Minimum edge length, squared.
+		double lmin  = *std::max_element(deltas.begin(), deltas.end()) * CLIPPER_OFFSET_SHORTEST_EDGE_FACTOR;
+		double l2min = lmin * lmin;
+		// Minimum angle to consider two edges to be parallel.
+		double sin_min_parallel = EPSILON + 1. / double(CLIPPER_OFFSET_SCALE);
+
+		// Find the last point further from pt by l2min.
+		Vec2d  pt     = contour.front().cast<double>();
+		size_t iprev  = contour.size() - 1;
+		Vec2d  ptprev;
+		for (; iprev > 0; -- iprev) {
+			ptprev = contour[iprev].cast<double>();
+			if ((ptprev - pt).squaredNorm() > l2min)
+				break;
+		}
+
+		if (iprev != 0) {
+			size_t ilast = iprev;
+			// Normal to the (pt - ptprev) segment.
+			Vec2d nprev = perp(pt - ptprev).normalized();
+			for (size_t i = 0; ; ) {
+				// Find the next point further from pt by l2min.
+				size_t j = i + 1;
+				Vec2d ptnext;
+				for (; j <= ilast; ++ j) {
+					ptnext = contour[j].cast<double>();
+					double l2 = (ptnext - pt).squaredNorm();
+					if (l2 > l2min)
+						break;
+				}
+				if (j > ilast)
+					ptnext = contour.front().cast<double>();
+
+				// Normal to the (ptnext - pt) segment.
+				Vec2d nnext  = perp(ptnext - pt).normalized();
+
+				double delta  = deltas[i];
+				double sin_a  = clamp(-1., 1., cross2(nprev, nnext));
+				double convex = sin_a * delta;
+				if (convex <= - sin_min_parallel) {
+					// Concave corner.
+					add_offset_point(pt + nprev * delta);
+					add_offset_point(pt);
+					add_offset_point(pt + nnext * delta);
+				} else if (convex < sin_min_parallel) {
+					// Nearly parallel.
+					add_offset_point((nprev.dot(nnext) > 0.) ? (pt + nprev * delta) : pt);
+				} else {
+					// Convex corner
+					double dot = nprev.dot(nnext);
+					double r = 1. + dot;
+				  	if (r >= miter_limit)
+						add_offset_point(pt + (nprev + nnext) * (delta / r));
+				  	else {
+						double dx = std::tan(std::atan2(sin_a, dot) / 4.);
+						Vec2d  newpt1 = pt + (nprev - perp(nprev) * dx) * delta;
+						Vec2d  newpt2 = pt + (nnext + perp(nnext) * dx) * delta;
+#ifndef NDEBUG
+						Vec2d vedge = 0.5 * (newpt1 + newpt2) - pt;
+						double dist_norm = vedge.norm();
+						assert(std::abs(dist_norm - delta) < EPSILON);
+#endif /* NDEBUG */
+						add_offset_point(newpt1);
+						add_offset_point(newpt2);
+				  	}
+				}
+
+				if (i == ilast)
+					break;
+
+				ptprev = pt;
+				nprev  = nnext;
+				pt     = ptnext;
+				i = j;
+			}
+		}
+	}
+
+	return out;
+}
+
+Polygons variable_offset_inner(const ExPolygon &expoly, const std::vector<std::vector<float>> &deltas, double miter_limit)
+{
+#ifndef NDEBUG
+	// Verify that the deltas are all non positive.
+	for (const std::vector<float> &ds : deltas)
+		for (float delta : ds)
+			assert(delta <= 0.);
+	assert(expoly.holes.size() + 1 == deltas.size());
+#endif /* NDEBUG */
+
+	// 1) Offset the outer contour.
+	ClipperLib::Paths contours = fix_after_inner_offset(mittered_offset_path_scaled(expoly.contour.points, deltas.front(), miter_limit), ClipperLib::pftNegative, true);
+
+	// 2) Offset the holes one by one, collect the results.
+	ClipperLib::Paths holes;
+	holes.reserve(expoly.holes.size());
+	for (const Polygon& hole : expoly.holes)
+		append(holes, fix_after_outer_offset(mittered_offset_path_scaled(hole, deltas[1 + &hole - expoly.holes.data()], miter_limit), ClipperLib::pftPositive, false));
+
+	// 3) Subtract holes from the contours.
+	ClipperLib::Paths output;
+	if (holes.empty())
+		output = std::move(contours);
+	else {
+		ClipperLib::Clipper clipper;
+		clipper.Clear();
+		clipper.AddPaths(contours, ClipperLib::ptSubject, true);
+		clipper.AddPaths(holes, ClipperLib::ptClip, true);
+		clipper.Execute(ClipperLib::ctDifference, output, ClipperLib::pftNonZero, ClipperLib::pftNonZero);
+	}
+
+	// 4) Unscale the output.
+	unscaleClipperPolygons(output);
+	return ClipperPaths_to_Slic3rPolygons(output);
+}
+
+Polygons variable_offset_outer(const ExPolygon &expoly, const std::vector<std::vector<float>> &deltas, double miter_limit)
+{
+#ifndef NDEBUG
+	// Verify that the deltas are all non positive.
+for (const std::vector<float>& ds : deltas)
+		for (float delta : ds)
+			assert(delta >= 0.);
+	assert(expoly.holes.size() + 1 == deltas.size());
+#endif /* NDEBUG */
+
+	// 1) Offset the outer contour.
+	ClipperLib::Paths contours = fix_after_outer_offset(mittered_offset_path_scaled(expoly.contour.points, deltas.front(), miter_limit), ClipperLib::pftPositive, false);
+
+	// 2) Offset the holes one by one, collect the results.
+	ClipperLib::Paths holes;
+	holes.reserve(expoly.holes.size());
+	for (const Polygon& hole : expoly.holes)
+		append(holes, fix_after_inner_offset(mittered_offset_path_scaled(hole, deltas[1 + &hole - expoly.holes.data()], miter_limit), ClipperLib::pftPositive, true));
+
+	// 3) Subtract holes from the contours.
+	ClipperLib::Paths output;
+	if (holes.empty())
+		output = std::move(contours);
+	else {
+		ClipperLib::Clipper clipper;
+		clipper.Clear();
+		clipper.AddPaths(contours, ClipperLib::ptSubject, true);
+		clipper.AddPaths(holes, ClipperLib::ptClip, true);
+		clipper.Execute(ClipperLib::ctDifference, output, ClipperLib::pftNonZero, ClipperLib::pftNonZero);
+	}
+
+	// 4) Unscale the output.
+	unscaleClipperPolygons(output);
+	return ClipperPaths_to_Slic3rPolygons(output);
+}
+
+ExPolygons variable_offset_outer_ex(const ExPolygon &expoly, const std::vector<std::vector<float>> &deltas, double miter_limit)
+{
+#ifndef NDEBUG
+	// Verify that the deltas are all non positive.
+for (const std::vector<float>& ds : deltas)
+		for (float delta : ds)
+			assert(delta >= 0.);
+	assert(expoly.holes.size() + 1 == deltas.size());
+#endif /* NDEBUG */
+
+	// 1) Offset the outer contour.
+	ClipperLib::Paths contours = fix_after_outer_offset(mittered_offset_path_scaled(expoly.contour.points, deltas.front(), miter_limit), ClipperLib::pftPositive, false);
+
+	// 2) Offset the holes one by one, collect the results.
+	ClipperLib::Paths holes;
+	holes.reserve(expoly.holes.size());
+	for (const Polygon& hole : expoly.holes)
+		append(holes, fix_after_inner_offset(mittered_offset_path_scaled(hole, deltas[1 + &hole - expoly.holes.data()], miter_limit), ClipperLib::pftPositive, true));
+
+	// 3) Subtract holes from the contours.
+	unscaleClipperPolygons(contours);
+	ExPolygons output;
+	if (holes.empty()) {
+		output.reserve(contours.size());
+		for (ClipperLib::Path &path : contours) 
+			output.emplace_back(ClipperPath_to_Slic3rPolygon(path));
+	} else {
+		ClipperLib::Clipper clipper;
+		unscaleClipperPolygons(holes);
+		clipper.AddPaths(contours, ClipperLib::ptSubject, true);
+		clipper.AddPaths(holes, ClipperLib::ptClip, true);
+	    ClipperLib::PolyTree polytree;
+		clipper.Execute(ClipperLib::ctDifference, polytree, ClipperLib::pftNonZero, ClipperLib::pftNonZero);
+	    output = PolyTreeToExPolygons(polytree);
+	}
+
+	return output;
+}
+
+
+ExPolygons variable_offset_inner_ex(const ExPolygon &expoly, const std::vector<std::vector<float>> &deltas, double miter_limit)
+{
+#ifndef NDEBUG
+	// Verify that the deltas are all non positive.
+for (const std::vector<float>& ds : deltas)
+		for (float delta : ds)
+			assert(delta <= 0.);
+	assert(expoly.holes.size() + 1 == deltas.size());
+#endif /* NDEBUG */
+
+	// 1) Offset the outer contour.
+	ClipperLib::Paths contours = fix_after_inner_offset(mittered_offset_path_scaled(expoly.contour.points, deltas.front(), miter_limit), ClipperLib::pftNegative, false);
+
+	// 2) Offset the holes one by one, collect the results.
+	ClipperLib::Paths holes;
+	holes.reserve(expoly.holes.size());
+	for (const Polygon& hole : expoly.holes)
+		append(holes, fix_after_outer_offset(mittered_offset_path_scaled(hole, deltas[1 + &hole - expoly.holes.data()], miter_limit), ClipperLib::pftNegative, true));
+
+	// 3) Subtract holes from the contours.
+	unscaleClipperPolygons(contours);
+	ExPolygons output;
+	if (holes.empty()) {
+		output.reserve(contours.size());
+		for (ClipperLib::Path &path : contours) 
+			output.emplace_back(ClipperPath_to_Slic3rPolygon(path));
+	} else {
+		ClipperLib::Clipper clipper;
+		unscaleClipperPolygons(holes);
+		clipper.AddPaths(contours, ClipperLib::ptSubject, true);
+		clipper.AddPaths(holes, ClipperLib::ptClip, true);
+	    ClipperLib::PolyTree polytree;
+		clipper.Execute(ClipperLib::ctDifference, polytree, ClipperLib::pftNonZero, ClipperLib::pftNonZero);
+	    output = PolyTreeToExPolygons(polytree);
+	}
+
+	return output;
+}
+
 }
diff --git a/src/libslic3r/ClipperUtils.hpp b/src/libslic3r/ClipperUtils.hpp
index d8f8a8f94..5a41a6a90 100644
--- a/src/libslic3r/ClipperUtils.hpp
+++ b/src/libslic3r/ClipperUtils.hpp
@@ -238,6 +238,11 @@ void safety_offset(ClipperLib::Paths* paths);
 
 Polygons top_level_islands(const Slic3r::Polygons &polygons);
 
+Polygons  variable_offset_inner(const ExPolygon &expoly, const std::vector<std::vector<float>> &deltas, double miter_limit = 2.);
+Polygons  variable_offset_outer(const ExPolygon &expoly, const std::vector<std::vector<float>> &deltas, double miter_limit = 2.);
+ExPolygons variable_offset_outer_ex(const ExPolygon &expoly, const std::vector<std::vector<float>> &deltas, double miter_limit = 2.);
+ExPolygons variable_offset_inner_ex(const ExPolygon &expoly, const std::vector<std::vector<float>> &deltas, double miter_limit = 2.);
+
 }
 
 #endif
diff --git a/src/libslic3r/ExPolygon.hpp b/src/libslic3r/ExPolygon.hpp
index c510b848f..08c4f7a07 100644
--- a/src/libslic3r/ExPolygon.hpp
+++ b/src/libslic3r/ExPolygon.hpp
@@ -301,6 +301,15 @@ inline bool expolygons_contain(ExPolygons &expolys, const Point &pt)
     return false;
 }
 
+inline ExPolygons expolygons_simplify(const ExPolygons &expolys, double tolerance)
+{
+	ExPolygons out;
+	out.reserve(expolys.size());
+	for (const ExPolygon &exp : expolys)
+		exp.simplify(tolerance, &out);
+	return out;
+}
+
 extern BoundingBox get_extents(const ExPolygon &expolygon);
 extern BoundingBox get_extents(const ExPolygons &expolygons);
 extern BoundingBox get_extents_rotated(const ExPolygon &poly, double angle);
diff --git a/src/libslic3r/Polygon.hpp b/src/libslic3r/Polygon.hpp
index 19be3068b..ad1b32feb 100644
--- a/src/libslic3r/Polygon.hpp
+++ b/src/libslic3r/Polygon.hpp
@@ -102,6 +102,15 @@ inline void        polygons_append(Polygons &dst, Polygons &&src)
     }
 }
 
+inline Polygons polygons_simplify(const Polygons &polys, double tolerance)
+{
+	Polygons out;
+	out.reserve(polys.size());
+	for (const Polygon &p : polys)
+		polygons_append(out, p.simplify(tolerance));
+	return out;
+}
+
 inline void polygons_rotate(Polygons &polys, double angle)
 {
     const double cos_angle = cos(angle);