#ifndef __JGEOMETRY2DLIBRARY__
#define __JGEOMETRY2DLIBRARY__

#include <utility>
#include <algorithm>

#include "JGeometry2D/JPoint2D.hh"


namespace JGEOMETRY2D {

  namespace {


    /**
     * Check sequence of three points.
     * 
     * \param  a       first  point
     * \param  b       second point
     * \param  c       third  point
     * \return         true if points a, b, c are counter clockwise; else false
     */
    inline bool ccw(const JPoint2D& a, 
		    const JPoint2D& b, 
		    const JPoint2D& c) 
    {
      const double A = a.getX() - b.getX();
      const double B = a.getY() - b.getY();
      const double C = c.getX() - b.getX();
      const double D = c.getY() - b.getY();
      
      return A*D - B*C <= 0;
    }
    
    
    /**
     * Sort criterion for lower hull.
     *
     * \param  first   first  point
     * \param  second  second point
     * \return         true if first before second; else false
     */
    inline bool cmpl(const JPoint2D& first, const JPoint2D& second)
    {
      if (first.getX() - second.getX() == 0) 
	return second.getY() > first.getY();
      else
	return first.getX() > second.getX();
    }
    
    
    /**
     * Sort criterion for upper hull.
     *
     * \param  first   first  point
     * \param  second  second point
     * \return         true if first before second; else false
     */
    inline bool cmph(const JPoint2D& first, const JPoint2D& second)
    { 
      return cmpl(second,first); 
    }
    

    /**
     * Partition half Hull.
     *
     * \param  __begin     begin of data
     * \param  __end       end   of data
     * \param  cmp         comparator
     * \return             end   of data
     */
    template<class T>
    inline T getConvexHull2D(T __begin, 
			     T __end, 
			     bool (*cmp)(const JPoint2D&,const JPoint2D&))
    {
      using namespace std;

      if (__begin != __end) {
	
	sort(__begin, __end, cmp);
	
	T __l = __begin;
	
	if (++__l != __end) {
	  
	  T __i = __begin; 
	  
	  ++(++__i);
	  
	  for (T __j, __k; __i != __end; ++__i) {
	    
	    for (__j = __k = __l; __j != __begin && ccw(*__i, *__j, *--__k); --__j) {}
	    
	    __l = __j;
	    ++__l;

	    iter_swap(__l,__i);
	  }
	}
	
	return __l;
	
      } else

	return __begin;
    }
  }    
  

  /**
   * Get convex Hull.
   *
   * \param  __begin     begin of data
   * \param  __end       end of data
   * \return             end of lower and upper Hull data
   */
  template<class T>
  inline std::pair<T,T> getConvexHull2D(T __begin, T __end)
  {
    using namespace std;

    T __p = getConvexHull2D(__begin, __end, cmpl);        // make lower hull
    
    if (__p == __begin || __p == __end) 
      return make_pair(__p, __p);
    
    // add first point of lower hull to upper hull
    
    reverse(__begin, __p);
    
    T __q = getConvexHull2D(--__p, __end, cmph);          // make upper hull
    
    ++__q;
    
    // re-store order and keep the extra point
      
    reverse(__p,     __q);
    reverse(__begin, __q);
    
    int n = distance(__p, __q);
    
    __p = __begin;

    advance(__p, n);
    
    return make_pair(__p, __q);
  }
  

  /**
   * Get area of a convex Hull polygon.
   *
   * \param   __begin   begin of data
   * \param   __end     end   of data
   * \return            area
   */
  template<class T>
  inline double getArea2D(T __begin, T __end)
  {
    using namespace std;

    double A = 0;

    if (distance(__begin,__end) >= 3) {

      T __i, __j, __k;
      
      __i = __j = __k = __begin;
      
      for (++__j, ++(++__k); __k != __end; ++__i, ++__j, ++__k)
	A += __j->getX() * (__k->getY()  -  __i->getY());
      
      __k = __begin;
      
      A +=  __j->getX() * (__k->getY()  -  __i->getY());
      
      ++__i;
      __j = __k; 
      ++__k;
      
      A +=  __j->getX() * (__k->getY()  -  __i->getY());
    }

    return 0.5 * fabs(A);
  }
}

#endif