+/*
+=============
+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, normalAccu;
+
+ // 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 are 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.
+
+ // Here is a further discussion regarding the choice of smallestEpsilonAllowed.
+ // In the 32 float world (we can assume vec_t is that), the "epsilon around 1.0" is
+ // 0.00000011921. In the 64 bit float world (we can assume vec_accu_t is that), the
+ // "epsilon around 1.0" is 0.00000000000000022204. (By the way these two epsilons
+ // are defined as VEC_SMALLEST_EPSILON_AROUND_ONE VEC_ACCU_SMALLEST_EPSILON_AROUND_ONE
+ // respectively.) If you divide the first by the second, you get approximately
+ // 536,885,246. Dividing that number by 200,000 (a typical base winding coordinate)
+ // gives 2684. So in other words, if our smallestEpsilonAllowed was chosen as exactly
+ // VEC_SMALLEST_EPSILON_AROUND_ONE, you would be guaranteed at least 2000 "ticks" in
+ // 64-bit land inside of the epsilon for all numbers we're dealing with.
+
+ static const vec_accu_t smallestEpsilonAllowed = ((vec_accu_t) VEC_SMALLEST_EPSILON_AROUND_ONE) * 0.5;
+ if (crudeEpsilon < smallestEpsilonAllowed) fineEpsilon = smallestEpsilonAllowed;
+ else fineEpsilon = (vec_accu_t) crudeEpsilon;
+
+ in = *inout;
+ counts[0] = counts[1] = counts[2] = 0;
+ VectorCopyRegularToAccu(normal, normalAccu);
+
+ for (i = 0; i < in->numpoints; i++)
+ {
+ dists[i] = DotProductAccu(in->p[i], normalAccu) - 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];
+
+ // I'm wondering if whatever code that handles duplicate planes is robust enough
+ // that we never get a case where two nearly equal planes result in 2 NULL windings
+ // due to the 'if' statement below. TODO: Investigate this.
+ if (!counts[SIDE_FRONT]) {
+ FreeWindingAccu(in);
+ *inout = NULL;
+ return;
+ }
+ if (!counts[SIDE_BACK]) {
+ return; // Winding is unmodified.
+ }
+
+ // NOTE: The least number of points that a winding can have at this point is 2.
+ // In that case, one point is SIDE_FRONT and the other is SIDE_BACK.
+
+ maxpts = counts[SIDE_FRONT] + 2; // We dynamically expand if this is too small.
+ f = AllocWindingAccu(maxpts);
+
+ for (i = 0; i < in->numpoints; i++)
+ {
+ p1 = in->p[i];
+
+ if (sides[i] == SIDE_ON || sides[i] == SIDE_FRONT)
+ {
+ if (f->numpoints >= MAX_POINTS_ON_WINDING)
+ Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
+ if (f->numpoints >= maxpts) // This will probably never happen.
+ {
+ Sys_FPrintf(SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n");
+ f = CopyWindingAccuIncreaseSizeAndFreeOld(f);
+ maxpts++;
+ }
+ VectorCopyAccu(p1, f->p[f->numpoints]);
+ f->numpoints++;
+ if (sides[i] == SIDE_ON) continue;
+ }
+ 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 >= MAX_POINTS_ON_WINDING)
+ Error("ChopWindingInPlaceAccu: MAX_POINTS_ON_WINDING");
+ if (f->numpoints >= maxpts) // This will probably never happen.
+ {
+ Sys_FPrintf(SYS_VRB, "WARNING: estimate on chopped winding size incorrect (no problem)\n");
+ f = CopyWindingAccuIncreaseSizeAndFreeOld(f);
+ maxpts++;
+ }
+ VectorCopyAccu(mid, f->p[f->numpoints]);
+ f->numpoints++;
+ }
+
+ FreeWindingAccu(in);
+ *inout = f;
+}
+