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
30 extern int numthreads;
32 // counters are only bumped when running single threaded,
33 // because they are an awefull coherence problem
34 int c_active_windings;
39 #define BOGUS_RANGE WORLD_SIZE
41 void pw( winding_t *w ){
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] );
53 winding_t *AllocWinding( int points ){
57 if ( points >= MAX_POINTS_ON_WINDING ) {
58 Error( "AllocWinding failed: MAX_POINTS_ON_WINDING exceeded" );
61 if ( numthreads == 1 ) {
63 c_winding_points += points;
65 if ( c_active_windings > c_peak_windings ) {
66 c_peak_windings = c_active_windings;
69 s = sizeof( vec_t ) * 3 * points + sizeof( int );
80 winding_accu_t *AllocWindingAccu( int points ){
84 if ( points >= MAX_POINTS_ON_WINDING ) {
85 Error( "AllocWindingAccu failed: MAX_POINTS_ON_WINDING exceeded" );
88 if ( numthreads == 1 ) {
89 // At the time of this writing, these statistics were not used in any way.
91 c_winding_points += points;
93 if ( c_active_windings > c_peak_windings ) {
94 c_peak_windings = c_active_windings;
97 s = sizeof( vec_accu_t ) * 3 * points + sizeof( int );
108 void FreeWinding( winding_t *w ){
110 Error( "FreeWinding: winding is NULL" );
113 if ( *(unsigned *)w == 0xdeaddead ) {
114 Error( "FreeWinding: freed a freed winding" );
116 *(unsigned *)w = 0xdeaddead;
118 if ( numthreads == 1 ) {
129 void FreeWindingAccu( winding_accu_t *w ){
131 Error( "FreeWindingAccu: winding is NULL" );
134 if ( *( (unsigned *) w ) == 0xdeaddead ) {
135 Error( "FreeWindingAccu: freed a freed winding" );
137 *( (unsigned *) w ) = 0xdeaddead;
139 if ( numthreads == 1 ) {
152 void RemoveColinearPoints( winding_t *w ){
156 vec3_t p[MAX_POINTS_ON_WINDING];
159 for ( i = 0 ; i < w->numpoints ; i++ )
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] );
173 if ( nump == w->numpoints ) {
177 if ( numthreads == 1 ) {
178 c_removed += w->numpoints - nump;
181 memcpy( w->p, p, nump * sizeof( p[0] ) );
189 void WindingPlane( winding_t *w, vec3_t normal, vec_t *dist ){
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 );
205 vec_t WindingArea( winding_t *w ){
207 vec3_t d1, d2, cross;
211 for ( i = 2 ; i < w->numpoints ; i++ )
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 );
221 void WindingBounds( winding_t *w, vec3_t mins, vec3_t maxs ){
225 mins[0] = mins[1] = mins[2] = 99999;
226 maxs[0] = maxs[1] = maxs[2] = -99999;
228 for ( i = 0 ; i < w->numpoints ; i++ )
230 for ( j = 0 ; j < 3 ; j++ )
248 void WindingCenter( winding_t *w, vec3_t center ){
252 VectorCopy( vec3_origin, center );
253 for ( i = 0 ; i < w->numpoints ; i++ )
254 VectorAdd( w->p[i], center, center );
256 scale = 1.0 / w->numpoints;
257 VectorScale( center, scale, center );
262 BaseWindingForPlaneAccu
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.
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.
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.
281 vec3_accu_t vright, vup, org, normalAccu;
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
290 for ( i = 0; i < 3; i++ ) {
291 v = (vec_t) fabs( normal[i] );
298 Error( "BaseWindingForPlaneAccu: no dominant axis found because normal is too short" );
302 case 0: // Fall through to next case.
304 vright[0] = (vec_accu_t) -normal[1];
305 vright[1] = (vec_accu_t) normal[0];
310 vright[1] = (vec_accu_t) -normal[2];
311 vright[2] = (vec_accu_t) normal[1];
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).
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.
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 );
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.
335 VectorCopyRegularToAccu( normal, normalAccu );
336 CrossProductAccu( normalAccu, vright, vup );
338 // vup now has length equal to that of vright.
340 VectorScaleAccu( normalAccu, (vec_accu_t) dist, org );
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
354 w = AllocWindingAccu( 4 );
356 VectorSubtractAccu( org, vright, w->p[0] );
357 VectorAddAccu( w->p[0], vup, w->p[0] );
359 VectorAddAccu( org, vright, w->p[1] );
360 VectorAddAccu( w->p[1], vup, w->p[1] );
362 VectorAddAccu( org, vright, w->p[2] );
363 VectorSubtractAccu( w->p[2], vup, w->p[2] );
365 VectorSubtractAccu( org, vright, w->p[3] );
366 VectorSubtractAccu( w->p[3], vup, w->p[3] );
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.
388 winding_t *BaseWindingForPlane( vec3_t normal, vec_t dist ){
391 vec3_t org, vright, vup;
394 // find the major axis
398 for ( i = 0 ; i < 3; i++ )
400 v = fabs( normal[i] );
407 Error( "BaseWindingForPlane: no axis found" );
410 VectorCopy( vec3_origin, vup );
422 v = DotProduct( vup, normal );
423 VectorMA( vup, -v, normal, vup );
424 VectorNormalize( vup, vup );
426 VectorScale( normal, dist, org );
428 CrossProduct( vup, normal, vright );
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 );
436 // project a really big axis aligned box onto the plane
437 w = AllocWinding( 4 );
439 VectorSubtract( org, vright, w->p[0] );
440 VectorAdd( w->p[0], vup, w->p[0] );
442 VectorAdd( org, vright, w->p[1] );
443 VectorAdd( w->p[1], vup, w->p[1] );
445 VectorAdd( org, vright, w->p[2] );
446 VectorSubtract( w->p[2], vup, w->p[2] );
448 VectorSubtract( org, vright, w->p[3] );
449 VectorSubtract( w->p[3], vup, w->p[3] );
461 winding_t *CopyWinding( winding_t *w ){
466 Error( "CopyWinding: winding is NULL" );
469 c = AllocWinding( w->numpoints );
470 size = (int)( (size_t)( (winding_t *)0 )->p[w->numpoints] );
471 memcpy( c, w, size );
477 CopyWindingAccuIncreaseSizeAndFreeOld
480 winding_accu_t *CopyWindingAccuIncreaseSizeAndFreeOld( winding_accu_t *w ){
485 Error( "CopyWindingAccuIncreaseSizeAndFreeOld: winding is NULL" );
488 c = AllocWindingAccu( w->numpoints + 1 );
489 c->numpoints = w->numpoints;
490 for ( i = 0; i < c->numpoints; i++ )
492 VectorCopyAccu( w->p[i], c->p[i] );
494 FreeWindingAccu( w );
500 CopyWindingAccuToRegular
503 winding_t *CopyWindingAccuToRegular( winding_accu_t *w ){
508 Error( "CopyWindingAccuToRegular: winding is NULL" );
511 c = AllocWinding( w->numpoints );
512 c->numpoints = w->numpoints;
513 for ( i = 0; i < c->numpoints; i++ )
515 VectorCopyAccuToRegular( w->p[i], c->p[i] );
525 winding_t *ReverseWinding( winding_t *w ){
529 c = AllocWinding( w->numpoints );
530 for ( i = 0 ; i < w->numpoints ; i++ )
532 VectorCopy( w->p[w->numpoints - 1 - i], c->p[i] );
534 c->numpoints = w->numpoints;
544 void ClipWindingEpsilon( 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];
549 static vec_t dot; // VC 4.2 optimizer bug if not static
556 counts[0] = counts[1] = counts[2] = 0;
558 // determine sides for each point
559 for ( i = 0 ; i < in->numpoints ; i++ )
562 dot = DotProduct( in->p[i], normal );
565 if ( dot > epsilon ) {
566 sides[i] = SIDE_FRONT;
568 else if ( dot < -epsilon ) {
569 sides[i] = SIDE_BACK;
580 *front = *back = NULL;
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" );
656 ChopWindingInPlaceAccu
659 void ChopWindingInPlaceAccu( winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon ){
660 vec_accu_t fineEpsilon;
664 vec_accu_t dists[MAX_POINTS_ON_WINDING + 1];
665 int sides[MAX_POINTS_ON_WINDING + 1];
670 vec3_accu_t mid, normalAccu;
672 // We require at least a very small epsilon. It's a good idea for several reasons.
673 // First, we will be dividing by a potentially very small distance below. We don't
674 // want that distance to be too small; otherwise, things "blow up" with little accuracy
675 // due to the division. (After a second look, the value w below is in range (0,1), but
676 // graininess problem remains.) Second, Having minimum epsilon also prevents the following
677 // situation. Say for example we have a perfect octagon defined by the input winding.
678 // Say our chopping plane (defined by normal and dist) is essentially the same plane
679 // that the octagon is sitting on. Well, due to rounding errors, it may be that point
680 // 1 of the octagon might be in front, point 2 might be in back, point 3 might be in
681 // front, point 4 might be in back, and so on. So we could end up with a very ugly-
682 // looking chopped winding, and this might be undesirable, and would at least lead to
683 // a possible exhaustion of MAX_POINTS_ON_WINDING. It's better to assume that points
684 // very very close to the plane are on the plane, using an infinitesimal epsilon amount.
686 // Now, the original ChopWindingInPlace() function used a vec_t-based winding_t.
687 // So this minimum epsilon is quite similar to casting the higher resolution numbers to
688 // the lower resolution and comparing them in the lower resolution mode. We explicitly
689 // choose the minimum epsilon as something around the vec_t epsilon of one because we
690 // want the resolution of vec_accu_t to have a large resolution around the epsilon.
691 // Some of that leftover resolution even goes away after we scale to points far away.
693 // Here is a further discussion regarding the choice of smallestEpsilonAllowed.
694 // In the 32 float world (we can assume vec_t is that), the "epsilon around 1.0" is
695 // 0.00000011921. In the 64 bit float world (we can assume vec_accu_t is that), the
696 // "epsilon around 1.0" is 0.00000000000000022204. (By the way these two epsilons
697 // are defined as VEC_SMALLEST_EPSILON_AROUND_ONE VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE
698 // respectively.) If you divide the first by the second, you get approximately
699 // 536,885,246. Dividing that number by 200,000 (a typical base winding coordinate)
700 // gives 2684. So in other words, if our smallestEpsilonAllowed was chosen as exactly
701 // VEC_SMALLEST_EPSILON_AROUND_ONE, you would be guaranteed at least 2000 "ticks" in
702 // 64-bit land inside of the epsilon for all numbers we're dealing with.
704 static const vec_accu_t smallestEpsilonAllowed = ( (vec_accu_t) VEC_SMALLEST_EPSILON_AROUND_ONE ) * 0.5;
705 if ( crudeEpsilon < smallestEpsilonAllowed ) {
706 fineEpsilon = smallestEpsilonAllowed;
708 else{fineEpsilon = (vec_accu_t) crudeEpsilon; }
711 counts[0] = counts[1] = counts[2] = 0;
712 VectorCopyRegularToAccu( normal, normalAccu );
714 for ( i = 0; i < in->numpoints; i++ )
716 dists[i] = DotProductAccu( in->p[i], normalAccu ) - dist;
717 if ( dists[i] > fineEpsilon ) {
718 sides[i] = SIDE_FRONT;
720 else if ( dists[i] < -fineEpsilon ) {
721 sides[i] = SIDE_BACK;
723 else{sides[i] = SIDE_ON; }
729 // I'm wondering if whatever code that handles duplicate planes is robust enough
730 // that we never get a case where two nearly equal planes result in 2 NULL windings
731 // due to the 'if' statement below. TODO: Investigate this.
732 if ( !counts[SIDE_FRONT] ) {
733 FreeWindingAccu( in );
737 if ( !counts[SIDE_BACK] ) {
738 return; // Winding is unmodified.
741 // NOTE: The least number of points that a winding can have at this point is 2.
742 // In that case, one point is SIDE_FRONT and the other is SIDE_BACK.
744 maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
745 f = AllocWindingAccu( maxpts );
747 for ( i = 0; i < in->numpoints; i++ )
751 if ( sides[i] == SIDE_ON || sides[i] == SIDE_FRONT ) {
752 if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
753 Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
755 if ( f->numpoints >= maxpts ) { // This will probably never happen.
756 Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
757 f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
760 VectorCopyAccu( p1, f->p[f->numpoints] );
762 if ( sides[i] == SIDE_ON ) {
766 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
770 // Generate a split point.
771 p2 = in->p[( ( i + 1 ) == in->numpoints ) ? 0 : ( i + 1 )];
773 // The divisor's absolute value is greater than the dividend's absolute value.
774 // w is in the range (0,1).
775 w = dists[i] / ( dists[i] - dists[i + 1] );
777 for ( j = 0; j < 3; j++ )
779 // Avoid round-off error when possible. Check axis-aligned normal.
780 if ( normal[j] == 1 ) {
783 else if ( normal[j] == -1 ) {
786 else{mid[j] = p1[j] + ( w * ( p2[j] - p1[j] ) ); }
788 if ( f->numpoints >= MAX_POINTS_ON_WINDING ) {
789 Error( "ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING" );
791 if ( f->numpoints >= maxpts ) { // This will probably never happen.
792 Sys_FPrintf( SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n" );
793 f = CopyWindingAccuIncreaseSizeAndFreeOld( f );
796 VectorCopyAccu( mid, f->p[f->numpoints] );
800 FreeWindingAccu( in );
809 void ChopWindingInPlace( winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon ){
811 vec_t dists[MAX_POINTS_ON_WINDING + 4];
812 int sides[MAX_POINTS_ON_WINDING + 4];
814 static vec_t dot; // VC 4.2 optimizer bug if not static
822 counts[0] = counts[1] = counts[2] = 0;
824 // determine sides for each point
825 for ( i = 0 ; i < in->numpoints ; i++ )
827 dot = DotProduct( in->p[i], normal );
830 if ( dot > epsilon ) {
831 sides[i] = SIDE_FRONT;
833 else if ( dot < -epsilon ) {
834 sides[i] = SIDE_BACK;
851 return; // inout stays the same
854 maxpts = in->numpoints + 4; // cant use counts[0]+2 because
855 // of fp grouping errors
857 f = AllocWinding( maxpts );
859 for ( i = 0 ; i < in->numpoints ; i++ )
863 if ( sides[i] == SIDE_ON ) {
864 VectorCopy( p1, f->p[f->numpoints] );
869 if ( sides[i] == SIDE_FRONT ) {
870 VectorCopy( p1, f->p[f->numpoints] );
874 if ( sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i] ) {
878 // generate a split point
879 p2 = in->p[( i + 1 ) % in->numpoints];
881 dot = dists[i] / ( dists[i] - dists[i + 1] );
882 for ( j = 0 ; j < 3 ; j++ )
883 { // avoid round off error when possible
884 if ( normal[j] == 1 ) {
887 else if ( normal[j] == -1 ) {
891 mid[j] = p1[j] + dot * ( p2[j] - p1[j] );
895 VectorCopy( mid, f->p[f->numpoints] );
899 if ( f->numpoints > maxpts ) {
900 Error( "ClipWinding: points exceeded estimate" );
902 if ( f->numpoints > MAX_POINTS_ON_WINDING ) {
903 Error( "ClipWinding: MAX_POINTS_ON_WINDING" );
915 Returns the fragment of in that is on the front side
916 of the cliping plane. The original is freed.
919 winding_t *ChopWinding( winding_t *in, vec3_t normal, vec_t dist ){
922 ClipWindingEpsilon( in, normal, dist, ON_EPSILON, &f, &b );
937 void CheckWinding( winding_t *w ){
941 vec3_t dir, edgenormal, facenormal;
945 if ( w->numpoints < 3 ) {
946 Error( "CheckWinding: %i points",w->numpoints );
949 area = WindingArea( w );
951 Error( "CheckWinding: %f area", area );
954 WindingPlane( w, facenormal, &facedist );
956 for ( i = 0 ; i < w->numpoints ; i++ )
960 for ( j = 0 ; j < 3 ; j++ )
961 if ( p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD ) {
962 Error( "CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j] );
965 j = i + 1 == w->numpoints ? 0 : i + 1;
967 // check the point is on the face plane
968 d = DotProduct( p1, facenormal ) - facedist;
969 if ( d < -ON_EPSILON || d > ON_EPSILON ) {
970 Error( "CheckWinding: point off plane" );
973 // check the edge isnt degenerate
975 VectorSubtract( p2, p1, dir );
977 if ( VectorLength( dir ) < ON_EPSILON ) {
978 Error( "CheckWinding: degenerate edge" );
981 CrossProduct( facenormal, dir, edgenormal );
982 VectorNormalize( edgenormal, edgenormal );
983 edgedist = DotProduct( p1, edgenormal );
984 edgedist += ON_EPSILON;
986 // all other points must be on front side
987 for ( j = 0 ; j < w->numpoints ; j++ )
992 d = DotProduct( w->p[j], edgenormal );
993 if ( d > edgedist ) {
994 Error( "CheckWinding: non-convex" );
1006 int WindingOnPlaneSide( winding_t *w, vec3_t normal, vec_t dist ){
1007 qboolean front, back;
1013 for ( i = 0 ; i < w->numpoints ; i++ )
1015 d = DotProduct( w->p[i], normal ) - dist;
1016 if ( d < -ON_EPSILON ) {
1023 if ( d > ON_EPSILON ) {
1044 AddWindingToConvexHull
1046 Both w and *hull are on the same plane
1049 #define MAX_HULL_POINTS 128
1050 void AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
1055 int numHullPoints, numNew;
1056 vec3_t hullPoints[MAX_HULL_POINTS];
1057 vec3_t newHullPoints[MAX_HULL_POINTS];
1058 vec3_t hullDirs[MAX_HULL_POINTS];
1059 qboolean hullSide[MAX_HULL_POINTS];
1063 *hull = CopyWinding( w );
1067 numHullPoints = ( *hull )->numpoints;
1068 memcpy( hullPoints, ( *hull )->p, numHullPoints * sizeof( vec3_t ) );
1070 for ( i = 0 ; i < w->numpoints ; i++ ) {
1073 // calculate hull side vectors
1074 for ( j = 0 ; j < numHullPoints ; j++ ) {
1075 k = ( j + 1 ) % numHullPoints;
1077 VectorSubtract( hullPoints[k], hullPoints[j], dir );
1078 VectorNormalize( dir, dir );
1079 CrossProduct( normal, dir, hullDirs[j] );
1083 for ( j = 0 ; j < numHullPoints ; j++ ) {
1084 VectorSubtract( p, hullPoints[j], dir );
1085 d = DotProduct( dir, hullDirs[j] );
1086 if ( d >= ON_EPSILON ) {
1089 if ( d >= -ON_EPSILON ) {
1090 hullSide[j] = qtrue;
1093 hullSide[j] = qfalse;
1097 // if the point is effectively inside, do nothing
1102 // find the back side to front side transition
1103 for ( j = 0 ; j < numHullPoints ; j++ ) {
1104 if ( !hullSide[ j % numHullPoints ] && hullSide[ ( j + 1 ) % numHullPoints ] ) {
1108 if ( j == numHullPoints ) {
1112 // insert the point here
1113 VectorCopy( p, newHullPoints[0] );
1116 // copy over all points that aren't double fronts
1117 j = ( j + 1 ) % numHullPoints;
1118 for ( k = 0 ; k < numHullPoints ; k++ ) {
1119 if ( hullSide[ ( j + k ) % numHullPoints ] && hullSide[ ( j + k + 1 ) % numHullPoints ] ) {
1122 copy = hullPoints[ ( j + k + 1 ) % numHullPoints ];
1123 VectorCopy( copy, newHullPoints[numNew] );
1127 numHullPoints = numNew;
1128 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof( vec3_t ) );
1131 FreeWinding( *hull );
1132 w = AllocWinding( numHullPoints );
1133 w->numpoints = numHullPoints;
1135 memcpy( w->p, hullPoints, numHullPoints * sizeof( vec3_t ) );