File indexing completed on 2024-04-21 14:46:23
0001 /* 0002 ** Author: Eric Veach, July 1994. 0003 ** 0004 */ 0005 0006 #include "gluos.h" 0007 #include <assert.h> 0008 #include "mesh.h" 0009 #include "geom.h" 0010 0011 int __gl_vertLeq(GLUvertex *u, GLUvertex *v) 0012 { 0013 /* Returns TRUE if u is lexicographically <= v. */ 0014 0015 return VertLeq(u, v); 0016 } 0017 0018 GLdouble __gl_edgeEval(GLUvertex *u, GLUvertex *v, GLUvertex *w) 0019 { 0020 /* Given three vertices u,v,w such that VertLeq(u,v) && VertLeq(v,w), 0021 * evaluates the t-coord of the edge uw at the s-coord of the vertex v. 0022 * Returns v->t - (uw)(v->s), ie. the signed distance from uw to v. 0023 * If uw is vertical (and thus passes thru v), the result is zero. 0024 * 0025 * The calculation is extremely accurate and stable, even when v 0026 * is very close to u or w. In particular if we set v->t = 0 and 0027 * let r be the negated result (this evaluates (uw)(v->s)), then 0028 * r is guaranteed to satisfy MIN(u->t,w->t) <= r <= MAX(u->t,w->t). 0029 */ 0030 GLdouble gapL, gapR; 0031 0032 assert(VertLeq(u, v) && VertLeq(v, w)); 0033 0034 gapL = v->s - u->s; 0035 gapR = w->s - v->s; 0036 0037 if (gapL + gapR > 0) 0038 { 0039 if (gapL < gapR) 0040 { 0041 return (v->t - u->t) + (u->t - w->t) * (gapL / (gapL + gapR)); 0042 } 0043 else 0044 { 0045 return (v->t - w->t) + (w->t - u->t) * (gapR / (gapL + gapR)); 0046 } 0047 } 0048 /* vertical line */ 0049 return 0; 0050 } 0051 0052 GLdouble __gl_edgeSign(GLUvertex *u, GLUvertex *v, GLUvertex *w) 0053 { 0054 /* Returns a number whose sign matches EdgeEval(u,v,w) but which 0055 * is cheaper to evaluate. Returns > 0, == 0 , or < 0 0056 * as v is above, on, or below the edge uw. 0057 */ 0058 GLdouble gapL, gapR; 0059 0060 assert(VertLeq(u, v) && VertLeq(v, w)); 0061 0062 gapL = v->s - u->s; 0063 gapR = w->s - v->s; 0064 0065 if (gapL + gapR > 0) 0066 { 0067 return (v->t - w->t) * gapL + (v->t - u->t) * gapR; 0068 } 0069 /* vertical line */ 0070 return 0; 0071 } 0072 0073 /*********************************************************************** 0074 * Define versions of EdgeSign, EdgeEval with s and t transposed. 0075 */ 0076 0077 GLdouble __gl_transEval(GLUvertex *u, GLUvertex *v, GLUvertex *w) 0078 { 0079 /* Given three vertices u,v,w such that TransLeq(u,v) && TransLeq(v,w), 0080 * evaluates the t-coord of the edge uw at the s-coord of the vertex v. 0081 * Returns v->s - (uw)(v->t), ie. the signed distance from uw to v. 0082 * If uw is vertical (and thus passes thru v), the result is zero. 0083 * 0084 * The calculation is extremely accurate and stable, even when v 0085 * is very close to u or w. In particular if we set v->s = 0 and 0086 * let r be the negated result (this evaluates (uw)(v->t)), then 0087 * r is guaranteed to satisfy MIN(u->s,w->s) <= r <= MAX(u->s,w->s). 0088 */ 0089 GLdouble gapL, gapR; 0090 0091 assert(TransLeq(u, v) && TransLeq(v, w)); 0092 0093 gapL = v->t - u->t; 0094 gapR = w->t - v->t; 0095 0096 if (gapL + gapR > 0) 0097 { 0098 if (gapL < gapR) 0099 { 0100 return (v->s - u->s) + (u->s - w->s) * (gapL / (gapL + gapR)); 0101 } 0102 else 0103 { 0104 return (v->s - w->s) + (w->s - u->s) * (gapR / (gapL + gapR)); 0105 } 0106 } 0107 /* vertical line */ 0108 return 0; 0109 } 0110 0111 GLdouble __gl_transSign(GLUvertex *u, GLUvertex *v, GLUvertex *w) 0112 { 0113 /* Returns a number whose sign matches TransEval(u,v,w) but which 0114 * is cheaper to evaluate. Returns > 0, == 0 , or < 0 0115 * as v is above, on, or below the edge uw. 0116 */ 0117 GLdouble gapL, gapR; 0118 0119 assert(TransLeq(u, v) && TransLeq(v, w)); 0120 0121 gapL = v->t - u->t; 0122 gapR = w->t - v->t; 0123 0124 if (gapL + gapR > 0) 0125 { 0126 return (v->s - w->s) * gapL + (v->s - u->s) * gapR; 0127 } 0128 /* vertical line */ 0129 return 0; 0130 } 0131 0132 int __gl_vertCCW(GLUvertex *u, GLUvertex *v, GLUvertex *w) 0133 { 0134 /* For almost-degenerate situations, the results are not reliable. 0135 * Unless the floating-point arithmetic can be performed without 0136 * rounding errors, *any* implementation will give incorrect results 0137 * on some degenerate inputs, so the client must have some way to 0138 * handle this situation. 0139 */ 0140 return (u->s * (v->t - w->t) + v->s * (w->t - u->t) + w->s * (u->t - v->t)) >= 0; 0141 } 0142 0143 /* Given parameters a,x,b,y returns the value (b*x+a*y)/(a+b), 0144 * or (x+y)/2 if a==b==0. It requires that a,b >= 0, and enforces 0145 * this in the rare case that one argument is slightly negative. 0146 * The implementation is extremely stable numerically. 0147 * In particular it guarantees that the result r satisfies 0148 * MIN(x,y) <= r <= MAX(x,y), and the results are very accurate 0149 * even when a and b differ greatly in magnitude. 0150 */ 0151 #define RealInterpolate(a, x, b, y) \ 0152 (a = (a < 0) ? 0 : a, b = (b < 0) ? 0 : b, \ 0153 ((a <= b) ? ((b == 0) ? ((x + y) / 2) : (x + (y - x) * (a / (a + b)))) : (y + (x - y) * (b / (a + b))))) 0154 0155 #ifndef FOR_TRITE_TEST_PROGRAM 0156 #define Interpolate(a, x, b, y) RealInterpolate(a, x, b, y) 0157 #else 0158 0159 /* Claim: the ONLY property the sweep algorithm relies on is that 0160 * MIN(x,y) <= r <= MAX(x,y). This is a nasty way to test that. 0161 */ 0162 #include <stdlib.h> 0163 extern int RandomInterpolate; 0164 0165 GLdouble Interpolate(GLdouble a, GLdouble x, GLdouble b, GLdouble y) 0166 { 0167 printf("*********************%d\n", RandomInterpolate); 0168 if (RandomInterpolate) 0169 { 0170 a = 1.2 * drand48() - 0.1; 0171 a = (a < 0) ? 0 : ((a > 1) ? 1 : a); 0172 b = 1.0 - a; 0173 } 0174 return RealInterpolate(a, x, b, y); 0175 } 0176 0177 #endif 0178 0179 #define Swap(a, b) \ 0180 do \ 0181 { \ 0182 GLUvertex *t = a; \ 0183 a = b; \ 0184 b = t; \ 0185 } while (0) 0186 0187 void __gl_edgeIntersect(GLUvertex *o1, GLUvertex *d1, GLUvertex *o2, GLUvertex *d2, GLUvertex *v) 0188 /* Given edges (o1,d1) and (o2,d2), compute their point of intersection. 0189 * The computed point is guaranteed to lie in the intersection of the 0190 * bounding rectangles defined by each edge. 0191 */ 0192 { 0193 GLdouble z1, z2; 0194 0195 /* This is certainly not the most efficient way to find the intersection 0196 * of two line segments, but it is very numerically stable. 0197 * 0198 * Strategy: find the two middle vertices in the VertLeq ordering, 0199 * and interpolate the intersection s-value from these. Then repeat 0200 * using the TransLeq ordering to find the intersection t-value. 0201 */ 0202 0203 if (!VertLeq(o1, d1)) 0204 { 0205 Swap(o1, d1); 0206 } 0207 if (!VertLeq(o2, d2)) 0208 { 0209 Swap(o2, d2); 0210 } 0211 if (!VertLeq(o1, o2)) 0212 { 0213 Swap(o1, o2); 0214 Swap(d1, d2); 0215 } 0216 0217 if (!VertLeq(o2, d1)) 0218 { 0219 /* Technically, no intersection -- do our best */ 0220 v->s = (o2->s + d1->s) / 2; 0221 } 0222 else if (VertLeq(d1, d2)) 0223 { 0224 /* Interpolate between o2 and d1 */ 0225 z1 = EdgeEval(o1, o2, d1); 0226 z2 = EdgeEval(o2, d1, d2); 0227 if (z1 + z2 < 0) 0228 { 0229 z1 = -z1; 0230 z2 = -z2; 0231 } 0232 v->s = Interpolate(z1, o2->s, z2, d1->s); 0233 } 0234 else 0235 { 0236 /* Interpolate between o2 and d2 */ 0237 z1 = EdgeSign(o1, o2, d1); 0238 z2 = -EdgeSign(o1, d2, d1); 0239 if (z1 + z2 < 0) 0240 { 0241 z1 = -z1; 0242 z2 = -z2; 0243 } 0244 v->s = Interpolate(z1, o2->s, z2, d2->s); 0245 } 0246 0247 /* Now repeat the process for t */ 0248 0249 if (!TransLeq(o1, d1)) 0250 { 0251 Swap(o1, d1); 0252 } 0253 if (!TransLeq(o2, d2)) 0254 { 0255 Swap(o2, d2); 0256 } 0257 if (!TransLeq(o1, o2)) 0258 { 0259 Swap(o1, o2); 0260 Swap(d1, d2); 0261 } 0262 0263 if (!TransLeq(o2, d1)) 0264 { 0265 /* Technically, no intersection -- do our best */ 0266 v->t = (o2->t + d1->t) / 2; 0267 } 0268 else if (TransLeq(d1, d2)) 0269 { 0270 /* Interpolate between o2 and d1 */ 0271 z1 = TransEval(o1, o2, d1); 0272 z2 = TransEval(o2, d1, d2); 0273 if (z1 + z2 < 0) 0274 { 0275 z1 = -z1; 0276 z2 = -z2; 0277 } 0278 v->t = Interpolate(z1, o2->t, z2, d1->t); 0279 } 0280 else 0281 { 0282 /* Interpolate between o2 and d2 */ 0283 z1 = TransSign(o1, o2, d1); 0284 z2 = -TransSign(o1, d2, d1); 0285 if (z1 + z2 < 0) 0286 { 0287 z1 = -z1; 0288 z2 = -z2; 0289 } 0290 v->t = Interpolate(z1, o2->t, z2, d2->t); 0291 } 0292 }