source: GTP/trunk/Lib/Vis/Preprocessing/src/Vector3.cpp @ 2575

#include "Matrix4x4.h"
#include "Vector3.h"
#include "Halton.h"
#include <float.h>

namespace GtpVisibilityPreprocessor {

// Given min a vector to minimize and a candidate vector, replace
// elements of min whose corresponding elements in Candidate are
// smaller.  This function is used for finding objects' bounds,
// among other things.
void
Minimize(Vector3 &min, const Vector3 &Candidate)
{
  if (Candidate.x < min.x)
    min.x = Candidate.x;
  if (Candidate.y < min.y)
    min.y = Candidate.y;
  if (Candidate.z < min.z)
    min.z = Candidate.z;
}

// Given min a vector to minimize and a candidate vector, replace
// elements of min whose corresponding elements in Candidate are
// larger.  This function is used for finding objects' bounds,
// among other things.
void
Maximize(Vector3 &max, const Vector3 &Candidate)
{
  if (Candidate.x > max.x)
    max.x = Candidate.x;
  if (Candidate.y > max.y)
    max.y = Candidate.y;
  if (Candidate.z > max.z)
    max.z = Candidate.z;
}

// Project the vector onto the YZ, XZ, or XY plane depending on which.
// which      Coordinate plane to project onto
// 0          YZ
// 1          XZ
// 2          XY
// This function is used by the polygon intersection code.

void
Vector3::ExtractVerts(float *px, float *py, int which) const
{
  switch (which) {
    case 0:
      *px = y;
      *py = z;
      break;
    case 1:
      *px = x;
      *py = z;
      break;
    case 2:
      *px = x;
      *py = y;
      break;
  }
}

// returns the axis, where the vector has the largest value
int
Vector3::DrivingAxis(void) const
{
  int axis = 0;
  float val = fabs(x);

  if (fabs(y) > val) {
    val = fabs(y);
    axis = 1;
  }

  if (fabs(z) > val)
    axis = 2;

  return axis;
}

// returns the axis, where the vector has the smallest value
int
Vector3::TinyAxis(void) const
{
  int axis = 0;
  float val = fabs(x);

  if (fabs(y) < val) {
    val = fabs(y);
    axis = 1;
  }

  if (fabs(z) < val)
    axis = 2;

  return axis;
}


// Construct a view vector ViewN, and the vector ViewU perpendicular
// to ViewN and lying in the plane given by ViewNxUpl
// the last vector of ortogonal system is ViewV, that is
// perpendicular to both ViewN and ViewU. 
// |ViewN| = |ViewU| = |ViewV| = 1
// The ViewN vector pierces the center of the synthesized image
//     ViewU vector goes from the center image rightwards
//     ViewV vector goes from the center image upwards
void
ViewVectors(const Vector3 &DirAt, const Vector3 &Viewer,
            const Vector3 &UpL, Vector3 &ViewV, Vector3 &ViewU,
            Vector3 &ViewN)
{
  Vector3 U, V, N;
  Vector3 Up = Normalize(UpL);

  N = -Normalize(DirAt);

  V = Normalize(Up - DirAt);
  V -= N * DotProd(V, N);
  V = Normalize(V);
  U = CrossProd(V, N);

  ViewU = U; // rightwards
  ViewV = V; // upwards
  ViewN = -N; // forwards
#ifdef GTP_DEBUG
  const float eps = 1e-3f;
  if (fabs(Magnitude(ViewU) - 1.0) > eps) {
    Debug << "ViewU magnitude error= " << Magnitude(ViewU) << "\n";
  }
  if (fabs(Magnitude(ViewV) - 1.0) > eps) {
    Debug << "ViewU magnitude error= " << Magnitude(ViewV) << "\n";
  }
  if (fabs(Magnitude(ViewN) - 1.0) > eps) {
    Debug << "ViewU magnitude error= " << Magnitude(ViewN) << "\n";
  }
#endif // GTP_DEBUG

  return;
}

// Given the intersection point `P', you have available normal `N'
// of unit length. Let us suppose the incoming ray has direction `D'.
// Then we can construct such two vectors `U' and `V' that
// `U',`N', and `D' are coplanar, and `V' is perpendicular
// to the vectors `N','D', and `V'. Then 'N', 'U', and 'V' create
// the orthonormal base in space R3.
void
TangentVectors(Vector3 &U,
               Vector3 &V, // output
               const Vector3 &normal, // input
               const Vector3 &dirIncoming)
{
#ifdef GTP_DEBUG
  float d = Magnitude(normal); 
  if ( (d < 0.99) ||
       (d > 1.01) ) {
    Debug << " The normal has not unit length = " << d << endl;
  }
  d = Magnitude(dirIncoming); 
  if ( (d < 0.99) ||
       (d > 1.01) ) {
    Debug << " The incoming dir has not unit length = " << d << endl;
  }
#endif
  
  V = CrossProd(normal, dirIncoming);

  if (SqrMagnitude(V) < 1e-3) {
    // the normal and dirIncoming are colinear
    // we can/have to generate arbitrary perpendicular vector to normal.
    if (fabs(normal.x) < 0.6)
      V.SetValue(0.0, -normal.z, normal.y);
    else {
      if (fabs(normal.y) < 0.6)
        V.SetValue(-normal.z, 0.0, normal.x);
      else
        V.SetValue(-normal.y, normal.x, 0.0);
    }
  }
  V = Normalize(V);

  U = CrossProd(normal, V);
#ifdef GTP_DEBUG
  d = SqrMagnitude(U);
  if ( (d < 0.99) ||
       (d > 1.01) ) {
    Debug << "The size of U vector incorrect\n";
  }
#endif
  return;  
}

#define USE_HALTON 0


Vector3
UniformRandomVector()
{
  //  float r1 = RandomValue(0.0f, 1.0f);
  //  float r2 = RandomValue(0.0f, 1.0f);
  float r1, r2;

#if USE_HALTON 
  halton2.GetNext(r1, r2);
#else
   r1 = RandomValue(0.0f, 1.0f);
   r2 = RandomValue(0.0f, 1.0f);
#endif

   return UniformRandomVector(r1, r2);
}

Vector3
UniformRandomVector(const float r1, const float r2)
{
   float cosTheta = 1.0f - 2*r1;
   float sinTheta = sqrt(1 - sqr(cosTheta));
   float fi = 2.0f*M_PI*r2;
   
   Vector3 dir(sinTheta*sin(fi),
                           cosTheta,
                           sinTheta*cos(fi));
   
   return dir;
}

Vector3
CosineRandomVector(const Vector3 &normal)
{
  float r1, r2;
  
  r1 = RandomValue(0.0f, 1.0f);
  r2 = RandomValue(0.0f, 1.0f);
  
  return CosineRandomVector(r1, r2, normal);
}

Vector3
CosineRandomVector(const float r1,
                                   const float r2,
                                   const Vector3 &normal)
{
  float theta = 2.0f*M_PI*r2;
  float radius = sqrt(r1);
  float x = radius*sin(theta);
  float z = radius*cos(theta);
  float y = sqrt(1.0f - x*x - z*z);
  
  Vector3 dir(x,
                          y,
                          z);
  
  //  return Normalize(dir);
   
  Matrix4x4 m = RotationVectorsMatrix(
                                                                          normal,
                                                                          Vector3(0,1,0)
                                                                          );
  Matrix4x4 mi = Invert(m);
  m = m*RotationVectorsMatrix(
                                                          Vector3(0,1,0),
                                                          Normalize(dir))*mi;
  
  return TransformNormal(m, normal);
}

Vector3
UniformRandomVector(const Vector3 &normal)
{
  //  float r1 = RandomValue(0.0f, 1.0f);
  //  float r2 = RandomValue(0.0f, 1.0f);
  float r1, r2;
  
#if USE_HALTON 
  halton2.GetNext(r1, r2);
#else
   r1 = RandomValue(0.0f, 1.0f);
   r2 = RandomValue(0.0f, 1.0f);
#endif
  
  
  float cosTheta = 1.0f - r1;
  float sinTheta = sqrt(1 - sqr(cosTheta));
  float fi = 2.0f*M_PI*r2;
  
  Vector3 dir(sinTheta*sin(fi),
                          cosTheta,
                          sinTheta*cos(fi));
  
  
  //  return Normalize(dir);
  
  Matrix4x4 m = RotationVectorsMatrix(
                                                                          normal,
                                                                          Vector3(0,1,0));
  Matrix4x4 mi = Invert(m);
  m = m*RotationVectorsMatrix(
                                                          Vector3(0,1,0),
                                                          Normalize(dir))*mi;
  
  return TransformNormal(m, normal);
  
  //    return TransformNormal(
  //                     RotationVectorsMatrix(
  //                                           Vector3(0,1,0),
  //                                           Normalize(dir)
  //                                           ),
  //                     normal
  //                     );
}


bool Vector3::CheckValidity() const
{
#ifdef _WIN32
  return !(_isnan(x) || _isnan(y) || _isnan(z));
#else
  return !(isnan(x) || isnan(y) || isnan(z));
#endif  
  //return ((x != x) || (y != y) || (z != z));
}

}