Merge commit '515673c08f8718a237e90c2130a1f5294f966d6a'
[xonotic/netradiant.git] / tools / quake2 / 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 #include "cmdlib.h"
23 #include "inout.h"
24 #include "mathlib.h"
25 #include "polylib.h"
26
27
28 extern int numthreads;
29
30 // counters are only bumped when running single threaded,
31 // because they are an awefull coherence problem
32 int     c_active_windings;
33 int     c_peak_windings;
34 int     c_winding_allocs;
35 int     c_winding_points;
36
37 #define BOGUS_RANGE     8192
38
39 void pw(winding_t *w)
40 {
41         int             i;
42         for (i=0 ; i<w->numpoints ; i++)
43                 printf ("(%5.1f, %5.1f, %5.1f)\n",w->p[i][0], w->p[i][1],w->p[i][2]);
44 }
45
46
47 /*
48 =============
49 AllocWinding
50 =============
51 */
52 winding_t       *AllocWinding (int points)
53 {
54         winding_t       *w;
55         int                     s;
56
57         if (numthreads == 1)
58         {
59                 c_winding_allocs++;
60                 c_winding_points += points;
61                 c_active_windings++;
62                 if (c_active_windings > c_peak_windings)
63                         c_peak_windings = c_active_windings;
64         }
65         s = sizeof(vec_t)*3*points + sizeof(int);
66         w = malloc (s);
67         memset (w, 0, s); 
68         return w;
69 }
70
71 void FreeWinding (winding_t *w)
72 {
73         if (*(unsigned *)w == 0xdeaddead)
74                 Error ("FreeWinding: freed a freed winding");
75         *(unsigned *)w = 0xdeaddead;
76
77         if (numthreads == 1)
78                 c_active_windings--;
79         free (w);
80 }
81
82 /*
83 ============
84 RemoveColinearPoints
85 ============
86 */
87 int     c_removed;
88
89 void    RemoveColinearPoints (winding_t *w)
90 {
91         int             i, j, k;
92         vec3_t  v1, v2;
93         int             nump;
94         vec3_t  p[MAX_POINTS_ON_WINDING];
95
96         nump = 0;
97         for (i=0 ; i<w->numpoints ; i++)
98         {
99                 j = (i+1)%w->numpoints;
100                 k = (i+w->numpoints-1)%w->numpoints;
101                 VectorSubtract (w->p[j], w->p[i], v1);
102                 VectorSubtract (w->p[i], w->p[k], v2);
103                 VectorNormalize(v1,v1);
104                 VectorNormalize(v2,v2);
105                 if (DotProduct(v1, v2) < 0.999)
106                 {
107                         VectorCopy (w->p[i], p[nump]);
108                         nump++;
109                 }
110         }
111
112         if (nump == w->numpoints)
113                 return;
114
115         if (numthreads == 1)
116                 c_removed += w->numpoints - nump;
117         w->numpoints = nump;
118         memcpy (w->p, p, nump*sizeof(p[0]));
119 }
120
121 /*
122 ============
123 WindingPlane
124 ============
125 */
126 void WindingPlane (winding_t *w, vec3_t normal, vec_t *dist)
127 {
128         vec3_t  v1, v2;
129
130         VectorSubtract (w->p[1], w->p[0], v1);
131         VectorSubtract (w->p[2], w->p[0], v2);
132         CrossProduct (v2, v1, normal);
133         VectorNormalize (normal, normal);
134         *dist = DotProduct (w->p[0], normal);
135
136 }
137
138 /*
139 =============
140 WindingArea
141 =============
142 */
143 vec_t   WindingArea (winding_t *w)
144 {
145         int             i;
146         vec3_t  d1, d2, cross;
147         vec_t   total;
148
149         total = 0;
150         for (i=2 ; i<w->numpoints ; i++)
151         {
152                 VectorSubtract (w->p[i-1], w->p[0], d1);
153                 VectorSubtract (w->p[i], w->p[0], d2);
154                 CrossProduct (d1, d2, cross);
155                 total += 0.5 * VectorLength ( cross );
156         }
157         return total;
158 }
159
160 void    WindingBounds (winding_t *w, vec3_t mins, vec3_t maxs)
161 {
162         vec_t   v;
163         int             i,j;
164
165         mins[0] = mins[1] = mins[2] = 99999;
166         maxs[0] = maxs[1] = maxs[2] = -99999;
167
168         for (i=0 ; i<w->numpoints ; i++)
169         {
170                 for (j=0 ; j<3 ; j++)
171                 {
172                         v = w->p[i][j];
173                         if (v < mins[j])
174                                 mins[j] = v;
175                         if (v > maxs[j])
176                                 maxs[j] = v;
177                 }
178         }
179 }
180
181 /*
182 =============
183 WindingCenter
184 =============
185 */
186 void    WindingCenter (winding_t *w, vec3_t center)
187 {
188         int             i;
189         float   scale;
190
191         VectorCopy (vec3_origin, center);
192         for (i=0 ; i<w->numpoints ; i++)
193                 VectorAdd (w->p[i], center, center);
194
195         scale = 1.0/w->numpoints;
196         VectorScale (center, scale, center);
197 }
198
199 /*
200 =================
201 BaseWindingForPlane
202 =================
203 */
204 winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
205 {
206         int             i, x;
207         vec_t   max, v;
208         vec3_t  org, vright, vup;
209         winding_t       *w;
210         
211 // find the major axis
212
213         max = -BOGUS_RANGE;
214         x = -1;
215         for (i=0 ; i<3; i++)
216         {
217                 v = fabs(normal[i]);
218                 if (v > max)
219                 {
220                         x = i;
221                         max = v;
222                 }
223         }
224         if (x==-1)
225                 Error ("BaseWindingForPlane: no axis found");
226                 
227         VectorCopy (vec3_origin, vup);  
228         switch (x)
229         {
230         case 0:
231         case 1:
232                 vup[2] = 1;
233                 break;          
234         case 2:
235                 vup[0] = 1;
236                 break;          
237         }
238
239         v = DotProduct (vup, normal);
240         VectorMA (vup, -v, normal, vup);
241         VectorNormalize (vup, vup);
242                 
243         VectorScale (normal, dist, org);
244         
245         CrossProduct (vup, normal, vright);
246         
247         VectorScale (vup, 8192, vup);
248         VectorScale (vright, 8192, vright);
249
250 // project a really big axis aligned box onto the plane
251         w = AllocWinding (4);
252         
253         VectorSubtract (org, vright, w->p[0]);
254         VectorAdd (w->p[0], vup, w->p[0]);
255         
256         VectorAdd (org, vright, w->p[1]);
257         VectorAdd (w->p[1], vup, w->p[1]);
258         
259         VectorAdd (org, vright, w->p[2]);
260         VectorSubtract (w->p[2], vup, w->p[2]);
261         
262         VectorSubtract (org, vright, w->p[3]);
263         VectorSubtract (w->p[3], vup, w->p[3]);
264         
265         w->numpoints = 4;
266         
267         return w;       
268 }
269
270 /*
271 ==================
272 CopyWinding
273 ==================
274 */
275 winding_t       *CopyWinding (winding_t *w)
276 {
277         int                     size;
278         winding_t       *c;
279
280         c = AllocWinding (w->numpoints);
281         size = (int)((winding_t *)0)->p[w->numpoints];
282         memcpy (c, w, size);
283         return c;
284 }
285
286 /*
287 ==================
288 ReverseWinding
289 ==================
290 */
291 winding_t       *ReverseWinding (winding_t *w)
292 {
293         int                     i;
294         winding_t       *c;
295
296         c = AllocWinding (w->numpoints);
297         for (i=0 ; i<w->numpoints ; i++)
298         {
299                 VectorCopy (w->p[w->numpoints-1-i], c->p[i]);
300         }
301         c->numpoints = w->numpoints;
302         return c;
303 }
304
305
306 /*
307 =============
308 ClipWindingEpsilon
309 =============
310 */
311 void    ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist, 
312                                 vec_t epsilon, winding_t **front, winding_t **back)
313 {
314         vec_t   dists[MAX_POINTS_ON_WINDING+4];
315         int             sides[MAX_POINTS_ON_WINDING+4];
316         int             counts[3];
317         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
318         int             i, j;
319         vec_t   *p1, *p2;
320         vec3_t  mid;
321         winding_t       *f, *b;
322         int             maxpts;
323         
324         counts[0] = counts[1] = counts[2] = 0;
325
326 // determine sides for each point
327         for (i=0 ; i<in->numpoints ; i++)
328         {
329                 dot = DotProduct (in->p[i], normal);
330                 dot -= dist;
331                 dists[i] = dot;
332                 if (dot > epsilon)
333                         sides[i] = SIDE_FRONT;
334                 else if (dot < -epsilon)
335                         sides[i] = SIDE_BACK;
336                 else
337                 {
338                         sides[i] = SIDE_ON;
339                 }
340                 counts[sides[i]]++;
341         }
342         sides[i] = sides[0];
343         dists[i] = dists[0];
344         
345         *front = *back = NULL;
346
347         if (!counts[0])
348         {
349                 *back = CopyWinding (in);
350                 return;
351         }
352         if (!counts[1])
353         {
354                 *front = CopyWinding (in);
355                 return;
356         }
357
358         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
359                                                                 // of fp grouping errors
360
361         *front = f = AllocWinding (maxpts);
362         *back = b = AllocWinding (maxpts);
363                 
364         for (i=0 ; i<in->numpoints ; i++)
365         {
366                 p1 = in->p[i];
367                 
368                 if (sides[i] == SIDE_ON)
369                 {
370                         VectorCopy (p1, f->p[f->numpoints]);
371                         f->numpoints++;
372                         VectorCopy (p1, b->p[b->numpoints]);
373                         b->numpoints++;
374                         continue;
375                 }
376         
377                 if (sides[i] == SIDE_FRONT)
378                 {
379                         VectorCopy (p1, f->p[f->numpoints]);
380                         f->numpoints++;
381                 }
382                 if (sides[i] == SIDE_BACK)
383                 {
384                         VectorCopy (p1, b->p[b->numpoints]);
385                         b->numpoints++;
386                 }
387
388                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
389                         continue;
390                         
391         // generate a split point
392                 p2 = in->p[(i+1)%in->numpoints];
393                 
394                 dot = dists[i] / (dists[i]-dists[i+1]);
395                 for (j=0 ; j<3 ; j++)
396                 {       // avoid round off error when possible
397                         if (normal[j] == 1)
398                                 mid[j] = dist;
399                         else if (normal[j] == -1)
400                                 mid[j] = -dist;
401                         else
402                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
403                 }
404                         
405                 VectorCopy (mid, f->p[f->numpoints]);
406                 f->numpoints++;
407                 VectorCopy (mid, b->p[b->numpoints]);
408                 b->numpoints++;
409         }
410         
411         if (f->numpoints > maxpts || b->numpoints > maxpts)
412                 Error ("ClipWinding: points exceeded estimate");
413         if (f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING)
414                 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
415 }
416
417
418 /*
419 =============
420 ChopWindingInPlace
421 =============
422 */
423 void ChopWindingInPlace (winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon)
424 {
425         winding_t       *in;
426         vec_t   dists[MAX_POINTS_ON_WINDING+4];
427         int             sides[MAX_POINTS_ON_WINDING+4];
428         int             counts[3];
429         static  vec_t   dot;            // VC 4.2 optimizer bug if not static
430         int             i, j;
431         vec_t   *p1, *p2;
432         vec3_t  mid;
433         winding_t       *f;
434         int             maxpts;
435
436         in = *inout;
437         counts[0] = counts[1] = counts[2] = 0;
438
439 // determine sides for each point
440         for (i=0 ; i<in->numpoints ; i++)
441         {
442                 dot = DotProduct (in->p[i], normal);
443                 dot -= dist;
444                 dists[i] = dot;
445                 if (dot > epsilon)
446                         sides[i] = SIDE_FRONT;
447                 else if (dot < -epsilon)
448                         sides[i] = SIDE_BACK;
449                 else
450                 {
451                         sides[i] = SIDE_ON;
452                 }
453                 counts[sides[i]]++;
454         }
455         sides[i] = sides[0];
456         dists[i] = dists[0];
457         
458         if (!counts[0])
459         {
460                 FreeWinding (in);
461                 *inout = NULL;
462                 return;
463         }
464         if (!counts[1])
465                 return;         // inout stays the same
466
467         maxpts = in->numpoints+4;       // cant use counts[0]+2 because
468                                                                 // of fp grouping errors
469
470         f = AllocWinding (maxpts);
471                 
472         for (i=0 ; i<in->numpoints ; i++)
473         {
474                 p1 = in->p[i];
475                 
476                 if (sides[i] == SIDE_ON)
477                 {
478                         VectorCopy (p1, f->p[f->numpoints]);
479                         f->numpoints++;
480                         continue;
481                 }
482         
483                 if (sides[i] == SIDE_FRONT)
484                 {
485                         VectorCopy (p1, f->p[f->numpoints]);
486                         f->numpoints++;
487                 }
488
489                 if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
490                         continue;
491                         
492         // generate a split point
493                 p2 = in->p[(i+1)%in->numpoints];
494                 
495                 dot = dists[i] / (dists[i]-dists[i+1]);
496                 for (j=0 ; j<3 ; j++)
497                 {       // avoid round off error when possible
498                         if (normal[j] == 1)
499                                 mid[j] = dist;
500                         else if (normal[j] == -1)
501                                 mid[j] = -dist;
502                         else
503                                 mid[j] = p1[j] + dot*(p2[j]-p1[j]);
504                 }
505                         
506                 VectorCopy (mid, f->p[f->numpoints]);
507                 f->numpoints++;
508         }
509         
510         if (f->numpoints > maxpts)
511                 Error ("ClipWinding: points exceeded estimate");
512         if (f->numpoints > MAX_POINTS_ON_WINDING)
513                 Error ("ClipWinding: MAX_POINTS_ON_WINDING");
514
515         FreeWinding (in);
516         *inout = f;
517 }
518
519
520 /*
521 =================
522 ChopWinding
523
524 Returns the fragment of in that is on the front side
525 of the cliping plane.  The original is freed.
526 =================
527 */
528 winding_t       *ChopWinding (winding_t *in, vec3_t normal, vec_t dist)
529 {
530         winding_t       *f, *b;
531
532         ClipWindingEpsilon (in, normal, dist, ON_EPSILON, &f, &b);
533         FreeWinding (in);
534         if (b)
535                 FreeWinding (b);
536         return f;
537 }
538
539
540 /*
541 =================
542 CheckWinding
543
544 =================
545 */
546 void CheckWinding (winding_t *w)
547 {
548         int             i, j;
549         vec_t   *p1, *p2;
550         vec_t   d, edgedist;
551         vec3_t  dir, edgenormal, facenormal;
552         vec_t   area;
553         vec_t   facedist;
554
555         if (w->numpoints < 3)
556                 Error ("CheckWinding: %i points",w->numpoints);
557         
558         area = WindingArea(w);
559         if (area < 1)
560                 Error ("CheckWinding: %f area", area);
561
562         WindingPlane (w, facenormal, &facedist);
563         
564         for (i=0 ; i<w->numpoints ; i++)
565         {
566                 p1 = w->p[i];
567
568                 for (j=0 ; j<3 ; j++)
569                         if (p1[j] > BOGUS_RANGE || p1[j] < -BOGUS_RANGE)
570                                 Error ("CheckFace: BUGUS_RANGE: %f",p1[j]);
571
572                 j = i+1 == w->numpoints ? 0 : i+1;
573                 
574         // check the point is on the face plane
575                 d = DotProduct (p1, facenormal) - facedist;
576                 if (d < -ON_EPSILON || d > ON_EPSILON)
577                         Error ("CheckWinding: point off plane");
578         
579         // check the edge isnt degenerate
580                 p2 = w->p[j];
581                 VectorSubtract (p2, p1, dir);
582                 
583                 if (VectorLength (dir) < ON_EPSILON)
584                         Error ("CheckWinding: degenerate edge");
585                         
586                 CrossProduct (facenormal, dir, edgenormal);
587                 VectorNormalize (edgenormal, edgenormal);
588                 edgedist = DotProduct (p1, edgenormal);
589                 edgedist += ON_EPSILON;
590                 
591         // all other points must be on front side
592                 for (j=0 ; j<w->numpoints ; j++)
593                 {
594                         if (j == i)
595                                 continue;
596                         d = DotProduct (w->p[j], edgenormal);
597                         if (d > edgedist)
598                                 Error ("CheckWinding: non-convex");
599                 }
600         }
601 }
602
603
604 /*
605 ============
606 WindingOnPlaneSide
607 ============
608 */
609 int             WindingOnPlaneSide (winding_t *w, vec3_t normal, vec_t dist)
610 {
611         qboolean        front, back;
612         int                     i;
613         vec_t           d;
614
615         front = false;
616         back = false;
617         for (i=0 ; i<w->numpoints ; i++)
618         {
619                 d = DotProduct (w->p[i], normal) - dist;
620                 if (d < -ON_EPSILON)
621                 {
622                         if (front)
623                                 return SIDE_CROSS;
624                         back = true;
625                         continue;
626                 }
627                 if (d > ON_EPSILON)
628                 {
629                         if (back)
630                                 return SIDE_CROSS;
631                         front = true;
632                         continue;
633                 }
634         }
635
636         if (back)
637                 return SIDE_BACK;
638         if (front)
639                 return SIDE_FRONT;
640         return SIDE_ON;
641 }
642