Angle between 3 points?

Question

Given points ABC, how could I find angle ABC? I'm making a feehand tool for a vector drawing application and to minimize the number of points it generates, I wont add a points unless the angle of the mouse position and the last 2 points is greater than a certain threshold. Thanks

what I had:

int CGlEngineFunctions::GetAngleABC( POINTFLOAT a, POINTFLOAT b, POINTFLOAT c )
{
    POINTFLOAT ab;
    POINTFLOAT ac;

    ab.x = b.x - a.x;
    ab.y = b.y - a.y;

    ac.x = b.x - c.x;
    ac.y = b.y - c.y;

    float dotabac = (ab.x * ab.y + ac.x * ac.y);
    float lenab = sqrt(ab.x * ab.x + ab.y * ab.y);
    float lenac = sqrt(ac.x * ac.x + ac.y * ac.y);

    float dacos = dotabac / lenab / lenac;

    float rslt = acos(dacos);
    float rs = (rslt * 180) / 3.141592;
     RoundNumber(rs);
     return (int)rs;


}

Answer 1

First suggestions regarding your method:

What you call ac is actually cb. But it's ok, this is what really needed. Next,

float dotabac = (ab.x * ab.y + ac.x * ac.y);

This is your first mistake. The real dot product of two vectors is:

float dotabac = (ab.x * ac.x + ab.y * ac.y);

Now,

float rslt = acos(dacos);

Here you should note that due to some precision loss during the calculation it's theoretically possible that dacos will become bigger than 1 (or lesser than -1). Hence - you should check this explicitly.

Plus a performance note: you call a heavy sqrt function twice for calculating the length of two vectors. Then you divide the dot product by those lengths. Instead you could call sqrt on the multiplication of squares of length of both vectors.

And lastly, you should note that your result is accurate up to the sign. That is, your method won't distinguish 20° and -20°, since the cosine of both are the same. Your method will yield the same angle for ABC and CBA.

One correct method for calculating the angle is as "oslvbo" suggests:

float angba = atan2(ab.y, ab.x);
float angbc = atan2(cb.y, cb.x);
float rslt = angba - angbc;
float rs = (rslt * 180) / 3.141592;

(I've just replaced atan by atan2).

It's the simplest method, which always yields the correct result. The drawback of this method is that you actually call a heavy trigonometry function atan2 twice.

I suggest the following method. It's a bit more complex (requires some trigonometry skills to understand), however it's superior from the performance point of view. It just calls once a trigonometry function atan2. And no square root calculations.

int CGlEngineFunctions::GetAngleABC( POINTFLOAT a, POINTFLOAT b, POINTFLOAT c )
{
    POINTFLOAT ab = { b.x - a.x, b.y - a.y };
    POINTFLOAT cb = { b.x - c.x, b.y - c.y };

    // dot product  
    float dot = (ab.x * cb.x + ab.y * cb.y);

    // length square of both vectors
    float abSqr = ab.x * ab.x + ab.y * ab.y;
    float cbSqr = cb.x * cb.x + cb.y * cb.y;

    // square of cosine of the needed angle    
    float cosSqr = dot * dot / abSqr / cbSqr;

    // this is a known trigonometric equality:
    // cos(alpha * 2) = [ cos(alpha) ]^2 * 2 - 1
    float cos2 = 2 * cosSqr - 1;

    // Here's the only invocation of the heavy function.
    // It's a good idea to check explicitly if cos2 is within [-1 .. 1] range

    const float pi = 3.141592f;

    float alpha2 =
        (cos2 <= -1) ? pi :
        (cos2 >= 1) ? 0 :
        acosf(cos2);

    float rslt = alpha2 / 2;

    float rs = rslt * 180. / pi;


    // Now revolve the ambiguities.
    // 1. If dot product of two vectors is negative - the angle is definitely
    // above 90 degrees. Still we have no information regarding the sign of the angle.

    // NOTE: This ambiguity is the consequence of our method: calculating the cosine
    // of the double angle. This allows us to get rid of calling sqrt.

    if (dot < 0)
        rs = 180 - rs;

    // 2. Determine the sign. For this we'll use the Determinant of two vectors.

    float det = (ab.x * cb.y - ab.y * cb.y);
    if (det < 0)
        rs = -rs;

    return (int) floor(rs + 0.5);


}

EDIT:

Recently I've been working on a related subject. And then I've realized there's a better way. It's actually more-or-less the same (behind the scenes). However it's more straightforward IMHO.

The idea is to rotate both vectors so that the first one is aligned to (positive) X-direction. Obviously rotating both vectors doesn't affect the angle between them. OTOH after such a rotation one just has to find out the angle of the 2nd vector relative to X-axis. And this is exactly what atan2 is for.

Rotation is achieved by multiplying a vector by the following matrix:

• a.x, a.y
• -a.y, a.x

Once may see that vector a multiplied by such a matrix indeed rotates toward positive X-axis.

Note: Strictly speaking the above matrix isn't just rotating, it's also scaling. But this is ok in our case, since the only thing that matters is the vector direction, not its length.

Rotated vector b becomes:

• a.x * b.x + a.y * b.y = a dot b
• -a.y * b.x + a.x * b.y = a cross b

Finally, the answer can be expressed as:

int CGlEngineFunctions::GetAngleABC( POINTFLOAT a, POINTFLOAT b, POINTFLOAT c )
{
    POINTFLOAT ab = { b.x - a.x, b.y - a.y };
    POINTFLOAT cb = { b.x - c.x, b.y - c.y };

    float dot = (ab.x * cb.x + ab.y * cb.y); // dot product
    float cross = (ab.x * cb.y - ab.y * cb.x); // cross product

    float alpha = atan2(cross, dot);

    return (int) floor(alpha * 180. / pi + 0.5);
}

Answer 2

β = arccos((a^2 + c^2 - b^2) / 2ac)

where a is the side opposite angle α, b is opposite angle β, and c is opposite angle γ. So β is what you called angle ABC.

Answer 3

Approach with arccos is dangerous, because we risk to have it's argument equal to, say, 1.0000001 and end up with EDOMAIN error. Even atan approach is dangerous, because it involves divisions, which may lead to division by zero error. Better use atan2, passing dx and dy values to it.

Answer 4

Here is a quick and correct way of calculating the right angle value:

double AngleBetweenThreePoints(POINTFLOAT pointA, POINTFLOAT pointB, POINTFLOAT pointC)
{
    float a = pointB.x - pointA.x;
    float b = pointB.y - pointA.y;
    float c = pointB.x - pointC.x;
    float d = pointB.y - pointC.y;

    float atanA = atan2(a, b);
    float atanB = atan2(c, d);

    return atanB - atanA;
} 

Answer 5

Off topic? But you can do it with the law of cosines:

Find the distance between A and B (call this x), and the distance between B and C (call this y), and the distance between A and C (call this z).

Then you know that z^2=x^2+y^2-2*xycos(ANGLE YOU WANT)

therefore, that angle is cos^-1((z^2-x^2-y^2)/(2xy))=ANGLE

Answer 6

float angba = atan((a.y - b.y) / (a.x - b.x));
float angbc = atan((c.y - b.y) / (c.x - b.y));
float rslt = angba - angbc;
float rs = (rslt * 180) / 3.141592;