2 Copyright (C) 1999-2006 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);
81 winding_accu_t *AllocWindingAccu(int points)
86 if (points >= MAX_POINTS_ON_WINDING)
87 Error("AllocWindingAccu failed: MAX_POINTS_ON_WINDING exceeded");
91 // At the time of this writing, these statistics were not used in any way.
93 c_winding_points += points;
95 if (c_active_windings > c_peak_windings)
96 c_peak_windings = c_active_windings;
98 s = sizeof(vec_accu_t) * 3 * points + sizeof(int);
109 void FreeWinding (winding_t *w)
111 if (!w) Error("FreeWinding: winding is NULL");
113 if (*(unsigned *)w == 0xdeaddead)
114 Error ("FreeWinding: freed a freed winding");
115 *(unsigned *)w = 0xdeaddead;
127 void FreeWindingAccu(winding_accu_t *w)
129 if (!w) Error("FreeWindingAccu: winding is NULL");
131 if (*((unsigned *) w) == 0xdeaddead)
132 Error("FreeWindingAccu: freed a freed winding");
133 *((unsigned *) w) = 0xdeaddead;
147 void RemoveColinearPoints (winding_t *w)
152 vec3_t p[MAX_POINTS_ON_WINDING];
155 for (i=0 ; i<w->numpoints ; i++)
157 j = (i+1)%w->numpoints;
158 k = (i+w->numpoints-1)%w->numpoints;
159 VectorSubtract (w->p[j], w->p[i], v1);
160 VectorSubtract (w->p[i], w->p[k], v2);
161 VectorNormalize(v1,v1);
162 VectorNormalize(v2,v2);
163 if (DotProduct(v1, v2) < 0.999)
165 VectorCopy (w->p[i], p[nump]);
170 if (nump == w->numpoints)
174 c_removed += w->numpoints - nump;
176 memcpy (w->p, p, nump*sizeof(p[0]));
184 void WindingPlane (winding_t *w, vec3_t normal, vec_t *dist)
188 VectorSubtract (w->p[1], w->p[0], v1);
189 VectorSubtract (w->p[2], w->p[0], v2);
190 CrossProduct (v2, v1, normal);
191 VectorNormalize (normal, normal);
192 *dist = DotProduct (w->p[0], normal);
201 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)
223 mins[0] = mins[1] = mins[2] = 99999;
224 maxs[0] = maxs[1] = maxs[2] = -99999;
226 for (i=0 ; i<w->numpoints ; i++)
228 for (j=0 ; j<3 ; j++)
244 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)
264 // The goal in this function is to replicate the behavior of the original BaseWindingForPlane()
265 // function (see below) but at the same time increasing accuracy substantially.
267 // The original code gave a preference for the vup vector to start out as (0, 0, 1), unless the
268 // normal had a dominant Z value, in which case vup started out as (1, 0, 0). After that, vup
269 // was "bent" [along the plane defined by normal and vup] to become perpendicular to normal.
270 // After that the vright vector was computed as the cross product of vup and normal.
272 // I'm constructing the winding polygon points in a fashion similar to the method used in the
273 // original function. Orientation is the same. The size of the winding polygon, however, is
274 // variable in this function (depending on the angle of normal), and is larger (by about a factor
275 // of 2) than the winding polygon in the original function.
279 vec3_accu_t vright, vup, org, normalAccu;
282 // One of the components of normal must have a magnitiude greater than this value,
283 // otherwise normal is not a unit vector. This is a little bit of inexpensive
284 // partial error checking we can do.
285 max = 0.56; // 1 / sqrt(1^2 + 1^2 + 1^2) = 0.577350269
288 for (i = 0; i < 3; i++) {
289 v = (vec_t) fabs(normal[i]);
295 if (x == -1) Error("BaseWindingForPlaneAccu: no dominant axis found because normal is too short");
298 case 0: // Fall through to next case.
300 vright[0] = (vec_accu_t) -normal[1];
301 vright[1] = (vec_accu_t) normal[0];
306 vright[1] = (vec_accu_t) -normal[2];
307 vright[2] = (vec_accu_t) normal[1];
311 // vright and normal are now perpendicular; you can prove this by taking their
312 // dot product and seeing that it's always exactly 0 (with no error).
314 // NOTE: vright is NOT a unit vector at this point. vright will have length
315 // not exceeding 1.0. The minimum length that vright can achieve happens when,
316 // for example, the Z and X components of the normal input vector are equal,
317 // and when normal's Y component is zero. In that case Z and X of the normal
318 // vector are both approximately 0.70711. The resulting vright vector in this
319 // case will have a length of 0.70711.
321 // We're relying on the fact that MAX_WORLD_COORD is a power of 2 to keep
322 // our calculation precise and relatively free of floating point error.
323 // [However, the code will still work fine if that's not the case.]
324 VectorScaleAccu(vright, ((vec_accu_t) MAX_WORLD_COORD) * 4.0, vright);
326 // At time time of this writing, MAX_WORLD_COORD was 65536 (2^16). Therefore
327 // the length of vright at this point is at least 185364. In comparison, a
328 // corner of the world at location (65536, 65536, 65536) is distance 113512
329 // away from the origin.
331 VectorCopyRegularToAccu(normal, normalAccu);
332 CrossProductAccu(normalAccu, vright, vup);
334 // vup now has length equal to that of vright.
336 VectorScaleAccu(normalAccu, (vec_accu_t) dist, org);
338 // org is now a point on the plane defined by normal and dist. Furthermore,
339 // org, vright, and vup are pairwise perpendicular. Now, the 4 vectors
340 // { (+-)vright + (+-)vup } have length that is at least sqrt(185364^2 + 185364^2),
341 // which is about 262144. That length lies outside the world, since the furthest
342 // point in the world has distance 113512 from the origin as mentioned above.
343 // Also, these 4 vectors are perpendicular to the org vector. So adding them
344 // to org will only increase their length. Therefore the 4 points defined below
345 // all lie outside of the world. Furthermore, it can be easily seen that the
346 // edges connecting these 4 points (in the winding_accu_t below) lie completely
347 // outside the world. sqrt(262144^2 + 262144^2)/2 = 185363, which is greater than
350 w = AllocWindingAccu(4);
352 VectorSubtractAccu(org, vright, w->p[0]);
353 VectorAddAccu(w->p[0], vup, w->p[0]);
355 VectorAddAccu(org, vright, w->p[1]);
356 VectorAddAccu(w->p[1], vup, w->p[1]);
358 VectorAddAccu(org, vright, w->p[2]);
359 VectorSubtractAccu(w->p[2], vup, w->p[2]);
361 VectorSubtractAccu(org, vright, w->p[3]);
362 VectorSubtractAccu(w->p[3], vup, w->p[3]);
373 Original BaseWindingForPlane() function that has serious accuracy problems. Here is why.
374 The base winding is computed as a rectangle with very large coordinates. These coordinates
375 are in the range 2^17 or 2^18. "Epsilon" (meaning the distance between two adjacent numbers)
376 at these magnitudes in 32 bit floating point world is about 0.02. So for example, if things
377 go badly (by bad luck), then the whole plane could be shifted by 0.02 units (its distance could
378 be off by that much). Then if we were to compute the winding of this plane and another of
379 the brush's planes met this winding at a very acute angle, that error could multiply to around
380 0.1 or more when computing the final vertex coordinates of the winding. 0.1 is a very large
381 error, and can lead to all sorts of disappearing triangle problems.
384 winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
388 vec3_t org, vright, vup;
391 // find the major axis
405 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 if (!w) Error("CopyWinding: winding is NULL");
465 c = AllocWinding (w->numpoints);
466 size = (size_t)((winding_t *)NULL)->p[w->numpoints];
473 CopyWindingAccuIncreaseSizeAndFreeOld
476 winding_accu_t *CopyWindingAccuIncreaseSizeAndFreeOld(winding_accu_t *w)
481 if (!w) Error("CopyWindingAccuIncreaseSizeAndFreeOld: winding is NULL");
483 c = AllocWindingAccu(w->numpoints + 1);
484 c->numpoints = w->numpoints;
485 for (i = 0; i < c->numpoints; i++)
487 VectorCopyAccu(w->p[i], c->p[i]);
495 CopyWindingAccuToRegular
498 winding_t *CopyWindingAccuToRegular(winding_accu_t *w)
503 if (!w) Error("CopyWindingAccuToRegular: winding is NULL");
505 c = AllocWinding(w->numpoints);
506 c->numpoints = w->numpoints;
507 for (i = 0; i < c->numpoints; i++)
509 VectorCopyAccuToRegular(w->p[i], c->p[i]);
519 winding_t *ReverseWinding (winding_t *w)
524 c = AllocWinding (w->numpoints);
525 for (i=0 ; i<w->numpoints ; i++)
527 VectorCopy (w->p[w->numpoints-1-i], c->p[i]);
529 c->numpoints = w->numpoints;
539 void ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist,
540 vec_t epsilon, winding_t **front, winding_t **back)
542 vec_t dists[MAX_POINTS_ON_WINDING+4];
543 int sides[MAX_POINTS_ON_WINDING+4];
545 static vec_t dot; // VC 4.2 optimizer bug if not static
552 counts[0] = counts[1] = counts[2] = 0;
554 // determine sides for each point
555 for (i=0 ; i<in->numpoints ; i++)
558 dot = DotProduct (in->p[i], normal);
562 sides[i] = SIDE_FRONT;
563 else if (dot < -epsilon)
564 sides[i] = SIDE_BACK;
574 *front = *back = NULL;
578 *back = CopyWinding (in);
583 *front = CopyWinding (in);
587 maxpts = in->numpoints+4; // cant use counts[0]+2 because
588 // of fp grouping errors
590 *front = f = AllocWinding (maxpts);
591 *back = b = AllocWinding (maxpts);
593 for (i=0 ; i<in->numpoints ; i++)
597 if (sides[i] == SIDE_ON)
599 VectorCopy (p1, f->p[f->numpoints]);
601 VectorCopy (p1, b->p[b->numpoints]);
606 if (sides[i] == SIDE_FRONT)
608 VectorCopy (p1, f->p[f->numpoints]);
611 if (sides[i] == SIDE_BACK)
613 VectorCopy (p1, b->p[b->numpoints]);
617 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
620 // generate a split point
621 p2 = in->p[(i+1)%in->numpoints];
623 dot = dists[i] / (dists[i]-dists[i+1]);
624 for (j=0 ; j<3 ; j++)
625 { // avoid round off error when possible
628 else if (normal[j] == -1)
631 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
634 VectorCopy (mid, f->p[f->numpoints]);
636 VectorCopy (mid, b->p[b->numpoints]);
640 if (f->numpoints > maxpts || b->numpoints > maxpts)
641 Error ("ClipWinding: points exceeded estimate");
642 if (f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING)
643 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
649 ChopWindingInPlaceAccu
652 void ChopWindingInPlaceAccu(winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon)
654 vec_accu_t fineEpsilon;
658 vec_accu_t dists[MAX_POINTS_ON_WINDING + 1];
659 int sides[MAX_POINTS_ON_WINDING + 1];
664 vec3_accu_t mid, normalAccu;
666 // We require at least a very small epsilon. It's a good idea for several reasons.
667 // First, we will be dividing by a potentially very small distance below. We don't
668 // want that distance to be too small; otherwise, things "blow up" with little accuracy
669 // due to the division. (After a second look, the value w below is in range (0,1), but
670 // graininess problem remains.) Second, Having minimum epsilon also prevents the following
671 // situation. Say for example we have a perfect octagon defined by the input winding.
672 // Say our chopping plane (defined by normal and dist) is essentially the same plane
673 // that the octagon is sitting on. Well, due to rounding errors, it may be that point
674 // 1 of the octagon might be in front, point 2 might be in back, point 3 might be in
675 // front, point 4 might be in back, and so on. So we could end up with a very ugly-
676 // looking chopped winding, and this might be undesirable, and would at least lead to
677 // a possible exhaustion of MAX_POINTS_ON_WINDING. It's better to assume that points
678 // very very close to the plane are on the plane, using an infinitesimal epsilon amount.
680 // Now, the original ChopWindingInPlace() function used a vec_t-based winding_t.
681 // So this minimum epsilon is quite similar to casting the higher resolution numbers to
682 // the lower resolution and comparing them in the lower resolution mode. We explicitly
683 // choose the minimum epsilon as something around the vec_t epsilon of one because we
684 // want the resolution of vec_accu_t to have a large resolution around the epsilon.
685 // Some of that leftover resolution even goes away after we scale to points far away.
687 // Here is a further discussion regarding the choice of smallestEpsilonAllowed.
688 // In the 32 float world (we can assume vec_t is that), the "epsilon around 1.0" is
689 // 0.00000011921. In the 64 bit float world (we can assume vec_accu_t is that), the
690 // "epsilon around 1.0" is 0.00000000000000022204. (By the way these two epsilons
691 // are defined as VEC_SMALLEST_EPSILON_AROUND_ONE VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE
692 // respectively.) If you divide the first by the second, you get approximately
693 // 536,885,246. Dividing that number by 200,000 (a typical base winding coordinate)
694 // gives 2684. So in other words, if our smallestEpsilonAllowed was chosen as exactly
695 // VEC_SMALLEST_EPSILON_AROUND_ONE, you would be guaranteed at least 2000 "ticks" in
696 // 64-bit land inside of the epsilon for all numbers we're dealing with.
698 static const vec_accu_t smallestEpsilonAllowed = ((vec_accu_t) VEC_SMALLEST_EPSILON_AROUND_ONE) * 0.5;
699 if (crudeEpsilon < smallestEpsilonAllowed) fineEpsilon = smallestEpsilonAllowed;
700 else fineEpsilon = (vec_accu_t) crudeEpsilon;
703 counts[0] = counts[1] = counts[2] = 0;
704 VectorCopyRegularToAccu(normal, normalAccu);
706 for (i = 0; i < in->numpoints; i++)
708 dists[i] = DotProductAccu(in->p[i], normalAccu) - dist;
709 if (dists[i] > fineEpsilon) sides[i] = SIDE_FRONT;
710 else if (dists[i] < -fineEpsilon) sides[i] = SIDE_BACK;
711 else sides[i] = SIDE_ON;
717 // I'm wondering if whatever code that handles duplicate planes is robust enough
718 // that we never get a case where two nearly equal planes result in 2 NULL windings
719 // due to the 'if' statement below. TODO: Investigate this.
720 if (!counts[SIDE_FRONT]) {
725 if (!counts[SIDE_BACK]) {
726 return; // Winding is unmodified.
729 // NOTE: The least number of points that a winding can have at this point is 2.
730 // In that case, one point is SIDE_FRONT and the other is SIDE_BACK.
732 maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
733 f = AllocWindingAccu(maxpts);
735 for (i = 0; i < in->numpoints; i++)
739 if (sides[i] == SIDE_ON || sides[i] == SIDE_FRONT)
741 if (f->numpoints >= MAX_POINTS_ON_WINDING)
742 Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
743 if (f->numpoints >= maxpts) // This will probably never happen.
745 Sys_FPrintf(SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n");
746 f = CopyWindingAccuIncreaseSizeAndFreeOld(f);
749 VectorCopyAccu(p1, f->p[f->numpoints]);
751 if (sides[i] == SIDE_ON) continue;
753 if (sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i])
758 // Generate a split point.
759 p2 = in->p[((i + 1) == in->numpoints) ? 0 : (i + 1)];
761 // The divisor's absolute value is greater than the dividend's absolute value.
762 // w is in the range (0,1).
763 w = dists[i] / (dists[i] - dists[i + 1]);
765 for (j = 0; j < 3; j++)
767 // Avoid round-off error when possible. Check axis-aligned normal.
768 if (normal[j] == 1) mid[j] = dist;
769 else if (normal[j] == -1) mid[j] = -dist;
770 else mid[j] = p1[j] + (w * (p2[j] - p1[j]));
772 if (f->numpoints >= MAX_POINTS_ON_WINDING)
773 Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
774 if (f->numpoints >= maxpts) // This will probably never happen.
776 Sys_FPrintf(SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n");
777 f = CopyWindingAccuIncreaseSizeAndFreeOld(f);
780 VectorCopyAccu(mid, f->p[f->numpoints]);
793 void ChopWindingInPlace (winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon)
796 vec_t dists[MAX_POINTS_ON_WINDING+4];
797 int sides[MAX_POINTS_ON_WINDING+4];
799 static vec_t dot; // VC 4.2 optimizer bug if not static
807 counts[0] = counts[1] = counts[2] = 0;
809 // determine sides for each point
810 for (i=0 ; i<in->numpoints ; i++)
812 dot = DotProduct (in->p[i], normal);
816 sides[i] = SIDE_FRONT;
817 else if (dot < -epsilon)
818 sides[i] = SIDE_BACK;
835 return; // inout stays the same
837 maxpts = in->numpoints+4; // cant use counts[0]+2 because
838 // of fp grouping errors
840 f = AllocWinding (maxpts);
842 for (i=0 ; i<in->numpoints ; i++)
846 if (sides[i] == SIDE_ON)
848 VectorCopy (p1, f->p[f->numpoints]);
853 if (sides[i] == SIDE_FRONT)
855 VectorCopy (p1, f->p[f->numpoints]);
859 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
862 // generate a split point
863 p2 = in->p[(i+1)%in->numpoints];
865 dot = dists[i] / (dists[i]-dists[i+1]);
866 for (j=0 ; j<3 ; j++)
867 { // avoid round off error when possible
870 else if (normal[j] == -1)
873 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
876 VectorCopy (mid, f->p[f->numpoints]);
880 if (f->numpoints > maxpts)
881 Error ("ClipWinding: points exceeded estimate");
882 if (f->numpoints > MAX_POINTS_ON_WINDING)
883 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
894 Returns the fragment of in that is on the front side
895 of the cliping plane. The original is freed.
898 winding_t *ChopWinding (winding_t *in, vec3_t normal, vec_t dist)
902 ClipWindingEpsilon (in, normal, dist, ON_EPSILON, &f, &b);
916 void CheckWinding (winding_t *w)
921 vec3_t dir, edgenormal, facenormal;
925 if (w->numpoints < 3)
926 Error ("CheckWinding: %i points",w->numpoints);
928 area = WindingArea(w);
930 Error ("CheckWinding: %f area", area);
932 WindingPlane (w, facenormal, &facedist);
934 for (i=0 ; i<w->numpoints ; i++)
938 for (j=0 ; j<3 ; j++)
939 if (p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD)
940 Error ("CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j]);
942 j = i+1 == w->numpoints ? 0 : i+1;
944 // check the point is on the face plane
945 d = DotProduct (p1, facenormal) - facedist;
946 if (d < -ON_EPSILON || d > ON_EPSILON)
947 Error ("CheckWinding: point off plane");
949 // check the edge isnt degenerate
951 VectorSubtract (p2, p1, dir);
953 if (VectorLength (dir) < ON_EPSILON)
954 Error ("CheckWinding: degenerate edge");
956 CrossProduct (facenormal, dir, edgenormal);
957 VectorNormalize (edgenormal, edgenormal);
958 edgedist = DotProduct (p1, edgenormal);
959 edgedist += ON_EPSILON;
961 // all other points must be on front side
962 for (j=0 ; j<w->numpoints ; j++)
966 d = DotProduct (w->p[j], edgenormal);
968 Error ("CheckWinding: non-convex");
979 int WindingOnPlaneSide (winding_t *w, vec3_t normal, vec_t dist)
981 qboolean front, back;
987 for (i=0 ; i<w->numpoints ; i++)
989 d = DotProduct (w->p[i], normal) - dist;
1016 AddWindingToConvexHull
1018 Both w and *hull are on the same plane
1021 #define MAX_HULL_POINTS 128
1022 void AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
1027 int numHullPoints, numNew;
1028 vec3_t hullPoints[MAX_HULL_POINTS];
1029 vec3_t newHullPoints[MAX_HULL_POINTS];
1030 vec3_t hullDirs[MAX_HULL_POINTS];
1031 qboolean hullSide[MAX_HULL_POINTS];
1035 *hull = CopyWinding( w );
1039 numHullPoints = (*hull)->numpoints;
1040 memcpy( hullPoints, (*hull)->p, numHullPoints * sizeof(vec3_t) );
1042 for ( i = 0 ; i < w->numpoints ; i++ ) {
1045 // calculate hull side vectors
1046 for ( j = 0 ; j < numHullPoints ; j++ ) {
1047 k = ( j + 1 ) % numHullPoints;
1049 VectorSubtract( hullPoints[k], hullPoints[j], dir );
1050 VectorNormalize( dir, dir );
1051 CrossProduct( normal, dir, hullDirs[j] );
1055 for ( j = 0 ; j < numHullPoints ; j++ ) {
1056 VectorSubtract( p, hullPoints[j], dir );
1057 d = DotProduct( dir, hullDirs[j] );
1058 if ( d >= ON_EPSILON ) {
1061 if ( d >= -ON_EPSILON ) {
1062 hullSide[j] = qtrue;
1064 hullSide[j] = qfalse;
1068 // if the point is effectively inside, do nothing
1073 // find the back side to front side transition
1074 for ( j = 0 ; j < numHullPoints ; j++ ) {
1075 if ( !hullSide[ j % numHullPoints ] && hullSide[ (j + 1) % numHullPoints ] ) {
1079 if ( j == numHullPoints ) {
1083 // insert the point here
1084 VectorCopy( p, newHullPoints[0] );
1087 // copy over all points that aren't double fronts
1088 j = (j+1)%numHullPoints;
1089 for ( k = 0 ; k < numHullPoints ; k++ ) {
1090 if ( hullSide[ (j+k) % numHullPoints ] && hullSide[ (j+k+1) % numHullPoints ] ) {
1093 copy = hullPoints[ (j+k+1) % numHullPoints ];
1094 VectorCopy( copy, newHullPoints[numNew] );
1098 numHullPoints = numNew;
1099 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof(vec3_t) );
1102 FreeWinding( *hull );
1103 w = AllocWinding( numHullPoints );
1104 w->numpoints = numHullPoints;
1106 memcpy( w->p, hullPoints, numHullPoints * sizeof(vec3_t) );