#ifndef _Vector3_h__
#define _Vector3_h__

#include <iostream>
using namespace std;
#include <math.h>
#include "common.h"


namespace GtpVisibilityPreprocessor {

// Forward-declare some other classes.
class Matrix4x4;
class Vector2;


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

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

  // for compatibility with pascal's code
  void  setX(float q) { x=q; }
  void  setY(float q) { y=q; }
  void  setZ(float q) { z=q; }
  float getX() const { return x; }
  float getY() const { return y; }
  float getZ() const { return z; }

  // constructors
  Vector3() { }

  Vector3(float X, float Y, float Z) { x = X; y = Y; z = Z; }
  Vector3(float X) { x = y = z = X; }
  Vector3(const Vector3 &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
    
  }

#ifdef __VECTOR_HACK
  operator const float*() const { return (const float*) this; }
#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 Vector3 Abs(void) const {
     return Vector3(fabs(x), fabs(y), fabs(z));
   }

  // normalizes the vector of unit size corresponding to given vector
  inline void Normalize();
  
  /** Returns false if this vector has a nan component.
  */
  bool CheckValidity() const;

  /**
    ===> Using ArbitraryNormal() for constructing coord systems 
    ===> is obsoleted by RightHandedBase() method (<JK> 12/20/03).   

     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 Vector3 ArbitraryNormal(const Vector3 &v); 

  /**
    Find a right handed coordinate system with (*this) being
    the z-axis. For a right-handed system, U x V = (*this) holds.
    This implementation is here to avoid inconsistence and confusion
    when construction coordinate systems using ArbitraryNormal():
    In fact:
      V = ArbitraryNormal(N);
      U = CrossProd(V,N);
    constructs a right-handed coordinate system as well, BUT:
    1) bugs can be introduced if one mistakenly constructs a 
       left handed sytems e.g. by doing
       U = ArbitraryNormal(N);
       V = CrossProd(U,N);
    2) this implementation gives non-negative base vectors
       for (*this)==(0,0,1) |  (0,1,0) | (1,0,0), which is
       good for debugging and is not the case with the implementation
       using ArbitraryNormal()

    ===> Using ArbitraryNormal() for constructing coord systems 
	     is obsoleted by this method (<JK> 12/20/03).   
  */
  void RightHandedBase(Vector3& U, Vector3& V) const;

  /// Transforms a vector to the global coordinate frame.
  /**
    Given a local coordinate frame (U,V,N) (i.e. U,V,N are 
    the x,y,z axes of the local coordinate system) and
    a vector 'loc' in the local coordiante system, this
    function returns a the coordinates of the same vector
    in global frame (i.e. frame (1,0,0), (0,1,0), (0,0,1).
  */
  friend inline Vector3 ToGlobalFrame(const Vector3& loc,
	  const Vector3& U,
	  const Vector3& V,
	  const Vector3& N);
  
  /// Transforms a vector to a local coordinate frame.
  /**
    Given a local coordinate frame (U,V,N) (i.e. U,V,N are 
    the x,y,z axes of the local coordinate system) and
    a vector 'loc' in the global coordiante system, this
    function returns a the coordinates of the same vector
    in the local frame.
  */
  friend inline Vector3 ToLocalFrame(const Vector3& loc,
	  const Vector3& U,
	  const Vector3& V,
	  const Vector3& N);

  /// the magnitude=size of the vector
  friend inline float Magnitude(const Vector3 &v);
  /// the squared magnitude of the vector .. for efficiency in some cases
  friend inline float SqrMagnitude(const Vector3 &v);
  /// Magnitude(v1-v2)
  friend inline float Distance(const Vector3 &v1, const Vector3 &v2);
  /// SqrMagnitude(v1-v2)
  friend inline float SqrDistance(const Vector3 &v1, const Vector3 &v2);

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

  // Rotate a normal vector.
  friend Vector3 PlaneRotate(const Matrix4x4 &, const Vector3 &);

  // construct view vectors .. DirAt is the main viewing direction
  // Viewer is the coordinates of viewer location, UpL is the vector.
  friend void ViewVectors(const Vector3 &DirAt, const Vector3 &Viewer,
						  const Vector3 &UpL, Vector3 &ViewV,
						  Vector3 &ViewU, Vector3 &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(Vector3 &U, Vector3 &V, // output
			     const Vector3 &normal, // input
			     const Vector3 &dirIncoming);
  // Unary operators
  Vector3 operator+ () const;
  Vector3 operator- () const;

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

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

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

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

  friend ostream& operator<< (ostream &s, const Vector3 &A);
  friend istream& operator>> (istream &s, Vector3 &A);
    
  friend void Minimize(Vector3 &min, const Vector3 &candidate);
  friend void Maximize(Vector3 &max, const Vector3 &candidate);

  friend inline int EpsilonEqualV3(const Vector3 &v1, const Vector3 &v2, float thr);
  friend inline int EpsilonEqualV3(const Vector3 &v1, const Vector3 &v2);

  friend Vector3 CosineRandomVector(const Vector3 &normal);
  friend Vector3 UniformRandomVector(const Vector3 &normal);
  friend Vector3 UniformRandomVector();

  friend Vector3
  UniformRandomVector(const float r1, const float r2);

  friend Vector3 CosineRandomVector(
									const float r1,
									const float r2,
									const Vector3 &normal
									);


};

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


inline void
Vector3::RightHandedBase(Vector3& U, Vector3& V) const
{
	// HACK
	V = ArbitraryNormal(*this);
	U = CrossProd(V, *this);
}


inline Vector3
ToGlobalFrame(const Vector3 &loc,
	      const Vector3 &U,
	      const Vector3 &V,
	      const Vector3 &N)
{
  return loc.x * U + loc.y * V + loc.z * N;
}

inline Vector3
ToLocalFrame(const Vector3 &loc,
	     const Vector3 &U,
	     const Vector3 &V,
	     const Vector3 &N)
{
  return Vector3( loc.x * U.x + loc.y * U.y + loc.z * U.z,
		    loc.x * V.x + loc.y * V.y + loc.z * V.z,
		    loc.x * N.x + loc.y * N.y + loc.z * N.z);
}

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

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

inline float
Distance(const Vector3 &v1, const Vector3 &v2)
{
  return sqrtf(sqr(v1.x-v2.x) + sqr(v1.y-v2.y) + sqr(v1.z-v2.z));
}

inline float
SqrDistance(const Vector3 &v1, const Vector3 &v2)
{
  return sqr(v1.x-v2.x)+sqr(v1.y-v2.y)+sqr(v1.z-v2.z);
}

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

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

// angle between two vectors with respect to a surface normal in the
// range [0 .. 2 * pi]
inline float
Angle(const Vector3 &A, const Vector3 &B, const Vector3 &norm)
{
	Vector3 cross = CrossProd(A, B);

	float signedAngle;

	if (DotProd(cross, norm) > 0)
		signedAngle = atan2(-Magnitude(CrossProd(A, B)), DotProd(A, B));
	else
		signedAngle = atan2(Magnitude(CrossProd(A, B)), DotProd(A, B));

	if (signedAngle < 0)
		return 2 * PI + signedAngle;

	return signedAngle;
}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

inline int
operator<= (const Vector3 &A, const Vector3 &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 Vector3 &A, const Vector3 &B)
{
  return (A.x == B.x) && (A.y == B.y) && (A.z == B.z);
}

inline Vector3
CrossProd (const Vector3 &A, const Vector3 &B)
{
  return 
	Vector3(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
Vector3::Normalize()
{
  float sqrmag = x * x + y * y + z * z;
  if (sqrmag > 0.0f)
    (*this) *= 1.0f / sqrtf(sqrmag);
}


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

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

inline int
EpsilonEqualV3(const Vector3 &v1, const Vector3 &v2, float thr)
{
  if ( fabsf(v1.x-v2.x) > thr )
    return false;
  if ( fabsf(v1.y-v2.y) > thr )
    return false;
  if ( fabsf(v1.z-v2.z) > thr )
    return false;
  return true;
}

inline int
EpsilonEqualV3(const Vector3 &v1, const Vector3 &v2)
{
  return EpsilonEqualV3(v1, v2, Limits::Small);
}



}

#endif