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 );
79 winding_accu_t *AllocWindingAccu( int points ){
83 if ( points >= MAX_POINTS_ON_WINDING ) {
84 Error( "AllocWindingAccu failed: MAX_POINTS_ON_WINDING exceeded" );
87 if ( numthreads == 1 ) {
88 // At the time of this writing, these statistics were not used in any way.
90 c_winding_points += points;
92 if ( c_active_windings > c_peak_windings ) {
93 c_peak_windings = c_active_windings;
96 s = sizeof( *w ) + ( points ? sizeof( w->p[0] ) * ( points - 1 ) : 0 );
107 void FreeWinding( winding_t *w ){
109 Error( "FreeWinding: winding is NULL" );
112 if ( *(unsigned *)w == 0xdeaddead ) {
113 Error( "FreeWinding: freed a freed winding" );
115 *(unsigned *)w = 0xdeaddead;
117 if ( numthreads == 1 ) {
128 void FreeWindingAccu( winding_accu_t *w ){
130 Error( "FreeWindingAccu: winding is NULL" );
133 if ( *( (unsigned *) w ) == 0xdeaddead ) {
134 Error( "FreeWindingAccu: freed a freed winding" );
136 *( (unsigned *) w ) = 0xdeaddead;
138 if ( numthreads == 1 ) {
151 void RemoveColinearPoints( winding_t *w ){
155 vec3_t p[MAX_POINTS_ON_WINDING];
158 for ( i = 0 ; i < w->numpoints ; i++ )
160 j = ( i + 1 ) % w->numpoints;
161 k = ( i + w->numpoints - 1 ) % w->numpoints;
162 VectorSubtract( w->p[j], w->p[i], v1 );
163 VectorSubtract( w->p[i], w->p[k], v2 );
164 VectorNormalize( v1,v1 );
165 VectorNormalize( v2,v2 );
166 if ( DotProduct( v1, v2 ) < 0.999 ) {
167 VectorCopy( w->p[i], p[nump] );
172 if ( nump == w->numpoints ) {
176 if ( numthreads == 1 ) {
177 c_removed += w->numpoints - nump;
180 memcpy( w->p, p, nump * sizeof( p[0] ) );
188 void WindingPlane( winding_t *w, vec3_t normal, vec_t *dist ){
191 VectorSubtract( w->p[1], w->p[0], v1 );
192 VectorSubtract( w->p[2], w->p[0], v2 );
193 CrossProduct( v2, v1, normal );
194 VectorNormalize( normal, normal );
195 *dist = DotProduct( w->p[0], normal );
204 vec_t WindingArea( winding_t *w ){
206 vec3_t d1, d2, cross;
210 for ( i = 2 ; i < w->numpoints ; i++ )
212 VectorSubtract( w->p[i - 1], w->p[0], d1 );
213 VectorSubtract( w->p[i], w->p[0], d2 );
214 CrossProduct( d1, d2, cross );
215 total += 0.5 * VectorLength( cross );
220 void WindingBounds( winding_t *w, vec3_t mins, vec3_t maxs ){
224 mins[0] = mins[1] = mins[2] = 99999;
225 maxs[0] = maxs[1] = maxs[2] = -99999;
227 for ( i = 0 ; i < w->numpoints ; i++ )
229 for ( j = 0 ; j < 3 ; j++ )
247 void WindingCenter( winding_t *w, vec3_t center ){
251 VectorCopy( vec3_origin, center );
252 for ( i = 0 ; i < w->numpoints ; i++ )
253 VectorAdd( w->p[i], center, center );
255 scale = 1.0 / w->numpoints;
256 VectorScale( center, scale, center );
261 BaseWindingForPlaneAccu
264 winding_accu_t *BaseWindingForPlaneAccu( vec3_t normal, vec_t dist ){
265 // The goal in this function is to replicate the behavior of the original BaseWindingForPlane()
266 // function (see below) but at the same time increasing accuracy substantially.
268 // The original code gave a preference for the vup vector to start out as (0, 0, 1), unless the
269 // normal had a dominant Z value, in which case vup started out as (1, 0, 0). After that, vup
270 // was "bent" [along the plane defined by normal and vup] to become perpendicular to normal.
271 // After that the vright vector was computed as the cross product of vup and normal.
273 // I'm constructing the winding polygon points in a fashion similar to the method used in the
274 // original function. Orientation is the same. The size of the winding polygon, however, is
275 // variable in this function (depending on the angle of normal), and is larger (by about a factor
276 // of 2) than the winding polygon in the original function.
280 vec3_accu_t vright, vup, org, normalAccu;
283 // One of the components of normal must have a magnitiude greater than this value,
284 // otherwise normal is not a unit vector. This is a little bit of inexpensive
285 // partial error checking we can do.
286 max = 0.56; // 1 / sqrt(1^2 + 1^2 + 1^2) = 0.577350269
289 for ( i = 0; i < 3; i++ ) {
290 v = (vec_t) fabs( normal[i] );
297 Error( "BaseWindingForPlaneAccu: no dominant axis found because normal is too short" );
301 case 0: // Fall through to next case.
303 vright[0] = (vec_accu_t) -normal[1];
304 vright[1] = (vec_accu_t) normal[0];
309 vright[1] = (vec_accu_t) -normal[2];
310 vright[2] = (vec_accu_t) normal[1];
314 // vright and normal are now perpendicular; you can prove this by taking their
315 // dot product and seeing that it's always exactly 0 (with no error).
317 // NOTE: vright is NOT a unit vector at this point. vright will have length
318 // not exceeding 1.0. The minimum length that vright can achieve happens when,
319 // for example, the Z and X components of the normal input vector are equal,
320 // and when normal's Y component is zero. In that case Z and X of the normal
321 // vector are both approximately 0.70711. The resulting vright vector in this
322 // case will have a length of 0.70711.
324 // We're relying on the fact that MAX_WORLD_COORD is a power of 2 to keep
325 // our calculation precise and relatively free of floating point error.
326 // [However, the code will still work fine if that's not the case.]
327 VectorScaleAccu( vright, ( (vec_accu_t) MAX_WORLD_COORD ) * 4.0, vright );
329 // At time time of this writing, MAX_WORLD_COORD was 65536 (2^16). Therefore
330 // the length of vright at this point is at least 185364. In comparison, a
331 // corner of the world at location (65536, 65536, 65536) is distance 113512
332 // away from the origin.
334 VectorCopyRegularToAccu( normal, normalAccu );
335 CrossProductAccu( normalAccu, vright, vup );
337 // vup now has length equal to that of vright.
339 VectorScaleAccu( normalAccu, (vec_accu_t) dist, org );
341 // org is now a point on the plane defined by normal and dist. Furthermore,
342 // org, vright, and vup are pairwise perpendicular. Now, the 4 vectors
343 // { (+-)vright + (+-)vup } have length that is at least sqrt(185364^2 + 185364^2),
344 // which is about 262144. That length lies outside the world, since the furthest
345 // point in the world has distance 113512 from the origin as mentioned above.
346 // Also, these 4 vectors are perpendicular to the org vector. So adding them
347 // to org will only increase their length. Therefore the 4 points defined below
348 // all lie outside of the world. Furthermore, it can be easily seen that the
349 // edges connecting these 4 points (in the winding_accu_t below) lie completely
350 // outside the world. sqrt(262144^2 + 262144^2)/2 = 185363, which is greater than
353 w = AllocWindingAccu( 4 );
355 VectorSubtractAccu( org, vright, w->p[0] );
356 VectorAddAccu( w->p[0], vup, w->p[0] );
358 VectorAddAccu( org, vright, w->p[1] );
359 VectorAddAccu( w->p[1], vup, w->p[1] );
361 VectorAddAccu( org, vright, w->p[2] );
362 VectorSubtractAccu( w->p[2], vup, w->p[2] );
364 VectorSubtractAccu( org, vright, w->p[3] );
365 VectorSubtractAccu( w->p[3], vup, w->p[3] );
376 Original BaseWindingForPlane() function that has serious accuracy problems. Here is why.
377 The base winding is computed as a rectangle with very large coordinates. These coordinates
378 are in the range 2^17 or 2^18. "Epsilon" (meaning the distance between two adjacent numbers)
379 at these magnitudes in 32 bit floating point world is about 0.02. So for example, if things
380 go badly (by bad luck), then the whole plane could be shifted by 0.02 units (its distance could
381 be off by that much). Then if we were to compute the winding of this plane and another of
382 the brush's planes met this winding at a very acute angle, that error could multiply to around
383 0.1 or more when computing the final vertex coordinates of the winding. 0.1 is a very large
384 error, and can lead to all sorts of disappearing triangle problems.
387 winding_t *BaseWindingForPlane( vec3_t normal, vec_t dist ){
390 vec3_t org, vright, vup;
393 // find the major axis
397 for ( i = 0 ; i < 3; i++ )
399 v = fabs( normal[i] );
406 Error( "BaseWindingForPlane: no axis found" );
409 VectorCopy( vec3_origin, vup );
421 v = DotProduct( vup, normal );
422 VectorMA( vup, -v, normal, vup );
423 VectorNormalize( vup, vup );
425 VectorScale( normal, dist, org );
427 CrossProduct( vup, normal, vright );
429 // LordHavoc: this has to use *2 because otherwise some created points may
430 // be inside the world (think of a diagonal case), and any brush with such
431 // points should be removed, failure to detect such cases is disasterous
432 VectorScale( vup, MAX_WORLD_COORD * 2, vup );
433 VectorScale( vright, MAX_WORLD_COORD * 2, vright );
435 // project a really big axis aligned box onto the plane
436 w = AllocWinding( 4 );
438 VectorSubtract( org, vright, w->p[0] );
439 VectorAdd( w->p[0], vup, w->p[0] );
441 VectorAdd( org, vright, w->p[1] );
442 VectorAdd( w->p[1], vup, w->p[1] );
444 VectorAdd( org, vright, w->p[2] );
445 VectorSubtract( w->p[2], vup, w->p[2] );
447 VectorSubtract( org, vright, w->p[3] );
448 VectorSubtract( w->p[3], vup, w->p[3] );
460 winding_t *CopyWinding( winding_t *w ){
465 Error( "CopyWinding: winding is NULL" );
468 c = AllocWinding( w->numpoints );
469 size = (size_t)( (winding_t *)NULL )->p[w->numpoints];
470 memcpy( c, w, size );
476 CopyWindingAccuIncreaseSizeAndFreeOld
479 winding_accu_t *CopyWindingAccuIncreaseSizeAndFreeOld( winding_accu_t *w ){
484 Error( "CopyWindingAccuIncreaseSizeAndFreeOld: winding is NULL" );
487 c = AllocWindingAccu( w->numpoints + 1 );
488 c->numpoints = w->numpoints;
489 for ( i = 0; i < c->numpoints; i++ )
491 VectorCopyAccu( w->p[i], c->p[i] );
493 FreeWindingAccu( w );
499 CopyWindingAccuToRegular
502 winding_t *CopyWindingAccuToRegular( winding_accu_t *w ){
507 Error( "CopyWindingAccuToRegular: winding is NULL" );
510 c = AllocWinding( w->numpoints );
511 c->numpoints = w->numpoints;
512 for ( i = 0; i < c->numpoints; i++ )
514 VectorCopyAccuToRegular( w->p[i], c->p[i] );
524 winding_t *ReverseWinding( winding_t *w ){
528 c = AllocWinding( w->numpoints );
529 for ( i = 0 ; i < w->numpoints ; i++ )
531 VectorCopy( w->p[w->numpoints - 1 - i], c->p[i] );
533 c->numpoints = w->numpoints;
543 void ClipWindingEpsilonStrict( winding_t *in, vec3_t normal, vec_t dist,
544 vec_t epsilon, winding_t **front, winding_t **back ){
545 vec_t dists[MAX_POINTS_ON_WINDING + 4];
546 int sides[MAX_POINTS_ON_WINDING + 4];
548 static vec_t dot; // VC 4.2 optimizer bug if not static
555 counts[0] = counts[1] = counts[2] = 0;
557 // determine sides for each point
558 for ( i = 0 ; i < in->numpoints ; i++ )
561 dot = DotProduct( in->p[i], normal );
564 if ( dot > epsilon ) {
565 sides[i] = SIDE_FRONT;
567 else if ( dot < -epsilon ) {
568 sides[i] = SIDE_BACK;
579 *front = *back = NULL;
581 if ( !counts[0] && !counts[1] ) {
585 *back = CopyWinding( in );
589 *front = CopyWinding( in );
593 maxpts = in->numpoints + 4; // cant use counts[0]+2 because
594 // of fp grouping errors
596 *front = f = AllocWinding( maxpts );
597 *back = b = AllocWinding( maxpts );
599 for ( i = 0 ; i < in->numpoints ; i++ )
603 if ( sides[i] == SIDE_ON ) {
604 VectorCopy( p1, f->p[f->numpoints] );
606 VectorCopy( p1, b->p[b->numpoints] );
611 if ( sides[i] == SIDE_FRONT ) {
612 VectorCopy( p1, f->p[f->numpoints] );
615 if ( sides[i] == SIDE_BACK ) {
616 VectorCopy( p1, b->p[b->numpoints] );
620 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
624 // generate a split point
625 p2 = in->p[( i + 1 ) % in->numpoints];
627 dot = dists[i] / ( dists[i] - dists[i + 1] );
628 for ( j = 0 ; j < 3 ; j++ )
629 { // avoid round off error when possible
630 if ( normal[j] == 1 ) {
633 else if ( normal[j] == -1 ) {
637 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
641 VectorCopy( mid, f->p[f->numpoints] );
643 VectorCopy( mid, b->p[b->numpoints] );
647 if ( f->numpoints > maxpts || b->numpoints > maxpts ) {
648 Error( "ClipWinding: points exceeded estimate" );
650 if ( f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING ) {
651 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
655 void ClipWindingEpsilon( winding_t *in, vec3_t normal, vec_t dist,
656 vec_t epsilon, winding_t **front, winding_t **back ){
657 ClipWindingEpsilonStrict( in, normal, dist, epsilon, front, back );
658 /* apparently most code expects that in the winding-on-plane case, the back winding is the original winding */
659 if ( !*front && !*back ) {
660 *back = CopyWinding( in );
667 ChopWindingInPlaceAccu
670 void ChopWindingInPlaceAccu( winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon ){
671 vec_accu_t fineEpsilon;
675 vec_accu_t dists[MAX_POINTS_ON_WINDING + 1];
676 int sides[MAX_POINTS_ON_WINDING + 1];
681 vec3_accu_t mid, normalAccu;
683 // We require at least a very small epsilon. It's a good idea for several reasons.
684 // First, we will be dividing by a potentially very small distance below. We don't
685 // want that distance to be too small; otherwise, things "blow up" with little accuracy
686 // due to the division. (After a second look, the value w below is in range (0,1), but
687 // graininess problem remains.) Second, Having minimum epsilon also prevents the following
688 // situation. Say for example we have a perfect octagon defined by the input winding.
689 // Say our chopping plane (defined by normal and dist) is essentially the same plane
690 // that the octagon is sitting on. Well, due to rounding errors, it may be that point
691 // 1 of the octagon might be in front, point 2 might be in back, point 3 might be in
692 // front, point 4 might be in back, and so on. So we could end up with a very ugly-
693 // looking chopped winding, and this might be undesirable, and would at least lead to
694 // a possible exhaustion of MAX_POINTS_ON_WINDING. It's better to assume that points
695 // very very close to the plane are on the plane, using an infinitesimal epsilon amount.
697 // Now, the original ChopWindingInPlace() function used a vec_t-based winding_t.
698 // So this minimum epsilon is quite similar to casting the higher resolution numbers to
699 // the lower resolution and comparing them in the lower resolution mode. We explicitly
700 // choose the minimum epsilon as something around the vec_t epsilon of one because we
701 // want the resolution of vec_accu_t to have a large resolution around the epsilon.
702 // Some of that leftover resolution even goes away after we scale to points far away.
704 // Here is a further discussion regarding the choice of smallestEpsilonAllowed.
705 // In the 32 float world (we can assume vec_t is that), the "epsilon around 1.0" is
706 // 0.00000011921. In the 64 bit float world (we can assume vec_accu_t is that), the
707 // "epsilon around 1.0" is 0.00000000000000022204. (By the way these two epsilons
708 // are defined as VEC_SMALLEST_EPSILON_AROUND_ONE VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE
709 // respectively.) If you divide the first by the second, you get approximately
710 // 536,885,246. Dividing that number by 200,000 (a typical base winding coordinate)
711 // gives 2684. So in other words, if our smallestEpsilonAllowed was chosen as exactly
712 // VEC_SMALLEST_EPSILON_AROUND_ONE, you would be guaranteed at least 2000 "ticks" in
713 // 64-bit land inside of the epsilon for all numbers we're dealing with.
715 static const vec_accu_t smallestEpsilonAllowed = ( (vec_accu_t) VEC_SMALLEST_EPSILON_AROUND_ONE ) * 0.5;
716 if ( crudeEpsilon < smallestEpsilonAllowed ) {
717 fineEpsilon = smallestEpsilonAllowed;
719 else{fineEpsilon = (vec_accu_t) crudeEpsilon; }
722 counts[0] = counts[1] = counts[2] = 0;
723 VectorCopyRegularToAccu( normal, normalAccu );
725 for ( i = 0; i < in->numpoints; i++ )
727 dists[i] = DotProductAccu( in->p[i], normalAccu ) - dist;
728 if ( dists[i] > fineEpsilon ) {
729 sides[i] = SIDE_FRONT;
731 else if ( dists[i] < -fineEpsilon ) {
732 sides[i] = SIDE_BACK;
734 else{sides[i] = SIDE_ON; }
740 // I'm wondering if whatever code that handles duplicate planes is robust enough
741 // that we never get a case where two nearly equal planes result in 2 NULL windings
742 // due to the 'if' statement below. TODO: Investigate this.
743 if ( !counts[SIDE_FRONT] ) {
744 FreeWindingAccu( in );
748 if ( !counts[SIDE_BACK] ) {
749 return; // Winding is unmodified.
752 // NOTE: The least number of points that a winding can have at this point is 2.
753 // In that case, one point is SIDE_FRONT and the other is SIDE_BACK.
755 maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
756 f = AllocWindingAccu( maxpts );
758 for ( i = 0; i < in->numpoints; i++ )
762 if ( sides[i] == SIDE_ON || sides[i] == SIDE_FRONT ) {
763 if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
764 Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
766 if ( f->numpoints >= maxpts ) { // This will probably never happen.
767 Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
768 f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
771 VectorCopyAccu( p1, f->p[f->numpoints] );
773 if ( sides[i] == SIDE_ON ) {
777 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
781 // Generate a split point.
782 p2 = in->p[( ( i + 1 ) == in->numpoints ) ? 0 : ( i + 1 )];
784 // The divisor's absolute value is greater than the dividend's absolute value.
785 // w is in the range (0,1).
786 w = dists[i] / ( dists[i] - dists[i + 1] );
788 for ( j = 0; j < 3; j++ )
790 // Avoid round-off error when possible. Check axis-aligned normal.
791 if ( normal[j] == 1 ) {
794 else if ( normal[j] == -1 ) {
797 else{mid[j] = p1[j] + ( w * ( p2[j] - p1[j] ) ); }
799 if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
800 Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
802 if ( f->numpoints >= maxpts ) { // This will probably never happen.
803 Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
804 f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
807 VectorCopyAccu( mid, f->p[f->numpoints] );
811 FreeWindingAccu( in );
820 void ChopWindingInPlace( winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon ){
822 vec_t dists[MAX_POINTS_ON_WINDING + 4];
823 int sides[MAX_POINTS_ON_WINDING + 4];
825 static vec_t dot; // VC 4.2 optimizer bug if not static
833 counts[0] = counts[1] = counts[2] = 0;
835 // determine sides for each point
836 for ( i = 0 ; i < in->numpoints ; i++ )
838 dot = DotProduct( in->p[i], normal );
841 if ( dot > epsilon ) {
842 sides[i] = SIDE_FRONT;
844 else if ( dot < -epsilon ) {
845 sides[i] = SIDE_BACK;
862 return; // inout stays the same
865 maxpts = in->numpoints + 4; // cant use counts[0]+2 because
866 // of fp grouping errors
868 f = AllocWinding( maxpts );
870 for ( i = 0 ; i < in->numpoints ; i++ )
874 if ( sides[i] == SIDE_ON ) {
875 VectorCopy( p1, f->p[f->numpoints] );
880 if ( sides[i] == SIDE_FRONT ) {
881 VectorCopy( p1, f->p[f->numpoints] );
885 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
889 // generate a split point
890 p2 = in->p[( i + 1 ) % in->numpoints];
892 dot = dists[i] / ( dists[i] - dists[i + 1] );
893 for ( j = 0 ; j < 3 ; j++ )
894 { // avoid round off error when possible
895 if ( normal[j] == 1 ) {
898 else if ( normal[j] == -1 ) {
902 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
906 VectorCopy( mid, f->p[f->numpoints] );
910 if ( f->numpoints > maxpts ) {
911 Error( "ClipWinding: points exceeded estimate" );
913 if ( f->numpoints > MAX_POINTS_ON_WINDING ) {
914 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
926 Returns the fragment of in that is on the front side
927 of the cliping plane. The original is freed.
930 winding_t *ChopWinding( winding_t *in, vec3_t normal, vec_t dist ){
933 ClipWindingEpsilon( in, normal, dist, ON_EPSILON, &f, &b );
948 void CheckWinding( winding_t *w ){
952 vec3_t dir, edgenormal, facenormal;
956 if ( w->numpoints < 3 ) {
957 Error( "CheckWinding: %i points",w->numpoints );
960 area = WindingArea( w );
962 Error( "CheckWinding: %f area", area );
965 WindingPlane( w, facenormal, &facedist );
967 for ( i = 0 ; i < w->numpoints ; i++ )
971 for ( j = 0 ; j < 3 ; j++ )
972 if ( p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD ) {
973 Error( "CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j] );
976 j = i + 1 == w->numpoints ? 0 : i + 1;
978 // check the point is on the face plane
979 d = DotProduct( p1, facenormal ) - facedist;
980 if ( d < -ON_EPSILON || d > ON_EPSILON ) {
981 Error( "CheckWinding: point off plane" );
984 // check the edge isnt degenerate
986 VectorSubtract( p2, p1, dir );
988 if ( VectorLength( dir ) < ON_EPSILON ) {
989 Error( "CheckWinding: degenerate edge" );
992 CrossProduct( facenormal, dir, edgenormal );
993 VectorNormalize( edgenormal, edgenormal );
994 edgedist = DotProduct( p1, edgenormal );
995 edgedist += ON_EPSILON;
997 // all other points must be on front side
998 for ( j = 0 ; j < w->numpoints ; j++ )
1003 d = DotProduct( w->p[j], edgenormal );
1004 if ( d > edgedist ) {
1005 Error( "CheckWinding: non-convex" );
1017 int WindingOnPlaneSide( winding_t *w, vec3_t normal, vec_t dist ){
1018 qboolean front, back;
1024 for ( i = 0 ; i < w->numpoints ; i++ )
1026 d = DotProduct( w->p[i], normal ) - dist;
1027 if ( d < -ON_EPSILON ) {
1034 if ( d > ON_EPSILON ) {
1055 AddWindingToConvexHull
1057 Both w and *hull are on the same plane
1060 #define MAX_HULL_POINTS 128
1061 void AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
1066 int numHullPoints, numNew;
1067 vec3_t hullPoints[MAX_HULL_POINTS];
1068 vec3_t newHullPoints[MAX_HULL_POINTS];
1069 vec3_t hullDirs[MAX_HULL_POINTS];
1070 qboolean hullSide[MAX_HULL_POINTS];
1074 *hull = CopyWinding( w );
1078 numHullPoints = ( *hull )->numpoints;
1079 memcpy( hullPoints, ( *hull )->p, numHullPoints * sizeof( vec3_t ) );
1081 for ( i = 0 ; i < w->numpoints ; i++ ) {
1084 // calculate hull side vectors
1085 for ( j = 0 ; j < numHullPoints ; j++ ) {
1086 k = ( j + 1 ) % numHullPoints;
1088 VectorSubtract( hullPoints[k], hullPoints[j], dir );
1089 VectorNormalize( dir, dir );
1090 CrossProduct( normal, dir, hullDirs[j] );
1094 for ( j = 0 ; j < numHullPoints ; j++ ) {
1095 VectorSubtract( p, hullPoints[j], dir );
1096 d = DotProduct( dir, hullDirs[j] );
1097 if ( d >= ON_EPSILON ) {
1100 if ( d >= -ON_EPSILON ) {
1101 hullSide[j] = qtrue;
1104 hullSide[j] = qfalse;
1108 // if the point is effectively inside, do nothing
1113 // find the back side to front side transition
1114 for ( j = 0 ; j < numHullPoints ; j++ ) {
1115 if ( !hullSide[ j % numHullPoints ] && hullSide[ ( j + 1 ) % numHullPoints ] ) {
1119 if ( j == numHullPoints ) {
1123 // insert the point here
1124 VectorCopy( p, newHullPoints[0] );
1127 // copy over all points that aren't double fronts
1128 j = ( j + 1 ) % numHullPoints;
1129 for ( k = 0 ; k < numHullPoints ; k++ ) {
1130 if ( hullSide[ ( j + k ) % numHullPoints ] && hullSide[ ( j + k + 1 ) % numHullPoints ] ) {
1133 copy = hullPoints[ ( j + k + 1 ) % numHullPoints ];
1134 VectorCopy( copy, newHullPoints[numNew] );
1138 numHullPoints = numNew;
1139 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof( vec3_t ) );
1142 FreeWinding( *hull );
1143 w = AllocWinding( numHullPoints );
1144 w->numpoints = numHullPoints;
1146 memcpy( w->p, hullPoints, numHullPoints * sizeof( vec3_t ) );