1 /*
2 ** Author: Eric Veach, July 1994.
3 **
4 */
5
6 #include "gluos.h"
7 #include <assert.h>
8 #include "mesh.h"
9 #include "geom.h"
10
__gl_vertLeq(GLUvertex * u,GLUvertex * v)11 int __gl_vertLeq(GLUvertex *u, GLUvertex *v)
12 {
13 /* Returns TRUE if u is lexicographically <= v. */
14
15 return VertLeq(u, v);
16 }
17
__gl_edgeEval(GLUvertex * u,GLUvertex * v,GLUvertex * w)18 GLdouble __gl_edgeEval(GLUvertex *u, GLUvertex *v, GLUvertex *w)
19 {
20 /* Given three vertices u,v,w such that VertLeq(u,v) && VertLeq(v,w),
21 * evaluates the t-coord of the edge uw at the s-coord of the vertex v.
22 * Returns v->t - (uw)(v->s), ie. the signed distance from uw to v.
23 * If uw is vertical (and thus passes thru v), the result is zero.
24 *
25 * The calculation is extremely accurate and stable, even when v
26 * is very close to u or w. In particular if we set v->t = 0 and
27 * let r be the negated result (this evaluates (uw)(v->s)), then
28 * r is guaranteed to satisfy MIN(u->t,w->t) <= r <= MAX(u->t,w->t).
29 */
30 GLdouble gapL, gapR;
31
32 assert(VertLeq(u, v) && VertLeq(v, w));
33
34 gapL = v->s - u->s;
35 gapR = w->s - v->s;
36
37 if (gapL + gapR > 0)
38 {
39 if (gapL < gapR)
40 {
41 return (v->t - u->t) + (u->t - w->t) * (gapL / (gapL + gapR));
42 }
43 else
44 {
45 return (v->t - w->t) + (w->t - u->t) * (gapR / (gapL + gapR));
46 }
47 }
48 /* vertical line */
49 return 0;
50 }
51
__gl_edgeSign(GLUvertex * u,GLUvertex * v,GLUvertex * w)52 GLdouble __gl_edgeSign(GLUvertex *u, GLUvertex *v, GLUvertex *w)
53 {
54 /* Returns a number whose sign matches EdgeEval(u,v,w) but which
55 * is cheaper to evaluate. Returns > 0, == 0 , or < 0
56 * as v is above, on, or below the edge uw.
57 */
58 GLdouble gapL, gapR;
59
60 assert(VertLeq(u, v) && VertLeq(v, w));
61
62 gapL = v->s - u->s;
63 gapR = w->s - v->s;
64
65 if (gapL + gapR > 0)
66 {
67 return (v->t - w->t) * gapL + (v->t - u->t) * gapR;
68 }
69 /* vertical line */
70 return 0;
71 }
72
73 /***********************************************************************
74 * Define versions of EdgeSign, EdgeEval with s and t transposed.
75 */
76
__gl_transEval(GLUvertex * u,GLUvertex * v,GLUvertex * w)77 GLdouble __gl_transEval(GLUvertex *u, GLUvertex *v, GLUvertex *w)
78 {
79 /* Given three vertices u,v,w such that TransLeq(u,v) && TransLeq(v,w),
80 * evaluates the t-coord of the edge uw at the s-coord of the vertex v.
81 * Returns v->s - (uw)(v->t), ie. the signed distance from uw to v.
82 * If uw is vertical (and thus passes thru v), the result is zero.
83 *
84 * The calculation is extremely accurate and stable, even when v
85 * is very close to u or w. In particular if we set v->s = 0 and
86 * let r be the negated result (this evaluates (uw)(v->t)), then
87 * r is guaranteed to satisfy MIN(u->s,w->s) <= r <= MAX(u->s,w->s).
88 */
89 GLdouble gapL, gapR;
90
91 assert(TransLeq(u, v) && TransLeq(v, w));
92
93 gapL = v->t - u->t;
94 gapR = w->t - v->t;
95
96 if (gapL + gapR > 0)
97 {
98 if (gapL < gapR)
99 {
100 return (v->s - u->s) + (u->s - w->s) * (gapL / (gapL + gapR));
101 }
102 else
103 {
104 return (v->s - w->s) + (w->s - u->s) * (gapR / (gapL + gapR));
105 }
106 }
107 /* vertical line */
108 return 0;
109 }
110
__gl_transSign(GLUvertex * u,GLUvertex * v,GLUvertex * w)111 GLdouble __gl_transSign(GLUvertex *u, GLUvertex *v, GLUvertex *w)
112 {
113 /* Returns a number whose sign matches TransEval(u,v,w) but which
114 * is cheaper to evaluate. Returns > 0, == 0 , or < 0
115 * as v is above, on, or below the edge uw.
116 */
117 GLdouble gapL, gapR;
118
119 assert(TransLeq(u, v) && TransLeq(v, w));
120
121 gapL = v->t - u->t;
122 gapR = w->t - v->t;
123
124 if (gapL + gapR > 0)
125 {
126 return (v->s - w->s) * gapL + (v->s - u->s) * gapR;
127 }
128 /* vertical line */
129 return 0;
130 }
131
__gl_vertCCW(GLUvertex * u,GLUvertex * v,GLUvertex * w)132 int __gl_vertCCW(GLUvertex *u, GLUvertex *v, GLUvertex *w)
133 {
134 /* For almost-degenerate situations, the results are not reliable.
135 * Unless the floating-point arithmetic can be performed without
136 * rounding errors, *any* implementation will give incorrect results
137 * on some degenerate inputs, so the client must have some way to
138 * handle this situation.
139 */
140 return (u->s * (v->t - w->t) + v->s * (w->t - u->t) + w->s * (u->t - v->t)) >= 0;
141 }
142
143 /* Given parameters a,x,b,y returns the value (b*x+a*y)/(a+b),
144 * or (x+y)/2 if a==b==0. It requires that a,b >= 0, and enforces
145 * this in the rare case that one argument is slightly negative.
146 * The implementation is extremely stable numerically.
147 * In particular it guarantees that the result r satisfies
148 * MIN(x,y) <= r <= MAX(x,y), and the results are very accurate
149 * even when a and b differ greatly in magnitude.
150 */
151 #define RealInterpolate(a, x, b, y) \
152 (a = (a < 0) ? 0 : a, b = (b < 0) ? 0 : b, \
153 ((a <= b) ? ((b == 0) ? ((x + y) / 2) : (x + (y - x) * (a / (a + b)))) : (y + (x - y) * (b / (a + b)))))
154
155 #ifndef FOR_TRITE_TEST_PROGRAM
156 #define Interpolate(a, x, b, y) RealInterpolate(a, x, b, y)
157 #else
158
159 /* Claim: the ONLY property the sweep algorithm relies on is that
160 * MIN(x,y) <= r <= MAX(x,y). This is a nasty way to test that.
161 */
162 #include <stdlib.h>
163 extern int RandomInterpolate;
164
Interpolate(GLdouble a,GLdouble x,GLdouble b,GLdouble y)165 GLdouble Interpolate(GLdouble a, GLdouble x, GLdouble b, GLdouble y)
166 {
167 printf("*********************%d\n", RandomInterpolate);
168 if (RandomInterpolate)
169 {
170 a = 1.2 * drand48() - 0.1;
171 a = (a < 0) ? 0 : ((a > 1) ? 1 : a);
172 b = 1.0 - a;
173 }
174 return RealInterpolate(a, x, b, y);
175 }
176
177 #endif
178
179 #define Swap(a, b) \
180 do \
181 { \
182 GLUvertex *t = a; \
183 a = b; \
184 b = t; \
185 } while (0)
186
__gl_edgeIntersect(GLUvertex * o1,GLUvertex * d1,GLUvertex * o2,GLUvertex * d2,GLUvertex * v)187 void __gl_edgeIntersect(GLUvertex *o1, GLUvertex *d1, GLUvertex *o2, GLUvertex *d2, GLUvertex *v)
188 /* Given edges (o1,d1) and (o2,d2), compute their point of intersection.
189 * The computed point is guaranteed to lie in the intersection of the
190 * bounding rectangles defined by each edge.
191 */
192 {
193 GLdouble z1, z2;
194
195 /* This is certainly not the most efficient way to find the intersection
196 * of two line segments, but it is very numerically stable.
197 *
198 * Strategy: find the two middle vertices in the VertLeq ordering,
199 * and interpolate the intersection s-value from these. Then repeat
200 * using the TransLeq ordering to find the intersection t-value.
201 */
202
203 if (!VertLeq(o1, d1))
204 {
205 Swap(o1, d1);
206 }
207 if (!VertLeq(o2, d2))
208 {
209 Swap(o2, d2);
210 }
211 if (!VertLeq(o1, o2))
212 {
213 Swap(o1, o2);
214 Swap(d1, d2);
215 }
216
217 if (!VertLeq(o2, d1))
218 {
219 /* Technically, no intersection -- do our best */
220 v->s = (o2->s + d1->s) / 2;
221 }
222 else if (VertLeq(d1, d2))
223 {
224 /* Interpolate between o2 and d1 */
225 z1 = EdgeEval(o1, o2, d1);
226 z2 = EdgeEval(o2, d1, d2);
227 if (z1 + z2 < 0)
228 {
229 z1 = -z1;
230 z2 = -z2;
231 }
232 v->s = Interpolate(z1, o2->s, z2, d1->s);
233 }
234 else
235 {
236 /* Interpolate between o2 and d2 */
237 z1 = EdgeSign(o1, o2, d1);
238 z2 = -EdgeSign(o1, d2, d1);
239 if (z1 + z2 < 0)
240 {
241 z1 = -z1;
242 z2 = -z2;
243 }
244 v->s = Interpolate(z1, o2->s, z2, d2->s);
245 }
246
247 /* Now repeat the process for t */
248
249 if (!TransLeq(o1, d1))
250 {
251 Swap(o1, d1);
252 }
253 if (!TransLeq(o2, d2))
254 {
255 Swap(o2, d2);
256 }
257 if (!TransLeq(o1, o2))
258 {
259 Swap(o1, o2);
260 Swap(d1, d2);
261 }
262
263 if (!TransLeq(o2, d1))
264 {
265 /* Technically, no intersection -- do our best */
266 v->t = (o2->t + d1->t) / 2;
267 }
268 else if (TransLeq(d1, d2))
269 {
270 /* Interpolate between o2 and d1 */
271 z1 = TransEval(o1, o2, d1);
272 z2 = TransEval(o2, d1, d2);
273 if (z1 + z2 < 0)
274 {
275 z1 = -z1;
276 z2 = -z2;
277 }
278 v->t = Interpolate(z1, o2->t, z2, d1->t);
279 }
280 else
281 {
282 /* Interpolate between o2 and d2 */
283 z1 = TransSign(o1, o2, d1);
284 z2 = -TransSign(o1, d2, d1);
285 if (z1 + z2 < 0)
286 {
287 z1 = -z1;
288 z2 = -z2;
289 }
290 v->t = Interpolate(z1, o2->t, z2, d2->t);
291 }
292 }
293