From 56e74492bcbbfadb24c4433f96c18249cf0d4302 Mon Sep 17 00:00:00 2001 From: rambetter Date: Thu, 30 Dec 2010 08:26:49 +0000 Subject: [PATCH] This is in branch Rambetter-math-fix-experiments. - Adding new function ChopWindingInPlaceAccu() in polylib.c. This function is just like ChopWindingInPlace() except that the winding input has the higher vect_accu_t resolution. I actually did a deep massage of this code as well (went over it with a fine-toothed comb). I still need to examine it at least one more time when I'm fresh, and fix an issue related to choosing the value for maxpts. - In mathlib.h, defining VEC_SMALLEST_EPSILON and VEC_ACCU_SMALLEST_EPSILON constants that will help determine suitable epsilons to use for certain operations. Comments are provided detailing what these constants do. - In mathlib.h, defining DotProductAccu() and VectorCopyAccu(). - Small comment in brush.c that has concerns about a certain use of BaseWindingForPlane() in WriteBSPBrushMap(). State of this code is: compiles on Linux. May not compile on Windows. Not tested (new functionality not used yet in any real code, only functions exist). git-svn-id: svn://svn.icculus.org/gtkradiant/GtkRadiant/branches/Rambetter-math-fix-experiments@382 8a3a26a2-13c4-0310-b231-cf6edde360e5 --- libs/mathlib.h | 20 +++++- tools/quake3/common/polylib.c | 132 ++++++++++++++++++++++++++++++++++ tools/quake3/common/polylib.h | 1 + tools/quake3/q3map2/brush.c | 2 + 4 files changed, 154 insertions(+), 1 deletion(-) diff --git a/libs/mathlib.h b/libs/mathlib.h index 8b474c0a..ae98071e 100644 --- a/libs/mathlib.h +++ b/libs/mathlib.h @@ -24,6 +24,7 @@ Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA // mathlib.h #include +#include #include "bytebool.h" @@ -37,6 +38,13 @@ typedef vec_t vec3_t[3]; typedef vec_t vec5_t[5]; typedef vec_t vec4_t[4]; +// Smallest positive value such that 1.0 + VEC_SMALLEST_EPSILON != 1.0 +// In the case of 32 bits (which is the case), it's 0.00000011921. +// Don't forget that your epsilons should depend on the possible range of values, +// because for example 1000.0 + VEC_SMALLEST_EPSILON will almost certainly be +// equal to 1000.0. +#define VEC_SMALLEST_EPSILON FLT_EPSILON + #define SIDE_FRONT 0 #define SIDE_ON 2 #define SIDE_BACK 1 @@ -308,13 +316,23 @@ vec_t ray_intersect_triangle(const ray_t *ray, qboolean bCullBack, const vec3_t typedef double vec_accu_t; typedef vec_accu_t vec3_accu_t[3]; +// Smallest positive value such that 1.0 + VEC_ACCU_SMALLEST_EPSILON != 1.0 +// In the case of 64 bits (which is the case), it's 0.00000000000000022204. +// Don't forget that your epsilons should depend on the possible range of values, +// because for example 1000.0 + VEC_ACCU_SMALLEST_EPSILON will almost certainly +// be equal to 1000.0. +#define VEC_ACCU_SMALLEST_EPSILON DBL_EPSILON + // 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. +// scaling by a large factor. That is why I'm copying and pasting the above #defines here, +// so that we can break these out into functions if need be. +#define DotProductAccu(x, y) ((x)[0] * (y)[0] + (x)[1] * (y)[1] + (x)[2] * (y)[2]) #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 VectorCopyAccu(a, b) ((b)[0] = (a)[0], (b)[1] = (a)[1], (b)[2] = (a)[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]) diff --git a/tools/quake3/common/polylib.c b/tools/quake3/common/polylib.c index 87bf79fe..93f15241 100644 --- a/tools/quake3/common/polylib.c +++ b/tools/quake3/common/polylib.c @@ -592,6 +592,138 @@ void ClipWindingEpsilon (winding_t *in, vec3_t normal, vec_t dist, } +/* +============= +ChopWindingInPlaceAccu +============= +*/ +void ChopWindingInPlaceAccu(winding_accu_t **inout, vec3_t normal, vec_t dist, vec_t crudeEpsilon) +{ + vec_accu_t fineEpsilon; + winding_accu_t *in; + int counts[3]; + int i, j; + vec_accu_t dists[MAX_POINTS_ON_WINDING + 1]; + int sides[MAX_POINTS_ON_WINDING + 1]; + int maxpts; + winding_accu_t *f; + vec_accu_t *p1, *p2; + vec_accu_t w; + vec3_accu_t mid; + + // We require at least a very small epsilon. It's a good idea for several reasons. + // First, we will be dividing by a potentially very small distance below. We don't + // want that distance to be too small; otherwise, things "blow up" with little accuracy + // due to the division. (After a second look, the value w below is in range (0,1), but + // graininess problem remains.) Second, Having minimum epsilon also prevents the following + // situation. Say for example we have a perfect octagon defined by the input winding. + // Say our chopping plane (defined by normal and dist) is essentially the same plane + // that the octagon is sitting on. Well, due to rounding errors, it may be that point + // 1 of the octagon might be in front, point 2 might be in back, point 3 might be in + // front, point 4 might be in back, and so on. So we could end up with a very ugly- + // looking chopped winding, and this might be undesirable, and would at least lead to + // a possible exhaustion of MAX_POINTS_ON_WINDING. It's better to assume that points + // very very close to the plane ore on the plane, using an infinitesimal epsilon amount. + + // Now, the original ChopWindingInPlace() function used a vec_t-based winding_t. + // So this minimum epsilon is quite similar to casting the higher resolution numbers to + // the lower resolution and comparing them in the lower resolution mode. We explicitly + // choose the minimum epsilon as something around the vec_t epsilon of one because we + // want the resolution of vec_accu_t to have a large resolution around the epsilon. + // Some of that leftover resolution even goes away after we scale to points far away. + + static const vec_accu_t smallestEpsilonAllowed = ((vec_accu_t) VEC_SMALLEST_EPSILON) * 0.5; + if (crudeEpsilon < smallestEpsilonAllowed) fineEpsilon = smallestEpsilonAllowed; + else fineEpsilon = (vec_accu_t) crudeEpsilon; + + in = *inout; + counts[0] = counts[1] = counts[2] = 0; + + for (i = 0; i < in->numpoints; i++) + { + dists[i] = DotProductAccu(in->p[i], normal) - dist; + if (dists[i] > fineEpsilon) sides[i] = SIDE_FRONT; + else if (dists[i] < -fineEpsilon) sides[i] = SIDE_BACK; + else sides[i] = SIDE_ON; + counts[sides[i]]++; + } + sides[i] = sides[0]; + dists[i] = dists[0]; + + if (!counts[SIDE_FRONT]) { + FreeWindingAccu(in); + *inout = NULL; + return; + } + if (!counts[SIDE_BACK]) { + return; // Winding is unmodified. + } + + // TODO: Why the hell are we using this number? The original comment says something along + // the lines of "Can't use counts[0] + 2 because of fp grouping errors", and I'm not + // entirely sure what that means. Is it a problem later on in how the caller of this function + // uses the returned winding afterwards? Is it a problem in that the algorithm below will + // bump into the max in some cases? Anyhow, investigate and fix this number to be a correct + // and predictable number. We can also dynamically grow inside of this function in case we + // do hit the max for some strange reason. + maxpts = in->numpoints + 4; + + f = AllocWindingAccu(maxpts); + + for (i = 0; i < in->numpoints; i++) + { + p1 = in->p[i]; + + if (sides[i] == SIDE_ON) + { + if (f->numpoints >= maxpts) + Error("ChopWindingInPlaceAccu: points exceeded estimate"); + if (f->numpoints >= MAX_POINTS_ON_WINDING) + Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING"); + VectorCopyAccu(p1, f->p[f->numpoints]); + f->numpoints++; + continue; + } + if (sides[i] == SIDE_FRONT) + { + if (f->numpoints >= maxpts) + Error("ChopWindingInPlaceAccu: points exceeded estimate"); + if (f->numpoints >= MAX_POINTS_ON_WINDING) + Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING"); + VectorCopyAccu(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) ? 0 : (i + 1)]; + + // The divisor's absolute value is greater than the dividend's absolute value. + // w is in the range (0,1). + w = dists[i] / (dists[i] - dists[i + 1]); + + for (j = 0; j < 3; j++) + { + // Avoid round-off error when possible. Check axis-aligned normal. + if (normal[j] == 1) mid[j] = dist; + else if (normal[j] == -1) mid[j] = -dist; + else mid[j] = p1[j] + (w * (p2[j] - p1[j])); + } + if (f->numpoints >= maxpts) + Error("ChopWindingInPlaceAccu: points exceeded estimate"); + if (f->numpoints >= MAX_POINTS_ON_WINDING) + Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING"); + VectorCopyAccu(mid, f->p[f->numpoints]); + f->numpoints++; + } + + FreeWindingAccu(in); + *inout = f; +} + /* ============= ChopWindingInPlace diff --git a/tools/quake3/common/polylib.h b/tools/quake3/common/polylib.h index 4eff5c0d..558c8a22 100644 --- a/tools/quake3/common/polylib.h +++ b/tools/quake3/common/polylib.h @@ -69,3 +69,4 @@ typedef struct } winding_accu_t; winding_accu_t *BaseWindingForPlaneAccu(vec3_t normal, vec_t dist); +void ChopWindingInPlaceAccu(winding_accu_t **w, vec3_t normal, vec_t dist, vec_t epsilon); diff --git a/tools/quake3/q3map2/brush.c b/tools/quake3/q3map2/brush.c index 2517ca63..abe3145c 100644 --- a/tools/quake3/q3map2/brush.c +++ b/tools/quake3/q3map2/brush.c @@ -506,6 +506,8 @@ void WriteBSPBrushMap( char *name, brush_t *list ) fprintf (f, "{\n"); for (i=0,s=list->sides ; inumsides ; i++,s++) { + // TODO: See if we can use a smaller winding to prevent resolution loss. + // Is WriteBSPBrushMap() used only to decompile maps? w = BaseWindingForPlane (mapplanes[s->planenum].normal, mapplanes[s->planenum].dist); fprintf (f,"( %i %i %i ) ", (int)w->p[0][0], (int)w->p[0][1], (int)w->p[0][2]);