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