// ===================================================================
// $Id: vector3.cpp,v 1.3 2002/12/13 10:34:15 havran Exp $
//
// vector.cpp
//     CVector3D support routines.
//
// Licence and usage conditions in November 2003:
// This file can be used and modified only with written permission of
// Vlastimil Havran, Pohodli 27, 57001 Litomysl, the Czech Republic.
// This library is granted to be used by Polina Kondratieva within
// her MPI HIWI job working on BRDFLIB library. The usage condition
// may change in the next version of the library. At the moment,
// the library development is supported by Max-Planck-Institut for
// Informatik, Stuhlsatzenhausweg 85, 66123 Saarbruecken Germany
// See WWW: http://www.mpi-sb.mpg.de
// REPLACEMENT_STRING
//
// Initial coding by Nicholas Wilt, recoded by Vlastimil Havran, 1998.
// Given to BRDFLIB project for the Polina Kondratieva's
// HIWI JOB assignment in November 2003.
// ===================================================================

// BRDFLIB headers
#include "config.h"
#include "basic.h"
#include "errorg.h"
#include "vector3.h"

__BEGIN_BRDFLIB_SOURCE

// 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(CVector3D &min, const CVector3D &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(CVector3D &max, const CVector3D &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
CVector3D::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
CVector3D::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
CVector3D::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 CVector3D &DirAt, const CVector3D &Viewer,
            const CVector3D &UpL, CVector3D &ViewV, CVector3D &ViewU,
            CVector3D &ViewN)
{
  CVector3D U, V, N;
  CVector3D 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 _DEBUG
  const float eps = 1e-3;
  if (fabs(Magnitude(ViewU) - 1.0) > eps) {
    FATAL << "ViewU magnitude error= " << Magnitude(ViewU) << "\n";
    FATAL_ABORT;
  }
  if (fabs(Magnitude(ViewV) - 1.0) > eps) {
    FATAL << "ViewU magnitude error= " << Magnitude(ViewV) << "\n";
    FATAL_ABORT;
  }
  if (fabs(Magnitude(ViewN) - 1.0) > eps) {
    FATAL << "ViewU magnitude error= " << Magnitude(ViewN) << "\n";
    FATAL_ABORT;
  }
#endif // _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(CVector3D &U, CVector3D &V, // output
	       const CVector3D &normal, // input
	       const CVector3D &dirIncoming)
{
#ifdef _DEBUG
  float d = Magnitude(normal); 
  if ( (d < 0.99) ||
       (d > 1.01) ) {
    FATAL << " The normal has not unit length = " << d << endl;
    FATAL_ABORT;
  }
  d = Magnitude(dirIncoming); 
  if ( (d < 0.99) ||
       (d > 1.01) ) {
    FATAL << " The incoming dir has not unit length = " << d << endl;
    FATAL_ABORT;
  }
#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 _DEBUG
  d = SqrMagnitude(U);
  if ( (d < 0.99) ||
       (d > 1.01) ) {
    FATAL << "The size of U vector incorrect\n";
    FATAL_EXIT;
  }
#endif
  return;  
}

// same as above, but assumes that 

__END_BRDFLIB_SOURCE

