Convex Hull Algorithm - Graham scan fastest compare function?

Question

I am already implemented Graham scan, but i see that the bottleneck of my program is the sorting (80% of the time). I want to improve in it, now I am doing the following:

 std::sort(intersections.begin() + 1, intersections.end(), [&minElement](Point const& a, Point const& b)
 {return angle (minElement - a, XAxis) < angle (minElement - b,XAxis);});

Which gives me the exact angle, but it's not cheap because my angle function looks like the following:

float angle (const Point& v1, const Point& v2) {
    return dot (v1, v2) / (v1.Length () * v2.Length ());
}

In the Length function I have to do a square root, which is one of the most expensive operations. But this way I get a good ordering.

I tried ordering the array by Slope, dot, ccw, and even only removing sqrt from my compare, but none of that provides me the same ordering. Can you give me any advice?

Answer 1

When you're sorting the points by their relative angles, you don't need to know the exact angle the two points make. Rather, you just need to know whether one point is to the left or to the right of the other.

Image, for example, you want to compare two points (x₁, y₁) and (x₂, y₂), assuming that the bottommost point is at (x_p, y_p). Look at the two vectors v₁ = (x₁ - x_p, y₁ - y_p) and v₂ = (x₂ - x_p, y₂ - y_p). To determine whether v₁ is to the left or to the right of v₂, which means you want to look at the sign of the angle made going from v₁ to v₂. If it's positive, then v₂ is on the left of v₁, and if it's negative, then v₁ is on the left of v₂.

The 2D cross product has the nice property that

v₁ × v₂ = |v₁| |v₂| sin(θ)

where θ is the angle made by going from v₁ to v₂. This means that θ > 0 if v₁ is to the right of v₂ and vice-versa, which is nice because that lets you compare which one goes where purely by taking a cross product!

In other words:

• If v₁ × v₂ > 0, then v₂ is to the left of v₁.
• If v₁ × v₂ = 0, then the points are collinear.
• If v₁ × v₂ < 0, then v₂ is to the right of v₁.

The 2D cross product formula is given by

(Δx₁, Δy₁) × (Δx₂, Δy₂) = (Δx₁ Δy₂ - Δx₂ Δy₁)

Where, here, Δx₁ represents x₁ - x_p, etc.

So you could compute the above quantity, then look at its sign to determine how the two points relate to one another. No square roots needed!

Answer 2

You can make ad hoc optimization through caching.

First, Length() can cache its result in member variable. You only need to invalidate this value in case point is changed (and probably it's not the case). You can make lazy calculation, so Length will check if it has value and calculate/store if it doesn't, then return stored value.

Second, for given minElement and xAxis, angle (minElement - a, XAxis) becomes function of 1 argument, so you can save (cache) values of it for each point before sorting and use comparator on prepared values.

Naive approach to these things: use subclass of Point in main algo, so every point already has necessary methods and place for cached values.