# Line segment intersection from [de Berg et al, 2000], computational geometry

I've implemented the following algorithm (From "Computational Geometry : Algorithms and Applications").

The code is below: (All the code actually is here: https://github.com/lukkio88/ComputationalGeometry/tree/master/Point, but the relevant bits for the algorithm are the .h and .cpp below)

The algorithm is quite complicated, this is why I'm asking for a review.

.h

#pragma once
#ifndef __LINE_SEG_INTERSECTION_H
#define __LINE_SEG_INTERSECTION_H
#include <segment.h>
#include <map>
#include <set>
#include <vector>

using std::map;
using std::set;
using std::pair;
using std::vector;

struct ComparePts {
bool operator()(const Point& p, const Point & q) const;
};

using PriorityQueue = map<Point, vector<Segment>, ComparePts>;
std::ostream& operator<<(std::ostream& os, const PriorityQueue& p);

Point getMin(const Point & p, const Point & q);
Point getMax(const Point & p, const Point & q);

Point getUp(const Segment &s);
Point getDown(const Segment & s);

struct SegmentComparator {
bool operator()(const Segment & s, const Segment & r) const;
bool * above;
Float * y;
};

using SweepLine = set<Segment, SegmentComparator>;
using SweepLineIter = SweepLine::iterator;
constexpr Float ths = 0.0001;

struct StatusStructure {

StatusStructure();

SweepLineIter getIncident(const Point& p);
bool findLeftNeighboor(Float x, Segment& sl) const;
bool findRightNeighboor(Float x, Segment& sr) const;
void findLeftmostAndRightmost(const Point& pt, SweepLineIter& it_l, SweepLineIter& it_r);

bool above;
Float y_line;
SegmentComparator segComp;
SweepLine sweepLine;
SweepLineIter nil;
};

std::ostream& operator<<(std::ostream& os, const StatusStructure& tau);

vector<Point> computeIntersection(vector<Segment> & S);

#endif


.cpp

#pragma once
#include <line_seg_intersection.h>

bool ComparePts::operator()(const Point& p, const Point & q) const {
return p.y - q.y >= ths || ((abs(p.y - q.y) < ths) && p.x < q.x);
}

Point getMin(const Point & p, const Point & q) {
ComparePts cmp;
if (cmp(p, q))
return p;
return q;
}

Point getMax(const Point & p, const Point & q) {
ComparePts cmp;
if (cmp(p, q))
return q;
return p;
}

Point getUp(const Segment &s) {
return getMin(s.p, s.q);
}

Point getDown(const Segment & s) {
return getMax(s.p, s.q);
}

bool SegmentComparator::operator()(const Segment& s, const Segment& r) const {
Float xs, xr;
bool
sIsHorizontal = s.isHorizontal(),
rIsHorizontal = r.isHorizontal();
if (sIsHorizontal)
xs = s.p.x;
else
s.getX(*y, xs);

if (rIsHorizontal)
xr = r.p.x;
else
r.getX(*y, xr);

if (xs != xr)
return xs < xr;
else {
Point u = (sIsHorizontal) ?
normalize(s.q - s.p):
normalize(getUp(s) - getDown(s));

Point v = (rIsHorizontal) ?
normalize(r.q - r.p):
normalize(getUp(r) - getDown(r));

Point ref{ 1.0,0.0 };

if (*above) {
return u * ref < v*ref;
}
else {
return u * ref > v*ref;
}

}
}

SweepLineIter StatusStructure::getIncident(const Point& p) {
return sweepLine.lower_bound(Segment{ p,p + Point{-1.0,0.0} });
}

std::ostream& operator<<(std::ostream& os, const PriorityQueue& p) {
for (auto el : p) {
std::cout << el.first << std::endl;
for (auto seg : el.second) {
std::cout << seg << std::endl;
}
}
return os;
}

static int size(const vector<Segment>& U, const vector<Segment>& C, const vector<Segment>& L) {
return U.size() + C.size() + L.size();
}

static int size(const vector<Segment>& U, const vector<Segment>& C) {
return U.size() + C.size();
}

static bool findEvent(const Segment& l, const Segment& r, Point& p) {
return l.intersect(r, p);
}

StatusStructure::StatusStructure() {
segComp.y = &y_line;
segComp.above = &above;
sweepLine = SweepLine(segComp);
nil = sweepLine.end();
}

bool StatusStructure::findLeftNeighboor(Float x, Segment& sl) const { //This assumes the flag "above" is false
Segment tmp{ Point{0.0,0.0},Point{1.0,0.0} };
SweepLineIter it = sweepLine.lower_bound(tmp);
while (it != nil && (--it) != nil) {
Float curr_x;
it->getX(y_line, curr_x);
if (curr_x != x) {
sl = *it;
return true;
}
}
return false;
}

bool StatusStructure::findRightNeighboor(Float x, Segment& sr) const {
Segment tmp{ Point{0.0,0.0},Point{1.0,0.0} };
SweepLineIter it = sweepLine.lower_bound(tmp);
while (it != nil && (++it) != nil) {
Float curr_x;
it->getX(y_line, curr_x);
if (curr_x != x) {
sr = *it;
return true;
}
}
return false;
}

/*
This function will find the leftMost and the rightMost in a star of segments
passing through pt
*/
void StatusStructure::findLeftmostAndRightmost(const Point& pt,SweepLineIter& it_l, SweepLineIter& it_r) {
Float x;
//Getting the segment whose dot product with e1 is >= to -1
it_l = sweepLine.lower_bound({ pt,pt + Point{-1.0,0.0} }); //this will return the actual iterator to the segment, this must exist
it_r = sweepLine.upper_bound({ pt,pt + Point{1.0,0.0} }); //this potentially might be nil
it_r--;
}

ostream& operator<<(ostream & os, const StatusStructure & tau)
{
std::string curr_str;
for (auto& seg : tau.sweepLine)
os << seg.label << " ";
return os;
}

vector<Point> computeIntersection(vector<Segment> & S) {
PriorityQueue queue;
while (!S.empty()) {
Segment s = S.back();
queue[getUp(s)].push_back(s);
queue[getDown(s)];
S.pop_back();
}

vector<Point> intersections;

//Init status structure
StatusStructure tau;
std::vector<Segment> C, L;
Float curr_x;
Segment sl, sr;

while (!queue.empty()) {

Point p = queue.begin()->first;
tau.y_line = p.y;
tau.above = true;
std::vector<Segment> U = queue.begin()->second;
queue.erase(queue.begin());

SweepLineIter it = tau.getIncident(p);

//populating L and C
while (it != tau.nil && it->getX(tau.y_line, curr_x) && (abs(curr_x - p.x) < 0.0001)) {
if (getDown(*it) == p)
L.push_back(*it);
else
C.push_back(*it);
it = tau.sweepLine.erase(it);
}

if (size(U, C, L) > 1)
intersections.push_back(p);

while (!L.empty())
L.pop_back();

tau.above = false;
int size_UC = size(U, C);

while (!U.empty()) {
tau.sweepLine.insert(U.back());
U.pop_back();
}
while (!C.empty()) {
tau.sweepLine.insert(C.back());
C.pop_back();
}

if (size_UC == 0) {
if (tau.findLeftNeighboor(p.x, sl) && tau.findRightNeighboor(p.x, sr)) {
Point new_event_pt;
if (findEvent(sl, sr, new_event_pt))
queue[new_event_pt];
}
}
else {
tau.above = true;
SweepLineIter it_l, it_r, it_ll, it_rr;
tau.findLeftmostAndRightmost(p, it_l, it_r);
it_ll = it_l;
it_rr = it_r;
--it_ll;
++it_rr;

if (it_ll != tau.nil) {
Point new_event_pt;
if (findEvent(*it_ll, *it_l, new_event_pt))
queue[new_event_pt];
}
if (it_rr != tau.nil) {
Point new_event_pt;
if (findEvent(*it_r, *it_rr, new_event_pt))
queue[new_event_pt];
}
}
}

return intersections;

}

• There's big chunks missing from this code - Point and Segment are not defined. Are they in the <segment.h> that isn't shown? I assume that's a library header (since it's <> rather than ""), but you didn't say which library it is! – Toby Speight Jun 10 '19 at 9:13
• @TobySpeight They are not existing libraries, but implemented in the referenced github repository. Perhaps it is better to include them here, a large part of my answer is actually about segment.cpp and point.cpp. – Discrete lizard Jun 10 '19 at 9:55

I'm not an expert in c++, so I will focus mostly on the correctness of your implementation w.r.t. the algorithm as described in the book.

## Numerical Robustness

Often, the main challenge in implementing geometric algorithms is the fact we cannot work with the real numbers and therefore implementing geometric primitives (is this point to the left of a line, find the intersection of two lines) is tricky. There are multiple ways to implement them, and combinations of them can be used as well:

• Floating point arithmetic, which is fast and has native language support (double), but imprecise.
• Fraction arithmetic, which is precise, but slow and often has no native language support. You also have to be very careful to avoid integer overflow (or use arbitrary precision integer arithmetic such as Java's BigInteger. I'm not aware of any comparable implementation in c++, however)
• Interval arithmetic. Is precise and relatively fast in the 'easy' cases, defers to exact arithmetic in the 'hard' cases. Doing this effectively, however, relies on the intricacies of the floating point standard and often some calculation of the error propagation within .

In general, numerical robustness is a hard problem to deal with efficiently. (These slides on some considerations of these issues within the CGAL library are nice to give an idea of the complexity involved)

You should decide on which approach to take depending on the purpose of your implementation. If your main purpose is to learn by yourselves, you might get away with using double (properly!) in the instances you test it on. If you are getting issues with the precision, consider using fraction arithmetic. If you want to use this algorithm in production, don't. Instead, use the naive O(n^2) line segment algorithm or use CGAL.

For the following, I will just assume that this implementation is supposed to be only educational.

## Implementing the geometric primitives

This is not in the class you posted, but it is important to do this correctly, so I will mention this anyway. I will first make some remarks on your current approach.

return (twiceArea(c, d, a)*twiceArea(c, d, b) <= 0.0) && (twiceArea(a, b, c)*twiceArea(a, b, d) <= 0.0);


You cannot rely on equality holding when the numbers are equal. Instead, you should consider everthing that is less than a small value as equal to zero (as you have done correctly in ComparePts):

double eps = 0.0001;
return (twiceArea(c, d, a)*twiceArea(c, d, b) < eps) && (twiceArea(a, b, c)*twiceArea(a, b, d) < eps);


But be aware that this will give an incorrect answer if there are true values in between 0 and eps or if the error in the calculation of twiceArea(c, d, a)*twiceArea(c, d, b) is larger than eps.

Also, twiceArea is not a good name. Areas should be positive, and it also isn't clear of what you are computing the area. Additionally, you only need the sign of this operation. I suggest you replace twiceArea by a function orientation that returns its sign, where you should set the sign to 0 if the value is less than some small constant and otherwise return the sign properly (1 if positive, -1 if negative)

To compute whether two segments intersect, I suggest the following approach:

Using orientation, you can easily test whether a segment intersects a line: this is the case if an only if one of the endpoints of the segment lies on the line or both endpoints lie on different sides of the line. With that, you can test whether two line segments intersect: test whether the first segment intersects the line extended from the second segment and the vice versa.

The advantage of this approach is that you replace the test sign(a * b) < 0 by sign(a) != sign(b) and avoid a floating point operation. It also places all the imprecision 'inside' the smallest function possible, which makes it easier to adapt it to more precise approaches or different ways of dealing with the imprecision.

As for Segment::intersect, I would instead implement a function to compute the intersection of lines extending from the segments. This means that you drop the check for checking whether the segments intersect, although you probably want to test for the case where the lines are parallel and do not intersect. One reason to do this is that you can then use this function to compute the intersection point of a segment with the sweepline, it it intersects the sweepline. So, you can replace getX with a call to the the line segment intersection function.

## Segment comparison

With these low-level implementation issues out of the way, let's look at the rest of the code.

In the function SegmentComparator:

if (xs != xr)
return xs < xr;


Two direct comparison of floats, do something like

if (xs - xr < ths)
return 1
if (xs - xr > ths)
return 0


instead. (The else is unnecessary.)

    Point u = (sIsHorizontal) ?
normalize(s.q - s.p):
normalize(getUp(s) - getDown(s));

Point v = (rIsHorizontal) ?
normalize(r.q - r.p):
normalize(getUp(r) - getDown(r));

Point ref{ 1.0,0.0 };

if (*above) {
return u * ref < v*ref;
}
else {
return u * ref > v*ref;
}


It seems you are comparing the inverse of the slope here. I would not normalize here and compute the inverse slope normally ((s.q.y-s.p.y)/(s.p.x-s.q-x)) which avoids having to determine which part is 'up', or at least not take the square root, that is not needed it makes it impossible to use fractions. (you can simply compare the squared Euclidean norm)

You are using above to switch the order of segments in the status at their intersection point. This is ok, but you have to be careful when you insert two segments with the same upper endpoint (when you put all eventpoints except the intersections in the queue). In that case, the above flag should be false. (it is not clear to me what the initial value of above is)

## Event handling

Set and Map use a red-black tree as the underlying data-structure, so they have the required logarithmic time for insertion, removal and lookup. So far so good.

    std::vector<Segment> U = queue.begin()->second;


This is incorrect, p does not have to be an upper point of all segments stored at the event point, the segments can either have this point as an intersection point or as lower point. You can check at this point to which part the points belong, but I'd recommend to store the type of point when inserting the event, as it is always known at that point. This also helps with visualizing your event queue for debugging.

static bool findEvent(const Segment& l, const Segment& r, Point& p) {
return l.intersect(r, p);
}


It is not enough to check whether the segments intersect, you have to verify that the intersection is not in the event queue, or any point that has been in the event queue. To do this, you have to check that the intersection point p is earlier than the current event point with ComparePts.

I'm not sure if there are more errors here, but I suggest you implement the simple case where an event-point can only be either 1. the top of a segment 2. the bottom of a segment or 3. the intersection of a pair of segments. I would only extend this to the general case if I test that that works.

## General recommendations

Implementing geometric algorithms properly is hard, so while it might be nice to try and implement the simpler ones such as the convex hull algorithms, you should keep in mind that implementing geometric algorithms is a completely different topic than what this textbook is trying to teach you.

Also, it helps to make some simplifying assumptions and test your implementation under those conditions first. Assuming that the points are in general position or that no points have the same y-coordinate simplifies the problem. (but be careful that this has to be the case for the intersection points as well)

• Thanks for your answer (finally I got one...). About the geometric primitives I'm planning to write better ones. About numerical robustness I also agree of course, but since this code is essentially not going to be used anywhere for now I'll just work with the assumption that double are good enough. I'm planning to implement few of the algorithms of that book. And later work on the numerical robustness. I'm more interested in manipulating the half edge data structure at very low level at the moment. – user8469759 Jun 9 '19 at 22:41
• About the simple case against the generic one, I'm only interested in the generic one since I'm currently implementing the map overlay algorithm, and I need the generic case. There're few questions I have about your answer, I'll collect them together and ask them later. – user8469759 Jun 9 '19 at 22:43
• @user8469759 Well, I understand you will at least need to support multiple endpoints of a segment at the same point. My main point is mostly that you should implement a simpler version first, test that, find bugs and only then implement the generic case. – Discrete lizard Jun 10 '19 at 6:39
• First question : You're pointing out my event handling is wrong, why is that? Line 1 of the pseudocode defines $U$ in the same way I did. – user8469759 Jun 10 '19 at 7:51
• @user8469759 No, or at least, I don't think your code selects the same U. The set U(p) is the set of segments that have p as their upper endpoint. Your code seems to assign the set of all segments that have p as an event-point to U. So, your code assigns some segments to U that should not be assigned to it, these are the segments that have p as intersection point or as lower endpoint. – Discrete lizard Jun 10 '19 at 8:08

Why is there #pragma once in the implementation file? Including it twice in the same translation unit is clearly an error that should be fixed, not a normal state of affairs as a header would be.

If the abs() is meant to be std::abs(), we need #include <cmath> and using std::abs;.

The getMin() and getMax functions seem to be exactly like std::min(p, q, ComparePts{}) and std::max(p, q, ComparePts{}) respectively, so provide very little value as functions.

I didn't get beyond this, as the rest seems to depend too much on the library that defines Float, Point and Segment, and there's no clue which of the many geometry libraries is used.

• There's no external library apart from the stl. I didn't feel like posting the whole, I was more looking potentially for people who have an undestanding of computational geometry to review that specific function. About Point and Segment the clue is in the link I gave. – user8469759 Jun 10 '19 at 10:21