Add Bentley-Ottman Algorithm #1168
Conversation
…uts a dictionary of vertex keys with segment values
…ntersections and ultimately calculates all intersections
|
@souma4 please let me know when the PR is ready for review. I have some review comments already, just waiting for your green light. |
|
@juliohm Thanks for reminding me! |
There was a problem hiding this comment.
Thank you @souma4 for working on this. It is really nice to see the amazing progress you've made already ❤️
Attached is my first round of review with a couple of suggestions to improve the PR. Let's work on it together slowly, possibly improving the BinaryTrees.jl dependency before coming back here.
Looking forward to the next round. It looks very promising!
…tests. Edited test to be more intensive. updated entire script to meet style and Julia best practices. Removed unneded functions. Shifted BinaryTrees relevant materials to PR JuliaGeometry#12 in BinaryTrees and adjusted code as needed. Reordered point processing to reduce the likelihood of a bug when processing start and intersection segments.
…d function. removed println calls used for debugging
|
Tests are failing because the utility function is not exported. We can either export it or qualify the calls in the test suite with the |
Co-authored-by: Júlio Hoffimann <[email protected]>
…eaner enum handling, and shifted the final dictionary to be unique values
…hes.jl into JRC_add_BentleyOttman
… that connect to the point
…status! (mostly has to do with the isless definition) but this works in all degenerate cases in FP precision.
…ty and reducing redundant computations
… needs more time or more brains. Also may have fixed FP32 robustness. trying to move away from snapping to grid until pushing to the output (maintaining true spatial relationships)
|
Thanks for pushing this @souma4 ❤️ I will try to review by the end of the week 🙏🏽 |
|
I'm honestly glad for the extra time, managed to fix some huge uneeded time sinks. I know there's still optimization there (just have to wrap my head around it) and FP32 still failing is aggravating. I need to research more robust checks. I saw one paper suggesting signed area is "second order" which makes FP64 work but not 32. A solution might be buried in there without requiring a whole new data structure. But you'll probably have some ideas to resolve type instabilities in the intersection calc. |
…orientation tests. also fixed tests for FP32. When the values of FP32 are correct, the algorithm works. Rotation isn't 100% FP32 robust and the intersection calculations for the tests at L127 don't end up near the correct points, thus rounding isn't working.
Benchmark Results (Julia vlts)Time benchmarks
Memory benchmarks
|
Benchmark Results (Julia v1)Time benchmarks
Memory benchmarks
|
There was a problem hiding this comment.
I know this PR is taking longer than what we expected, but I think the only way to speed things up is to make sure we have a version that is ready for review. My understanding is that you are still exploring ideas and variations of the code, and that it is not polished for final review yet. Could you please try to clean the implementation as much as possible, follow our code style and variable names, and then add comments explaining what is happening?
If you feel that Float32 is an unattainable goal, we can focus on a working version that is robust with Float64, and only consider tests in that context. The Float32 case could be useful in the context of GPUs, but we need a bullet-proof CPU implementation first.
…wrappers to keep things simpler. Updated test documentation to be more clear about intent
src/Meshes.jl
Outdated
| @@ -29,6 +30,7 @@ using DelaunayTriangulation: get_polygon_points | |||
| using StatsBase: AbstractWeights, Weights, quantile | |||
| using TiledIteration: TileIterator | |||
| using ScopedValues: ScopedValue | |||
| using AdaptivePredicates: orient2, orient2fast | |||
There was a problem hiding this comment.
orient2fast is not used in the PR
Also, could you please confirm that Meshes.signarea was not enough and that we really need orient2 to make this work?
There was a problem hiding this comment.
Thank you for asking me to confirm because signarea actually was sufficient. I added AdaptivePredicates.jl as a way to have maximum robustness and to somewhat get around FP32 lack of precision. Plus, the overhead for it was minimal. But signarea is giving the same passing tests. I've just left a breadcrumb comment in case we ever need to revisit this.
test/utils.jl
Outdated
| segs = Segment.([(cart(0, 0), cart(2, 2)), (cart(0, 2), cart(2, 0)), (cart(0, 1), cart(0.5, 1))]) | ||
| points, seginds = Meshes.bentleyottmann(segs) | ||
| @test length(points) == 7 | ||
| @test length(seginds) == 7 |
There was a problem hiding this comment.
I believe we discussed this issue before. The algorithm should only return the intersection points. Right now it is returning all points, leading to unnecessary allocations.
…e input data is ordered in a useful way. updated documentation to note FP32 concerns. split start and end initialization into one function that can handle either (maintainability)
|
@souma4 I tried to cleanup the code further, but I am finding it difficult to do so. I believe you are used to other programming languages where types appear in many places. Can you please try to reduce the amount of explicit type handling? Many functions like Also, try to follow our code style with comments in lowercase, avoid unnecessary underscores (e.g., I still believe that we can simplify this implementation further to a point where we can clearly read it and know what is happening without too much back and forth with low-level information passed around. For example, sometimes you define auxiliary functions in terms of I really think we have something promising here, but we will only reach a point where we feel safe about the implementation when these iterations and reviews start to converge into a cohesive, consistent code. Thanks for the amazing effort you already put into this 🙏🏽 ❤️ |
…roughout function barriers to easily track what kind of object is pass to which function. Fixed the tests to account for ordering of outputs because hashtables change upon new sessions, leading to variable ordering of outputs. This simplifies tests.
| seg = first(segments) | ||
| ℒ = lentype(seg) | ||
| τ = ustrip(eps(ℒ)) | ||
| round(Int, 0.8 * (-log10(τ))) # 0.8 is a heuristic to avoid numerical issues |
There was a problem hiding this comment.
This 0.8 was not here before. Can we simply use tau = atol(numtype(lentype(seg)))?
There was a problem hiding this comment.
That does work, however it does seem to make one more Float32 test fail (easy to skip it). The 0.8 arose because I had switched to using machine epsilon. If I remember right that was largely to more specifically reflect the precision of the numtype. But this atol switch is a bit more forgiving and seems to work. I'm cool with tau = atol(numtype(lentype(seg)) rather than being reliant on a tuned parameter from eps.
There was a problem hiding this comment.
I don't remember the code that was there before, but the 0.8 factor is very ad-hoc. We shouldn't have that.
| _keyseg(ns) = _segment(BinaryTrees.key(ns)) | ||
|
|
||
| # handles the degenerate case of trivial (0-length) segments | ||
| _istrivial(seg) = seg[1] == seg[2] |
There was a problem hiding this comment.
I am concerned that I am still reviewing dead code. This function isn't used anywhere.
see issue #1126 and PR #1125
I needed to refactor and I messed up with the last draft #1143.
This produces a dictionary of vertices with a vector of segment values. Future needs: A test that the outputs are correct, and general method for rebuilding polygons with this structure.
@juliohm sorry about messing up the prior draft. Live and learn. I'll let you look this over. I think a polygon build method for this output is a new PR. This outputs vertices with segment info. If someone just wants the vertices they can grab the keys. If someone wants a polygon they can use a build method.