From: sloan@cis.uab.edu (Kenneth Sloan)
Subject: Re: Point inside a triangle
Organization: CIS, University of Alabama at Birmingham
Date: Tue, 17 Nov 1992 04:12:16 GMT

In article <1e9c3cINN4hi@uniwa.uwa.edu.au> watson@maths.uwa.oz.au (David Watson) writes:
>...
>A fail-proof method is to compute the barycentric coordinates.  For a
>triangle {(x1,y1), (x2,y2), (x3,y3)} and some point (x0,y0), calculate
>
>b0 = (x2 - x1) * (y3 - y1) - (x3 - x1) * (y2 - y1)
>b1 = ((x2 - x0) * (y3 - y0) - (x3 - x0) * (y2 - y0)) / b0 
>b2 = ((x3 - x0) * (y1 - y0) - (x1 - x0) * (y3 - y0)) / b0
>b3 = ((x1 - x0) * (y2 - y0) - (x2 - x0) * (y1 - y0)) / b0 
>
>Then if b1, b2, and b3 are all > 0, (x0,y0) is strictly inside the
>triangle; if bi = 0 and the other two coordinates are positive, (x0,y0)
>lies on the edge opposite (xi,yi); if bi and bj = 0, (x0,y0) lies on
>(xk,yk); if bi < 0, (x0,y0) lies outside the edge opposite (xi,yi);
>if all three coordinates are negative, something else is wrong.  This
>method does not depend on the cyclic order of the vertices.

So far, so good.  If the above tests are all you require, I might take
the trouble to get the order "right" so that the areas are positive, and
eliminate the "/ b0"'s.

Even better is to replace the last line by:

b3 = 1.0 - b1 - b2

Oh yes - be sure to check that:

       |b0| > EPSILON

or all bets are off.

Here's a piece of code that demonstrates the sort of trouble you can get
into if you aren't careful.  Bugs to me, please.  Please don't bother to
"optimize" it - all of the code which uses this function has been
extensively profiled - trust me, this is not even close to the hot spot.

And, beware - there are a few assumptions hidden in the code which
happen to be enforced by the surrounding environment.  I suppose I ought
to fix that. Just lazy, I guess.

================================================================
 /*
   InsideTriangle decides if a point P is Inside, or OnEdge
   of the triangle defined by A, B, C.

   If P is on an edge, V1 and V2 identify that edge
 */
void InsideTriangle (A, B, C, P, Inside, OnEdge, V1, V2)
 int A, B, C, P;
 int *Inside, *OnEdge;
 int *V1, *V2;
 {
  double ax, ay, bx, by, cx, cy, apx, apy, bpx, bpy, cpx, cpy;
  double cCROSSap, bCROSScp, aCROSSbp, DOT;
  int OnEdgeAB, OnEdgeBC, OnEdgeCA;

  *Inside = false; *OnEdge = false;
  OnEdgeAB = false; OnEdgeBC = false; OnEdgeCA = false;
  ax = Points[C].x - Points[B].x;  ay = Points[C].y - Points[B].y;
  bx = Points[A].x - Points[C].x;  by = Points[A].y - Points[C].y;
  cx = Points[B].x - Points[A].x;  cy = Points[B].y - Points[A].y;
  apx= Points[P].x - Points[A].x;  apy= Points[P].y - Points[A].y;
  bpx= Points[P].x - Points[B].x;  bpy= Points[P].y - Points[B].y;
  cpx= Points[P].x - Points[C].x;  cpy= Points[P].y - Points[C].y;

  aCROSSbp = ax*bpy - ay*bpx;
  cCROSSap = cx*apy - cy*apx;
  bCROSScp = bx*cpy - by*cpx;

  *Inside = ((aCROSSbp > 0.0) && (bCROSScp > 0.0) && (cCROSSap > 0.0));

  if (!(*Inside))
   {
    if ((EPSILON > aCROSSbp) && (aCROSSbp > -EPSILON))
     {
      DOT = (ax*bpx)+(ay*bpy);
      OnEdgeBC = ((0.0 <= DOT) && (DOT <= (ax*ax)+(ay*ay)));
     }

    if ((EPSILON > bCROSScp) && (bCROSScp > -EPSILON))
     {
      DOT  = (bx*cpx)+(by*cpy);
      OnEdgeCA = ((0.0 <= DOT) && (DOT <= (bx*bx)+(by*by)));
     }

    if ((EPSILON > cCROSSap) && (cCROSSap > -EPSILON))
     {
      DOT = (cx*apx)+(cy*apy);
      OnEdgeAB = ((0.0 <= DOT) && (DOT <= (cx*cx)+(cy*cy)));
     }
    if      (OnEdgeAB) { *V1 = A; *V2 = B; *OnEdge = true; }
    else if (OnEdgeBC) { *V1 = B; *V2 = C; *OnEdge = true; }
    else if (OnEdgeCA) { *V1 = C; *V2 = A; *OnEdge = true; }
   }
 } /* InsideTriangle */
================================================================
-- 
Kenneth Sloan                   Computer and Information Sciences
sloan@cis.uab.edu               University of Alabama at Birmingham
(205) 934-2213                  115A Campbell Hall, UAB Station 
(205) 934-5473 FAX              Birmingham, AL 35294-1170