Search code examples
performancesortingmathgeometryangle

Fastest way to sort vectors by angle without actually computing that angle


Many algorithms (e.g. Graham scan) require points or vectors to be sorted by their angle (perhaps as seen from some other point, i.e. using difference vectors). This order is inherently cyclic, and where this cycle is broken to compute linear values often doesn't matter that much. But the real angle value doesn't matter much either, as long as cyclic order is maintained. So doing an atan2 call for every point might be wasteful. What faster methods are there to compute a value which is strictly monotonic in the angle, the way atan2 is? Such functions apparently have been called “pseudoangle” by some.


Solution

  • I started to play around with this and realised that the spec is kind of incomplete. atan2 has a discontinuity, because as dx and dy are varied, there's a point where atan2 will jump between -pi and +pi. The graph below shows the two formulas suggested by @MvG, and in fact they both have the discontinuity in a different place compared to atan2. (NB: I added 3 to the first formula and 4 to the alternative so that the lines don't overlap on the graph). If I added atan2 to that graph then it would be the straight line y=x. So it seems to me that there could be various answers, depending on where one wants to put the discontinuity. If one really wants to replicate atan2, the answer (in this genre) would be

    # Input:  dx, dy: coordinates of a (difference) vector.
    # Output: a number from the range [-2 .. 2] which is monotonic
    #         in the angle this vector makes against the x axis.
    #         and with the same discontinuity as atan2
    def pseudoangle(dx, dy):
        p = dx/(abs(dx)+abs(dy)) # -1 .. 1 increasing with x
        if dy < 0: return p - 1  # -2 .. 0 increasing with x
        else:      return 1 - p  #  0 .. 2 decreasing with x
    

    This means that if the language that you're using has a sign function, you could avoid branching by returning sign(dy)(1-p), which has the effect of putting an answer of 0 at the discontinuity between returning -2 and +2. And the same trick would work with @MvG's original methodology, one could return sign(dx)(p-1).

    Update In a comment below, @MvG suggests a one-line C implementation of this, namely

    pseudoangle = copysign(1. - dx/(fabs(dx)+fabs(dy)),dy)
    

    @MvG says it works well, and it looks good to me :-).

    enter image description here