l_poly.c

Go to the documentation of this file.

/*
===========================================================================
Copyright (C) 1999-2005 Id Software, Inc.

This file is part of Quake III Arena source code.

Quake III Arena source code 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.

Quake III Arena source code 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 Foobar; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
===========================================================================
*/

#include <malloc.h>
#include "l_cmd.h"
#include "l_math.h"
#include "l_poly.h"
#include "l_log.h"
#include "l_mem.h"

#define BOGUS_RANGE     65535

extern int numthreads;

// counters are only bumped when running single threaded,
// because they are an awefull coherence problem
int c_active_windings;
int c_peak_windings;
int c_winding_allocs;
int c_winding_points;
int c_windingmemory;
int c_peak_windingmemory;

char windingerror[1024];

void pw(winding_t *w)
{
    int     i;
    for (i=0 ; i<w->numpoints ; i++)
        printf ("(%5.3f, %5.3f, %5.3f)\n",w->p[i][0], w->p[i][1],w->p[i][2]);
}


void ResetWindings(void)
{
    c_active_windings = 0;
    c_peak_windings = 0;
    c_winding_allocs = 0;
    c_winding_points = 0;
    c_windingmemory = 0;
    c_peak_windingmemory = 0;

    strcpy(windingerror, "");
} //end of the function ResetWindings
/*
=============
AllocWinding
=============
*/
winding_t *AllocWinding (int points)
{
    winding_t   *w;
    int         s;

    s = sizeof(vec_t)*3*points + sizeof(int);
    w = GetMemory(s);
    memset(w, 0, s);

    if (numthreads == 1)
    {
        c_winding_allocs++;
        c_winding_points += points;
        c_active_windings++;
        if (c_active_windings > c_peak_windings)
            c_peak_windings = c_active_windings;
        c_windingmemory += MemorySize(w);
        if (c_windingmemory > c_peak_windingmemory)
            c_peak_windingmemory = c_windingmemory;
    } //end if
    return w;
} //end of the function AllocWinding

void FreeWinding (winding_t *w)
{
    if (*(unsigned *)w == 0xdeaddead)
        Error ("FreeWinding: freed a freed winding");

    if (numthreads == 1)
    {
        c_active_windings--;
        c_windingmemory -= MemorySize(w);
    } //end if

    *(unsigned *)w = 0xdeaddead;

    FreeMemory(w);
} //end of the function FreeWinding

int WindingMemory(void)
{
    return c_windingmemory;
} //end of the function WindingMemory

int WindingPeakMemory(void)
{
    return c_peak_windingmemory;
} //end of the function WindingPeakMemory

int ActiveWindings(void)
{
    return c_active_windings;
} //end of the function ActiveWindings
/*
============
RemoveColinearPoints
============
*/
int c_removed;

void RemoveColinearPoints (winding_t *w)
{
    int     i, j, k;
    vec3_t  v1, v2;
    int     nump;
    vec3_t  p[MAX_POINTS_ON_WINDING];

    nump = 0;
    for (i=0 ; i<w->numpoints ; i++)
    {
        j = (i+1)%w->numpoints;
        k = (i+w->numpoints-1)%w->numpoints;
        VectorSubtract (w->p[j], w->p[i], v1);
        VectorSubtract (w->p[i], w->p[k], v2);
        VectorNormalize(v1);
        VectorNormalize(v2);
        if (DotProduct(v1, v2) < 0.999)
        {
            if (nump >= MAX_POINTS_ON_WINDING)
                Error("RemoveColinearPoints: MAX_POINTS_ON_WINDING");
            VectorCopy (w->p[i], p[nump]);
            nump++;
        }
    }

    if (nump == w->numpoints)
        return;

    if (numthreads == 1)
        c_removed += w->numpoints - nump;
    w->numpoints = nump;
    memcpy (w->p, p, nump*sizeof(p[0]));
}

/*
============
WindingPlane
============
*/
void WindingPlane (winding_t *w, vec3_t normal, vec_t *dist)
{
    vec3_t v1, v2;
    int i;

    //find two vectors each longer than 0.5 units
    for (i = 0; i < w->numpoints; i++)
    {
        VectorSubtract(w->p[(i+1) % w->numpoints], w->p[i], v1);
        VectorSubtract(w->p[(i+2) % w->numpoints], w->p[i], v2);
        if (VectorLength(v1) > 0.5 && VectorLength(v2) > 0.5) break;
    } //end for
    CrossProduct(v2, v1, normal);
    VectorNormalize(normal);
    *dist = DotProduct(w->p[0], normal);
} //end of the function WindingPlane

/*
=============
WindingArea
=============
*/
vec_t   WindingArea (winding_t *w)
{
    int     i;
    vec3_t  d1, d2, cross;
    vec_t   total;

    total = 0;
    for (i=2 ; i<w->numpoints ; i++)
    {
        VectorSubtract (w->p[i-1], w->p[0], d1);
        VectorSubtract (w->p[i], w->p[0], d2);
        CrossProduct (d1, d2, cross);
        total += 0.5 * VectorLength ( cross );
    }
    return total;
}

void WindingBounds (winding_t *w, vec3_t mins, vec3_t maxs)
{
    vec_t   v;
    int     i,j;

    mins[0] = mins[1] = mins[2] = 99999;
    maxs[0] = maxs[1] = maxs[2] = -99999;

    for (i=0 ; i<w->numpoints ; i++)
    {
        for (j=0 ; j<3 ; j++)
        {
            v = w->p[i][j];
            if (v < mins[j])
                mins[j] = v;
            if (v > maxs[j])
                maxs[j] = v;
        }
    }
}

/*
=============
WindingCenter
=============
*/
void    WindingCenter (winding_t *w, vec3_t center)
{
    int     i;
    float   scale;

    VectorCopy (vec3_origin, center);
    for (i=0 ; i<w->numpoints ; i++)
        VectorAdd (w->p[i], center, center);

    scale = 1.0/w->numpoints;
    VectorScale (center, scale, center);
}

/*
=================
BaseWindingForPlane
=================
*/
winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
{
    int     i, x;
    vec_t   max, v;
    vec3_t  org, vright, vup;
    winding_t   *w;
    
// find the major axis

    max = -BOGUS_RANGE;
    x = -1;
    for (i=0 ; i<3; i++)
    {
        v = fabs(normal[i]);
        if (v > max)
        {
            x = i;
            max = v;
        }
    }
    if (x==-1)
        Error ("BaseWindingForPlane: no axis found");
        
    VectorCopy (vec3_origin, vup);  
    switch (x)
    {
    case 0:
    case 1:
        vup[2] = 1;
        break;      
    case 2:
        vup[0] = 1;
        break;      
    }

    v = DotProduct (vup, normal);
    VectorMA (vup, -v, normal, vup);
    VectorNormalize (vup);
        
    VectorScale (normal, dist, org);
    
    CrossProduct (vup, normal, vright);
    
    VectorScale (vup, BOGUS_RANGE, vup);
    VectorScale (vright, BOGUS_RANGE, vright);

// project a really big axis aligned box onto the plane
    w = AllocWinding (4);
    
    VectorSubtract (org, vright, w->p[0]);
    VectorAdd (w->p[0], vup, w->p[0]);
    
    VectorAdd (org, vright, w->p[1]);
    VectorAdd (w->p[1], vup, w->p[1]);
    
    VectorAdd (org, vright, w->p[2]);
    VectorSubtract (w->p[2], vup, w->p[2]);
    
    VectorSubtract (org, vright, w->p[3]);
    VectorSubtract (w->p[3], vup, w->p[3]);
    
    w->numpoints = 4;
    
    return w;   
}

/*
==================
CopyWinding
==================
*/
winding_t *CopyWinding (winding_t *w)
{
    int         size;
    winding_t   *c;

    c = AllocWinding (w->numpoints);
    size = (int)((winding_t *)0)->p[w->numpoints];
    memcpy (c, w, size);
    return c;
}

/*
==================
ReverseWinding
==================
*/
winding_t   *ReverseWinding (winding_t *w)
{
    int         i;
    winding_t   *c;

    c = AllocWinding (w->numpoints);
    for (i=0 ; i<w->numpoints ; i++)
    {
        VectorCopy (w->p[w->numpoints-1-i], c->p[i]);
    }
    c->numpoints = w->numpoints;
    return c;
}


/*
=============
ClipWindingEpsilon
=============
*/
void ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist, 
                vec_t epsilon, winding_t **front, winding_t **back)
{
    vec_t   dists[MAX_POINTS_ON_WINDING+4];
    int     sides[MAX_POINTS_ON_WINDING+4];
    int     counts[3];
    //MrElusive: DOH can't use statics when unsing multithreading!!!
    vec_t dot;      // VC 4.2 optimizer bug if not static
    int     i, j;
    vec_t   *p1, *p2;
    vec3_t  mid;
    winding_t   *f, *b;
    int     maxpts;
    
    counts[0] = counts[1] = counts[2] = 0;

// determine sides for each point
    for (i=0 ; i<in->numpoints ; i++)
    {
        dot = DotProduct (in->p[i], normal);
        dot -= dist;
        dists[i] = dot;
        if (dot > epsilon)
            sides[i] = SIDE_FRONT;
        else if (dot < -epsilon)
            sides[i] = SIDE_BACK;
        else
        {
            sides[i] = SIDE_ON;
        }
        counts[sides[i]]++;
    }
    sides[i] = sides[0];
    dists[i] = dists[0];
    
    *front = *back = NULL;

    if (!counts[0])
    {
        *back = CopyWinding (in);
        return;
    }
    if (!counts[1])
    {
        *front = CopyWinding (in);
        return;
    }

    maxpts = in->numpoints+4;   // cant use counts[0]+2 because
                                // of fp grouping errors

    *front = f = AllocWinding (maxpts);
    *back = b = AllocWinding (maxpts);
        
    for (i=0 ; i<in->numpoints ; i++)
    {
        p1 = in->p[i];
        
        if (sides[i] == SIDE_ON)
        {
            VectorCopy (p1, f->p[f->numpoints]);
            f->numpoints++;
            VectorCopy (p1, b->p[b->numpoints]);
            b->numpoints++;
            continue;
        }
    
        if (sides[i] == SIDE_FRONT)
        {
            VectorCopy (p1, f->p[f->numpoints]);
            f->numpoints++;
        }
        if (sides[i] == SIDE_BACK)
        {
            VectorCopy (p1, b->p[b->numpoints]);
            b->numpoints++;
        }

        if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
            continue;
            
    // generate a split point
        p2 = in->p[(i+1)%in->numpoints];
        
        dot = dists[i] / (dists[i]-dists[i+1]);
        for (j=0 ; j<3 ; j++)
        {   // avoid round off error when possible
            if (normal[j] == 1)
                mid[j] = dist;
            else if (normal[j] == -1)
                mid[j] = -dist;
            else
                mid[j] = p1[j] + dot*(p2[j]-p1[j]);
        }
            
        VectorCopy (mid, f->p[f->numpoints]);
        f->numpoints++;
        VectorCopy (mid, b->p[b->numpoints]);
        b->numpoints++;
    }
    
    if (f->numpoints > maxpts || b->numpoints > maxpts)
        Error ("ClipWinding: points exceeded estimate");
    if (f->numpoints > MAX_POINTS_ON_WINDING || b->numpoints > MAX_POINTS_ON_WINDING)
        Error ("ClipWinding: MAX_POINTS_ON_WINDING");
}


/*
=============
ChopWindingInPlace
=============
*/
void ChopWindingInPlace (winding_t **inout, vec3_t normal, vec_t dist, vec_t epsilon)
{
    winding_t *in;
    vec_t   dists[MAX_POINTS_ON_WINDING+4];
    int sides[MAX_POINTS_ON_WINDING+4];
    int counts[3];
    //MrElusive: DOH can't use statics when unsing multithreading!!!
    vec_t dot;      // VC 4.2 optimizer bug if not static
    int i, j;
    vec_t *p1, *p2;
    vec3_t mid;
    winding_t *f;
    int maxpts;

    in = *inout;
    counts[0] = counts[1] = counts[2] = 0;

// determine sides for each point
    for (i=0 ; i<in->numpoints ; i++)
    {
        dot = DotProduct (in->p[i], normal);
        dot -= dist;
        dists[i] = dot;
        if (dot > epsilon)
            sides[i] = SIDE_FRONT;
        else if (dot < -epsilon)
            sides[i] = SIDE_BACK;
        else
        {
            sides[i] = SIDE_ON;
        }
        counts[sides[i]]++;
    }
    sides[i] = sides[0];
    dists[i] = dists[0];
    
    if (!counts[0])
    {
        FreeWinding (in);
        *inout = NULL;
        return;
    }
    if (!counts[1])
        return;     // inout stays the same

    maxpts = in->numpoints+4;   // cant use counts[0]+2 because
                                // of fp grouping errors

    f = AllocWinding (maxpts);
        
    for (i=0 ; i<in->numpoints ; i++)
    {
        p1 = in->p[i];
        
        if (sides[i] == SIDE_ON)
        {
            VectorCopy (p1, f->p[f->numpoints]);
            f->numpoints++;
            continue;
        }
    
        if (sides[i] == SIDE_FRONT)
        {
            VectorCopy (p1, f->p[f->numpoints]);
            f->numpoints++;
        }

        if (sides[i+1] == SIDE_ON || sides[i+1] == sides[i])
            continue;
            
    // generate a split point
        p2 = in->p[(i+1)%in->numpoints];
        
        dot = dists[i] / (dists[i]-dists[i+1]);
        for (j=0 ; j<3 ; j++)
        {   // avoid round off error when possible
            if (normal[j] == 1)
                mid[j] = dist;
            else if (normal[j] == -1)
                mid[j] = -dist;
            else
                mid[j] = p1[j] + dot*(p2[j]-p1[j]);
        }
            
        VectorCopy (mid, f->p[f->numpoints]);
        f->numpoints++;
    }
    
    if (f->numpoints > maxpts)
        Error ("ClipWinding: points exceeded estimate");
    if (f->numpoints > MAX_POINTS_ON_WINDING)
        Error ("ClipWinding: MAX_POINTS_ON_WINDING");

    FreeWinding (in);
    *inout = f;
}


/*
=================
ChopWinding

Returns the fragment of in that is on the front side
of the cliping plane.  The original is freed.
=================
*/
winding_t   *ChopWinding (winding_t *in, vec3_t normal, vec_t dist)
{
    winding_t   *f, *b;

    ClipWindingEpsilon (in, normal, dist, ON_EPSILON, &f, &b);
    FreeWinding (in);
    if (b)
        FreeWinding (b);
    return f;
}


/*
=================
CheckWinding

=================
*/
void CheckWinding (winding_t *w)
{
    int     i, j;
    vec_t   *p1, *p2;
    vec_t   d, edgedist;
    vec3_t  dir, edgenormal, facenormal;
    vec_t   area;
    vec_t   facedist;

    if (w->numpoints < 3)
        Error ("CheckWinding: %i points",w->numpoints);
    
    area = WindingArea(w);
    if (area < 1)
        Error ("CheckWinding: %f area", area);

    WindingPlane (w, facenormal, &facedist);
    
    for (i=0 ; i<w->numpoints ; i++)
    {
        p1 = w->p[i];

        for (j=0 ; j<3 ; j++)
            if (p1[j] > BOGUS_RANGE || p1[j] < -BOGUS_RANGE)
                Error ("CheckWinding: BUGUS_RANGE: %f",p1[j]);

        j = i+1 == w->numpoints ? 0 : i+1;
        
    // check the point is on the face plane
        d = DotProduct (p1, facenormal) - facedist;
        if (d < -ON_EPSILON || d > ON_EPSILON)
            Error ("CheckWinding: point off plane");
    
    // check the edge isnt degenerate
        p2 = w->p[j];
        VectorSubtract (p2, p1, dir);
        
        if (VectorLength (dir) < ON_EPSILON)
            Error ("CheckWinding: degenerate edge");
            
        CrossProduct (facenormal, dir, edgenormal);
        VectorNormalize (edgenormal);
        edgedist = DotProduct (p1, edgenormal);
        edgedist += ON_EPSILON;
        
    // all other points must be on front side
        for (j=0 ; j<w->numpoints ; j++)
        {
            if (j == i)
                continue;
            d = DotProduct (w->p[j], edgenormal);
            if (d > edgedist)
                Error ("CheckWinding: non-convex");
        }
    }
}


/*
============
WindingOnPlaneSide
============
*/
int WindingOnPlaneSide (winding_t *w, vec3_t normal, vec_t dist)
{
    qboolean    front, back;
    int         i;
    vec_t       d;

    front = false;
    back = false;
    for (i=0 ; i<w->numpoints ; i++)
    {
        d = DotProduct (w->p[i], normal) - dist;
        if (d < -ON_EPSILON)
        {
            if (front)
                return SIDE_CROSS;
            back = true;
            continue;
        }
        if (d > ON_EPSILON)
        {
            if (back)
                return SIDE_CROSS;
            front = true;
            continue;
        }
    }

    if (back)
        return SIDE_BACK;
    if (front)
        return SIDE_FRONT;
    return SIDE_ON;
}

//#ifdef ME
    #define CONTINUOUS_EPSILON  0.005
//#else
//  #define CONTINUOUS_EPSILON  0.001
//#endif

/*
=============
TryMergeWinding

If two polygons share a common edge and the edges that meet at the
common points are both inside the other polygons, merge them

Returns NULL if the faces couldn't be merged, or the new face.
The originals will NOT be freed.
=============
*/

winding_t *TryMergeWinding (winding_t *f1, winding_t *f2, vec3_t planenormal)
{
    vec_t       *p1, *p2, *p3, *p4, *back;
    winding_t   *newf;
    int         i, j, k, l;
    vec3_t      normal, delta;
    vec_t       dot;
    qboolean    keep1, keep2;
    

    //
    // find a common edge
    //  
    p1 = p2 = NULL; // stop compiler warning
    j = 0;          // 
    
    for (i = 0; i < f1->numpoints; i++)
    {
        p1 = f1->p[i];
        p2 = f1->p[(i+1) % f1->numpoints];
        for (j = 0; j < f2->numpoints; j++)
        {
            p3 = f2->p[j];
            p4 = f2->p[(j+1) % f2->numpoints];
            for (k = 0; k < 3; k++)
            {
                if (fabs(p1[k] - p4[k]) > 0.1)//EQUAL_EPSILON) //ME
                    break;
                if (fabs(p2[k] - p3[k]) > 0.1)//EQUAL_EPSILON) //ME
                    break;
            } //end for
            if (k==3)
                break;
        } //end for
        if (j < f2->numpoints)
            break;
    } //end for
    
    if (i == f1->numpoints)
        return NULL;            // no matching edges

    //
    // check slope of connected lines
    // if the slopes are colinear, the point can be removed
    //
    back = f1->p[(i+f1->numpoints-1)%f1->numpoints];
    VectorSubtract (p1, back, delta);
    CrossProduct (planenormal, delta, normal);
    VectorNormalize (normal);
    
    back = f2->p[(j+2)%f2->numpoints];
    VectorSubtract (back, p1, delta);
    dot = DotProduct (delta, normal);
    if (dot > CONTINUOUS_EPSILON)
        return NULL;            // not a convex polygon
    keep1 = (qboolean)(dot < -CONTINUOUS_EPSILON);
    
    back = f1->p[(i+2)%f1->numpoints];
    VectorSubtract (back, p2, delta);
    CrossProduct (planenormal, delta, normal);
    VectorNormalize (normal);

    back = f2->p[(j+f2->numpoints-1)%f2->numpoints];
    VectorSubtract (back, p2, delta);
    dot = DotProduct (delta, normal);
    if (dot > CONTINUOUS_EPSILON)
        return NULL;            // not a convex polygon
    keep2 = (qboolean)(dot < -CONTINUOUS_EPSILON);

    //
    // build the new polygon
    //
    newf = AllocWinding (f1->numpoints + f2->numpoints);
    
    // copy first polygon
    for (k=(i+1)%f1->numpoints ; k != i ; k=(k+1)%f1->numpoints)
    {
        if (k==(i+1)%f1->numpoints && !keep2)
            continue;
        
        VectorCopy (f1->p[k], newf->p[newf->numpoints]);
        newf->numpoints++;
    }
    
    // copy second polygon
    for (l= (j+1)%f2->numpoints ; l != j ; l=(l+1)%f2->numpoints)
    {
        if (l==(j+1)%f2->numpoints && !keep1)
            continue;
        VectorCopy (f2->p[l], newf->p[newf->numpoints]);
        newf->numpoints++;
    }

    return newf;
}

//#ifdef ME
//===========================================================================
//
// Parameter:               -
// Returns:                 -
// Changes Globals:     -
//===========================================================================
winding_t *MergeWindings(winding_t *w1, winding_t *w2, vec3_t planenormal)
{
    winding_t *neww;
    float dist;
    int i, j, n, found, insertafter;
    int sides[MAX_POINTS_ON_WINDING+4];
    vec3_t newp[MAX_POINTS_ON_WINDING+4];
    int numpoints;
    vec3_t edgevec, sepnormal, v;

    RemoveEqualPoints(w1, 0.2);
    numpoints = w1->numpoints;
    memcpy(newp, w1->p, w1->numpoints * sizeof(vec3_t));
    //
    for (i = 0; i < w2->numpoints; i++)
    {
        VectorCopy(w2->p[i], v);
        for (j = 0; j < numpoints; j++)
        {
            VectorSubtract(newp[(j+1)%numpoints],
                            newp[(j)%numpoints], edgevec);
            CrossProduct(edgevec, planenormal, sepnormal);
            VectorNormalize(sepnormal);
            if (VectorLength(sepnormal) < 0.9)
            {
                //remove the point from the new winding
                for (n = j; n < numpoints-1; n++)
                {
                    VectorCopy(newp[n+1], newp[n]);
                    sides[n] = sides[n+1];
                } //end for
                numpoints--;
                j--;
                Log_Print("MergeWindings: degenerate edge on winding %f %f %f\n", sepnormal[0],
                                                                                sepnormal[1],
                                                                                sepnormal[2]);
                continue;
            } //end if
            dist = DotProduct(newp[(j)%numpoints], sepnormal);
            if (DotProduct(v, sepnormal) - dist < -0.1) sides[j] = SIDE_BACK;
            else sides[j] = SIDE_FRONT;
        } //end for
        //remove all unnecesary points
        for (j = 0; j < numpoints;)
        {
            if (sides[j] == SIDE_BACK
                && sides[(j+1)%numpoints] == SIDE_BACK)
            {
                //remove the point from the new winding
                for (n = (j+1)%numpoints; n < numpoints-1; n++)
                {
                    VectorCopy(newp[n+1], newp[n]);
                    sides[n] = sides[n+1];
                } //end for
                numpoints--;
            } //end if
            else
            {
                j++;
            } //end else
        } //end for
        //
        found = false;
        for (j = 0; j < numpoints; j++)
        {
            if (sides[j] == SIDE_FRONT
                && sides[(j+1)%numpoints] == SIDE_BACK)
            {
                if (found) Log_Print("Warning: MergeWindings: front to back found twice\n");
                found = true;
            } //end if
        } //end for
        //
        for (j = 0; j < numpoints; j++)
        {
            if (sides[j] == SIDE_FRONT
                && sides[(j+1)%numpoints] == SIDE_BACK)
            {
                insertafter = (j+1)%numpoints;
                //insert the new point after j+1
                for (n = numpoints-1; n > insertafter; n--)
                {
                    VectorCopy(newp[n], newp[n+1]);
                } //end for
                numpoints++;
                VectorCopy(v, newp[(insertafter+1)%numpoints]);
                break;
            } //end if
        } //end for
    } //end for
    neww = AllocWinding(numpoints);
    neww->numpoints = numpoints;
    memcpy(neww->p, newp, numpoints * sizeof(vec3_t));
    RemoveColinearPoints(neww);
    return neww;
} //end of the function MergeWindings
//===========================================================================
//
// Parameter:               -
// Returns:                 -
// Changes Globals:     -
//===========================================================================
char *WindingErrorString(void)
{
    return windingerror;
} //end of the function WindingErrorString
//===========================================================================
//
// Parameter:               -
// Returns:                 -
// Changes Globals:     -
//===========================================================================
int WindingError(winding_t *w)
{
    int     i, j;
    vec_t   *p1, *p2;
    vec_t   d, edgedist;
    vec3_t  dir, edgenormal, facenormal;
    vec_t   area;
    vec_t   facedist;

    if (w->numpoints < 3)
    {
        sprintf(windingerror, "winding %i points", w->numpoints);
        return WE_NOTENOUGHPOINTS;
    } //end if
    
    area = WindingArea(w);
    if (area < 1)
    {
        sprintf(windingerror, "winding %f area", area);
        return WE_SMALLAREA;
    } //end if

    WindingPlane (w, facenormal, &facedist);
    
    for (i=0 ; i<w->numpoints ; i++)
    {
        p1 = w->p[i];

        for (j=0 ; j<3 ; j++)
        {
            if (p1[j] > BOGUS_RANGE || p1[j] < -BOGUS_RANGE)
            {
                sprintf(windingerror, "winding point %d BUGUS_RANGE \'%f %f %f\'", j, p1[0], p1[1], p1[2]);
                return WE_POINTBOGUSRANGE;
            } //end if
        } //end for

        j = i+1 == w->numpoints ? 0 : i+1;
        
    // check the point is on the face plane
        d = DotProduct (p1, facenormal) - facedist;
        if (d < -ON_EPSILON || d > ON_EPSILON)
        {
            sprintf(windingerror, "winding point %d off plane", i);
            return WE_POINTOFFPLANE;
        } //end if
    
    // check the edge isnt degenerate
        p2 = w->p[j];
        VectorSubtract (p2, p1, dir);
        
        if (VectorLength (dir) < ON_EPSILON)
        {
            sprintf(windingerror, "winding degenerate edge %d-%d", i, j);
            return WE_DEGENERATEEDGE;
        } //end if
            
        CrossProduct (facenormal, dir, edgenormal);
        VectorNormalize (edgenormal);
        edgedist = DotProduct (p1, edgenormal);
        edgedist += ON_EPSILON;
        
    // all other points must be on front side
        for (j=0 ; j<w->numpoints ; j++)
        {
            if (j == i)
                continue;
            d = DotProduct (w->p[j], edgenormal);
            if (d > edgedist)
            {
                sprintf(windingerror, "winding non-convex");
                return WE_NONCONVEX;
            } //end if
        } //end for
    } //end for
    return WE_NONE;
} //end of the function WindingError
//===========================================================================
//
// Parameter:               -
// Returns:                 -
// Changes Globals:     -
//===========================================================================
void RemoveEqualPoints(winding_t *w, float epsilon)
{
    int i, nump;
    vec3_t v;
    vec3_t p[MAX_POINTS_ON_WINDING];

    VectorCopy(w->p[0], p[0]);
    nump = 1;
    for (i = 1; i < w->numpoints; i++)
    {
        VectorSubtract(w->p[i], p[nump-1], v);
        if (VectorLength(v) > epsilon)
        {
            if (nump >= MAX_POINTS_ON_WINDING)
                Error("RemoveColinearPoints: MAX_POINTS_ON_WINDING");
            VectorCopy (w->p[i], p[nump]);
            nump++;
        } //end if
    } //end for

    if (nump == w->numpoints)
        return;

    w->numpoints = nump;
    memcpy(w->p, p, nump * sizeof(p[0]));
} //end of the function RemoveEqualPoints
//===========================================================================
// adds the given point to a winding at the given spot
// (for instance when spot is zero then the point is added at position zero)
// the original winding is NOT freed
//
// Parameter:               -
// Returns:                 the new winding with the added point
// Changes Globals:     -
//===========================================================================
winding_t *AddWindingPoint(winding_t *w, vec3_t point, int spot)
{
    int i, j;
    winding_t *neww;

    if (spot > w->numpoints)
    {
        Error("AddWindingPoint: num > w->numpoints");
    } //end if
    if (spot < 0)
    {
        Error("AddWindingPoint: num < 0");
    } //end if
    neww = AllocWinding(w->numpoints + 1);
    neww->numpoints = w->numpoints + 1;
    for (i = 0, j = 0; i < neww->numpoints; i++)
    {
        if (i == spot)
        {
            VectorCopy(point, neww->p[i]);
        } //end if
        else
        {
            VectorCopy(w->p[j], neww->p[i]);
            j++;
        } //end else
    } //end for
    return neww;
} //end of the function AddWindingPoint
//===========================================================================
// the position where the new point should be added in the winding is
// stored in *spot
//
// Parameter:               -
// Returns:                 true if the point is on the winding
// Changes Globals:     -
//===========================================================================
#define MELT_ON_EPSILON     0.2

int PointOnWinding(winding_t *w, vec3_t normal, float dist, vec3_t point, int *spot)
{
    int i, j;
    vec3_t v1, v2;
    vec3_t edgenormal, edgevec;
    float edgedist, dot;

    *spot = 0;
    //the point must be on the winding plane
    dot = DotProduct(point, normal) - dist;
    if (dot < -MELT_ON_EPSILON || dot > MELT_ON_EPSILON) return false;
    //
    for (i = 0; i < w->numpoints; i++)
    {
        j = (i+1) % w->numpoints;
        //get a plane orthogonal to the winding plane through the edge
        VectorSubtract(w->p[j], w->p[i], edgevec);
        CrossProduct(normal, edgevec, edgenormal);
        VectorNormalize(edgenormal);
        edgedist = DotProduct(edgenormal, w->p[i]);
        //point must be not too far from the plane
        dot = DotProduct(point, edgenormal) - edgedist;
        if (dot < -MELT_ON_EPSILON || dot > MELT_ON_EPSILON) continue;
        //vector from first point of winding to the point to test
        VectorSubtract(point, w->p[i], v1);
        //vector from second point of winding to the point to test
        VectorSubtract(point, w->p[j], v2);
        //if the length of the vector is not larger than 0.5 units then
        //the point is assumend to be the same as one of the winding points
        if (VectorNormalize(v1) < 0.5) return false;
        if (VectorNormalize(v2) < 0.5) return false;
        //point must be between the two winding points
        //(the two vectors must be directed towards each other, and on the
        //same straight line)
        if (DotProduct(v1, v2) < -0.99)
        {
            *spot = i + 1;
            return true;
        } //end if
    } //end for
    return false;
} //end of the function PointOnWinding
//===========================================================================
//
// Parameter:               -
// Returns:                 -
// Changes Globals:     -
//===========================================================================
int FindPlaneSeperatingWindings(winding_t *w1, winding_t *w2, vec3_t dir,
                                            vec3_t normal, float *dist)
{
    int i, i2, j, j2, n;
    int sides1[3], sides2[3];
    float dist1, dist2, dot, diff;
    vec3_t normal1, normal2;
    vec3_t v1, v2;

    for (i = 0; i < w1->numpoints; i++)
    {
        i2 = (i+1) % w1->numpoints;
        //
        VectorSubtract(w1->p[i2], w1->p[i], v1);
        if (VectorLength(v1) < 0.1)
        {
            //Log_Write("FindPlaneSeperatingWindings: winding1 with degenerate edge\r\n");
            continue;
        } //end if
        CrossProduct(v1, dir, normal1);
        VectorNormalize(normal1);
        dist1 = DotProduct(normal1, w1->p[i]);
        //
        for (j = 0; j < w2->numpoints; j++)
        {
            j2 = (j+1) % w2->numpoints;
            //
            VectorSubtract(w2->p[j2], w2->p[j], v2);
            if (VectorLength(v2) < 0.1)
            {
                //Log_Write("FindPlaneSeperatingWindings: winding2 with degenerate edge\r\n");
                continue;
            } //end if
            CrossProduct(v2, dir, normal2);
            VectorNormalize(normal2);
            dist2 = DotProduct(normal2, w2->p[j]);
            //
            diff = dist1 - dist2;
            if (diff < -0.1 || diff > 0.1)
            {
                dist2 = -dist2;
                VectorNegate(normal2, normal2);
                diff = dist1 - dist2;
                if (diff < -0.1 || diff > 0.1) continue;
            } //end if
            //check if the normal vectors are equal
            for (n = 0; n < 3; n++)
            {
                diff = normal1[n] - normal2[n];
                if (diff < -0.0001 || diff > 0.0001) break;
            } //end for
            if (n != 3) continue;
            //check on which side of the seperating plane the points of
            //the first winding are
            sides1[0] = sides1[1] = sides1[2] = 0;
            for (n = 0; n < w1->numpoints; n++)
            {
                dot = DotProduct(w1->p[n], normal1) - dist1;
                if (dot > 0.1) sides1[0]++;
                else if (dot < -0.1) sides1[1]++;
                else sides1[2]++;
            } //end for
            //check on which side of the seperating plane the points of
            //the second winding are
            sides2[0] = sides2[1] = sides2[2] = 0;
            for (n = 0; n < w2->numpoints; n++)
            {
                //used normal1 and dist1 (they are equal to normal2 and dist2)
                dot = DotProduct(w2->p[n], normal1) - dist1;
                if (dot > 0.1) sides2[0]++;
                else if (dot < -0.1) sides2[1]++;
                else sides2[2]++;
            } //end for
            //if the first winding has points at both sides
            if (sides1[0] && sides1[1])
            {
                Log_Write("FindPlaneSeperatingWindings: winding1 non-convex\r\n");
                continue;
            } //end if
            //if the second winding has points at both sides
            if (sides2[0] && sides2[1])
            {
                Log_Write("FindPlaneSeperatingWindings: winding2 non-convex\r\n");
                continue;
            } //end if
            //
            if ((!sides1[0] && !sides1[1]) || (!sides2[0] && !sides2[1]))
            {
                //don't use one of the winding planes as the seperating plane
                continue;
            } //end if
            //the windings must be at different sides of the seperating plane
            if ((!sides1[0] && !sides2[1]) || (!sides1[1] && !sides2[0]))
            {
                VectorCopy(normal1, normal);
                *dist = dist1;
                return true;
            } //end if
        } //end for
    } //end for
    return false;
} //end of the function FindPlaneSeperatingWindings
//===========================================================================
//
// Parameter:               -
// Returns:                 -
// Changes Globals:     -
//===========================================================================
#define WCONVEX_EPSILON     0.2

int WindingsNonConvex(winding_t *w1, winding_t *w2,
                             vec3_t normal1, vec3_t normal2,
                             float dist1, float dist2)
{
    int i;

    if (!w1 || !w2) return false;

    //check if one of the points of face1 is at the back of the plane of face2
    for (i = 0; i < w1->numpoints; i++)
    {
        if (DotProduct(normal2, w1->p[i]) - dist2 > WCONVEX_EPSILON) return true;
    } //end for
    //check if one of the points of face2 is at the back of the plane of face1
    for (i = 0; i < w2->numpoints; i++)
    {
        if (DotProduct(normal1, w2->p[i]) - dist1 > WCONVEX_EPSILON) return true;
    } //end for

    return false;
} //end of the function WindingsNonConvex
//===========================================================================
//
// Parameter:               -
// Returns:                 -
// Changes Globals:     -
//===========================================================================
/*
#define VERTEX_EPSILON      0.5

qboolean EqualVertexes(vec3_t v1, vec3_t v2)
{
    float diff;

    diff = v1[0] - v2[0];
    if (diff > -VERTEX_EPSILON && diff < VERTEX_EPSILON)
    {
        diff = v1[1] - v2[1];
        if (diff > -VERTEX_EPSILON && diff < VERTEX_EPSILON)
        {
            diff = v1[2] - v2[2];
            if (diff > -VERTEX_EPSILON && diff < VERTEX_EPSILON)
            {
                return true;
            } //end if
        } //end if
    } //end if
    return false;
} //end of the function EqualVertexes

#define CONTINUOUS_EPSILON  0.001

winding_t *AAS_MergeWindings(winding_t *w1, winding_t *w2, vec3_t windingnormal)
{
    int n, i, k;
    vec3_t normal, delta;
    winding_t *winding, *neww;
    float dist, dot;
    int p1, p2;
    int points[2][64];
    int numpoints[2] = {0, 0};
    int newnumpoints;
    int keep[2];

    if (!FindPlaneSeperatingWindings(w1, w2, windingnormal, normal, &dist)) return NULL;

    //for both windings
    for (n = 0; n < 2; n++)
    {
        if (n == 0) winding = w1;
        else winding = w2;
        //get the points of the winding which are on the seperating plane
        for (i = 0; i < winding->numpoints; i++)
        {
            dot = DotProduct(winding->p[i], normal) - dist;
            if (dot > -ON_EPSILON && dot < ON_EPSILON)
            {
                //don't allow more than 64 points on the seperating plane
                if (numpoints[n] >= 64) Error("AAS_MergeWindings: more than 64 points on seperating plane\n");
                points[n][numpoints[n]++] = i;
            } //end if
        } //end for
        //there must be at least two points of each winding on the seperating plane
        if (numpoints[n] < 2) return NULL;
    } //end for

    //if the first point of winding1 (which is on the seperating plane) is unequal
    //to the last point of winding2 (which is on the seperating plane)
    if (!EqualVertexes(w1->p[points[0][0]], w2->p[points[1][numpoints[1]-1]]))
    {
        return NULL;
    } //end if
    //if the last point of winding1 (which is on the seperating plane) is unequal
    //to the first point of winding2 (which is on the seperating plane)
    if (!EqualVertexes(w1->p[points[0][numpoints[0]-1]], w2->p[points[1][0]]))
    {
        return NULL;
    } //end if
    //
    // check slope of connected lines
    // if the slopes are colinear, the point can be removed
    //