1 #ifndef lint 2 static char *sccsid ="inter.c (CWI) 1.1 85/03/01"; 3 #endif 4 #include "ideal.h" 5 #include "y.tab.h" 6 7 boolean llinter (x1,y1, x2,y2, z1,w1, z2,w2, alpha, beta, collinear) 8 /* if this function returns TRUE, 9 /* then alpha[(x1,y1),(x2,y2)] = beta[(z1,w1),(z2,w2)] */ 10 float x1,y1, x2,y2, z1,w1, z2,w2; 11 float *alpha; 12 float *beta; 13 boolean *collinear; 14 { 15 float A1, B1, C1, A2, B2, C2, D, x, y; 16 dprintf "(%f,%f) -- (%f,%f)\n", x1,y1, x2,y2); 17 dprintf "(%f,%f) -- (%f,%f)\n", z1,w1, z2,w2); 18 A1 = y1 - y2; 19 B1 = x2 - x1; 20 C1 = -B1*y1 - A1*x1; 21 A2 = w1 - w2; 22 B2 = z2 - z1; 23 C2 = -B2*w1 - A2*z1; 24 D = A1*B2 - A2*B1; 25 if (fabs(D) < EPSILON) { 26 *collinear = arecollinear(x1,y1,x2,y2,z1,w1); 27 dprintf "%s\n", (*collinear)?"coincident":"disjoint"); 28 return (FALSE); 29 } 30 *collinear = FALSE; 31 x = (B1*C2 - B2*C1)/D; 32 y = (A2*C1 - A1*C2)/D; 33 if (fabs(x2 - x1) > EPSILON) { 34 *alpha = (x - x1)/(x2 - x1); 35 } else if (fabs(y2 - y1) > EPSILON) { 36 *alpha = (y - y1)/(y2 - y1); 37 } else fprintf (stderr, "ideal: llinter: can't happen\n"); 38 if (fabs(z2 - z1) > EPSILON) { 39 *beta = (x - z1)/(z2 - z1); 40 } else if (fabs(w2 - w1) > EPSILON) { 41 *beta = (y - w1)/(w2 - w1); 42 } else fprintf (stderr, "ideal: llinter: can't happen\n"); 43 dprintf "intersection alpha = %f; beta = %f\n", *alpha, *beta); 44 return (TRUE); 45 } /* llinter */ 46 47 boolean lcinter (x1,y1, x2,y2, x0,y0, r, alpha1,theta1, alpha2,theta2) 48 /* if this function returns TRUE, 49 /* then alpha1[(x1,y1),(x2,y2)] = (x0,y0) + r*cis(theta1) 50 /* and alpha2[(x1,y1),(x2,y2)] = (x0,y0) + r*cis(theta2) */ 51 float x1,y1, x2,y2, x0,y0, r; 52 float *alpha1; 53 float *theta1; 54 float *alpha2; 55 float *theta2; 56 { 57 float dx1, dx2, dy1, dy2; 58 float A, B, C, D; 59 60 dprintf "intersection parameters:\n"); 61 dprintf "%f, %f -- %f, %f\n", x1, y1, x2, y2); 62 dprintf "%f, %f (%f)\n", x0, y0, r); 63 r = fabs(r); 64 dx1 = x1 - x0; 65 dx2 = x2 - x1; 66 dy1 = y1 - y0; 67 dy2 = y2 - y1; 68 A = dx2*dx2 + dy2*dy2; 69 dprintf "A=%f\n", A); 70 if (A < EPSILON) { 71 if (fabs (hypot (dx1, dy1) - r) < EPSILON) { 72 *alpha1 = *alpha2 = 0.0; 73 *theta1 = atan2 (dy1, dx1); 74 *theta2 = *theta1 = rprin (*theta1); 75 dprintf "alpha1 = alpha2 = %f theta1 = theta2 = %f\n", 76 *alpha1, *theta1); 77 return (TRUE); 78 } 79 else 80 return (FALSE); 81 } 82 B = 2*(dx1*dx2 + dy1*dy2); 83 C = dx1*dx1 + dy1*dy1 - r*r; 84 D = B*B - 4*A*C; 85 dprintf "B=%f C=%f D=%f\n", B, C, D); 86 if (D < -EPSILON) 87 return (FALSE); 88 if (fabs(D) < EPSILON) 89 D = 0.0; 90 D = sqrt(D); 91 *alpha1 = (-B + D)/(2.0*A); 92 *theta1 = rprin (atan2 (dy1 + *alpha1*dy2, dx1 + *alpha1*dx2)); 93 *alpha2 = (-B - D)/(2.0*A); 94 *theta2 = rprin (atan2 (dy1 + *alpha2*dy2, dx1 + *alpha2*dx2)); 95 dprintf "intersection alpha1 = %f, theta1 = %f\n", *alpha1, *theta1); 96 dprintf "intersection alpha2 = %f, theta2 = %f\n", *alpha2, *theta2); 97 return (TRUE); 98 } 99 100 boolean ccinter (x0,y0,r0, x1,y1,r1, theta1,phi1, theta2,phi2) 101 /* if this function returns TRUE, 102 /* then (x0,y0) + r0*cis(theta1) = (x1,y1) + r1*cis(phi1) 103 /* and (x0,y0) + r0*cis(theta2) = (x1,y1) + r1*cis(phi2) */ 104 float x0,y0,r0; 105 float x1,y1,r1; 106 float *theta1; 107 float *phi1; 108 float *theta2; 109 float *phi2; 110 { 111 float xcoeff, ycoeff, const; 112 float u1, v1, u2, v2; 113 boolean lncrc; 114 115 dprintf "intersection parameters\n"); 116 dprintf "%f %f (%f)\n", x0, y0, r0); 117 dprintf "%f %f (%f)\n", x1, y1, r1); 118 119 r0 = fabs(r0); 120 r1 = fabs(r1); 121 xcoeff = 2*(x1 - x0); 122 ycoeff = 2*(y1 - y0); 123 const = r0*r0 - x0*x0 - y0*y0 - r1*r1 + x1*x1 + y1*y1; 124 if (fabs(xcoeff) < EPSILON && fabs(ycoeff) < EPSILON) 125 return (FALSE); 126 if (fabs(xcoeff) < EPSILON) { 127 u1 = 0.0; 128 u2 = 1.0; 129 v1 = v2 = const/ycoeff; 130 } else if (fabs(ycoeff) < EPSILON) { 131 v1 = 0.0; 132 v2 = 1.0; 133 u1 = u2 = const/xcoeff; 134 } else if (fabs(const) < EPSILON) { 135 u1 = 0.0; 136 v1 = 0.0; 137 u2 = 1.0; 138 v2 = (const - 1.0/xcoeff)/ycoeff; 139 } else { 140 u1 = 0.0; 141 v1 = const/ycoeff; 142 u2 = const/xcoeff; 143 v2 = 0.0; 144 } 145 lncrc = lcinter (u1,v1, u2,v2, x1,y1,r1, theta1,phi1, theta2,phi2); 146 if (lncrc) { 147 *phi1 = rprin (*phi1); 148 *phi2 = rprin (*phi2); 149 *theta1 = atan2 (y1 + r1*sin(*phi1) - y0, x1 + r1*cos(*phi1) - x0); 150 *theta2 = atan2 (y1 + r1*sin(*phi2) - y0, x1 + r1*cos(*phi2) - x0); 151 *theta1 = rprin (*theta1); 152 *theta2 = rprin (*theta2); 153 dprintf "intersection theta1 = %f phi1 = %f\n", *theta1, *phi1); 154 dprintf "intersection theta2 = %f phi2 = %f\n", *theta2, *phi2); 155 return (TRUE); 156 } else 157 return (FALSE); 158 } 159