ClipperLib: Optimized PointInPolygon() to calculate cross products

with int64s instead of doubles.
Polygon / ExPolygon: contains() reworked to use ClipperLib::PointInPolygon().
	The Slic3r own implementation was not robust.
Fixed test_perimeters after recent refactoring (sorting of extrusions
	into LayerIslands)

src/clipper/clipper.cpp

@@ -225,7 +225,7 @@ int PointInPolygon(const IntPoint &pt, const Path &path)
         if (ipNext.x() > pt.x()) result = 1 - result;
         else
         {
-          double d = (double)(ip.x() - pt.x()) * (ipNext.y() - pt.y()) - (double)(ipNext.x() - pt.x()) * (ip.y() - pt.y());
+          auto d = CrossProductType(ip.x() - pt.x()) * CrossProductType(ipNext.y() - pt.y()) - CrossProductType(ipNext.x() - pt.x()) * CrossProductType(ip.y() - pt.y());
           if (!d) return -1;
           if ((d > 0) == (ipNext.y() > ip.y())) result = 1 - result;
         }
@@ -233,7 +233,7 @@ int PointInPolygon(const IntPoint &pt, const Path &path)
       {
         if (ipNext.x() > pt.x())
         {
-          double d = (double)(ip.x() - pt.x()) * (ipNext.y() - pt.y()) - (double)(ipNext.x() - pt.x()) * (ip.y() - pt.y());
+          auto d = CrossProductType(ip.x() - pt.x()) * CrossProductType(ipNext.y() - pt.y()) - CrossProductType(ipNext.x() - pt.x()) * CrossProductType(ip.y() - pt.y());
           if (!d) return -1;
           if ((d > 0) == (ipNext.y() > ip.y())) result = 1 - result;
         }
@@ -246,7 +246,7 @@ int PointInPolygon(const IntPoint &pt, const Path &path)
 //------------------------------------------------------------------------------
 
 // Called by Poly2ContainsPoly1()
-int PointInPolygon (const IntPoint &pt, OutPt *op)
+int PointInPolygon(const IntPoint &pt, OutPt *op)
 {
   //returns 0 if false, +1 if true, -1 if pt ON polygon boundary
   int result = 0;
@@ -265,7 +265,7 @@ int PointInPolygon (const IntPoint &pt, OutPt *op)
         if (op->Next->Pt.x() > pt.x()) result = 1 - result;
         else
         {
-          double d = (double)(op->Pt.x() - pt.x()) * (op->Next->Pt.y() - pt.y()) - (double)(op->Next->Pt.x() - pt.x()) * (op->Pt.y() - pt.y());
+          auto d = CrossProductType(op->Pt.x() - pt.x()) * CrossProductType(op->Next->Pt.y() - pt.y()) - CrossProductType(op->Next->Pt.x() - pt.x()) * CrossProductType(op->Pt.y() - pt.y());
           if (!d) return -1;
           if ((d > 0) == (op->Next->Pt.y() > op->Pt.y())) result = 1 - result;
         }
@@ -273,7 +273,7 @@ int PointInPolygon (const IntPoint &pt, OutPt *op)
       {
         if (op->Next->Pt.x() > pt.x())
         {
-          double d = (double)(op->Pt.x() - pt.x()) * (op->Next->Pt.y() - pt.y()) - (double)(op->Next->Pt.x() - pt.x()) * (op->Pt.y() - pt.y());
+          auto d = CrossProductType(op->Pt.x() - pt.x()) * CrossProductType(op->Next->Pt.y() - pt.y()) - CrossProductType(op->Next->Pt.x() - pt.x()) * CrossProductType(op->Pt.y() - pt.y());
           if (!d) return -1;
           if ((d > 0) == (op->Next->Pt.y() > op->Pt.y())) result = 1 - result;
         }

src/clipper/clipper.hpp

@@ -85,9 +85,11 @@ enum PolyFillType { pftEvenOdd, pftNonZero, pftPositive, pftNegative };
 // Point coordinate type
 #ifdef CLIPPERLIB_INT32
   // Coordinates and their differences (vectors of the coordinates) have to fit int32_t.
-  typedef int32_t cInt;
+  using cInt = int32_t;
+  using CrossProductType = int64_t;
 #else
-  typedef int64_t cInt;
+  using cInt = int64_t;
+  using CrossProductType = double;
   // Maximum cInt value to allow a cross product calculation using 32bit expressions.
   static constexpr cInt const loRange = 0x3FFFFFFF; // 0x3FFFFFFF = 1 073 741 823
   // Maximum allowed cInt value.

src/libslic3r/ExPolygon.cpp

@@ -100,10 +100,12 @@ bool ExPolygon::contains(const Polylines &polylines) const
 
 bool ExPolygon::contains(const Point &point) const
 {
-    if (! this->contour.contains(point))
+    if (! Slic3r::contains(contour, point, true))
+        // Outside the outer contour, not on the contour boundary.
         return false;
     for (const Polygon &hole : this->holes)
-        if (hole.contains(point))
+        if (Slic3r::contains(hole, point, false))
+            // Inside a hole, not on the hole boundary.
             return false;
     return true;
 }
@@ -114,13 +116,13 @@ bool ExPolygon::contains_b(const Point &point) const
     return this->contains(point) || this->has_boundary_point(point);
 }
 
-bool
+bool ExPolygon::has_boundary_point(const Point &point) const
-ExPolygon::has_boundary_point(const Point &point) const
 {
-    if (this->contour.has_boundary_point(point)) return true;
+    if (this->contour.has_boundary_point(point))
-    for (Polygons::const_iterator h = this->holes.begin(); h != this->holes.end(); ++h) {
+        return true;
-        if (h->has_boundary_point(point)) return true;
+    for (const Polygon &hole : this->holes)
-    }
+        if (hole.has_boundary_point(point))
+            return true;
     return false;
 }
 

src/libslic3r/ExPolygon.hpp

@@ -373,14 +373,14 @@ inline void expolygons_append(ExPolygons &dst, ExPolygons &&src)
 
 inline void expolygons_rotate(ExPolygons &expolys, double angle)
 {
-    for (ExPolygons::iterator p = expolys.begin(); p != expolys.end(); ++p)
+    for (ExPolygon &expoly : expolys)
-        p->rotate(angle);
+        expoly.rotate(angle);
 }
 
 inline bool expolygons_contain(ExPolygons &expolys, const Point &pt)
 {
-    for (ExPolygons::iterator p = expolys.begin(); p != expolys.end(); ++p)
+    for (const ExPolygon &expoly : expolys)
-        if (p->contains(pt))
+        if (expoly.contains(pt))
             return true;
     return false;
 }

src/libslic3r/Layer.cpp

@@ -691,7 +691,7 @@ void Layer::sort_perimeters_into_islands(
     for (auto it_source_slice  = perimeter_slices_queue.begin(); it_source_slice != perimeter_slices_queue.end(); ++ it_source_slice) {
         double d2min = std::numeric_limits<double>::max();
         int    lslice_idx_min = -1;
-        for (int lslice_idx = int(lslices_ex.size()) - 1; lslice_idx >= 0 && ! perimeter_slices_queue.empty(); -- lslice_idx)
+        for (int lslice_idx = int(lslices_ex.size()) - 1; lslice_idx >= 0; -- lslice_idx)
             if (double d2 = point_inside_surface_dist2(lslice_idx, it_source_slice->second); d2 < d2min) {
                 d2min = d2;
                 lslice_idx_min = lslice_idx;

src/libslic3r/Polygon.cpp

@@ -89,26 +89,9 @@ void Polygon::douglas_peucker(double tolerance)
 }
 
 // Does an unoriented polygon contain a point?
-// Tested by counting intersections along a horizontal line.
 bool Polygon::contains(const Point &p) const
 {
-    // http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
+    return Slic3r::contains(*this, p, true);
-    bool result = false;
-    Points::const_iterator i = this->points.begin();
-    Points::const_iterator j = this->points.end() - 1;
-    for (; i != this->points.end(); j = i ++)
-        if (i->y() > p.y() != j->y() > p.y())
-#if 1
-            if (Vec2d v = (*j - *i).cast<double>();
-                // p.x() is below the line
-                p.x() - i->x() < double(p.y() - i->y()) * v.x() / v.y())
-#else
-            // Orientation predicated relative to i-th point.
-            if (double orient = (double)(p.x() - i->x()) * (double)(j->y() - i->y()) - (double)(p.y() - i->y()) * (double)(j->x() - i->x());
-                (i->y() > j->y()) ? (orient > 0.) : (orient < 0.))
-#endif
-                result = !result;
-    return result;
 }
 
 // this only works on CCW polygons as CW will be ripped out by Clipper's simplify_polygons()
@@ -536,6 +519,15 @@ bool polygons_match(const Polygon &l, const Polygon &r)
     return true;
 }
 
+bool contains(const Polygon &polygon, const Point &p, bool border_result)
+{
+    if (const int poly_count_inside = ClipperLib::PointInPolygon(p, polygon.points); 
+        poly_count_inside == -1)
+        return border_result;
+    else
+        return (poly_count_inside % 2) == 1;
+}
+
 bool contains(const Polygons &polygons, const Point &p, bool border_result)
 {
     int poly_count_inside = 0;

src/libslic3r/Polygon.hpp

@@ -255,6 +255,7 @@ inline Polygons to_polygons(std::vector<Points> &&paths)
 bool polygons_match(const Polygon &l, const Polygon &r);
 
 // Returns true if inside. Returns border_result if on boundary.
+bool contains(const Polygon& polygon, const Point& p, bool border_result = true);
 bool contains(const Polygons& polygons, const Point& p, bool border_result = true);
 
 Polygon make_circle(double radius, double error);

tests/fff_print/test_perimeters.cpp

@@ -470,6 +470,8 @@ SCENARIO("Some weird coverage test", "[Perimeters]")
     LayerRegion *layerm = layer->get_region(0);
     layerm->m_slices.clear();
     layerm->m_slices.append({ expolygon }, stInternal);
+    layer->lslices = { expolygon };
+    layer->lslices_ex = { { get_extents(expolygon) } };
     
     // make perimeters
     layer->make_perimeters();