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
44 for (i=0 ; i<w->numpoints ; i++)
45 Sys_Printf ("(%5.1f, %5.1f, %5.1f)\n",w->p[i][0], w->p[i][1],w->p[i][2]);
54 winding_t *AllocWinding (int points)
59 if (points >= MAX_POINTS_ON_WINDING)
60 Error ("AllocWinding failed: MAX_POINTS_ON_WINDING exceeded");
65 c_winding_points += points;
67 if (c_active_windings > c_peak_windings)
68 c_peak_windings = c_active_windings;
70 s = sizeof(vec_t)*3*points + sizeof(int);
76 void FreeWinding (winding_t *w)
78 if (*(unsigned *)w == 0xdeaddead)
79 Error ("FreeWinding: freed a freed winding");
80 *(unsigned *)w = 0xdeaddead;
94 void RemoveColinearPoints (winding_t *w)
99 vec3_t p[MAX_POINTS_ON_WINDING];
102 for (i=0 ; i<w->numpoints ; i++)
104 j = (i+1)%w->numpoints;
105 k = (i+w->numpoints-1)%w->numpoints;
106 VectorSubtract (w->p[j], w->p[i], v1);
107 VectorSubtract (w->p[i], w->p[k], v2);
108 VectorNormalize(v1,v1);
109 VectorNormalize(v2,v2);
110 if (DotProduct(v1, v2) < 0.999)
112 VectorCopy (w->p[i], p[nump]);
117 if (nump == w->numpoints)
121 c_removed += w->numpoints - nump;
123 memcpy (w->p, p, nump*sizeof(p[0]));
131 void WindingPlane (winding_t *w, vec3_t normal, vec_t *dist)
135 VectorSubtract (w->p[1], w->p[0], v1);
136 VectorSubtract (w->p[2], w->p[0], v2);
137 CrossProduct (v2, v1, normal);
138 VectorNormalize (normal, normal);
139 *dist = DotProduct (w->p[0], normal);
148 vec_t WindingArea (winding_t *w)
151 vec3_t d1, d2, cross;
155 for (i=2 ; i<w->numpoints ; i++)
157 VectorSubtract (w->p[i-1], w->p[0], d1);
158 VectorSubtract (w->p[i], w->p[0], d2);
159 CrossProduct (d1, d2, cross);
160 total += 0.5 * VectorLength ( cross );
165 void WindingBounds (winding_t *w, vec3_t mins, vec3_t maxs)
170 mins[0] = mins[1] = mins[2] = 99999;
171 maxs[0] = maxs[1] = maxs[2] = -99999;
173 for (i=0 ; i<w->numpoints ; i++)
175 for (j=0 ; j<3 ; j++)
191 void WindingCenter (winding_t *w, vec3_t center)
196 VectorCopy (vec3_origin, center);
197 for (i=0 ; i<w->numpoints ; i++)
198 VectorAdd (w->p[i], center, center);
200 scale = 1.0/w->numpoints;
201 VectorScale (center, scale, center);
209 winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
211 // The goal in this function is to replicate the exact behavior that was in the original
212 // BaseWindingForPlane() function (see below). The only thing we're going to change is the
213 // accuracy of the operation. The original code gave a preference for the vup vector to start
214 // out as (0, 0, 1), unless the normal had a dominant Z value, in which case vup started out
215 // as (1, 0, 0). After that, vup was "bent" [along the plane defined by normal and vup] to
216 // become perpendicular to normal. After that the vright vector was computed as the cross
217 // product of vup and normal.
219 // Once these vectors are calculated, I'm constructing the winding points in exactly the same
220 // way as was done in the original function. Orientation is the same.
222 // EDIT: We're also changing the size of the winding polygon; this is a side effect of
223 // eliminating a VectorNormalize() call. The new winding polygon is actually bigger.
227 vec3_t vright, vup, org;
232 for (i = 0; i < 3; i++) {
239 if (x == -1) Error("BaseWindingForPlane: no axis found");
242 case 0: // Fall through to next case.
244 vright[0] = -normal[1];
245 vright[1] = normal[0];
250 vright[1] = -normal[2];
251 vright[2] = normal[1];
254 // vright and normal are now perpendicular, you can prove this by taking their
255 // dot product and seeing that it's always exactly 0 (with no error).
257 // NOTE: vright is NOT a unit vector at this point. vright will have length
258 // not exceeding 1.0. The minimum length that vright can achieve happens when,
259 // for example, the Z and X components of the normal input vector are equal,
260 // and when its Y component is zero. In that case Z and X of the normal vector
261 // are both approximately 0.70711. The resulting vright vector in this case
262 // will have a length of 0.70711.
264 // We're relying on the fact that MAX_WORLD_COORD is a power of 2 to keep
265 // our calculation precise and relatively free of floating point error.
266 // The code will work if that's not the case, but not as well.
267 VectorScale(vright, MAX_WORLD_COORD * 4, vright);
269 // At time time of this writing, MAX_WORLD_COORD was 65536 (2^16). Therefore
270 // the length of vright at this point is at least 185364. A corner of the world
271 // at location (65536, 65536, 65536) is distance 113512 away from the origin.
273 CrossProduct(normal, vright, vup);
275 // vup now has length equal to that of vright.
277 VectorScale(normal, dist, org);
279 // org is now a point on the plane defined by normal and dist. Furthermore,
280 // org, vright, and vup are pairwise perpendicular. Now, the 4 vectors
281 // (+-)vright + (+-)vup have length that is at least sqrt(185364^2 + 185364^2),
282 // which is about 262144. That length lies outside the world, since the furthest
283 // point in the world has distance 113512 from the origin as mentioned above.
284 // Also, these 4 vectors are perpendicular to the org vector. So adding them
285 // to org will only increase their length. Therefore the 4 points defined below
286 // all lie outside of the world. Furthermore, it can be easily seen that the
287 // edges connecting these 4 points (in the winding_t below) lie completely outside
288 // the world. sqrt(262144^2 + 262144^2)/2 = 185363, which is greater than 113512.
292 VectorSubtract(org, vright, w->p[0]);
293 VectorAdd(w->p[0], vup, w->p[0]);
295 VectorAdd(org, vright, w->p[1]);
296 VectorAdd(w->p[1], vup, w->p[1]);
298 VectorAdd(org, vright, w->p[2]);
299 VectorSubtract(w->p[2], vup, w->p[2]);
301 VectorSubtract(org, vright, w->p[3]);
302 VectorSubtract(w->p[3], vup, w->p[3]);
309 // Old function, not used but here for reference. Please do not modify it.
310 // (You may remove it at some point.)
311 winding_t *_BaseWindingForPlane_orig_(vec3_t normal, vec_t dist)
315 vec3_t org, vright, vup;
318 // find the major axis
332 Error ("BaseWindingForPlane: no axis found");
334 VectorCopy (vec3_origin, vup);
346 v = DotProduct (vup, normal);
347 VectorMA (vup, -v, normal, vup);
348 VectorNormalize (vup, vup);
350 VectorScale (normal, dist, org);
352 CrossProduct (vup, normal, vright);
354 // LordHavoc: this has to use *2 because otherwise some created points may
355 // be inside the world (think of a diagonal case), and any brush with such
356 // points should be removed, failure to detect such cases is disasterous
357 VectorScale (vup, MAX_WORLD_COORD*2, vup);
358 VectorScale (vright, MAX_WORLD_COORD*2, vright);
360 // project a really big axis aligned box onto the plane
361 w = AllocWinding (4);
363 VectorSubtract (org, vright, w->p[0]);
364 VectorAdd (w->p[0], vup, w->p[0]);
366 VectorAdd (org, vright, w->p[1]);
367 VectorAdd (w->p[1], vup, w->p[1]);
369 VectorAdd (org, vright, w->p[2]);
370 VectorSubtract (w->p[2], vup, w->p[2]);
372 VectorSubtract (org, vright, w->p[3]);
373 VectorSubtract (w->p[3], vup, w->p[3]);
385 winding_t *CopyWinding (winding_t *w)
390 c = AllocWinding (w->numpoints);
391 size = (int)((size_t)((winding_t *)0)->p[w->numpoints]);
401 winding_t *ReverseWinding (winding_t *w)
406 c = AllocWinding (w->numpoints);
407 for (i=0 ; i<w->numpoints ; i++)
409 VectorCopy (w->p[w->numpoints-1-i], c->p[i]);
411 c->numpoints = w->numpoints;
421 void ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist,
422 vec_t epsilon, winding_t **front, winding_t **back)
424 vec_t dists[MAX_POINTS_ON_WINDING+4];
425 int sides[MAX_POINTS_ON_WINDING+4];
427 static vec_t dot; // VC 4.2 optimizer bug if not static
434 counts[0] = counts[1] = counts[2] = 0;
436 // determine sides for each point
437 for (i=0 ; i<in->numpoints ; i++)
440 dot = DotProduct (in->p[i], normal);
444 sides[i] = SIDE_FRONT;
445 else if (dot < -epsilon)
446 sides[i] = SIDE_BACK;
456 *front = *back = NULL;
460 *back = CopyWinding (in);
465 *front = CopyWinding (in);
469 maxpts = in->numpoints+4; // cant use counts[0]+2 because
470 // of fp grouping errors
472 *front = f = AllocWinding (maxpts);
473 *back = b = AllocWinding (maxpts);
475 for (i=0 ; i<in->numpoints ; i++)
479 if (sides[i] == SIDE_ON)
481 VectorCopy (p1, f->p[f->numpoints]);
483 VectorCopy (p1, b->p[b->numpoints]);
488 if (sides[i] == SIDE_FRONT)
490 VectorCopy (p1, f->p[f->numpoints]);
493 if (sides[i] == SIDE_BACK)
495 VectorCopy (p1, b->p[b->numpoints]);
499 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
502 // generate a split point
503 p2 = in->p[(i+1)%in->numpoints];
505 dot = dists[i] / (dists[i]-dists[i+1]);
506 for (j=0 ; j<3 ; j++)
507 { // avoid round off error when possible
510 else if (normal[j] == -1)
513 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
516 VectorCopy (mid, f->p[f->numpoints]);
518 VectorCopy (mid, b->p[b->numpoints]);
522 if (f->numpoints > maxpts || b->numpoints > maxpts)
523 Error ("ClipWinding: points exceeded estimate");
524 if (f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING)
525 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
534 void ChopWindingInPlace (winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon)
537 vec_t dists[MAX_POINTS_ON_WINDING+4];
538 int sides[MAX_POINTS_ON_WINDING+4];
540 static vec_t dot; // VC 4.2 optimizer bug if not static
548 counts[0] = counts[1] = counts[2] = 0;
550 // determine sides for each point
551 for (i=0 ; i<in->numpoints ; i++)
553 dot = DotProduct (in->p[i], normal);
557 sides[i] = SIDE_FRONT;
558 else if (dot < -epsilon)
559 sides[i] = SIDE_BACK;
576 return; // inout stays the same
578 maxpts = in->numpoints+4; // cant use counts[0]+2 because
579 // of fp grouping errors
581 f = AllocWinding (maxpts);
583 for (i=0 ; i<in->numpoints ; i++)
587 if (sides[i] == SIDE_ON)
589 VectorCopy (p1, f->p[f->numpoints]);
594 if (sides[i] == SIDE_FRONT)
596 VectorCopy (p1, f->p[f->numpoints]);
600 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
603 // generate a split point
604 p2 = in->p[(i+1)%in->numpoints];
606 dot = dists[i] / (dists[i]-dists[i+1]);
607 for (j=0 ; j<3 ; j++)
608 { // avoid round off error when possible
611 else if (normal[j] == -1)
614 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
617 VectorCopy (mid, f->p[f->numpoints]);
621 if (f->numpoints > maxpts)
622 Error ("ClipWinding: points exceeded estimate");
623 if (f->numpoints > MAX_POINTS_ON_WINDING)
624 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
635 Returns the fragment of in that is on the front side
636 of the cliping plane. The original is freed.
639 winding_t *ChopWinding (winding_t *in, vec3_t normal, vec_t dist)
643 ClipWindingEpsilon (in, normal, dist, ON_EPSILON, &f, &b);
657 void CheckWinding (winding_t *w)
662 vec3_t dir, edgenormal, facenormal;
666 if (w->numpoints < 3)
667 Error ("CheckWinding: %i points",w->numpoints);
669 area = WindingArea(w);
671 Error ("CheckWinding: %f area", area);
673 WindingPlane (w, facenormal, &facedist);
675 for (i=0 ; i<w->numpoints ; i++)
679 for (j=0 ; j<3 ; j++)
680 if (p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD)
681 Error ("CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j]);
683 j = i+1 == w->numpoints ? 0 : i+1;
685 // check the point is on the face plane
686 d = DotProduct (p1, facenormal) - facedist;
687 if (d < -ON_EPSILON || d > ON_EPSILON)
688 Error ("CheckWinding: point off plane");
690 // check the edge isnt degenerate
692 VectorSubtract (p2, p1, dir);
694 if (VectorLength (dir) < ON_EPSILON)
695 Error ("CheckWinding: degenerate edge");
697 CrossProduct (facenormal, dir, edgenormal);
698 VectorNormalize (edgenormal, edgenormal);
699 edgedist = DotProduct (p1, edgenormal);
700 edgedist += ON_EPSILON;
702 // all other points must be on front side
703 for (j=0 ; j<w->numpoints ; j++)
707 d = DotProduct (w->p[j], edgenormal);
709 Error ("CheckWinding: non-convex");
720 int WindingOnPlaneSide (winding_t *w, vec3_t normal, vec_t dist)
722 qboolean front, back;
728 for (i=0 ; i<w->numpoints ; i++)
730 d = DotProduct (w->p[i], normal) - dist;
757 AddWindingToConvexHull
759 Both w and *hull are on the same plane
762 #define MAX_HULL_POINTS 128
763 void AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
768 int numHullPoints, numNew;
769 vec3_t hullPoints[MAX_HULL_POINTS];
770 vec3_t newHullPoints[MAX_HULL_POINTS];
771 vec3_t hullDirs[MAX_HULL_POINTS];
772 qboolean hullSide[MAX_HULL_POINTS];
776 *hull = CopyWinding( w );
780 numHullPoints = (*hull)->numpoints;
781 memcpy( hullPoints, (*hull)->p, numHullPoints * sizeof(vec3_t) );
783 for ( i = 0 ; i < w->numpoints ; i++ ) {
786 // calculate hull side vectors
787 for ( j = 0 ; j < numHullPoints ; j++ ) {
788 k = ( j + 1 ) % numHullPoints;
790 VectorSubtract( hullPoints[k], hullPoints[j], dir );
791 VectorNormalize( dir, dir );
792 CrossProduct( normal, dir, hullDirs[j] );
796 for ( j = 0 ; j < numHullPoints ; j++ ) {
797 VectorSubtract( p, hullPoints[j], dir );
798 d = DotProduct( dir, hullDirs[j] );
799 if ( d >= ON_EPSILON ) {
802 if ( d >= -ON_EPSILON ) {
805 hullSide[j] = qfalse;
809 // if the point is effectively inside, do nothing
814 // find the back side to front side transition
815 for ( j = 0 ; j < numHullPoints ; j++ ) {
816 if ( !hullSide[ j % numHullPoints ] && hullSide[ (j + 1) % numHullPoints ] ) {
820 if ( j == numHullPoints ) {
824 // insert the point here
825 VectorCopy( p, newHullPoints[0] );
828 // copy over all points that aren't double fronts
829 j = (j+1)%numHullPoints;
830 for ( k = 0 ; k < numHullPoints ; k++ ) {
831 if ( hullSide[ (j+k) % numHullPoints ] && hullSide[ (j+k+1) % numHullPoints ] ) {
834 copy = hullPoints[ (j+k+1) % numHullPoints ];
835 VectorCopy( copy, newHullPoints[numNew] );
839 numHullPoints = numNew;
840 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof(vec3_t) );
843 FreeWinding( *hull );
844 w = AllocWinding( numHullPoints );
845 w->numpoints = numHullPoints;
847 memcpy( w->p, hullPoints, numHullPoints * sizeof(vec3_t) );