C++

Forum

 
A little geometry

kanesoban (39)
Helo,

I have the following problem:

I have points p1,p2,p3 (x and y coordinates). I need to determine, that if i draw a line between p1 and p2, drawing a line to p3 from p2 is taking a "left turn" or taking a "right turn". Some fast and easy formula would be nice.
Mitsakos (343)
You can determine the angle between them by using the vector property:

(x1 * x2) + (y1 * y2)
cos( φ ) = -------------------------------------------
sqrt( x1^2 + y1^2 ) * sqrt( x2^2 + y2^2 )

Then depending on the angle you can say if it is left or right.
exception (323)
But using trigonometry, square roots and division is suboptimal. *Very* suboptimal for several reasons.
Hint: you can determine the solution of every geometric predicate with a determinant. Often this solution is not optimal either, but at least it is "only" a polynomial.
In this case, substract one point from the others and write the result in a 2x2 matrix. Then the sign of the determinant is your answer. For problems regarding rounding errors see here: http://www.cplusplus.com/forum/articles/3827/
There is an exact, lazy adaptive solution for this problem by Shewchuk: http://www.cs.berkeley.edu/~jrs/papers/robustr.ps
It describes floating point expansions and then uses them for the orientation tests in 2D and 3D.
If you use the formula I described above, you get (for the points p,q,r to avoid confusing them with the coordinates):
(qx-px)*(ry-py)-(rx-px)*(qy-py)
...which has an error bound of 22b-48, where b is the number of bits you require to store ceil(max(px,py,qx,qy,rx,ry)).
So, basically, if you don't want Shewchuk's algorithm (it's a bit complicated), you can do:
1
2
3
4
5
6
7
8
 
int bitsize = calc_bitsize(ceil(max(px,py,qx,qy,rx,ry)));
double errorbound = exp(2, bitsize-48);
double result = (qx-px)*(ry-py)-(rx-px)*(qy-py);
if(result < errorbound)
  return true;
if(result > errorbound)
  return false;
return <calculate result using exact arithmetic, or, if you want a potentially wrong answer, return "p, q and r are on a line">

This, of course, implies that neither overflow nor underflow occurs. And if you know the largest number at compile time, you should perhaps precompute it (it is usually better to use, say, 10 for b instead of computing it every time, getting, say, 4-8 depending on the case (if you don't make a case study about degenerate cases))
kanesoban (39)
How does the (qx-px)*(ry-py)-(rx-px)*(qy-py) formula change, if i say that the y coordinate 'grows' as you go down ?
xabnu (72)
Why don't you just take the 'slope' between two points, then between two other points? You can determine enough from that, if you're only working with 2 dimensions.

If you are afraid of an undefined slope, you can always check for divide-by-zero, and have a way to distinguish between positive and negative infinity, if your computer won't do that for you.

1
2
3
4
5
6
7
8
9
10
11
12
 

double findslope(point p0, point p1){
  return (p0.y - p1.y) / (p0.x - p1.y);
}

void whichwaydiditurn(point p0, point p1, point p2){
  if( findslope(p0,p1) > findslope(p0,p2) )
    printf("you turned right!\n");
  else
    printf("you turned left or went straight!\n");
}

Last edited on
exception (323)
How does the (qx-px)*(ry-py)-(rx-px)*(qy-py) formula change, if i say that the y coordinate 'grows' as you go down ?

If I understand correctly, you want to invert the y axis. Well, the sign will switch in this case. Since it does so for all cases, it isn't an issue.

Why don't you just take the 'slope' between two points, then between two other points?

You can, but it is (a) not a general solution (extra check, which is expensive) and it requires (b) a division. You should try to (you always can!) avoid divisions in geometric decision making (like this predicate), since they are introduce (I) a large error and are (II) slow (and(III) many error bound calculations can't deal with them). As I stated, determinants can always be used for geometric predicates.
kanesoban (39)
If I understand correctly, you want to invert the y axis. Well, the sign will switch in this case. Since it does so for all cases, it isn't an issue.


You mean (rx-px)*(qy-py) - (qx-px)*(ry-py) instead of (qx-px)*(ry-py)-(rx-px)*(qy-py) ? I tried that, but it doesn't seem to work.
exception (323)
Then your code doesn't seem to be correct. Do you understand the solution? If not, perhaps you should alctually use one of the alternatives described here. Or learn enough about linear algebra and projective geometry to understand the given formula. (Hint: the formula is an abbreviation for

| px qx rx |
det(| py qy ry |)
| 1 1 1 |

(or the same transposed). Due to the rules applicable to determinants, this is equal to the formula I gave (where one column is substructed from the other two, so that the bottom row becomes 0), but requires more calculations -> has a larger error bound)
Topic archived. No new replies allowed.