::zerowing-base=422
[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 /*
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(vec_accu_t) * 3 * points + sizeof(int);
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    ClipWindingEpsilon (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])
577         {
578                 *back = CopyWinding (in);
579                 return;
580         }
581         if (!counts[1])
582         {
583                 *front = CopyWinding (in);
584                 return;
585         }
586
587         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
588                                                                 // of fp grouping errors
589
590         *front = f = AllocWinding (maxpts);
591         *back = b = AllocWinding (maxpts);
592                 
593         for (i=0 ; i<in->numpoints ; i++)
594         {
595                 p1 = in->p[i];
596                 
597                 if (sides[i] == SIDE_ON)
598                 {
599                         VectorCopy (p1, f->p[f->numpoints]);
600                         f->numpoints++;
601                         VectorCopy (p1, b->p[b->numpoints]);
602                         b->numpoints++;
603                         continue;
604                 }
605         
606                 if (sides[i] == SIDE_FRONT)
607                 {
608                         VectorCopy (p1, f->p[f->numpoints]);
609                         f->numpoints++;
610                 }
611                 if (sides[i] == SIDE_BACK)
612                 {
613                         VectorCopy (p1, b->p[b->numpoints]);
614                         b->numpoints++;
615                 }
616
617                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
618                         continue;
619                         
620         // generate a split point
621                 p2 = in->p[(i+1)%in->numpoints];
622                 
623                 dot = dists[i] / (dists[i]-dists[i+1]);
624                 for (j=0 ; j<3 ; j++)
625                 {       // avoid round off error when possible
626                         if (normal[j] == 1)
627                                 mid[j] = dist;
628                         else if (normal[j] == -1)
629                                 mid[j] = -dist;
630                         else
631                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
632                 }
633                         
634                 VectorCopy (mid, f->p[f->numpoints]);
635                 f->numpoints++;
636                 VectorCopy (mid, b->p[b->numpoints]);
637                 b->numpoints++;
638         }
639         
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");
644 }
645
646
647 /*
648 =============
649 ChopWindingInPlaceAccu
650 =============
651 */
652 void ChopWindingInPlaceAccu(winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon)
653 {
654         vec_accu_t      fineEpsilon;
655         winding_accu_t  *in;
656         int             counts[3];
657         int             i, j;
658         vec_accu_t      dists[MAX_POINTS_ON_WINDING + 1];
659         int             sides[MAX_POINTS_ON_WINDING + 1];
660         int             maxpts;
661         winding_accu_t  *f;
662         vec_accu_t      *p1, *p2;
663         vec_accu_t      w;
664         vec3_accu_t     mid, normalAccu;
665
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.
679
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.
686
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.
697
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;
701
702         in = *inout;
703         counts[0] = counts[1] = counts[2] = 0;
704         VectorCopyRegularToAccu(normal, normalAccu);
705
706         for (i = 0; i < in->numpoints; i++)
707         {
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;
712                 counts[sides[i]]++;
713         }
714         sides[i] = sides[0];
715         dists[i] = dists[0];
716         
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]) {
721                 FreeWindingAccu(in);
722                 *inout = NULL;
723                 return;
724         }
725         if (!counts[SIDE_BACK]) {
726                 return; // Winding is unmodified.
727         }
728
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.
731
732         maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
733         f = AllocWindingAccu(maxpts);
734
735         for (i = 0; i < in->numpoints; i++)
736         {
737                 p1 = in->p[i];
738
739                 if (sides[i] == SIDE_ON || sides[i] == SIDE_FRONT)
740                 {
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.
744                         {
745                                 Sys_FPrintf(SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n");
746                                 f = CopyWindingAccuIncreaseSizeAndFreeOld(f);
747                                 maxpts++;
748                         }
749                         VectorCopyAccu(p1, f->p[f->numpoints]);
750                         f->numpoints++;
751                         if (sides[i] == SIDE_ON) continue;
752                 }
753                 if (sides[i + 1] == SIDE_ON || sides[i + 1] == sides[i])
754                 {
755                         continue;
756                 }
757                         
758                 // Generate a split point.
759                 p2 = in->p[((i + 1) == in->numpoints) ? 0 : (i + 1)];
760
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]);
764
765                 for (j = 0; j < 3; j++)
766                 {
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]));
771                 }
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.
775                 {
776                         Sys_FPrintf(SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n");
777                         f = CopyWindingAccuIncreaseSizeAndFreeOld(f);
778                         maxpts++;
779                 }
780                 VectorCopyAccu(mid, f->p[f->numpoints]);
781                 f->numpoints++;
782         }
783
784         FreeWindingAccu(in);
785         *inout = f;
786 }
787
788 /*
789 =============
790 ChopWindingInPlace
791 =============
792 */
793 void ChopWindingInPlace (winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon)
794 {
795         winding_t       *in;
796         vec_t   dists[MAX_POINTS_ON_WINDING+4];
797         int             sides[MAX_POINTS_ON_WINDING+4];
798         int             counts[3];
799         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
800         int             i, j;
801         vec_t   *p1, *p2;
802         vec3_t  mid;
803         winding_t       *f;
804         int             maxpts;
805
806         in = *inout;
807         counts[0] = counts[1] = counts[2] = 0;
808
809 // determine sides for each point
810         for (i=0 ; i<in->numpoints ; i++)
811         {
812                 dot = DotProduct (in->p[i], normal);
813                 dot -= dist;
814                 dists[i] = dot;
815                 if (dot > epsilon)
816                         sides[i] = SIDE_FRONT;
817                 else if (dot < -epsilon)
818                         sides[i] = SIDE_BACK;
819                 else
820                 {
821                         sides[i] = SIDE_ON;
822                 }
823                 counts[sides[i]]++;
824         }
825         sides[i] = sides[0];
826         dists[i] = dists[0];
827         
828         if (!counts[0])
829         {
830                 FreeWinding (in);
831                 *inout = NULL;
832                 return;
833         }
834         if (!counts[1])
835                 return;         // inout stays the same
836
837         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
838                                                                 // of fp grouping errors
839
840         f = AllocWinding (maxpts);
841                 
842         for (i=0 ; i<in->numpoints ; i++)
843         {
844                 p1 = in->p[i];
845                 
846                 if (sides[i] == SIDE_ON)
847                 {
848                         VectorCopy (p1, f->p[f->numpoints]);
849                         f->numpoints++;
850                         continue;
851                 }
852         
853                 if (sides[i] == SIDE_FRONT)
854                 {
855                         VectorCopy (p1, f->p[f->numpoints]);
856                         f->numpoints++;
857                 }
858
859                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
860                         continue;
861                         
862         // generate a split point
863                 p2 = in->p[(i+1)%in->numpoints];
864                 
865                 dot = dists[i] / (dists[i]-dists[i+1]);
866                 for (j=0 ; j<3 ; j++)
867                 {       // avoid round off error when possible
868                         if (normal[j] == 1)
869                                 mid[j] = dist;
870                         else if (normal[j] == -1)
871                                 mid[j] = -dist;
872                         else
873                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
874                 }
875                         
876                 VectorCopy (mid, f->p[f->numpoints]);
877                 f->numpoints++;
878         }
879         
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");
884
885         FreeWinding (in);
886         *inout = f;
887 }
888
889
890 /*
891 =================
892 ChopWinding
893
894 Returns the fragment of in that is on the front side
895 of the cliping plane.  The original is freed.
896 =================
897 */
898 winding_t       *ChopWinding (winding_t *in, vec3_t normal, vec_t dist)
899 {
900         winding_t       *f, *b;
901
902         ClipWindingEpsilon (in, normal, dist, ON_EPSILON, &f, &b);
903         FreeWinding (in);
904         if (b)
905                 FreeWinding (b);
906         return f;
907 }
908
909
910 /*
911 =================
912 CheckWinding
913
914 =================
915 */
916 void CheckWinding (winding_t *w)
917 {
918         int             i, j;
919         vec_t   *p1, *p2;
920         vec_t   d, edgedist;
921         vec3_t  dir, edgenormal, facenormal;
922         vec_t   area;
923         vec_t   facedist;
924
925         if (w->numpoints < 3)
926                 Error ("CheckWinding: %i points",w->numpoints);
927         
928         area = WindingArea(w);
929         if (area < 1)
930                 Error ("CheckWinding: %f area", area);
931
932         WindingPlane (w, facenormal, &facedist);
933         
934         for (i=0 ; i<w->numpoints ; i++)
935         {
936                 p1 = w->p[i];
937
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]);
941
942                 j = i+1 == w->numpoints ? 0 : i+1;
943                 
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");
948         
949         // check the edge isnt degenerate
950                 p2 = w->p[j];
951                 VectorSubtract (p2, p1, dir);
952                 
953                 if (VectorLength (dir) < ON_EPSILON)
954                         Error ("CheckWinding: degenerate edge");
955                         
956                 CrossProduct (facenormal, dir, edgenormal);
957                 VectorNormalize (edgenormal, edgenormal);
958                 edgedist = DotProduct (p1, edgenormal);
959                 edgedist += ON_EPSILON;
960                 
961         // all other points must be on front side
962                 for (j=0 ; j<w->numpoints ; j++)
963                 {
964                         if (j == i)
965                                 continue;
966                         d = DotProduct (w->p[j], edgenormal);
967                         if (d > edgedist)
968                                 Error ("CheckWinding: non-convex");
969                 }
970         }
971 }
972
973
974 /*
975 ============
976 WindingOnPlaneSide
977 ============
978 */
979 int             WindingOnPlaneSide (winding_t *w, vec3_t normal, vec_t dist)
980 {
981         qboolean        front, back;
982         int                     i;
983         vec_t           d;
984
985         front = qfalse;
986         back = qfalse;
987         for (i=0 ; i<w->numpoints ; i++)
988         {
989                 d = DotProduct (w->p[i], normal) - dist;
990                 if (d < -ON_EPSILON)
991                 {
992                         if (front)
993                                 return SIDE_CROSS;
994                         back = qtrue;
995                         continue;
996                 }
997                 if (d > ON_EPSILON)
998                 {
999                         if (back)
1000                                 return SIDE_CROSS;
1001                         front = qtrue;
1002                         continue;
1003                 }
1004         }
1005
1006         if (back)
1007                 return SIDE_BACK;
1008         if (front)
1009                 return SIDE_FRONT;
1010         return SIDE_ON;
1011 }
1012
1013
1014 /*
1015 =================
1016 AddWindingToConvexHull
1017
1018 Both w and *hull are on the same plane
1019 =================
1020 */
1021 #define MAX_HULL_POINTS         128
1022 void    AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
1023         int                     i, j, k;
1024         float           *p, *copy;
1025         vec3_t          dir;
1026         float           d;
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];
1032         qboolean        outside;
1033
1034         if ( !*hull ) {
1035                 *hull = CopyWinding( w );
1036                 return;
1037         }
1038
1039         numHullPoints = (*hull)->numpoints;
1040         memcpy( hullPoints, (*hull)->p, numHullPoints * sizeof(vec3_t) );
1041
1042         for ( i = 0 ; i < w->numpoints ; i++ ) {
1043                 p = w->p[i];
1044
1045                 // calculate hull side vectors
1046                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1047                         k = ( j + 1 ) % numHullPoints;
1048
1049                         VectorSubtract( hullPoints[k], hullPoints[j], dir );
1050                         VectorNormalize( dir, dir );
1051                         CrossProduct( normal, dir, hullDirs[j] );
1052                 }
1053
1054                 outside = qfalse;
1055                 for ( j = 0 ; j < numHullPoints ; j++ ) {
1056                         VectorSubtract( p, hullPoints[j], dir );
1057                         d = DotProduct( dir, hullDirs[j] );
1058                         if ( d >= ON_EPSILON ) {
1059                                 outside = qtrue;
1060                         }
1061                         if ( d >= -ON_EPSILON ) {
1062                                 hullSide[j] = qtrue;
1063                         } else {
1064                                 hullSide[j] = qfalse;
1065                         }
1066                 }
1067
1068                 // if the point is effectively inside, do nothing
1069                 if ( !outside ) {
1070                         continue;
1071                 }
1072
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 ] ) {
1076                                 break;
1077                         }
1078                 }
1079                 if ( j == numHullPoints ) {
1080                         continue;
1081                 }
1082
1083                 // insert the point here
1084                 VectorCopy( p, newHullPoints[0] );
1085                 numNew = 1;
1086
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 ] ) {
1091                                 continue;
1092                         }
1093                         copy = hullPoints[ (j+k+1) % numHullPoints ];
1094                         VectorCopy( copy, newHullPoints[numNew] );
1095                         numNew++;
1096                 }
1097
1098                 numHullPoints = numNew;
1099                 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof(vec3_t) );
1100         }
1101
1102         FreeWinding( *hull );
1103         w = AllocWinding( numHullPoints );
1104         w->numpoints = numHullPoints;
1105         *hull = w;
1106         memcpy( w->p, hullPoints, numHullPoints * sizeof(vec3_t) );
1107 }
1108
1109