I have a 2D plane, partitioned into n-sided, convex polygons. I'm using WRF's PNPOLY algorithm for polygon inclusion to ensure that a point belongs inside one and only one polygon.
Is there an algorithm I can use to clip a line segment PO to a given polygon in the plane, assuming that pnpoly(O) == true, such that pnpoly(P') will always be true?
My current implementation of clipToPoly does a line-line intersection test with the segment and each edge of the polygon, then uses the intersection point (as detailed in this SO answer), but this does not always yield a point that satisfies PNPOLY.
function clipPointToPoly(p, o, poly) {
var i, j, n = poly.length,
q, r = {}, s = {}, pq = {},
rxs, t, u;
function cross2(v, w) {
return v.x * w.y - v.y * w.x;
}
for (i = 0, j = n - 1; i < n; j = i++) {
q = poly[i];
s.x = poly[j].x - q.x;
s.y = poly[j].y - q.y;
r.x = o.x - p.x;
r.y = o.y - p.y;
rxs = cross2(r, s);
if (rxs !== 0) {
pq.x = q.x - p.x;
pq.y = q.y - p.y;
t = cross2(pq, s) / rxs;
u = cross2(pq, r) / rxs;
if (0 < u && u < 1 && 0 < t && t < 1) {
p.x = p.x + t * r.x;
p.y = p.y + t * r.y;
return true;
}
}
}
return false;
};
Here is my implementation of PNPOLY:
function pnpoly(p, poly) {
var i, j, e0, e1,
n = polygon.length,
inside = false;
for (i = 0, j = n - 1; i < n; j = i++) {
e0 = poly[i];
e1 = poly[j];
if ( ((p0.y > p.y) !== (p1.y > p.y)) &&
((p.x < (p1.x - p0.x) * (p.y - p0.y) / (p1.y - p0.y) + p0.x)) ) {
inside = !inside;
}
}
return inside;
};
I don't think I understand enough about the Simulation of Simplicity, or how to deal with PNPOLY's half-open sets when using floating point numbers to handle the edge cases properly.
For example:
poly: [(1,1), (-1,1), (-1,-1), (1,-1)]
p: (5,5)
o: (0,0)
p' = (1,1)
This fails because (1,1) is not included according to PNPOLY (it's on the open side of the set), but clipToPoly does not take that into account. I suppose I could nudge it by an epsilon if I knew it was on an open end of the set, but I'd prefer a more stable solution.
Another example:
poly: [-995.9592341908675, -88.48705014724577
-1040.5031753180106, -176.53192722405026
-549.9211095905894, -330.8462151682281
-653.7143990581328, -211.59193148034612]
p: -1032.3773586525654, -208.3586379393678
o: -957.4172402148379, -202.6668958854324
In this case, clipToPoly fails because O is so close to the edge of the polygon, it doesn't even detect an intersection due to floating point imprecision.
t: 1.0000000000000002 u: 0.8306380503739466
Is there a way to get clipToPoly's floating point imprecision to match PNPOLY's, so that both are consistent?
