Add Bentley-Ottman Algorithm

[souma4]

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.

[juliohm]

@souma4 please let me know when the PR is ready for review. I have some review comments already, just waiting for your green light.

[souma4]

@juliohm Thanks for reminding me!

[juliohm]

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!

[juliohm]

Attaching a few minor comments. Tomorrow I will try to find some time to run tests locally.

The good news is that the output of @code_warntype is all green.

[juliohm]

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 Meshes prefix.

[juliohm]

More basic fixes before the actual review of the algorithm.

[juliohm]

Tests are failing due to the latest review changes.

[juliohm]

I introduced a helper _digits function to extract the digits based on the atol of the underlying floating point type. Let's see if tests pass. If they still fail in Float32, we can relax the heuristic a bit more.

[juliohm]

The heuristic for the number of digits worked well. Assuming a tolerance of 1e-10, we get 10-1 digits by default. The same for 1e-5, which becomes 5-1 digits.

I noticed that the output of the algorithm is including the original points besides the intersections. For instance, here is what I get in the second test:

viz(segs)
viz!(points, color="red")

Could you please take a look to make sure we are only returning intersection points?

[souma4]

Oh if we only want intersections (not all points) then this is an easy change. I'll commit fixes and update tests to reflect. Although I am finding some variance in tests, so I'll need to debug and see if it's real or not.

Outside of that, the facetgrid test returns nothing because the facet generator only generates cornerintersections, which are start/endpoints and not included.

[juliohm]

Outside of that, the facetgrid test returns nothing because the facet generator only generates cornerintersections, which are start/endpoints and not included.

Very good point. I think we need to refactor this test by constructing segments manually in a grid-like arrangement. We could do the following instead:

horizontal = [Segment((1, i), (n, i)) for i in 1:n]
vertical = [Segment((i, 1), (i, n)) for i in 1:n]
segs = [horizontal; vertical]

[juliohm]

Thank you @souma4 for looking into it ❤️ Tests with random segments show that the algorithm is not there yet. Perhaps we could try to follow the algorithmic structure in Chapter 2 of Berg et al. Computational Geometry - Applications and Algorithms instead. This structure is discussed step-by-step in the video I linked in the issue. More specifically, it is discussed in this second video: https://www.youtube.com/watch?v=qkhUNzCGDt0

(If you need access to the book, I have the PDF and can send it to you directly)

We are almost there!

[souma4]

agreed. My plan was to pull up De Berg to see exactly what they do. I'm fairly busy this week, but, it'll keep moving haha. I think we are pretty close though.

[juliohm]

Did you find an issue with our built-in intersection method? It is not clear to me yet why we need to handle coordinates manually within a custom intersection algorithm. From what I recall, the Bentley-Ottmann algorithm assumes an atom function intersect(s1, s2) with the usual definition.

[souma4]

I did. The only consistent way of calculating the x intersect (which now that I think about this again, maybe x intersect and all derivates should be y intersect, because x=# intersect segment can just be seen as a "y intercept", but my intial logic was "this is the x at the intersect"). Is to either do this cross product solution or to set a Segment((x,-Inf),(x,Inf)) intersection, but in the intersection script, it uses a method that can't compute the solution due to the infinities. Another potential solution is to use a y extent determined in the initialization and send that to find intersection. Right now I'm finding the current intersect methods struggle with the curved lines, so I felt fine doing the linear thing.

[juliohm]

Right now I'm finding the current intersect methods struggle with the curved lines, so I felt fine doing the linear thing.

Do you mean our intersect atom has bugs that we need to tackle first?

Is to either do this cross product solution or to set a Segment((x,-Inf),(x,Inf)) intersection, but in the intersection script, it uses a method that can't compute the solution due to the infinities.

Handing coordinates at infinity is indeed tricky. I wonder what the code would look like with the alternative solution

Another potential solution is to use a y extent determined in the initialization and send that to find intersection.

I understand it would just call our built-in intersect with the original segments and 4 additional segments beyond the boundingbox. For example, segs |> boundingbox |> boundary |> segments |> collect. We could also add an extra room to the bounding box with Stretch(factor).

[juliohm]

@souma4 thank you very much for your resilience! I wish Bentley-Ottmann had shared a precise pseudo-code in their original article. That would have saved us a huge amount of time handing these tricky corner cases.

Could you please remind what was the issue with that original vanilla implementation? Was it purely caused by floating point arithmetic or there was something else?

[souma4]

The issue with the original is it struggled with verticals and multiple intersections. Floating point was mainly an implementation detail that didn't transfer over nicely from the good math. But hey, it works, it seems to have good time complexity, woohoo.

[juliohm]

The issue with the original is it struggled with verticals and multiple intersections. Floating point was mainly an implementation detail that didn't transfer over nicely from the good math. But hey, it works, it seems to have good time complexity, woohoo.

Let's investigate these two issues separately: verticals and multiple intersections.

Regarding the verticals, what exactly is the issue? I remember the solution in that video involved a small rotation of the sweep line, which was implemented implicitly in the order relation when we compare "x" and then "y" coordinate in order (lexicographical comparison). Our isless operation used in the binary tree should already take care of that?

Regarding the multiple intersections, I believe this is a more serious issue. The algorithm presented in the video (and in the book) seemed ok with multiple intersections. The algorithm implemented in the PR may be doing something slightly different, even though you adapted it from the original paper by Bentley-Ottmann.

[souma4]

Right now, the algorithm in this PR follows the version from the video and de Berg et al.—essentially modified Bentley-Ottmann.

Handling Multiple Intersections

I'll start with handling multiple segments intersecting at a point, since that ties into vertical segment handling later. When several segments meet at a point p, we need to:

1. Find which segments intersect at p (_findintersections)

2. Remove all those segments from the status.

3. Reinsert them in order from bottom to top based on where they intersect the sweepline (I used _nudge to make the order correct)

4. Identify the lowest and highest segments that pass through p.

5. Check for new intersections just below the lowest and above the highest segments.

To do this ordering properly, I compute the xintersection in the _SweepSegment struct. That allows comparison using a custom is.less that stores the segment based on the y-coordinate where each segment intersects the sweep line. This makes the status structure behave correctly—it’s ordered by y-position at the current sweep x.

I tried other methods (hardcoding order, swapping axes, using midpoints, holding segments as values instead of keys), but none handled edge cases well. Luckily, computing the intersection is cheap. Once calculated, comparisons are simple: go right if you're below the segment, left if you're above.

Status Structure and Staleness

The status tree is ordered by these computed xintersect values. But I only recompute _SweepSegments when they're involved in an intersection. That means the structure can get stale: the relative order is still valid, but the actual intersection points might not be.

So, when checking "is this segment active?" we need to recalculate the intersection to be sure. That’s what _findintersections does—it performs a guided search using the same relative order logic as before. It stays O(log n), unlike previous O(n) brute-force checks.

I considered trying to do what I'll call a "Meta Tree" (there might be a real name for it) where all nodes are recomputed upon a change of some "meta state," for example, a change in the sweepline. But I opted against that because this would cause redundant calculations since the relative order of the tree is still unchanged. That's why I do the recompute on the fly strategy.

Verticals

Vertical segments are a special case. They must always be considered the "top" segment at any point they pass through. But I couldn’t just set their xintersect to their top endpoint because it would miss every intersection between p and the endpoint. Instead, I make it relative to the current sweep position at p.y, while ensuring it’s still above everything else. This effectively "pushes" the vertical segment to the top of the tree as needed.

[juliohm]

Thank you @souma4 for the detailed comments. I understand that your _SweepSegment struct is an auxiliary type with custom isless for the binary tree. And that this is only needed because of the +/-Inf coordinates. Is that correct?

Could we try that other approach you suggested where we add two additional sentinel segments in the beginning of the algorithm using the boundingbox of the original set of segments? If you do box = segs |> boundingbox |> Stretch(factor), you get a bounding box that covers the original set of segments with some margin. You can then take the left and right segments of the boundary as the sentinels, which will certainly be positioned in the first and last position of the sweep line. Any iteration of the algorithm will update the active segments intersecting the sweep line without ever going beyond these sentinel segments. Makes sense? Please let me know if I am missing something.

Ideally, we could achieve the desired results with two tricks besides the built-in intersect predicate:

1. Quantization of coordinates in the beginning of the function and throughout the intersections (our digits keyword argument) to avoid floating point issues.
2. Addition of two sentinel segments at "infinity" using the boundingbox of the original set of segments.

[souma4]

So the auxiliary type is there because we need to order the status structure by where the segments intersect the sweepline. The +/-Inf does not affect that. I'd say the _SweepSegment structure is necessary for correct ordering, regardless of whether it uses intersect or not.

Where the +/-Inf comes into play is in the computation of intersections. That would allow us to use the intersect predicate. My thinking is to use the ymin and ymax from the bounding box (plus stretch factor) to construct the sweepline and modify x to update the sweepline position.