1 /* 2 * SGI FREE SOFTWARE LICENSE B (Version 2.0, Sept. 18, 2008) 3 * Copyright (C) 1991-2000 Silicon Graphics, Inc. All Rights Reserved. 4 * 5 * Permission is hereby granted, free of charge, to any person obtaining a 6 * copy of this software and associated documentation files (the "Software"), 7 * to deal in the Software without restriction, including without limitation 8 * the rights to use, copy, modify, merge, publish, distribute, sublicense, 9 * and/or sell copies of the Software, and to permit persons to whom the 10 * Software is furnished to do so, subject to the following conditions: 11 * 12 * The above copyright notice including the dates of first publication and 13 * either this permission notice or a reference to 14 * http://oss.sgi.com/projects/FreeB/ 15 * shall be included in all copies or substantial portions of the Software. 16 * 17 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS 18 * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, 19 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 20 * SILICON GRAPHICS, INC. BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, 21 * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF 22 * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE 23 * SOFTWARE. 24 * 25 * Except as contained in this notice, the name of Silicon Graphics, Inc. 26 * shall not be used in advertising or otherwise to promote the sale, use or 27 * other dealings in this Software without prior written authorization from 28 * Silicon Graphics, Inc. 29 */ 30 /* 31 ** Author: Eric Veach, July 1994. 32 ** 33 */ 34 35 #include "gluos.h" 36 #include <assert.h> 37 //#include "mesh.h" 38 #include "geom.h" 39 40 int __gl_vertLeq( GLUvertex *u, GLUvertex *v ) 41 { 42 /* Returns TRUE if u is lexicographically <= v. */ 43 44 return VertLeq( u, v ); 45 } 46 47 GLdouble __gl_edgeEval( GLUvertex *u, GLUvertex *v, GLUvertex *w ) 48 { 49 /* Given three vertices u,v,w such that VertLeq(u,v) && VertLeq(v,w), 50 * evaluates the t-coord of the edge uw at the s-coord of the vertex v. 51 * Returns v->t - (uw)(v->s), ie. the signed distance from uw to v. 52 * If uw is vertical (and thus passes thru v), the result is zero. 53 * 54 * The calculation is extremely accurate and stable, even when v 55 * is very close to u or w. In particular if we set v->t = 0 and 56 * let r be the negated result (this evaluates (uw)(v->s)), then 57 * r is guaranteed to satisfy MIN(u->t,w->t) <= r <= MAX(u->t,w->t). 58 */ 59 GLdouble gapL, gapR; 60 61 assert( VertLeq( u, v ) && VertLeq( v, w )); 62 63 gapL = v->s - u->s; 64 gapR = w->s - v->s; 65 66 if( gapL + gapR > 0 ) { 67 if( gapL < gapR ) { 68 return (v->t - u->t) + (u->t - w->t) * (gapL / (gapL + gapR)); 69 } else { 70 return (v->t - w->t) + (w->t - u->t) * (gapR / (gapL + gapR)); 71 } 72 } 73 /* vertical line */ 74 return 0; 75 } 76 77 GLdouble __gl_edgeSign( GLUvertex *u, GLUvertex *v, GLUvertex *w ) 78 { 79 /* Returns a number whose sign matches EdgeEval(u,v,w) but which 80 * is cheaper to evaluate. Returns > 0, == 0 , or < 0 81 * as v is above, on, or below the edge uw. 82 */ 83 GLdouble gapL, gapR; 84 85 assert( VertLeq( u, v ) && VertLeq( v, w )); 86 87 gapL = v->s - u->s; 88 gapR = w->s - v->s; 89 90 if( gapL + gapR > 0 ) { 91 return (v->t - w->t) * gapL + (v->t - u->t) * gapR; 92 } 93 /* vertical line */ 94 return 0; 95 } 96 97 98 /*********************************************************************** 99 * Define versions of EdgeSign, EdgeEval with s and t transposed. 100 */ 101 102 GLdouble __gl_transEval( GLUvertex *u, GLUvertex *v, GLUvertex *w ) 103 { 104 /* Given three vertices u,v,w such that TransLeq(u,v) && TransLeq(v,w), 105 * evaluates the t-coord of the edge uw at the s-coord of the vertex v. 106 * Returns v->s - (uw)(v->t), ie. the signed distance from uw to v. 107 * If uw is vertical (and thus passes thru v), the result is zero. 108 * 109 * The calculation is extremely accurate and stable, even when v 110 * is very close to u or w. In particular if we set v->s = 0 and 111 * let r be the negated result (this evaluates (uw)(v->t)), then 112 * r is guaranteed to satisfy MIN(u->s,w->s) <= r <= MAX(u->s,w->s). 113 */ 114 GLdouble gapL, gapR; 115 116 assert( TransLeq( u, v ) && TransLeq( v, w )); 117 118 gapL = v->t - u->t; 119 gapR = w->t - v->t; 120 121 if( gapL + gapR > 0 ) { 122 if( gapL < gapR ) { 123 return (v->s - u->s) + (u->s - w->s) * (gapL / (gapL + gapR)); 124 } else { 125 return (v->s - w->s) + (w->s - u->s) * (gapR / (gapL + gapR)); 126 } 127 } 128 /* vertical line */ 129 return 0; 130 } 131 132 GLdouble __gl_transSign( GLUvertex *u, GLUvertex *v, GLUvertex *w ) 133 { 134 /* Returns a number whose sign matches TransEval(u,v,w) but which 135 * is cheaper to evaluate. Returns > 0, == 0 , or < 0 136 * as v is above, on, or below the edge uw. 137 */ 138 GLdouble gapL, gapR; 139 140 assert( TransLeq( u, v ) && TransLeq( v, w )); 141 142 gapL = v->t - u->t; 143 gapR = w->t - v->t; 144 145 if( gapL + gapR > 0 ) { 146 return (v->s - w->s) * gapL + (v->s - u->s) * gapR; 147 } 148 /* vertical line */ 149 return 0; 150 } 151 152 153 int __gl_vertCCW( GLUvertex *u, GLUvertex *v, GLUvertex *w ) 154 { 155 /* For almost-degenerate situations, the results are not reliable. 156 * Unless the floating-point arithmetic can be performed without 157 * rounding errors, *any* implementation will give incorrect results 158 * on some degenerate inputs, so the client must have some way to 159 * handle this situation. 160 */ 161 return (u->s*(v->t - w->t) + v->s*(w->t - u->t) + w->s*(u->t - v->t)) >= 0; 162 } 163 164 /* Given parameters a,x,b,y returns the value (b*x+a*y)/(a+b), 165 * or (x+y)/2 if a==b==0. It requires that a,b >= 0, and enforces 166 * this in the rare case that one argument is slightly negative. 167 * The implementation is extremely stable numerically. 168 * In particular it guarantees that the result r satisfies 169 * MIN(x,y) <= r <= MAX(x,y), and the results are very accurate 170 * even when a and b differ greatly in magnitude. 171 */ 172 #define RealInterpolate(a,x,b,y) \ 173 (a = (a < 0) ? 0 : a, b = (b < 0) ? 0 : b, \ 174 ((a <= b) ? ((b == 0) ? ((x+y) / 2) \ 175 : (x + (y-x) * (a/(a+b)))) \ 176 : (y + (x-y) * (b/(a+b))))) 177 178 #ifndef FOR_TRITE_TEST_PROGRAM 179 #define Interpolate(a,x,b,y) RealInterpolate(a,x,b,y) 180 #else 181 182 /* Claim: the ONLY property the sweep algorithm relies on is that 183 * MIN(x,y) <= r <= MAX(x,y). This is a nasty way to test that. 184 */ 185 #include <stdlib.h> 186 extern int RandomInterpolate; 187 188 GLdouble Interpolate( GLdouble a, GLdouble x, GLdouble b, GLdouble y) 189 { 190 printf("*********************%d\n",RandomInterpolate); 191 if( RandomInterpolate ) { 192 a = 1.2 * drand48() - 0.1; 193 a = (a < 0) ? 0 : ((a > 1) ? 1 : a); 194 b = 1.0 - a; 195 } 196 return RealInterpolate(a,x,b,y); 197 } 198 199 #endif 200 201 #define Swap(a,b) do { GLUvertex *t = a; a = b; b = t; } while (0) 202 203 void __gl_edgeIntersect( GLUvertex *o1, GLUvertex *d1, 204 GLUvertex *o2, GLUvertex *d2, 205 GLUvertex *v ) 206 /* Given edges (o1,d1) and (o2,d2), compute their point of intersection. 207 * The computed point is guaranteed to lie in the intersection of the 208 * bounding rectangles defined by each edge. 209 */ 210 { 211 GLdouble z1, z2; 212 213 /* This is certainly not the most efficient way to find the intersection 214 * of two line segments, but it is very numerically stable. 215 * 216 * Strategy: find the two middle vertices in the VertLeq ordering, 217 * and interpolate the intersection s-value from these. Then repeat 218 * using the TransLeq ordering to find the intersection t-value. 219 */ 220 221 if( ! VertLeq( o1, d1 )) { Swap( o1, d1 ); } 222 if( ! VertLeq( o2, d2 )) { Swap( o2, d2 ); } 223 if( ! VertLeq( o1, o2 )) { Swap( o1, o2 ); Swap( d1, d2 ); } 224 225 if( ! VertLeq( o2, d1 )) { 226 /* Technically, no intersection -- do our best */ 227 v->s = (o2->s + d1->s) / 2; 228 } else if( VertLeq( d1, d2 )) { 229 /* Interpolate between o2 and d1 */ 230 z1 = EdgeEval( o1, o2, d1 ); 231 z2 = EdgeEval( o2, d1, d2 ); 232 if( z1+z2 < 0 ) { z1 = -z1; z2 = -z2; } 233 v->s = Interpolate( z1, o2->s, z2, d1->s ); 234 } else { 235 /* Interpolate between o2 and d2 */ 236 z1 = EdgeSign( o1, o2, d1 ); 237 z2 = -EdgeSign( o1, d2, d1 ); 238 if( z1+z2 < 0 ) { z1 = -z1; z2 = -z2; } 239 v->s = Interpolate( z1, o2->s, z2, d2->s ); 240 } 241 242 /* Now repeat the process for t */ 243 244 if( ! TransLeq( o1, d1 )) { Swap( o1, d1 ); } 245 if( ! TransLeq( o2, d2 )) { Swap( o2, d2 ); } 246 if( ! TransLeq( o1, o2 )) { Swap( o1, o2 ); Swap( d1, d2 ); } 247 248 if( ! TransLeq( o2, d1 )) { 249 /* Technically, no intersection -- do our best */ 250 v->t = (o2->t + d1->t) / 2; 251 } else if( TransLeq( d1, d2 )) { 252 /* Interpolate between o2 and d1 */ 253 z1 = TransEval( o1, o2, d1 ); 254 z2 = TransEval( o2, d1, d2 ); 255 if( z1+z2 < 0 ) { z1 = -z1; z2 = -z2; } 256 v->t = Interpolate( z1, o2->t, z2, d1->t ); 257 } else { 258 /* Interpolate between o2 and d2 */ 259 z1 = TransSign( o1, o2, d1 ); 260 z2 = -TransSign( o1, d2, d1 ); 261 if( z1+z2 < 0 ) { z1 = -z1; z2 = -z2; } 262 v->t = Interpolate( z1, o2->t, z2, d2->t ); 263 } 264 } 265