Polygon triangulation is broken for US states

[emackey]

I'll submit the US states GeoJSON on a branch. It's public domain, from Natural Earth.

[emackey]

Added the new file on branch us-states in 3e7a34d. Please use the same branch to develop the fix. Scott and I would like the new GeoJSON file added to master as an example.

[emackey]

To reproduce the bug, drag-and-drop the new file ne_10m_us_states.json onto Cesium Viewer.

message: "Tried 50 times to find a valid cut and couldn't."

[pjcozzi]

@hpinkos want to give this a shot? @DavidHudlow any chance you can advise?

[DavidHudlow]

Some of these polygons seem to be complex. I think that's where the problem is coming from.

[DavidHudlow]

The complex polygon I ran across was very close, geometrically, to being simple. Is it possible that a projection subtly pushed it over the edge? If necessary, I can add complex polygon handling to the triangulation algorithm. It already isolates the "complex" part of the polygon (usually 4-6 vertices) before failing.

[pjcozzi]

@DavidHudlow it's possible I suppose but we've been using the same tangent plane projection for years with a variety of use cases. It is the call to projectPointsOntoPlane from PolygonGeometry.js. I'm pretty sure I tested this same dataset with my C# implementation of this projection and our old ear clipping algorithm.

@emackey can you reproduce this with b19?

[DavidHudlow]

If you view that json file, there's a complex area on the southeast coast of Texas just south of Sweeny.

[emackey]

The file was working on the dataSourceBrowser branch up until a recent
merge from master. By "working" I mean displaying state outlines and not
crashing (no idea if the triangulation itself was correct).

On Fri, Sep 6, 2013 at 2:00 PM, David Hudlow notifications@gh.neting.ccwrote:

If you view that json file, there's a complex area on the southeast coast
of Texas just south of Sweeny.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/1121#issuecomment-23958109
.

[pjcozzi]

@DavidHudlow I think you are referring to here:

Zoomed out:

I agree that this is not a polygon we can - or should have to - triangulate. @emackey can you try tweaking the dataset? Could be fun to use http://geojson.io/

[shunter]

I generated the TopoJSON file in question using the tool at https://github.com/mbostock/world-atlas. Since complex polygons are presumably valid in TopoJSON/GeoJSON, I guess the question is what level should deal with this. Should we simplify polygons all the way up in the GeoJSON data source layer?

[DavidHudlow]

I suspect it'll be easier to handle this at the triangulation level. Detecting complex polygons requires a somewhat complex algorithm to accomplish in O(n log n) rather than O(n^2) time. Because the current triangulation algorithm already isolates the "complex" portion of the polygon, the analysis/simplification could be done at that point without much concern for performance.

[pjcozzi]

Here's another file with a similar problem: https://raw.github.com/turban/N2000-Kartdata/master/NO_Admin_latlng.topojson

Just drag and drop onto Cesium viewer.

[DavidHudlow]

So do y'all want me to take a stab at getting these to triangulate?

[emackey]

Yes please! 👍

I'm sure we could fix these files by hand, but I don't think we can claim full GeoJSON support without them.

[pjcozzi]

@DavidHudlow that would be awesome 👍

[DavidHudlow]

Just giving an update on this. I've made some progress on this, and I think I have a decent solution. Unfortunately, it does require adding intersection points to the polygon in question (which means some changes to PolygonGeometry).

The easiest way to implement this is using EllipsoidTangentPlane's projectPointsOntoEllipsoid method (using the same plane used to project to 2D). It seems, however, that this method provides increasingly erroneous results as the size of the ellipsoid increases from the Unit Sphere.

I added the following unit test to EllipsoidTangentPlane, and it fails. I haven't been able to track down why, and I'd love some support on that.

    it('projectPointsOntoEllipsoid works with an arbitrary ellipsoid using fromPoints', function () {
        var ellipsoid = new Ellipsoid(2,2,2); // Passes if unit circle, fails otherwise

        var points = ellipsoid.cartographicArrayToCartesianArray([
            Cartographic.fromDegrees(-72.0, 40.0),
            Cartographic.fromDegrees(-68.0, 35.0),
            Cartographic.fromDegrees(-75.0, 30.0),
            Cartographic.fromDegrees(-70.0, 30.0),
            Cartographic.fromDegrees(-68.0, 40.0)
        ]);

        var tangentPlane = EllipsoidTangentPlane.fromPoints(points, ellipsoid);
        var points2D = tangentPlane.projectPointsOntoPlane(points);
        var positionsBack = tangentPlane.projectPointsOntoEllipsoid(points2D);

        expect(positionsBack[0].x).toBeCloseTo(points[0].x);
        expect(positionsBack[0].y).toBeCloseTo(points[0].y);
        expect(positionsBack[0].z).toBeCloseTo(points[0].z);
    });

[pjcozzi]

@bagnell can you help @DavidHudlow with this?

[bagnell]

@DavidHudlow I'm looking into the problem, but there is a workaround you can use until its fixed. Here's the modified test:

    it('projectPointsOntoEllipsoid works with an arbitrary ellipsoid using fromPoints', function () {
        var ellipsoid = Ellipsoid.WGS84;
        var oneOverRadii = ellipsoid.getOneOverRadii();
        var radii = ellipsoid.getRadii();

        var points = ellipsoid.cartographicArrayToCartesianArray([
            Cartographic.fromDegrees(-72.0, 40.0),
            Cartographic.fromDegrees(-68.0, 35.0),
            Cartographic.fromDegrees(-75.0, 30.0),
            Cartographic.fromDegrees(-70.0, 30.0),
            Cartographic.fromDegrees(-68.0, 40.0)
        ]);

        var i, point;

        var scaledPoints = [];
        for (i = 0; i < points.length; ++i) {
            point = points[i];
            scaledPoints.push(Cartesian3.multiplyComponents(oneOverRadii, point));
        }

        var tangentPlane = EllipsoidTangentPlane.fromPoints(scaledPoints, Ellipsoid.UNIT_SPHERE);
        var points2D = tangentPlane.projectPointsOntoPlane(scaledPoints);
        var positionsBack = tangentPlane.projectPointsOntoEllipsoid(points2D);

        for (i = 0; i < points.length; ++i) {
            point = positionsBack[i];
            Cartesian3.multiplyComponents(radii, point, point);
        }

        expect(positionsBack[0].x).toBeCloseTo(points[0].x);
        expect(positionsBack[0].y).toBeCloseTo(points[0].y);
        expect(positionsBack[0].z).toBeCloseTo(points[0].z);
    });

[bagnell]

@DavidHudlow I fixed the tangent plane problem in the above pull request link.

[DavidHudlow]

Thanks, Dan! I'll get back to work on this then!
On Sep 13, 2013 5:15 PM, "Dan Bagnell" notifications@github.com wrote:

@DavidHudlow https://github.com/DavidHudlow I fixed the tangent plane
problem in the above pull request link.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/1121#issuecomment-24428359
.

[pjcozzi]

@DavidHudlow I merged @bagnell's changes into master. Any luck with this?

[DavidHudlow]

@pjcozzi Thanks, Patrick! I've been busy with some other stuff, but I'll take another look at this today.

[DavidHudlow]

Update: I feel like I'm almost there, but can't get the complex polygons to render. I have to update the 3D positions as well as come up with indices because I have to add the intersection points. I have PolygonGeometry passing unit tests building a complex polygon, but when I try to render the same polygon, it fails. (A zero Cartesian3 popped up and caused a dereferencing error when it was converted to a Cartographic).

I have my latest attempt here: https://github.com/DavidHudlow/cesium/commits/master

I overrode the Polygons Sandcastle example with the case I'm studying. I'll try continue on it tomorrow, but if anyone has some immediate insight, it would be appreciated.

[DavidHudlow]

I've got the algorithm working, but I've been able to find some cases that break it. I may not get back to working on this until next week. If anyone wants to pick up what I have and run with it, feel free.

https://github.com/DavidHudlow/cesium/commits/master

[pjcozzi]

Looks good now:

[pjcozzi]

http://gis.civilsolutions.biz/geojson/HammontonOnlyGeo2.topojson is also fixed.

[pjcozzi]

https://raw.github.com/turban/N2000-Kartdata/master/NO_Admin_latlng.topojson also works now. @mramato this is really slow now so it will be a nice test case for when you add geometry and appearances. It should get 10-100x faster.

[pjcozzi]

https://raw.github.com/turban/N2000-Kartdata/master/NO_Admin_latlng.topojson also works now.

@DavidHudlow actually this one is hit or miss. Sometimes it works:

Sometimes it doesn't:

[mramato]

90% of the geojson files I have fail with this bug, so I've marked it as a showstopper for next release.

[pjcozzi]

@DavidHudlow have you had any luck reproducing the issue? Or switching to a consistent random seed?

I agree with @mramato that this is a pretty significant problem for the b21 release on Tuesday. We have a few options.

1. Ship b21 with a full fix with consistent behavior between runs.
2. Ship b21 with the partial fix in Support Complex Polygons #1163.
3. Ship b21 without a fix since we already shipped this in b20.
4. Roll back to ear clipping.

Ideally, I'd like to do (1). I'm OK with (2) if we have to. (3) doesn't feel like a real option given how this issue comes up so often in practice and on the forum. (4) is probably more trouble than it's worth, and it moves us backwards, not forward.

@DavidHudlow any thoughts here?

[DavidHudlow]

@pjcozzi I've been covered up with other stuff, but let me take a shot at #1 today and let you know how it goes.

[pjcozzi]

@DavidHudlow thanks so much. This will be a much appreciated fix.

[DavidHudlow]

@pjcozzi Can't get fully complex polygon support in time for release, so how about option 5) Render non-complex portions of complex polygons.