gtkradiant/tools/quake3/common/polylib.c
rambetter c5449b8cc5 Rewriting BaseWindingForPlane() in polylib.c from the ground up. The behavior
is totally unchanged, and I verified this by running extensive tests.  The
only difference is that the math precision is much much better now.
Performance should be better as well (but that is not tested).
This is a major milestone because it fixes two regression tests:
disappearing_sliver2 and sparkly_seam.
Improvements to math precision is ongoing and more improvements can probably
be made even after this patch.
I will update the README.txt files in the regression tests in a separate
commit.  This commit only includes the actual fixed code.


git-svn-id: svn://svn.icculus.org/gtkradiant/GtkRadiant/trunk@371 8a3a26a2-13c4-0310-b231-cf6edde360e5
2010-12-28 10:02:11 +00:00

825 lines
17 KiB
C

/*
Copyright (C) 1999-2007 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 "cmdlib.h"
#include "mathlib.h"
#include "inout.h"
#include "polylib.h"
#include "qfiles.h"
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;
#define BOGUS_RANGE WORLD_SIZE
void pw(winding_t *w)
{
int i;
for (i=0 ; i<w->numpoints ; i++)
Sys_Printf ("(%5.1f, %5.1f, %5.1f)\n",w->p[i][0], w->p[i][1],w->p[i][2]);
}
/*
=============
AllocWinding
=============
*/
winding_t *AllocWinding (int points)
{
winding_t *w;
int s;
if (points >= MAX_POINTS_ON_WINDING)
Error ("AllocWinding failed: MAX_POINTS_ON_WINDING exceeded");
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;
}
s = sizeof(vec_t)*3*points + sizeof(int);
w = safe_malloc (s);
memset (w, 0, s);
return w;
}
void FreeWinding (winding_t *w)
{
if (*(unsigned *)w == 0xdeaddead)
Error ("FreeWinding: freed a freed winding");
*(unsigned *)w = 0xdeaddead;
if (numthreads == 1)
c_active_windings--;
free (w);
}
/*
============
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,v1);
VectorNormalize(v2,v2);
if (DotProduct(v1, v2) < 0.999)
{
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;
VectorSubtract (w->p[1], w->p[0], v1);
VectorSubtract (w->p[2], w->p[0], v2);
CrossProduct (v2, v1, normal);
VectorNormalize (normal, normal);
*dist = DotProduct (w->p[0], normal);
}
/*
=============
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)
{
// The goal in this function is to replicate the exact behavior that was in the original
// BaseWindingForPlane() function (see below). The only thing we're going to change is the
// accuracy of the operation. The original code gave a preference for the vup vector to start
// out as (0, 0, 1), unless the normal had a dominant Z value, in which case vup started out
// as (1, 0, 0). After that, vup was "bent" [along the plane defined by normal and vup] to
// become perpendicular to normal. After that the vright vector was computed as the cross
// product of vup and normal.
// Once these vectors are calculated, I'm constructing the winding points in exactly the same
// way as was done in the original function. Orientation is the same.
// Note that the 4 points in the returned winding_t may actually not be necessary (3 might
// be enough). However, I want to minimize the chance of ANY bugs popping up due to any
// change in behavior of this function. Therefore, behavior stays exactly the same, except
// for precision of math. Performance might be better in the new function as well.
int x, i;
vec_t max, v;
vec3_t vright, vup, org;
winding_t *w;
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");
switch (x) {
case 0: // Fall through to next case.
case 1:
vright[0] = -normal[1];
vright[1] = normal[0];
vright[2] = 0;
break;
case 2:
vright[0] = 0;
vright[1] = -normal[2];
vright[2] = normal[1];
break;
}
CrossProduct(normal, vright, vup);
// IMPORTANT NOTE: vright and vup are NOT unit vectors at this point.
// However, normal, vup, and vright are pairwise perpendicular.
VectorSetLength(vup, MAX_WORLD_COORD * 2, vup);
VectorSetLength(vright, MAX_WORLD_COORD * 2, vright);
VectorScale(normal, dist, org);
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;
}
// Old function, not used but here for reference. Please do not modify it.
// (You may remove it at some point.)
winding_t *_BaseWindingForPlane_orig_(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, vup);
VectorScale (normal, dist, org);
CrossProduct (vup, normal, vright);
// LordHavoc: this has to use *2 because otherwise some created points may
// be inside the world (think of a diagonal case), and any brush with such
// points should be removed, failure to detect such cases is disasterous
VectorScale (vup, MAX_WORLD_COORD*2, vup);
VectorScale (vright, MAX_WORLD_COORD*2, 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)((size_t)((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];
static 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];
static 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] > MAX_WORLD_COORD || p1[j] < MIN_WORLD_COORD)
Error ("CheckFace: MAX_WORLD_COORD exceeded: %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, 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 = qfalse;
back = qfalse;
for (i=0 ; i<w->numpoints ; i++)
{
d = DotProduct (w->p[i], normal) - dist;
if (d < -ON_EPSILON)
{
if (front)
return SIDE_CROSS;
back = qtrue;
continue;
}
if (d > ON_EPSILON)
{
if (back)
return SIDE_CROSS;
front = qtrue;
continue;
}
}
if (back)
return SIDE_BACK;
if (front)
return SIDE_FRONT;
return SIDE_ON;
}
/*
=================
AddWindingToConvexHull
Both w and *hull are on the same plane
=================
*/
#define MAX_HULL_POINTS 128
void AddWindingToConvexHull( winding_t *w, winding_t **hull, vec3_t normal ) {
int i, j, k;
float *p, *copy;
vec3_t dir;
float d;
int numHullPoints, numNew;
vec3_t hullPoints[MAX_HULL_POINTS];
vec3_t newHullPoints[MAX_HULL_POINTS];
vec3_t hullDirs[MAX_HULL_POINTS];
qboolean hullSide[MAX_HULL_POINTS];
qboolean outside;
if ( !*hull ) {
*hull = CopyWinding( w );
return;
}
numHullPoints = (*hull)->numpoints;
memcpy( hullPoints, (*hull)->p, numHullPoints * sizeof(vec3_t) );
for ( i = 0 ; i < w->numpoints ; i++ ) {
p = w->p[i];
// calculate hull side vectors
for ( j = 0 ; j < numHullPoints ; j++ ) {
k = ( j + 1 ) % numHullPoints;
VectorSubtract( hullPoints[k], hullPoints[j], dir );
VectorNormalize( dir, dir );
CrossProduct( normal, dir, hullDirs[j] );
}
outside = qfalse;
for ( j = 0 ; j < numHullPoints ; j++ ) {
VectorSubtract( p, hullPoints[j], dir );
d = DotProduct( dir, hullDirs[j] );
if ( d >= ON_EPSILON ) {
outside = qtrue;
}
if ( d >= -ON_EPSILON ) {
hullSide[j] = qtrue;
} else {
hullSide[j] = qfalse;
}
}
// if the point is effectively inside, do nothing
if ( !outside ) {
continue;
}
// find the back side to front side transition
for ( j = 0 ; j < numHullPoints ; j++ ) {
if ( !hullSide[ j % numHullPoints ] && hullSide[ (j + 1) % numHullPoints ] ) {
break;
}
}
if ( j == numHullPoints ) {
continue;
}
// insert the point here
VectorCopy( p, newHullPoints[0] );
numNew = 1;
// copy over all points that aren't double fronts
j = (j+1)%numHullPoints;
for ( k = 0 ; k < numHullPoints ; k++ ) {
if ( hullSide[ (j+k) % numHullPoints ] && hullSide[ (j+k+1) % numHullPoints ] ) {
continue;
}
copy = hullPoints[ (j+k+1) % numHullPoints ];
VectorCopy( copy, newHullPoints[numNew] );
numNew++;
}
numHullPoints = numNew;
memcpy( hullPoints, newHullPoints, numHullPoints * sizeof(vec3_t) );
}
FreeWinding( *hull );
w = AllocWinding( numHullPoints );
w->numpoints = numHullPoints;
*hull = w;
memcpy( w->p, hullPoints, numHullPoints * sizeof(vec3_t) );
}