svn r377 by Rambetter:
[divverent/netradiant.git] / tools / quake3 / common / polylib.c
1 /*
2 Copyright (C) 1999-2006 Id Software, Inc. and contributors.
3 For a list of contributors, see the accompanying CONTRIBUTORS file.
4
5 This file is part of GtkRadiant.
6
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.
11
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.
16
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
20 */
21
22
23 #include "cmdlib.h"
24 #include "mathlib.h"
25 #include "inout.h"
26 #include "polylib.h"
27 #include "qfiles.h"
28
29
30 extern int numthreads;
31
32 // counters are only bumped when running single threaded,
33 // because they are an awefull coherence problem
34 int     c_active_windings;
35 int     c_peak_windings;
36 int     c_winding_allocs;
37 int     c_winding_points;
38
39 #define BOGUS_RANGE     WORLD_SIZE
40
41 void pw(winding_t *w)
42 {
43         int             i;
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]);
46 }
47
48
49 /*
50 =============
51 AllocWinding
52 =============
53 */
54 winding_t       *AllocWinding (int points)
55 {
56         winding_t       *w;
57         int                     s;
58
59   if (points >= MAX_POINTS_ON_WINDING)
60     Error ("AllocWinding failed: MAX_POINTS_ON_WINDING exceeded");
61
62         if (numthreads == 1)
63         {
64                 c_winding_allocs++;
65                 c_winding_points += points;
66                 c_active_windings++;
67                 if (c_active_windings > c_peak_windings)
68                         c_peak_windings = c_active_windings;
69         }
70         s = sizeof(vec_t)*3*points + sizeof(int);
71         w = safe_malloc (s);
72         memset (w, 0, s); 
73         return w;
74 }
75
76 void FreeWinding (winding_t *w)
77 {
78         if (*(unsigned *)w == 0xdeaddead)
79                 Error ("FreeWinding: freed a freed winding");
80         *(unsigned *)w = 0xdeaddead;
81
82         if (numthreads == 1)
83                 c_active_windings--;
84         free (w);
85 }
86
87 /*
88 ============
89 RemoveColinearPoints
90 ============
91 */
92 int     c_removed;
93
94 void    RemoveColinearPoints (winding_t *w)
95 {
96         int             i, j, k;
97         vec3_t  v1, v2;
98         int             nump;
99         vec3_t  p[MAX_POINTS_ON_WINDING];
100
101         nump = 0;
102         for (i=0 ; i<w->numpoints ; i++)
103         {
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)
111                 {
112                         VectorCopy (w->p[i], p[nump]);
113                         nump++;
114                 }
115         }
116
117         if (nump == w->numpoints)
118                 return;
119
120         if (numthreads == 1)
121                 c_removed += w->numpoints - nump;
122         w->numpoints = nump;
123         memcpy (w->p, p, nump*sizeof(p[0]));
124 }
125
126 /*
127 ============
128 WindingPlane
129 ============
130 */
131 void WindingPlane (winding_t *w, vec3_t normal, vec_t *dist)
132 {
133         vec3_t  v1, v2;
134
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);
140
141 }
142
143 /*
144 =============
145 WindingArea
146 =============
147 */
148 vec_t   WindingArea (winding_t *w)
149 {
150         int             i;
151         vec3_t  d1, d2, cross;
152         vec_t   total;
153
154         total = 0;
155         for (i=2 ; i<w->numpoints ; i++)
156         {
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 );
161         }
162         return total;
163 }
164
165 void    WindingBounds (winding_t *w, vec3_t mins, vec3_t maxs)
166 {
167         vec_t   v;
168         int             i,j;
169
170         mins[0] = mins[1] = mins[2] = 99999;
171         maxs[0] = maxs[1] = maxs[2] = -99999;
172
173         for (i=0 ; i<w->numpoints ; i++)
174         {
175                 for (j=0 ; j<3 ; j++)
176                 {
177                         v = w->p[i][j];
178                         if (v < mins[j])
179                                 mins[j] = v;
180                         if (v > maxs[j])
181                                 maxs[j] = v;
182                 }
183         }
184 }
185
186 /*
187 =============
188 WindingCenter
189 =============
190 */
191 void    WindingCenter (winding_t *w, vec3_t center)
192 {
193         int             i;
194         float   scale;
195
196         VectorCopy (vec3_origin, center);
197         for (i=0 ; i<w->numpoints ; i++)
198                 VectorAdd (w->p[i], center, center);
199
200         scale = 1.0/w->numpoints;
201         VectorScale (center, scale, center);
202 }
203
204 /*
205 =================
206 BaseWindingForPlane
207 =================
208 */
209 winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
210 {
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.
218
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.
221
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.
224
225         int             x, i;
226         vec_t           max, v;
227         vec3_t          vright, vup, org;
228         winding_t       *w;
229
230         max = -BOGUS_RANGE;
231         x = -1;
232         for (i = 0; i < 3; i++) {
233                 v = fabs(normal[i]);
234                 if (v > max) {
235                         x = i;
236                         max = v;
237                 }
238         }
239         if (x == -1) Error("BaseWindingForPlane: no axis found");
240
241         switch (x) {
242                 case 0: // Fall through to next case.
243                 case 1:
244                         vright[0] = -normal[1];
245                         vright[1] = normal[0];
246                         vright[2] = 0;
247                         break;
248                 case 2:
249                         vright[0] = 0;
250                         vright[1] = -normal[2];
251                         vright[2] = normal[1];
252                         break;
253         }
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).
256
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.
263
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);
268
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.
272
273         CrossProduct(normal, vright, vup);
274
275         // vup now has length equal to that of vright.
276
277         VectorScale(normal, dist, org);
278
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.
289
290         w = AllocWinding(4);
291
292         VectorSubtract(org, vright, w->p[0]);
293         VectorAdd(w->p[0], vup, w->p[0]);
294
295         VectorAdd(org, vright, w->p[1]);
296         VectorAdd(w->p[1], vup, w->p[1]);
297
298         VectorAdd(org, vright, w->p[2]);
299         VectorSubtract(w->p[2], vup, w->p[2]);
300
301         VectorSubtract(org, vright, w->p[3]);
302         VectorSubtract(w->p[3], vup, w->p[3]);
303
304         w->numpoints = 4;
305
306         return w;
307 }
308
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)
312 {
313         int             i, x;
314         vec_t   max, v;
315         vec3_t  org, vright, vup;
316         winding_t       *w;
317         
318 // find the major axis
319
320         max = -BOGUS_RANGE;
321         x = -1;
322         for (i=0 ; i<3; i++)
323         {
324                 v = fabs(normal[i]);
325                 if (v > max)
326                 {
327                         x = i;
328                         max = v;
329                 }
330         }
331         if (x==-1)
332                 Error ("BaseWindingForPlane: no axis found");
333                 
334         VectorCopy (vec3_origin, vup);  
335         switch (x)
336         {
337         case 0:
338         case 1:
339                 vup[2] = 1;
340                 break;          
341         case 2:
342                 vup[0] = 1;
343                 break;          
344         }
345
346         v = DotProduct (vup, normal);
347         VectorMA (vup, -v, normal, vup);
348         VectorNormalize (vup, vup);
349                 
350         VectorScale (normal, dist, org);
351         
352         CrossProduct (vup, normal, vright);
353         
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);
359
360   // project a really big       axis aligned box onto the plane
361         w = AllocWinding (4);
362         
363         VectorSubtract (org, vright, w->p[0]);
364         VectorAdd (w->p[0], vup, w->p[0]);
365         
366         VectorAdd (org, vright, w->p[1]);
367         VectorAdd (w->p[1], vup, w->p[1]);
368         
369         VectorAdd (org, vright, w->p[2]);
370         VectorSubtract (w->p[2], vup, w->p[2]);
371         
372         VectorSubtract (org, vright, w->p[3]);
373         VectorSubtract (w->p[3], vup, w->p[3]);
374         
375         w->numpoints = 4;
376         
377         return w;       
378 }
379
380 /*
381 ==================
382 CopyWinding
383 ==================
384 */
385 winding_t       *CopyWinding (winding_t *w)
386 {
387         size_t                  size;
388         winding_t       *c;
389
390         c = AllocWinding (w->numpoints);
391         size = (size_t)((winding_t *)NULL)->p[w->numpoints];
392         memcpy (c, w, size);
393         return c;
394 }
395
396 /*
397 ==================
398 ReverseWinding
399 ==================
400 */
401 winding_t       *ReverseWinding (winding_t *w)
402 {
403         int                     i;
404         winding_t       *c;
405
406         c = AllocWinding (w->numpoints);
407         for (i=0 ; i<w->numpoints ; i++)
408         {
409                 VectorCopy (w->p[w->numpoints-1-i], c->p[i]);
410         }
411         c->numpoints = w->numpoints;
412         return c;
413 }
414
415
416 /*
417 =============
418 ClipWindingEpsilon
419 =============
420 */
421 void    ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist, 
422                                 vec_t epsilon, winding_t **front, winding_t **back)
423 {
424         vec_t   dists[MAX_POINTS_ON_WINDING+4];
425         int             sides[MAX_POINTS_ON_WINDING+4];
426         int             counts[3];
427         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
428         int             i, j;
429         vec_t   *p1, *p2;
430         vec3_t  mid;
431         winding_t       *f, *b;
432         int             maxpts;
433         
434         counts[0] = counts[1] = counts[2] = 0;
435
436 // determine sides for each point
437         for (i=0 ; i<in->numpoints ; i++)
438   {
439
440                 dot = DotProduct (in->p[i], normal);
441                 dot -= dist;
442                 dists[i] = dot;
443                 if (dot > epsilon)
444                         sides[i] = SIDE_FRONT;
445                 else if (dot < -epsilon)
446                         sides[i] = SIDE_BACK;
447                 else
448                 {
449                         sides[i] = SIDE_ON;
450                 }
451                 counts[sides[i]]++;
452         }
453         sides[i] = sides[0];
454         dists[i] = dists[0];
455         
456         *front = *back = NULL;
457
458         if (!counts[0])
459         {
460                 *back = CopyWinding (in);
461                 return;
462         }
463         if (!counts[1])
464         {
465                 *front = CopyWinding (in);
466                 return;
467         }
468
469         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
470                                                                 // of fp grouping errors
471
472         *front = f = AllocWinding (maxpts);
473         *back = b = AllocWinding (maxpts);
474                 
475         for (i=0 ; i<in->numpoints ; i++)
476         {
477                 p1 = in->p[i];
478                 
479                 if (sides[i] == SIDE_ON)
480                 {
481                         VectorCopy (p1, f->p[f->numpoints]);
482                         f->numpoints++;
483                         VectorCopy (p1, b->p[b->numpoints]);
484                         b->numpoints++;
485                         continue;
486                 }
487         
488                 if (sides[i] == SIDE_FRONT)
489                 {
490                         VectorCopy (p1, f->p[f->numpoints]);
491                         f->numpoints++;
492                 }
493                 if (sides[i] == SIDE_BACK)
494                 {
495                         VectorCopy (p1, b->p[b->numpoints]);
496                         b->numpoints++;
497                 }
498
499                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
500                         continue;
501                         
502         // generate a split point
503                 p2 = in->p[(i+1)%in->numpoints];
504                 
505                 dot = dists[i] / (dists[i]-dists[i+1]);
506                 for (j=0 ; j<3 ; j++)
507                 {       // avoid round off error when possible
508                         if (normal[j] == 1)
509                                 mid[j] = dist;
510                         else if (normal[j] == -1)
511                                 mid[j] = -dist;
512                         else
513                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
514                 }
515                         
516                 VectorCopy (mid, f->p[f->numpoints]);
517                 f->numpoints++;
518                 VectorCopy (mid, b->p[b->numpoints]);
519                 b->numpoints++;
520         }
521         
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");
526 }
527
528
529 /*
530 =============
531 ChopWindingInPlace
532 =============
533 */
534 void ChopWindingInPlace (winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon)
535 {
536         winding_t       *in;
537         vec_t   dists[MAX_POINTS_ON_WINDING+4];
538         int             sides[MAX_POINTS_ON_WINDING+4];
539         int             counts[3];
540         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
541         int             i, j;
542         vec_t   *p1, *p2;
543         vec3_t  mid;
544         winding_t       *f;
545         int             maxpts;
546
547         in = *inout;
548         counts[0] = counts[1] = counts[2] = 0;
549
550 // determine sides for each point
551         for (i=0 ; i<in->numpoints ; i++)
552         {
553                 dot = DotProduct (in->p[i], normal);
554                 dot -= dist;
555                 dists[i] = dot;
556                 if (dot > epsilon)
557                         sides[i] = SIDE_FRONT;
558                 else if (dot < -epsilon)
559                         sides[i] = SIDE_BACK;
560                 else
561                 {
562                         sides[i] = SIDE_ON;
563                 }
564                 counts[sides[i]]++;
565         }
566         sides[i] = sides[0];
567         dists[i] = dists[0];
568         
569         if (!counts[0])
570         {
571                 FreeWinding (in);
572                 *inout = NULL;
573                 return;
574         }
575         if (!counts[1])
576                 return;         // inout stays the same
577
578         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
579                                                                 // of fp grouping errors
580
581         f = AllocWinding (maxpts);
582                 
583         for (i=0 ; i<in->numpoints ; i++)
584         {
585                 p1 = in->p[i];
586                 
587                 if (sides[i] == SIDE_ON)
588                 {
589                         VectorCopy (p1, f->p[f->numpoints]);
590                         f->numpoints++;
591                         continue;
592                 }
593         
594                 if (sides[i] == SIDE_FRONT)
595                 {
596                         VectorCopy (p1, f->p[f->numpoints]);
597                         f->numpoints++;
598                 }
599
600                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
601                         continue;
602                         
603         // generate a split point
604                 p2 = in->p[(i+1)%in->numpoints];
605                 
606                 dot = dists[i] / (dists[i]-dists[i+1]);
607                 for (j=0 ; j<3 ; j++)
608                 {       // avoid round off error when possible
609                         if (normal[j] == 1)
610                                 mid[j] = dist;
611                         else if (normal[j] == -1)
612                                 mid[j] = -dist;
613                         else
614                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
615                 }
616                         
617                 VectorCopy (mid, f->p[f->numpoints]);
618                 f->numpoints++;
619         }
620         
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");
625
626         FreeWinding (in);
627         *inout = f;
628 }
629
630
631 /*
632 =================
633 ChopWinding
634
635 Returns the fragment of in that is on the front side
636 of the cliping plane.  The original is freed.
637 =================
638 */
639 winding_t       *ChopWinding (winding_t *in, vec3_t normal, vec_t dist)
640 {
641         winding_t       *f, *b;
642
643         ClipWindingEpsilon (in, normal, dist, ON_EPSILON, &f, &b);
644         FreeWinding (in);
645         if (b)
646                 FreeWinding (b);
647         return f;
648 }
649
650
651 /*
652 =================
653 CheckWinding
654
655 =================
656 */
657 void CheckWinding (winding_t *w)
658 {
659         int             i, j;
660         vec_t   *p1, *p2;
661         vec_t   d, edgedist;
662         vec3_t  dir, edgenormal, facenormal;
663         vec_t   area;
664         vec_t   facedist;
665
666         if (w->numpoints < 3)
667                 Error ("CheckWinding: %i points",w->numpoints);
668         
669         area = WindingArea(w);
670         if (area < 1)
671                 Error ("CheckWinding: %f area", area);
672
673         WindingPlane (w, facenormal, &facedist);
674         
675         for (i=0 ; i<w->numpoints ; i++)
676         {
677                 p1 = w->p[i];
678
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]);
682
683                 j = i+1 == w->numpoints ? 0 : i+1;
684                 
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");
689         
690         // check the edge isnt degenerate
691                 p2 = w->p[j];
692                 VectorSubtract (p2, p1, dir);
693                 
694                 if (VectorLength (dir) < ON_EPSILON)
695                         Error ("CheckWinding: degenerate edge");
696                         
697                 CrossProduct (facenormal, dir, edgenormal);
698                 VectorNormalize (edgenormal, edgenormal);
699                 edgedist = DotProduct (p1, edgenormal);
700                 edgedist += ON_EPSILON;
701                 
702         // all other points must be on front side
703                 for (j=0 ; j<w->numpoints ; j++)
704                 {
705                         if (j == i)
706                                 continue;
707                         d = DotProduct (w->p[j], edgenormal);
708                         if (d > edgedist)
709                                 Error ("CheckWinding: non-convex");
710                 }
711         }
712 }
713
714
715 /*
716 ============
717 WindingOnPlaneSide
718 ============
719 */
720 int             WindingOnPlaneSide (winding_t *w, vec3_t normal, vec_t dist)
721 {
722         qboolean        front, back;
723         int                     i;
724         vec_t           d;
725
726         front = qfalse;
727         back = qfalse;
728         for (i=0 ; i<w->numpoints ; i++)
729         {
730                 d = DotProduct (w->p[i], normal) - dist;
731                 if (d < -ON_EPSILON)
732                 {
733                         if (front)
734                                 return SIDE_CROSS;
735                         back = qtrue;
736                         continue;
737                 }
738                 if (d > ON_EPSILON)
739                 {
740                         if (back)
741                                 return SIDE_CROSS;
742                         front = qtrue;
743                         continue;
744                 }
745         }
746
747         if (back)
748                 return SIDE_BACK;
749         if (front)
750                 return SIDE_FRONT;
751         return SIDE_ON;
752 }
753
754
755 /*
756 =================
757 AddWindingToConvexHull
758
759 Both w and *hull are on the same plane
760 =================
761 */
762 #define MAX_HULL_POINTS         128
763 void    AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
764         int                     i, j, k;
765         float           *p, *copy;
766         vec3_t          dir;
767         float           d;
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];
773         qboolean        outside;
774
775         if ( !*hull ) {
776                 *hull = CopyWinding( w );
777                 return;
778         }
779
780         numHullPoints = (*hull)->numpoints;
781         memcpy( hullPoints, (*hull)->p, numHullPoints * sizeof(vec3_t) );
782
783         for ( i = 0 ; i < w->numpoints ; i++ ) {
784                 p = w->p[i];
785
786                 // calculate hull side vectors
787                 for ( j = 0 ; j < numHullPoints ; j++ ) {
788                         k = ( j + 1 ) % numHullPoints;
789
790                         VectorSubtract( hullPoints[k], hullPoints[j], dir );
791                         VectorNormalize( dir, dir );
792                         CrossProduct( normal, dir, hullDirs[j] );
793                 }
794
795                 outside = qfalse;
796                 for ( j = 0 ; j < numHullPoints ; j++ ) {
797                         VectorSubtract( p, hullPoints[j], dir );
798                         d = DotProduct( dir, hullDirs[j] );
799                         if ( d >= ON_EPSILON ) {
800                                 outside = qtrue;
801                         }
802                         if ( d >= -ON_EPSILON ) {
803                                 hullSide[j] = qtrue;
804                         } else {
805                                 hullSide[j] = qfalse;
806                         }
807                 }
808
809                 // if the point is effectively inside, do nothing
810                 if ( !outside ) {
811                         continue;
812                 }
813
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 ] ) {
817                                 break;
818                         }
819                 }
820                 if ( j == numHullPoints ) {
821                         continue;
822                 }
823
824                 // insert the point here
825                 VectorCopy( p, newHullPoints[0] );
826                 numNew = 1;
827
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 ] ) {
832                                 continue;
833                         }
834                         copy = hullPoints[ (j+k+1) % numHullPoints ];
835                         VectorCopy( copy, newHullPoints[numNew] );
836                         numNew++;
837                 }
838
839                 numHullPoints = numNew;
840                 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof(vec3_t) );
841         }
842
843         FreeWinding( *hull );
844         w = AllocWinding( numHullPoints );
845         w->numpoints = numHullPoints;
846         *hull = w;
847         memcpy( w->p, hullPoints, numHullPoints * sizeof(vec3_t) );
848 }
849
850