JuliaGeo / GeometryOps.jl

GeoInterface-based geometry operations
https://juliageo.org/GeometryOps.jl/
MIT License
29 stars 4 forks source link

Predicates #4

Closed DanielVandH closed 10 months ago

DanielVandH commented 1 year ago

Following on here https://github.com/asinghvi17/GeometryOps.jl/pull/3

Using exact predicates will be something important to consider, making use of ExactPredicates.jl being the best choice. There is an issue of having to convert to Tuples, but this isn't an issue unless points are represented as (non-static) vectors - you already access individual coordinates anyway, so there is no difference in allocations. Also importantly, the performance loss of considering exact predicates is nothing major - for coordinates far from degeneracy, the fast path in ExactPredicates.jl makes it the same as using standard determinant definitions. For point sets where you end up needing to use the more complicated paths in ExactPredicates.jl, the performance does slow down - but that's exactly where you need it!

As an example, my original implementations in DelaunayTriangulation.jl just used determinant definitions for the predicates. You quickly run into infinite loops or e.g. triangles with the incorrect orientation. Some good examples demonstrating this importance here https://www.sciencedirect.com/science/article/pii/S0925772107000697 if you are curious. See also the first chapter here.

For the function is_point_on_segment I mentioned in the PR, where you need an exact predicate is in the computation of cross - this is just the definition of the orient predicate. It's really difficult to detect points on line segments without exact predicates, and when it matters it really matters that we detect collinear points in certain geometric algorithms (at least the ones I use anyway). You can see how I determine whether a point is on a line segment with ExactPredicates.jl here (which assumes already that the points are collinear, detected via point_position_relative_to_line (same as orient).

Lastly: You might be interested in considering three-valued logic rather than Bools, as argued e.g. here. I implement this with a Certificate system in DelaunayTriangulation.jl, but that's rather specific; ExactPredicates.jl uses (-1, 0, 1). It just turns out to be nice to be able to distinguish between cases like e.g. on, left, right. (Certificates were needed since some predicates even needed > 3 outcomes! like classifying the type of an intersection or specifying how a line intersects a triangle).

Anyway, don't have to take all of this on but thought you might appreciate some input. I had to think a lot about this at the time of initially implementing DelaunayTriangulation.jl.

rafaqz commented 1 year ago

This all sounds great, and like something we want. I think the allocations are actually fine, I didn't understand the scope of your package at first.

But the main problem is that - I don't really understand your package. The docs are terse and the code generation and macros means its hard to just drop down to reading the code.

We will likely need someone with your expertise to actually hook this up.

For now at least I will just bang out the structure of the package, as far as I can tell we can plug in better predicates at any time?

DanielVandH commented 1 year ago

(I've nothing to do with the code in that package, FYI - My contributions there are just some typo fixes. It's just a very useful tool / impressive piece of work.)

I'm not sure the main difficulties with the package's code/documentation. The start of the docs here https://danielvandh.github.io/DelaunayTriangulation.jl/dev/predicates/, after Certificate and up to the meet_predicate docstring, might help clarify some things, the functions in that part are all from ExactPredicates.jl. The aim is to, essentially, compute the value of the predicate and use that value to determine whether more precision is needed, adapting until the sign of the predicate is known exactly. (Codegen is very important for this unfortunately, hence the need for so much of it in the code.) Anyway, yes I could probably look into it myself eventually if needed.

Indeed, you can just hook into it anytime. That's something I considered here https://github.com/DanielVandH/DelaunayTriangulation.jl/blob/main/src/predicates/general.jl#L1-L110 also, just redefining orient_predicate = orient and incircle_predicate = incircle, etc., so that the definitions can be easily changed later if needed.

rafaqz commented 10 months ago

@DanielVandH we now use ExactPredicates.jl :)