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
llinter(x1,y1,x2,y2,z1,w1,z2,w2,alpha,beta,collinear)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
lcinter(x1,y1,x2,y2,x0,y0,r,alpha1,theta1,alpha2,theta2)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
ccinter(x0,y0,r0,x1,y1,r1,theta1,phi1,theta2,phi2)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