mapbox / wagyu

A general library for geometry operations of union, intersections, difference, and xor
Other
166 stars 33 forks source link

Drops small polygons #94

Closed Algunenano closed 5 years ago

Algunenano commented 5 years ago

I have a geometry that is producing invalid outputs when intersecting it with others or itself.

Here is an example of an intersection that works. Since I'm intersecting a polygon (POLYGON((0 0, 0 10, 10 10, 10 0, 0 0)) with itself, I'm expecting the result to be the same polygon):

Click to expand

```c++ TEST_CASE("Self intersection (OK)") { mapbox::geometry::linear_ring ring0_0; ring0_0.push_back(mapbox::geometry::point( 0, 0)); ring0_0.push_back(mapbox::geometry::point( 0, 10)); ring0_0.push_back(mapbox::geometry::point(10, 10)); ring0_0.push_back(mapbox::geometry::point(10, 0)); ring0_0.push_back(mapbox::geometry::point( 0, 0)); mapbox::geometry::polygon pol0; pol0.push_back(ring0_0); mapbox::geometry::linear_ring ring1_0; ring1_0.push_back(mapbox::geometry::point( 0, 0)); ring1_0.push_back(mapbox::geometry::point( 0, 10)); ring1_0.push_back(mapbox::geometry::point(10, 10)); ring1_0.push_back(mapbox::geometry::point(10, 0)); ring1_0.push_back(mapbox::geometry::point( 0, 0)); mapbox::geometry::polygon pol1; pol1.push_back(ring1_0); mapbox::geometry::wagyu::wagyu wagyu; wagyu.add_polygon(pol0, mapbox::geometry::wagyu::polygon_type::polygon_type_subject); wagyu.add_polygon(pol1, mapbox::geometry::wagyu::polygon_type::polygon_type_clip); mapbox::geometry::multi_polygon solution; wagyu.execute(mapbox::geometry::wagyu::clip_type_intersection, solution, mapbox::geometry::wagyu::fill_type_even_odd, mapbox::geometry::wagyu::fill_type_even_odd); REQUIRE(solution.size() == 1); REQUIRE(solution[0].size() == 1); REQUIRE(solution[0][0].size() == 5); } ``` ```bash $ ./cmake-build/unit-tests "Self intersection (OK)" =============================================================================== All tests passed (3 assertions in 1 test case) ```

On the other hand, when I change the geometry to this one: `POLYGON((1252 1904,1253 1905,1253 1906,1251 1904,1252 1904))` what I don't get anything out of the intersection with itself:
Click to expand

```c++ TEST_CASE("Self intersection (KO)") { mapbox::geometry::linear_ring ring0_0; ring0_0.push_back(mapbox::geometry::point(1252, 1904)); ring0_0.push_back(mapbox::geometry::point(1253, 1905)); ring0_0.push_back(mapbox::geometry::point(1253, 1906)); ring0_0.push_back(mapbox::geometry::point(1251, 1904)); ring0_0.push_back(mapbox::geometry::point(1252, 1904)); mapbox::geometry::polygon pol0; pol0.push_back(ring0_0); mapbox::geometry::linear_ring ring1_0; ring1_0.push_back(mapbox::geometry::point(1252, 1904)); ring1_0.push_back(mapbox::geometry::point(1253, 1905)); ring1_0.push_back(mapbox::geometry::point(1253, 1906)); ring1_0.push_back(mapbox::geometry::point(1251, 1904)); ring1_0.push_back(mapbox::geometry::point(1252, 1904)); mapbox::geometry::polygon pol1; pol1.push_back(ring1_0); mapbox::geometry::wagyu::wagyu wagyu; wagyu.add_polygon(pol0, mapbox::geometry::wagyu::polygon_type::polygon_type_subject); wagyu.add_polygon(pol1, mapbox::geometry::wagyu::polygon_type::polygon_type_clip); mapbox::geometry::multi_polygon solution; wagyu.execute(mapbox::geometry::wagyu::clip_type_intersection, solution, mapbox::geometry::wagyu::fill_type_even_odd, mapbox::geometry::wagyu::fill_type_even_odd); REQUIRE(solution.size() == 1); REQUIRE(solution[0].size() == 1); REQUIRE(solution[0][0].size() == 5); } ``` ```bash ./cmake-build/unit-tests "Self intersection (KO)" ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ unit-tests is a Catch v1.9.6 host application. Run with -? for options ------------------------------------------------------------------------------- Self intersection (KO) ------------------------------------------------------------------------------- /home/raul/dev/public/wagyu/tests/unit/vatti.cpp:443 ............................................................................... /home/raul/dev/public/wagyu/tests/unit/vatti.cpp:475: FAILED: REQUIRE( solution.size() == 1 ) with expansion: 0 == 1 =============================================================================== test cases: 1 | 1 failed assertions: 1 | 1 failed ```

In both cases the geometries are identical, but in one the intersection works as expected while in the other one it doesn't. In the second case, if I remove `correct_topology` the resulting polygon has the following points: ``` 1252 1904 1253 1905 1253 1906 1253 1905 1252 1904 1251 1904 1252 1904 ``` So at some point, the algorithm has decided to turn back around (1253 1905 => 1253 1906 => 1253 1905) instead of continuing its path. A similar thing happens if I intersect the bad polygon with its bounding box or with a bigger box (in my case an standard mvt tile (0 0, 4096 4096). Is the code to calculate the intersection correct? I'm investigating what's causing the issue so any pointers are greatly welcomed.
flippmoke commented 5 years ago

@Algunenano sorry for the late responses here as I have been on holiday. I am going to dig into this right now

flippmoke commented 5 years ago

Alright, I do not think this is a bug -- this has to do with the way that snap rounding occurs and the danger of using wagyu with very small differences between your integers. The snap rounding within Wagyu attempts to snap any possible crossing of a line within a 0.5x0.5 square of any point being processed. After iterating through it snaps down this entire shape into a single line, which is then discarded during processing.

Here is a screenshot of the original shape:

screen shot 2019-01-02 at 3 37 41 pm

Here is a screenshot of the result after the vatti step (before it is then passed to correct_topology)

screen shot 2019-01-02 at 3 41 22 pm
flippmoke commented 5 years ago

I do not like that the algorithm does this, but it is one of the few ways that it can fix really nasty issues with data. The solution might be to consider how we could scale up and then down the values here based on a parameter before entering and exiting the algorithm. Something like "snap_rounding_precision", where if you set it to a value of 0.5 for example, it would scale all values up and then back down before performing the entire operation.

However, in this case the results would then need to be in a floating point format as a result, because you could have points at 0.5 etc.

Algunenano commented 5 years ago

@Algunenano sorry for the late responses here as I have been on holiday. I am going to dig into this right now

No worries; happy new year!

Alright, I do not think this is a bug -- this has to do with the way that snap rounding occurs and the danger of using wagyu with very small differences between your integers. The snap rounding within Wagyu attempts to snap any possible crossing of a line within a 0.5x0.5 square of any point being processed. After iterating through it snaps down this entire shape into a single line, which is then discarded during processing.

Are you sure that that snapping is the issue? If I scale the polygon (I've tried both x2 and x10) I get nothing out of the Vatti execution, and if I translate the polygon (x += 10000, y stays the same) I also don't get anything out of that step.

However, in this case the results would then need to be in a floating point format as a result, because you could have points at 0.5 etc.

I guess that even if this worked I'd need to be able to do the scale down before validation, as rounding to int values might produce invalid polygons again.

Algunenano commented 5 years ago

Are you sure that that snapping is the issue? If I scale the polygon (I've tried both x2 and x10) I get nothing out of the Vatti execution, and if I translate the polygon (x += 10000, y stays the same) I also don't get anything out of that step.

There was an error in these tests, this indeed works but it has other effects (see above).

So I've tried 2 things to workaround this small polygons disappearing:

In both cases I end up seeing more small parts that before were removed (like the ones in the issue I opened) but I also see self-intersections and spikes; for example: Before: image

After: image

Obviously it'd be great to have the cake and eat it, but I don't see an easy way to keep the small polygons and the validity. Taking into account that dropping small polygons in MVT isn't that bad (it's just that we are dropping some extra ;D), I'll stick with having the cake and maintain validity.

Closing this for now, maybe I'll revisit it in the future.