20

I'm working on a problem that the professor assigned, and I'm having a problem looking for a way to detect if the angle between 3 points is more than 180 degrees, e.g:

I want to detect if alpha is more than 180 degrees. Anyways, my professor has a code that solves the problem, but he has a function called zcross, but I don't exactly know how it works. Could anyone tell me? His code is here:

#include <fstream.h>
#include <math.h>
#include <stdlib.h>

struct point {
    double  x;
    double  y;
    double  angle;
};

struct vector {
    double  i;
    double  j;
};

point   P[10000];
int     hull[10000];

int 
zcross (vector * u, vector * v)
{
    double  p = u->i * v->j - v->i * u->j;
    if (p > 0)
    return 1;
    if (p < 0)
    return -1;
    return 0;
}

int 
cmpP (const void *a, const void *b)
{
    if (((point *) a)->angle < ((point *) b)->angle)
    return -1;
    if (((point *) a)->angle > ((point *) b)->angle)
    return 1;
    return 0;
}

void 
main ()
{
    int     N, i, hullstart, hullend, a, b;
    double  midx, midy, length;
    vector  v1, v2;

    ifstream fin ("fc.in");
    fin >> N;
    midx = 0, midy = 0;
    for (i = 0; i < N; i++) {
        fin >> P[i].x >> P[i].y;
        midx += P[i].x;
        midy += P[i].y;
    }
    fin.close ();
    midx = (double) midx / N;
    midy = (double) midy / N;
    for (i = 0; i < N; i++)
        P[i].angle = atan2 (P[i].y - midy, P[i].x - midx);
    qsort (P, N, sizeof (P[0]), cmpP);

    hull[0] = 0;
    hull[1] = 1;
    hullend = 2;
    for (i = 2; i < N - 1; i++) {
        while (hullend > 1) {
            v1.i = P[hull[hullend - 2]].x - P[hull[hullend - 1]].x;
            v1.j = P[hull[hullend - 2]].y - P[hull[hullend - 1]].y;
            v2.i = P[i].x - P[hull[hullend - 1]].x;
            v2.j = P[i].y - P[hull[hullend - 1]].y;
            if (zcross (&v1, &v2) < 0)
                break;
            hullend--;
        }
        hull[hullend] = i;
        hullend++;
    }

    while (hullend > 1) {
        v1.i = P[hull[hullend - 2]].x - P[hull[hullend - 1]].x;
        v1.j = P[hull[hullend - 2]].y - P[hull[hullend - 1]].y;
        v2.i = P[i].x - P[hull[hullend - 1]].x;
        v2.j = P[i].y - P[hull[hullend - 1]].y;
        if (zcross (&v1, &v2) < 0)
            break;
        hullend--;
    }
    hull[hullend] = i;

    hullstart = 0;
    while (true) {
        v1.i = P[hull[hullend - 1]].x - P[hull[hullend]].x;
        v1.j = P[hull[hullend - 1]].y - P[hull[hullend]].y;
        v2.i = P[hull[hullstart]].x - P[hull[hullend]].x;
        v2.j = P[hull[hullstart]].y - P[hull[hullend]].y;
        if (hullend - hullstart > 1 && zcross (&v1, &v2) >= 0) {
            hullend--;
            continue;
        }
        v1.i = P[hull[hullend]].x - P[hull[hullstart]].x;
        v1.j = P[hull[hullend]].y - P[hull[hullstart]].y;
        v2.i = P[hull[hullstart + 1]].x - P[hull[hullstart]].x;
        v2.j = P[hull[hullstart + 1]].y - P[hull[hullstart]].y;
        if (hullend - hullstart > 1 && zcross (&v1, &v2) >= 0) {
            hullstart++;
            continue;
        }
        break;
    }

    length = 0;
    for (i = hullstart; i <= hullend; i++) {
        a = hull[i];
        if (i == hullend)
            b = hull[hullstart];
        else
            b = hull[i + 1];
        length += sqrt ((P[a].x - P[b].x) * (P[a].x - P[b].x) + (P[a].y - P[b].y) * (P[a].y - P[b].y));
    }

    ofstream fout ("fc.out");
    fout.setf (ios: :fixed);
    fout.precision (2);
    fout << length << '\n';
    fout.close ();
}
Bill the Lizard
  • 398,270
  • 210
  • 566
  • 880
Richie Li
  • 1,257
  • 2
  • 13
  • 20

4 Answers4

40

First, we know that if sin(a) is negative, then the angle is more than 180 degrees.

How do we find the sign of sin(a)? Here is where cross product comes into play.

First, let's define two vectors:

v1 = p1-p2
v2 = p3-p2

This means that the two vectors start at p2 and one points to p1 and the other points to p3.

Cross product is defined as:

(x1, y1, z1) x (x2, y2, z2) = (y1z2-y2z1, z1x2-z2x1, x1y2-x2y1)

Since your vectors are in 2d, then z1 and z2 are 0 and hence:

(x1, y1, 0) x (x2, y2, 0) = (0, 0, x1y2-x2y1)

That is why they call it zcross because only the z element of the product has a value other than 0.

Now, on the other hand, we know that:

||v1 x v2|| = ||v1|| * ||v2|| * abs(sin(a))

where ||v|| is the norm (size) of vector v. Also, we know that if the angle a is less than 180, then v1 x v2 will point upwards (right hand rule), while if it is larger than 180 it will point down. So in your special case:

(v1 x v2).z = ||v1|| * ||v2|| * sin(a)

Simply put, if the z value of v1 x v2 is positive, then a is smaller than 180. If it is negative, then it's bigger (The z value was x1y2-x2y1). If the cross product is 0, then the two vectors are parallel and the angle is either 0 or 180, depending on whether the two vectors have respectively same or opposite direction.

Shahbaz
  • 46,337
  • 19
  • 116
  • 182
  • 2
    In 2D, what you are really doing is computing the "outer product", which is a more general concept than cross product and works in any number of dimensions. They don't teach it in introductory linear algebra classes, which is a shame. (The formula is mostly the same, just with no mention of "z" coordinates, so it's simpler.) – Dietrich Epp Nov 03 '11 at 04:07
4

zcross is using the sign of the vector cross product (plus or minus in the z direction) to determine if the angle is more or less than 180 degrees, as you've put it.

Justin Peel
  • 46,722
  • 6
  • 58
  • 80
1

In 3D, find the cross product of the vectors, find the minimum length for the cross product which is basically just finding the smallest number of x, y and z.

If the smallest value is smaller than 0, the angle of the vectors is negative.

So in code:

float Vector3::Angle(const Vector3 &v) const
{
    float a = SquareLength();
    float b = v.SquareLength();
    if (a > 0.0f && b > 0.0f)
    {
        float sign = (CrossProduct(v)).MinLength();
        if (sign < 0.0f)
            return -acos(DotProduct(v) / sqrtf(a * b));
        else
            return acos(DotProduct(v) / sqrtf(a * b));
    }
    return 0.0f;
}
  • I think its important to mention, that the function returns a angle between [-180°;180°] - not a angle between [0;360°] - works perfect! – Vertexwahn Apr 26 '13 at 06:57
0

Another way to do it would be as follows:

calculate vector v1=p2-p1, v2 = p2 -p3. Then, use the cross-product formula : u.v = ||u|| ||v|| cos(theta)

grdvnl
  • 636
  • 6
  • 9