计算两个向量顺时针夹角的直接方法

我想找出2个向量(2D,3D)之间的顺时针角度。

点乘的经典方法给出了内角(0-180度) ,我需要使用一些 if 语句来确定结果是我需要的角度还是它的补数。

你知道一种计算顺时针角度的直接方法吗?

99092 次浏览

Scalar (dot) product of two vectors lets you get the cosinus of the angle between them. To get the 'direction' of the angle, you should also calculate the cross product, it will let you check (via z coordinate) is angle is clockwise or not (i.e. should you extract it from 360 degrees or not).

To compute angle you just need to call atan2(v1.s_cross(v2), v1.dot(v2)) for 2D case. Where s_cross is scalar analogue of cross production (signed area of parallelogram). For 2D case that would be wedge production. For 3D case you need to define clockwise rotation because from one side of plane clockwise is one direction, from other side of plane is another direction =)

Edit: this is counter clockwise angle, clockwise angle is just opposite

If by "direct way" you mean avoiding the if statement, then I don't think there is a really general solution.

However, if your specific problem would allow loosing some precision in angle discretization and you are ok with loosing some time in type conversions, you can map the [-pi,pi) allowed range of phi angle onto the allowed range of some signed integer type. Then you would get the complementarity for free. However, I didn't really use this trick in practice. Most likely, the expense of float-to-integer and integer-to-float conversions would outweigh any benefit of the directness. It's better to set your priorities on writing autovectorizable or parallelizable code when this angle computation is done a lot.

Also, if your problem details are such that there is a definite more likely outcome for the angle direction, then you can use compilers' builtin functions to supply this information to the compiler, so it can optimize the branching more efficiently. E.g., in case of gcc, that's __builtin_expect function. It's somewhat more handy to use when you wrap it into such likely and unlikely macros (like in linux kernel):

#define likely(x)      __builtin_expect(!!(x), 1)
#define unlikely(x)    __builtin_expect(!!(x), 0)

2D case

Just like the dot product is proportional to the cosine of the angle, the determinant is proprortional to its sine. So you can compute the angle like this:

dot = x1*x2 + y1*y2      # dot product between [x1, y1] and [x2, y2]
det = x1*y2 - y1*x2      # determinant
angle = atan2(det, dot)  # atan2(y, x) or atan2(sin, cos)

The orientation of this angle matches that of the coordinate system. In a left-handed coordinate system, i.e. x pointing right and y down as is common for computer graphics, this will mean you get a positive sign for clockwise angles. If the orientation of the coordinate system is mathematical with y up, you get counter-clockwise angles as is the convention in mathematics. Changing the order of the inputs will change the sign, so if you are unhappy with the signs just swap the inputs.

3D case

In 3D, two arbitrarily placed vectors define their own axis of rotation, perpendicular to both. That axis of rotation does not come with a fixed orientation, which means that you cannot uniquely fix the direction of the angle of rotation either. One common convention is to let angles be always positive, and to orient the axis in such a way that it fits a positive angle. In this case, the dot product of the normalized vectors is enough to compute angles.

dot = x1*x2 + y1*y2 + z1*z2    #between [x1, y1, z1] and [x2, y2, z2]
lenSq1 = x1*x1 + y1*y1 + z1*z1
lenSq2 = x2*x2 + y2*y2 + z2*z2
angle = acos(dot/sqrt(lenSq1 * lenSq2))

Edit: Note that some comments and alternate answers advise against the use of acos for numeric reasons, in particular if the angles to be measured are small.

Plane embedded in 3D

One special case is the case where your vectors are not placed arbitrarily, but lie within a plane with a known normal vector n. Then the axis of rotation will be in direction n as well, and the orientation of n will fix an orientation for that axis. In this case, you can adapt the 2D computation above, including n into the determinant to make its size 3×3.

dot = x1*x2 + y1*y2 + z1*z2
det = x1*y2*zn + x2*yn*z1 + xn*y1*z2 - z1*y2*xn - z2*yn*x1 - zn*y1*x2
angle = atan2(det, dot)

One condition for this to work is that the normal vector n has unit length. If not, you'll have to normalize it.

As triple product

This determinant could also be expressed as the triple product, as @Excrubulent pointed out in a suggested edit.

det = n · (v1 × v2)

This might be easier to implement in some APIs, and gives a different perspective on what's going on here: The cross product is proportional to the sine of the angle, and will lie perpendicular to the plane, hence be a multiple of n. The dot product will therefore basically measure the length of that vector, but with the correct sign attached to it.

For a 2D method, you could use the law of cosines and the "direction" method.

To calculate the angle of segment P3:P1 sweeping clockwise to segment P3:P2.

P1     P2


P3
    double d = direction(x3, y3, x2, y2, x1, y1);


// c
int d1d3 = distanceSqEucl(x1, y1, x3, y3);


// b
int d2d3 = distanceSqEucl(x2, y2, x3, y3);


// a
int d1d2 = distanceSqEucl(x1, y1, x2, y2);


//cosine A = (b^2 + c^2 - a^2)/2bc
double cosA = (d1d3 + d2d3 - d1d2)
/ (2 * Math.sqrt(d1d3 * d2d3));


double angleA = Math.acos(cosA);


if (d > 0) {
angleA = 2.*Math.PI - angleA;
}


This has the same number of transcendental

operations as suggestions above and only one more or so floating point operation.

the methods it uses are:

 public int distanceSqEucl(int x1, int y1,
int x2, int y2) {


int diffX = x1 - x2;
int diffY = y1 - y2;
return (diffX * diffX + diffY * diffY);
}


public int direction(int x1, int y1, int x2, int y2,
int x3, int y3) {


int d = ((x2 - x1)*(y3 - y1)) - ((y2 - y1)*(x3 - x1));


return d;
}

This answer is the same as MvG's, but explains it differently (it's the result of my efforts in trying to understand why MvG's solution works). I'm posting it on the off chance that others find it helpful.

The anti-clockwise angle theta from x to y, with respect to the viewpoint of their given normal n (||n|| = 1), is given by

atan2( dot(n, cross(x,y)), dot(x,y) )

(1) = atan2( ||x|| ||y|| sin(theta),  ||x|| ||y|| cos(theta) )

(2) = atan2( sin(theta), cos(theta) )

(3) = anti-clockwise angle between x axis and the vector (cos(theta), sin(theta))

(4) = theta

where ||x|| denotes the magnitude of x.

Step (1) follows by noting that

cross(x,y) = ||x|| ||y|| sin(theta) n,

and so

dot(n, cross(x,y))

= dot(n, ||x|| ||y|| sin(theta) n)

= ||x|| ||y|| sin(theta) dot(n, n)

which equals

||x|| ||y|| sin(theta)

if ||n|| = 1.

Step (2) follows from the definition of atan2, noting that atan2(cy, cx) = atan2(y,x), where c is a scalar. Step (3) follows from the definition of atan2. Step (4) follows from the geometric definitions of cos and sin.

A formula for clockwise angle,2D case, between 2 vectors, xa,ya and xb,yb.

Angle(vec.a-vec,b)=
pi()/2*((1+sign(ya))*
(1-sign(xa^2))-(1+sign(yb))*
(1-sign(xb^2))) +pi()/4*
((2+sign(ya))*sign(xa)-(2+sign(yb))*
sign(xb)) +sign(xa*ya)*
atan((abs(ya)-abs(xa))/(abs(ya)+abs(xa)))-sign(xb*yb)*
atan((abs(yb)-abs(xb))/(abs(yb)+abs(xb)))




                     

  

just copy & paste this.

angle = (acos((v1.x * v2.x + v1.y * v2.y)/((sqrt(v1.x*v1.x + v1.y*v1.y) * sqrt(v2.x*v2.x + v2.y*v2.y))))/pi*180);

you're welcome ;-)

Since one of the simplest and most elegant solutions is hidden in one the comments, I think it might be useful to post it as a separate answer. acos can cause inaccuracies for very small angles, so atan2 is usually preferred. For the 3D case:

dot = x1 * x2 + y1 * y2 + z1 * z2
cross_x = (y1 * z2 – z1 * y2)
cross_y = (z1 * x2 – x1 * z2)
cross_z = (x1 * y2 – y1 * x2)
det = sqrt(cross_x * cross_x + cross_y * cross_y + cross_z * cross_z)
angle = atan2(det, dot)

For 2D case atan2 can easily calculate angle between (1, 0) vector (X-axis) and one of your vectors. Formula is:

Atan2(y, x)

So you can easily calculate difference of two angles relatively X-axis

angle = -(atan2(y2, x2) - atan2(y1, x1))

Why is it not used as default solution? atan2 is not efficient enough. Solution from the top answer is better. Tests on C# showed that this method has 19.6% less performance (100 000 000 iterations). It's not critical but unpleasant.

So, another info that can be useful:

The smallest angle between outer and inner in degrees:

abs(angle * 180 / PI)

Full angle in degrees:

angle = angle * 180 / PI
angle = angle > 0 ? angle : 360 - angle

or

angle = angle * 180 / PI
if (angle < 0) angle = 360 - angle;