1 /*
2 Copyright (C) 1999-2006 Id Software, Inc. and contributors.
3 For a list of contributors, see the accompanying CONTRIBUTORS file.
4
5 This file is part of GtkRadiant.
6
7 GtkRadiant is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 GtkRadiant is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with GtkRadiant; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20 */
21
22
23 #include "cmdlib.h"
24 #include "mathlib.h"
25 #include "inout.h"
26 #include "polylib.h"
27 #include "qfiles.h"
28
29
30 extern int numthreads;
31
32 // counters are only bumped when running single threaded,
33 // because they are an awefull coherence problem
34 int c_active_windings;
35 int c_peak_windings;
36 int c_winding_allocs;
37 int c_winding_points;
38
39 #define BOGUS_RANGE WORLD_SIZE
40
pw(winding_t * w)41 void pw( winding_t *w ){
42 int i;
43 for ( i = 0 ; i < w->numpoints ; i++ )
44 Sys_Printf( "(%5.1f, %5.1f, %5.1f)\n",w->p[i][0], w->p[i][1],w->p[i][2] );
45 }
46
47
48 /*
49 =============
50 AllocWinding
51 =============
52 */
AllocWinding(int points)53 winding_t *AllocWinding( int points ){
54 winding_t *w;
55 int s;
56
57 if ( points >= MAX_POINTS_ON_WINDING ) {
58 Error( "AllocWinding failed: MAX_POINTS_ON_WINDING exceeded" );
59 }
60
61 if ( numthreads == 1 ) {
62 c_winding_allocs++;
63 c_winding_points += points;
64 c_active_windings++;
65 if ( c_active_windings > c_peak_windings ) {
66 c_peak_windings = c_active_windings;
67 }
68 }
69 s = sizeof( *w ) + ( points ? sizeof( w->p[0] ) * ( points - 1 ) : 0 );
70 w = safe_malloc( s );
71 memset( w, 0, s );
72 return w;
73 }
74
75 /*
76 =============
77 AllocWindingAccu
78 =============
79 */
AllocWindingAccu(int points)80 winding_accu_t *AllocWindingAccu( int points ){
81 winding_accu_t *w;
82 int s;
83
84 if ( points >= MAX_POINTS_ON_WINDING ) {
85 Error( "AllocWindingAccu failed: MAX_POINTS_ON_WINDING exceeded" );
86 }
87
88 if ( numthreads == 1 ) {
89 // At the time of this writing, these statistics were not used in any way.
90 c_winding_allocs++;
91 c_winding_points += points;
92 c_active_windings++;
93 if ( c_active_windings > c_peak_windings ) {
94 c_peak_windings = c_active_windings;
95 }
96 }
97 s = sizeof( *w ) + ( points ? sizeof( w->p[0] ) * ( points - 1 ) : 0 );
98 w = safe_malloc( s );
99 memset( w, 0, s );
100 return w;
101 }
102
103 /*
104 =============
105 FreeWinding
106 =============
107 */
FreeWinding(winding_t * w)108 void FreeWinding( winding_t *w ){
109 if ( !w ) {
110 Error( "FreeWinding: winding is NULL" );
111 }
112
113 if ( *(unsigned *)w == 0xdeaddead ) {
114 Error( "FreeWinding: freed a freed winding" );
115 }
116 *(unsigned *)w = 0xdeaddead;
117
118 if ( numthreads == 1 ) {
119 c_active_windings--;
120 }
121 free( w );
122 }
123
124 /*
125 =============
126 FreeWindingAccu
127 =============
128 */
FreeWindingAccu(winding_accu_t * w)129 void FreeWindingAccu( winding_accu_t *w ){
130 if ( !w ) {
131 Error( "FreeWindingAccu: winding is NULL" );
132 }
133
134 if ( *( (unsigned *) w ) == 0xdeaddead ) {
135 Error( "FreeWindingAccu: freed a freed winding" );
136 }
137 *( (unsigned *) w ) = 0xdeaddead;
138
139 if ( numthreads == 1 ) {
140 c_active_windings--;
141 }
142 free( w );
143 }
144
145 /*
146 ============
147 RemoveColinearPoints
148 ============
149 */
150 int c_removed;
151
RemoveColinearPoints(winding_t * w)152 void RemoveColinearPoints( winding_t *w ){
153 int i, j, k;
154 vec3_t v1, v2;
155 int nump;
156 vec3_t p[MAX_POINTS_ON_WINDING];
157
158 nump = 0;
159 for ( i = 0 ; i < w->numpoints ; i++ )
160 {
161 j = ( i + 1 ) % w->numpoints;
162 k = ( i + w->numpoints - 1 ) % w->numpoints;
163 VectorSubtract( w->p[j], w->p[i], v1 );
164 VectorSubtract( w->p[i], w->p[k], v2 );
165 VectorNormalize( v1,v1 );
166 VectorNormalize( v2,v2 );
167 if ( DotProduct( v1, v2 ) < 0.999 ) {
168 VectorCopy( w->p[i], p[nump] );
169 nump++;
170 }
171 }
172
173 if ( nump == w->numpoints ) {
174 return;
175 }
176
177 if ( numthreads == 1 ) {
178 c_removed += w->numpoints - nump;
179 }
180 w->numpoints = nump;
181 memcpy( w->p, p, nump * sizeof( p[0] ) );
182 }
183
184 /*
185 ============
186 WindingPlane
187 ============
188 */
WindingPlane(winding_t * w,vec3_t normal,vec_t * dist)189 void WindingPlane( winding_t *w, vec3_t normal, vec_t *dist ){
190 vec3_t v1, v2;
191
192 VectorSubtract( w->p[1], w->p[0], v1 );
193 VectorSubtract( w->p[2], w->p[0], v2 );
194 CrossProduct( v2, v1, normal );
195 VectorNormalize( normal, normal );
196 *dist = DotProduct( w->p[0], normal );
197
198 }
199
200 /*
201 =============
202 WindingArea
203 =============
204 */
WindingArea(winding_t * w)205 vec_t WindingArea( winding_t *w ){
206 int i;
207 vec3_t d1, d2, cross;
208 vec_t total;
209
210 total = 0;
211 for ( i = 2 ; i < w->numpoints ; i++ )
212 {
213 VectorSubtract( w->p[i - 1], w->p[0], d1 );
214 VectorSubtract( w->p[i], w->p[0], d2 );
215 CrossProduct( d1, d2, cross );
216 total += 0.5 * VectorLength( cross );
217 }
218 return total;
219 }
220
WindingBounds(winding_t * w,vec3_t mins,vec3_t maxs)221 void WindingBounds( winding_t *w, vec3_t mins, vec3_t maxs ){
222 vec_t v;
223 int i,j;
224
225 mins[0] = mins[1] = mins[2] = 99999;
226 maxs[0] = maxs[1] = maxs[2] = -99999;
227
228 for ( i = 0 ; i < w->numpoints ; i++ )
229 {
230 for ( j = 0 ; j < 3 ; j++ )
231 {
232 v = w->p[i][j];
233 if ( v < mins[j] ) {
234 mins[j] = v;
235 }
236 if ( v > maxs[j] ) {
237 maxs[j] = v;
238 }
239 }
240 }
241 }
242
243 /*
244 =============
245 WindingCenter
246 =============
247 */
WindingCenter(winding_t * w,vec3_t center)248 void WindingCenter( winding_t *w, vec3_t center ){
249 int i;
250 float scale;
251
252 VectorCopy( vec3_origin, center );
253 for ( i = 0 ; i < w->numpoints ; i++ )
254 VectorAdd( w->p[i], center, center );
255
256 scale = 1.0 / w->numpoints;
257 VectorScale( center, scale, center );
258 }
259
260 /*
261 =================
262 BaseWindingForPlaneAccu
263 =================
264 */
BaseWindingForPlaneAccu(vec3_t normal,vec_t dist)265 winding_accu_t *BaseWindingForPlaneAccu( vec3_t normal, vec_t dist ){
266 // The goal in this function is to replicate the behavior of the original BaseWindingForPlane()
267 // function (see below) but at the same time increasing accuracy substantially.
268
269 // The original code gave a preference for the vup vector to start out as (0, 0, 1), unless the
270 // normal had a dominant Z value, in which case vup started out as (1, 0, 0). After that, vup
271 // was "bent" [along the plane defined by normal and vup] to become perpendicular to normal.
272 // After that the vright vector was computed as the cross product of vup and normal.
273
274 // I'm constructing the winding polygon points in a fashion similar to the method used in the
275 // original function. Orientation is the same. The size of the winding polygon, however, is
276 // variable in this function (depending on the angle of normal), and is larger (by about a factor
277 // of 2) than the winding polygon in the original function.
278
279 int x, i;
280 vec_t max, v;
281 vec3_accu_t vright, vup, org, normalAccu;
282 winding_accu_t *w;
283
284 // One of the components of normal must have a magnitiude greater than this value,
285 // otherwise normal is not a unit vector. This is a little bit of inexpensive
286 // partial error checking we can do.
287 max = 0.56; // 1 / sqrt(1^2 + 1^2 + 1^2) = 0.577350269
288
289 x = -1;
290 for ( i = 0; i < 3; i++ ) {
291 v = (vec_t) fabs( normal[i] );
292 if ( v > max ) {
293 x = i;
294 max = v;
295 }
296 }
297 if ( x == -1 ) {
298 Error( "BaseWindingForPlaneAccu: no dominant axis found because normal is too short" );
299 }
300
301 switch ( x ) {
302 case 0: // Fall through to next case.
303 case 1:
304 vright[0] = (vec_accu_t) -normal[1];
305 vright[1] = (vec_accu_t) normal[0];
306 vright[2] = 0;
307 break;
308 case 2:
309 vright[0] = 0;
310 vright[1] = (vec_accu_t) -normal[2];
311 vright[2] = (vec_accu_t) normal[1];
312 break;
313 }
314
315 // vright and normal are now perpendicular; you can prove this by taking their
316 // dot product and seeing that it's always exactly 0 (with no error).
317
318 // NOTE: vright is NOT a unit vector at this point. vright will have length
319 // not exceeding 1.0. The minimum length that vright can achieve happens when,
320 // for example, the Z and X components of the normal input vector are equal,
321 // and when normal's Y component is zero. In that case Z and X of the normal
322 // vector are both approximately 0.70711. The resulting vright vector in this
323 // case will have a length of 0.70711.
324
325 // We're relying on the fact that MAX_WORLD_COORD is a power of 2 to keep
326 // our calculation precise and relatively free of floating point error.
327 // [However, the code will still work fine if that's not the case.]
328 VectorScaleAccu( vright, ( (vec_accu_t) MAX_WORLD_COORD ) * 4.0, vright );
329
330 // At time time of this writing, MAX_WORLD_COORD was 65536 (2^16). Therefore
331 // the length of vright at this point is at least 185364. In comparison, a
332 // corner of the world at location (65536, 65536, 65536) is distance 113512
333 // away from the origin.
334
335 VectorCopyRegularToAccu( normal, normalAccu );
336 CrossProductAccu( normalAccu, vright, vup );
337
338 // vup now has length equal to that of vright.
339
340 VectorScaleAccu( normalAccu, (vec_accu_t) dist, org );
341
342 // org is now a point on the plane defined by normal and dist. Furthermore,
343 // org, vright, and vup are pairwise perpendicular. Now, the 4 vectors
344 // { (+-)vright + (+-)vup } have length that is at least sqrt(185364^2 + 185364^2),
345 // which is about 262144. That length lies outside the world, since the furthest
346 // point in the world has distance 113512 from the origin as mentioned above.
347 // Also, these 4 vectors are perpendicular to the org vector. So adding them
348 // to org will only increase their length. Therefore the 4 points defined below
349 // all lie outside of the world. Furthermore, it can be easily seen that the
350 // edges connecting these 4 points (in the winding_accu_t below) lie completely
351 // outside the world. sqrt(262144^2 + 262144^2)/2 = 185363, which is greater than
352 // 113512.
353
354 w = AllocWindingAccu( 4 );
355
356 VectorSubtractAccu( org, vright, w->p[0] );
357 VectorAddAccu( w->p[0], vup, w->p[0] );
358
359 VectorAddAccu( org, vright, w->p[1] );
360 VectorAddAccu( w->p[1], vup, w->p[1] );
361
362 VectorAddAccu( org, vright, w->p[2] );
363 VectorSubtractAccu( w->p[2], vup, w->p[2] );
364
365 VectorSubtractAccu( org, vright, w->p[3] );
366 VectorSubtractAccu( w->p[3], vup, w->p[3] );
367
368 w->numpoints = 4;
369
370 return w;
371 }
372
373 /*
374 =================
375 BaseWindingForPlane
376
377 Original BaseWindingForPlane() function that has serious accuracy problems. Here is why.
378 The base winding is computed as a rectangle with very large coordinates. These coordinates
379 are in the range 2^17 or 2^18. "Epsilon" (meaning the distance between two adjacent numbers)
380 at these magnitudes in 32 bit floating point world is about 0.02. So for example, if things
381 go badly (by bad luck), then the whole plane could be shifted by 0.02 units (its distance could
382 be off by that much). Then if we were to compute the winding of this plane and another of
383 the brush's planes met this winding at a very acute angle, that error could multiply to around
384 0.1 or more when computing the final vertex coordinates of the winding. 0.1 is a very large
385 error, and can lead to all sorts of disappearing triangle problems.
386 =================
387 */
BaseWindingForPlane(vec3_t normal,vec_t dist)388 winding_t *BaseWindingForPlane( vec3_t normal, vec_t dist ){
389 int i, x;
390 vec_t max, v;
391 vec3_t org, vright, vup;
392 winding_t *w;
393
394 // find the major axis
395
396 max = -BOGUS_RANGE;
397 x = -1;
398 for ( i = 0 ; i < 3; i++ )
399 {
400 v = fabs( normal[i] );
401 if ( v > max ) {
402 x = i;
403 max = v;
404 }
405 }
406 if ( x == -1 ) {
407 Error( "BaseWindingForPlane: no axis found" );
408 }
409
410 VectorCopy( vec3_origin, vup );
411 switch ( x )
412 {
413 case 0:
414 case 1:
415 vup[2] = 1;
416 break;
417 case 2:
418 vup[0] = 1;
419 break;
420 }
421
422 v = DotProduct( vup, normal );
423 VectorMA( vup, -v, normal, vup );
424 VectorNormalize( vup, vup );
425
426 VectorScale( normal, dist, org );
427
428 CrossProduct( vup, normal, vright );
429
430 // LordHavoc: this has to use *2 because otherwise some created points may
431 // be inside the world (think of a diagonal case), and any brush with such
432 // points should be removed, failure to detect such cases is disasterous
433 VectorScale( vup, MAX_WORLD_COORD * 2, vup );
434 VectorScale( vright, MAX_WORLD_COORD * 2, vright );
435
436 // project a really big axis aligned box onto the plane
437 w = AllocWinding( 4 );
438
439 VectorSubtract( org, vright, w->p[0] );
440 VectorAdd( w->p[0], vup, w->p[0] );
441
442 VectorAdd( org, vright, w->p[1] );
443 VectorAdd( w->p[1], vup, w->p[1] );
444
445 VectorAdd( org, vright, w->p[2] );
446 VectorSubtract( w->p[2], vup, w->p[2] );
447
448 VectorSubtract( org, vright, w->p[3] );
449 VectorSubtract( w->p[3], vup, w->p[3] );
450
451 w->numpoints = 4;
452
453 return w;
454 }
455
456 /*
457 ==================
458 CopyWinding
459 ==================
460 */
CopyWinding(winding_t * w)461 winding_t *CopyWinding( winding_t *w ){
462 size_t size;
463 winding_t *c;
464
465 if ( !w ) {
466 Error( "CopyWinding: winding is NULL" );
467 }
468
469 c = AllocWinding( w->numpoints );
470 size = (size_t)( (winding_t *)NULL )->p[w->numpoints];
471 memcpy( c, w, size );
472 return c;
473 }
474
475 /*
476 ==================
477 CopyWindingAccuIncreaseSizeAndFreeOld
478 ==================
479 */
CopyWindingAccuIncreaseSizeAndFreeOld(winding_accu_t * w)480 winding_accu_t *CopyWindingAccuIncreaseSizeAndFreeOld( winding_accu_t *w ){
481 int i;
482 winding_accu_t *c;
483
484 if ( !w ) {
485 Error( "CopyWindingAccuIncreaseSizeAndFreeOld: winding is NULL" );
486 }
487
488 c = AllocWindingAccu( w->numpoints + 1 );
489 c->numpoints = w->numpoints;
490 for ( i = 0; i < c->numpoints; i++ )
491 {
492 VectorCopyAccu( w->p[i], c->p[i] );
493 }
494 FreeWindingAccu( w );
495 return c;
496 }
497
498 /*
499 ==================
500 CopyWindingAccuToRegular
501 ==================
502 */
CopyWindingAccuToRegular(winding_accu_t * w)503 winding_t *CopyWindingAccuToRegular( winding_accu_t *w ){
504 int i;
505 winding_t *c;
506
507 if ( !w ) {
508 Error( "CopyWindingAccuToRegular: winding is NULL" );
509 }
510
511 c = AllocWinding( w->numpoints );
512 c->numpoints = w->numpoints;
513 for ( i = 0; i < c->numpoints; i++ )
514 {
515 VectorCopyAccuToRegular( w->p[i], c->p[i] );
516 }
517 return c;
518 }
519
520 /*
521 ==================
522 ReverseWinding
523 ==================
524 */
ReverseWinding(winding_t * w)525 winding_t *ReverseWinding( winding_t *w ){
526 int i;
527 winding_t *c;
528
529 c = AllocWinding( w->numpoints );
530 for ( i = 0 ; i < w->numpoints ; i++ )
531 {
532 VectorCopy( w->p[w->numpoints - 1 - i], c->p[i] );
533 }
534 c->numpoints = w->numpoints;
535 return c;
536 }
537
538
539 /*
540 =============
541 ClipWindingEpsilon
542 =============
543 */
ClipWindingEpsilonStrict(winding_t * in,vec3_t normal,vec_t dist,vec_t epsilon,winding_t ** front,winding_t ** back)544 void ClipWindingEpsilonStrict( winding_t *in, vec3_t normal, vec_t dist,
545 vec_t epsilon, winding_t **front, winding_t **back ){
546 vec_t dists[MAX_POINTS_ON_WINDING + 4];
547 int sides[MAX_POINTS_ON_WINDING + 4];
548 int counts[3];
549 static vec_t dot; // VC 4.2 optimizer bug if not static
550 int i, j;
551 vec_t *p1, *p2;
552 vec3_t mid;
553 winding_t *f, *b;
554 int maxpts;
555
556 counts[0] = counts[1] = counts[2] = 0;
557
558 // determine sides for each point
559 for ( i = 0 ; i < in->numpoints ; i++ )
560 {
561
562 dot = DotProduct( in->p[i], normal );
563 dot -= dist;
564 dists[i] = dot;
565 if ( dot > epsilon ) {
566 sides[i] = SIDE_FRONT;
567 }
568 else if ( dot < -epsilon ) {
569 sides[i] = SIDE_BACK;
570 }
571 else
572 {
573 sides[i] = SIDE_ON;
574 }
575 counts[sides[i]]++;
576 }
577 sides[i] = sides[0];
578 dists[i] = dists[0];
579
580 *front = *back = NULL;
581
582 if ( !counts[0] && !counts[1] ) {
583 return;
584 }
585 if ( !counts[0] ) {
586 *back = CopyWinding( in );
587 return;
588 }
589 if ( !counts[1] ) {
590 *front = CopyWinding( in );
591 return;
592 }
593
594 maxpts = in->numpoints + 4; // cant use counts[0]+2 because
595 // of fp grouping errors
596
597 *front = f = AllocWinding( maxpts );
598 *back = b = AllocWinding( maxpts );
599
600 for ( i = 0 ; i < in->numpoints ; i++ )
601 {
602 p1 = in->p[i];
603
604 if ( sides[i] == SIDE_ON ) {
605 VectorCopy( p1, f->p[f->numpoints] );
606 f->numpoints++;
607 VectorCopy( p1, b->p[b->numpoints] );
608 b->numpoints++;
609 continue;
610 }
611
612 if ( sides[i] == SIDE_FRONT ) {
613 VectorCopy( p1, f->p[f->numpoints] );
614 f->numpoints++;
615 }
616 if ( sides[i] == SIDE_BACK ) {
617 VectorCopy( p1, b->p[b->numpoints] );
618 b->numpoints++;
619 }
620
621 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
622 continue;
623 }
624
625 // generate a split point
626 p2 = in->p[( i + 1 ) % in->numpoints];
627
628 dot = dists[i] / ( dists[i] - dists[i + 1] );
629 for ( j = 0 ; j < 3 ; j++ )
630 { // avoid round off error when possible
631 if ( normal[j] == 1 ) {
632 mid[j] = dist;
633 }
634 else if ( normal[j] == -1 ) {
635 mid[j] = -dist;
636 }
637 else{
638 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
639 }
640 }
641
642 VectorCopy( mid, f->p[f->numpoints] );
643 f->numpoints++;
644 VectorCopy( mid, b->p[b->numpoints] );
645 b->numpoints++;
646 }
647
648 if ( f->numpoints > maxpts || b->numpoints > maxpts ) {
649 Error( "ClipWinding: points exceeded estimate" );
650 }
651 if ( f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING ) {
652 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
653 }
654 }
655
ClipWindingEpsilon(winding_t * in,vec3_t normal,vec_t dist,vec_t epsilon,winding_t ** front,winding_t ** back)656 void ClipWindingEpsilon( winding_t *in, vec3_t normal, vec_t dist,
657 vec_t epsilon, winding_t **front, winding_t **back ){
658 ClipWindingEpsilonStrict( in, normal, dist, epsilon, front, back );
659 /* apparently most code expects that in the winding-on-plane case, the back winding is the original winding */
660 if ( !*front && !*back ) {
661 *back = CopyWinding( in );
662 }
663 }
664
665
666 /*
667 =============
668 ChopWindingInPlaceAccu
669 =============
670 */
ChopWindingInPlaceAccu(winding_accu_t ** inout,vec3_t normal,vec_t dist,vec_t crudeEpsilon)671 void ChopWindingInPlaceAccu( winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon ){
672 vec_accu_t fineEpsilon;
673 winding_accu_t *in;
674 int counts[3];
675 int i, j;
676 vec_accu_t dists[MAX_POINTS_ON_WINDING + 1];
677 int sides[MAX_POINTS_ON_WINDING + 1];
678 int maxpts;
679 winding_accu_t *f;
680 vec_accu_t *p1, *p2;
681 vec_accu_t w;
682 vec3_accu_t mid, normalAccu;
683
684 // We require at least a very small epsilon. It's a good idea for several reasons.
685 // First, we will be dividing by a potentially very small distance below. We don't
686 // want that distance to be too small; otherwise, things "blow up" with little accuracy
687 // due to the division. (After a second look, the value w below is in range (0,1), but
688 // graininess problem remains.) Second, Having minimum epsilon also prevents the following
689 // situation. Say for example we have a perfect octagon defined by the input winding.
690 // Say our chopping plane (defined by normal and dist) is essentially the same plane
691 // that the octagon is sitting on. Well, due to rounding errors, it may be that point
692 // 1 of the octagon might be in front, point 2 might be in back, point 3 might be in
693 // front, point 4 might be in back, and so on. So we could end up with a very ugly-
694 // looking chopped winding, and this might be undesirable, and would at least lead to
695 // a possible exhaustion of MAX_POINTS_ON_WINDING. It's better to assume that points
696 // very very close to the plane are on the plane, using an infinitesimal epsilon amount.
697
698 // Now, the original ChopWindingInPlace() function used a vec_t-based winding_t.
699 // So this minimum epsilon is quite similar to casting the higher resolution numbers to
700 // the lower resolution and comparing them in the lower resolution mode. We explicitly
701 // choose the minimum epsilon as something around the vec_t epsilon of one because we
702 // want the resolution of vec_accu_t to have a large resolution around the epsilon.
703 // Some of that leftover resolution even goes away after we scale to points far away.
704
705 // Here is a further discussion regarding the choice of smallestEpsilonAllowed.
706 // In the 32 float world (we can assume vec_t is that), the "epsilon around 1.0" is
707 // 0.00000011921. In the 64 bit float world (we can assume vec_accu_t is that), the
708 // "epsilon around 1.0" is 0.00000000000000022204. (By the way these two epsilons
709 // are defined as VEC_SMALLEST_EPSILON_AROUND_ONE VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE
710 // respectively.) If you divide the first by the second, you get approximately
711 // 536,885,246. Dividing that number by 200,000 (a typical base winding coordinate)
712 // gives 2684. So in other words, if our smallestEpsilonAllowed was chosen as exactly
713 // VEC_SMALLEST_EPSILON_AROUND_ONE, you would be guaranteed at least 2000 "ticks" in
714 // 64-bit land inside of the epsilon for all numbers we're dealing with.
715
716 static const vec_accu_t smallestEpsilonAllowed = ( (vec_accu_t) VEC_SMALLEST_EPSILON_AROUND_ONE ) * 0.5;
717 if ( crudeEpsilon < smallestEpsilonAllowed ) {
718 fineEpsilon = smallestEpsilonAllowed;
719 }
720 else{fineEpsilon = (vec_accu_t) crudeEpsilon; }
721
722 in = *inout;
723 counts[0] = counts[1] = counts[2] = 0;
724 VectorCopyRegularToAccu( normal, normalAccu );
725
726 for ( i = 0; i < in->numpoints; i++ )
727 {
728 dists[i] = DotProductAccu( in->p[i], normalAccu ) - dist;
729 if ( dists[i] > fineEpsilon ) {
730 sides[i] = SIDE_FRONT;
731 }
732 else if ( dists[i] < -fineEpsilon ) {
733 sides[i] = SIDE_BACK;
734 }
735 else{sides[i] = SIDE_ON; }
736 counts[sides[i]]++;
737 }
738 sides[i] = sides[0];
739 dists[i] = dists[0];
740
741 // I'm wondering if whatever code that handles duplicate planes is robust enough
742 // that we never get a case where two nearly equal planes result in 2 NULL windings
743 // due to the 'if' statement below. TODO: Investigate this.
744 if ( !counts[SIDE_FRONT] ) {
745 FreeWindingAccu( in );
746 *inout = NULL;
747 return;
748 }
749 if ( !counts[SIDE_BACK] ) {
750 return; // Winding is unmodified.
751 }
752
753 // NOTE: The least number of points that a winding can have at this point is 2.
754 // In that case, one point is SIDE_FRONT and the other is SIDE_BACK.
755
756 maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
757 f = AllocWindingAccu( maxpts );
758
759 for ( i = 0; i < in->numpoints; i++ )
760 {
761 p1 = in->p[i];
762
763 if ( sides[i] == SIDE_ON || sides[i] == SIDE_FRONT ) {
764 if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
765 Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
766 }
767 if ( f->numpoints >= maxpts ) { // This will probably never happen.
768 Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
769 f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
770 maxpts++;
771 }
772 VectorCopyAccu( p1, f->p[f->numpoints] );
773 f->numpoints++;
774 if ( sides[i] == SIDE_ON ) {
775 continue;
776 }
777 }
778 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
779 continue;
780 }
781
782 // Generate a split point.
783 p2 = in->p[( ( i + 1 ) == in->numpoints ) ? 0 : ( i + 1 )];
784
785 // The divisor's absolute value is greater than the dividend's absolute value.
786 // w is in the range (0,1).
787 w = dists[i] / ( dists[i] - dists[i + 1] );
788
789 for ( j = 0; j < 3; j++ )
790 {
791 // Avoid round-off error when possible. Check axis-aligned normal.
792 if ( normal[j] == 1 ) {
793 mid[j] = dist;
794 }
795 else if ( normal[j] == -1 ) {
796 mid[j] = -dist;
797 }
798 else{mid[j] = p1[j] + ( w * ( p2[j] - p1[j] ) ); }
799 }
800 if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
801 Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
802 }
803 if ( f->numpoints >= maxpts ) { // This will probably never happen.
804 Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
805 f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
806 maxpts++;
807 }
808 VectorCopyAccu( mid, f->p[f->numpoints] );
809 f->numpoints++;
810 }
811
812 FreeWindingAccu( in );
813 *inout = f;
814 }
815
816 /*
817 =============
818 ChopWindingInPlace
819 =============
820 */
ChopWindingInPlace(winding_t ** inout,vec3_t normal,vec_t dist,vec_t epsilon)821 void ChopWindingInPlace( winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon ){
822 winding_t *in;
823 vec_t dists[MAX_POINTS_ON_WINDING + 4];
824 int sides[MAX_POINTS_ON_WINDING + 4];
825 int counts[3];
826 static vec_t dot; // VC 4.2 optimizer bug if not static
827 int i, j;
828 vec_t *p1, *p2;
829 vec3_t mid;
830 winding_t *f;
831 int maxpts;
832
833 in = *inout;
834 counts[0] = counts[1] = counts[2] = 0;
835
836 // determine sides for each point
837 for ( i = 0 ; i < in->numpoints ; i++ )
838 {
839 dot = DotProduct( in->p[i], normal );
840 dot -= dist;
841 dists[i] = dot;
842 if ( dot > epsilon ) {
843 sides[i] = SIDE_FRONT;
844 }
845 else if ( dot < -epsilon ) {
846 sides[i] = SIDE_BACK;
847 }
848 else
849 {
850 sides[i] = SIDE_ON;
851 }
852 counts[sides[i]]++;
853 }
854 sides[i] = sides[0];
855 dists[i] = dists[0];
856
857 if ( !counts[0] ) {
858 FreeWinding( in );
859 *inout = NULL;
860 return;
861 }
862 if ( !counts[1] ) {
863 return; // inout stays the same
864
865 }
866 maxpts = in->numpoints + 4; // cant use counts[0]+2 because
867 // of fp grouping errors
868
869 f = AllocWinding( maxpts );
870
871 for ( i = 0 ; i < in->numpoints ; i++ )
872 {
873 p1 = in->p[i];
874
875 if ( sides[i] == SIDE_ON ) {
876 VectorCopy( p1, f->p[f->numpoints] );
877 f->numpoints++;
878 continue;
879 }
880
881 if ( sides[i] == SIDE_FRONT ) {
882 VectorCopy( p1, f->p[f->numpoints] );
883 f->numpoints++;
884 }
885
886 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
887 continue;
888 }
889
890 // generate a split point
891 p2 = in->p[( i + 1 ) % in->numpoints];
892
893 dot = dists[i] / ( dists[i] - dists[i + 1] );
894 for ( j = 0 ; j < 3 ; j++ )
895 { // avoid round off error when possible
896 if ( normal[j] == 1 ) {
897 mid[j] = dist;
898 }
899 else if ( normal[j] == -1 ) {
900 mid[j] = -dist;
901 }
902 else{
903 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
904 }
905 }
906
907 VectorCopy( mid, f->p[f->numpoints] );
908 f->numpoints++;
909 }
910
911 if ( f->numpoints > maxpts ) {
912 Error( "ClipWinding: points exceeded estimate" );
913 }
914 if ( f->numpoints > MAX_POINTS_ON_WINDING ) {
915 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
916 }
917
918 FreeWinding( in );
919 *inout = f;
920 }
921
922
923 /*
924 =================
925 ChopWinding
926
927 Returns the fragment of in that is on the front side
928 of the cliping plane. The original is freed.
929 =================
930 */
ChopWinding(winding_t * in,vec3_t normal,vec_t dist)931 winding_t *ChopWinding( winding_t *in, vec3_t normal, vec_t dist ){
932 winding_t *f, *b;
933
934 ClipWindingEpsilon( in, normal, dist, ON_EPSILON, &f, &b );
935 FreeWinding( in );
936 if ( b ) {
937 FreeWinding( b );
938 }
939 return f;
940 }
941
942
943 /*
944 =================
945 CheckWinding
946
947 =================
948 */
CheckWinding(winding_t * w)949 void CheckWinding( winding_t *w ){
950 int i, j;
951 vec_t *p1, *p2;
952 vec_t d, edgedist;
953 vec3_t dir, edgenormal, facenormal;
954 vec_t area;
955 vec_t facedist;
956
957 if ( w->numpoints < 3 ) {
958 Error( "CheckWinding: %i points",w->numpoints );
959 }
960
961 area = WindingArea( w );
962 if ( area < 1 ) {
963 Error( "CheckWinding: %f area", area );
964 }
965
966 WindingPlane( w, facenormal, &facedist );
967
968 for ( i = 0 ; i < w->numpoints ; i++ )
969 {
970 p1 = w->p[i];
971
972 for ( j = 0 ; j < 3 ; j++ )
973 if ( p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD ) {
974 Error( "CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j] );
975 }
976
977 j = i + 1 == w->numpoints ? 0 : i + 1;
978
979 // check the point is on the face plane
980 d = DotProduct( p1, facenormal ) - facedist;
981 if ( d < -ON_EPSILON || d > ON_EPSILON ) {
982 Error( "CheckWinding: point off plane" );
983 }
984
985 // check the edge isnt degenerate
986 p2 = w->p[j];
987 VectorSubtract( p2, p1, dir );
988
989 if ( VectorLength( dir ) < ON_EPSILON ) {
990 Error( "CheckWinding: degenerate edge" );
991 }
992
993 CrossProduct( facenormal, dir, edgenormal );
994 VectorNormalize( edgenormal, edgenormal );
995 edgedist = DotProduct( p1, edgenormal );
996 edgedist += ON_EPSILON;
997
998 // all other points must be on front side
999 for ( j = 0 ; j < w->numpoints ; j++ )
1000 {
1001 if ( j == i ) {
1002 continue;
1003 }
1004 d = DotProduct( w->p[j], edgenormal );
1005 if ( d > edgedist ) {
1006 Error( "CheckWinding: non-convex" );
1007 }
1008 }
1009 }
1010 }
1011
1012
1013 /*
1014 ============
1015 WindingOnPlaneSide
1016 ============
1017 */
WindingOnPlaneSide(winding_t * w,vec3_t normal,vec_t dist)1018 int WindingOnPlaneSide( winding_t *w, vec3_t normal, vec_t dist ){
1019 qboolean front, back;
1020 int i;
1021 vec_t d;
1022
1023 front = qfalse;
1024 back = qfalse;
1025 for ( i = 0 ; i < w->numpoints ; i++ )
1026 {
1027 d = DotProduct( w->p[i], normal ) - dist;
1028 if ( d < -ON_EPSILON ) {
1029 if ( front ) {
1030 return SIDE_CROSS;
1031 }
1032 back = qtrue;
1033 continue;
1034 }
1035 if ( d > ON_EPSILON ) {
1036 if ( back ) {
1037 return SIDE_CROSS;
1038 }
1039 front = qtrue;
1040 continue;
1041 }
1042 }
1043
1044 if ( back ) {
1045 return SIDE_BACK;
1046 }
1047 if ( front ) {
1048 return SIDE_FRONT;
1049 }
1050 return SIDE_ON;
1051 }
1052
1053
1054 /*
1055 =================
1056 AddWindingToConvexHull
1057
1058 Both w and *hull are on the same plane
1059 =================
1060 */
1061 #define MAX_HULL_POINTS 128
AddWindingToConvexHull(winding_t * w,winding_t ** hull,vec3_t normal)1062 void AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
1063 int i, j, k;
1064 float *p, *copy;
1065 vec3_t dir;
1066 float d;
1067 int numHullPoints, numNew;
1068 vec3_t hullPoints[MAX_HULL_POINTS];
1069 vec3_t newHullPoints[MAX_HULL_POINTS];
1070 vec3_t hullDirs[MAX_HULL_POINTS];
1071 qboolean hullSide[MAX_HULL_POINTS];
1072 qboolean outside;
1073
1074 if ( !*hull ) {
1075 *hull = CopyWinding( w );
1076 return;
1077 }
1078
1079 numHullPoints = ( *hull )->numpoints;
1080 memcpy( hullPoints, ( *hull )->p, numHullPoints * sizeof( vec3_t ) );
1081
1082 for ( i = 0 ; i < w->numpoints ; i++ ) {
1083 p = w->p[i];
1084
1085 // calculate hull side vectors
1086 for ( j = 0 ; j < numHullPoints ; j++ ) {
1087 k = ( j + 1 ) % numHullPoints;
1088
1089 VectorSubtract( hullPoints[k], hullPoints[j], dir );
1090 VectorNormalize( dir, dir );
1091 CrossProduct( normal, dir, hullDirs[j] );
1092 }
1093
1094 outside = qfalse;
1095 for ( j = 0 ; j < numHullPoints ; j++ ) {
1096 VectorSubtract( p, hullPoints[j], dir );
1097 d = DotProduct( dir, hullDirs[j] );
1098 if ( d >= ON_EPSILON ) {
1099 outside = qtrue;
1100 }
1101 if ( d >= -ON_EPSILON ) {
1102 hullSide[j] = qtrue;
1103 }
1104 else {
1105 hullSide[j] = qfalse;
1106 }
1107 }
1108
1109 // if the point is effectively inside, do nothing
1110 if ( !outside ) {
1111 continue;
1112 }
1113
1114 // find the back side to front side transition
1115 for ( j = 0 ; j < numHullPoints ; j++ ) {
1116 if ( !hullSide[ j % numHullPoints ] && hullSide[ ( j + 1 ) % numHullPoints ] ) {
1117 break;
1118 }
1119 }
1120 if ( j == numHullPoints ) {
1121 continue;
1122 }
1123
1124 // insert the point here
1125 VectorCopy( p, newHullPoints[0] );
1126 numNew = 1;
1127
1128 // copy over all points that aren't double fronts
1129 j = ( j + 1 ) % numHullPoints;
1130 for ( k = 0 ; k < numHullPoints ; k++ ) {
1131 if ( hullSide[ ( j + k ) % numHullPoints ] && hullSide[ ( j + k + 1 ) % numHullPoints ] ) {
1132 continue;
1133 }
1134 copy = hullPoints[ ( j + k + 1 ) % numHullPoints ];
1135 VectorCopy( copy, newHullPoints[numNew] );
1136 numNew++;
1137 }
1138
1139 numHullPoints = numNew;
1140 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof( vec3_t ) );
1141 }
1142
1143 FreeWinding( *hull );
1144 w = AllocWinding( numHullPoints );
1145 w->numpoints = numHullPoints;
1146 *hull = w;
1147 memcpy( w->p, hullPoints, numHullPoints * sizeof( vec3_t ) );
1148 }
1149