Bentley-Ottmann Algorithm in Haskell?

Question

So I've been writing a computational geometry library in Haskell because I couldn't find one on Hackage, and I figured it would be fun to do anyway. However, I've been stuck for nearly a week on one particular algorithm that I just can't seem to get into a nice 'haskell-like' form. That algorithm is the Bentley-Ottmann algorithm for finding the intersections in a set of line segments. If you are familiar with the algorithm, you can skip down to the last paragraph for my begging :)

The way I choose to implement this function is as a function that takes a list of line segments and returns a list of points, and the line segments that intersect at that point. This lets us deal with the case where multiple segments intersect at the same point.

bentleyOttmann :: [Segment] -> [(Point, [Segment])]

The algorithm is a sweep-line algorithm. We imagine a line sweeping the plane, doing algorithmic work at various even points. The event points in the Bentley-Ottmann algorithm are:

• The beginning end-point of a line segment.
• The ending end-point of a line segment.
• The intersection point of a bunch of segments.

Note that an event point can be associated with more than one line segment in more than one way. In order to keep track of which segments correspond to which end points, I use a map from the containers package. The keys of this map are points, and the values are lists of segments, labeled by whether they start at that point, end at that point, or intersect at that point.

The sweep-line determines the order of points. Imagine a vertical line sweeping the plane, stopping at event points in order to do work. Event points are ordered first by their x value, with smaller points being processed first. Generically, this is all we need. In degenerate cases, event points may have the same x-coordinate. We also order by their y coordinates, event points with smaller y coordinates are processed first if there is an x-coordinate tie.

So the structure I use, naturally, is a priority queue. The one I use is from the heaps package from Hackage.

What is the work that we do at each event point? Well, first we check which segments are associated to the event point. If there is more than one, it is an intersection point. We can add it to the list of intersections we have found so far.

Here comes the tricky part. While we are sweeping across the plane, we keep track of a set of segments, ordered with respect to the point where they intersect the sweep line. When we process an event point, we first remove all the segments that ended at that event point. Then, all the segments that intersected at that point are reversed in order. Finally, we add segments that begin at that event point to the ordered set. Note, that since these segments all intersect at the event point, they have to be ordered with respect to a sweep line that is slightly perturbed ahead.

At each event point, we must add any new event points, new intersections that occur. Because we keep track of the relative order of segments intersecting the sweep line, we do one of two things:

• If we swapped two segments or added a new segment, we find the bottommost (with respect to the sweep line) modified segment, the topmost modified segment, and test them for intersections with their immediate non-modified neighbors.

• If we did not swap or add new segments, then we at the very least removed a segment, thus making its former neighbors now adjacent. We test its these new neighbors for intersection.

This is the key to the Bentley-Ottmann algorithm, as we sweep across the plane, we only test new candidate segments with their neighbors. This means that we beat the naive O(n^2) algorithm when there are relatively few intersections.

My problem (finally, I'm sorry this is so long-winded) is this: I have no idea how to implement this ordering logic. I cannot use Data.Set because the ordering changes as we sweep. I am trying to implement my own data structure to keep track of the information, but it's grungy, buggy, probably inefficient, and also ugly! I hate ugly code.

I know Haskell is all about pretty code. I also believe that if I cannot implement an algorithm in a pretty way, it means I don't really understand it. Can anyone give me an insight into cleanly implementing this algorithm?

Edit: I have a 'working' implementation now. I intended it to work with generic input, as well as multiple segments intersecting at the same point, and vertical segments. It seems to work on those inputs with the paltry tests I did. It does not work when segments are overlapping. I do not know how to deal with those yet. I would appreciate input on how to accommodate them. Currently, my sweep line structure keeps track of them in the same node, but it will only use one of them in intersection tests, and can give inconsistent results.

I use Data.Set for my event queue, Data.Map for lookup, and an implementation of red-black trees with zippers that I based off of Okasaki's in his book. If my snippet is not enough context, I can add more.

I would appreciate tips on restructuring the implementation so it's... less ugly. I can't tell how correct it is and it makes me nervous.

The code can be found on hpaste here

Answer 1

The ordering if the segments change only at intersection points, and only for the segments that do intersect at the given point. This can be implemented by removing the intersecting segments and inserting them again.

The ordering function is by y coordinate, and when ys are equal, by the slope. The intersecting segments will then be inserted in the correct order. As the sweep progresses, the actual y coordinates of the segments' intersections of the sweep line will change. It doesn't matter as the order will remain the same (until we swap, that is, remove and insert again, intersecting segments). The actual y coordinate need not be stored anyway. It should be calculated dynamically for any given position of the sweep line, as we insert or remove segments.

The data structure in question should not be called Set, it's a Map or, more precisely, an Ordered Map. The operation of finding neighbours of a given element is essential here.