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.

Answers


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 :-).


I know one possible such function, which I will describe here.

# Input:  dx, dy: coordinates of a (difference) vector.
# Output: a number from the range [-1 .. 3] (or [0 .. 4] with the comment enabled)
#         which is monotonic in the angle this vector makes against the x axis.
def pseudoangle(dx, dy):
    ax = abs(dx)
    ay = abs(dy)
    p = dy/(ax+ay)
    if dx < 0: p = 2 - p
    # elif dy < 0: p = 4 + p
    return p

So why does this work? One thing to note is that scaling all input lengths will not affect the ouput. So the length of the vector (dx, dy) is irrelevant, only its direction matters. Concentrating on the first quadrant, we may for the moment assume dx == 1. Then dy/(1+dy) grows monotonically from zero for dy == 0 to one for infinite dy (i.e. for dx == 0). Now the other quadrants have to be handled as well. If dy is negative, then so is the initial p. So for positive dx we already have a range -1 <= p <= 1 monotonic in the angle. For dx < 0 we change the sign and add two. That gives a range 1 <= p <= 3 for dx < 0, and a range of -1 <= p <= 3 on the whole. If negative numbers are for some reason undesirable, the elif comment line can be included, which will shift the 4th quadrant from -1…0 to 3…4.

I don't know if the above function has an established name, and who might have published it first. I've gotten it quite a while ago and copied it from one project to the next. I have however found occurrences of this on the web, so I'd consider this snipped public enough for re-use.

There is a way to obtain the range [0 … 4] (for real angles [0 … 2π]) without introducing a further case distinction:

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

I kinda like trigonometry, so I know the best way of mapping an angle to some values we usually have is a tangent. Of course, if we want a finite number in order to not have the hassle of comparing {sign(x),y/x}, it gets a bit more confusing.

But there is a function that maps [1,+inf[ to [1,0[ known as inverse, that will allow us to have a finite range to which we will map angles. The inverse of the tangent is the well known cotangent, thus x/y (yes, it's as simple as that).

A little illustration, showing the values of tangent and cotangent on a unit circle :

You see the values are the same when |x| = |y|, and you see also that if we color the parts that output a value between [-1,1] on both circles, we manage to color a full circle. To have this mapping of values be continuous and monotonous, we can do two this :

  • use the opposite of the cotangent to have the same monotony as tangent
  • add 2 to -cotan, to have the values coincide where tan=1
  • add 4 to one half of the circle (say, below the x=-y diagonal) to have values fit on the one of the discontinuities.

That gives the following piecewise function, which is a continuous and monotonous function of the angles, with only one discontinuity (which is the minimum) :

double pseudoangle(double dy, double dx)
{
    if( dy == 0 )
        return (double)(4 * (dx < 0));

    // 1 for above, 0 for below the diagonal/anti-diagonal
    int diag = dx > dy;
    int adiag = dx > -dy;

    double r;
    if( diag ^ adiag )
        r = dy / dx;
    else
        r = 2 - dx / dy;

    return 4 * (!adiag) + reduced;
}

Note that this is very close to Fowler angles, with the same properties. Formally, pseudoangle(dx,dy) + 1 % 8 == Fowler(dx,dy)

To talk performance, it's much less branchy than Fowler's (old) code (and generally less complicated imo). The minimum you need to do is one for 0 checking of either x or y, and another one which you could trade in for a floating-point division, though that doesn't sound great.

This doesn't rely on abs, which I find nice, but a good compiler wouldn't throw in a branch for that. However other solutions that forget to check whether abs(dx) + abs(dy) == 0 before dividing by it, would be equivalent to this code when removing the if( dy == 0 ).

Finally, I would argue this version is more precise, since it only uses mantissa preserving operations (until shifting the result to the right interval). This should be especially visible when |x| << |y| or |y| >> |x|, then the operation |x| + |y| looses quite some precision.


nice.. here is a varient that returns -Pi , Pi like many arctan2 functions.

edit note: changed my pseudoscode to proper python.. arg order changed for compatibility with pythons math module atan2(). Edit2 bother more code to catch the case dx=0.

def pseudoangle( dy , dx ):
  """ returns approximation to math.atan2(dy,dx)*2/pi"""
  if dx == 0 :
      s = cmp(dy,0)
  else::
      s = cmp(dx*dy,0)  # cmp == "sign" in many other languages.
  if s == 0 : return 0 # doesnt hurt performance much.but can omit if 0,0 never happens
  p = dy/(dx+s*dy)
  if dx < 0: return p-2*s
  return  p

In this form the max error is only ~0.07 radian for all angles. (of course leave out the Pi/2 if you don't care about the magnitude.)

Now for the bad news -- on my system using python math.atan2 is about 25% faster Obviously replacing a simple interpreted code doesnt beat a compiled intrisic.


Just use a cross-product function. The direction you rotate one segment relative to the other will give either a positive or negative number. No trig functions and no division. Fast and simple. Just Google it.


Need Your Help

How to decode Norwegian or any special character in a xml to display it in TextView

ios xml ipad nsstring special-characters

I am dealing with strings in my project where in i send some data (string) to a service and the service responses me XML with some data (String).Earlier i used to face probs while sending special

Golang crypto multiple calls have different responses

cryptography go hmac sha256 scrypt

I'm having a problem with some Go code I've written for a password authentication library. The general idea is to provide 2 functions, Check() and New(), which are both provided a password and a 25...

About UNIX Resources Network

Original, collect and organize Developers related documents, information and materials, contains jQuery, Html, CSS, MySQL, .NET, ASP.NET, SQL, objective-c, iPhone, Ruby on Rails, C, SQL Server, Ruby, Arrays, Regex, ASP.NET MVC, WPF, XML, Ajax, DataBase, and so on.