Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!samsung!usc!ucsd!nosc!crash!jcs
From: jcs@crash.cts.com (John Schultz)
Newsgroups: comp.graphics
Subject: Re: Algorithm needed to determine if a point is inside a polygon
Summary: Here's some code to try.
Keywords: polygon
Message-ID: <3752@crash.cts.com>
Date: 30 Jul 90 17:40:50 GMT
References: <1990Jul27.201939.17816@cbnewsh.att.com] <46093@brunix.UUCP> <6951@helios.TAMU.EDU>
Distribution: na
Organization: Crash TimeSharing, El Cajon, CA
Lines: 242


/*
  Here's a line/convex polygon intersection routine. I haven't fully tested
it, so if you spot a bug, let me know.
This code is set up for 32 bit integer fixed point, but can easily be
modified for floats by removing the shifts and declaring floats.
Another trick is to rotate and translate the *ray* (two points) and not
the polygons/normals.  The ray transform matrix is the inverse of the
polygon (local) transform. See the Siggraph article for more info.

  John

*/
/* Collision.c */
/* Created by John C. Schultz 18-July-90 */
/* Idea from Siggraph Computer Graphics Vol 21, Number 4, July '87,
   page 126 */

#include <math.h>
#include <stdio.h>
#define sgn(x) ((x)>0 ? 1: ((x)<0 ? -1: (0)))
#define FBITS 14

typedef struct point3d  {long defx,defy,defz,x,y,z;} point3d;

/* Prototypes */

short linepolyint(point3d * ip,point3d * b,point3d * d,point3d * n,
                  point3d * poly,short count);

/* Functions */

short linepolyint(point3d * ip,point3d * b,point3d * d,point3d * n,
                  point3d * poly,short count) {

long k,t;   /* k = Plane constant, t = parametric var */ 
long sn;   /* sign max normal */
long absx,absy,absz; /* abs x,y,z of normal */
long nd; /* N dot D */
long beta; /* Result of cross product */
short i,i1;
point3d p,q; /* temps */

/* Compute:
  t = parametric variable
  k = plane constant.
  N = normal (plane equation).
  B,D = line base, direction.
  
                  k - (N dot B)
             t = -----------
                    (N dot D)
*/

/* Compute plane contant k */
  k = poly[0].x*n->x + poly[0].y*n->y + poly[0].z*n->z; /* Keep upshifted */

  nd = (n->x*d->x + n->y*d->y + n->z*d->z) >> FBITS;

  if (nd) {

/* Compute parametric variable */

    t = (k - (n->x*b->x + n->y*b->y + n->z*b->z)) / nd;

  } else {

    return 0; /* Line is parallel to plane */

  } /* if nd */

/* Now see if point is within polygon */
  absx = abs(n->x);
  absy = abs(n->y);
  absz = abs(n->z);

/* Find max of normal{x,y,z}, do specific case */

  if (absx > absy) {
    if (absx > absz) { /************ X CASE ************/

      sn = sgn(n->x);

/* Compute only intersection points being checked */
      ip->y = b->y + ((t*d->y) >> FBITS);
      ip->z = b->z + ((t*d->z) >> FBITS);

      for (i=0; i<count; i++) {

        if ( (i1 = i+1) >= count) i1 = 0; /* Modulo count */     

        p.y = poly[i1].y - poly[i].y;
        p.z = poly[i1].z - poly[i].z;

        q.y = ip->y - poly[i].y;
        q.z = ip->z - poly[i].z;

        beta = (p.y*q.z - p.z*q.y); /* cross product: makes an x */

/* beta == 0 when on polygon boundary */

        if ((beta) && (sgn(beta) != sn)) return 0;
/* The above line may need to be modified to implement:
             0 <= beta <= 1, if n->x is positive and
            -1 <= beta <= 0, if n->x is negative.
 I simplified it for speed.
 The above code occurs twice more below.
*/

      }/* for i */

/* Success, now compute last point */
      ip->x = b->x + ((t*d->x) >> FBITS);

    } else {

      goto zmax; /* For speed and non-redunancy */

    } /* if x */

  } else {

    if (absy > absz) { /************ Y CASE ************/

      sn = sgn(n->y);

/* Compute only intersection points being checked */
      ip->x = b->x + ((t*d->x) >> FBITS);
      ip->z = b->z + ((t*d->z) >> FBITS);

      for (i=0; i<count; i++) {

        if ( (i1 = i+1) >= count) i1 = 0;      

        p.x = poly[i1].x - poly[i].x;
        p.z = poly[i1].z - poly[i].z;

        q.x = ip->x - poly[i].x;
        q.z = ip->z - poly[i].z;

        beta = (p.z*q.x - p.x*q.z); /* cross product: makes a y */

/* beta == 0 when on polygon boundary */

        if ((beta) && (sgn(beta) != sn)) return 0;

      }/* for i */

/* Success, now compute last point */
      ip->y = b->y + ((t*d->y) >> FBITS);

    } else {
zmax:                  /************ Z CASE ************/
      sn = sgn(n->z);

/* Compute only intersection points being checked */
      ip->x = b->x + ((t*d->x) >> FBITS);
      ip->y = b->y + ((t*d->y) >> FBITS);

      for (i=0; i<count; i++) {

        if ( (i1 = i+1) >= count) i1 = 0;      

        p.x = poly[i1].x - poly[i].x;
        p.y = poly[i1].y - poly[i].y;

        q.x = ip->x - poly[i].x;
        q.y = ip->y - poly[i].y;

        beta = (p.x*q.y - p.y*q.x); /* cross product: makes a z */

/* beta == 0 when on polygon boundary */

        if ((beta) && (sgn(beta) != sn)) return 0;

      }/* for i */

/* Success, now compute last point */
      ip->z = b->z + ((t*d->z) >> FBITS);

    } /* if y/z */
  } /* if x/y */

  return 1; /* Line intersected polygon */

} /* end raypolyint */

void main(long argc,char ** argv) {
point3d b,d,n,ip,p1,p2;
point3d poly[3];
short count = 3;
short hit;

/* Initialize intersection point */
  ip.x = 0;
  ip.y = 0;
  ip.z = 0;

/* Start of line */
  p1.x = 0;
  p1.y = -1000;
  p1.z =-9992;

/* End of line */
  p2.x = 0;
  p2.y = 1000;
  p2.z = -10002;

/* Normal */
  n.x = 0;
  n.y = -((1<<FBITS)-1); /* -1.0 */
  n.z = 0;

/* Line base is p1 */
  b = p1;

/* Line direction is p2 - p1 */
  d.x = p2.x - p1.x;
  d.y = p2.y - p1.y;
  d.z = p2.z - p1.z;

/* Define polygon */
  poly[0].x =  0;
  poly[0].y =  500;
  poly[0].z =  10000;

  poly[2].x =  10000;
  poly[2].y =  500;
  poly[2].z = -10000;

  poly[1].x = -10000;
  poly[1].y =  500;
  poly[1].z = -10000;

/* Test the code */
  hit = linepolyint(&ip,&b,&d,&n,poly,count);

  printf("hit %d, x %d y %d z %d\n",hit,ip.x,ip.y,ip.z);

} /* end main */

/* end collision.c */