The top image contains the banding. The lower one does not.
http://home.comcast.net/~exspacea/banding_artifacts.png

The solution follows.

It uses precisely one, (I repeat for emphasis ONE) trig function call.  It
is tan() and it is contained in randomcone()  (Not included here).  I have
no idea how many trig functions Axis_Rotate() uses.  No doubt it is
probably something like 4 or 8 trig function calls.

C_Vector Vector_Machine::jittervector( C_Vector *N, double steradians )
{

 //Given some vector, N, deviate it in a random angular direction
 //within a solid cone that is 'steradians' wide.  Return the new
 //vector.   Despite the name of the variable I chose here,
 //steradians should be in DEGREES, not radians.  (..let the groaning
 //commence..)


 C_Vector bU, bV, bW, projIy;
 C_Vector Iy;

 // Get a random vector that deviates around the positive Y axis.
 Iy = randomcone(steradians);

 // Create an axis system where N is the "height" or "Up" direction
 getorthobasisvectors(N, &bU, &bV, &bW );

 // Project our random vector "as if" N were the Y axis.
 projIy = changeaxissystem( &Iy, &bU, &bV, &bW );

 return projIy;
}


void Vector_Machine::getorthobasisvectors( C_Vector *N,C_Vector *U,
C_Vector *V, C_Vector *W )
{
 /*
  Given some vector in space, N, derive an axis system of
  unit orthogonal basis vectors such that V=N.

   U is considered the "right" direction.
   V is considered the "up" direction.
   W is considered the "into" direction.
 */

 C_Vector rightV, intoV, uN;
 C_Vector eastpole, westpole;
 double singularity;

 uN = (*N);
// We will assume that the right/U direction is always X+ ie <1,0,0>
//  This assumption breaks down at two singularities, namely the east and
west
//  poles of the unit hemisphere.
//  For the east singularity, assume right/U is Z+ (<0,0,1>) and into/W is
Y+
// (<0,1,0>)
//  For the west singularity, assume right/U is Y+ and into/W is Z+

 eastpole.x=1.0;  eastpole.y=0; eastpole.z=0;
 westpole.x=-1.0;  westpole.y=0; westpole.z=0;

 singularity = uN * eastpole;
 if( singularity > 0.0 ) {
  singularity = 1.0 - singularity;
  if( singularity < EPSILON ) {
   rightV.x = 0; rightV.y = 0; rightV.z = 1; // Assumed right.
   goto sendreturn;
  }
 }

 singularity = uN * westpole;
 if( singularity > 0.0 ) {
  singularity = 1.0 - singularity;
  if( singularity < EPSILON ) {
   rightV.x = 0; rightV.y =1; rightV.z = 0; // Assumed right.
   goto sendreturn;
  }
 }

 //  ...otherwise, everyone else gets an X+.
 rightV.x = 1; rightV.y =0; rightV.z = 0;

sendreturn:
 intoV = rightV % uN;  // Give me an orthogonal vector. I dont carewhere.
 intoV.Unitize();
 rightV = uN % intoV;  // The 'true' right.
 rightV.Unitize();

 (*U) = rightV;
 (*V) = uN;
 (*W) = intoV;
}


C_Vector Vector_Machine::changeaxissystem(C_Vector *E, C_Vector *U, C_Vector
*V, C_Vector *W )
{
/*
Given some unit axis system U-V-W, project the point E into this
system.  Return this new point.
*/

 C_Vector P;

 P.x = (U->x)*(E->x)  + (V->x)*(E->y) + (W->x)*(E->z);
 P.y = (U->y)*(E->x)  + (V->y)*(E->y) + (W->y)*(E->z);
 P.z = (U->z)*(E->x)  + (V->z)*(E->y) + (W->z)*(E->z);

 return P;
}

Post a reply to this message