Author: 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         int             i, x;
212         vec_t   max, v;
213         vec3_t  org, vright, vup;
214         winding_t       *w;
215         
216 // find the major axis
217
218         max = -BOGUS_RANGE;
219         x = -1;
220         for (i=0 ; i<3; i++)
221         {
222                 v = fabs(normal[i]);
223                 if (v > max)
224                 {
225                         x = i;
226                         max = v;
227                 }
228         }
229         if (x==-1)
230                 Error ("BaseWindingForPlane: no axis found");
231                 
232         VectorCopy (vec3_origin, vup);  
233         switch (x)
234         {
235         case 0:
236         case 1:
237                 vup[2] = 1;
238                 break;          
239         case 2:
240                 vup[0] = 1;
241                 break;          
242         }
243
244         v = DotProduct (vup, normal);
245         VectorMA (vup, -v, normal, vup);
246         VectorNormalize (vup, vup);
247                 
248         VectorScale (normal, dist, org);
249         
250         CrossProduct (vup, normal, vright);
251         
252         // LordHavoc: this has to use *2 because otherwise some created points may
253         // be inside the world (think of a diagonal case), and any brush with such
254         // points should be removed, failure to detect such cases is disasterous
255         VectorScale (vup, MAX_WORLD_COORD*2, vup);
256         VectorScale (vright, MAX_WORLD_COORD*2, vright);
257
258   // project a really big       axis aligned box onto the plane
259         w = AllocWinding (4);
260         
261         VectorSubtract (org, vright, w->p[0]);
262         VectorAdd (w->p[0], vup, w->p[0]);
263         
264         VectorAdd (org, vright, w->p[1]);
265         VectorAdd (w->p[1], vup, w->p[1]);
266         
267         VectorAdd (org, vright, w->p[2]);
268         VectorSubtract (w->p[2], vup, w->p[2]);
269         
270         VectorSubtract (org, vright, w->p[3]);
271         VectorSubtract (w->p[3], vup, w->p[3]);
272         
273         w->numpoints = 4;
274         
275         return w;       
276 }
277
278 /*
279 ==================
280 CopyWinding
281 ==================
282 */
283 winding_t       *CopyWinding (winding_t *w)
284 {
285         size_t                  size;
286         winding_t       *c;
287
288         c = AllocWinding (w->numpoints);
289         size = (size_t)((winding_t *)NULL)->p[w->numpoints];
290         memcpy (c, w, size);
291         return c;
292 }
293
294 /*
295 ==================
296 ReverseWinding
297 ==================
298 */
299 winding_t       *ReverseWinding (winding_t *w)
300 {
301         int                     i;
302         winding_t       *c;
303
304         c = AllocWinding (w->numpoints);
305         for (i=0 ; i<w->numpoints ; i++)
306         {
307                 VectorCopy (w->p[w->numpoints-1-i], c->p[i]);
308         }
309         c->numpoints = w->numpoints;
310         return c;
311 }
312
313
314 /*
315 =============
316 ClipWindingEpsilon
317 =============
318 */
319 void    ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist, 
320                                 vec_t epsilon, winding_t **front, winding_t **back)
321 {
322         vec_t   dists[MAX_POINTS_ON_WINDING+4];
323         int             sides[MAX_POINTS_ON_WINDING+4];
324         int             counts[3];
325         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
326         int             i, j;
327         vec_t   *p1, *p2;
328         vec3_t  mid;
329         winding_t       *f, *b;
330         int             maxpts;
331         
332         counts[0] = counts[1] = counts[2] = 0;
333
334 // determine sides for each point
335         for (i=0 ; i<in->numpoints ; i++)
336   {
337
338                 dot = DotProduct (in->p[i], normal);
339                 dot -= dist;
340                 dists[i] = dot;
341                 if (dot > epsilon)
342                         sides[i] = SIDE_FRONT;
343                 else if (dot < -epsilon)
344                         sides[i] = SIDE_BACK;
345                 else
346                 {
347                         sides[i] = SIDE_ON;
348                 }
349                 counts[sides[i]]++;
350         }
351         sides[i] = sides[0];
352         dists[i] = dists[0];
353         
354         *front = *back = NULL;
355
356         if (!counts[0])
357         {
358                 *back = CopyWinding (in);
359                 return;
360         }
361         if (!counts[1])
362         {
363                 *front = CopyWinding (in);
364                 return;
365         }
366
367         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
368                                                                 // of fp grouping errors
369
370         *front = f = AllocWinding (maxpts);
371         *back = b = AllocWinding (maxpts);
372                 
373         for (i=0 ; i<in->numpoints ; i++)
374         {
375                 p1 = in->p[i];
376                 
377                 if (sides[i] == SIDE_ON)
378                 {
379                         VectorCopy (p1, f->p[f->numpoints]);
380                         f->numpoints++;
381                         VectorCopy (p1, b->p[b->numpoints]);
382                         b->numpoints++;
383                         continue;
384                 }
385         
386                 if (sides[i] == SIDE_FRONT)
387                 {
388                         VectorCopy (p1, f->p[f->numpoints]);
389                         f->numpoints++;
390                 }
391                 if (sides[i] == SIDE_BACK)
392                 {
393                         VectorCopy (p1, b->p[b->numpoints]);
394                         b->numpoints++;
395                 }
396
397                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
398                         continue;
399                         
400         // generate a split point
401                 p2 = in->p[(i+1)%in->numpoints];
402                 
403                 dot = dists[i] / (dists[i]-dists[i+1]);
404                 for (j=0 ; j<3 ; j++)
405                 {       // avoid round off error when possible
406                         if (normal[j] == 1)
407                                 mid[j] = dist;
408                         else if (normal[j] == -1)
409                                 mid[j] = -dist;
410                         else
411                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
412                 }
413                         
414                 VectorCopy (mid, f->p[f->numpoints]);
415                 f->numpoints++;
416                 VectorCopy (mid, b->p[b->numpoints]);
417                 b->numpoints++;
418         }
419         
420         if (f->numpoints > maxpts || b->numpoints > maxpts)
421                 Error ("ClipWinding: points exceeded estimate");
422         if (f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING)
423                 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
424 }
425
426
427 /*
428 =============
429 ChopWindingInPlace
430 =============
431 */
432 void ChopWindingInPlace (winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon)
433 {
434         winding_t       *in;
435         vec_t   dists[MAX_POINTS_ON_WINDING+4];
436         int             sides[MAX_POINTS_ON_WINDING+4];
437         int             counts[3];
438         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
439         int             i, j;
440         vec_t   *p1, *p2;
441         vec3_t  mid;
442         winding_t       *f;
443         int             maxpts;
444
445         in = *inout;
446         counts[0] = counts[1] = counts[2] = 0;
447
448 // determine sides for each point
449         for (i=0 ; i<in->numpoints ; i++)
450         {
451                 dot = DotProduct (in->p[i], normal);
452                 dot -= dist;
453                 dists[i] = dot;
454                 if (dot > epsilon)
455                         sides[i] = SIDE_FRONT;
456                 else if (dot < -epsilon)
457                         sides[i] = SIDE_BACK;
458                 else
459                 {
460                         sides[i] = SIDE_ON;
461                 }
462                 counts[sides[i]]++;
463         }
464         sides[i] = sides[0];
465         dists[i] = dists[0];
466         
467         if (!counts[0])
468         {
469                 FreeWinding (in);
470                 *inout = NULL;
471                 return;
472         }
473         if (!counts[1])
474                 return;         // inout stays the same
475
476         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
477                                                                 // of fp grouping errors
478
479         f = AllocWinding (maxpts);
480                 
481         for (i=0 ; i<in->numpoints ; i++)
482         {
483                 p1 = in->p[i];
484                 
485                 if (sides[i] == SIDE_ON)
486                 {
487                         VectorCopy (p1, f->p[f->numpoints]);
488                         f->numpoints++;
489                         continue;
490                 }
491         
492                 if (sides[i] == SIDE_FRONT)
493                 {
494                         VectorCopy (p1, f->p[f->numpoints]);
495                         f->numpoints++;
496                 }
497
498                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
499                         continue;
500                         
501         // generate a split point
502                 p2 = in->p[(i+1)%in->numpoints];
503                 
504                 dot = dists[i] / (dists[i]-dists[i+1]);
505                 for (j=0 ; j<3 ; j++)
506                 {       // avoid round off error when possible
507                         if (normal[j] == 1)
508                                 mid[j] = dist;
509                         else if (normal[j] == -1)
510                                 mid[j] = -dist;
511                         else
512                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
513                 }
514                         
515                 VectorCopy (mid, f->p[f->numpoints]);
516                 f->numpoints++;
517         }
518         
519         if (f->numpoints > maxpts)
520                 Error ("ClipWinding: points exceeded estimate");
521         if (f->numpoints > MAX_POINTS_ON_WINDING)
522                 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
523
524         FreeWinding (in);
525         *inout = f;
526 }
527
528
529 /*
530 =================
531 ChopWinding
532
533 Returns the fragment of in that is on the front side
534 of the cliping plane.  The original is freed.
535 =================
536 */
537 winding_t       *ChopWinding (winding_t *in, vec3_t normal, vec_t dist)
538 {
539         winding_t       *f, *b;
540
541         ClipWindingEpsilon (in, normal, dist, ON_EPSILON, &f, &b);
542         FreeWinding (in);
543         if (b)
544                 FreeWinding (b);
545         return f;
546 }
547
548
549 /*
550 =================
551 CheckWinding
552
553 =================
554 */
555 void CheckWinding (winding_t *w)
556 {
557         int             i, j;
558         vec_t   *p1, *p2;
559         vec_t   d, edgedist;
560         vec3_t  dir, edgenormal, facenormal;
561         vec_t   area;
562         vec_t   facedist;
563
564         if (w->numpoints < 3)
565                 Error ("CheckWinding: %i points",w->numpoints);
566         
567         area = WindingArea(w);
568         if (area < 1)
569                 Error ("CheckWinding: %f area", area);
570
571         WindingPlane (w, facenormal, &facedist);
572         
573         for (i=0 ; i<w->numpoints ; i++)
574         {
575                 p1 = w->p[i];
576
577                 for (j=0 ; j<3 ; j++)
578                         if (p1[j] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD)
579                                 Error ("CheckFace: MAX_WORLD_COORD exceeded: %f",p1[j]);
580
581                 j = i+1 == w->numpoints ? 0 : i+1;
582                 
583         // check the point is on the face plane
584                 d = DotProduct (p1, facenormal) - facedist;
585                 if (d < -ON_EPSILON || d > ON_EPSILON)
586                         Error ("CheckWinding: point off plane");
587         
588         // check the edge isnt degenerate
589                 p2 = w->p[j];
590                 VectorSubtract (p2, p1, dir);
591                 
592                 if (VectorLength (dir) < ON_EPSILON)
593                         Error ("CheckWinding: degenerate edge");
594                         
595                 CrossProduct (facenormal, dir, edgenormal);
596                 VectorNormalize (edgenormal, edgenormal);
597                 edgedist = DotProduct (p1, edgenormal);
598                 edgedist += ON_EPSILON;
599                 
600         // all other points must be on front side
601                 for (j=0 ; j<w->numpoints ; j++)
602                 {
603                         if (j == i)
604                                 continue;
605                         d = DotProduct (w->p[j], edgenormal);
606                         if (d > edgedist)
607                                 Error ("CheckWinding: non-convex");
608                 }
609         }
610 }
611
612
613 /*
614 ============
615 WindingOnPlaneSide
616 ============
617 */
618 int             WindingOnPlaneSide (winding_t *w, vec3_t normal, vec_t dist)
619 {
620         qboolean        front, back;
621         int                     i;
622         vec_t           d;
623
624         front = qfalse;
625         back = qfalse;
626         for (i=0 ; i<w->numpoints ; i++)
627         {
628                 d = DotProduct (w->p[i], normal) - dist;
629                 if (d < -ON_EPSILON)
630                 {
631                         if (front)
632                                 return SIDE_CROSS;
633                         back = qtrue;
634                         continue;
635                 }
636                 if (d > ON_EPSILON)
637                 {
638                         if (back)
639                                 return SIDE_CROSS;
640                         front = qtrue;
641                         continue;
642                 }
643         }
644
645         if (back)
646                 return SIDE_BACK;
647         if (front)
648                 return SIDE_FRONT;
649         return SIDE_ON;
650 }
651
652
653 /*
654 =================
655 AddWindingToConvexHull
656
657 Both w and *hull are on the same plane
658 =================
659 */
660 #define MAX_HULL_POINTS         128
661 void    AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
662         int                     i, j, k;
663         float           *p, *copy;
664         vec3_t          dir;
665         float           d;
666         int                     numHullPoints, numNew;
667         vec3_t          hullPoints[MAX_HULL_POINTS];
668         vec3_t          newHullPoints[MAX_HULL_POINTS];
669         vec3_t          hullDirs[MAX_HULL_POINTS];
670         qboolean        hullSide[MAX_HULL_POINTS];
671         qboolean        outside;
672
673         if ( !*hull ) {
674                 *hull = CopyWinding( w );
675                 return;
676         }
677
678         numHullPoints = (*hull)->numpoints;
679         memcpy( hullPoints, (*hull)->p, numHullPoints * sizeof(vec3_t) );
680
681         for ( i = 0 ; i < w->numpoints ; i++ ) {
682                 p = w->p[i];
683
684                 // calculate hull side vectors
685                 for ( j = 0 ; j < numHullPoints ; j++ ) {
686                         k = ( j + 1 ) % numHullPoints;
687
688                         VectorSubtract( hullPoints[k], hullPoints[j], dir );
689                         VectorNormalize( dir, dir );
690                         CrossProduct( normal, dir, hullDirs[j] );
691                 }
692
693                 outside = qfalse;
694                 for ( j = 0 ; j < numHullPoints ; j++ ) {
695                         VectorSubtract( p, hullPoints[j], dir );
696                         d = DotProduct( dir, hullDirs[j] );
697                         if ( d >= ON_EPSILON ) {
698                                 outside = qtrue;
699                         }
700                         if ( d >= -ON_EPSILON ) {
701                                 hullSide[j] = qtrue;
702                         } else {
703                                 hullSide[j] = qfalse;
704                         }
705                 }
706
707                 // if the point is effectively inside, do nothing
708                 if ( !outside ) {
709                         continue;
710                 }
711
712                 // find the back side to front side transition
713                 for ( j = 0 ; j < numHullPoints ; j++ ) {
714                         if ( !hullSide[ j % numHullPoints ] && hullSide[ (j + 1) % numHullPoints ] ) {
715                                 break;
716                         }
717                 }
718                 if ( j == numHullPoints ) {
719                         continue;
720                 }
721
722                 // insert the point here
723                 VectorCopy( p, newHullPoints[0] );
724                 numNew = 1;
725
726                 // copy over all points that aren't double fronts
727                 j = (j+1)%numHullPoints;
728                 for ( k = 0 ; k < numHullPoints ; k++ ) {
729                         if ( hullSide[ (j+k) % numHullPoints ] && hullSide[ (j+k+1) % numHullPoints ] ) {
730                                 continue;
731                         }
732                         copy = hullPoints[ (j+k+1) % numHullPoints ];
733                         VectorCopy( copy, newHullPoints[numNew] );
734                         numNew++;
735                 }
736
737                 numHullPoints = numNew;
738                 memcpy( hullPoints, newHullPoints, numHullPoints * sizeof(vec3_t) );
739         }
740
741         FreeWinding( *hull );
742         w = AllocWinding( numHullPoints );
743         w->numpoints = numHullPoints;
744         *hull = w;
745         memcpy( w->p, hullPoints, numHullPoints * sizeof(vec3_t) );
746 }
747
748