Curve.cpp

// Curve.cpp
// Copyright 2011, Dan Heeks
// This program is released under the BSD license. See the file COPYING for details.

#include "Curve.h"
#include "Circle.h"
#include "Arc.h"
#include "Area.h"
#include "kurve/geometry.h"

const Point operator*(const double &d, const Point &p){ return p * d;}
double Point::tolerance = 0.001;

//static const double PI = 3.1415926535897932; duplicated in kurve/geometry.h

double Point::length()const
{
    return sqrt( x*x + y*y );
}

double Point::normalize()
{
	double len = length();
	if(fabs(len)> 0.000000000000001)
		*this = (*this) / len;
	return len;
}

Line::Line(const Point& P0, const Point& V):p0(P0), v(V)
{
}

double Line::Dist(const Point& p)const
{
	Point vn = v;
	vn.normalize();
	double d1 = p0 * vn;
	double d2 = p * vn;
	Point pn = p0 + vn * (d2 - d1);

	return pn.dist(p);
}

CVertex::CVertex(int type, const Point& p, const Point& c, int user_data):m_type(type), m_p(p), m_c(c), m_user_data(user_data)
{
}

CVertex::CVertex(const Point& p, int user_data):m_type(0), m_p(p), m_c(0.0, 0.0), m_user_data(user_data)
{
}

void CCurve::append(const CVertex& vertex)
{
	m_vertices.push_back(vertex);
}

bool CCurve::CheckForArc(const CVertex& prev_vt, std::list<const CVertex*>& might_be_an_arc, Arc &arc_returned)
{
	// this examines the vertices in might_be_an_arc
	// if they do fit an arc, set arc to be the arc that they fit and return true
	// returns true, if arc added
	if(might_be_an_arc.size() < 2)return false;

	// find middle point
	int num = might_be_an_arc.size();
	int i = 0;
	const CVertex* mid_vt = NULL;
	int mid_i = (num-1)/2;
	for(std::list<const CVertex*>::iterator It = might_be_an_arc.begin(); It != might_be_an_arc.end(); It++, i++)
	{
		if(i == mid_i)
		{
			mid_vt = *It;
			break;
		}
	}

	// create a circle to test
	Point p0(prev_vt.m_p);
	Point p1(mid_vt->m_p);
	Point p2(might_be_an_arc.back()->m_p);
	Circle c(p0, p1, p2);

	const CVertex* current_vt = &prev_vt;
	double accuracy = CArea::m_accuracy * 1.4 / CArea::m_units;
	for(std::list<const CVertex*>::iterator It = might_be_an_arc.begin(); It != might_be_an_arc.end(); It++)
	{
		const CVertex* vt = *It;

		if(!c.LineIsOn(current_vt->m_p, vt->m_p, accuracy))
			return false;
		current_vt = vt;
	}

	Arc arc;
	arc.m_c = c.m_c;
	arc.m_s = prev_vt.m_p;
	arc.m_e = might_be_an_arc.back()->m_p;
	arc.SetDirWithPoint(might_be_an_arc.front()->m_p);
	arc.m_user_data = might_be_an_arc.back()->m_user_data;

	double angs = atan2(arc.m_s.y - arc.m_c.y, arc.m_s.x - arc.m_c.x);
	double ange = atan2(arc.m_e.y - arc.m_c.y, arc.m_e.x - arc.m_c.x);
	if(arc.m_dir)
	{
		// make sure ange > angs
		if(ange < angs)ange += 6.2831853071795864;
	}
	else
	{
		// make sure angs > ange
		if(angs < ange)angs += 6.2831853071795864;
	}

	if(arc.IncludedAngle() >= 3.15)return false; // We don't want full arcs, so limit to about 180 degrees

	for(std::list<const CVertex*>::iterator It = might_be_an_arc.begin(); It != might_be_an_arc.end(); It++)
	{
		const CVertex* vt = *It;
		double angp = atan2(vt->m_p.y - arc.m_c.y, vt->m_p.x - arc.m_c.x);
		if(arc.m_dir)
		{
			// make sure angp > angs
			if(angp < angs)angp += 6.2831853071795864;
			if(angp > ange)return false;
		}
		else
		{
			// make sure angp > ange
			if(angp < ange)angp += 6.2831853071795864;
			if(angp > angs)return false;
		}
	}

	arc_returned = arc;
	return true;
}

void CCurve::AddArcOrLines(bool check_for_arc, std::list<CVertex> &new_vertices, std::list<const CVertex*>& might_be_an_arc, Arc &arc, bool &arc_found, bool &arc_added)
{
	if(check_for_arc && CheckForArc(new_vertices.back(), might_be_an_arc, arc))
	{
		arc_found = true;
	}
	else
	{
		if(arc_found)
		{
			if(arc.AlmostALine())
			{
				new_vertices.push_back(CVertex(arc.m_e, arc.m_user_data));
			}
			else
			{
				new_vertices.push_back(CVertex(arc.m_dir ? 1:-1, arc.m_e, arc.m_c, arc.m_user_data));
			}

			arc_added = true;
			arc_found = false;
			const CVertex* back_vt = might_be_an_arc.back();
			might_be_an_arc.clear();
			if(check_for_arc)might_be_an_arc.push_back(back_vt);
		}
		else
		{
			const CVertex* back_vt = might_be_an_arc.back();
			if(check_for_arc)might_be_an_arc.pop_back();
			for(std::list<const CVertex*>::iterator It = might_be_an_arc.begin(); It != might_be_an_arc.end(); It++)
			{
				const CVertex* v = *It;
				if(It != might_be_an_arc.begin() || (new_vertices.size() == 0) || (new_vertices.back().m_p != v->m_p))
				{
					new_vertices.push_back(*v);
				}
			}
			might_be_an_arc.clear();
			if(check_for_arc)might_be_an_arc.push_back(back_vt);
		}
	}
}

void CCurve::FitArcs()
{
	std::list<CVertex> new_vertices;

	std::list<const CVertex*> might_be_an_arc;
	Arc arc;
	bool arc_found = false;
	bool arc_added = false;

	int i = 0;
	for(std::list<CVertex>::iterator It = m_vertices.begin(); It != m_vertices.end(); It++, i++)
	{
		CVertex& vt = *It;
		if(vt.m_type || i == 0)
			new_vertices.push_back(vt);
		else
		{
			might_be_an_arc.push_back(&vt);

			if(might_be_an_arc.size() == 1)
			{
			}
			else
			{
				AddArcOrLines(true, new_vertices, might_be_an_arc, arc, arc_found, arc_added);
			}
		}
	}

	if(might_be_an_arc.size() > 0)AddArcOrLines(false, new_vertices, might_be_an_arc, arc, arc_found, arc_added);

	if(arc_added)
	{
		m_vertices.clear();
		for(std::list<CVertex>::iterator It = new_vertices.begin(); It != new_vertices.end(); It++)m_vertices.push_back(*It);
		for(std::list<const CVertex*>::iterator It = might_be_an_arc.begin(); It != might_be_an_arc.end(); It++)m_vertices.push_back(*(*It));
	}
}

void CCurve::UnFitArcs()
{
	std::list<Point> new_pts;

	const CVertex* prev_vertex = NULL;
	for(std::list<CVertex>::const_iterator It2 = m_vertices.begin(); It2 != m_vertices.end(); It2++)
	{
		const CVertex& vertex = *It2;
		if(vertex.m_type == 0 || prev_vertex == NULL)
		{
			new_pts.push_back(vertex.m_p * CArea::m_units);
		}
		else
		{
			if(vertex.m_p != prev_vertex->m_p)
			{
				double phi,dphi,dx,dy;
				int Segments;
				int i;
				double ang1,ang2,phit;

				dx = (prev_vertex->m_p.x - vertex.m_c.x) * CArea::m_units;
				dy = (prev_vertex->m_p.y - vertex.m_c.y) * CArea::m_units;

				ang1=atan2(dy,dx);
				if (ang1<0) ang1+=2.0*PI;
				dx = (vertex.m_p.x - vertex.m_c.x) * CArea::m_units;
				dy = (vertex.m_p.y - vertex.m_c.y) * CArea::m_units;
				ang2=atan2(dy,dx);
				if (ang2<0) ang2+=2.0*PI;

				if (vertex.m_type == -1)
				{ //clockwise
					if (ang2 > ang1)
						phit=2.0*PI-ang2+ ang1;
					else
						phit=ang1-ang2;
				}
				else
				{ //counter_clockwise
					if (ang1 > ang2)
						phit=-(2.0*PI-ang1+ ang2);
					else
						phit=-(ang2-ang1);
				}

				//what is the delta phi to get an accurancy of aber
				double radius = sqrt(dx*dx + dy*dy);
				dphi=2*acos((radius-CArea::m_accuracy)/radius);

				//set the number of segments
				if (phit > 0)
					Segments=(int)ceil(phit/dphi);
				else
					Segments=(int)ceil(-phit/dphi);

				if (Segments < 1)
					Segments=1;
				if (Segments > 100)
					Segments=100;

				dphi=phit/(Segments);

				double px = prev_vertex->m_p.x * CArea::m_units;
				double py = prev_vertex->m_p.y * CArea::m_units;

				for (i=1; i<=Segments; i++)
				{
					dx = px - vertex.m_c.x * CArea::m_units;
					dy = py - vertex.m_c.y * CArea::m_units;
					phi=atan2(dy,dx);

					double nx = vertex.m_c.x * CArea::m_units + radius * cos(phi-dphi);
					double ny = vertex.m_c.y * CArea::m_units + radius * sin(phi-dphi);

					new_pts.push_back(Point(nx, ny));

					px = nx;
					py = ny;
				}
			}
		}
		prev_vertex = &vertex;
	}

	m_vertices.clear();

	for(std::list<Point>::iterator It = new_pts.begin(); It != new_pts.end(); It++)
	{
		Point &pt = *It;
		CVertex vertex(0, pt / CArea::m_units, Point(0.0, 0.0));
		m_vertices.push_back(vertex);
	}
}

Point CCurve::NearestPoint(const Point& p)const
{
	double best_dist = 0.0;
	Point best_point = Point(0, 0);
	bool best_point_valid = false;
	Point prev_p = Point(0, 0);
	bool prev_p_valid = false;
	bool first_span = true;
	for(std::list<CVertex>::const_iterator It = m_vertices.begin(); It != m_vertices.end(); It++)
	{
		const CVertex& vertex = *It;
		if(prev_p_valid)
		{
			Point near_point = Span(prev_p, vertex, first_span).NearestPoint(p);
			first_span = false;
			double dist = near_point.dist(p);
			if(!best_point_valid || dist < best_dist)
			{
				best_dist = dist;
				best_point = near_point;
				best_point_valid = true;
			}
		}
		prev_p = vertex.m_p;
		prev_p_valid = true;
	}
	return best_point;
}

Point CCurve::NearestPoint(const CCurve& c, double *d)const
{
	double best_dist = 0.0;
	Point best_point = Point(0, 0);
	bool best_point_valid = false;
	Point prev_p = Point(0, 0);
	bool prev_p_valid = false;
	bool first_span = true;
	for(std::list<CVertex>::const_iterator It = c.m_vertices.begin(); It != c.m_vertices.end(); It++)
	{
		const CVertex& vertex = *It;
		if(prev_p_valid)
		{
			double dist;
			Point near_point = NearestPoint(Span(prev_p, vertex, first_span), &dist);
			first_span = false;
			if(!best_point_valid || dist < best_dist)
			{
				best_dist = dist;
				best_point = near_point;
				best_point_valid = true;
			}
		}
		prev_p = vertex.m_p;
		prev_p_valid = true;
	}
	if(d)*d = best_dist;
	return best_point;
}

void CCurve::GetBox(CAreaBox &box)
{
	Point prev_p = Point(0, 0);
	bool prev_p_valid = false;
	for(std::list<CVertex>::iterator It = m_vertices.begin(); It != m_vertices.end(); It++)
	{
		CVertex& vertex = *It;
		if(prev_p_valid)
		{
			Span(prev_p, vertex).GetBox(box);
		}
		prev_p = vertex.m_p;
		prev_p_valid = true;
	}
}

void CCurve::Reverse()
{
	std::list<CVertex> new_vertices;

	CVertex* prev_v = NULL;

	for(std::list<CVertex>::reverse_iterator It = m_vertices.rbegin(); It != m_vertices.rend(); It++)
	{
		CVertex &v = *It;
		int type = 0;
		Point cp(0.0, 0.0);
		if(prev_v)
		{
			type = -prev_v->m_type;
			cp = prev_v->m_c;
		}
		CVertex new_v(type, v.m_p, cp);
		new_vertices.push_back(new_v);
		prev_v = &v;
	}

	m_vertices = new_vertices;
}

double CCurve::GetArea()const
{
	double area = 0.0;
	Point prev_p = Point(0, 0);
	bool prev_p_valid = false;
	for(std::list<CVertex>::const_iterator It = m_vertices.begin(); It != m_vertices.end(); It++)
	{
		const CVertex& vertex = *It;
		if(prev_p_valid)
		{
			area += Span(prev_p, vertex).GetArea();
		}
		prev_p = vertex.m_p;
		prev_p_valid = true;
	}
	return area;
}

bool CCurve::IsClosed()const
{
	if(m_vertices.size() == 0)return false;
	return m_vertices.front().m_p == m_vertices.back().m_p;
}

void CCurve::ChangeStart(const Point &p) {
	CCurve new_curve;

	bool started = false;
	bool finished = false;
	int start_span;

	for(int i = 0; i < 2; i++)
	{
		const Point *prev_p = NULL;

		int span_index = 0;
		for(std::list<CVertex>::const_iterator VIt = m_vertices.begin(); VIt != m_vertices.end() && !finished; VIt++)
		{
			const CVertex& vertex = *VIt;

			if(prev_p)
			{
				Span span(*prev_p, vertex);
				if(span.On(p))
				{
					if(started)
					{
						if(p == *prev_p || span_index != start_span)
						{
							new_curve.m_vertices.push_back(vertex);
						}
						else
						{
							if(p == vertex.m_p)new_curve.m_vertices.push_back(vertex);
							else
							{
								CVertex v(vertex);
								v.m_p = p;
								new_curve.m_vertices.push_back(v);
							}
							finished = true;
						}
					}
					else
					{
						new_curve.m_vertices.push_back(CVertex(p));
						started = true;
						start_span = span_index;
						if(p != vertex.m_p)new_curve.m_vertices.push_back(vertex);
					}
				}
				else
				{
					if(started)
					{
						new_curve.m_vertices.push_back(vertex);
					}
				}
				span_index++;
			}
			prev_p = &(vertex.m_p);
		}
	}

	if(finished)
	{
		*this = new_curve;
	}
}

void CCurve::Break(const Point &p) {
	// inserts a point, if it lies on the curve
	const Point *prev_p = NULL;

	for(std::list<CVertex>::iterator VIt = m_vertices.begin(); VIt != m_vertices.end(); VIt++)
	{
		CVertex& vertex = *VIt;

		if(p == vertex.m_p)break; // point is already on a vertex

		if(prev_p)
		{
			Span span(*prev_p, vertex);
			if(span.On(p))
			{
				CVertex v(vertex);
				v.m_p = p;
				m_vertices.insert(VIt, v);
				break;
			}
		}
		prev_p = &(vertex.m_p);
	}
}

void CCurve::RemoveTinySpans() {
	CCurve new_curve;

	std::list<CVertex>::const_iterator VIt = m_vertices.begin(); 
	new_curve.m_vertices.push_back(*VIt);
	VIt++;

	for(; VIt != m_vertices.end(); VIt++)
	{
		const CVertex& vertex = *VIt;

		if(vertex.m_type != 0 || new_curve.m_vertices.back().m_p.dist(vertex.m_p) > Point::tolerance)
		{
			new_curve.m_vertices.push_back(vertex);
		}
	}
    new_curve.m_recur_depth = m_recur_depth;
	*this = new_curve;
}

void CCurve::ChangeEnd(const Point &p) {
	// changes the end position of the Kurve, doesn't keep closed kurves closed
	CCurve new_curve;

	const Point *prev_p = NULL;

	for(std::list<CVertex>::const_iterator VIt = m_vertices.begin(); VIt != m_vertices.end(); VIt++)
	{
		const CVertex& vertex = *VIt;

		if(prev_p)
		{
			Span span(*prev_p, vertex);
			if(span.On(p))
			{
				CVertex v(vertex);
				v.m_p = p;
				new_curve.m_vertices.push_back(v);
				break;
			}
			else
			{
				if(p != vertex.m_p)new_curve.m_vertices.push_back(vertex);
			}
		}
		else
		{
			new_curve.m_vertices.push_back(vertex);
		}
		prev_p = &(vertex.m_p);
	}

	*this = new_curve;
}

Point CCurve::NearestPoint(const Span& p, double *d)const
{
	double best_dist = 0.0;
	Point best_point = Point(0, 0);
	bool best_point_valid = false;
	Point prev_p = Point(0, 0);
	bool prev_p_valid = false;
	bool first_span = true;
	for(std::list<CVertex>::const_iterator It = m_vertices.begin(); It != m_vertices.end(); It++)
	{
		const CVertex& vertex = *It;
		if(prev_p_valid)
		{
			double dist;
			Point near_point = Span(prev_p, vertex, first_span).NearestPoint(p, &dist);
			first_span = false;
			if(!best_point_valid || dist < best_dist)
			{
				best_dist = dist;
				best_point = near_point;
				best_point_valid = true;
			}
		}
		prev_p = vertex.m_p;
		prev_p_valid = true;
	}
	if(d)*d = best_dist;
	return best_point;
}

static geoff_geometry::Kurve MakeKurve(const CCurve& curve)
{
	geoff_geometry::Kurve k;
	for(std::list<CVertex>::const_iterator It = curve.m_vertices.begin(); It != curve.m_vertices.end(); It++)
	{
		const CVertex& v = *It;
		k.Add(geoff_geometry::spVertex(v.m_type, geoff_geometry::Point(v.m_p.x, v.m_p.y), geoff_geometry::Point(v.m_c.x, v.m_c.y)));
	}
	return k;
}

static CCurve MakeCCurve(const geoff_geometry::Kurve& k)
{
	CCurve c;
	int n = k.nSpans();
	for(int i = 0; i<= n; i++)
	{
		geoff_geometry::spVertex spv;
		k.Get(i, spv);
		c.append(CVertex(spv.type, Point(spv.p.x, spv.p.y), Point(spv.pc.x, spv.pc.y)));
	}
	return c;
}

bool CCurve::Offset(double leftwards_value)
{
	// use the kurve code donated by Geoff Hawkesford, to offset the curve as an open curve
	// returns true for success, false for failure
	geoff_geometry::Kurve k = MakeKurve(*this);
	geoff_geometry::Kurve kOffset;
	bool success = true;

	try
	{
		int ret = 0;
		k.OffsetMethod1(kOffset, fabs(leftwards_value), (leftwards_value > 0) ? 1:-1, 1, ret);
		success = (ret == 0);
		if(success)*this = MakeCCurve(kOffset);
	}
	catch(...)
	{
		success = false;
	}

	return success;
}

void CCurve::GetSpans(std::list<Span> &spans)const
{
	const Point *prev_p = NULL;
	for(std::list<CVertex>::const_iterator It = m_vertices.begin(); It != m_vertices.end(); It++)
	{
		const CVertex& vertex = *It;
		if(prev_p)
		{
			spans.push_back(Span(*prev_p, vertex));
		}
		prev_p = &(vertex.m_p);
	}
}

void CCurve::OffsetForward(double forwards_value, bool refit_arcs)
{
	// for drag-knife compensation

	// replace arcs with lines
	UnFitArcs();

	std::list<Span> spans;
	GetSpans(spans);

	m_vertices.clear();

	// shift all the spans
	for(std::list<Span>::iterator It = spans.begin(); It != spans.end(); It++)
	{
		Span &span = *It;
		Point v = span.GetVector(0.0);
		v.normalize();
		Point shift = v * forwards_value;
		span.m_p = span.m_p + shift;
		span.m_v.m_p = span.m_v.m_p + shift;
	}

	// loop through the shifted spans
	for(std::list<Span>::iterator It = spans.begin(); It != spans.end();)
	{
		Span &span = *It;
		Point v = span.GetVector(0.0);
		v.normalize();

		// add the span
		if(It == spans.begin())m_vertices.push_back(span.m_p);
		m_vertices.push_back(span.m_v.m_p);

		It++;
		if(It != spans.end())
		{
			Span &next_span = *It;
			Point nv = next_span.GetVector(0.0);
			nv.normalize();
			double sin_angle = v ^ nv;
			bool sharp_corner = ( fabs(sin_angle) > 0.5 ); // angle > 30 degrees

			if(sharp_corner)
			{
				// add an arc to the start of the next span
				int arc_type = ((sin_angle > 0) ? 1 : (-1));
				Point centre = span.m_v.m_p - v * forwards_value;
				m_vertices.push_back(CVertex(arc_type, next_span.m_p, centre));
			}
		}
	}

	if(refit_arcs)
		FitArcs(); // find the arcs again
	else
		UnFitArcs(); // convert those little arcs added to lines
}

double CCurve::Perim()const
{
	const Point *prev_p = NULL;
	double perim = 0.0;
	for(std::list<CVertex>::const_iterator It = m_vertices.begin(); It != m_vertices.end(); It++)
	{
		const CVertex& vertex = *It;
		if(prev_p)
		{
			Span span(*prev_p, vertex);
			perim += span.Length();
		}
		prev_p = &(vertex.m_p);
	}

	return perim;
}

Point CCurve::PerimToPoint(double perim)const
{
	if(m_vertices.size() == 0)return Point(0, 0);

	const Point *prev_p = NULL;
	double kperim = 0.0;
	for(std::list<CVertex>::const_iterator It = m_vertices.begin(); It != m_vertices.end(); It++)
	{
		const CVertex& vertex = *It;
		if(prev_p)
		{
			Span span(*prev_p, vertex);
			double length = span.Length();
			if(perim < kperim + length)
			{
				Point p = span.MidPerim(perim - kperim);
				return p;
			}
			kperim += length;
		}
		prev_p = &(vertex.m_p);
	}

	return m_vertices.back().m_p;
}

double CCurve::PointToPerim(const Point& p)const
{
	double best_dist = 0.0;
	double perim_at_best_dist = 0.0;
	Point best_point = Point(0, 0);
	bool best_dist_found = false;

	double perim = 0.0;

	const Point *prev_p = NULL;
	bool first_span = true;
	for(std::list<CVertex>::const_iterator It = m_vertices.begin(); It != m_vertices.end(); It++)
	{
		const CVertex& vertex = *It;
		if(prev_p)
		{
			Span span(*prev_p, vertex, first_span);
			Point near_point = span.NearestPoint(p);
			first_span = false;
			double dist = near_point.dist(p);
			if(!best_dist_found || dist < best_dist)
			{
				best_dist = dist;
				Span span_to_point(*prev_p, CVertex(span.m_v.m_type, near_point, span.m_v.m_c));
				perim_at_best_dist = perim + span_to_point.Length();
				best_dist_found = true;
			}
			perim += span.Length();
		}
		prev_p = &(vertex.m_p);
	}
	return perim_at_best_dist;
}

void CCurve::operator+=(const CCurve& curve)
{
	for(std::list<CVertex>::const_iterator It = curve.m_vertices.begin(); It != curve.m_vertices.end(); It++)
	{
		const CVertex &vt = *It;
		if(It == curve.m_vertices.begin())
		{
			if((m_vertices.size() == 0) || (It->m_p != m_vertices.back().m_p))
			{
				m_vertices.push_back(CVertex(It->m_p));
			}
		}
		else
		{
			m_vertices.push_back(vt);
		}
	}
}

const Point Span::null_point = Point(0, 0);
const CVertex Span::null_vertex = CVertex(Point(0, 0));

Span::Span():m_start_span(false), m_p(null_point), m_v(null_vertex){}

Point Span::NearestPointNotOnSpan(const Point& p)const
{
	if(m_v.m_type == 0)
	{
		Point Vs = (m_v.m_p - m_p);
		Vs.normalize();
		double dp = (p - m_p) * Vs;		
		return (Vs * dp) + m_p;
	}
	else
	{
		double radius = m_p.dist(m_v.m_c);
		double r = p.dist(m_v.m_c);
		if(r < Point::tolerance)return m_p;
		Point vc = (m_v.m_c - p);
		return p + vc * ((r - radius) / r);
	}
}

Point Span::NearestPoint(const Point& p)const
{
	Point np = NearestPointNotOnSpan(p);
	double t = Parameter(np);
	if(t >= 0.0 && t <= 1.0)return np;

	double d1 = p.dist(this->m_p);
	double d2 = p.dist(this->m_v.m_p);

	if(d1 < d2)return this->m_p;
	else return m_v.m_p;
}

Point Span::MidPerim(double d)const {
	/// returns a point which is 0-d along span
	Point p;
	if(m_v.m_type == 0) {
		Point vs = m_v.m_p - m_p;
		vs.normalize();
		p = vs * d + m_p;
	}
	else {
		Point v = m_p - m_v.m_c;
		double radius = v.length();
		v.Rotate(d * m_v.m_type / radius);
		p = v + m_v.m_c;
	}
	return p;
}

Point Span::MidParam(double param)const {
	/// returns a point which is 0-1 along span
	if(fabs(param) < 0.00000000000001)return m_p;
	if(fabs(param - 1.0) < 0.00000000000001)return m_v.m_p;

	Point p;
	if(m_v.m_type == 0) {
		Point vs = m_v.m_p - m_p;
		p = vs * param + m_p;
	}
	else {
		Point v = m_p - m_v.m_c;
		v.Rotate(param * IncludedAngle());
		p = v + m_v.m_c;
	}
	return p;
}

Point Span::NearestPointToSpan(const Span& p, double &d)const
{
	Point midpoint = MidParam(0.5);
	Point np = p.NearestPoint(m_p);
	Point best_point = m_p;
	double dist = np.dist(m_p);
	if(p.m_start_span)dist -= (CArea::m_accuracy * 2); // give start of curve most priority
	Point npm = p.NearestPoint(midpoint);
	double dm = npm.dist(midpoint) - CArea::m_accuracy; // lie about midpoint distance to give midpoints priority
	if(dm < dist){dist = dm; best_point = midpoint;}
	Point np2 = p.NearestPoint(m_v.m_p);
	double dp2 = np2.dist(m_v.m_p);
	if(dp2 < dist){dist = dp2; best_point = m_v.m_p;}
	d = dist;
	return best_point;
}

Point Span::NearestPoint(const Span& p, double *d)const
{
	double best_dist;
	Point best_point = this->NearestPointToSpan(p, best_dist);

	// try the other way round too
	double best_dist2;
	Point best_point2 = p.NearestPointToSpan(*this, best_dist2);
	if(best_dist2 < best_dist)
	{
		best_point = NearestPoint(best_point2);
		best_dist = best_dist2;
	}

	if(d)*d = best_dist;
	return best_point;
}

static int GetQuadrant(const Point& v){
	// 0 = [+,+], 1 = [-,+], 2 = [-,-], 3 = [+,-]
	if(v.x > 0)
	{
		if(v.y > 0)
			return 0;
		return 3;
	}
	if(v.y > 0)
		return 1;
	return 2;
}

static Point QuadrantEndPoint(int i)
{
	if(i >3)i-=4;
	switch(i)
	{
	case 0:
		return Point(0.0,1.0);
	case 1:
		return Point(-1.0,0.0);
	case 2:
		return Point(0.0,-1.0);
	default:
		return Point(1.0,0.0);
	}
}

void Span::GetBox(CAreaBox &box)
{
	box.Insert(m_p);
	box.Insert(m_v.m_p);

	if(this->m_v.m_type)
	{
		// arc, add quadrant points
		Point vs = m_p - m_v.m_c;
		Point ve = m_v.m_p - m_v.m_c;
		int qs = GetQuadrant(vs);
		int qe = GetQuadrant(ve);
		if(m_v.m_type == -1)
		{
			// swap qs and qe
			int t=qs;
			qs = qe;
			qe = t;
		}

		if(qe<qs)qe = qe + 4;

		double rad = m_v.m_p.dist(m_v.m_c);

		for(int i = qs; i<qe; i++)
		{
			box.Insert(m_v.m_c + QuadrantEndPoint(i) * rad);
		}
	}
}

double IncludedAngle(const Point& v0, const Point& v1, int dir) {
	// returns the absolute included angle between 2 vectors in the direction of dir ( 1=acw  -1=cw)
	double inc_ang = v0 * v1;
	if(inc_ang > 1. - 1.0e-10) return 0;
	if(inc_ang < -1. + 1.0e-10)
		inc_ang = PI;  
	else {									// dot product,   v1 . v2  =  cos ang
		if(inc_ang > 1.0) inc_ang = 1.0;
		inc_ang = acos(inc_ang);									// 0 to pi radians

		if(dir * (v0 ^ v1) < 0) inc_ang = 2 * PI - inc_ang ;		// cp
	}
	return dir * inc_ang;
}

double Span::IncludedAngle()const
{
	if(m_v.m_type)
	{
		Point vs = ~(m_p - m_v.m_c);
		Point ve = ~(m_v.m_p - m_v.m_c);
		if(m_v.m_type == -1)
		{
			vs = -vs;
			ve = -ve;
		}
		vs.normalize();
		ve.normalize();

		return ::IncludedAngle(vs, ve, m_v.m_type);
	}

	return 0.0;
}

double Span::GetArea()const
{
	if(m_v.m_type)
	{
		double angle = IncludedAngle();
		double radius = m_p.dist(m_v.m_c);
		return ( 0.5 * ((m_v.m_c.x - m_p.x) * (m_v.m_c.y + m_p.y) - (m_v.m_c.x - m_v.m_p.x) * (m_v.m_c.y + m_v.m_p.y) - angle * radius * radius));
	}

	return 0.5 * (m_v.m_p.x - m_p.x) * (m_p.y + m_v.m_p.y);
}

double Span::Parameter(const Point& p)const
{
	double t;
	if(m_v.m_type == 0) {
		Point v0 = p - m_p;
		Point vs = m_v.m_p - m_p;
		double length = vs.length();
		vs.normalize();
		t = vs * v0;
		t = t / length;
	}
	else
	{
		// true if p lies on arc span sp (p must be on circle of span)
		Point vs = ~(m_p - m_v.m_c);
		Point v = ~(p - m_v.m_c);
		vs.normalize();
		v.normalize();
		if(m_v.m_type == -1){
			vs = -vs;
			v = -v;
		}
		double ang = ::IncludedAngle(vs, v, m_v.m_type);
		double angle = IncludedAngle();
		t = ang / angle;
	}
	return t;
}

bool Span::On(const Point& p, double* t)const
{
	if(p != NearestPoint(p))return false;
	if(t)*t = Parameter(p);
	return true;
}

double Span::Length()const
{
	if(m_v.m_type) {
		double radius = m_p.dist(m_v.m_c);
		return fabs(IncludedAngle()) * radius;
	}

	return m_p.dist(m_v.m_p);
}

Point Span::GetVector(double fraction)const
{
	/// returns the direction vector at point which is 0-1 along span
	if(m_v.m_type == 0){
		Point v(m_p, m_v.m_p);
		v.normalize();
		return v;
	}

	Point p= MidParam(fraction);
	Point v(m_v.m_c, p);
	v.normalize();
	if(m_v.m_type == 1)
	{
		return Point(-v.y, v.x);
	}
	else
	{
		return Point(v.y, -v.x);
	}
}

void tangential_arc(const Point &p0, const Point &p1, const Point &v0, Point &c, int &dir)
{
	geoff_geometry::Point gp0(p0.x, p0.y);
	geoff_geometry::Point gp1(p1.x, p1.y);
	geoff_geometry::Vector2d gv0(v0.x, v0.y);
	geoff_geometry::Point gc;
	geoff_geometry::tangential_arc(gp0, gp1, gv0, gc, dir);
	c = Point(gc.x, gc.y);
}