home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / gems / gemsiii / polyintr.c < prev    next >
C/C++ Source or Header  |  1992-03-16  |  11KB  |  300 lines
  1.  
  2.  
  3.  
  4.  
  5.  
  6.  
  7.  
  8.  
  9.  
  10.  
  11.  
  12.  
  13.  
  14.  
  15.  
  16.  
  17.  
  18.  
  19. #include "GraphicsGems.h"
  20. #undef ON
  21.  
  22. /* Comparison macros */
  23. #define LEFT  -1  /* value to left of (less than) another */
  24. #define ON     0  /* two values equal */
  25. #define RIGHT  1  /* value to right of (greater than) another */
  26.  
  27. typedef struct endpoint
  28. { /* An a priori endpoint */
  29.   short x, y;
  30. } ENDPOINT;
  31.  
  32. typedef struct segment
  33. { /* A priori line segment with integer endpoints */
  34.   ENDPOINT first, last; /* defined to be ordered from "first" to "last" */
  35. } SEGMENT;
  36.  
  37. typedef struct param
  38. { /* Parameterized description of an intersection along a SEGMENT */
  39.   long num, denom;
  40. } PARAM;
  41.  
  42. typedef struct subsegment
  43. {
  44.   SEGMENT apriori;     /* The a priori segment this subsegment falls on */
  45.   PARAM ParOne, ParTwo;/* Parameterized description of intersection points */
  46. } SUBSEGMENT;
  47.  
  48. typedef struct intpoint
  49. { /* Intersection point returned by SegIntersect */
  50.   PARAM par[2];  /* par[0] is on the first segment, par[1] on the second */
  51.   long a, b, c, d; /* storing these allows fast computation for direction of
  52.               crossing */
  53. } INTPOINT;
  54.  
  55. /* All a priori endpoints must have coordinates falling in this range */
  56. #define PointRange(X)  (((X) > -16384) && ((X) < 16383))
  57.  
  58. #define SHORTMASK 0xffff  /* Used by mult64 */
  59.  
  60. /* TRUE iff A and B have same signs. */
  61. #define SAME_SIGNS(A, B) (((long)((unsigned long)A ^ (unsigned long)B)) >= 0)
  62.  
  63. /* Return the max value, storing the minimum value in min */
  64. #define  maxmin(x1, x2, min) (x1 >= x2 ? (min = x2, x1) : (min = x1, x2))
  65.  
  66. /* *********************************************************************
  67.    Below are two utility functions for implementing exact intersection
  68.    calculation: SubsegIntersect and SideOfPoint.  SubsegIntersect uses
  69.    SegIntersect to do the bulk of its work.
  70. ********************************************************************* */
  71.  
  72. /* *********************************************************************
  73.    Compute the intersection points between two a priori line segments.
  74.    Return 0 if no intersection, 1 if intersects in a point,
  75.    and 2 if intersecting lines are colinear. Parameter values for
  76.    intersection point are returned in ipt.
  77.  
  78.    Entry: s1, s2: the line segments
  79.    Exit:  ipt: the intersection point.
  80. ********************************************************************* */
  81. int SegIntersect(SEGMENT *s1, SEGMENT *s2,INTPOINT *ipt)
  82. {
  83.   long a, b, c, d, tdet, sdet, det;  /* parameter calculation variables */
  84.   short max1, max2, min1, min2; /* bounding box check variables */
  85.   ENDPOINT p1, p2, q1, q2; /* dereference a priori endpoints */
  86.  
  87.   p1 = s1->first;    p2 = s1->last;
  88.   q1 = s2->first;    q2 = s2->last;
  89.  
  90.   /*  First make the bounding box test. */
  91.   max1 = maxmin(p1.x, p2.x, min1);
  92.   max2 = maxmin(q1.x, q2.x, min2);
  93.   if((max1 < min2) || (min1 > max2)) return(0); /* no intersection */
  94.   max1 = maxmin(p1.y, p2.y, min1);
  95.   max2 = maxmin(q1.y, q2.y, min2);
  96.   if((max1 < min2) || (min1 > max2)) return(0); /* no intersection */
  97.  
  98.   /* See if the endpoints of the second segment lie on the opposite
  99.      sides of the first.  If not, return 0. */
  100.   a = (long)(q1.x - p1.x) * (long)(p2.y - p1.y) -
  101.       (long)(q1.y - p1.y) * (long)(p2.x - p1.x);
  102.   b = (long)(q2.x - p1.x) * (long)(p2.y - p1.y) -
  103.       (long)(q2.y - p1.y) * (long)(p2.x - p1.x);
  104.   if(a!=0 && b!=0 && SAME_SIGNS(a, b)) return(0);
  105.  
  106.   /* See if the endpoints of the first segment lie on the opposite
  107.      sides of the second.  If not, return 0.  */
  108.   c = (long)(p1.x - q1.x) * (long)(q2.y - q1.y) -
  109.       (long)(p1.y - q1.y) * (long)(q2.x - q1.x);
  110.   d = (long)(p2.x - q1.x) * (long)(q2.y - q1.y) -
  111.       (long)(p2.y - q1.y) * (long)(q2.x - q1.x);
  112.   if(c!=0 && d!=0 && SAME_SIGNS(c, d) ) return(0);
  113.  
  114.   /* At this point each segment meets the line of the other. */
  115.   det = a - b;
  116.   if(det == 0) return(2); /* The segments are colinear.  Determining
  117.      colinear intersection parameters would be tedious and not instructive. */
  118.  
  119.   /* The segments intersect since each segment crosses the other's line;
  120.      however, since the lines are not parallel, either a or b is not
  121.      zero.  Similarly either c or d is not zero. */
  122.   tdet = -c;   sdet =   a;
  123.   if(det < 0)  /* The denominator of the parameter must be positive. */
  124.     {  det = -det;  sdet = -sdet;  tdet = -tdet;  }
  125.   ipt->a = a; ipt->b = b;
  126.   ipt->c = c; ipt->d = d;
  127.   ipt->par[0].num = tdet;
  128.   ipt->par[0].denom = det;
  129.   ipt->par[1].num = sdet;
  130.   ipt->par[1].denom = det;
  131.   return(1);
  132. }
  133.  
  134. /* *********************************************************************
  135.    Returns the number of intersection points (0, 1 or 2) between line
  136.    subsegments ss1 and ss2.  Note that 2 intersection points would
  137.    represent the endpoints of overlap between two colinear line segments.
  138.    If there is a single intersection point, it is returned in ipt.
  139.  
  140.    Entry: ss1, ss2: the line subsegments.
  141.    Exit:  ipt: the intersection point (if exactly one).
  142. ********************************************************************* */
  143. int SubsegIntersect(SUBSEGMENT *ss1, SUBSEGMENT *ss2, INTPOINT *ipt)
  144. {
  145.   int i;  /* Number of intersection points */
  146.  
  147.   /* intersect a priori segments */
  148.   i = SegIntersect(&ss1->apriori, &ss2->apriori, ipt);
  149.   if(i != 1) return(i);  /* either no intersection or colinear */
  150.   /* one intersection point - check if it falls outside of subsegments */
  151.   if(LessThan(&(ipt->par[0]), &(ss1->ParOne)) ||
  152.      GreaterThan(&ipt->par[0], &(ss1->ParTwo)))
  153.     return(0);  /* A priori segments intersect, but subsegements do not */
  154.   if(LessThan(&ipt->par[1], &ss2->ParOne) ||
  155.      GreaterThan(&ipt->par[1], &ss2->ParTwo))
  156.     return(0);  /* A priori segments intersect, but subsegements do not */
  157.   return(1);
  158. }
  159.  
  160.  
  161. /* *********************************************************************
  162.    64 bit multiplication function.
  163.    Multiply two long integers in1 and in2, returning the result in out.
  164. ********************************************************************* */
  165. void mult64(long in1, long in2, unsigned long out[2])
  166. {
  167.   unsigned short *x, *y, *z;
  168.   unsigned long temp;
  169.  
  170.   x = (unsigned short *) &in1;  y = (unsigned short *) &in2;
  171.   z = (unsigned short *) out;
  172.   temp = x[0] * y[0];
  173.   z[1] = temp >> 16;   z[0] = temp & SHORTMASK;
  174.   temp = x[0] * y[1];
  175.   z[2] = temp >> 16;   z[1] += temp & SHORTMASK;
  176.   z[2] += z[1] >>16;   z[1] = z[1] & SHORTMASK;
  177.   temp = x[1] * y[0];
  178.   z[2] += temp >> 16;  z[1] += temp & SHORTMASK;
  179.   z[2] += z[1] >>16;  z[1] = z[1] & SHORTMASK;
  180.   z[3] = z[2] >>16;  z[2] = z[2] & SHORTMASK;
  181.   temp = x[1] * y[1];
  182.   z[3] += temp >> 16;   z[2] += temp & SHORTMASK;
  183.   z[3] += z[2] >>16;   z[2] = z[2] & SHORTMASK;
  184. }
  185.  
  186.  
  187. /* *********************************************************************
  188. Comparison primitive to test par1 <= par2.
  189. Return LEFT if par1 < par2; return ON if par1 = par2;
  190. return RIGHT if par1 > par2.
  191. ********************************************************************* */
  192. int CompPrim(PARAM *par1, PARAM *par2)
  193. {
  194.   unsigned long r1[2],r2[2];
  195.  
  196.   mult64(par1->num, par2->denom, r1);
  197.   mult64(par2->num, par1->denom, r2);
  198.   if(r1[1] != r2[1]) { if(r1[1] < r2[1]) return(LEFT); else return(RIGHT);}
  199.   if(r1[0] != r2[0]) { if(r1[0] < r2[0]) return(LEFT); else return(RIGHT);}
  200.   return(ON);
  201. }
  202.  
  203.  
  204. /* *********************************************************************
  205.   Helper function for all parameter comparisons.
  206.   Returns LEFT if par1 < par2, ON if par1 = par2, and RIGHT otherwise.
  207.   Denominators must be positive.
  208. ********************************************************************* */
  209. int CompHelp(PARAM *par1, PARAM *par2)
  210. {
  211.   PARAM tpar1, tpar2;
  212.  
  213.   tpar1 = *par1;   tpar2 = *par2;
  214.   if(tpar1.num < 0)
  215.      {
  216.        if(tpar2.num >= 0) return(LEFT);
  217.        tpar1.num = -tpar1.num;  tpar2.num = -tpar2.num;
  218.        return(CompPrim(&tpar2, &tpar1));
  219.      }
  220.   if(tpar2.num < 0) return(RIGHT);
  221.   return(CompPrim(&tpar1, &tpar2));
  222. }
  223.  
  224.  
  225. /* *********************************************************************
  226.   Returns TRUE if par1 < par2 and FALSE otherwise;
  227. ********************************************************************* */
  228. boolean LessThan(PARAM *par1, PARAM *par2)
  229. {
  230.   double x1, x2;
  231.  
  232.   x1 = ((double)par1->num) * par2->denom;
  233.   x2 = ((double)par2->num) * par1->denom;
  234.   if(x1 != x2)
  235.     if(x1 < x2) return(TRUE);
  236.     else return(FALSE);
  237.   return(CompHelp(par1, par2) == LEFT);
  238. }
  239.  
  240.  
  241. /* *********************************************************************
  242.   Returns TRUE if par1 > par2 and FALSE otherwise;
  243. ********************************************************************* */
  244. boolean GreaterThan(PARAM *par1, PARAM *par2)
  245. {
  246.   double x1, x2;
  247.  
  248.   x1 = ((double)par1->num) * par2->denom;
  249.   x2 = ((double)par2->num) * par1->denom;
  250.   if(x1 != x2)
  251.     if(x1 > x2) return(TRUE);
  252.     else return(FALSE);
  253.   return(CompHelp(par1, par2) == RIGHT);
  254. }
  255.  
  256.  
  257. /* *********************************************************************
  258.   Returns TRUE if par1 = par2 and FALSE otherwise;
  259. ********************************************************************* */
  260. boolean Equal(PARAM *par1, PARAM *par2)
  261. {
  262.   double x1, x2;
  263.  
  264.   x1 = ((double)par1->num) * par2->denom;
  265.   x2 = ((double)par2->num) * par1->denom;
  266.   if(x1 != x2) return(FALSE);
  267.   return(CompHelp(par1, par2) == ON);
  268. }
  269.  
  270.  
  271. /* *********************************************************************
  272.    Determine if a given point is to the left of, right of,
  273.    or on segment s1.  The point is defined as the position (inpt)
  274.    along a line segment (in).  Returns LEFT, RIGHT or ON respectively.
  275.  
  276.    Entry:
  277.      in: line segment defining intersection point
  278.      inpt: intersection point parameter along "in"
  279.      s1: segment to check on which side of
  280. ********************************************************************* */
  281. int SideOfPoint(SEGMENT *in, PARAM *inpar, SEGMENT *s1)
  282. {
  283.   PARAM par;  /* build a fraction for use in comparison against inpar */
  284.   long delx, dely;  /* multiplications must be done as longs */
  285.   short *x1, *x2, *y1, *y2, *x3, *x4, *y3, *y4; /* alias names */
  286.  
  287.   x1 = &s1->first.x; y1 = &s1->first.y; x2 = &s1->last.x; y2 = &s1->last.y;
  288.   x3 = &in->first.x; y3 = &in->first.y; x4 = &in->last.x; y4 = &in->last.y;
  289.   delx = *x2 - *x1;  dely = *y2 - *y1;
  290.   par.denom =  (*y4 - *y3) * delx - (*x4 - *x3) * dely;
  291.   par.num   =  (*x3 - *x1) * dely - (*y3 - *y1) * delx;
  292.   if(par.denom > 0) return(CompHelp(&par, inpar));
  293.   else if(par.denom < 0) /* switch signs and compare */
  294.     { par.denom = -par.denom; par.num = -par.num;
  295.       return(CompHelp(inpar, &par)); }
  296.   else if(par.num > 0) return(RIGHT);
  297.   else if(par.num < 0) return(LEFT);
  298.   else return(ON);
  299. }
  300.