// ===================================================================
// $Id: vector3.h,v 1.3 2002/11/25 20:31:12 havran Exp $
//
// vector3.h
//     Header file for CVector3D class in BRDFLIB.
//
// Class: CVector3D
//
// 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 in OORT package.
// Several functions by Vlastimil Havran and Jiri Bittner.
// ArbitraryNormal by Jaroslav Krivanek.
//
// Given for the use to BRDFLIB project for the Polina Kondratieva's
// HIWI JOB assignment in November 2003.
// ===================================================================


#ifndef __VECTOR3_H__
#define __VECTOR3_H__

// BRDFLIB headers
#include "configh.h"

//standard headers
#include <cmath>

__BEGIN_BRDFLIB_HEADER

// Forward-declare some other classes.
class CMatrix4;

// HACK of returning vector components as array fields.
// NOT guarrantied to work with some strange variable allignment !
#define __VECTOR_HACK

class CVector3D
{
public:
  float x, y, z;

  // constructors
  CVector3D() { }
  CVector3D(float X, float Y, float Z) { x = X; y = Y; z = Z; }
  CVector3D(float X) { x = y = z = X; }
  CVector3D(const CVector3D &v) { x = v.x; y = v.y; z = v.z; }

  // Functions to get at the vector components
  float& operator[] (int inx) {
#ifndef __VECTOR_HACK
    if (inx == 0)
      return x;
    else
      if (inx == 1)
	return y;
    else
      return z;
#else
    return (&x)[inx];
#endif
    
  }

  const float& operator[] (int inx) const {
#ifndef __VECTOR_HACK
    if (inx == 0)
      return x;
    else
      if (inx == 1)
	return y;
      else
	return z;
#else
    return *(&x+inx);
#endif
  }
  void ExtractVerts(float *px, float *py, int which) const;
  void SetValue(const float &a, const float &b, const float &c)
  {  x=a; y=b; z=c; }

  void SetValue(const float a) { x = y = z = a; }
  
  // returns the axis, where the vector has the largest value
  int DrivingAxis(void) const;

  // returns the axis, where the vector has the smallest value
  int TinyAxis(void) const;

  inline float MaxComponent(void) const {
     // return (x > y && x > z) ? x : ((y > z) ? y : z);
     return (x > y) ? ( (x > z) ? x : z) : ( (y > z) ? y : z);
   }
 
   inline CVector3D Abs(void) const {
     return CVector3D(fabs(x), fabs(y), fabs(z));
   }

  // normalizes the vector of unit size corresponding to given vector
  inline void Normalize();
  
  /**
     Return an arbitrary normal to `v'.
     In fact it tries v x (0,0,1) an if the result is too small,
     it definitely does v x (0,1,0). It will always work for
     non-degenareted vector and is much faster than to use
     TangentVectors.

     @param v(in) The vector we want to find normal for.
     @return The normal vector to v.
  */
  friend inline CVector3D ArbitraryNormal(const CVector3D &v); 
  
  // the magnitude=size of the vector
  friend inline float Magnitude(const CVector3D &v);

  // the squared magnitude of the vector .. for efficiency in some cases
  friend inline float SqrMagnitude(const CVector3D &v);

  // creates the vector of unit size corresponding to given vector
  friend inline CVector3D Normalize(const CVector3D &A);

  // Rotate a normal vector.
  friend CVector3D PlaneRotate(const CMatrix4 &, const CVector3D &);

  // construct view vectors .. DirAt is the main viewing direction
  // Viewer is the coordinates of viewer location, UpL is the vector.
  friend void ViewVectors(const CVector3D &DirAt, const CVector3D &Viewer,
	                  const CVector3D &UpL, CVector3D &ViewV,
			  CVector3D &ViewU, CVector3D &ViewN );

  // 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.
  friend void TangentVectors(CVector3D &U, CVector3D &V, // output
			     const CVector3D &normal, // input
			     const CVector3D &dirIncoming);
  // Unary operators
  CVector3D operator+ () const;
  CVector3D operator- () const;

  // Assignment operators
  CVector3D& operator+= (const CVector3D &A);
  CVector3D& operator-= (const CVector3D &A);
  CVector3D& operator*= (const CVector3D &A);
  CVector3D& operator*= (float A);
  CVector3D& operator/= (float A);

  // Binary operators
  friend inline CVector3D operator+ (const CVector3D &A, const CVector3D &B);
  friend inline CVector3D operator- (const CVector3D &A, const CVector3D &B);
  friend inline CVector3D operator* (const CVector3D &A, const CVector3D &B);
  friend inline CVector3D operator* (const CVector3D &A, float B);
  friend inline CVector3D operator* (float A, const CVector3D &B);
  friend CVector3D operator* (const CMatrix4 &, const CVector3D &);
  friend inline CVector3D operator/ (const CVector3D &A, const CVector3D &B);

  friend inline int operator< (const CVector3D &A, const CVector3D &B);
  friend inline int operator<= (const CVector3D &A, const CVector3D &B);

  friend inline CVector3D operator/ (const CVector3D &A, float B);
  friend inline int operator== (const CVector3D &A, const CVector3D &B);
  friend inline float DotProd(const CVector3D &A, const CVector3D &B);
  friend inline CVector3D CrossProd (const CVector3D &A, const CVector3D &B);

  friend ostream& operator<< (ostream &s, const CVector3D &A);
  friend istream& operator>> (istream &s, CVector3D &A);
    

  friend void Minimize(CVector3D &min, const CVector3D &candidate);
  friend void Maximize(CVector3D &max, const CVector3D &candidate);
};

inline CVector3D ArbitraryNormal(const CVector3D &N) {
  float dist2 = N.x * N.x + N.y * N.y;
  if (dist2 > 0.0001) {
    float inv_size = 1.0f/sqrt(dist2);
    return CVector3D(N.y * inv_size, -N.x * inv_size, 0.0f); // N x (0,0,1)
  }
  else {
    float inv_size = 1.0f/sqrt(N.z * N.z + N.x * N.x);
    return CVector3D(-N.z * inv_size, 0.0f, N.x * inv_size); // N x (0,1,0)
  }
}


inline float Magnitude(const CVector3D& v) 
{
  return (float)sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
}

inline float SqrMagnitude(const CVector3D& v) 
{
  return v.x * v.x + v.y * v.y + v.z * v.z;
}

inline CVector3D Normalize(const CVector3D& A)
{
  return A * (1.0f/Magnitude(A));
}

inline float DotProd (const CVector3D& A, const CVector3D& B)
{
  return A.x * B.x + A.y * B.y + A.z * B.z;
}

inline CVector3D CVector3D::operator+ () const
{
  return *this;
}

inline CVector3D CVector3D::operator- () const
{
  return CVector3D(-x, -y, -z);
}

inline CVector3D& CVector3D::operator+= (const CVector3D& A)
{
  x += A.x;  y += A.y;  z += A.z;
  return *this;
}

inline CVector3D& CVector3D::operator-= (const CVector3D& A)
{
  x -= A.x;  y -= A.y;  z -= A.z;
  return *this;
}

inline CVector3D& CVector3D::operator*= (float A)
{
  x *= A;  y *= A;  z *= A;
  return *this;
}

inline CVector3D& CVector3D::operator/= (float A)
{
  float a = 1.0/A;
  x *= a;  y *= a;  z *= a;
  return *this;
}

inline CVector3D& CVector3D::operator*= (const CVector3D &A)
{
  x *= A.x;  y *= A.y;  z *= A.z;
  return *this;
}

inline CVector3D operator+ (const CVector3D &A, const CVector3D &B)
{
  return CVector3D(A.x + B.x, A.y + B.y, A.z + B.z);
}

inline CVector3D operator- (const CVector3D &A, const CVector3D &B)
{
  return CVector3D(A.x - B.x, A.y - B.y, A.z - B.z);
}

inline CVector3D operator* (const CVector3D &A, const CVector3D &B)
{
  return CVector3D(A.x * B.x, A.y * B.y, A.z * B.z);
}

inline CVector3D operator* (const CVector3D &A, float B)
{
  return CVector3D(A.x * B, A.y * B, A.z * B);
}

inline CVector3D operator* (float A, const CVector3D &B)
{
  return CVector3D(B.x * A, B.y * A, B.z * A);
}

inline CVector3D operator/ (const CVector3D &A, const CVector3D &B)
{
  return CVector3D(A.x / B.x, A.y / B.y, A.z / B.z);
}

inline CVector3D operator/ (const CVector3D &A, float B)
{
  float b = 1.0 / B;
  return CVector3D(A.x * b, A.y * b, A.z * b);
}

inline int operator< (const CVector3D &A, const CVector3D& B)
{
  return A.x < B.x && A.y < B.y && A.z < B.z;
}

inline int operator<= (const CVector3D& A, const CVector3D& B)
{
  return A.x <= B.x && A.y <= B.y && A.z <= B.z;
}

// Might replace floating-point == with comparisons of
// magnitudes, if needed.
inline int operator== (const CVector3D& A, const CVector3D& B)
{
  return (A.x == B.x) && (A.y == B.y) && (A.z == B.z);
}

inline CVector3D CrossProd (const CVector3D& A, const CVector3D& B)
{
  return CVector3D(A.y * B.z - A.z * B.y,
		   A.z * B.x - A.x * B.z,
		   A.x * B.y - A.y * B.x);
}

inline void CVector3D::Normalize()
{
  float sqrmag = x * x + y * y + z * z;
  if (sqrmag > 0.0)
    (*this) *= 1.0/sqrt(sqrmag);
}

// Overload << operator for C++-style output
inline ostream&
operator<< (ostream &s, const CVector3D &A)
{
  return s << "(" << A.x << ", " << A.y << ", " << A.z << ")";
}

// Overload >> operator for C++-style input
inline istream&
operator>> (istream &s, CVector3D &A)
{
  char a;
  // read "(x, y, z)"
  return s >> a >> A.x >> a >> A.y >> a >> A.z >> a;
}

__END_BRDFLIB_HEADER

#endif // __VECTOR3_H__
