Intersection may have badly ordered points.

[untoldwind]

const martinez = require('martinez-polygon-clipping');
const geojson2svg = require('geojson-to-svg'); 

const gj1 = {
    "type": "Feature", 
    "properties": {},
    "geometry": { "type": "Polygon", "coordinates": [[[ -530, -530], [-530, 530], [530, 530], [530, -530], [ -530, -530]]] }
}
const gj2 = {
    "type": "Feature", 
    "properties": {},
    "geometry": { "type": "Polygon", "coordinates": [[[1.2500125250252, -531],[-98, -531],[-98, 531],[1.250012525025, 531],[1.2500125250252, -531]]] }
}


const intersection = {
    "type": "Feature",
    "geometry": {
      "type": "Polygon",
      "coordinates": martinez.intersection(gj1.geometry.coordinates, gj2.geometry.coordinates)
    }
};

console.log(JSON.stringify(intersection, null, "  "));

console.log(geojson2svg()
  .data(gj1)
  .render());

  console.log(geojson2svg()
  .data(gj2)
  .render());

Will produce:

{
  "type": "Feature",
  "geometry": {
    "type": "Polygon",
    "coordinates": [
      [
        [
          [
            -98,
            -530
          ],
          [
            -98,
            -530
          ],
          [
            -98,
            530
          ],
          [
            1.2500125250249994,
            530
          ],
          [
            1.2500125250249994,
            530
          ]
        ]
      ]
    ]
  }
}

This seems to be a rounding-ish issue (ordering is correct when using 631 instead of 531)

[rowanwins]

Resolved in 0.6.0 by #111

Now produces

[-98, -530],
[1.250013, -530],
[1.250013, 530],
[-98, 530],
[-98, -530]

[rowanwins]

Not resolved by #111 , turns out the demo site had some issues :( Needs further investigation

[bluenote10]

I had a look into that and figured out what is happening. It is indeed a floating point stability problem. When exaggerating the coordinates a little bit, the example is structurally equivalent to this:

Note that the last edge of the smaller polygon is perfectly straight, but goes from a slightly smaller x on the top (A) to a slightly bigger x on the bottom (B). Now when computing the intersection points I1 and I2 the we get numerically unstable results. These are the involved x values:

x of A:  1.250012525025
x of I1: 1.2500125250249994 (actually slightly to the left of A)
x of I2: 1.2500125250252125 (actually slightly to the right of B)
x of B:  1.2500125250252

The problem here is the x value of I2 in combination with this optimization. In order to produce a proper result polygon we have to process the edges involving I2. However, the optimization breaks out of the while loop for any coordinate that falls outside the range of one of the polygons, because that area isn't relevant for processing the intersection in theory. So because I2.x > B.x it cannot be processed. The algorithm works correctly when disabled the optimization.

This particular example has a simple work-around: We could add a check after the intersection computation that ensures that intersection points don't fall out of the possible range. In the example, the valid x-values of an intersection points are [A.x, B.x]. If the intersection coordinates fall outside the range, we "reset" them to the border of the range. This works fine for this case, but I still have to think about it if it is the best solution. Maybe a robust intersection point computation is what we need in general.

Note: This test case is the first example that goes into this branch, which feels a bit like a hack to account for out-of-order intersection points, and currently doesn't seem to be covered by tests. It feels like we could get rid of it if the intersection points would never be out-of-order.

What would be an interesting test case: Similar scenario, but A.x and B.x are just one ULP (unit of least precision) apart from each other, and a bunch of intersection points between them. With the proposed work-around they would probably end up randomly on A.x or B.x -- not sure what would happen.

[rowanwins]

Good research @bluenote10 - I wonder if we could apply add an epsilon to the comparison and whether that would be enough to catch the edge cases...

[bluenote10]

@w8r / @rowanwins For the Rust version I'd propose this solution now: 21re/rust-geo-booleanop#11.

[rowanwins]

This is now resolved in 0.7.3