#ifndef TETRAHEDRONPARTITION_H
#define TETRAHEDRONPARTITION_H

#include <halfspaceD3.h>
#include <partitionspace.h>
#include <typedefs.h>

/*!  
\brief Test if a point is within a tetrahedron.
*/
template< typename T, typename D >
class tetrahedronpartition : public partitionspace<T>
{
public:

  /** The four points of the tetrahedron. */
  T pi[4];

  /** The halfspaces are opposite the vertex with 
      the same index, and point inwards. */
  halfspaceD3<T,D> hi[4];

  /** The first three points form a anticlockwise triangle
      facing inwards. ie its half space sees the fourth point. */
  tetrahedronpartition
  (
    T const & p0_,
    T const & p1_,
    T const & p2_,
    T const & p3_
  );
  /** Unconstructed partition. */
  tetrahedronpartition() {}
  /** The first three points form a anticlockwise triangle
      facing inwards. ie its half space sees the fourth point. */
  void construct
  (
    T const & p0_,
    T const & p1_,
    T const & p2_,
    T const & p3_
  );

  /** Is the point inside the tetrahedron? */
  boolc isInside( T const & x ) const; 

};


// --------------------------------------------------------- 
// Implementation


template< typename T, typename D >
tetrahedronpartition<T,D>::tetrahedronpartition
(
  T const & p0_,
  T const & p1_,
  T const & p2_,
  T const & p3_
)
{
  construct(p0_,p1_,p2_,p3_);
}


template< typename T, typename D >
void tetrahedronpartition<T,D>::construct
(
  T const & p0_,
  T const & p1_,
  T const & p2_,
  T const & p3_
)
{
  pi[0] = p0_;
  pi[1] = p1_;
  pi[2] = p2_;
  pi[3] = p3_;

  hi[3] = halfspaceD3<T,D>(p0_,p1_,p2_);
  assert( hi[3].isInside(pi[3])==true );
  hi[1] = halfspaceD3<T,D>(p0_,p2_,p3_);
  assert( hi[1].isInside(pi[1])==true );
  hi[2] = halfspaceD3<T,D>(p3_,p1_,p0_);
  assert( hi[2].isInside(pi[2])==true );
  hi[0] = halfspaceD3<T,D>(p2_,p1_,p3_);
  assert( hi[0].isInside(pi[0])==true );
}

template< typename T, typename D >
boolc tetrahedronpartition<T,D>::isInside( T const & x ) const
{
  for (uint i=0; i<4; ++i)
    if (!hi[i].isInside(x))
      return false;

  return true;
}

#endif