ioquake/gtkradiant
0
0
Fork 0
mirror of https://github.com/TTimo/GtkRadiant.git synced 2024-11-10 07:11:54 +00:00
Code Issues Projects Releases Packages Wiki Activity

This is in branch Rambetter-math-fix-experiments.

Browse source 
- Adding new function BaseWindingForPlaneAccu() in polylib.c.  Analogous to
original BaseWindingForPlane() only that the return value is a
winding_accu_t (new type w/ double precision). 	This function is not being
used yet.

- Restoring the original BaseWindingForPlane() function (pre-r371).  Will keep
this as a reference.  What I did exactly:
  * Renamed BaseWindingForPlane() to BaseWindingForPlaneAccu(), changed stuff.
  * Renamed _BaseWindingForPlane_orig_() to BaseWindingForPlane(), preserved.

- Adding things like vec_accu_t, vec3_accu_t, VectorSubtractAccu(),
VectorAddAccu(), VectorScaleAccu(), CrossProductAccu(), and
winding_accu_t.	 Also AllocWindingAccu() and FreeWindingAccu().

- Removing the recently	added VectorSetLength() function because it is no
longer needed.

State of this code is: compiles, not tested (new functionality not used	yet).


git-svn-id: svn://svn.icculus.org/gtkradiant/GtkRadiant/branches/Rambetter-math-fix-experiments@381 8a3a26a2-13c4-0310-b231-cf6edde360e5
This commit is contained in:
rambetter 2010-12-30 01:25:37 +00:00
parent 40dab90db9
commit a28d27640c
4 changed files with 166 additions and 57 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

21
libs/mathlib.h
View file

@ -81,7 +81,6 @@ void VectorMA( const vec3_t va, vec_t scale, const vec3_t vb, vec3_t vc );
void _CrossProduct (vec3_t v1, vec3_t v2, vec3_t cross);
vec_t VectorNormalize (const vec3_t in, vec3_t out);
vec_t VectorSetLength (const vec3_t in, vec_t length, vec3_t out);
vec_t ColorNormalize( const vec3_t in, vec3_t out );
void VectorInverse (vec3_t v);
void VectorPolar(vec3_t v, float radius, float theta, float phi);
@ -299,6 +298,26 @@ vec_t ray_intersect_point(const ray_t *ray, const vec3_t point, vec_t epsilon, v
/*! return true if triangle intersects ray... dist = dist from intersection point to ray-origin */
vec_t ray_intersect_triangle(const ray_t *ray, qboolean bCullBack, const vec3_t vert0, const vec3_t vert1, const vec3_t vert2);
////////////////////////////////////////////////////////////////////////////////////////////
// Below is double-precision math stuff. This was initially needed by the base winding code
// in q3map2 brush processing. These definitions can be used wherever extra precision is an
// absolute must.
////////////////////////////////////////////////////////////////////////////////////////////
typedef double vec_accu_t;
typedef vec_accu_t vec3_accu_t[3];
// TODO: I have a feeling it may be safer to break these function out into actual functions
// in order to avoid accidental loss of precision. For example, say you call
// VectorScaleAccu(vec3_t, vec_t, vec3_accu_t). The scale would take place in 32 bit land
// and the result would be cast to 64 bit, which would cause total loss of precision when
// scaling by a large factor.
#define VectorSubtractAccu(a, b, c) ((c)[0] = (a)[0] - (b)[0], (c)[1] = (a)[1] - (b)[1], (c)[2] = (a)[2] - (b)[2])
#define VectorAddAccu(a, b, c) ((c)[0] = (a)[0] + (b)[0], (c)[1] = (a)[1] + (b)[1], (c)[2] = (a)[2] + (b)[2])
#define VectorScaleAccu(a, b, c) ((c)[0] = (b) * (a)[0], (c)[1] = (b) * (a)[1], (c)[2] = (b) * (a)[2])
#define CrossProductAccu(a, b, c) ((c)[0] = (a)[1] * (b)[2] - (a)[2] * (b)[1], (c)[1] = (a)[2] * (b)[0] - (a)[0] * (b)[2], (c)[2] = (a)[0] * (b)[1] - (a)[1] * (b)[0])
#ifdef __cplusplus
}
#endif

15
libs/mathlib/mathlib.c
View file

@ -143,21 +143,6 @@ vec_t VectorNormalize( const vec3_t in, vec3_t out ) {
return length;
}
vec_t VectorSetLength(const vec3_t in, vec_t length, vec3_t out) {
vec_t origLength;
origLength = (vec_t) sqrt((in[0] * in[0]) + (in[1] * in[1]) + (in[2] * in[2]));
if (origLength == 0)
{
VectorClear(out);
return 0;
}
VectorScale(in, length / origLength, out);
return origLength;
}
vec_t ColorNormalize( const vec3_t in, vec3_t out ) {
float max, scale;

173
tools/quake3/common/polylib.c
View file

@ -73,6 +73,39 @@ winding_t *AllocWinding (int points)
return w;
}
/*
=============
AllocWindingAccu
=============
*/
winding_accu_t *AllocWindingAccu(int points)
{
winding_accu_t *w;
int s;
if (points >= MAX_POINTS_ON_WINDING)
Error("AllocWindingAccu failed: MAX_POINTS_ON_WINDING exceeded");
if (numthreads == 1)
{
// At the time of this writing, these statistics were not used in any way.
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_accu_t) * 3 * points + sizeof(int);
w = safe_malloc(s);
memset(w, 0, s);
return w;
}
/*
=============
FreeWinding
=============
*/
void FreeWinding (winding_t *w)
{
if (*(unsigned *)w == 0xdeaddead)
@ -84,6 +117,22 @@ void FreeWinding (winding_t *w)
free (w);
}
/*
=============
FreeWindingAccu
=============
*/
void FreeWindingAccu(winding_accu_t *w)
{
if (*((unsigned *) w) == 0xdeaddead)
Error("FreeWindingAccu: freed a freed winding");
*((unsigned *) w) = 0xdeaddead;
if (numthreads == 1)
c_active_windings--;
free(w);
}
/*
============
RemoveColinearPoints
@ -203,33 +252,34 @@ void WindingCenter (winding_t *w, vec3_t center)
/*
=================
BaseWindingForPlane
BaseWindingForPlaneAccu
=================
*/
winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
winding_accu_t *BaseWindingForPlaneAccu(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.
// The goal in this function is to replicate the behavior of the original BaseWindingForPlane()
// function (see below) but at the same time increasing accuracy substantially.
// 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.
// 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.
// 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.
// I'm constructing the winding polygon points in a fashion similar to the method used in the
// original function. Orientation is the same. The size of the winding polygon, however, is
// variable in this function (depending on the angle of normal), and is larger (by about a factor
// of 2) than the winding polygon in the original function.
int x, i;
vec_t max, v;
vec3_t vright, vup, org;
winding_t *w;
vec3_accu_t vright, vup, org;
winding_accu_t *w;
// One of the components of normal must have a magnitiude greater than this value,
// otherwise normal is not a unit vector. This is a little bit of inexpensive
// partial error checking we can do.
max = 0.56; // 1 / sqrt(1^2 + 1^2 + 1^2) = 0.577350269
max = -BOGUS_RANGE;
x = -1;
for (i = 0; i < 3; i++) {
v = fabs(normal[i]);
@ -238,52 +288,93 @@ winding_t *BaseWindingForPlane (vec3_t normal, vec_t dist)
max = v;
}
}
if (x == -1) Error("BaseWindingForPlane: no axis found");
if (x == -1) Error("BaseWindingForPlaneAccu: no dominant axis found because normal is too short");
switch (x) {
case 0: // Fall through to next case.
case 1:
vright[0] = -normal[1];
vright[1] = normal[0];
vright[0] = (vec_accu_t) -normal[1];
vright[1] = (vec_accu_t) normal[0];
vright[2] = 0;
break;
case 2:
vright[0] = 0;
vright[1] = -normal[2];
vright[2] = normal[1];
vright[1] = (vec_accu_t) -normal[2];
vright[2] = (vec_accu_t) 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.
// vright and normal are now perpendicular; you can prove this by taking their
// dot product and seeing that it's always exactly 0 (with no error).
VectorSetLength(vup, MAX_WORLD_COORD * 2, vup);
VectorSetLength(vright, MAX_WORLD_COORD * 2, vright);
VectorScale(normal, dist, org);
// NOTE: vright is NOT a unit vector at this point. vright will have length
// not exceeding 1.0. The minimum length that vright can achieve happens when,
// for example, the Z and X components of the normal input vector are equal,
// and when normal's Y component is zero. In that case Z and X of the normal
// vector are both approximately 0.70711. The resulting vright vector in this
// case will have a length of 0.70711.
w = AllocWinding(4);
// We're relying on the fact that MAX_WORLD_COORD is a power of 2 to keep
// our calculation precise and relatively free of floating point error.
// [However, the code will still work fine if that's not the case.]
VectorScaleAccu(vright, ((vec_accu_t) MAX_WORLD_COORD) * 4.0, vright);
VectorSubtract(org, vright, w->p[0]);
VectorAdd(w->p[0], vup, w->p[0]);
// At time time of this writing, MAX_WORLD_COORD was 65536 (2^16). Therefore
// the length of vright at this point is at least 185364. In comparison, a
// corner of the world at location (65536, 65536, 65536) is distance 113512
// away from the origin.
VectorAdd(org, vright, w->p[1]);
VectorAdd(w->p[1], vup, w->p[1]);
CrossProductAccu(normal, vright, vup); // NOTE: normal is NOT a vec3_accu_t.
VectorAdd(org, vright, w->p[2]);
VectorSubtract(w->p[2], vup, w->p[2]);
// vup now has length equal to that of vright.
VectorSubtract(org, vright, w->p[3]);
VectorSubtract(w->p[3], vup, w->p[3]);
// Here, we're casting dist to the more accurate data type just in case
// VectorScaleAccu() is a #define and not a true function, since normal
// has the lower resolution too. We want the multiply operations to take
// place in the higher resolution data types to prevent loss of accuracy,
// especially considering that we may be scaling by a very large amount.
VectorScaleAccu(normal, (vec_accu_t) dist, org);
// org is now a point on the plane defined by normal and dist. Furthermore,
// org, vright, and vup are pairwise perpendicular. Now, the 4 vectors
// { (+-)vright + (+-)vup } have length that is at least sqrt(185364^2 + 185364^2),
// which is about 262144. That length lies outside the world, since the furthest
// point in the world has distance 113512 from the origin as mentioned above.
// Also, these 4 vectors are perpendicular to the org vector. So adding them
// to org will only increase their length. Therefore the 4 points defined below
// all lie outside of the world. Furthermore, it can be easily seen that the
// edges connecting these 4 points (in the winding_accu_t below) lie completely
// outside the world. sqrt(262144^2 + 262144^2)/2 = 185363, which is greater than
// 113512.
w = AllocWindingAccu(4);
VectorSubtractAccu(org, vright, w->p[0]);
VectorAddAccu(w->p[0], vup, w->p[0]);
VectorAddAccu(org, vright, w->p[1]);
VectorAddAccu(w->p[1], vup, w->p[1]);
VectorAddAccu(org, vright, w->p[2]);
VectorSubtractAccu(w->p[2], vup, w->p[2]);
VectorSubtractAccu(org, vright, w->p[3]);
VectorSubtractAccu(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)
/*
=================
BaseWindingForPlane
Original BaseWindingForPlane() function that has serious accuracy problems. Here is why.
TODO: Explain why.
=================
*/
winding_t *BaseWindingForPlane(vec3_t normal, vec_t dist)
{
int i, x;
vec_t max, v;

14
tools/quake3/common/polylib.h
View file

@ -55,3 +55,17 @@ void ChopWindingInPlace (winding_t **w, vec3_t normal, vec_t dist, vec_t epsilon
// frees the original if clipped
void pw(winding_t *w);
///////////////////////////////////////////////////////////////////////////////////////
// Below is double-precision stuff. This was initially needed by the base winding code
// in q3map2 brush processing.
///////////////////////////////////////////////////////////////////////////////////////
typedef struct
{
int numpoints;
vec3_accu_t p[4]; // variable sized
} winding_accu_t;
winding_accu_t *BaseWindingForPlaneAccu(vec3_t normal, vec_t dist);