#ifndef __JCIRCLE2D__
#define __JCIRCLE2D__

#include <istream>
#include <ostream>
#include <iterator>
#include <cmath>
#include <limits>

#include "JLang/JManip.hh"
#include "JLang/JVectorize.hh"
#include "JIO/JSerialisable.hh"
#include "JGeometry2D/JPosition2D.hh"


/**
 * \author mdejong
 */

namespace JGEOMETRY2D {}
namespace JPP { using namespace JGEOMETRY2D; }

namespace JGEOMETRY2D {

  using JLANG::array_type;
  using JIO::JReader;
  using JIO::JWriter;


  /**
   * Data structure for circle in two dimensions.
   */
  class JCircle2D :
    public JPosition2D
  {
  public:
    /**
     * Default constructor.
     */
    JCircle2D() :
      JPosition2D(),
      __r(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  point    point
     * \param  r        radius
     */
    JCircle2D(const JVector2D& point,
	      const double     r) :
      JPosition2D(point),
      __r(r)
    {}


    /**
     * Constructor.
     * Determines circle through two points.
     *
     * \param  p0         first  point
     * \param  p1         second point
     */
    JCircle2D(const JVector2D& p0,
	      const JVector2D& p1) :
      JPosition2D(),
      __r(0.0)
    {
      set(p0, p1);
    }


    /**
     * Constructor.
     * Determines circle through three points.
     *
     * \param  p0         first  point
     * \param  p1         second point
     * \param  p2         third  point
     * \param  precision  precision
     */
    JCircle2D(const JVector2D& p0,
	      const JVector2D& p1,
	      const JVector2D& p2,
	      const double precision = std::numeric_limits<double>::epsilon()) :
      JPosition2D(),
      __r(0.0)
    {
      set(p0, p1, p2, precision);
    }


    /**
     * Constructor.
     * Determines smallest enclosing circle around set of points.\n
     * The type of data should have the following member methods:
     * <pre>
     *     double getX();   // x coordinate
     *     double getY();   // y coordinate
     * </pre> 
     *
     * Reference:\n
     * Computational Geometry
     * Algorithms and Applications\n
     * Authors: de Berg, M., Cheong, O., van Kreveld, M., Overmars, M.
     *
     * \param  __begin    begin of data
     * \param  __end      end   of data
     * \param  precision  precision
     */
    template<class T>
    JCircle2D(T __begin,
	      T __end,
	      const double precision = std::numeric_limits<double>::epsilon()) :
      JPosition2D(),
      __r(0.0)
    {
      set(__begin, __end, precision);
    }


    /**
     * Constructor.
     * \param  buffer     input data
     * \param  precision  precision
     */
    template<class JElement_t, class JAllocator_t>
    JCircle2D(const array_type<JElement_t, JAllocator_t>& buffer,
	      const double precision = std::numeric_limits<double>::epsilon()) :
      JPosition2D(),
      __r(0.0)
    {
      set(buffer.begin(), buffer.end(), precision);
    }


    /**
     * Get radius.
     *
     * \return            radius
     */
    double getRadius() const 
    { 
      return __r;
    }


    /**
     * Set circle.
     * Determines circle through two points.
     *
     * \param  p0         first  point
     * \param  p1         second point
     */
    void set(const JVector2D& p0,
	     const JVector2D& p1)
    {
      __x = 0.5 * (p0.getX() + p1.getX());
      __y = 0.5 * (p0.getY() + p1.getY());	
      __r = this->getDistance(p0);
    }


    /**
     * Set circle.
     * Determines circle through three points.
     *
     * \param  p0         first  point
     * \param  p1         second point
     * \param  p2         third  point
     * \param  precision  precision
     */
    void set(const JVector2D& p0,
	     const JVector2D& p1,
	     const JVector2D& p2,
	     const double precision = std::numeric_limits<double>::epsilon())
    {
      const double x0 = p2.getX() - p1.getX();
      const double x1 = p0.getX() - p2.getX();
      const double x2 = p1.getX() - p0.getX();

      const double y0 = p1.getY() - p2.getY();
      const double y1 = p2.getY() - p0.getY();
      const double y2 = p0.getY() - p1.getY();

      const double D  = 2.0 * (p0.getX()*y0 + p1.getX()*y1 + p2.getX()*y2);

      if (fabs(D) > precision) {

	const double a = p0.getLengthSquared();
	const double b = p1.getLengthSquared();
	const double c = p2.getLengthSquared();

	__x = (a*y0 + b*y1 + c*y2) / D;
	__y = (a*x0 + b*x1 + c*x2) / D;
	__r = this->getDistance(p0);

      } else {

	set(p0, p1);

	const JCircle2D c1(p1, p2);
	const JCircle2D c2(p0, p2);

	if (c1.getRadius() > this->getRadius()) { *this = c1; }
	if (c2.getRadius() > this->getRadius()) { *this = c2; }
      }
    }


    /**
     * Set circle.
     * Determines smallest enclosing circle around set of points.
     * <pre>
     *     double getX();   // x coordinate
     *     double getY();   // y coordinate
     * </pre> 
     *
     * \param  __begin    begin of data
     * \param  __end      end   of data
     * \param  precision  precision
     */
    template<class T>
    void set(T __begin,
	     T __end,
	     const double precision = std::numeric_limits<double>::epsilon())
    {
      if (__begin != __end) {

	__x = __begin->getX();
	__y = __begin->getY();
	__r = 0.0;

	T i = __begin;
	
	const JVector2D p0(i->getX(), i->getY());

	while (++i != __end && p0.getDistance(JVector2D(i->getX(), i->getY())) <= precision) {}
	
	if (i != __end) {

	  const JVector2D p1(i->getX(), i->getY());

	  set(p0, p1);

	  while (++i != __end) {

	    const JVector2D p2(i->getX(), i->getY());

	    if (this->getDistance(p2) > this->getRadius() + precision  &&
		p0.getDistance(p2)    > precision                      &&
		p1.getDistance(p2)    > precision) {
	      configure(__begin, i, p2, precision);
	    }
	  }
	}
      }
    }


    /**
     * Check whether given point is inside circle.
     *
     * \param  pos         position
     * \param  precision   precision
     * \return             true if inside; else false
     */
    inline bool is_inside(const JVector2D& pos,
			  const double     precision = std::numeric_limits<double>::min()) const
    {
      return (this->getDistance(pos) <= this->getRadius() + precision);
    }


    /**
     * Read circle from input stream.
     *
     * \param  in         input stream
     * \param  circle     circle
     * \return            input stream
     */
    friend inline std::istream& operator>>(std::istream& in, JCircle2D& circle)
    {
      in >> static_cast<JPosition2D&>(circle);
      in >> circle.__r;

      return in;
    }


    /**
     * Write circle to output stream.
     *
     * \param  out        output stream
     * \param  circle     circle
     * \return            output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JCircle2D& circle)
    {
      const JFormat format(out, getFormat<JPosition2D>(JFormat_t(9, 3, std::ios::fixed | std::ios::showpos)));

      out << static_cast<const JPosition2D&>(circle);
      out << ' ';
      out << format << circle.getRadius();

      return out;
    }


    /**
     * Read circle from input.
     *
     * \param  in         reader
     * \param  circle     circle
     * \return            reader
     */
    friend inline JReader& operator>>(JReader& in, JCircle2D& circle)
    {
      in >> static_cast<JPosition2D&>(circle);
      in >> circle.__r;

      return in;
    }


    /**
     * Write circle to output.
     *
     * \param  out        writer
     * \param  circle     circle
     * \return            writer
     */
    friend inline JWriter& operator<<(JWriter& out, const JCircle2D& circle)
    {
      out << static_cast<const JPosition2D&>(circle);
      out << circle.__r;

      return out;
    }


  protected:
    double __r;


  private:
    /**
     * Determine smallest enclosing circle.
     *
     * \param  __begin    begin of data
     * \param  __end      end   of data
     * \param  p0         point on circle
     * \param  precision  precision
     */
    template<class T>
    void configure(T __begin, 
		   T __end, 
		   const JVector2D& p0,
		   const double precision)
    {
      this->set(JVector2D(__begin->getX(), __begin->getY()), p0);
      
      for (T i = __begin; ++i != __end; ) {

	const JVector2D p1(i->getX(), i->getY());

	if (this->getDistance(p1)     > this->getRadius() + precision  &&
	    p0.getDistanceSquared(p1) > precision) {
	  configure(__begin, i, p0, p1, precision);
	}
      }
    }


    /**
     * Determine smallest enclosing circle.
     *
     * \param  __begin    begin of data
     * \param  __end      end   of data
     * \param  p0         point on circle
     * \param  p1         point on circle
     * \param  precision  precision
     */
    template<class T>
    void configure(T __begin, 
		   T __end, 
		   const JVector2D& p0, 
		   const JVector2D& p1,
		   const double precision)
    {
      this->set(p0, p1);

      for (T i = __begin; i != __end; ++i) {

	const JVector2D p2(i->getX(), i->getY());

	if (this->getDistance(p2)     > this->getRadius() + precision  &&
	    p0.getDistanceSquared(p2) > precision                      &&
	    p1.getDistanceSquared(p2) > precision) {
	  this->set(p0, p1, p2, precision);
	}
      }
    }
  };
}

#endif