#ifndef TRIANGLESPACE_H
#define TRIANGLESPACE_H

#include <halfspaceD2.h>

/*!
  \brief Not inside the triangle partition.

  The boundary is not included, so inverting the isInside operator 
  would return true if a point was inside the triangle or on the
  triangles boundary.

*/
template< typename T, typename D >
class trianglespace 
{
public:

  /** Anti clockwise ordering of points. */
  T pi[3];

  /** Half spaces pointing outside of triangle.
      The half space is opposite the point with 
      the same local index. */
  halfspaceD2< T, D > hi[3];

  /** Anticlockwise ordering of points expected. */
  trianglespace(T const & p0, T const & p1, T const & p2)
  {
    pi[0] = p0;
    pi[1] = p1;
    pi[2] = p2;

    hi[0].set(pi[2],pi[1]);
    hi[1].set(pi[0],pi[2]);
    hi[2].set(pi[1],pi[0]);
  }

  /** Is the point outside the triangle boundary? */
  bool const isInside( T const & x ) const
  {
    if (hi[0].isInside(x))
      return true;

    if (hi[1].isInside(x))
      return true;

    if (hi[2].isInside(x))
      return true;

    return false;
  }

  bool const isOnEdge(T const & a, T const & b, T const & w, D const zero) const
  {
    assert( (b.dot()!=0) );
    D t = (w-a).dot(b);
    t /= b.dot();
    if (t<0)
       return false;
    if (t>1)
      return false;

    T p(a+b*t-w);
    D d = p.dot();
 
    if (d+zero<0)
      return false;

    if (d<zero)
      return true;

    return false;
  }

  bool const isOnBoundary(T const & x, D const zero ) const
  {
    if (isInside(x))
      return false;

    if (isOnEdge(pi[0],pi[1]-pi[0],x,zero))
      return true;
    
    if (isOnEdge(pi[0],pi[2]-pi[0],x,zero))
      return true;

    if (isOnEdge(pi[1],pi[2]-pi[1],x,zero))
      return true;

    return false;
  }
};


#endif