#ifndef APPROX_CURVE_SEGMENT_H
#define APPROX_CURVE_SEGMENT_H 

#include <math.h>
#include <CGAL/basic.h>
#include <CGAL/Cartesian.h>
#include <CGAL/Point_2.h>
#include <CGAL/Vector_2.h>
#include <CGAL/Line_2.h>
#include <CGAL/Conic_2.h>
#include <CGAL/Segment_circle_2.h>
#include <CGAL/utils.h>
#include <CGAL/number_utils.h>
#include <CGAL/user_classes.h>
#include <CGAL/squared_distance_2.h>
#include <CGAL/enum.h>
#include <fstream>
#include <list>
#include <deque>
#include <iterator>


CGAL_BEGIN_NAMESPACE

template <class R >
class Arc_2
{

 public:


  typedef Point_2<R>           Point;

  typedef Circle_2<R>          Circle;
  typedef Segment_2<R>         Segment;

  private:
  int      _deg;                // The degree of the conic (either 1 or 2).
  Point    _source;             // The source of the arc. 
  Point    _target;             // The target of the arc.
  Point    _center;             // The center of the basic circle of the arc.
  typename R::FT _radius;       // The squared radius of the basic circle.
  Orientation    _orient;       // The orientation of the basic circle.

 public:
  // Default constructor.
  Arc_2 () :
    _deg(0)
        {}

  // Copy constructor.
  Arc_2 (const Arc_2<R>& arc) :
   
    _deg(arc._deg),
    _source(arc._source),
    _target(arc._target),
    _center(arc._center),
    _radius(arc._radius),
    _orient(arc._orient)
     {}

  // Constuct a segment arc from a segment.
  Arc_2 (const Segment& segment) :
    _deg(1),
    _source(segment.source()),
    _target(segment.target()) 
  {}

  // Construct a circular arc from a circle and two points on that circle

  Arc_2 (const Circle& circle,
	 const Point& source, const Point& target) :
    _deg(2),
    _source(source),
    _target(target),
    _center(circle.center()),
    _radius(circle.squared_radius()),
    _orient(circle.orientation())
  {
    
    CGAL_precondition(_source != _target);       
   
  }



   // Destructor.
  virtual ~Arc_2 ()
  {}

 

  // Get the arc's source.
  Point source () const
  {
    return (_source);
  }

  // Get the arc's target.
  Point target () const
  {
    return (_target);
  }

  // Check whether the curve is a segment.
  bool is_segment() const
  {
    return (_deg == 1);
  } 

 Segment segment() const
  {
    CGAL_precondition(is_segment());

    return (Segment (_source, _target));
  }

// Check whether the curve is a circle.
  bool is_circle() const
  {
    return (_deg == 2);
  }

  // Get a circle if the arc is indeed a circular arc.
  Circle circle() const
  {
    return  Circle (_center,_radius,_orient);
    
  }

  //Split the arc at the splitpoint
  Arc_2<R>*  split( Arc_2<R>& arc1,const Point& splitpoint)
    {
      Circle* circle1 = new Circle(arc1.source(),splitpoint,arc1.target());    
      //Construct the new arc, which will be inserted before arc1
      Arc_2<R>* Arc1 = new Arc_2<R>(*circle1,arc1.source(),splitpoint);
      Arc_2<R> Arc2(*circle1,splitpoint,arc1.target());
      //Change the source of arc1
      arc1 = Arc2;
      
      return Arc1;
    }

 typename R::FT
   find_farthest( const Arc_2<R> &arc, const Line_2<R> &line, Point &point)
   {
     typename R::FT  dist1(0),dist2(0),dist3(0),r(0);
     typedef typename R::Vector_2    Vector;
     typedef typename R::RT          Coefficient;
     bool                            max1on, max2on;
    
     if(arc.is_segment())
       { 
	 dist1 = squared_distance(arc.source(),line);
	 dist2 = squared_distance(arc.target(),line);
	 if (dist1 > dist2)
	   {
	     point = arc.source();
	     return dist1;
	   }
	 else
	   {
	     point = arc.target();
	     return dist2;
	   }
       }
     else
       {
	 
	 Circle circle = arc.circle();
	 Point m = circle.center();
	 
	 Coefficient a = line.a();         //compute the normal vector nl to the line
	 Coefficient b = line.b();
	 Vector nl(a,b);
	 
	 nl = nl/(sqrt(b*b+a*a));
	 
	 r = sqrt(circle.squared_radius());
	 
	 
	 nl = nl*r;                    //compute the first extreme point
	 Vector v1(ORIGIN,m);
	 v1 = v1+nl;
	 Point max1(v1.x(),v1.y());
	 
	 nl = -nl;                    //compute the secound extreme point
	 Vector v2(ORIGIN,m);
	 v2 = v2+nl;
	 Point max2(v2.x(),v2.y());
	 
	 if (circle.orientation() == COUNTERCLOCKWISE)  //check if extreme points are on the arc
	   {
	     
	     max1on = leftturn<R>(arc.source(),max1,arc.target());
	     
	     max2on = leftturn<R>(arc.source(),max2,arc.target());
	   }
	 else
	   {
	     max1on = rightturn<R>(arc.source(),max1,arc.target());
	     max2on = rightturn<R>(arc.source(),max2,arc.target());
	   }
	 
	 
	 if ((max1on==0) && (max2on==0))//The extreme points are not on the arc, source
	   {                           // or target has maximal distance to line
	     dist1 = squared_distance(arc.source(),line);
	     dist2 = squared_distance(arc.target(),line);
	     if (dist1 > dist2)
	       {
		 point = arc.source();
		 return dist1;
	       }
	     else
	       {
		 point = arc.target();
		 return dist2;
	       }
	   }
	 else                        
	   {
	     if (max1on && max2on)    //The  extreme points are on the arc,one of them  
	       {                     // is the point with max distance to the line
		 dist1 = squared_distance(max1,line);
		 dist2 = squared_distance(max2,line);
		 
		 if (dist1 > dist2)       // Get the extremum with the max distance to the line
		   {
		     point = max1;
		     return dist1;
		   }
		 else
		   {
		     point = max2;
		     return dist2;
		   }
	       }
	     else                  //Only one extreme point is on the arc
	       {
		 dist1 = squared_distance(arc.source(),line);
		 dist2 = squared_distance(arc.target(),line);
		 
		 if (max1on)
		   {
		     dist3 = squared_distance(max1,line);
		     if (dist1 >= dist2 && dist1 >dist3)
		       {
			 point = arc.source();
			 return dist1;
		       }
		     if (dist2 > dist1 && dist2 > dist3)
		       {
			 point = arc.target();
			 return dist2;
		       }
		     if (dist3 > dist1 && dist3 > dist2)
		       {
			 point = max1;
			 return dist3;
		       }
		     
		   }
		 else
		   {	    
		     
		     dist3 = squared_distance(max2,line);
		     
		     if ((dist1 >= dist2) && (dist1 >dist3))
		       {
			 point = arc.source();
			 return dist1;
		       }
		     if ((dist2 > dist1) && (dist2 > dist3))
		       {
			 point = arc.target();
			 
			 return dist2;
		       }
		     if ((dist3 > dist1) && (dist3 > dist2))
		       {
			 point = max2;
			 return dist3;
		       }
		   }
		 
	       }
	   }
       }
   }
};

template <class InputIterator, class OutputIterator, class Rep,class Iterator> 
OutputIterator approximate_curve_p(InputIterator first, InputIterator past,Iterator element, OutputIterator result,double eps, Rep param )
{
  double eps_2 = square(eps);
  std::list<Iterator>     list;  
  InputIterator h;
  
  //Create a list as a copy of the original container, because "split" needs to insert
  //elements 
  h=first;
  list.push_back(*h);
  int i=distance(first,past);
  int j;
  for (j=1;j!=i;j++)
    {    
      h++;
      list.push_back(*h);
    }
  
 
   Point_2<Rep> p;
   p=(**first).source();
   
  
  *result++ = (**first).source();
  
  //start the approximation 
  result =  approx_p(list.begin(), list.end(), result, eps_2, list, param);
 
  return result;
}




template <class InputIterator, class Container, class Rep,class Iterator> 
void approximate_curve_c(InputIterator first, InputIterator past, Iterator element,Container& result,double eps, Rep param )
{
  double eps_2 = square(eps);
  
  std::list<Iterator>     list;  
  InputIterator h;
  
  //Create a list as a copy of the original container, because "split" needs to insert
  //elements 
  h=first;
  list.push_back(*h);
  int i=distance(first,past);
  int j;
  for (j=1;j!=i;j++)
    {    
      h++;
      list.push_back(*h);
    }
  
 
   Point_2<Rep> p;
   p=(**first).source();
   
  
  result.push_back((**first).source());
  
  //start the approximation 
  approx_c(list.begin(), list.end(), result, eps_2, list, param);
 
 
}


template <class InputIterator, class OutputIterator, class R, class List>
OutputIterator approx_p(InputIterator first, InputIterator last, OutputIterator result,double  eps_2,List& list, R param)
{
  double maxdist, dist;
  typedef Point_2<R>   Point;
  typedef Arc_2<R>     Arc;
  typedef Circle_2<R> Circle;
  Point maxpoint,splitpoint;
  InputIterator   i,j,m,p;

  int d=distance(first,last);
 

  if (d==1) //Check if only one object has to be approximated 
    {
      Line_2<R> Line((**first).source(),(**first).target());
      
     maxdist = find_farthest(**first, Line, maxpoint);
      
     if (maxdist < eps_2)//if the error is smaller then the error bound - done
	{	
	  

	  *result++ = (**first).target();
	  return result;
	}
      else
	{

	  p = m =  first; 
	  splitpoint = maxpoint;

	  split(**p,splitpoint,p,list);
	 
	  //because of the insertion of one element (from split) the first iterator 
	  //needs re-arangement
	  first=p;
	  --first;
	}
    }
  else //more then one object 
    {
      i=last;
      --i;
     
      Line_2<R> Line((**first).source(),(**i).target());
      

      maxdist = find_farthest(**first, Line, maxpoint);
      p = first;
      splitpoint = maxpoint;

      j=first;
      ++j;  
      for ( ;j!=last; j++) //find the object with the farthest point to line
	                  //and compute the farthest point
	{

	  dist = find_farthest(**j, Line, maxpoint);
	  if (dist > maxdist)
	    {
	      maxdist = dist;
	      p = j;
	      splitpoint = maxpoint;

	    }
	}
        
      
      if (maxdist < eps_2) // the error is smaller then the error bound
	{
	  
	  *result++ = (**i).target();//the target of the last object is one of 
	  //vertices of the approximationg polygonal line
	
	  return result;
	}
      else
	{ 
	  if (splitpoint == (**p).target())
	    {
	      
	      ++p;	      

	    }
	  else
	    {
	      if (splitpoint != (**p).source())//the splitpoint is on the arc
		{
		  
		  split(**p,splitpoint,p,list);//the arc needs to be splited at the 
		  //splitpoint into two arcs
		  if (p==first)
		    first--;
		 
		}
		 
	    }
	}
      
      
    }
  
  approx_p(first,p,result,eps_2,list, param);//approximate the curve of the left interval
  approx_p(p,last,result,eps_2,list, param);//approximate the curve of the right interval
        
}


template <class InputIterator, class Container, class R, class List>
void approx_c(InputIterator first, InputIterator last, Container& result,double  eps_2,List& list, R param)
{
  double maxdist, dist;
  typedef Point_2<R>   Point;
  typedef Arc_2<R>     Arc;
  typedef Circle_2<R> Circle;
  Point maxpoint,splitpoint;
  InputIterator   i,j,m,p;

  int d=distance(first,last);
 

  if (d==1) //Check if only one object has to be approximated 
    {
      Line_2<R> Line((**first).source(),(**first).target());
      
     maxdist = find_farthest(**first, Line, maxpoint);
      
     if (maxdist < eps_2)//if the error is smaller then the error bound - done
	{	
	  

	  result.push_back((**first).target());
	  return;
	}
      else
	{

	  p = m =  first; 
	  splitpoint = maxpoint;

	  split(**p,splitpoint,p,list);
	 
	  //because of the insertion of one element (from split) the first iterator 
	  //needs re-arangement
	  first=p;
	  --first;
	}
    }
  else //more then one object 
    {
      i=last;
      --i;
     
      Line_2<R> Line((**first).source(),(**i).target());
      

      maxdist = find_farthest(**first, Line, maxpoint);
      p = first;
      splitpoint = maxpoint;

      j=first;
      ++j;  
      for ( ;j!=last; j++) //find the object with the farthest point to line
	                  //and compute the farthest point
	{

	  dist = find_farthest(**j, Line, maxpoint);
	  if (dist > maxdist)
	    {
	      maxdist = dist;
	      p = j;
	      splitpoint = maxpoint;

	    }
	}
        
      
      if (maxdist < eps_2) // the error is smaller then the error bound
	{
	  
	  result.push_back((**i).target());//the target of the last object is one of 
	  //vertices of the approximationg polygonal line
	
	  return;
	}
      else
	{ 
	  if (splitpoint == (**p).target())
	    {
	      
	      ++p;	      

	    }
	  else
	    {
	      if (splitpoint != (**p).source())//the splitpoint is on the arc
		{
		  
		  split(**p,splitpoint,p,list);//the arc needs to be splited at the 
		  //splitpoint into two arcs
		  if (p==first)
		    first--;
		 
		}
		 
	    }
	}
      
      
    }
  
  approx_c(first,p,result,eps_2,list, param);//approximate the curve of the left interval
  approx_c(p,last,result,eps_2,list, param);//approximate the curve of the right interval
}



template <class R, class Object,class List,class Iterator>
	       void split(Object &object,const Point_2<R> &splitpoint,const Iterator p, List &list) 
{
  Object* nobject;
 
  nobject = object.split(object,splitpoint);//Call the split function of the object

  list.insert(p,nobject);//insert the pointer of the new object in the list
  
}


template <class R,class Object>
typename R::FT
find_farthest( Object &object, const Line_2<R> &line, Point_2<R> &point)
{
  return object.find_farthest(object,line,point);//Call the find_farthest function of the object
 
}
CGAL_END_NAMESPACE
#endif