fix(boolean-contains): reject inner polygon whose edge crosses a conc…

[aloktomarr]

Fixes #2242.

Problem

booleanContains(polygon, polygon) only checked that every vertex of the inner polygon lies inside the outer polygon. When the outer polygon is concave, an inner edge can pass through a notch (the exterior) while all of its vertices remain inside — so a polygon that is not actually contained was reported as contained.

Fix

In isPolyInPoly, after the existing all-vertices-inside check, split each edge of the inner ring on the outer polygon's boundary and reject containment if any resulting sub-segment's midpoint lies strictly in the outer polygon's exterior. This mirrors the midpoint-based approach already used by lineInPolyStatus. Edges that run along the outer boundary are still allowed (checked via a small-epsilon boundary test).

Tests

Added test/false/PolygonIsNotContainedByConcavePolygon.geojson: a concave outer polygon (a square with a V-notch) and an inner rectangle whose vertices are all inside but whose top edge crosses the notch. It returns false after the fix; the full @turf/boolean-contains fixture suite still passes.

[mfedderly]

Thanks for taking this on! I think the same problem that applies to the vertices can also apply to the midpoints though. So you'll still have false positives, just fewer than before adding the midpoint check. I think the final step in this check probably has to be a more complicated version of intersection checking using something like a sweepline.

The boolean-* methods are intended to be faster versions because they aren't required to find all intersections, but I'm not sure we have a sweeplines implementation yet that allows us to exit on the first intersection that allows us to return false. @turf/line-intersect might at least give us a correct implementation to start with? We could add an option to return after the first intersecting point was discovered to retain as much speed as possible.

[aloktomarr]

Thanks for the quick look!

It's actually not checking a single midpoint per edge. It splits each of feature2's edges at
the points where they cross feature1's boundary (using @turf/line-split, which is lineIntersect
under the hood) and then checks the midpoint of each resulting sub-segment. A sub-segment has no
crossing inside it, so its midpoint is a fair stand-in for the whole piece, and that's what lets
it catch the concave notch in the fixture (the piece sitting in the notch has its midpoint outside).

That said, you're right that it isn't complete. It leans on line-split catching every crossing,
so the awkward cases still get through: mainly an edge that crosses the boundary right at one of
feature1's vertices (where the intersection can get missed or deduped), plus tangents and
collinear overlaps. In those the midpoint can land inside even when part of the segment is
outside, so you'd still get a false positive. So it's a solid improvement over what's there today
(it fixes the reported case), just not the whole story.

Happy to take it further either way. I can leave this as the incremental fix and open a follow-up
for the rigorous version, or do the rigorous one right here: collect all the inner-edge/outer-boundary
intersections with @turf/line-intersect so the vertex hits are handled explicitly, split on those,
and bail out with false on the first sub-segment whose midpoint is outside, so the boolean-* fast
path stays fast. Do you have a preference?

[mfedderly]

Ah sorry I missed that yesterday. I was bouncing around on a bunch of different tickets and didn't fully understand what you meant at first. This implementation is basically what I had in my head after a much closer review. The only thing that makes me a little uncomfortable is the hardcoded epsilon in there.

In a look around I noticed that de9im has a pretty different implementation involving triangulation. I didn't reason about its correctness, but it may be speedy. JTS (or JSTS) might also have a reference implementation we can look at.
https://github.com/dpmcmlxxvi/de9im/blob/master/src/util/polygon.js#L68-L69
https://github.com/dpmcmlxxvi/de9im/blob/master/src/util/polygon.js#L100-L138

Another thing that might work is just doing a full polygon/polygon difference operation, but polyclip-ts that powers @turf/difference at the moment is known to be quite slow (it uses BigDecimal under the hood). I have geoclipper2 ready to handle buffer, union, difference, etc but we unfortunately have to hold off on using it until we do a v8 (it introduces BigInt usage for the first time, and the buffer operation needs an API break for some parameter changes). I'm currently trying to work through all of the open PRs for a v7.4 before making the major break though. I don't actually know if using difference would be faster than this implementation.

A few thoughts on the 'boolean-* as faster' thing. This can be really hard to achieve if you're still maintaining correctness. At some point the work required to return the correct answer is just really expensive (like for polygon pairs here). I'm pretty sure its possible if you're doing the difference operation to interrupt the sweepline process as soon as you know there's something, but I don't know that I've actually seen an implementation that had that out of the box. I think its probably possible to thread something into geoclipper2 that would do it though.

So I think there's a few paths forward:

1. merge this PR, maybe with some thought around the epsilon hardcoding
2. similar PR, but it just uses @turf/difference here and we accept that its really slow today, but faster in v8
3. we like the de9im version and choose that instead of any of this

Thoughts?

[aloktomarr]

Thanks for taking another look, and no worries at all. Glad the approach lines up with what you had in mind after the closer review.

On the epsilon

You're right to be uneasy about it. It gets passed straight into turf's existing booleanPointOnLine epsilon option, where it's compared against the collinearity cross product |(pt − a) × (b − a)|. That quantity scales with edge length and coordinate magnitude, so a fixed 1e-9 isn't scale invariant. It's fine for lng/lat, but for large or projected coordinates it gets far too strict and starts rejecting valid containment. So yes, the magic number is a real wart.

I'd like to keep this approach (path 1) but swap the constant for a scale relative tolerance derived from the geometry, so it behaves the same no matter what coordinate units we're in. Two ways to do that: normalize the cross product by the edge length, which turns it into an actual distance tolerance, or compute point to segment distance directly and compare it to REL * scale. Either way the hardcoded value goes away. If you'd rather keep the change minimal, I'm just as happy to lift it into a single documented constant, or to run the fixture suite with the tolerance removed so we can see exactly which fixtures lean on it.

One thing that holds regardless of which we pick: the tolerance only ever widens the "this midpoint lies on the shared boundary" acceptance. A real exterior crossing sits strictly outside and fails both the point in polygon test and the bounded booleanPointOnLine check, so no small tolerance can turn a real exterior edge into a false accept. It only stops floating point false negatives on boundary coincident edges.

On @turf/difference (path 2)

Agreed this is the cleanest long term formulation, but like you said it's slow on polyclip-ts today and the geoclipper2 swap is gated on v8. Also worth flagging that without an interruptible sweepline, a full difference is strictly more work than this PR for the common "not contained" case. This implementation already bails early: it returns false on the first subsegment whose midpoint lands in the exterior, with no global operation. So it feels like a natural v8 follow up once geoclipper2 can early exit, rather than something to change now.

On de9im (path 3)

I went and read the lines you referenced. isInPolygon (lines 62 to 70) starts with the same "all of #1's vertices inside #2" check we do, then hands off to relate (lines 114 to 138), which triangulates both polygons, partitions #1 against #2 (lines 129 to 131) so each triangle is wholly in or out, and tests each triangle's centroid (lines 134 to 136). It's basically this PR's strategy one dimension up: triangle centroid instead of edge subsegment midpoint. Two things jump out though. It triangulates and partitions both polygons upfront on every call with no early exit, and it doesn't actually get rid of the float sensitivity. relate has to special case identical polygons (lines 119 to 124) because the clipping tools it uses for partitioning break on identical high precision triangulations. So it swaps one hardcoded epsilon for a dependency on triangulation plus the clipping layer, plus its own structural workaround. I don't think it's clearly more correct or faster here, but I'm happy to benchmark it against this once the v8 clipping story lands. (JTS/JSTS would be the authoritative reference, but that's RelateOp and PreparedPolygon.contains via robust DE-9IM edge intersection rather than triangulation, and porting that is well beyond the scope of this fix.)

My vote

Merge path 1 with the scale relative epsilon, and treat difference based containment as a v8 follow up once geoclipper2 supports an interruptible sweepline. Happy to push the epsilon change to this branch now if you're good with that.