Inexact calculations effects on clippings

[skygering]

Here is an example of inexact calculations causing issues:

scale_fac = 1e5
c1 = [[-28083.868447876892, -58059.13401805979], [-9833.052704767595, -48001.726711609794], [-16111.439295815226, -2.856614689791036e-11], [-76085.95770326033, -2.856614689791036e-11], [-28083.868447876892, -58059.13401805979]]
c1 = [c ./ scale_fac for c in c1]
c2 = [[-53333.333333333336, 0.0], [0.0, 0.0], [0.0, -80000.0], [-60000.0, -80000.0], [-53333.333333333336, 0.0]]
c2 = [c ./ scale_fac for c in c2]
p1 = GI.Polygon([c1])
p2 = GI.Polygon([c2])
GO.intersection(p1, p2; target = GI.PolygonTrait(), fix_multipoly = nothing)

When scale_fac = 1 this doesn't work, but when scale_fac = 1e5 it does work...

Basically, there are two edges over one another on top of the x-axis, but one of the lines is slightly below the x-axis. I am working on a "patch" for this type of issue, but we might want a stronger fix at some point.

[asinghvi17]

We could always add ExactPredicates as a "backend", and use regular floating point computations as the default? Or the other way around. It looks like the refactor from earlier should allow for that if I understand correctly?

[skygering]

It is the math here though... not comparisons like ExactPredicates would be used for. Maybe we can do an exact comparison if the numbers are really close to zero... Or somehow use the backed package ? I don't really understand how ExactPredicates works yet...

Also, how do we know when to use it? It shouldn't be used all of the time... it is much to expensive and I was just running my simulation tests and I am getting these errors. It isn't some crazy edge case for me at least.

What do we think about adding an atol parameter that the user can pass in that allows for a bit of fudging to help with this problem for now? This is still problematic and should be fixed long term, but only happens when the polygons are of a large scale and then things are super close to zero... So for example if my polygons above are in meters, I am okay with fudging of up to like 0.5 meters personally for my applications, which would get rid of this problem I think.

I am drafting a PR with this to give the user more flexibility.

[asinghvi17]

Would these issues have happened before #72? More specifically, is this what you're looking for?

I got that result without scaling on the commit immediately before #72, so if that's the case then it shows that the exact arithmetic was probably helping here...

[skygering]

Yeah that should be it. Ugh! I just don't think we should be calling it all the time. The performance is SO bad on the geometries in #72.

But this is an example where it would actually help... it is just that Exact Predicates doesn't actually do the math. It just tells us a bit about what type of intersection it is. That doesn't actually help fix everything I don't think... but it would help haha

[asinghvi17]

ExactPredicates probably also needs some help becoming faster (though there's only so much type stability you can have in that case :D). Maybe playing with turning down the threshold for where ExactPredicates falls into adaptive math might help?

In general this is why I was advocating for a FastMath()/ExactMath() (or maybe even ; exact = _True()) - though I guess you would need to have it on all the time unless you are 100% sure about how your polygons are conditioned...

[skygering]

Yeah... It is hard to know if you are going to need it until you need it. And you only know because the polygons come out weird.

[rafaqz]

I think the majority of our exact predicates calls would be removed with judicious extents checks, and sorts or other preparations .

We should try to optimise everything else hard so its less painful to turn it back on.

[asinghvi17]

Yeah, preparation should definitely speed this up a lot if LibGEOS is any indicator!

[skygering]

Another example:

p1 = GI.Polygon([[[-57725.80869813739, -52709.704377648755], [-53333.333333333336, 0.0], [-41878.01362848005, 0.0], [-36022.23699059147, -43787.61366192682], [-48268.44121252392, -52521.18593721105], [-57725.80869813739, -52709.704377648755]]])
p2 = GI.Polygon([[[-60000.0, 80000.0], [0.0, 80000.0], [0.0, 0.0], [-53333.33333333333, 0.0], [-50000.0, 40000.0], [-60000.0, 80000.0]]])

[skygering]

I also did find another bug that maybe this one sets off... but even after fixing that bug this one still doesn't work. So if it errors rather than just breaks, let me know.

[asinghvi17]

Are those polygons meant to be unioned or intersected?

[skygering]

Intersected. Basically, the problem is that the left-most intersection point is right right next to the orange one's endpoint... by like 1e-11.

[asinghvi17]

LibGEOS just returns a linestring for this second scenario - LINESTRING (-53333.33333333333 0, -41878.01362848005 0)

[asinghvi17]

A more faithful overview of different predicate schemes:

ExactPredicates is actually >1000x slower, not 50x slower! Adaptive methods are 10x slower but that itself should speed up polygon operations significantly over ExactPredicates