Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!sun-barr!newstop!exodus!texsun!csccat!ncmicro!chris
From: chris@ncmicro.lonestar.org (Chris Arps)
Newsgroups: comp.graphics
Subject: Re: Intersections
Summary: I have a 'C' double precision algorithm that I use in NC
	 work.  I hope that it will suffice as it is used for
	 accuracy not for integer or display type calculations.
Message-ID: <194@ncmicro.lonestar.org>
Date: 12 Oct 90 15:41:30 GMT
References: <31960@nigel.ee.udel.edu> <10669@goofy.Apple.COM>
Organization: NC Microproducts, Richardson, Tx.
Lines: 70



/* Intersection of two 2d lines - Chris Arps, NC Microproducts Inc. */


/* $Func  - interltl() intersect two lines */
/* $Input - x1,y1 start point of line #1   */
/*          x2,y2 end point of line #1     */
/*          x3,y3 start point of line #2   */
/*          x4,y4 end point of line #2     */
/* $Output - x5,y5 intersection point of two lines */
/* $Returns - < 0 no intersection due to parallel,too small, etc. */
/*                    x5,y5 are NOT computed  			  */
/*              0 x5,y5 computed , line segments do not touch */
/*              1 x5,y5 computed , line segments touch */

#define ABSZERO 1.0E-10    /* tolerance for zero, works on PC/Sun */

int interltl(x1,y1,x2,y2,x3,y3,x4,y4,x5,y5)
double x1,y1,x2,y2,x3,y3,x4,y4;
double *x5,*y5;
{
   extern double fabs(),sqrt();

   double a1,b1,c1,a2,b2,c2;
   double x,y,t0,t1;

   a1 = y1 - y2;
   b1 = x2 - x1;
   c1 = x1*y2 - x2*y1;

   a2 = y3 - y4;
   b2 = x4 - x3;
   c2 = x3*y4 - x4*y3;

   y = a1*b2 - a2*b1;
   x = a1 - a2;

   if ( fabs(y) < ABSZERO) {  /* parallel */
      x = x1*y2 + x3*y1 + x2*y3 - x3*y2 - x1*y3 - x2*y1;
      if ( fabs(x) < ABSZERO ) return -2;  /* collinear */
      return -1;
   }

   y = (a2*c1 - a1*c2) / y;

   if ( fabs(x) < ABSZERO) return -3; /* both horizontal */

   x = (b2*y - b1*y + c2 - c1) / x;

   if ( fabs(x2-x1) < ABSZERO ) { 
      if (fabs(y2-y1) < ABSZERO ) return -4; /* line too small */
      else t0 = (y - y1) / (y2 - y1);
   } 
   else t0 = (x - x1) / (x2 - x1);

   if ( fabs(x4-x3) < ABSZERO ) {
      if (fabs(y4-y3) < ABSZERO ) return -4; /* line too small */
      else t1 = (y - y3) / (y4 - y3);
   } 
   else t1 = (x - x3) / (x4 - x3);

   *x5 = x;
   *y5 = y;
  
   if ( t0 < ABSZERO || t0 > (1.0 - ABSZERO) || 
        t1 < ABSZERO || t1 > (1.0 - ABSZERO) ) return 0;

   return 1;
}