37

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

danharaj
  • 1,613
  • 14
  • 19

1 Answers1

4

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.

n. m. could be an AI
  • 112,515
  • 14
  • 128
  • 243
  • Most haskell maps (e.g. `Data.Map`, `unordered-containers`) don't have a function to return neighbors of a given element, but you can use `split` and then `minView` or `maxView` as needed. – John L Jun 27 '11 at 07:55
  • 1
    @John L: One can also drop the set/map abstraction and use some kind of balanced search tree under its own name. – n. m. could be an AI Jun 27 '11 at 07:59
  • 1
    @n.m. - true, but that's how `Data.Map.Map` is implemented anyway, as I understand it. That's why map keys have an `Ord` constraint. Of course you may be able to find a more efficient implementation for any given purpose, plus then you wouldn't need to perform a split+rebalance (presuming Data.Map does so), which is probably costly. – John L Jun 27 '11 at 09:33
  • @John L: Yes, both `Map` and `Set` are balanced binary trees. `IntMap` and `IntSet` are not, which is why they're distinct from just using `Int` as a key. – C. A. McCann Jun 27 '11 at 13:04
  • I understand the ordering a bit more now. How do I deal with the special case of segments that are collinear with the sweep line? – danharaj Jun 27 '11 at 19:18
  • 1
    @danharaj: Special cases like this are hard. This particular case can be handled by perturbing the input ever so slightly, such that no segment is vertical. There are other cases, like multiple collinear segments, that Bentley-Ottmann doesn't handle very well. There's an entire separate science about handling degenerate cases in comp-geom. Oh, and I hope you know that one generally needs to use rational numbers, not any floating point with rounding, to avoid nasty topological problems. There are ways around this, though, see e.g. [this](http://ect.bell-labs.com/who/hobby/93_2-27.pdf) (pdf). – n. m. could be an AI Jun 27 '11 at 19:47
  • @n.m. : Tell me about it! I think it's the degenerate cases that are tripping me up here. As for precision issues: My intent is not to be precise, but to provide some comp. geometry for interactive real-time graphical programs like video games, so it's important to be fast and 'good enough' instead of exact. Are there any other segment intersection algorithms I should look at? I looked at space partitioning algorithms, but as far as I could tell, they weren't output sensitive like line sweeps. – danharaj Jun 27 '11 at 21:11
  • 1
    @danharaj: The problem with floating point is not that it's inexact. The problem is that the algorithms simply stop working. Intersection points appear out of nowhere, the paper that I linked has an example. This is more dangerous than it sounds, as each intersection involves swapping segment order. A false intersection then leads to wrong Y-ordering of the segments, and consequently, makes the algorithm miss more intersections, which makes the Y-ordering of more segments wrong, etc. – n. m. could be an AI Jun 27 '11 at 22:15
  • @n.m.: Ah I see. I am reading the paper now. Is there a way I can side-step the issue by using extra precision in the calculation? Like storing data as `Float`s and computing with `Double`s? Oi, this stuff takes a lot of thought! – danharaj Jun 27 '11 at 23:02
  • 1
    @danharaj: Why not use `Rational`, Haskell already has it? It's slower than `Double` but at least the algo will work. I think just having the source coords in `Float` is not enough, you need an algo specially designed for that; I remember reading. a paper on it, but can't seem to find it right now. – n. m. could be an AI Jun 28 '11 at 02:36
  • @n.m.: I could use rational, I guess. But it would be slightly inconveniencing. I am using `Float` because that is how I represent vertices to be rendered by OpenGL. It would be convenient to use the same representation for rendering as for computation. I looked for papers, and I found [this one](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.91.1578&rep=rep1&type=pdf). It looks promising. I'll continue to explore this avenue. I want to conclude this question when I have a working robust implementation, hopefully soon. – danharaj Jun 28 '11 at 04:26