#include "GraphicsGems.h" #undef ON /* Comparison macros */ #define LEFT -1 /* value to left of (less than) another */ #define ON 0 /* two values equal */ #define RIGHT 1 /* value to right of (greater than) another */ typedef struct endpoint { /* An a priori endpoint */ short x, y; } ENDPOINT; typedef struct segment { /* A priori line segment with integer endpoints */ ENDPOINT first, last; /* defined to be ordered from "first" to "last" */ } SEGMENT; typedef struct param { /* Parameterized description of an intersection along a SEGMENT */ long num, denom; } PARAM; typedef struct subsegment { SEGMENT apriori; /* The a priori segment this subsegment falls on */ PARAM ParOne, ParTwo;/* Parameterized description of intersection points */ } SUBSEGMENT; typedef struct intpoint { /* Intersection point returned by SegIntersect */ PARAM par[2]; /* par[0] is on the first segment, par[1] on the second */ long a, b, c, d; /* storing these allows fast computation for direction of crossing */ } INTPOINT; /* All a priori endpoints must have coordinates falling in this range */ #define PointRange(X) (((X) > -16384) && ((X) < 16383)) #define SHORTMASK 0xffff /* Used by mult64 */ /* TRUE iff A and B have same signs. */ #define SAME_SIGNS(A, B) (((long)((unsigned long)A ^ (unsigned long)B)) >= 0) /* Return the max value, storing the minimum value in min */ #define maxmin(x1, x2, min) (x1 >= x2 ? (min = x2, x1) : (min = x1, x2)) /* ********************************************************************* Below are two utility functions for implementing exact intersection calculation: SubsegIntersect and SideOfPoint. SubsegIntersect uses SegIntersect to do the bulk of its work. ********************************************************************* */ /* ********************************************************************* Compute the intersection points between two a priori line segments. Return 0 if no intersection, 1 if intersects in a point, and 2 if intersecting lines are colinear. Parameter values for intersection point are returned in ipt. Entry: s1, s2: the line segments Exit: ipt: the intersection point. ********************************************************************* */ int SegIntersect(SEGMENT *s1, SEGMENT *s2,INTPOINT *ipt) { long a, b, c, d, tdet, sdet, det; /* parameter calculation variables */ short max1, max2, min1, min2; /* bounding box check variables */ ENDPOINT p1, p2, q1, q2; /* dereference a priori endpoints */ p1 = s1->first; p2 = s1->last; q1 = s2->first; q2 = s2->last; /* First make the bounding box test. */ max1 = maxmin(p1.x, p2.x, min1); max2 = maxmin(q1.x, q2.x, min2); if((max1 < min2) || (min1 > max2)) return(0); /* no intersection */ max1 = maxmin(p1.y, p2.y, min1); max2 = maxmin(q1.y, q2.y, min2); if((max1 < min2) || (min1 > max2)) return(0); /* no intersection */ /* See if the endpoints of the second segment lie on the opposite sides of the first. If not, return 0. */ a = (long)(q1.x - p1.x) * (long)(p2.y - p1.y) - (long)(q1.y - p1.y) * (long)(p2.x - p1.x); b = (long)(q2.x - p1.x) * (long)(p2.y - p1.y) - (long)(q2.y - p1.y) * (long)(p2.x - p1.x); if(a!=0 && b!=0 && SAME_SIGNS(a, b)) return(0); /* See if the endpoints of the first segment lie on the opposite sides of the second. If not, return 0. */ c = (long)(p1.x - q1.x) * (long)(q2.y - q1.y) - (long)(p1.y - q1.y) * (long)(q2.x - q1.x); d = (long)(p2.x - q1.x) * (long)(q2.y - q1.y) - (long)(p2.y - q1.y) * (long)(q2.x - q1.x); if(c!=0 && d!=0 && SAME_SIGNS(c, d) ) return(0); /* At this point each segment meets the line of the other. */ det = a - b; if(det == 0) return(2); /* The segments are colinear. Determining colinear intersection parameters would be tedious and not instructive. */ /* The segments intersect since each segment crosses the other's line; however, since the lines are not parallel, either a or b is not zero. Similarly either c or d is not zero. */ tdet = -c; sdet = a; if(det < 0) /* The denominator of the parameter must be positive. */ { det = -det; sdet = -sdet; tdet = -tdet; } ipt->a = a; ipt->b = b; ipt->c = c; ipt->d = d; ipt->par[0].num = tdet; ipt->par[0].denom = det; ipt->par[1].num = sdet; ipt->par[1].denom = det; return(1); } /* ********************************************************************* Returns the number of intersection points (0, 1 or 2) between line subsegments ss1 and ss2. Note that 2 intersection points would represent the endpoints of overlap between two colinear line segments. If there is a single intersection point, it is returned in ipt. Entry: ss1, ss2: the line subsegments. Exit: ipt: the intersection point (if exactly one). ********************************************************************* */ int SubsegIntersect(SUBSEGMENT *ss1, SUBSEGMENT *ss2, INTPOINT *ipt) { int i; /* Number of intersection points */ /* intersect a priori segments */ i = SegIntersect(&ss1->apriori, &ss2->apriori, ipt); if(i != 1) return(i); /* either no intersection or colinear */ /* one intersection point - check if it falls outside of subsegments */ if(LessThan(&(ipt->par[0]), &(ss1->ParOne)) || GreaterThan(&ipt->par[0], &(ss1->ParTwo))) return(0); /* A priori segments intersect, but subsegments do not */ if(LessThan(&ipt->par[1], &ss2->ParOne) || GreaterThan(&ipt->par[1], &ss2->ParTwo)) return(0); /* A priori segments intersect, but subsegments do not */ return(1); } /* ********************************************************************* 64 bit multiplication function. Multiply two long integers in1 and in2, returning the result in out. ********************************************************************* */ void mult64(long in1, long in2, unsigned long out[2]) { unsigned short *x, *y, *z; unsigned long temp; x = (unsigned short *) &in1; y = (unsigned short *) &in2; z = (unsigned short *) out; temp = x[0] * y[0]; z[1] = temp >> 16; z[0] = temp & SHORTMASK; temp = x[0] * y[1]; z[2] = temp >> 16; z[1] += temp & SHORTMASK; z[2] += z[1] >>16; z[1] = z[1] & SHORTMASK; temp = x[1] * y[0]; z[2] += temp >> 16; z[1] += temp & SHORTMASK; z[2] += z[1] >>16; z[1] = z[1] & SHORTMASK; z[3] = z[2] >>16; z[2] = z[2] & SHORTMASK; temp = x[1] * y[1]; z[3] += temp >> 16; z[2] += temp & SHORTMASK; z[3] += z[2] >>16; z[2] = z[2] & SHORTMASK; } /* ********************************************************************* Comparison primitive to test par1 <= par2. Return LEFT if par1 < par2; return ON if par1 = par2; return RIGHT if par1 > par2. ********************************************************************* */ int CompPrim(PARAM *par1, PARAM *par2) { unsigned long r1[2],r2[2]; mult64(par1->num, par2->denom, r1); mult64(par2->num, par1->denom, r2); if(r1[1] != r2[1]) { if(r1[1] < r2[1]) return(LEFT); else return(RIGHT);} if(r1[0] != r2[0]) { if(r1[0] < r2[0]) return(LEFT); else return(RIGHT);} return(ON); } /* ********************************************************************* Helper function for all parameter comparisons. Returns LEFT if par1 < par2, ON if par1 = par2, and RIGHT otherwise. Denominators must be positive. ********************************************************************* */ int CompHelp(PARAM *par1, PARAM *par2) { PARAM tpar1, tpar2; tpar1 = *par1; tpar2 = *par2; if(tpar1.num < 0) { if(tpar2.num >= 0) return(LEFT); tpar1.num = -tpar1.num; tpar2.num = -tpar2.num; return(CompPrim(&tpar2, &tpar1)); } if(tpar2.num < 0) return(RIGHT); return(CompPrim(&tpar1, &tpar2)); } /* ********************************************************************* Returns TRUE if par1 < par2 and FALSE otherwise; ********************************************************************* */ boolean LessThan(PARAM *par1, PARAM *par2) { double x1, x2; x1 = ((double)par1->num) * par2->denom; x2 = ((double)par2->num) * par1->denom; if(x1 != x2) if(x1 < x2) return(TRUE); else return(FALSE); return(CompHelp(par1, par2) == LEFT); } /* ********************************************************************* Returns TRUE if par1 > par2 and FALSE otherwise; ********************************************************************* */ boolean GreaterThan(PARAM *par1, PARAM *par2) { double x1, x2; x1 = ((double)par1->num) * par2->denom; x2 = ((double)par2->num) * par1->denom; if(x1 != x2) if(x1 > x2) return(TRUE); else return(FALSE); return(CompHelp(par1, par2) == RIGHT); } /* ********************************************************************* Returns TRUE if par1 = par2 and FALSE otherwise; ********************************************************************* */ boolean Equal(PARAM *par1, PARAM *par2) { double x1, x2; x1 = ((double)par1->num) * par2->denom; x2 = ((double)par2->num) * par1->denom; if(x1 != x2) return(FALSE); return(CompHelp(par1, par2) == ON); } /* ********************************************************************* Determine if a given point is to the left of, right of, or on segment s1. The point is defined as the position (inpt) along a line segment (in). Returns LEFT, RIGHT or ON respectively. Entry: in: line segment defining intersection point inpt: intersection point parameter along "in" s1: segment to check on which side of ********************************************************************* */ int SideOfPoint(SEGMENT *in, PARAM *inpar, SEGMENT *s1) { PARAM par; /* build a fraction for use in comparison against inpar */ long delx, dely; /* multiplications must be done as longs */ short *x1, *x2, *y1, *y2, *x3, *x4, *y3, *y4; /* alias names */ x1 = &s1->first.x; y1 = &s1->first.y; x2 = &s1->last.x; y2 = &s1->last.y; x3 = &in->first.x; y3 = &in->first.y; x4 = &in->last.x; y4 = &in->last.y; delx = *x2 - *x1; dely = *y2 - *y1; par.denom = (*y4 - *y3) * delx - (*x4 - *x3) * dely; par.num = (*x3 - *x1) * dely - (*y3 - *y1) * delx; if(par.denom > 0) return(CompHelp(&par, inpar)); else if(par.denom < 0) /* switch signs and compare */ { par.denom = -par.denom; par.num = -par.num; return(CompHelp(inpar, &par)); } else if(par.num > 0) return(RIGHT); else if(par.num < 0) return(LEFT); else return(ON); }