NOW I do it right: #woxblox#
[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(*w) + (points ? sizeof(w->p[0])*(points-1) : 0);
71         w = safe_malloc (s);
72         memset (w, 0, s); 
73         return w;
74 }
75
76 /*
77 =============
78 AllocWindingAccu
79 =============
80 */
81 winding_accu_t *AllocWindingAccu(int points)
82 {
83         winding_accu_t  *w;
84         int             s;
85
86         if (points >= MAX_POINTS_ON_WINDING)
87                 Error("AllocWindingAccu failed: MAX_POINTS_ON_WINDING exceeded");
88
89         if (numthreads == 1)
90         {
91                 // At the time of this writing, these statistics were not used in any way.
92                 c_winding_allocs++;
93                 c_winding_points += points;
94                 c_active_windings++;
95                 if (c_active_windings > c_peak_windings)
96                         c_peak_windings = c_active_windings;
97         }
98         s = sizeof(*w) + (points ? sizeof(w->p[0])*(points-1) : 0);
99         w = safe_malloc(s);
100         memset(w, 0, s); 
101         return w;
102 }
103
104 /*
105 =============
106 FreeWinding
107 =============
108 */
109 void FreeWinding (winding_t *w)
110 {
111         if (!w) Error("FreeWinding: winding is NULL");
112
113         if (*(unsigned *)w == 0xdeaddead)
114                 Error ("FreeWinding: freed a freed winding");
115         *(unsigned *)w = 0xdeaddead;
116
117         if (numthreads == 1)
118                 c_active_windings--;
119         free (w);
120 }
121
122 /*
123 =============
124 FreeWindingAccu
125 =============
126 */
127 void FreeWindingAccu(winding_accu_t *w)
128 {
129         if (!w) Error("FreeWindingAccu: winding is NULL");
130
131         if (*((unsigned *) w) == 0xdeaddead)
132                 Error("FreeWindingAccu: freed a freed winding");
133         *((unsigned *) w) = 0xdeaddead;
134
135         if (numthreads == 1)
136                 c_active_windings--;
137         free(w);
138 }
139
140 /*
141 ============
142 RemoveColinearPoints
143 ============
144 */
145 int     c_removed;
146
147 void    RemoveColinearPoints (winding_t *w)
148 {
149         int             i, j, k;
150         vec3_t  v1, v2;
151         int             nump;
152         vec3_t  p[MAX_POINTS_ON_WINDING];
153
154         nump = 0;
155         for (i=0 ; i<w->numpoints ; i++)
156         {
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)
164                 {
165                         VectorCopy (w->p[i], p[nump]);
166                         nump++;
167                 }
168         }
169
170         if (nump == w->numpoints)
171                 return;
172
173         if (numthreads == 1)
174                 c_removed += w->numpoints - nump;
175         w->numpoints = nump;
176         memcpy (w->p, p, nump*sizeof(p[0]));
177 }
178
179 /*
180 ============
181 WindingPlane
182 ============
183 */
184 void WindingPlane (winding_t *w, vec3_t normal, vec_t *dist)
185 {
186         vec3_t  v1, v2;
187
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);
193
194 }
195
196 /*
197 =============
198 WindingArea
199 =============
200 */
201 vec_t   WindingArea (winding_t *w)
202 {
203         int             i;
204         vec3_t  d1, d2, cross;
205         vec_t   total;
206
207         total = 0;
208         for (i=2 ; i<w->numpoints ; i++)
209         {
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 );
214         }
215         return total;
216 }
217
218 void    WindingBounds (winding_t *w, vec3_t mins, vec3_t maxs)
219 {
220         vec_t   v;
221         int             i,j;
222
223         mins[0] = mins[1] = mins[2] = 99999;
224         maxs[0] = maxs[1] = maxs[2] = -99999;
225
226         for (i=0 ; i<w->numpoints ; i++)
227         {
228                 for (j=0 ; j<3 ; j++)
229                 {
230                         v = w->p[i][j];
231                         if (v < mins[j])
232                                 mins[j] = v;
233                         if (v > maxs[j])
234                                 maxs[j] = v;
235                 }
236         }
237 }
238
239 /*
240 =============
241 WindingCenter
242 =============
243 */
244 void    WindingCenter (winding_t *w, vec3_t center)
245 {
246         int             i;
247         float   scale;
248
249         VectorCopy (vec3_origin, center);
250         for (i=0 ; i<w->numpoints ; i++)
251                 VectorAdd (w->p[i], center, center);
252
253         scale = 1.0/w->numpoints;
254         VectorScale (center, scale, center);
255 }
256
257 /*
258 =================
259 BaseWindingForPlaneAccu
260 =================
261 */
262 winding_accu_t *BaseWindingForPlaneAccu(vec3_t normal, vec_t dist)
263 {
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.
266
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.
271
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.
276
277         int             x, i;
278         vec_t           max, v;
279         vec3_accu_t     vright, vup, org, normalAccu;
280         winding_accu_t  *w;
281
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
286
287         x = -1;
288         for (i = 0; i < 3; i++) {
289                 v = (vec_t) fabs(normal[i]);
290                 if (v > max) {
291                         x = i;
292                         max = v;
293                 }
294         }
295         if (x == -1) Error("BaseWindingForPlaneAccu: no dominant axis found because normal is too short");
296
297         switch (x) {
298                 case 0: // Fall through to next case.
299                 case 1:
300                         vright[0] = (vec_accu_t) -normal[1];
301                         vright[1] = (vec_accu_t) normal[0];
302                         vright[2] = 0;
303                         break;
304                 case 2:
305                         vright[0] = 0;
306                         vright[1] = (vec_accu_t) -normal[2];
307                         vright[2] = (vec_accu_t) normal[1];
308                         break;
309         }
310
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).
313
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.
320
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);
325
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.
330
331         VectorCopyRegularToAccu(normal, normalAccu);
332         CrossProductAccu(normalAccu, vright, vup);
333
334         // vup now has length equal to that of vright.
335
336         VectorScaleAccu(normalAccu, (vec_accu_t) dist, org);
337
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
348         // 113512.
349
350         w = AllocWindingAccu(4);
351
352         VectorSubtractAccu(org, vright, w->p[0]);
353         VectorAddAccu(w->p[0], vup, w->p[0]);
354
355         VectorAddAccu(org, vright, w->p[1]);
356         VectorAddAccu(w->p[1], vup, w->p[1]);
357
358         VectorAddAccu(org, vright, w->p[2]);
359         VectorSubtractAccu(w->p[2], vup, w->p[2]);
360
361         VectorSubtractAccu(org, vright, w->p[3]);
362         VectorSubtractAccu(w->p[3], vup, w->p[3]);
363
364         w->numpoints = 4;
365
366         return w;
367 }
368
369 /*
370 =================
371 BaseWindingForPlane
372
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.
382 =================
383 */
384 winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
385 {
386         int             i, x;
387         vec_t   max, v;
388         vec3_t  org, vright, vup;
389         winding_t       *w;
390         
391 // find the major axis
392
393         max = -BOGUS_RANGE;
394         x = -1;
395         for (i=0 ; i<3; i++)
396         {
397                 v = fabs(normal[i]);
398                 if (v > max)
399                 {
400                         x = i;
401                         max = v;
402                 }
403         }
404         if (x==-1)
405                 Error ("BaseWindingForPlane: no axis found");
406                 
407         VectorCopy (vec3_origin, vup);  
408         switch (x)
409         {
410         case 0:
411         case 1:
412                 vup[2] = 1;
413                 break;          
414         case 2:
415                 vup[0] = 1;
416                 break;          
417         }
418
419         v = DotProduct (vup, normal);
420         VectorMA (vup, -v, normal, vup);
421         VectorNormalize (vup, vup);
422                 
423         VectorScale (normal, dist, org);
424         
425         CrossProduct (vup, normal, vright);
426         
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);
432
433   // project a really big       axis aligned box onto the plane
434         w = AllocWinding (4);
435         
436         VectorSubtract (org, vright, w->p[0]);
437         VectorAdd (w->p[0], vup, w->p[0]);
438         
439         VectorAdd (org, vright, w->p[1]);
440         VectorAdd (w->p[1], vup, w->p[1]);
441         
442         VectorAdd (org, vright, w->p[2]);
443         VectorSubtract (w->p[2], vup, w->p[2]);
444         
445         VectorSubtract (org, vright, w->p[3]);
446         VectorSubtract (w->p[3], vup, w->p[3]);
447         
448         w->numpoints = 4;
449         
450         return w;       
451 }
452
453 /*
454 ==================
455 CopyWinding
456 ==================
457 */
458 winding_t       *CopyWinding (winding_t *w)
459 {
460         size_t                  size;
461         winding_t       *c;
462
463         if (!w) Error("CopyWinding: winding is NULL");
464
465         c = AllocWinding (w->numpoints);
466         size = (size_t)((winding_t *)NULL)->p[w->numpoints];
467         memcpy (c, w, size);
468         return c;
469 }
470
471 /*
472 ==================
473 CopyWindingAccuIncreaseSizeAndFreeOld
474 ==================
475 */
476 winding_accu_t *CopyWindingAccuIncreaseSizeAndFreeOld(winding_accu_t *w)
477 {
478         int             i;
479         winding_accu_t  *c;
480
481         if (!w) Error("CopyWindingAccuIncreaseSizeAndFreeOld: winding is NULL");
482
483         c = AllocWindingAccu(w->numpoints + 1);
484         c->numpoints = w->numpoints;
485         for (i = 0; i < c->numpoints; i++)
486         {
487                 VectorCopyAccu(w->p[i], c->p[i]);
488         }
489         FreeWindingAccu(w);
490         return c;
491 }
492
493 /*
494 ==================
495 CopyWindingAccuToRegular
496 ==================
497 */
498 winding_t       *CopyWindingAccuToRegular(winding_accu_t *w)
499 {
500         int             i;
501         winding_t       *c;
502
503         if (!w) Error("CopyWindingAccuToRegular: winding is NULL");
504
505         c = AllocWinding(w->numpoints);
506         c->numpoints = w->numpoints;
507         for (i = 0; i < c->numpoints; i++)
508         {
509                 VectorCopyAccuToRegular(w->p[i], c->p[i]);
510         }
511         return c;
512 }
513
514 /*
515 ==================
516 ReverseWinding
517 ==================
518 */
519 winding_t       *ReverseWinding (winding_t *w)
520 {
521         int                     i;
522         winding_t       *c;
523
524         c = AllocWinding (w->numpoints);
525         for (i=0 ; i<w->numpoints ; i++)
526         {
527                 VectorCopy (w->p[w->numpoints-1-i], c->p[i]);
528         }
529         c->numpoints = w->numpoints;
530         return c;
531 }
532
533
534 /*
535 =============
536 ClipWindingEpsilon
537 =============
538 */
539 void    ClipWindingEpsilonStrict (winding_t *in, vec3_t normal, vec_t dist, 
540                                 vec_t epsilon, winding_t **front, winding_t **back)
541 {
542         vec_t   dists[MAX_POINTS_ON_WINDING+4];
543         int             sides[MAX_POINTS_ON_WINDING+4];
544         int             counts[3];
545         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
546         int             i, j;
547         vec_t   *p1, *p2;
548         vec3_t  mid;
549         winding_t       *f, *b;
550         int             maxpts;
551         
552         counts[0] = counts[1] = counts[2] = 0;
553
554 // determine sides for each point
555         for (i=0 ; i<in->numpoints ; i++)
556   {
557
558                 dot = DotProduct (in->p[i], normal);
559                 dot -= dist;
560                 dists[i] = dot;
561                 if (dot > epsilon)
562                         sides[i] = SIDE_FRONT;
563                 else if (dot < -epsilon)
564                         sides[i] = SIDE_BACK;
565                 else
566                 {
567                         sides[i] = SIDE_ON;
568                 }
569                 counts[sides[i]]++;
570         }
571         sides[i] = sides[0];
572         dists[i] = dists[0];
573         
574         *front = *back = NULL;
575
576         if (!counts[0] && !counts[1])
577         {
578                 return;
579         }
580         if (!counts[0])
581         {
582                 *back = CopyWinding (in);
583                 return;
584         }
585         if (!counts[1])
586         {
587                 *front = CopyWinding (in);
588                 return;
589         }
590
591         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
592                                                                 // of fp grouping errors
593
594         *front = f = AllocWinding (maxpts);
595         *back = b = AllocWinding (maxpts);
596                 
597         for (i=0 ; i<in->numpoints ; i++)
598         {
599                 p1 = in->p[i];
600                 
601                 if (sides[i] == SIDE_ON)
602                 {
603                         VectorCopy (p1, f->p[f->numpoints]);
604                         f->numpoints++;
605                         VectorCopy (p1, b->p[b->numpoints]);
606                         b->numpoints++;
607                         continue;
608                 }
609         
610                 if (sides[i] == SIDE_FRONT)
611                 {
612                         VectorCopy (p1, f->p[f->numpoints]);
613                         f->numpoints++;
614                 }
615                 if (sides[i] == SIDE_BACK)
616                 {
617                         VectorCopy (p1, b->p[b->numpoints]);
618                         b->numpoints++;
619                 }
620
621                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
622                         continue;
623                         
624         // generate a split point
625                 p2 = in->p[(i+1)%in->numpoints];
626                 
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)
631                                 mid[j] = dist;
632                         else if (normal[j] == -1)
633                                 mid[j] = -dist;
634                         else
635                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
636                 }
637                         
638                 VectorCopy (mid, f->p[f->numpoints]);
639                 f->numpoints++;
640                 VectorCopy (mid, b->p[b->numpoints]);
641                 b->numpoints++;
642         }
643         
644         if (f->numpoints > maxpts || b->numpoints > maxpts)
645                 Error ("ClipWinding: points exceeded estimate");
646         if (f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING)
647                 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
648 }
649
650 void    ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist, 
651                                 vec_t epsilon, winding_t **front, winding_t **back)
652 {
653         ClipWindingEpsilonStrict(in, normal, dist, epsilon, front, back);
654         /* apparently most code expects that in the winding-on-plane case, the back winding is the original winding */
655         if(!*front && !*back)
656                 *back = CopyWinding(in);
657 }
658
659
660 /*
661 =============
662 ChopWindingInPlaceAccu
663 =============
664 */
665 void ChopWindingInPlaceAccu(winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon)
666 {
667         vec_accu_t      fineEpsilon;
668         winding_accu_t  *in;
669         int             counts[3];
670         int             i, j;
671         vec_accu_t      dists[MAX_POINTS_ON_WINDING + 1];
672         int             sides[MAX_POINTS_ON_WINDING + 1];
673         int             maxpts;
674         winding_accu_t  *f;
675         vec_accu_t      *p1, *p2;
676         vec_accu_t      w;
677         vec3_accu_t     mid, normalAccu;
678
679         // We require at least a very small epsilon.  It's a good idea for several reasons.
680         // First, we will be dividing by a potentially very small distance below.  We don't
681         // want that distance to be too small; otherwise, things "blow up" with little accuracy
682         // due to the division.  (After a second look, the value w below is in range (0,1), but
683         // graininess problem remains.)  Second, Having minimum epsilon also prevents the following
684         // situation.  Say for example we have a perfect octagon defined by the input winding.
685         // Say our chopping plane (defined by normal and dist) is essentially the same plane
686         // that the octagon is sitting on.  Well, due to rounding errors, it may be that point
687         // 1 of the octagon might be in front, point 2 might be in back, point 3 might be in
688         // front, point 4 might be in back, and so on.  So we could end up with a very ugly-
689         // looking chopped winding, and this might be undesirable, and would at least lead to
690         // a possible exhaustion of MAX_POINTS_ON_WINDING.  It's better to assume that points
691         // very very close to the plane are on the plane, using an infinitesimal epsilon amount.
692
693         // Now, the original ChopWindingInPlace() function used a vec_t-based winding_t.
694         // So this minimum epsilon is quite similar to casting the higher resolution numbers to
695         // the lower resolution and comparing them in the lower resolution mode.  We explicitly
696         // choose the minimum epsilon as something around the vec_t epsilon of one because we
697         // want the resolution of vec_accu_t to have a large resolution around the epsilon.
698         // Some of that leftover resolution even goes away after we scale to points far away.
699
700         // Here is a further discussion regarding the choice of smallestEpsilonAllowed.
701         // In the 32 float world (we can assume vec_t is that), the "epsilon around 1.0" is
702         // 0.00000011921.  In the 64 bit float world (we can assume vec_accu_t is that), the
703         // "epsilon around 1.0" is 0.00000000000000022204.  (By the way these two epsilons
704         // are defined as VEC_SMALLEST_EPSILON_AROUND_ONE VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE
705         // respectively.)  If you divide the first by the second, you get approximately
706         // 536,885,246.  Dividing that number by 200,000 (a typical base winding coordinate)
707         // gives 2684.  So in other words, if our smallestEpsilonAllowed was chosen as exactly
708         // VEC_SMALLEST_EPSILON_AROUND_ONE, you would be guaranteed at least 2000 "ticks" in
709         // 64-bit land inside of the epsilon for all numbers we're dealing with.
710
711         static const vec_accu_t smallestEpsilonAllowed = ((vec_accu_t) VEC_SMALLEST_EPSILON_AROUND_ONE) * 0.5;
712         if (crudeEpsilon < smallestEpsilonAllowed)      fineEpsilon = smallestEpsilonAllowed;
713         else                                            fineEpsilon = (vec_accu_t) crudeEpsilon;
714
715         in = *inout;
716         counts[0] = counts[1] = counts[2] = 0;
717         VectorCopyRegularToAccu(normal, normalAccu);
718
719         for (i = 0; i < in->numpoints; i++)
720         {
721                 dists[i] = DotProductAccu(in->p[i], normalAccu) - dist;
722                 if (dists[i] > fineEpsilon)             sides[i] = SIDE_FRONT;
723                 else if (dists[i] < -fineEpsilon)       sides[i] = SIDE_BACK;
724                 else                                    sides[i] = SIDE_ON;
725                 counts[sides[i]]++;
726         }
727         sides[i] = sides[0];
728         dists[i] = dists[0];
729         
730         // I'm wondering if whatever code that handles duplicate planes is robust enough
731         // that we never get a case where two nearly equal planes result in 2 NULL windings
732         // due to the 'if' statement below.  TODO: Investigate this.
733         if (!counts[SIDE_FRONT]) {
734                 FreeWindingAccu(in);
735                 *inout = NULL;
736                 return;
737         }
738         if (!counts[SIDE_BACK]) {
739                 return; // Winding is unmodified.
740         }
741
742         // NOTE: The least number of points that a winding can have at this point is 2.
743         // In that case, one point is SIDE_FRONT and the other is SIDE_BACK.
744
745         maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
746         f = AllocWindingAccu(maxpts);
747
748         for (i = 0; i < in->numpoints; i++)
749         {
750                 p1 = in->p[i];
751
752                 if (sides[i] == SIDE_ON || sides[i] == SIDE_FRONT)
753                 {
754                         if (f->numpoints >= MAX_POINTS_ON_WINDING)
755                                 Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
756                         if (f->numpoints >= maxpts) // This will probably never happen.
757                         {
758                                 Sys_FPrintf(SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n");
759                                 f = CopyWindingAccuIncreaseSizeAndFreeOld(f);
760                                 maxpts++;
761                         }
762                         VectorCopyAccu(p1, f->p[f->numpoints]);
763                         f->numpoints++;
764                         if (sides[i] == SIDE_ON) continue;
765                 }
766                 if (sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i])
767                 {
768                         continue;
769                 }
770                         
771                 // Generate a split point.
772                 p2 = in->p[((i + 1) == in->numpoints) ? 0 : (i + 1)];
773
774                 // The divisor's absolute value is greater than the dividend's absolute value.
775                 // w is in the range (0,1).
776                 w = dists[i] / (dists[i] - dists[i + 1]);
777
778                 for (j = 0; j < 3; j++)
779                 {
780                         // Avoid round-off error when possible.  Check axis-aligned normal.
781                         if (normal[j] == 1)             mid[j] = dist;
782                         else if (normal[j] == -1)       mid[j] = -dist;
783                         else                            mid[j] = p1[j] + (w * (p2[j] - p1[j]));
784                 }
785                 if (f->numpoints >= MAX_POINTS_ON_WINDING)
786                         Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
787                 if (f->numpoints >= maxpts) // This will probably never happen.
788                 {
789                         Sys_FPrintf(SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n");
790                         f = CopyWindingAccuIncreaseSizeAndFreeOld(f);
791                         maxpts++;
792                 }
793                 VectorCopyAccu(mid, f->p[f->numpoints]);
794                 f->numpoints++;
795         }
796
797         FreeWindingAccu(in);
798         *inout = f;
799 }
800
801 /*
802 =============
803 ChopWindingInPlace
804 =============
805 */
806 void ChopWindingInPlace (winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon)
807 {
808         winding_t       *in;
809         vec_t   dists[MAX_POINTS_ON_WINDING+4];
810         int             sides[MAX_POINTS_ON_WINDING+4];
811         int             counts[3];
812         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
813         int             i, j;
814         vec_t   *p1, *p2;
815         vec3_t  mid;
816         winding_t       *f;
817         int             maxpts;
818
819         in = *inout;
820         counts[0] = counts[1] = counts[2] = 0;
821
822 // determine sides for each point
823         for (i=0 ; i<in->numpoints ; i++)
824         {
825                 dot = DotProduct (in->p[i], normal);
826                 dot -= dist;
827                 dists[i] = dot;
828                 if (dot > epsilon)
829                         sides[i] = SIDE_FRONT;
830                 else if (dot < -epsilon)
831                         sides[i] = SIDE_BACK;
832                 else
833                 {
834                         sides[i] = SIDE_ON;
835                 }
836                 counts[sides[i]]++;
837         }
838         sides[i] = sides[0];
839         dists[i] = dists[0];
840         
841         if (!counts[0])
842         {
843                 FreeWinding (in);
844                 *inout = NULL;
845                 return;
846         }
847         if (!counts[1])
848                 return;         // inout stays the same
849
850         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
851                                                                 // of fp grouping errors
852
853         f = AllocWinding (maxpts);
854                 
855         for (i=0 ; i<in->numpoints ; i++)
856         {
857                 p1 = in->p[i];
858                 
859                 if (sides[i] == SIDE_ON)
860                 {
861                         VectorCopy (p1, f->p[f->numpoints]);
862                         f->numpoints++;
863                         continue;
864                 }
865         
866                 if (sides[i] == SIDE_FRONT)
867                 {
868                         VectorCopy (p1, f->p[f->numpoints]);
869                         f->numpoints++;
870                 }
871
872                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
873                         continue;
874                         
875         // generate a split point
876                 p2 = in->p[(i+1)%in->numpoints];
877                 
878                 dot = dists[i] / (dists[i]-dists[i+1]);
879                 for (j=0 ; j<3 ; j++)
880                 {       // avoid round off error when possible
881                         if (normal[j] == 1)
882                                 mid[j] = dist;
883                         else if (normal[j] == -1)
884                                 mid[j] = -dist;
885                         else
886                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
887                 }
888                         
889                 VectorCopy (mid, f->p[f->numpoints]);
890                 f->numpoints++;
891         }
892         
893         if (f->numpoints > maxpts)
894                 Error ("ClipWinding: points exceeded estimate");
895         if (f->numpoints > MAX_POINTS_ON_WINDING)
896                 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
897
898         FreeWinding (in);
899         *inout = f;
900 }
901
902
903 /*
904 =================
905 ChopWinding
906
907 Returns the fragment of in that is on the front side
908 of the cliping plane.  The original is freed.
909 =================
910 */
911 winding_t       *ChopWinding (winding_t *in, vec3_t normal, vec_t dist)
912 {
913         winding_t       *f, *b;
914
915         ClipWindingEpsilon (in, normal, dist, ON_EPSILON, &f, &b);
916         FreeWinding (in);
917         if (b)
918                 FreeWinding (b);
919         return f;
920 }
921
922
923 /*
924 =================
925 CheckWinding
926
927 =================
928 */
929 void CheckWinding (winding_t *w)
930 {
931         int             i, j;
932         vec_t   *p1, *p2;
933         vec_t   d, edgedist;
934         vec3_t  dir, edgenormal, facenormal;
935         vec_t   area;
936         vec_t   facedist;
937
938         if (w->numpoints < 3)
939                 Error ("CheckWinding: %i points",w->numpoints);
940         
941         area = WindingArea(w);
942         if (area < 1)
943                 Error ("CheckWinding: %f area", area);
944
945         WindingPlane (w, facenormal, &facedist);
946         
947         for (i=0 ; i<w->numpoints ; i++)
948         {
949                 p1 = w->p[i];
950
951                 for (j=0 ; j<3 ; j++)
952                         if (p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD)
953                                 Error ("CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j]);
954
955                 j = i+1 == w->numpoints ? 0 : i+1;
956                 
957         // check the point is on the face plane
958                 d = DotProduct (p1, facenormal) - facedist;
959                 if (d < -ON_EPSILON || d > ON_EPSILON)
960                         Error ("CheckWinding: point off plane");
961         
962         // check the edge isnt degenerate
963                 p2 = w->p[j];
964                 VectorSubtract (p2, p1, dir);
965                 
966                 if (VectorLength (dir) < ON_EPSILON)
967                         Error ("CheckWinding: degenerate edge");
968                         
969                 CrossProduct (facenormal, dir, edgenormal);
970                 VectorNormalize (edgenormal, edgenormal);
971                 edgedist = DotProduct (p1, edgenormal);
972                 edgedist += ON_EPSILON;
973                 
974         // all other points must be on front side
975                 for (j=0 ; j<w->numpoints ; j++)
976                 {
977                         if (j == i)
978                                 continue;
979                         d = DotProduct (w->p[j], edgenormal);
980                         if (d > edgedist)
981                                 Error ("CheckWinding: non-convex");
982                 }
983         }
984 }
985
986
987 /*
988 ============
989 WindingOnPlaneSide
990 ============
991 */
992 int             WindingOnPlaneSide (winding_t *w, vec3_t normal, vec_t dist)
993 {
994         qboolean        front, back;
995         int                     i;
996         vec_t           d;
997
998         front = qfalse;
999         back = qfalse;
1000         for (i=0 ; i<w->numpoints ; i++)
1001         {
1002                 d = DotProduct (w->p[i], normal) - dist;
1003                 if (d < -ON_EPSILON)
1004                 {
1005                         if (front)
1006                                 return SIDE_CROSS;
1007                         back = qtrue;
1008                         continue;
1009                 }
1010                 if (d > ON_EPSILON)
1011                 {
1012                         if (back)
1013                                 return SIDE_CROSS;
1014                         front = qtrue;
1015                         continue;
1016                 }
1017         }
1018
1019         if (back)
1020                 return SIDE_BACK;
1021         if (front)
1022                 return SIDE_FRONT;
1023         return SIDE_ON;
1024 }
1025
1026
1027 /*
1028 =================
1029 AddWindingToConvexHull
1030
1031 Both w and *hull are on the same plane
1032 =================
1033 */
1034 #define MAX_HULL_POINTS         128
1035 void    AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
1036         int                     i, j, k;
1037         float           *p, *copy;
1038         vec3_t          dir;
1039         float           d;
1040         int                     numHullPoints, numNew;
1041         vec3_t          hullPoints[MAX_HULL_POINTS];
1042         vec3_t          newHullPoints[MAX_HULL_POINTS];
1043         vec3_t          hullDirs[MAX_HULL_POINTS];
1044         qboolean        hullSide[MAX_HULL_POINTS];
1045         qboolean        outside;
1046
1047         if ( !*hull ) {
1048                 *hull = CopyWinding( w );
1049                 return;
1050         }
1051
1052         numHullPoints = (*hull)->numpoints;
1053         memcpy( hullPoints, (*hull)->p, numHullPoints * sizeof(vec3_t) );
1054
1055         for ( i = 0 ; i < w->numpoints ; i++ ) {
1056                 p = w->p[i];
1057
1058                 // calculate hull side vectors
1059                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1060                         k = ( j + 1 ) % numHullPoints;
1061
1062                         VectorSubtract( hullPoints[k], hullPoints[j], dir );
1063                         VectorNormalize( dir, dir );
1064                         CrossProduct( normal, dir, hullDirs[j] );
1065                 }
1066
1067                 outside = qfalse;
1068                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1069                         VectorSubtract( p, hullPoints[j], dir );
1070                         d = DotProduct( dir, hullDirs[j] );
1071                         if ( d >= ON_EPSILON ) {
1072                                 outside = qtrue;
1073                         }
1074                         if ( d >= -ON_EPSILON ) {
1075                                 hullSide[j] = qtrue;
1076                         } else {
1077                                 hullSide[j] = qfalse;
1078                         }
1079                 }
1080
1081                 // if the point is effectively inside, do nothing
1082                 if ( !outside ) {
1083                         continue;
1084                 }
1085
1086                 // find the back side to front side transition
1087                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1088                         if ( !hullSide[ j % numHullPoints ] && hullSide[ (j + 1) % numHullPoints ] ) {
1089                                 break;
1090                         }
1091                 }
1092                 if ( j == numHullPoints ) {
1093                         continue;
1094                 }
1095
1096                 // insert the point here
1097                 VectorCopy( p, newHullPoints[0] );
1098                 numNew = 1;
1099
1100                 // copy over all points that aren't double fronts
1101                 j = (j+1)%numHullPoints;
1102                 for ( k = 0 ; k < numHullPoints ; k++ ) {
1103                         if ( hullSide[ (j+k) % numHullPoints ] && hullSide[ (j+k+1) % numHullPoints ] ) {
1104                                 continue;
1105                         }
1106                         copy = hullPoints[ (j+k+1) % numHullPoints ];
1107                         VectorCopy( copy, newHullPoints[numNew] );
1108                         numNew++;
1109                 }
1110
1111                 numHullPoints = numNew;
1112                 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof(vec3_t) );
1113         }
1114
1115         FreeWinding( *hull );
1116         w = AllocWinding( numHullPoints );
1117         w->numpoints = numHullPoints;
1118         *hull = w;
1119         memcpy( w->p, hullPoints, numHullPoints * sizeof(vec3_t) );
1120 }
1121
1122