+ =============
+ 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;
+}
+
+/*
+ =============
+ 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;