]> git.xonotic.org Git - xonotic/netradiant.git/blobdiff - radiant/winding.cpp
Merge commit '020d0244e4239b21dc804d630edff926386ea34f' into master-merge
[xonotic/netradiant.git] / radiant / winding.cpp
index 65393ce908dbd31f1ec6dd4fe682726284562b9b..9651af2c2c027dc665fbda6d9f91dfa7ddb721ec 100644 (file)
-/*\r
-Copyright (C) 1999-2007 id Software, Inc. and contributors.\r
-For a list of contributors, see the accompanying CONTRIBUTORS file.\r
-\r
-This file is part of GtkRadiant.\r
-\r
-GtkRadiant is free software; you can redistribute it and/or modify\r
-it under the terms of the GNU General Public License as published by\r
-the Free Software Foundation; either version 2 of the License, or\r
-(at your option) any later version.\r
-\r
-GtkRadiant is distributed in the hope that it will be useful,\r
-but WITHOUT ANY WARRANTY; without even the implied warranty of\r
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\r
-GNU General Public License for more details.\r
-\r
-You should have received a copy of the GNU General Public License\r
-along with GtkRadiant; if not, write to the Free Software\r
-Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\r
-*/\r
-\r
-\r
-\r
-#include "stdafx.h"\r
-#include <assert.h>\r
-#include "winding.h"\r
-\r
-#define BOGUS_RANGE (g_MaxWorldCoord+1)\r
-\r
-/*\r
-=============\r
-Plane_Equal\r
-=============\r
-*/\r
-#define        NORMAL_EPSILON  0.0001\r
-#define        DIST_EPSILON    0.02\r
-\r
-int Plane_Equal(plane_t *a, plane_t *b, int flip)\r
-{\r
-       vec3_t normal;\r
-       float dist;\r
-\r
-       if (flip) {\r
-               normal[0] = - b->normal[0];\r
-               normal[1] = - b->normal[1];\r
-               normal[2] = - b->normal[2];\r
-               dist = - b->dist;\r
-       }\r
-       else {\r
-               normal[0] = b->normal[0];\r
-               normal[1] = b->normal[1];\r
-               normal[2] = b->normal[2];\r
-               dist = b->dist;\r
-       }\r
-       if (\r
-          fabs(a->normal[0] - normal[0]) < NORMAL_EPSILON\r
-       && fabs(a->normal[1] - normal[1]) < NORMAL_EPSILON\r
-       && fabs(a->normal[2] - normal[2]) < NORMAL_EPSILON\r
-       && fabs(a->dist - dist) < DIST_EPSILON )\r
-               return true;\r
-       return false;\r
-}\r
-\r
-/*\r
-============\r
-Plane_FromPoints\r
-============\r
-*/\r
-int Plane_FromPoints(vec3_t p1, vec3_t p2, vec3_t p3, plane_t *plane)\r
-{\r
-       vec3_t v1, v2;\r
-\r
-       VectorSubtract(p2, p1, v1);\r
-       VectorSubtract(p3, p1, v2);\r
-       //CrossProduct(v2, v1, plane->normal);\r
-       CrossProduct(v1, v2, plane->normal);\r
-       if (VectorNormalize(plane->normal, plane->normal) < 0.1) return false;\r
-       plane->dist = DotProduct(p1, plane->normal);\r
-       return true;\r
-}\r
-\r
-/*\r
-=================\r
-Point_Equal\r
-=================\r
-*/\r
-int Point_Equal(vec3_t p1, vec3_t p2, float epsilon)\r
-{\r
-       int i;\r
-\r
-       for (i = 0; i < 3; i++)\r
-       {\r
-               if (fabs(p1[i] - p2[i]) > epsilon) return false;\r
-       }\r
-       return true;\r
-}\r
-\r
-\r
-/*\r
-=================\r
-Winding_BaseForPlane\r
-=================\r
-*/\r
-//#define DBG_WNDG\r
-winding_t *Winding_BaseForPlane (plane_t *p)\r
-{\r
-       int             i, x;\r
-       vec_t   max, v;\r
-       vec3_t  org, vright, vup;\r
-       winding_t       *w;\r
-       \r
-       // find the major axis\r
-#ifdef DBG_WNDG\r
-  Sys_Printf("Winding_BaseForPlane %p\n",p);\r
-#endif\r
-\r
-  max = -BOGUS_RANGE;\r
-       x = -1;\r
-       for (i=0 ; i<3; i++)\r
-       {\r
-               v = fabs(p->normal[i]);\r
-               if (v > max)\r
-               {\r
-                       x = i;\r
-                       max = v;\r
-               }\r
-       }\r
-       if (x==-1)\r
-               Error ("Winding_BaseForPlane: no axis found");\r
-               \r
-       VectorCopy (vec3_origin, vup);  \r
-       switch (x)\r
-       {\r
-       case 0:\r
-       case 1:\r
-               vup[2] = 1;\r
-               break;          \r
-       case 2:\r
-               vup[0] = 1;\r
-               break;          \r
-       }\r
-\r
-\r
-       v = DotProduct (vup, p->normal);\r
-       VectorMA (vup, -v, p->normal, vup);\r
-       VectorNormalize (vup, vup);\r
-               \r
-       VectorScale (p->normal, p->dist, org);\r
-       \r
-       CrossProduct (vup, p->normal, vright);\r
-       \r
-       VectorScale (vup, BOGUS_RANGE, vup);\r
-       VectorScale (vright, BOGUS_RANGE, vright);\r
-\r
-  // project a really big      axis aligned box onto the plane\r
-       w = Winding_Alloc (4);\r
-       \r
-       VectorSubtract (org, vright, w->points[0]);\r
-       VectorAdd (w->points[0], vup, w->points[0]);\r
-       \r
-       VectorAdd (org, vright, w->points[1]);\r
-       VectorAdd (w->points[1], vup, w->points[1]);\r
-       \r
-       VectorAdd (org, vright, w->points[2]);\r
-       VectorSubtract (w->points[2], vup, w->points[2]);\r
-       \r
-       VectorSubtract (org, vright, w->points[3]);\r
-       VectorSubtract (w->points[3], vup, w->points[3]);\r
-       \r
-       w->numpoints = 4;\r
-\r
-       return w;       \r
-}\r
-\r
-// macro to compute winding size\r
-#define WINDING_SIZE(pt) (sizeof(int)*2+sizeof(float)*5*(pt))\r
-\r
-/*\r
-==================\r
-Winding_Alloc\r
-==================\r
-*/\r
-winding_t *Winding_Alloc (int points)\r
-{\r
-       winding_t       *w;\r
-       int                     size;\r
-       \r
-       if (points > MAX_POINTS_ON_WINDING)\r
-               Error ("Winding_Alloc: %i points", points);\r
-       \r
-//     size = (int)((winding_t *)0)->points[points];\r
-  size = WINDING_SIZE(points);\r
-       w = (winding_t*) malloc (size);\r
-       memset (w, 0, size);\r
-       w->maxpoints = points;\r
-       \r
-       return w;\r
-}\r
-\r
-void Winding_Free (winding_t *w)\r
-{\r
-       free(w);\r
-}\r
-\r
-/*\r
-==================\r
-Winding_Clone\r
-==================\r
-*/\r
-winding_t *Winding_Clone(winding_t *w)\r
-{\r
-       int                     size;\r
-       winding_t       *c;\r
-       \r
-//     size = (int)((winding_t *)0)->points[w->numpoints];\r
-  size = WINDING_SIZE(w->numpoints);\r
-       c = (winding_t*)qmalloc (size);\r
-       memcpy (c, w, size);\r
-       return c;\r
-}\r
-\r
-/*\r
-==================\r
-ReverseWinding\r
-==================\r
-*/\r
-winding_t *Winding_Reverse(winding_t *w)\r
-{\r
-       int                     i;\r
-       winding_t       *c;\r
-\r
-       c = Winding_Alloc(w->numpoints);\r
-       for (i = 0; i < w->numpoints; i++)\r
-       {\r
-               VectorCopy (w->points[w->numpoints-1-i], c->points[i]);\r
-       }\r
-       c->numpoints = w->numpoints;\r
-       return c;\r
-}\r
-\r
-/*\r
-==============\r
-Winding_RemovePoint\r
-==============\r
-*/\r
-void Winding_RemovePoint(winding_t *w, int point)\r
-{\r
-       if (point < 0 || point >= w->numpoints)\r
-               Error("Winding_RemovePoint: point out of range");\r
-\r
-       if (point < w->numpoints-1)\r
-       {\r
-               memmove(&w->points[point], &w->points[point+1], (int)((winding_t *)0)->points[w->numpoints - point - 1]);\r
-       }\r
-       w->numpoints--;\r
-}\r
-\r
-/*\r
-=============\r
-Winding_InsertPoint\r
-=============\r
-*/\r
-winding_t *Winding_InsertPoint(winding_t *w, vec3_t point, int spot)\r
-{\r
-       int i, j;\r
-       winding_t *neww;\r
-\r
-       if (spot > w->numpoints)\r
-       {\r
-               Error("Winding_InsertPoint: spot > w->numpoints");\r
-       } //end if\r
-       if (spot < 0)\r
-       {\r
-               Error("Winding_InsertPoint: spot < 0");\r
-       } //end if\r
-       neww = Winding_Alloc(w->numpoints + 1);\r
-       neww->numpoints = w->numpoints + 1;\r
-       for (i = 0, j = 0; i < neww->numpoints; i++)\r
-       {\r
-               if (i == spot)\r
-               {\r
-                       VectorCopy(point, neww->points[i]);\r
-               }\r
-               else\r
-               {\r
-                       VectorCopy(w->points[j], neww->points[i]);\r
-                       j++;\r
-               }\r
-       }\r
-       return neww;\r
-}\r
-\r
-/*\r
-==============\r
-Winding_IsTiny\r
-==============\r
-*/\r
-#define        EDGE_LENGTH     0.2\r
-\r
-int Winding_IsTiny (winding_t *w)\r
-{\r
-       int             i, j;\r
-       vec_t   len;\r
-       vec3_t  delta;\r
-       int             edges;\r
-\r
-       edges = 0;\r
-       for (i=0 ; i<w->numpoints ; i++)\r
-       {\r
-               j = i == w->numpoints - 1 ? 0 : i+1;\r
-               VectorSubtract (w->points[j], w->points[i], delta);\r
-               len = VectorLength (delta);\r
-               if (len > EDGE_LENGTH)\r
-               {\r
-                       if (++edges == 3)\r
-                               return false;\r
-               }\r
-       }\r
-       return true;\r
-}\r
-\r
-/*\r
-==============\r
-Winding_IsHuge\r
-==============\r
-*/\r
-int Winding_IsHuge(winding_t *w)\r
-{\r
-       int             i, j;\r
-\r
-       for (i=0 ; i<w->numpoints ; i++)\r
-       {\r
-               for (j=0 ; j<3 ; j++)\r
-                       if (w->points[i][j] < -BOGUS_RANGE+1 || w->points[i][j] > BOGUS_RANGE-1)\r
-                               return true;\r
-       }\r
-       return false;\r
-}\r
-\r
-/*\r
-=============\r
-Winding_PlanesConcave\r
-=============\r
-*/\r
-#define WCONVEX_EPSILON                0.2\r
-\r
-int Winding_PlanesConcave(winding_t *w1, winding_t *w2,\r
-                                                        vec3_t normal1, vec3_t normal2,\r
-                                                        float dist1, float dist2)\r
-{\r
-       int i;\r
-\r
-       if (!w1 || !w2) return false;\r
-\r
-       // check if one of the points of winding 1 is at the back of the plane of winding 2\r
-       for (i = 0; i < w1->numpoints; i++)\r
-       {\r
-               if (DotProduct(normal2, w1->points[i]) - dist2 > WCONVEX_EPSILON) return true;\r
-       }\r
-       // check if one of the points of winding 2 is at the back of the plane of winding 1\r
-       for (i = 0; i < w2->numpoints; i++)\r
-       {\r
-               if (DotProduct(normal1, w2->points[i]) - dist1 > WCONVEX_EPSILON) return true;\r
-       }\r
-\r
-       return false;\r
-}\r
-\r
-/*\r
-==================\r
-Winding_Clip\r
-\r
-Clips the winding to the plane, returning the new winding on the positive side\r
-Frees the input winding.\r
-If keepon is true, an exactly on-plane winding will be saved, otherwise\r
-it will be clipped away.\r
-==================\r
-*/\r
-winding_t *Winding_Clip (winding_t *in, plane_t *split, qboolean keepon)\r
-{\r
-       vec_t   dists[MAX_POINTS_ON_WINDING];\r
-       int             sides[MAX_POINTS_ON_WINDING];\r
-       int             counts[3];\r
-       vec_t   dot;\r
-       int             i, j;\r
-       vec_t   *p1, *p2;\r
-       vec3_t  mid;\r
-       winding_t       *neww;\r
-       int             maxpts;\r
-       \r
-       counts[0] = counts[1] = counts[2] = 0;\r
-\r
-       // determine sides for each point\r
-       for (i=0 ; i<in->numpoints ; i++)\r
-       {\r
-               dot = DotProduct (in->points[i], split->normal);\r
-               dot -= split->dist;\r
-               dists[i] = dot;\r
-               if (dot > ON_EPSILON)\r
-                       sides[i] = SIDE_FRONT;\r
-               else if (dot < -ON_EPSILON)\r
-                       sides[i] = SIDE_BACK;\r
-               else\r
-               {\r
-                       sides[i] = SIDE_ON;\r
-               }\r
-               counts[sides[i]]++;\r
-       }\r
-       sides[i] = sides[0];\r
-       dists[i] = dists[0];\r
-       \r
-       if (keepon && !counts[0] && !counts[1])\r
-               return in;\r
-               \r
-       if (!counts[0])\r
-       {\r
-               Winding_Free (in);\r
-               return NULL;\r
-       }\r
-       if (!counts[1])\r
-               return in;\r
-       \r
-       maxpts = in->numpoints+4;       // can't use counts[0]+2 because\r
-                                                               // of fp grouping errors\r
-       neww = Winding_Alloc (maxpts);\r
-               \r
-       for (i=0 ; i<in->numpoints ; i++)\r
-       {\r
-               p1 = in->points[i];\r
-               \r
-               if (sides[i] == SIDE_ON)\r
-               {\r
-                       VectorCopy (p1, neww->points[neww->numpoints]);\r
-                       neww->numpoints++;\r
-                       continue;\r
-               }\r
-       \r
-               if (sides[i] == SIDE_FRONT)\r
-               {\r
-                       VectorCopy (p1, neww->points[neww->numpoints]);\r
-                       neww->numpoints++;\r
-               }\r
-               \r
-               if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])\r
-                       continue;\r
-                       \r
-               // generate a split point\r
-               p2 = in->points[(i+1)%in->numpoints];\r
-               \r
-               dot = dists[i] / (dists[i]-dists[i+1]);\r
-               for (j=0 ; j<3 ; j++)\r
-               {       // avoid round off error when possible\r
-                       if (split->normal[j] == 1)\r
-                               mid[j] = split->dist;\r
-                       else if (split->normal[j] == -1)\r
-                               mid[j] = -split->dist;\r
-                       else\r
-                               mid[j] = p1[j] + dot*(p2[j]-p1[j]);\r
-               }\r
-                       \r
-               VectorCopy (mid, neww->points[neww->numpoints]);\r
-               neww->numpoints++;\r
-       }\r
-       \r
-       if (neww->numpoints > maxpts)\r
-               Error ("Winding_Clip: points exceeded estimate");\r
-               \r
-       // free the original winding\r
-       Winding_Free (in);\r
-       \r
-       return neww;\r
-}\r
-\r
-/*\r
-=============\r
-Winding_SplitEpsilon\r
-\r
-  split the input winding with the plane\r
-  the input winding stays untouched\r
-=============\r
-*/\r
-void Winding_SplitEpsilon (winding_t *in, vec3_t normal, double dist, \r
-                               vec_t epsilon, winding_t **front, winding_t **back)\r
-{\r
-       vec_t   dists[MAX_POINTS_ON_WINDING+4];\r
-       int             sides[MAX_POINTS_ON_WINDING+4];\r
-       int             counts[3];\r
-       vec_t   dot;\r
-       int             i, j;\r
-       vec_t   *p1, *p2;\r
-       vec3_t  mid;\r
-       winding_t       *f, *b;\r
-       int             maxpts;\r
-       \r
-       counts[0] = counts[1] = counts[2] = 0;\r
-\r
-       // determine sides for each point\r
-       for (i = 0; i < in->numpoints; i++)\r
-       {\r
-               dot = DotProduct (in->points[i], normal);\r
-               dot -= dist;\r
-               dists[i] = dot;\r
-               if (dot > epsilon)\r
-                       sides[i] = SIDE_FRONT;\r
-               else if (dot < -epsilon)\r
-                       sides[i] = SIDE_BACK;\r
-               else\r
-               {\r
-                       sides[i] = SIDE_ON;\r
-               }\r
-               counts[sides[i]]++;\r
-       }\r
-       sides[i] = sides[0];\r
-       dists[i] = dists[0];\r
-       \r
-       *front = *back = NULL;\r
-\r
-       if (!counts[0])\r
-       {\r
-               *back = Winding_Clone(in);\r
-               return;\r
-       }\r
-       if (!counts[1])\r
-       {\r
-               *front = Winding_Clone(in);\r
-               return;\r
-       }\r
-\r
-       maxpts = in->numpoints+4;       // cant use counts[0]+2 because\r
-                                                               // of fp grouping errors\r
-\r
-       *front = f = Winding_Alloc (maxpts);\r
-       *back = b = Winding_Alloc (maxpts);\r
-               \r
-       for (i = 0; i < in->numpoints; i++)\r
-       {\r
-               p1 = in->points[i];\r
-               \r
-               if (sides[i] == SIDE_ON)\r
-               {\r
-                       VectorCopy (p1, f->points[f->numpoints]);\r
-                       f->numpoints++;\r
-                       VectorCopy (p1, b->points[b->numpoints]);\r
-                       b->numpoints++;\r
-                       continue;\r
-               }\r
-       \r
-               if (sides[i] == SIDE_FRONT)\r
-               {\r
-                       VectorCopy (p1, f->points[f->numpoints]);\r
-                       f->numpoints++;\r
-               }\r
-               if (sides[i] == SIDE_BACK)\r
-               {\r
-                       VectorCopy (p1, b->points[b->numpoints]);\r
-                       b->numpoints++;\r
-               }\r
-\r
-               if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])\r
-                       continue;\r
-                       \r
-               // generate a split point\r
-               p2 = in->points[(i+1)%in->numpoints];\r
-               \r
-               dot = dists[i] / (dists[i]-dists[i+1]);\r
-               for (j = 0; j < 3; j++)\r
-               {\r
-                       // avoid round off error when possible\r
-                       if (normal[j] == 1)\r
-                               mid[j] = dist;\r
-                       else if (normal[j] == -1)\r
-                               mid[j] = -dist;\r
-                       else\r
-                               mid[j] = p1[j] + dot*(p2[j]-p1[j]);\r
-               }\r
-                       \r
-               VectorCopy (mid, f->points[f->numpoints]);\r
-               f->numpoints++;\r
-               VectorCopy (mid, b->points[b->numpoints]);\r
-               b->numpoints++;\r
-       }\r
-       \r
-       if (f->numpoints > maxpts || b->numpoints > maxpts)\r
-               Error ("Winding_Clip: points exceeded estimate");\r
-       if (f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING)\r
-               Error ("Winding_Clip: MAX_POINTS_ON_WINDING");\r
-}\r
-\r
-/*\r
-=============\r
-Winding_TryMerge\r
-\r
-If two windings share a common edge and the edges that meet at the\r
-common points are both inside the other polygons, merge them\r
-\r
-Returns NULL if the windings couldn't be merged, or the new winding.\r
-The originals will NOT be freed.\r
-\r
-if keep is true no points are ever removed\r
-=============\r
-*/\r
-#define        CONTINUOUS_EPSILON      0.005\r
-\r
-winding_t *Winding_TryMerge(winding_t *f1, winding_t *f2, vec3_t planenormal, int keep)\r
-{\r
-       vec_t           *p1, *p2, *p3, *p4, *back;\r
-       winding_t       *newf;\r
-       int                     i, j, k, l;\r
-       vec3_t          normal, delta;\r
-       vec_t           dot;\r
-       qboolean        keep1, keep2;\r
-       \r
-\r
-       //\r
-       // find a common edge\r
-       //      \r
-       p1 = p2 = NULL; // stop compiler warning\r
-       j = 0;                  // \r
-       \r
-       for (i = 0; i < f1->numpoints; i++)\r
-       {\r
-               p1 = f1->points[i];\r
-               p2 = f1->points[(i+1) % f1->numpoints];\r
-               for (j = 0; j < f2->numpoints; j++)\r
-               {\r
-                       p3 = f2->points[j];\r
-                       p4 = f2->points[(j+1) % f2->numpoints];\r
-                       for (k = 0; k < 3; k++)\r
-                       {\r
-                               if (fabs(p1[k] - p4[k]) > 0.1)//EQUAL_EPSILON) //ME\r
-                                       break;\r
-                               if (fabs(p2[k] - p3[k]) > 0.1)//EQUAL_EPSILON) //ME\r
-                                       break;\r
-                       } //end for\r
-                       if (k==3)\r
-                               break;\r
-               } //end for\r
-               if (j < f2->numpoints)\r
-                       break;\r
-       } //end for\r
-       \r
-       if (i == f1->numpoints)\r
-               return NULL;                    // no matching edges\r
-\r
-       //\r
-       // check slope of connected lines\r
-       // if the slopes are colinear, the point can be removed\r
-       //\r
-       back = f1->points[(i+f1->numpoints-1)%f1->numpoints];\r
-       VectorSubtract (p1, back, delta);\r
-       CrossProduct (planenormal, delta, normal);\r
-       VectorNormalize (normal, normal);\r
-       \r
-       back = f2->points[(j+2)%f2->numpoints];\r
-       VectorSubtract (back, p1, delta);\r
-       dot = DotProduct (delta, normal);\r
-       if (dot > CONTINUOUS_EPSILON)\r
-               return NULL;                    // not a convex polygon\r
-       keep1 = (qboolean)(dot < -CONTINUOUS_EPSILON);\r
-       \r
-       back = f1->points[(i+2)%f1->numpoints];\r
-       VectorSubtract (back, p2, delta);\r
-       CrossProduct (planenormal, delta, normal);\r
-       VectorNormalize (normal, normal);\r
-\r
-       back = f2->points[(j+f2->numpoints-1)%f2->numpoints];\r
-       VectorSubtract (back, p2, delta);\r
-       dot = DotProduct (delta, normal);\r
-       if (dot > CONTINUOUS_EPSILON)\r
-               return NULL;                    // not a convex polygon\r
-       keep2 = (qboolean)(dot < -CONTINUOUS_EPSILON);\r
-\r
-       //\r
-       // build the new polygon\r
-       //\r
-       newf = Winding_Alloc (f1->numpoints + f2->numpoints);\r
-       \r
-       // copy first polygon\r
-       for (k=(i+1)%f1->numpoints ; k != i ; k=(k+1)%f1->numpoints)\r
-       {\r
-               if (!keep && k==(i+1)%f1->numpoints && !keep2)\r
-                       continue;\r
-               \r
-               VectorCopy (f1->points[k], newf->points[newf->numpoints]);\r
-               newf->numpoints++;\r
-       }\r
-       \r
-       // copy second polygon\r
-       for (l= (j+1)%f2->numpoints ; l != j ; l=(l+1)%f2->numpoints)\r
-       {\r
-               if (!keep && l==(j+1)%f2->numpoints && !keep1)\r
-                       continue;\r
-               VectorCopy (f2->points[l], newf->points[newf->numpoints]);\r
-               newf->numpoints++;\r
-       }\r
-\r
-       return newf;\r
-}\r
-\r
-/*\r
-============\r
-Winding_Plane\r
-============\r
-*/\r
-void Winding_Plane (winding_t *w, vec3_t normal, double *dist)\r
-{\r
-       vec3_t v1, v2;\r
-       int i;\r
-\r
-       //find two vectors each longer than 0.5 units\r
-       for (i = 0; i < w->numpoints; i++)\r
-       {\r
-               VectorSubtract(w->points[(i+1) % w->numpoints], w->points[i], v1);\r
-               VectorSubtract(w->points[(i+2) % w->numpoints], w->points[i], v2);\r
-               if (VectorLength(v1) > 0.5 && VectorLength(v2) > 0.5) break;\r
-       }\r
-       CrossProduct(v2, v1, normal);\r
-       VectorNormalize(normal, normal);\r
-       *dist = DotProduct(w->points[0], normal);\r
-}\r
-\r
-/*\r
-=============\r
-Winding_Area\r
-=============\r
-*/\r
-float Winding_Area (winding_t *w)\r
-{\r
-       int             i;\r
-       vec3_t  d1, d2, cross;\r
-       float   total;\r
-\r
-       total = 0;\r
-       for (i=2 ; i<w->numpoints ; i++)\r
-       {\r
-               VectorSubtract (w->points[i-1], w->points[0], d1);\r
-               VectorSubtract (w->points[i], w->points[0], d2);\r
-               CrossProduct (d1, d2, cross);\r
-               total += 0.5 * VectorLength ( cross );\r
-       }\r
-       return total;\r
-}\r
-\r
-/*\r
-=============\r
-Winding_Bounds\r
-=============\r
-*/\r
-void Winding_Bounds (winding_t *w, vec3_t mins, vec3_t maxs)\r
-{\r
-       vec_t   v;\r
-       int             i,j;\r
-\r
-       mins[0] = mins[1] = mins[2] = 99999;\r
-       maxs[0] = maxs[1] = maxs[2] = -99999;\r
-\r
-       for (i=0 ; i<w->numpoints ; i++)\r
-       {\r
-               for (j=0 ; j<3 ; j++)\r
-               {\r
-                       v = w->points[i][j];\r
-                       if (v < mins[j])\r
-                               mins[j] = v;\r
-                       if (v > maxs[j])\r
-                               maxs[j] = v;\r
-               }\r
-       }\r
-}\r
-\r
-\r
-/*\r
-=================\r
-Winding_PointInside\r
-=================\r
-*/\r
-int Winding_PointInside(winding_t *w, plane_t *plane, vec3_t point, float epsilon)\r
-{\r
-       int i;\r
-       vec3_t dir, normal, pointvec;\r
-\r
-       for (i = 0; i < w->numpoints; i++)\r
-       {\r
-               VectorSubtract(w->points[(i+1) % w->numpoints], w->points[i], dir);\r
-               VectorSubtract(point, w->points[i], pointvec);\r
-               //\r
-               CrossProduct(dir, plane->normal, normal);\r
-               //\r
-               if (DotProduct(pointvec, normal) < -epsilon) return false;\r
-       }\r
-       return true;\r
-}\r
-\r
-/*\r
-=================\r
-Winding_VectorIntersect\r
-=================\r
-*/\r
-int Winding_VectorIntersect(winding_t *w, plane_t *plane, vec3_t p1, vec3_t p2, float epsilon)\r
-{\r
-       float front, back, frac;\r
-       vec3_t mid;\r
-\r
-       front = DotProduct(p1, plane->normal) - plane->dist;\r
-       back = DotProduct(p2, plane->normal) - plane->dist;\r
-       //if both points at the same side of the plane\r
-       if (front < -epsilon && back < -epsilon) return false;\r
-       if (front > epsilon && back > epsilon) return false;\r
-       //get point of intersection with winding plane\r
-       if (fabs(front-back) < 0.001)\r
-       {\r
-               VectorCopy(p2, mid);\r
-       }\r
-       else\r
-       {\r
-               frac = front/(front-back);\r
-               mid[0] = p1[0] + (p2[0] - p1[0]) * frac;\r
-               mid[1] = p1[1] + (p2[1] - p1[1]) * frac;\r
-               mid[2] = p1[2] + (p2[2] - p1[2]) * frac;\r
-       }\r
-       return Winding_PointInside(w, plane, mid, epsilon);\r
-}\r
-\r
+/*
+   Copyright (C) 1999-2006 Id Software, Inc. and contributors.
+   For a list of contributors, see the accompanying CONTRIBUTORS file.
+
+   This file is part of GtkRadiant.
+
+   GtkRadiant is free software; you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation; either version 2 of the License, or
+   (at your option) any later version.
+
+   GtkRadiant is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with GtkRadiant; if not, write to the Free Software
+   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+#include "winding.h"
+
+#include <algorithm>
+
+#include "math/line.h"
+
+
+inline double plane3_distance_to_point( const Plane3& plane, const DoubleVector3& point ){
+       return vector3_dot( point, plane.normal() ) - plane.dist();
+}
+
+inline double plane3_distance_to_point( const Plane3& plane, const Vector3& point ){
+       return vector3_dot( point, plane.normal() ) - plane.dist();
+}
+
+/// \brief Returns the point at which \p line intersects \p plane, or an undefined value if there is no intersection.
+inline DoubleVector3 line_intersect_plane( const DoubleLine& line, const Plane3& plane ){
+       return line.origin + vector3_scaled(
+                          line.direction,
+                          -plane3_distance_to_point( plane, line.origin )
+                          / vector3_dot( line.direction, plane.normal() )
+                          );
+}
+
+inline bool float_is_largest_absolute( double axis, double other ){
+       return fabs( axis ) > fabs( other );
+}
+
+/// \brief Returns the index of the component of \p v that has the largest absolute value.
+inline int vector3_largest_absolute_component_index( const DoubleVector3& v ){
+       return ( float_is_largest_absolute( v[1], v[0] ) )
+                  ? ( float_is_largest_absolute( v[1], v[2] ) )
+                  ? 1
+                  : 2
+                  : ( float_is_largest_absolute( v[0], v[2] ) )
+                  ? 0
+                  : 2;
+}
+
+/// \brief Returns the infinite line that is the intersection of \p plane and \p other.
+inline DoubleLine plane3_intersect_plane3( const Plane3& plane, const Plane3& other ){
+       DoubleLine line;
+       line.direction = vector3_cross( plane.normal(), other.normal() );
+       switch ( vector3_largest_absolute_component_index( line.direction ) )
+       {
+       case 0:
+               line.origin.x() = 0;
+               line.origin.y() = ( -other.dist() * plane.normal().z() - -plane.dist() * other.normal().z() ) / line.direction.x();
+               line.origin.z() = ( -plane.dist() * other.normal().y() - -other.dist() * plane.normal().y() ) / line.direction.x();
+               break;
+       case 1:
+               line.origin.x() = ( -plane.dist() * other.normal().z() - -other.dist() * plane.normal().z() ) / line.direction.y();
+               line.origin.y() = 0;
+               line.origin.z() = ( -other.dist() * plane.normal().x() - -plane.dist() * other.normal().x() ) / line.direction.y();
+               break;
+       case 2:
+               line.origin.x() = ( -other.dist() * plane.normal().y() - -plane.dist() * other.normal().y() ) / line.direction.z();
+               line.origin.y() = ( -plane.dist() * other.normal().x() - -other.dist() * plane.normal().x() ) / line.direction.z();
+               line.origin.z() = 0;
+               break;
+       default:
+               break;
+       }
+
+       return line;
+}
+
+
+/// \brief Keep the value of \p infinity as small as possible to improve precision in Winding_Clip.
+void Winding_createInfinite( FixedWinding& winding, const Plane3& plane, double infinity ){
+       double max = -infinity;
+       int x = -1;
+       for ( int i = 0 ; i < 3; i++ )
+       {
+               double d = fabs( plane.normal()[i] );
+               if ( d > max ) {
+                       x = i;
+                       max = d;
+               }
+       }
+       if ( x == -1 ) {
+               globalErrorStream() << "invalid plane\n";
+               return;
+       }
+
+       DoubleVector3 vup = g_vector3_identity;
+       switch ( x )
+       {
+       case 0:
+       case 1:
+               vup[2] = 1;
+               break;
+       case 2:
+               vup[0] = 1;
+               break;
+       }
+
+
+       vector3_add( vup, vector3_scaled( plane.normal(), -vector3_dot( vup, plane.normal() ) ) );
+       vector3_normalise( vup );
+
+       DoubleVector3 org = vector3_scaled( plane.normal(), plane.dist() );
+
+       DoubleVector3 vright = vector3_cross( vup, plane.normal() );
+
+       vector3_scale( vup, infinity );
+       vector3_scale( vright, infinity );
+
+       // project a really big  axis aligned box onto the plane
+
+       DoubleLine r1, r2, r3, r4;
+       r1.origin = vector3_added( vector3_subtracted( org, vright ), vup );
+       r1.direction = vector3_normalised( vright );
+       winding.push_back( FixedWindingVertex( r1.origin, r1, c_brush_maxFaces ) );
+       r2.origin = vector3_added( vector3_added( org, vright ), vup );
+       r2.direction = vector3_normalised( vector3_negated( vup ) );
+       winding.push_back( FixedWindingVertex( r2.origin, r2, c_brush_maxFaces ) );
+       r3.origin = vector3_subtracted( vector3_added( org, vright ), vup );
+       r3.direction = vector3_normalised( vector3_negated( vright ) );
+       winding.push_back( FixedWindingVertex( r3.origin, r3, c_brush_maxFaces ) );
+       r4.origin = vector3_subtracted( vector3_subtracted( org, vright ), vup );
+       r4.direction = vector3_normalised( vup );
+       winding.push_back( FixedWindingVertex( r4.origin, r4, c_brush_maxFaces ) );
+}
+
+
+inline PlaneClassification Winding_ClassifyDistance( const double distance, const double epsilon ){
+       if ( distance > epsilon ) {
+               return ePlaneFront;
+       }
+       if ( distance < -epsilon ) {
+               return ePlaneBack;
+       }
+       return ePlaneOn;
+}
+
+/// \brief Returns true if
+/// !flipped && winding is completely BACK or ON
+/// or flipped && winding is completely FRONT or ON
+bool Winding_TestPlane( const Winding& winding, const Plane3& plane, bool flipped ){
+       const int test = ( flipped ) ? ePlaneBack : ePlaneFront;
+       for ( Winding::const_iterator i = winding.begin(); i != winding.end(); ++i )
+       {
+               if ( test == Winding_ClassifyDistance( plane3_distance_to_point( plane, ( *i ).vertex ), ON_EPSILON ) ) {
+                       return false;
+               }
+       }
+       return true;
+}
+
+/// \brief Returns true if any point in \p w1 is in front of plane2, or any point in \p w2 is in front of plane1
+bool Winding_PlanesConcave( const Winding& w1, const Winding& w2, const Plane3& plane1, const Plane3& plane2 ){
+       return !Winding_TestPlane( w1, plane2, false ) || !Winding_TestPlane( w2, plane1, false );
+}
+
+brushsplit_t Winding_ClassifyPlane( const Winding& winding, const Plane3& plane ){
+       brushsplit_t split;
+       for ( Winding::const_iterator i = winding.begin(); i != winding.end(); ++i )
+       {
+               ++split.counts[Winding_ClassifyDistance( plane3_distance_to_point( plane, ( *i ).vertex ), ON_EPSILON )];
+       }
+       return split;
+}
+
+
+#define DEBUG_EPSILON ON_EPSILON
+const double DEBUG_EPSILON_SQUARED = DEBUG_EPSILON * DEBUG_EPSILON;
+
+#define WINDING_DEBUG 0
+
+/// \brief Clip \p winding which lies on \p plane by \p clipPlane, resulting in \p clipped.
+/// If \p winding is completely in front of the plane, \p clipped will be identical to \p winding.
+/// If \p winding is completely in back of the plane, \p clipped will be empty.
+/// If \p winding intersects the plane, the edge of \p clipped which lies on \p clipPlane will store the value of \p adjacent.
+void Winding_Clip( const FixedWinding& winding, const Plane3& plane, const Plane3& clipPlane, std::size_t adjacent, FixedWinding& clipped ){
+       PlaneClassification classification = Winding_ClassifyDistance( plane3_distance_to_point( clipPlane, winding.back().vertex ), ON_EPSILON );
+       PlaneClassification nextClassification;
+       // for each edge
+       for ( std::size_t next = 0, i = winding.size() - 1; next != winding.size(); i = next, ++next, classification = nextClassification )
+       {
+               nextClassification = Winding_ClassifyDistance( plane3_distance_to_point( clipPlane, winding[next].vertex ), ON_EPSILON );
+               const FixedWindingVertex& vertex = winding[i];
+
+               // if first vertex of edge is ON
+               if ( classification == ePlaneOn ) {
+                       // append first vertex to output winding
+                       if ( nextClassification == ePlaneBack ) {
+                               // this edge lies on the clip plane
+                               clipped.push_back( FixedWindingVertex( vertex.vertex, plane3_intersect_plane3( plane, clipPlane ), adjacent ) );
+                       }
+                       else
+                       {
+                               clipped.push_back( vertex );
+                       }
+                       continue;
+               }
+
+               // if first vertex of edge is FRONT
+               if ( classification == ePlaneFront ) {
+                       // add first vertex to output winding
+                       clipped.push_back( vertex );
+               }
+               // if second vertex of edge is ON
+               if ( nextClassification == ePlaneOn ) {
+                       continue;
+               }
+               // else if second vertex of edge is same as first
+               else if ( nextClassification == classification ) {
+                       continue;
+               }
+               // else if first vertex of edge is FRONT and there are only two edges
+               else if ( classification == ePlaneFront && winding.size() == 2 ) {
+                       continue;
+               }
+               // else first vertex is FRONT and second is BACK or vice versa
+               else
+               {
+                       // append intersection point of line and plane to output winding
+                       DoubleVector3 mid( line_intersect_plane( vertex.edge, clipPlane ) );
+
+                       if ( classification == ePlaneFront ) {
+                               // this edge lies on the clip plane
+                               clipped.push_back( FixedWindingVertex( mid, plane3_intersect_plane3( plane, clipPlane ), adjacent ) );
+                       }
+                       else
+                       {
+                               clipped.push_back( FixedWindingVertex( mid, vertex.edge, vertex.adjacent ) );
+                       }
+               }
+       }
+}
+
+std::size_t Winding_FindAdjacent( const Winding& winding, std::size_t face ){
+       for ( std::size_t i = 0; i < winding.numpoints; ++i )
+       {
+               ASSERT_MESSAGE( winding[i].adjacent != c_brush_maxFaces, "edge connectivity data is invalid" );
+               if ( winding[i].adjacent == face ) {
+                       return i;
+               }
+       }
+       return c_brush_maxFaces;
+}
+
+std::size_t Winding_Opposite( const Winding& winding, const std::size_t index, const std::size_t other ){
+       ASSERT_MESSAGE( index < winding.numpoints && other < winding.numpoints, "Winding_Opposite: index out of range" );
+
+       double dist_best = 0;
+       std::size_t index_best = c_brush_maxFaces;
+
+       Ray edge( ray_for_points( winding[index].vertex, winding[other].vertex ) );
+
+       for ( std::size_t i = 0; i < winding.numpoints; ++i )
+       {
+               if ( i == index || i == other ) {
+                       continue;
+               }
+
+               double dist_squared = ray_squared_distance_to_point( edge, winding[i].vertex );
+
+               if ( dist_squared > dist_best ) {
+                       dist_best = dist_squared;
+                       index_best = i;
+               }
+       }
+       return index_best;
+}
+
+std::size_t Winding_Opposite( const Winding& winding, const std::size_t index ){
+       return Winding_Opposite( winding, index, Winding_next( winding, index ) );
+}
+
+/// \brief Calculate the \p centroid of the polygon defined by \p winding which lies on plane \p plane.
+void Winding_Centroid( const Winding& winding, const Plane3& plane, Vector3& centroid ){
+       double area2 = 0, x_sum = 0, y_sum = 0;
+       const ProjectionAxis axis = projectionaxis_for_normal( plane.normal() );
+       const indexremap_t remap = indexremap_for_projectionaxis( axis );
+       for ( std::size_t i = winding.numpoints - 1, j = 0; j < winding.numpoints; i = j, ++j )
+       {
+               const double ai = winding[i].vertex[remap.x] * winding[j].vertex[remap.y] - winding[j].vertex[remap.x] * winding[i].vertex[remap.y];
+               area2 += ai;
+               x_sum += ( winding[j].vertex[remap.x] + winding[i].vertex[remap.x] ) * ai;
+               y_sum += ( winding[j].vertex[remap.y] + winding[i].vertex[remap.y] ) * ai;
+       }
+
+       centroid[remap.x] = static_cast<float>( x_sum / ( 3 * area2 ) );
+       centroid[remap.y] = static_cast<float>( y_sum / ( 3 * area2 ) );
+       {
+               Ray ray( Vector3( 0, 0, 0 ), Vector3( 0, 0, 0 ) );
+               ray.origin[remap.x] = centroid[remap.x];
+               ray.origin[remap.y] = centroid[remap.y];
+               ray.direction[remap.z] = 1;
+               centroid[remap.z] = static_cast<float>( ray_distance_to_plane( ray, plane ) );
+       }
+}