//----------------------------------------------------------------------------
// William Baxter III's Ray Tracer
//
//     Project for Comp 238, Raster Graphics
//     University of North Carolina at Chapel Hill
//     
// $Id:$
//----------------------------------------------------------------------------

#include "wb3Sampling.hpp"
#include <stdlib.h>

namespace wb3Sampling 
{


static const float PI = 3.14159265358979f;
static const float TWOPI = 2*PI;
#define ABS(x) (((x)>0)?(x):(-x))

//----------------------------------------------------------------------------
void ComputePerpendiculars(const Vec3f& L, Vec3f& U, Vec3f& V)
{
  // Construct a pair of vectors U,V orthogonal to L (this is probably 
  // Graham-Schmidt) 

  // Find smallest component, zero it, 
  // negate one of the other two and swap them.  Then normalize.
  // Ex. if x1 is the smallest, assign (x2,y2,z2):=(0,z1,-y1)
  // Clearly now v1 dot v2 = x1*0 + y1*z1 + z1*-y1 = 0;
  U = L;
  float *s1 = &(U.x), *s2 = &(U.y), *s3 = &(U.z);
  float *tmpa;
  if (ABS(*s1)>ABS(*s2))
  {
    tmpa = s1;
    s1 = s2;
    s2 = tmpa;
  }
  // xy min in s1 now.
  if (ABS(*s1)>ABS(*s3)) {
    tmpa = s3;
    s3 = s1;
    s1 = tmpa;
  }
  // xyz min in s1 now
  *s1 = 0;
  float tmp = -*s2;
  *s2 = *s3;
  *s3 = tmp;

  // Voila U should be perpendicular now
  U.Normalize();

  // Check that it really is
  //float f1 = L*U;
  //assert(ABS(f1) < 0.001f);

  // And so it's easy to find a v that is too.
  V = (L ^ U);
  V.Normalize();

  // Check again that we're right
  //float f2 = L*v;
  //float f3 = U*v;
  //assert(ABS(f2) < 0.001f);
  //assert(ABS(f3) < 0.001f);

}

//----------------------------------------------------------------------------
void PerturbSampleRadially(Vec3f &pos, const Vec3f& L, float radius)
{
  Vec3f U, V;
  ComputePerpendiculars(L, U, V);
  PerturbSampleRadially(pos, U, V, radius);
}
//----------------------------------------------------------------------------
void PerturbSampleRadially(
    Vec3f &pos, const Vec3f& U, const Vec3f& V, float radius)
{
  // We want to perturb up to 'radius' units away from pos
  // and in the plane perpendicular defined by u,v

#if 0
  // Uniform sampling in R,Theta
  // PROBLEM: This gives more weight to the middle of the 
  // disk.  We could weight the r distribution linearly towards
  // the outer edge...
  float r = (float(rand())/float(RAND_MAX))*radius;
  float t = (float(rand())/float(RAND_MAX))*TWOPI;
  pos += ( r * (cos(t)*U + sin(t)*V) );
#else
  // Uniform area sampling. Throw darts at a square of side 2*radius;
  // just throw another if you find you missed the inscribed circle.
  float r2 = radius*radius;
  float su,sv;
  do {
    su = ((float(rand())/float(RAND_MAX))*2 - 1.0)*radius;
    sv = ((float(rand())/float(RAND_MAX))*2 - 1.0)*radius;
  } while(su*su + sv*sv > r2);
  pos += su*U + sv*V;
#endif
}

//----------------------------------------------------------------------------





} // end namespace