2 Copyright (C) 1999-2007 id Software, Inc. and contributors.
3 For a list of contributors, see the accompanying CONTRIBUTORS file.
5 This file is part of GtkRadiant.
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.
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.
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
29 extern int numthreads;
31 // counters are only bumped when running single threaded,
32 // because they are an awefull coherence problem
33 int c_active_windings;
38 #define BOGUS_RANGE WORLD_SIZE
40 void pw( winding_t *w ){
42 for ( i = 0 ; i < w->numpoints ; i++ )
43 Sys_Printf( "(%5.1f, %5.1f, %5.1f)\n",w->p[i][0], w->p[i][1],w->p[i][2] );
52 winding_t *AllocWinding( int points ){
56 if ( points >= MAX_POINTS_ON_WINDING ) {
57 Error( "AllocWinding failed: MAX_POINTS_ON_WINDING exceeded" );
60 if ( numthreads == 1 ) {
62 c_winding_points += points;
64 if ( c_active_windings > c_peak_windings ) {
65 c_peak_windings = c_active_windings;
68 s = sizeof( *w ) + ( points ? sizeof( w->p[0] ) * ( points - 1 ) : 0 );
69 w = safe_malloc0( s );
78 winding_accu_t *AllocWindingAccu( int points ){
82 if ( points >= MAX_POINTS_ON_WINDING ) {
83 Error( "AllocWindingAccu failed: MAX_POINTS_ON_WINDING exceeded" );
86 if ( numthreads == 1 ) {
87 // At the time of this writing, these statistics were not used in any way.
89 c_winding_points += points;
91 if ( c_active_windings > c_peak_windings ) {
92 c_peak_windings = c_active_windings;
95 s = sizeof( *w ) + ( points ? sizeof( w->p[0] ) * ( points - 1 ) : 0 );
96 w = safe_malloc0( s );
105 void FreeWinding( winding_t *w ){
107 Error( "FreeWinding: winding is NULL" );
110 if ( *(unsigned *)w == 0xdeaddead ) {
111 Error( "FreeWinding: freed a freed winding" );
113 *(unsigned *)w = 0xdeaddead;
115 if ( numthreads == 1 ) {
126 void FreeWindingAccu( winding_accu_t *w ){
128 Error( "FreeWindingAccu: winding is NULL" );
131 if ( *( (unsigned *) w ) == 0xdeaddead ) {
132 Error( "FreeWindingAccu: freed a freed winding" );
134 *( (unsigned *) w ) = 0xdeaddead;
136 if ( numthreads == 1 ) {
149 void RemoveColinearPoints( winding_t *w ){
153 vec3_t p[MAX_POINTS_ON_WINDING];
156 for ( i = 0 ; i < w->numpoints ; i++ )
158 j = ( i + 1 ) % w->numpoints;
159 k = ( i + w->numpoints - 1 ) % w->numpoints;
160 VectorSubtract( w->p[j], w->p[i], v1 );
161 VectorSubtract( w->p[i], w->p[k], v2 );
162 VectorNormalize( v1,v1 );
163 VectorNormalize( v2,v2 );
164 if ( DotProduct( v1, v2 ) < 0.999 ) {
165 VectorCopy( w->p[i], p[nump] );
170 if ( nump == w->numpoints ) {
174 if ( numthreads == 1 ) {
175 c_removed += w->numpoints - nump;
178 memcpy( w->p, p, nump * sizeof( p[0] ) );
186 void WindingPlane( winding_t *w, vec3_t normal, vec_t *dist ){
189 VectorSubtract( w->p[1], w->p[0], v1 );
190 VectorSubtract( w->p[2], w->p[0], v2 );
191 CrossProduct( v2, v1, normal );
192 VectorNormalize( normal, normal );
193 *dist = DotProduct( w->p[0], normal );
202 vec_t WindingArea( winding_t *w ){
204 vec3_t d1, d2, cross;
208 for ( i = 2 ; i < w->numpoints ; i++ )
210 VectorSubtract( w->p[i - 1], w->p[0], d1 );
211 VectorSubtract( w->p[i], w->p[0], d2 );
212 CrossProduct( d1, d2, cross );
213 total += 0.5 * VectorLength( cross );
218 void WindingBounds( winding_t *w, vec3_t mins, vec3_t maxs ){
222 mins[0] = mins[1] = mins[2] = 99999;
223 maxs[0] = maxs[1] = maxs[2] = -99999;
225 for ( i = 0 ; i < w->numpoints ; i++ )
227 for ( j = 0 ; j < 3 ; j++ )
245 void WindingCenter( winding_t *w, vec3_t center ){
249 VectorCopy( vec3_origin, center );
250 for ( i = 0 ; i < w->numpoints ; i++ )
251 VectorAdd( w->p[i], center, center );
253 scale = 1.0 / w->numpoints;
254 VectorScale( center, scale, center );
259 BaseWindingForPlaneAccu
262 winding_accu_t *BaseWindingForPlaneAccu( vec3_t normal, vec_t dist ){
263 // The goal in this function is to replicate the behavior of the original BaseWindingForPlane()
264 // function (see below) but at the same time increasing accuracy substantially.
266 // The original code gave a preference for the vup vector to start out as (0, 0, 1), unless the
267 // normal had a dominant Z value, in which case vup started out as (1, 0, 0). After that, vup
268 // was "bent" [along the plane defined by normal and vup] to become perpendicular to normal.
269 // After that the vright vector was computed as the cross product of vup and normal.
271 // I'm constructing the winding polygon points in a fashion similar to the method used in the
272 // original function. Orientation is the same. The size of the winding polygon, however, is
273 // variable in this function (depending on the angle of normal), and is larger (by about a factor
274 // of 2) than the winding polygon in the original function.
278 vec3_accu_t vright, vup, org, normalAccu;
281 // One of the components of normal must have a magnitiude greater than this value,
282 // otherwise normal is not a unit vector. This is a little bit of inexpensive
283 // partial error checking we can do.
284 max = 0.56; // 1 / sqrt(1^2 + 1^2 + 1^2) = 0.577350269
287 for ( i = 0; i < 3; i++ ) {
288 v = (vec_t) fabs( normal[i] );
295 Error( "BaseWindingForPlaneAccu: no dominant axis found because normal is too short" );
299 case 0: // Fall through to next case.
301 vright[0] = (vec_accu_t) -normal[1];
302 vright[1] = (vec_accu_t) normal[0];
307 vright[1] = (vec_accu_t) -normal[2];
308 vright[2] = (vec_accu_t) normal[1];
312 // vright and normal are now perpendicular; you can prove this by taking their
313 // dot product and seeing that it's always exactly 0 (with no error).
315 // NOTE: vright is NOT a unit vector at this point. vright will have length
316 // not exceeding 1.0. The minimum length that vright can achieve happens when,
317 // for example, the Z and X components of the normal input vector are equal,
318 // and when normal's Y component is zero. In that case Z and X of the normal
319 // vector are both approximately 0.70711. The resulting vright vector in this
320 // case will have a length of 0.70711.
322 // We're relying on the fact that MAX_WORLD_COORD is a power of 2 to keep
323 // our calculation precise and relatively free of floating point error.
324 // [However, the code will still work fine if that's not the case.]
325 VectorScaleAccu( vright, ( (vec_accu_t) MAX_WORLD_COORD ) * 4.0, vright );
327 // At time time of this writing, MAX_WORLD_COORD was 65536 (2^16). Therefore
328 // the length of vright at this point is at least 185364. In comparison, a
329 // corner of the world at location (65536, 65536, 65536) is distance 113512
330 // away from the origin.
332 VectorCopyRegularToAccu( normal, normalAccu );
333 CrossProductAccu( normalAccu, vright, vup );
335 // vup now has length equal to that of vright.
337 VectorScaleAccu( normalAccu, (vec_accu_t) dist, org );
339 // org is now a point on the plane defined by normal and dist. Furthermore,
340 // org, vright, and vup are pairwise perpendicular. Now, the 4 vectors
341 // { (+-)vright + (+-)vup } have length that is at least sqrt(185364^2 + 185364^2),
342 // which is about 262144. That length lies outside the world, since the furthest
343 // point in the world has distance 113512 from the origin as mentioned above.
344 // Also, these 4 vectors are perpendicular to the org vector. So adding them
345 // to org will only increase their length. Therefore the 4 points defined below
346 // all lie outside of the world. Furthermore, it can be easily seen that the
347 // edges connecting these 4 points (in the winding_accu_t below) lie completely
348 // outside the world. sqrt(262144^2 + 262144^2)/2 = 185363, which is greater than
351 w = AllocWindingAccu( 4 );
353 VectorSubtractAccu( org, vright, w->p[0] );
354 VectorAddAccu( w->p[0], vup, w->p[0] );
356 VectorAddAccu( org, vright, w->p[1] );
357 VectorAddAccu( w->p[1], vup, w->p[1] );
359 VectorAddAccu( org, vright, w->p[2] );
360 VectorSubtractAccu( w->p[2], vup, w->p[2] );
362 VectorSubtractAccu( org, vright, w->p[3] );
363 VectorSubtractAccu( w->p[3], vup, w->p[3] );
374 Original BaseWindingForPlane() function that has serious accuracy problems. Here is why.
375 The base winding is computed as a rectangle with very large coordinates. These coordinates
376 are in the range 2^17 or 2^18. "Epsilon" (meaning the distance between two adjacent numbers)
377 at these magnitudes in 32 bit floating point world is about 0.02. So for example, if things
378 go badly (by bad luck), then the whole plane could be shifted by 0.02 units (its distance could
379 be off by that much). Then if we were to compute the winding of this plane and another of
380 the brush's planes met this winding at a very acute angle, that error could multiply to around
381 0.1 or more when computing the final vertex coordinates of the winding. 0.1 is a very large
382 error, and can lead to all sorts of disappearing triangle problems.
385 winding_t *BaseWindingForPlane( vec3_t normal, vec_t dist ){
388 vec3_t org, vright, vup;
391 // find the major axis
395 for ( i = 0 ; i < 3; i++ )
397 v = fabs( normal[i] );
404 Error( "BaseWindingForPlane: no axis found" );
407 VectorCopy( vec3_origin, vup );
419 v = DotProduct( vup, normal );
420 VectorMA( vup, -v, normal, vup );
421 VectorNormalize( vup, vup );
423 VectorScale( normal, dist, org );
425 CrossProduct( vup, normal, vright );
427 // LordHavoc: this has to use *2 because otherwise some created points may
428 // be inside the world (think of a diagonal case), and any brush with such
429 // points should be removed, failure to detect such cases is disasterous
430 VectorScale( vup, MAX_WORLD_COORD * 2, vup );
431 VectorScale( vright, MAX_WORLD_COORD * 2, vright );
433 // project a really big axis aligned box onto the plane
434 w = AllocWinding( 4 );
436 VectorSubtract( org, vright, w->p[0] );
437 VectorAdd( w->p[0], vup, w->p[0] );
439 VectorAdd( org, vright, w->p[1] );
440 VectorAdd( w->p[1], vup, w->p[1] );
442 VectorAdd( org, vright, w->p[2] );
443 VectorSubtract( w->p[2], vup, w->p[2] );
445 VectorSubtract( org, vright, w->p[3] );
446 VectorSubtract( w->p[3], vup, w->p[3] );
458 winding_t *CopyWinding( winding_t *w ){
463 Error( "CopyWinding: winding is NULL" );
466 c = AllocWinding( w->numpoints );
467 size = (size_t)( (winding_t *)NULL )->p[w->numpoints];
468 memcpy( c, w, size );
474 CopyWindingAccuIncreaseSizeAndFreeOld
477 winding_accu_t *CopyWindingAccuIncreaseSizeAndFreeOld( winding_accu_t *w ){
482 Error( "CopyWindingAccuIncreaseSizeAndFreeOld: winding is NULL" );
485 c = AllocWindingAccu( w->numpoints + 1 );
486 c->numpoints = w->numpoints;
487 for ( i = 0; i < c->numpoints; i++ )
489 VectorCopyAccu( w->p[i], c->p[i] );
491 FreeWindingAccu( w );
497 CopyWindingAccuToRegular
500 winding_t *CopyWindingAccuToRegular( winding_accu_t *w ){
505 Error( "CopyWindingAccuToRegular: winding is NULL" );
508 c = AllocWinding( w->numpoints );
509 c->numpoints = w->numpoints;
510 for ( i = 0; i < c->numpoints; i++ )
512 VectorCopyAccuToRegular( w->p[i], c->p[i] );
522 winding_t *ReverseWinding( winding_t *w ){
526 c = AllocWinding( w->numpoints );
527 for ( i = 0 ; i < w->numpoints ; i++ )
529 VectorCopy( w->p[w->numpoints - 1 - i], c->p[i] );
531 c->numpoints = w->numpoints;
541 void ClipWindingEpsilonStrict( winding_t *in, vec3_t normal, vec_t dist,
542 vec_t epsilon, winding_t **front, winding_t **back ){
543 vec_t dists[MAX_POINTS_ON_WINDING + 4];
544 int sides[MAX_POINTS_ON_WINDING + 4];
546 static vec_t dot; // VC 4.2 optimizer bug if not static
553 counts[0] = counts[1] = counts[2] = 0;
555 // determine sides for each point
556 for ( i = 0 ; i < in->numpoints ; i++ )
559 dot = DotProduct( in->p[i], normal );
562 if ( dot > epsilon ) {
563 sides[i] = SIDE_FRONT;
565 else if ( dot < -epsilon ) {
566 sides[i] = SIDE_BACK;
577 *front = *back = NULL;
579 if ( !counts[0] && !counts[1] ) {
583 *back = CopyWinding( in );
587 *front = CopyWinding( in );
591 maxpts = in->numpoints + 4; // cant use counts[0]+2 because
592 // of fp grouping errors
594 *front = f = AllocWinding( maxpts );
595 *back = b = AllocWinding( maxpts );
597 for ( i = 0 ; i < in->numpoints ; i++ )
601 if ( sides[i] == SIDE_ON ) {
602 VectorCopy( p1, f->p[f->numpoints] );
604 VectorCopy( p1, b->p[b->numpoints] );
609 if ( sides[i] == SIDE_FRONT ) {
610 VectorCopy( p1, f->p[f->numpoints] );
613 if ( sides[i] == SIDE_BACK ) {
614 VectorCopy( p1, b->p[b->numpoints] );
618 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
622 // generate a split point
623 p2 = in->p[( i + 1 ) % in->numpoints];
625 dot = dists[i] / ( dists[i] - dists[i + 1] );
626 for ( j = 0 ; j < 3 ; j++ )
627 { // avoid round off error when possible
628 if ( normal[j] == 1 ) {
631 else if ( normal[j] == -1 ) {
635 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
639 VectorCopy( mid, f->p[f->numpoints] );
641 VectorCopy( mid, b->p[b->numpoints] );
645 if ( f->numpoints > maxpts || b->numpoints > maxpts ) {
646 Error( "ClipWinding: points exceeded estimate" );
648 if ( f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING ) {
649 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
653 void ClipWindingEpsilon( winding_t *in, vec3_t normal, vec_t dist,
654 vec_t epsilon, winding_t **front, winding_t **back ){
655 ClipWindingEpsilonStrict( in, normal, dist, epsilon, front, back );
656 /* apparently most code expects that in the winding-on-plane case, the back winding is the original winding */
657 if ( !*front && !*back ) {
658 *back = CopyWinding( in );
665 ChopWindingInPlaceAccu
668 void ChopWindingInPlaceAccu( winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon ){
669 vec_accu_t fineEpsilon;
673 vec_accu_t dists[MAX_POINTS_ON_WINDING + 1];
674 int sides[MAX_POINTS_ON_WINDING + 1];
679 vec3_accu_t mid, normalAccu;
681 // We require at least a very small epsilon. It's a good idea for several reasons.
682 // First, we will be dividing by a potentially very small distance below. We don't
683 // want that distance to be too small; otherwise, things "blow up" with little accuracy
684 // due to the division. (After a second look, the value w below is in range (0,1), but
685 // graininess problem remains.) Second, Having minimum epsilon also prevents the following
686 // situation. Say for example we have a perfect octagon defined by the input winding.
687 // Say our chopping plane (defined by normal and dist) is essentially the same plane
688 // that the octagon is sitting on. Well, due to rounding errors, it may be that point
689 // 1 of the octagon might be in front, point 2 might be in back, point 3 might be in
690 // front, point 4 might be in back, and so on. So we could end up with a very ugly-
691 // looking chopped winding, and this might be undesirable, and would at least lead to
692 // a possible exhaustion of MAX_POINTS_ON_WINDING. It's better to assume that points
693 // very very close to the plane are on the plane, using an infinitesimal epsilon amount.
695 // Now, the original ChopWindingInPlace() function used a vec_t-based winding_t.
696 // So this minimum epsilon is quite similar to casting the higher resolution numbers to
697 // the lower resolution and comparing them in the lower resolution mode. We explicitly
698 // choose the minimum epsilon as something around the vec_t epsilon of one because we
699 // want the resolution of vec_accu_t to have a large resolution around the epsilon.
700 // Some of that leftover resolution even goes away after we scale to points far away.
702 // Here is a further discussion regarding the choice of smallestEpsilonAllowed.
703 // In the 32 float world (we can assume vec_t is that), the "epsilon around 1.0" is
704 // 0.00000011921. In the 64 bit float world (we can assume vec_accu_t is that), the
705 // "epsilon around 1.0" is 0.00000000000000022204. (By the way these two epsilons
706 // are defined as VEC_SMALLEST_EPSILON_AROUND_ONE VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE
707 // respectively.) If you divide the first by the second, you get approximately
708 // 536,885,246. Dividing that number by 200,000 (a typical base winding coordinate)
709 // gives 2684. So in other words, if our smallestEpsilonAllowed was chosen as exactly
710 // VEC_SMALLEST_EPSILON_AROUND_ONE, you would be guaranteed at least 2000 "ticks" in
711 // 64-bit land inside of the epsilon for all numbers we're dealing with.
713 static const vec_accu_t smallestEpsilonAllowed = ( (vec_accu_t) VEC_SMALLEST_EPSILON_AROUND_ONE ) * 0.5;
714 if ( crudeEpsilon < smallestEpsilonAllowed ) {
715 fineEpsilon = smallestEpsilonAllowed;
717 else{fineEpsilon = (vec_accu_t) crudeEpsilon; }
720 counts[0] = counts[1] = counts[2] = 0;
721 VectorCopyRegularToAccu( normal, normalAccu );
723 for ( i = 0; i < in->numpoints; i++ )
725 dists[i] = DotProductAccu( in->p[i], normalAccu ) - dist;
726 if ( dists[i] > fineEpsilon ) {
727 sides[i] = SIDE_FRONT;
729 else if ( dists[i] < -fineEpsilon ) {
730 sides[i] = SIDE_BACK;
732 else{sides[i] = SIDE_ON; }
738 // I'm wondering if whatever code that handles duplicate planes is robust enough
739 // that we never get a case where two nearly equal planes result in 2 NULL windings
740 // due to the 'if' statement below. TODO: Investigate this.
741 if ( !counts[SIDE_FRONT] ) {
742 FreeWindingAccu( in );
746 if ( !counts[SIDE_BACK] ) {
747 return; // Winding is unmodified.
750 // NOTE: The least number of points that a winding can have at this point is 2.
751 // In that case, one point is SIDE_FRONT and the other is SIDE_BACK.
753 maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
754 f = AllocWindingAccu( maxpts );
756 for ( i = 0; i < in->numpoints; i++ )
760 if ( sides[i] == SIDE_ON || sides[i] == SIDE_FRONT ) {
761 if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
762 Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
764 if ( f->numpoints >= maxpts ) { // This will probably never happen.
765 Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
766 f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
769 VectorCopyAccu( p1, f->p[f->numpoints] );
771 if ( sides[i] == SIDE_ON ) {
775 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
779 // Generate a split point.
780 p2 = in->p[( ( i + 1 ) == in->numpoints ) ? 0 : ( i + 1 )];
782 // The divisor's absolute value is greater than the dividend's absolute value.
783 // w is in the range (0,1).
784 w = dists[i] / ( dists[i] - dists[i + 1] );
786 for ( j = 0; j < 3; j++ )
788 // Avoid round-off error when possible. Check axis-aligned normal.
789 if ( normal[j] == 1 ) {
792 else if ( normal[j] == -1 ) {
795 else{mid[j] = p1[j] + ( w * ( p2[j] - p1[j] ) ); }
797 if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
798 Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
800 if ( f->numpoints >= maxpts ) { // This will probably never happen.
801 Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
802 f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
805 VectorCopyAccu( mid, f->p[f->numpoints] );
809 FreeWindingAccu( in );
818 void ChopWindingInPlace( winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon ){
820 vec_t dists[MAX_POINTS_ON_WINDING + 4];
821 int sides[MAX_POINTS_ON_WINDING + 4];
823 static vec_t dot; // VC 4.2 optimizer bug if not static
831 counts[0] = counts[1] = counts[2] = 0;
833 // determine sides for each point
834 for ( i = 0 ; i < in->numpoints ; i++ )
836 dot = DotProduct( in->p[i], normal );
839 if ( dot > epsilon ) {
840 sides[i] = SIDE_FRONT;
842 else if ( dot < -epsilon ) {
843 sides[i] = SIDE_BACK;
860 return; // inout stays the same
863 maxpts = in->numpoints + 4; // cant use counts[0]+2 because
864 // of fp grouping errors
866 f = AllocWinding( maxpts );
868 for ( i = 0 ; i < in->numpoints ; i++ )
872 if ( sides[i] == SIDE_ON ) {
873 VectorCopy( p1, f->p[f->numpoints] );
878 if ( sides[i] == SIDE_FRONT ) {
879 VectorCopy( p1, f->p[f->numpoints] );
883 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
887 // generate a split point
888 p2 = in->p[( i + 1 ) % in->numpoints];
890 dot = dists[i] / ( dists[i] - dists[i + 1] );
891 for ( j = 0 ; j < 3 ; j++ )
892 { // avoid round off error when possible
893 if ( normal[j] == 1 ) {
896 else if ( normal[j] == -1 ) {
900 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
904 VectorCopy( mid, f->p[f->numpoints] );
908 if ( f->numpoints > maxpts ) {
909 Error( "ClipWinding: points exceeded estimate" );
911 if ( f->numpoints > MAX_POINTS_ON_WINDING ) {
912 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
924 Returns the fragment of in that is on the front side
925 of the cliping plane. The original is freed.
928 winding_t *ChopWinding( winding_t *in, vec3_t normal, vec_t dist ){
931 ClipWindingEpsilon( in, normal, dist, ON_EPSILON, &f, &b );
946 void CheckWinding( winding_t *w ){
950 vec3_t dir, edgenormal, facenormal;
954 if ( w->numpoints < 3 ) {
955 Error( "CheckWinding: %i points",w->numpoints );
958 area = WindingArea( w );
960 Error( "CheckWinding: %f area", area );
963 WindingPlane( w, facenormal, &facedist );
965 for ( i = 0 ; i < w->numpoints ; i++ )
969 for ( j = 0 ; j < 3 ; j++ )
970 if ( p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD ) {
971 Error( "CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j] );
974 j = i + 1 == w->numpoints ? 0 : i + 1;
976 // check the point is on the face plane
977 d = DotProduct( p1, facenormal ) - facedist;
978 if ( d < -ON_EPSILON || d > ON_EPSILON ) {
979 Error( "CheckWinding: point off plane" );
982 // check the edge isnt degenerate
984 VectorSubtract( p2, p1, dir );
986 if ( VectorLength( dir ) < ON_EPSILON ) {
987 Error( "CheckWinding: degenerate edge" );
990 CrossProduct( facenormal, dir, edgenormal );
991 VectorNormalize( edgenormal, edgenormal );
992 edgedist = DotProduct( p1, edgenormal );
993 edgedist += ON_EPSILON;
995 // all other points must be on front side
996 for ( j = 0 ; j < w->numpoints ; j++ )
1001 d = DotProduct( w->p[j], edgenormal );
1002 if ( d > edgedist ) {
1003 Error( "CheckWinding: non-convex" );
1015 int WindingOnPlaneSide( winding_t *w, vec3_t normal, vec_t dist ){
1016 qboolean front, back;
1022 for ( i = 0 ; i < w->numpoints ; i++ )
1024 d = DotProduct( w->p[i], normal ) - dist;
1025 if ( d < -ON_EPSILON ) {
1032 if ( d > ON_EPSILON ) {
1053 AddWindingToConvexHull
1055 Both w and *hull are on the same plane
1058 #define MAX_HULL_POINTS 128
1059 void AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
1064 int numHullPoints, numNew;
1065 vec3_t hullPoints[MAX_HULL_POINTS];
1066 vec3_t newHullPoints[MAX_HULL_POINTS];
1067 vec3_t hullDirs[MAX_HULL_POINTS];
1068 qboolean hullSide[MAX_HULL_POINTS];
1072 *hull = CopyWinding( w );
1076 numHullPoints = ( *hull )->numpoints;
1077 memcpy( hullPoints, ( *hull )->p, numHullPoints * sizeof( vec3_t ) );
1079 for ( i = 0 ; i < w->numpoints ; i++ ) {
1082 // calculate hull side vectors
1083 for ( j = 0 ; j < numHullPoints ; j++ ) {
1084 k = ( j + 1 ) % numHullPoints;
1086 VectorSubtract( hullPoints[k], hullPoints[j], dir );
1087 VectorNormalize( dir, dir );
1088 CrossProduct( normal, dir, hullDirs[j] );
1092 for ( j = 0 ; j < numHullPoints ; j++ ) {
1093 VectorSubtract( p, hullPoints[j], dir );
1094 d = DotProduct( dir, hullDirs[j] );
1095 if ( d >= ON_EPSILON ) {
1098 if ( d >= -ON_EPSILON ) {
1099 hullSide[j] = qtrue;
1102 hullSide[j] = qfalse;
1106 // if the point is effectively inside, do nothing
1111 // find the back side to front side transition
1112 for ( j = 0 ; j < numHullPoints ; j++ ) {
1113 if ( !hullSide[ j % numHullPoints ] && hullSide[ ( j + 1 ) % numHullPoints ] ) {
1117 if ( j == numHullPoints ) {
1121 // insert the point here
1122 VectorCopy( p, newHullPoints[0] );
1125 // copy over all points that aren't double fronts
1126 j = ( j + 1 ) % numHullPoints;
1127 for ( k = 0 ; k < numHullPoints ; k++ ) {
1128 if ( hullSide[ ( j + k ) % numHullPoints ] && hullSide[ ( j + k + 1 ) % numHullPoints ] ) {
1131 copy = hullPoints[ ( j + k + 1 ) % numHullPoints ];
1132 VectorCopy( copy, newHullPoints[numNew] );
1136 numHullPoints = numNew;
1137 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof( vec3_t ) );
1140 FreeWinding( *hull );
1141 w = AllocWinding( numHullPoints );
1142 w->numpoints = numHullPoints;
1144 memcpy( w->p, hullPoints, numHullPoints * sizeof( vec3_t ) );