#ifndef _Plane3_H__
#define _Plane3_H__

#include "Vector3.h"

  
/** 3D Plane */
class Plane3 {
public:

  Vector3 mNormal;
  float mD;
  
  Plane3() {}
  
  Plane3(const Vector3 &a,
	 const Vector3 &b,
	 const Vector3 &c
	 ) {
    Vector3 v1=a-b, v2=c-b;
    mNormal = Normalize(CrossProd(v2,v1));
    mD = -DotProd(b, mNormal);
  }

  Plane3(const Vector3 &normal,
	 const Vector3 &point
	 ):mNormal(normal)
  {
    mD = -DotProd(normal, point);
  }


  float Distance(const Vector3 &v) const {
    return DotProd(v, mNormal) + mD;
  }

  int Side(const Vector3 &v, const float threshold = 1e-6) const {
    return signum(Distance(v), threshold);
  }

  Vector3 FindIntersection(const Vector3 &a,
			   const Vector3 &b,
			   float *t = NULL,
			   bool *coplanar = NULL
			   ) const 
  {
    const Vector3 v = b - a; // line from A to B
    float dv = DotProd(v, mNormal);
    
    if (signum(dv) == 0)
	{
		if (coplanar) (*coplanar) = true;	
		if (t) (*t) = 1;
		return a;
	}
	
    float u = - Distance(a) / dv; // TODO: could be done more efficiently
    
	if (coplanar) (*coplanar) = false;
	if (t) (*t) = u;
	
	//Debug << "t: " << u << ", b - a: " << v << ", norm: " << mNormal << ", dist(a): " << - Distance(a) << ", dv: " << dv << endl;
    //return a - Distance(a) * b / dv + Distance(a) * a / dv; // NOTE: gives better precision than calclulating a + u * v
    return a + u * v;
  }
  
  friend bool
  PlaneIntersection(const Plane3 &a, const Plane3 &b, const Plane3 &c, Vector3 &result);
  
  friend ostream &operator<<(ostream &s, const Plane3 p) {
    s<<p.mNormal<<" "<<p.mD;
    return s;
  }

  
};
  



#endif