initial add of libarea files

src/Mod/Path/libarea/Arc.cpp

@@ -0,0 +1,122 @@
+// Arc.cpp
+
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#include "Arc.h"
+#include "Curve.h"
+
+void CArc::SetDirWithPoint(const Point& p)
+{
+	double angs = atan2(m_s.y - m_c.y, m_s.x - m_c.x);
+	double ange = atan2(m_e.y - m_c.y, m_e.x - m_c.x);
+	double angp = atan2(p.y - m_c.y, p.x - m_c.x);
+	if(ange < angs)ange += 6.2831853071795864;
+	if(angp < angs - 0.0000000000001)angp += 6.2831853071795864;
+	if(angp > ange + 0.0000000000001)m_dir = false;
+	else m_dir = true;
+}
+
+double CArc::IncludedAngle()const
+{
+	double angs = atan2(m_s.y - m_c.y, m_s.x - m_c.x);
+	double ange = atan2(m_e.y - m_c.y, m_e.x - m_c.x);
+	if(m_dir)
+	{
+		// make sure ange > angs
+		if(ange < angs)ange += 6.2831853071795864;
+	}
+	else
+	{
+		// make sure angs > ange
+		if(angs < ange)angs += 6.2831853071795864;
+	}
+
+	return fabs(ange - angs);
+}
+
+bool CArc::AlmostALine()const
+{
+	Point mid_point = MidParam(0.5);
+	if(Line(m_s, m_e - m_s).Dist(mid_point) <= Point::tolerance)
+		return true;
+
+	const double max_arc_radius = 1.0 / Point::tolerance;
+	double radius = m_c.dist(m_s);
+	if (radius > max_arc_radius)
+	{
+		return true;	// We don't want to produce an arc whose radius is too large.
+	}
+
+	return false;
+}
+
+Point CArc::MidParam(double param)const {
+	/// returns a point which is 0-1 along arc
+	if(fabs(param) < 0.00000000000001)return m_s;
+	if(fabs(param - 1.0) < 0.00000000000001)return m_e;
+
+	Point p;
+	Point v = m_s - m_c;
+	v.Rotate(param * IncludedAngle());
+	p = v + m_c;
+
+	return p;
+}
+
+//segments - number of segments per full revolution!
+//d_angle - determines the direction and the ammount of the arc to draw
+void CArc::GetSegments(void(*callbackfunc)(const double *p), double pixels_per_mm, bool want_start_point)const
+{
+	if(m_s == m_e)
+		return;
+
+	Point Va = m_s - m_c;
+	Point Vb = m_e - m_c;
+
+	double start_angle = atan2(Va.y, Va.x);
+	double end_angle = atan2(Vb.y, Vb.x);
+
+	if(m_dir)
+	{
+		if(start_angle > end_angle)end_angle += 6.28318530717958;
+	}
+	else
+	{
+		if(start_angle < end_angle)end_angle -= 6.28318530717958;
+	}
+
+	double radius = m_c.dist(m_s);
+	double d_angle = end_angle - start_angle;
+	int segments = (int)(fabs(pixels_per_mm * radius * d_angle / 6.28318530717958 + 1));
+
+    double theta = d_angle / (double)segments;
+	while(theta>1.0){segments*=2;theta = d_angle / (double)segments;}
+    double tangetial_factor = tan(theta);
+    double radial_factor = 1 - cos(theta);
+
+    double x = radius * cos(start_angle);
+    double y = radius * sin(start_angle);
+
+	double pp[3] = {0.0, 0.0, 0.0};
+
+   for(int i = 0; i < segments + 1; i++)
+    {
+		Point p = m_c + Point(x, y);
+		pp[0] = p.x;
+		pp[1] = p.y;
+		(*callbackfunc)(pp);
+
+        double tx = -y;
+        double ty = x;
+
+        x += tx * tangetial_factor;
+        y += ty * tangetial_factor;
+
+        double rx = - x;
+        double ry = - y;
+
+        x += rx * radial_factor;
+        y += ry * radial_factor;
+    }
+}

src/Mod/Path/libarea/Arc.h

@@ -0,0 +1,25 @@
+// Arc.h
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#pragma once
+
+#include "Point.h"
+
+class CArc{
+public:
+	Point m_s;
+	Point m_e;
+	Point m_c;
+	bool m_dir; // true - anti-clockwise, false - clockwise
+	int m_user_data;
+
+	CArc():m_dir(true), m_user_data(0){}
+	CArc(const Point& s, const Point& e, const Point& c, bool dir, int user_data):m_s(s), m_e(e), m_c(c), m_dir(dir), m_user_data(user_data){}
+
+	void SetDirWithPoint(const Point& p); // set m_dir, such that this point lies between m_s and m_e
+	double IncludedAngle()const; // always > 0
+	bool AlmostALine()const;
+	Point MidParam(double param)const;
+	void GetSegments(void(*callbackfunc)(const double *p), double pixels_per_mm, bool want_start_point = true)const;
+};

src/Mod/Path/libarea/Area.cpp

@@ -0,0 +1,727 @@
+// Area.cpp
+
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#include "Area.h"
+#include "AreaOrderer.h"
+
+#include <map>
+
+double CArea::m_accuracy = 0.01;
+double CArea::m_units = 1.0;
+bool CArea::m_fit_arcs = true;
+double CArea::m_single_area_processing_length = 0.0;
+double CArea::m_processing_done = 0.0;
+bool CArea::m_please_abort = false;
+double CArea::m_MakeOffsets_increment = 0.0;
+double CArea::m_split_processing_length = 0.0;
+bool CArea::m_set_processing_length_in_split = false;
+double CArea::m_after_MakeOffsets_length = 0.0;
+//static const double PI = 3.1415926535897932;
+
+void CArea::append(const CCurve& curve)
+{
+	m_curves.push_back(curve);
+}
+
+void CArea::FitArcs(){
+	for(std::list<CCurve>::iterator It = m_curves.begin(); It != m_curves.end(); It++)
+	{
+		CCurve& curve = *It;
+		curve.FitArcs();
+	}
+}
+
+Point CArea::NearestPoint(const Point& p)const
+{
+	double best_dist = 0.0;
+	Point best_point = Point(0, 0);
+	for(std::list<CCurve>::const_iterator It = m_curves.begin(); It != m_curves.end(); It++)
+	{
+		const CCurve& curve = *It;
+		Point near_point = curve.NearestPoint(p);
+		double dist = near_point.dist(p);
+		if(It == m_curves.begin() || dist < best_dist)
+		{
+			best_dist = dist;
+			best_point = near_point;
+		}
+	}
+	return best_point;
+}
+
+void CArea::GetBox(CBox2D &box)
+{
+	for(std::list<CCurve>::iterator It = m_curves.begin(); It != m_curves.end(); It++)
+	{
+		CCurve& curve = *It;
+		curve.GetBox(box);
+	}
+}
+
+void CArea::Reorder()
+{
+	// curves may have been added with wrong directions
+	// test all kurves to see which one are outsides and which are insides and 
+	// make sure outsides are anti-clockwise and insides are clockwise
+
+	// returns 0, if the curves are OK
+	// returns 1, if the curves are overlapping
+
+	CAreaOrderer ao;
+	for(std::list<CCurve>::iterator It = m_curves.begin(); It != m_curves.end(); It++)
+	{
+		CCurve& curve = *It;
+		ao.Insert(&curve);
+		if(m_set_processing_length_in_split)
+		{
+			CArea::m_processing_done += (m_split_processing_length / m_curves.size());
+		}
+	}
+
+	*this = ao.ResultArea();
+}
+
+class ZigZag
+{
+public:
+	CCurve zig;
+	CCurve zag;
+	ZigZag(const CCurve& Zig, const CCurve& Zag):zig(Zig), zag(Zag){}
+};
+
+static double stepover_for_pocket = 0.0;
+static std::list<ZigZag> zigzag_list_for_zigs;
+static std::list<CCurve> *curve_list_for_zigs = NULL;
+static bool rightward_for_zigs = true;
+static double sin_angle_for_zigs = 0.0;
+static double cos_angle_for_zigs = 0.0;
+static double sin_minus_angle_for_zigs = 0.0;
+static double cos_minus_angle_for_zigs = 0.0;
+static double one_over_units = 0.0;
+
+static Point rotated_point(const Point &p)
+{
+	return Point(p.x * cos_angle_for_zigs - p.y * sin_angle_for_zigs, p.x * sin_angle_for_zigs + p.y * cos_angle_for_zigs);
+}
+    
+static Point unrotated_point(const Point &p)
+{
+    return Point(p.x * cos_minus_angle_for_zigs - p.y * sin_minus_angle_for_zigs, p.x * sin_minus_angle_for_zigs + p.y * cos_minus_angle_for_zigs);
+}
+
+static CVertex rotated_vertex(const CVertex &v)
+{
+	if(v.m_type)
+	{
+		return CVertex(v.m_type, rotated_point(v.m_p), rotated_point(v.m_c));
+	}
+    return CVertex(v.m_type, rotated_point(v.m_p), Point(0, 0));
+}
+
+static CVertex unrotated_vertex(const CVertex &v)
+{
+	if(v.m_type)
+	{
+		return CVertex(v.m_type, unrotated_point(v.m_p), unrotated_point(v.m_c));
+	}
+	return CVertex(v.m_type, unrotated_point(v.m_p), Point(0, 0));
+}
+
+static void rotate_area(CArea &a)
+{
+	for(std::list<CCurve>::iterator It = a.m_curves.begin(); It != a.m_curves.end(); It++)
+	{
+		CCurve& curve = *It;
+		for(std::list<CVertex>::iterator CIt = curve.m_vertices.begin(); CIt != curve.m_vertices.end(); CIt++)
+		{
+			CVertex& vt = *CIt;
+			vt = rotated_vertex(vt);
+		}
+	}
+}
+
+void test_y_point(int i, const Point& p, Point& best_p, bool &found, int &best_index, double y, bool left_not_right)
+{
+	// only consider points at y
+	if(fabs(p.y - y) < 0.002 * one_over_units)
+	{
+		if(found)
+		{
+			// equal high point
+			if(left_not_right)
+			{
+				// use the furthest left point
+				if(p.x < best_p.x)
+				{
+					best_p = p;
+					best_index = i;
+				}
+			}
+			else
+			{
+				// use the furthest right point
+				if(p.x > best_p.x)
+				{
+					best_p = p;
+					best_index = i;
+				}
+			}
+		}
+		else
+		{
+			best_p = p;
+			best_index = i;
+			found = true;
+		}
+	}
+}
+
+static void make_zig_curve(const CCurve& input_curve, double y0, double y)
+{
+	CCurve curve(input_curve);
+
+	if(rightward_for_zigs)
+	{
+		if(curve.IsClockwise())
+			curve.Reverse();
+	}
+	else
+	{
+		if(!curve.IsClockwise())
+			curve.Reverse();
+	}
+
+    // find a high point to start looking from
+	Point top_left;
+	int top_left_index;
+	bool top_left_found = false;
+	Point top_right;
+	int top_right_index;
+	bool top_right_found = false;
+	Point bottom_left;
+	int bottom_left_index;
+	bool bottom_left_found = false;
+
+	int i =0;
+	for(std::list<CVertex>::const_iterator VIt = curve.m_vertices.begin(); VIt != curve.m_vertices.end(); VIt++, i++)
+	{
+		const CVertex& vertex = *VIt;
+
+		test_y_point(i, vertex.m_p, top_right, top_right_found, top_right_index, y, !rightward_for_zigs);
+		test_y_point(i, vertex.m_p, top_left, top_left_found, top_left_index, y, rightward_for_zigs);
+		test_y_point(i, vertex.m_p, bottom_left, bottom_left_found, bottom_left_index, y0, rightward_for_zigs);
+	}
+
+	int start_index = 0;
+	int end_index = 0;
+	int zag_end_index = 0;
+
+	if(bottom_left_found)start_index = bottom_left_index;
+	else if(top_left_found)start_index = top_left_index;
+
+	if(top_right_found)
+	{
+		end_index = top_right_index;
+		zag_end_index = top_left_index;
+	}
+	else
+	{
+		end_index = bottom_left_index;
+		zag_end_index =  bottom_left_index;
+	}
+	if(end_index <= start_index)end_index += (i-1);
+	if(zag_end_index <= start_index)zag_end_index += (i-1);
+
+    CCurve zig, zag;
+    
+    bool zig_started = false;
+    bool zig_finished = false;
+    bool zag_finished = false;
+    
+	int v_index = 0;
+	for(int i = 0; i < 2; i++)
+	{
+		// process the curve twice because we don't know where it will start
+		if(zag_finished)
+			break;
+		for(std::list<CVertex>::const_iterator VIt = curve.m_vertices.begin(); VIt != curve.m_vertices.end(); VIt++)
+		{
+			if(i == 1 && VIt == curve.m_vertices.begin())
+			{
+				continue;
+			}
+
+			const CVertex& vertex = *VIt;
+
+			if(zig_finished)
+			{
+				zag.m_vertices.push_back(unrotated_vertex(vertex));
+				if(v_index == zag_end_index)
+				{
+					zag_finished = true;
+					break;
+				}
+			}
+			else if(zig_started)
+			{
+				zig.m_vertices.push_back(unrotated_vertex(vertex));
+				if(v_index == end_index)
+				{
+					zig_finished = true;
+					if(v_index == zag_end_index)
+					{
+						zag_finished = true;
+						break;
+					}
+					zag.m_vertices.push_back(unrotated_vertex(vertex));
+				}
+			}
+			else
+			{
+				if(v_index == start_index)
+				{
+					zig.m_vertices.push_back(unrotated_vertex(vertex));
+					zig_started = true;
+				}
+			}
+			v_index++;
+		}
+	}
+        
+    if(zig_finished)
+		zigzag_list_for_zigs.push_back(ZigZag(zig, zag));
+}
+
+void make_zig(const CArea &a, double y0, double y)
+{
+	for(std::list<CCurve>::const_iterator It = a.m_curves.begin(); It != a.m_curves.end(); It++)
+	{
+		const CCurve &curve = *It;
+		make_zig_curve(curve, y0, y);
+	}
+}
+        
+std::list< std::list<ZigZag> > reorder_zig_list_list;
+        
+void add_reorder_zig(ZigZag &zigzag)
+{
+    // look in existing lists
+
+	// see if the zag is part of an existing zig
+	if(zigzag.zag.m_vertices.size() > 1)
+	{
+		const Point& zag_e = zigzag.zag.m_vertices.front().m_p;
+		bool zag_removed = false;
+		for(std::list< std::list<ZigZag> >::iterator It = reorder_zig_list_list.begin(); It != reorder_zig_list_list.end() && !zag_removed; It++)
+		{
+			std::list<ZigZag> &zigzag_list = *It;
+			for(std::list<ZigZag>::iterator It2 = zigzag_list.begin(); It2 != zigzag_list.end() && !zag_removed; It2++)
+			{
+				const ZigZag& z = *It2;
+				for(std::list<CVertex>::const_iterator It3 = z.zig.m_vertices.begin(); It3 != z.zig.m_vertices.end() && !zag_removed; It3++)
+				{
+					const CVertex &v = *It3;
+					if((fabs(zag_e.x - v.m_p.x) < (0.002 * one_over_units)) && (fabs(zag_e.y - v.m_p.y) < (0.002 * one_over_units)))
+					{
+						// remove zag from zigzag
+						zigzag.zag.m_vertices.clear();
+						zag_removed = true;
+					}
+				}
+			}
+		}
+	}
+
+	// see if the zigzag can join the end of an existing list
+	const Point& zig_s = zigzag.zig.m_vertices.front().m_p;
+	for(std::list< std::list<ZigZag> >::iterator It = reorder_zig_list_list.begin(); It != reorder_zig_list_list.end(); It++)
+	{
+		std::list<ZigZag> &zigzag_list = *It;
+		const ZigZag& last_zigzag = zigzag_list.back();
+        const Point& e = last_zigzag.zig.m_vertices.back().m_p;
+        if((fabs(zig_s.x - e.x) < (0.002 * one_over_units)) && (fabs(zig_s.y - e.y) < (0.002 * one_over_units)))
+		{
+            zigzag_list.push_back(zigzag);
+			return;
+		}
+	}
+        
+    // else add a new list
+    std::list<ZigZag> zigzag_list;
+    zigzag_list.push_back(zigzag);
+    reorder_zig_list_list.push_back(zigzag_list);
+}
+
+void reorder_zigs()
+{
+	for(std::list<ZigZag>::iterator It = zigzag_list_for_zigs.begin(); It != zigzag_list_for_zigs.end(); It++)
+	{
+		ZigZag &zigzag = *It;
+        add_reorder_zig(zigzag);
+	}
+        
+	zigzag_list_for_zigs.clear();
+
+	for(std::list< std::list<ZigZag> >::iterator It = reorder_zig_list_list.begin(); It != reorder_zig_list_list.end(); It++)
+	{
+		std::list<ZigZag> &zigzag_list = *It;
+		if(zigzag_list.size() == 0)continue;
+
+		curve_list_for_zigs->push_back(CCurve());
+		for(std::list<ZigZag>::const_iterator It = zigzag_list.begin(); It != zigzag_list.end();)
+		{
+			const ZigZag &zigzag = *It;
+			for(std::list<CVertex>::const_iterator It2 = zigzag.zig.m_vertices.begin(); It2 != zigzag.zig.m_vertices.end(); It2++)
+			{
+				if(It2 == zigzag.zig.m_vertices.begin() && It != zigzag_list.begin())continue; // only add the first vertex if doing the first zig
+				const CVertex &v = *It2;
+				curve_list_for_zigs->back().m_vertices.push_back(v);
+			}
+
+			It++;
+			if(It == zigzag_list.end())
+			{
+				for(std::list<CVertex>::const_iterator It2 = zigzag.zag.m_vertices.begin(); It2 != zigzag.zag.m_vertices.end(); It2++)
+				{
+					if(It2 == zigzag.zag.m_vertices.begin())continue; // don't add the first vertex of the zag
+					const CVertex &v = *It2;
+					curve_list_for_zigs->back().m_vertices.push_back(v);
+				}
+			}
+		}
+	}
+	reorder_zig_list_list.clear();
+}
+
+static void zigzag(const CArea &input_a)
+{
+	if(input_a.m_curves.size() == 0)
+	{
+		CArea::m_processing_done += CArea::m_single_area_processing_length;
+		return;
+	}
+    
+    one_over_units = 1 / CArea::m_units;
+    
+	CArea a(input_a);
+    rotate_area(a);
+    
+    CBox2D b;
+	a.GetBox(b);
+    
+    double x0 = b.MinX() - 1.0;
+    double x1 = b.MaxX() + 1.0;
+
+    double height = b.MaxY() - b.MinY();
+    int num_steps = int(height / stepover_for_pocket + 1);
+    double y = b.MinY();// + 0.1 * one_over_units;
+    Point null_point(0, 0);
+	rightward_for_zigs = true;
+
+	if(CArea::m_please_abort)return;
+
+	double step_percent_increment = 0.8 * CArea::m_single_area_processing_length / num_steps;
+
+	for(int i = 0; i<num_steps; i++)
+	{
+		double y0 = y;
+		y = y + stepover_for_pocket;
+		Point p0(x0, y0);
+		Point p1(x0, y);
+		Point p2(x1, y);
+		Point p3(x1, y0);
+		CCurve c;
+		c.m_vertices.push_back(CVertex(0, p0, null_point, 0));
+		c.m_vertices.push_back(CVertex(0, p1, null_point, 0));
+		c.m_vertices.push_back(CVertex(0, p2, null_point, 1));
+		c.m_vertices.push_back(CVertex(0, p3, null_point, 0));
+		c.m_vertices.push_back(CVertex(0, p0, null_point, 1));
+		CArea a2;
+		a2.m_curves.push_back(c);
+		a2.Intersect(a);
+		make_zig(a2, y0, y);
+		rightward_for_zigs = !rightward_for_zigs;
+		if(CArea::m_please_abort)return;
+		CArea::m_processing_done += step_percent_increment;
+	}
+
+	reorder_zigs();
+	CArea::m_processing_done += 0.2 * CArea::m_single_area_processing_length;
+}
+
+void CArea::SplitAndMakePocketToolpath(std::list<CCurve> &curve_list, const CAreaPocketParams &params)const
+{
+	CArea::m_processing_done = 0.0;
+
+	double save_units = CArea::m_units;
+	CArea::m_units = 1.0;
+	std::list<CArea> areas;
+	m_split_processing_length = 50.0; // jump to 50 percent after split
+	m_set_processing_length_in_split = true;
+	Split(areas);
+	m_set_processing_length_in_split = false;
+	CArea::m_processing_done = m_split_processing_length;
+	CArea::m_units = save_units;
+
+	if(areas.size() == 0)return;
+
+	double single_area_length = 50.0 / areas.size();
+
+	for(std::list<CArea>::iterator It = areas.begin(); It != areas.end(); It++)
+	{
+		CArea::m_single_area_processing_length = single_area_length;
+		CArea &ar = *It;
+		ar.MakePocketToolpath(curve_list, params);
+	}
+}
+
+void CArea::MakePocketToolpath(std::list<CCurve> &curve_list, const CAreaPocketParams &params)const
+{
+	double radians_angle = params.zig_angle * PI / 180;
+	sin_angle_for_zigs = sin(-radians_angle);
+	cos_angle_for_zigs = cos(-radians_angle);
+	sin_minus_angle_for_zigs = sin(radians_angle);
+	cos_minus_angle_for_zigs = cos(radians_angle);
+	stepover_for_pocket = params.stepover;
+
+	CArea a_offset = *this;
+	double current_offset = params.tool_radius + params.extra_offset;
+
+	a_offset.Offset(current_offset);
+
+	if(params.mode == ZigZagPocketMode || params.mode == ZigZagThenSingleOffsetPocketMode)
+	{
+		curve_list_for_zigs = &curve_list;
+		zigzag(a_offset);
+	}
+	else if(params.mode == SpiralPocketMode)
+	{
+		std::list<CArea> m_areas;
+		a_offset.Split(m_areas);
+		if(CArea::m_please_abort)return;
+		if(m_areas.size() == 0)
+		{
+			CArea::m_processing_done += CArea::m_single_area_processing_length;
+			return;
+		}
+
+		CArea::m_single_area_processing_length /= m_areas.size();
+
+		for(std::list<CArea>::iterator It = m_areas.begin(); It != m_areas.end(); It++)
+		{
+			CArea &a2 = *It;
+			a2.MakeOnePocketCurve(curve_list, params);
+		}
+	}
+
+	if(params.mode == SingleOffsetPocketMode || params.mode == ZigZagThenSingleOffsetPocketMode)
+	{
+		// add the single offset too
+		for(std::list<CCurve>::iterator It = a_offset.m_curves.begin(); It != a_offset.m_curves.end(); It++)
+		{
+			CCurve& curve = *It;
+			curve_list.push_back(curve);
+		}
+	}
+}
+
+void CArea::Split(std::list<CArea> &m_areas)const
+{
+	if(HolesLinked())
+	{
+		for(std::list<CCurve>::const_iterator It = m_curves.begin(); It != m_curves.end(); It++)
+		{
+			const CCurve& curve = *It;
+			m_areas.push_back(CArea());
+			m_areas.back().m_curves.push_back(curve);
+		}
+	}
+	else
+	{
+		CArea a = *this;
+		a.Reorder();
+
+		if(CArea::m_please_abort)return;
+
+		for(std::list<CCurve>::const_iterator It = a.m_curves.begin(); It != a.m_curves.end(); It++)
+		{
+			const CCurve& curve = *It;
+			if(curve.IsClockwise())
+			{
+				if(m_areas.size() > 0)
+					m_areas.back().m_curves.push_back(curve);
+			}
+			else
+			{
+				m_areas.push_back(CArea());
+				m_areas.back().m_curves.push_back(curve);
+			}
+		}
+	}
+}
+
+double CArea::GetArea(bool always_add)const
+{
+	// returns the area of the area
+	double area = 0.0;
+	for(std::list<CCurve>::const_iterator It = m_curves.begin(); It != m_curves.end(); It++)
+	{
+		const CCurve& curve = *It;
+		double a = curve.GetArea();
+		if(always_add)area += fabs(a);
+		else area += a;
+	}
+	return area;
+}
+
+eOverlapType GetOverlapType(const CCurve& c1, const CCurve& c2)
+{
+	CArea a1;
+	a1.m_curves.push_back(c1);
+	CArea a2;
+	a2.m_curves.push_back(c2);
+
+	return GetOverlapType(a1, a2);
+}
+
+eOverlapType GetOverlapType(const CArea& a1, const CArea& a2)
+{
+	CArea A1(a1);
+
+	A1.Subtract(a2);
+	if(A1.m_curves.size() == 0)
+	{
+		return eInside;
+	}
+
+	CArea A2(a2);
+	A2.Subtract(a1);
+	if(A2.m_curves.size() == 0)
+	{
+		return eOutside;
+	}
+
+	A1 = a1;
+	A1.Intersect(a2);
+	if(A1.m_curves.size() == 0)
+	{
+		return eSiblings;
+	}
+
+	return eCrossing;
+}
+
+bool IsInside(const Point& p, const CCurve& c)
+{
+	CArea a;
+	a.m_curves.push_back(c);
+	return IsInside(p, a);
+}
+
+bool IsInside(const Point& p, const CArea& a)
+{
+	CArea a2;
+	CCurve c;
+	c.m_vertices.push_back(CVertex(Point(p.x - 0.01, p.y - 0.01)));
+	c.m_vertices.push_back(CVertex(Point(p.x + 0.01, p.y - 0.01)));
+	c.m_vertices.push_back(CVertex(Point(p.x + 0.01, p.y + 0.01)));
+	c.m_vertices.push_back(CVertex(Point(p.x - 0.01, p.y + 0.01)));
+	c.m_vertices.push_back(CVertex(Point(p.x - 0.01, p.y - 0.01)));
+	a2.m_curves.push_back(c);
+	a2.Intersect(a);
+	if(fabs(a2.GetArea()) < 0.0004)return false;
+	return true;
+}
+
+void CArea::SpanIntersections(const Span& span, std::list<Point> &pts)const
+{
+	// this returns all the intersections of this area with the given span, ordered along the span
+
+	// get all points where this area's curves intersect the span
+	std::list<Point> pts2;
+	for(std::list<CCurve>::const_iterator It = m_curves.begin(); It != m_curves.end(); It++)
+	{
+		const CCurve &c = *It;
+		c.SpanIntersections(span, pts2);
+	}
+
+	// order them along the span
+	std::multimap<double, Point> ordered_points;
+	for(std::list<Point>::iterator It = pts2.begin(); It != pts2.end(); It++)
+	{
+		Point &p = *It;
+		double t;
+		if(span.On(p, &t))
+		{
+			ordered_points.insert(std::make_pair(t, p));
+		}
+	}
+
+	// add them to the given list of points
+	for(std::multimap<double, Point>::iterator It = ordered_points.begin(); It != ordered_points.end(); It++)
+	{
+		Point p = It->second;
+		pts.push_back(p);
+	}
+}
+
+void CArea::CurveIntersections(const CCurve& curve, std::list<Point> &pts)const
+{
+	// this returns all the intersections of this area with the given curve, ordered along the curve
+	std::list<Span> spans;
+	curve.GetSpans(spans);
+	for(std::list<Span>::iterator It = spans.begin(); It != spans.end(); It++)
+	{
+		Span& span = *It;
+		std::list<Point> pts2;
+		SpanIntersections(span, pts2);
+		for(std::list<Point>::iterator It = pts2.begin(); It != pts2.end(); It++)
+		{
+			Point &pt = *It;
+			if(pts.size() == 0)
+			{
+				pts.push_back(pt);
+			}
+			else
+			{
+				if(pt != pts.back())pts.push_back(pt);
+			}
+		}
+	}
+}
+
+class ThickLine
+{
+public:
+	CArea m_area;
+	CCurve m_curve;
+
+	ThickLine(const CCurve& curve)
+	{
+		m_curve = curve;
+		m_area.append(curve);
+		m_area.Thicken(0.001);
+	}
+};
+
+void CArea::InsideCurves(const CCurve& curve, std::list<CCurve> &curves_inside)const
+{
+	//1. find the intersectionpoints between these two curves.
+	std::list<Point> pts;
+	CurveIntersections(curve, pts);
+
+	//2.seperate curve2 in multiple curves between these intersections.
+	std::list<CCurve> separate_curves;
+	curve.ExtractSeparateCurves(pts, separate_curves);
+
+	//3. if the midpoint of a seperate curve lies in a1, then we return it.
+	for(std::list<CCurve>::iterator It = separate_curves.begin(); It != separate_curves.end(); It++)
+	{
+		CCurve &curve = *It;
+		double length = curve.Perim();
+		Point mid_point = curve.PerimToPoint(length * 0.5);
+		if(IsInside(mid_point, *this))curves_inside.push_back(curve);
+	}
+}

src/Mod/Path/libarea/Area.h

@@ -0,0 +1,89 @@
+// Area.h
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#ifndef AREA_HEADER
+#define AREA_HEADER
+
+#include "Curve.h"
+
+enum PocketMode
+{
+	SpiralPocketMode,
+	ZigZagPocketMode,
+	SingleOffsetPocketMode,
+	ZigZagThenSingleOffsetPocketMode,
+};
+
+struct CAreaPocketParams
+{
+	double tool_radius;
+	double extra_offset;
+	double stepover;
+	bool from_center;
+	PocketMode mode;
+	double zig_angle;
+	bool only_cut_first_offset;
+	CAreaPocketParams(double Tool_radius, double Extra_offset, double Stepover, bool From_center, PocketMode Mode, double Zig_angle)
+	{
+		tool_radius = Tool_radius;
+		extra_offset = Extra_offset;
+		stepover = Stepover;
+		from_center = From_center;
+		mode = Mode;
+		zig_angle = Zig_angle;
+	}
+};
+
+class CArea
+{
+public:
+	std::list<CCurve> m_curves;
+	static double m_accuracy;
+	static double m_units; // 1.0 for mm, 25.4 for inches. All points are multiplied by this before going to the engine
+	static bool m_fit_arcs;
+	static double m_processing_done; // 0.0 to 100.0, set inside MakeOnePocketCurve
+	static double m_single_area_processing_length;
+	static double m_after_MakeOffsets_length;
+	static double m_MakeOffsets_increment;
+	static double m_split_processing_length;
+	static bool m_set_processing_length_in_split;
+	static bool m_please_abort; // the user sets this from another thread, to tell MakeOnePocketCurve to finish with no result.
+
+	void append(const CCurve& curve);
+	void Subtract(const CArea& a2);
+	void Intersect(const CArea& a2);
+	void Union(const CArea& a2);
+	void Xor(const CArea& a2);
+	void Offset(double inwards_value);
+	void Thicken(double value);
+	void FitArcs();
+	unsigned int num_curves(){return m_curves.size();}
+	Point NearestPoint(const Point& p)const;
+	void GetBox(CBox2D &box);
+	void Reorder();
+	void MakePocketToolpath(std::list<CCurve> &toolpath, const CAreaPocketParams &params)const;
+	void SplitAndMakePocketToolpath(std::list<CCurve> &toolpath, const CAreaPocketParams &params)const;
+	void MakeOnePocketCurve(std::list<CCurve> &curve_list, const CAreaPocketParams &params)const;
+	static bool HolesLinked();
+	void Split(std::list<CArea> &m_areas)const;
+	double GetArea(bool always_add = false)const;
+	void SpanIntersections(const Span& span, std::list<Point> &pts)const; 
+	void CurveIntersections(const CCurve& curve, std::list<Point> &pts)const; 
+	void InsideCurves(const CCurve& curve, std::list<CCurve> &curves_inside)const;
+};
+
+enum eOverlapType
+{
+	eOutside,
+	eInside,
+	eSiblings,
+	eCrossing,
+};
+
+eOverlapType GetOverlapType(const CCurve& c1, const CCurve& c2);
+eOverlapType GetOverlapType(const CArea& a1, const CArea& a2);
+bool IsInside(const Point& p, const CCurve& c);
+bool IsInside(const Point& p, const CArea& a);
+
+#endif // #define AREA_HEADER

src/Mod/Path/libarea/AreaClipper.cpp

@@ -0,0 +1,479 @@
+// AreaClipper.cpp
+
+// implements CArea methods using Angus Johnson's "Clipper"
+
+#include "Area.h"
+#include "clipper.hpp"
+using namespace ClipperLib;
+
+#define TPolygon Path
+#define TPolyPolygon Paths
+
+bool CArea::HolesLinked(){ return false; }
+
+//static const double PI = 3.1415926535897932;
+static double Clipper4Factor = 10000.0;
+
+class DoubleAreaPoint
+{
+public:
+	double X, Y;
+
+	DoubleAreaPoint(double x, double y){X = x; Y = y;}
+	DoubleAreaPoint(const IntPoint& p){X = (double)(p.X) / Clipper4Factor; Y = (double)(p.Y) / Clipper4Factor;}
+	IntPoint int_point(){return IntPoint((long64)(X * Clipper4Factor), (long64)(Y * Clipper4Factor));}
+};
+
+static std::list<DoubleAreaPoint> pts_for_AddVertex;
+
+static void AddPoint(const DoubleAreaPoint& p)
+{
+	pts_for_AddVertex.push_back(p);
+}
+
+static void AddVertex(const CVertex& vertex, const CVertex* prev_vertex)
+{
+	if(vertex.m_type == 0 || prev_vertex == NULL)
+	{
+		AddPoint(DoubleAreaPoint(vertex.m_p.x * CArea::m_units, vertex.m_p.y * 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);
+
+			AddPoint(DoubleAreaPoint(nx, ny));
+
+			px = nx;
+			py = ny;
+		}
+		}
+	}
+}
+
+static bool IsPolygonClockwise(const TPolygon& p)
+{
+#if 1
+	double area = 0.0;
+	unsigned int s = p.size();
+	for(unsigned int i = 0; i<s; i++)
+	{
+		int im1 = i-1;
+		if(im1 < 0)im1 += s;
+
+		DoubleAreaPoint pt0(p[im1]);
+		DoubleAreaPoint pt1(p[i]);
+
+		area += 0.5 * (pt1.X - pt0.X) * (pt0.Y + pt1.Y);
+	}
+
+	return area > 0.0;
+#else
+	return IsClockwise(p);
+#endif
+}
+
+static void MakeLoop(const DoubleAreaPoint &pt0, const DoubleAreaPoint &pt1, const DoubleAreaPoint &pt2, double radius)
+{
+	Point p0(pt0.X, pt0.Y);
+	Point p1(pt1.X, pt1.Y);
+	Point p2(pt2.X, pt2.Y);
+	Point forward0 = p1 - p0;
+	Point right0(forward0.y, -forward0.x);
+	right0.normalize();
+	Point forward1 = p2 - p1;
+	Point right1(forward1.y, -forward1.x);
+	right1.normalize();
+
+	int arc_dir = (radius > 0) ? 1 : -1;
+
+	CVertex v0(0, p1 + right0 * radius, Point(0, 0));
+	CVertex v1(arc_dir, p1 + right1 * radius, p1);
+	CVertex v2(0, p2 + right1 * radius, Point(0, 0));
+
+	double save_units = CArea::m_units;
+	CArea::m_units = 1.0;
+
+	AddVertex(v1, &v0);
+	AddVertex(v2, &v1);
+
+	CArea::m_units = save_units;
+}
+
+static void OffsetWithLoops(const TPolyPolygon &pp, TPolyPolygon &pp_new, double inwards_value)
+{
+	Clipper c;
+
+	bool inwards = (inwards_value > 0);
+	bool reverse = false;
+	double radius = -fabs(inwards_value);
+
+	if(inwards)
+	{
+		// add a large square on the outside, to be removed later
+		TPolygon p;
+		p.push_back(DoubleAreaPoint(-10000.0, -10000.0).int_point());
+		p.push_back(DoubleAreaPoint(-10000.0, 10000.0).int_point());
+		p.push_back(DoubleAreaPoint(10000.0, 10000.0).int_point());
+		p.push_back(DoubleAreaPoint(10000.0, -10000.0).int_point());
+		c.AddPath(p, ptSubject, true);
+	}
+	else
+	{
+		reverse = true;
+	}
+
+	for(unsigned int i = 0; i < pp.size(); i++)
+	{
+		const TPolygon& p = pp[i];
+
+		pts_for_AddVertex.clear();
+
+		if(p.size() > 2)
+		{
+			if(reverse)
+			{
+				for(unsigned int j = p.size()-1; j > 1; j--)MakeLoop(p[j], p[j-1], p[j-2], radius);
+				MakeLoop(p[1], p[0], p[p.size()-1], radius);
+				MakeLoop(p[0], p[p.size()-1], p[p.size()-2], radius);
+			}
+			else
+			{
+				MakeLoop(p[p.size()-2], p[p.size()-1], p[0], radius);
+				MakeLoop(p[p.size()-1], p[0], p[1], radius);
+				for(unsigned int j = 2; j < p.size(); j++)MakeLoop(p[j-2], p[j-1], p[j], radius);
+			}
+
+			TPolygon loopy_polygon;
+			loopy_polygon.reserve(pts_for_AddVertex.size());
+			for(std::list<DoubleAreaPoint>::iterator It = pts_for_AddVertex.begin(); It != pts_for_AddVertex.end(); It++)
+			{
+				loopy_polygon.push_back(It->int_point());
+			}
+			c.AddPath(loopy_polygon, ptSubject, true);
+			pts_for_AddVertex.clear();
+		}
+	}
+
+	//c.ForceOrientation(false);
+	c.Execute(ctUnion, pp_new, pftNonZero, pftNonZero);
+
+	if(inwards)
+	{
+		// remove the large square
+		if(pp_new.size() > 0)
+		{
+			pp_new.erase(pp_new.begin());
+		}
+	}
+	else
+	{
+		// reverse all the resulting polygons
+		TPolyPolygon copy = pp_new;
+		pp_new.clear();
+		pp_new.resize(copy.size());
+		for(unsigned int i = 0; i < copy.size(); i++)
+		{
+			const TPolygon& p = copy[i];
+			TPolygon p_new;
+			p_new.resize(p.size());
+			int size_minus_one = p.size() - 1;
+			for(unsigned int j = 0; j < p.size(); j++)p_new[j] = p[size_minus_one - j];
+			pp_new[i] = p_new;
+		}
+	}
+}
+
+static void MakeObround(const Point &pt0, const CVertex &vt1, double radius)
+{
+	Span span(pt0, vt1);
+	Point forward0 = span.GetVector(0.0);
+	Point forward1 = span.GetVector(1.0);
+	Point right0(forward0.y, -forward0.x);
+	Point right1(forward1.y, -forward1.x);
+	right0.normalize();
+	right1.normalize();
+
+	CVertex v0(pt0 + right0 * radius);
+	CVertex v1(vt1.m_type, vt1.m_p + right1 * radius, vt1.m_c);
+	CVertex v2(1, vt1.m_p + right1 * -radius, vt1.m_p);
+	CVertex v3(-vt1.m_type, pt0 + right0 * -radius, vt1.m_c);
+	CVertex v4(1, pt0 + right0 * radius, pt0);
+
+	double save_units = CArea::m_units;
+	CArea::m_units = 1.0;
+
+	AddVertex(v0, NULL);
+	AddVertex(v1, &v0);
+	AddVertex(v2, &v1);
+	AddVertex(v3, &v2);
+	AddVertex(v4, &v3);
+
+	CArea::m_units = save_units;
+}
+
+static void OffsetSpansWithObrounds(const CArea& area, TPolyPolygon &pp_new, double radius)
+{
+	Clipper c;
+
+
+	for(std::list<CCurve>::const_iterator It = area.m_curves.begin(); It != area.m_curves.end(); It++)
+	{
+		pts_for_AddVertex.clear();
+		const CCurve& curve = *It;
+		const CVertex* prev_vertex = NULL;
+		for(std::list<CVertex>::const_iterator It2 = curve.m_vertices.begin(); It2 != curve.m_vertices.end(); It2++)
+		{
+			const CVertex& vertex = *It2;
+			if(prev_vertex)
+			{
+				MakeObround(prev_vertex->m_p, vertex, radius);
+
+				TPolygon loopy_polygon;
+				loopy_polygon.reserve(pts_for_AddVertex.size());
+				for(std::list<DoubleAreaPoint>::iterator It = pts_for_AddVertex.begin(); It != pts_for_AddVertex.end(); It++)
+				{
+					loopy_polygon.push_back(It->int_point());
+				}
+				c.AddPath(loopy_polygon, ptSubject, true);
+				pts_for_AddVertex.clear();
+			}
+			prev_vertex = &vertex;
+		}
+	}
+
+	pp_new.clear();
+	c.Execute(ctUnion, pp_new, pftNonZero, pftNonZero);
+
+	// reverse all the resulting polygons
+	TPolyPolygon copy = pp_new;
+	pp_new.clear();
+	pp_new.resize(copy.size());
+	for(unsigned int i = 0; i < copy.size(); i++)
+	{
+		const TPolygon& p = copy[i];
+		TPolygon p_new;
+		p_new.resize(p.size());
+		int size_minus_one = p.size() - 1;
+		for(unsigned int j = 0; j < p.size(); j++)p_new[j] = p[size_minus_one - j];
+		pp_new[i] = p_new;
+	}
+}
+
+static void MakePolyPoly( const CArea& area, TPolyPolygon &pp, bool reverse = true ){
+	pp.clear();
+
+	for(std::list<CCurve>::const_iterator It = area.m_curves.begin(); It != area.m_curves.end(); It++)
+	{
+		pts_for_AddVertex.clear();
+		const CCurve& curve = *It;
+		const CVertex* prev_vertex = NULL;
+		for(std::list<CVertex>::const_iterator It2 = curve.m_vertices.begin(); It2 != curve.m_vertices.end(); It2++)
+		{
+			const CVertex& vertex = *It2;
+			if(prev_vertex)AddVertex(vertex, prev_vertex);
+			prev_vertex = &vertex;
+		}
+
+		TPolygon p;
+		p.resize(pts_for_AddVertex.size());
+		if(reverse)
+		{
+			unsigned int i = pts_for_AddVertex.size() - 1;// clipper wants them the opposite way to CArea
+			for(std::list<DoubleAreaPoint>::iterator It = pts_for_AddVertex.begin(); It != pts_for_AddVertex.end(); It++, i--)
+			{
+				p[i] = It->int_point();
+			}
+		}
+		else
+		{
+			unsigned int i = 0;
+			for(std::list<DoubleAreaPoint>::iterator It = pts_for_AddVertex.begin(); It != pts_for_AddVertex.end(); It++, i++)
+			{
+				p[i] = It->int_point();
+			}
+		}
+
+		pp.push_back(p);
+	}
+}
+
+static void SetFromResult( CCurve& curve, const TPolygon& p, bool reverse = true )
+{
+	for(unsigned int j = 0; j < p.size(); j++)
+	{
+		const IntPoint &pt = p[j];
+		DoubleAreaPoint dp(pt);
+		CVertex vertex(0, Point(dp.X / CArea::m_units, dp.Y / CArea::m_units), Point(0.0, 0.0));
+		if(reverse)curve.m_vertices.push_front(vertex);
+		else curve.m_vertices.push_back(vertex);
+	}
+	// make a copy of the first point at the end
+	if(reverse)curve.m_vertices.push_front(curve.m_vertices.back());
+	else curve.m_vertices.push_back(curve.m_vertices.front());
+
+	if(CArea::m_fit_arcs)curve.FitArcs();
+}
+
+static void SetFromResult( CArea& area, const TPolyPolygon& pp, bool reverse = true )
+{
+	// delete existing geometry
+	area.m_curves.clear();
+
+	for(unsigned int i = 0; i < pp.size(); i++)
+	{
+		const TPolygon& p = pp[i];
+
+		area.m_curves.push_back(CCurve());
+		CCurve &curve = area.m_curves.back();
+		SetFromResult(curve, p, reverse);
+    }
+}
+
+void CArea::Subtract(const CArea& a2)
+{
+	Clipper c;
+	TPolyPolygon pp1, pp2;
+	MakePolyPoly(*this, pp1);
+	MakePolyPoly(a2, pp2);
+	c.AddPaths(pp1, ptSubject, true);
+	c.AddPaths(pp2, ptClip, true);
+	TPolyPolygon solution;
+	c.Execute(ctDifference, solution);
+	SetFromResult(*this, solution);
+}
+
+void CArea::Intersect(const CArea& a2)
+{
+	Clipper c;
+	TPolyPolygon pp1, pp2;
+	MakePolyPoly(*this, pp1);
+	MakePolyPoly(a2, pp2);
+	c.AddPaths(pp1, ptSubject, true);
+	c.AddPaths(pp2, ptClip, true);
+	TPolyPolygon solution;
+	c.Execute(ctIntersection, solution);
+	SetFromResult(*this, solution);
+}
+
+void CArea::Union(const CArea& a2)
+{
+	Clipper c;
+	TPolyPolygon pp1, pp2;
+	MakePolyPoly(*this, pp1);
+	MakePolyPoly(a2, pp2);
+	c.AddPaths(pp1, ptSubject, true);
+	c.AddPaths(pp2, ptClip, true);
+	TPolyPolygon solution;
+	c.Execute(ctUnion, solution);
+	SetFromResult(*this, solution);
+}
+
+void CArea::Xor(const CArea& a2)
+{
+	Clipper c;
+	TPolyPolygon pp1, pp2;
+	MakePolyPoly(*this, pp1);
+	MakePolyPoly(a2, pp2);
+	c.AddPaths(pp1, ptSubject, true);
+	c.AddPaths(pp2, ptClip, true);
+	TPolyPolygon solution;
+	c.Execute(ctXor, solution);
+	SetFromResult(*this, solution);
+}
+
+void CArea::Offset(double inwards_value)
+{
+	TPolyPolygon pp, pp2;
+	MakePolyPoly(*this, pp, false);
+	OffsetWithLoops(pp, pp2, inwards_value * m_units);
+	SetFromResult(*this, pp2, false);
+	this->Reorder();
+}
+
+void CArea::Thicken(double value)
+{
+	TPolyPolygon pp;
+	OffsetSpansWithObrounds(*this, pp, value * m_units);
+	SetFromResult(*this, pp, false);
+	this->Reorder();
+}
+
+void UnFitArcs(CCurve &curve)
+{
+	pts_for_AddVertex.clear();
+	const CVertex* prev_vertex = NULL;
+	for(std::list<CVertex>::const_iterator It2 = curve.m_vertices.begin(); It2 != curve.m_vertices.end(); It2++)
+	{
+		const CVertex& vertex = *It2;
+		AddVertex(vertex, prev_vertex);
+		prev_vertex = &vertex;
+	}
+
+	curve.m_vertices.clear();
+
+	for(std::list<DoubleAreaPoint>::iterator It = pts_for_AddVertex.begin(); It != pts_for_AddVertex.end(); It++)
+	{
+		DoubleAreaPoint &pt = *It;
+		CVertex vertex(0, Point(pt.X / CArea::m_units, pt.Y / CArea::m_units), Point(0.0, 0.0));
+		curve.m_vertices.push_back(vertex);
+	}
+}

src/Mod/Path/libarea/AreaDxf.cpp

@@ -0,0 +1,31 @@
+// AreaDxf.cpp
+// Copyright (c) 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#include "AreaDxf.h"
+#include "Area.h"
+
+AreaDxfRead::AreaDxfRead(CArea* area, const char* filepath):CDxfRead(filepath), m_area(area){}
+
+void AreaDxfRead::StartCurveIfNecessary(const double* s)
+{
+	Point ps(s);
+	if((m_area->m_curves.size() == 0) || (m_area->m_curves.back().m_vertices.size() == 0) || (m_area->m_curves.back().m_vertices.back().m_p != ps))
+	{
+		// start a new curve
+		m_area->m_curves.push_back(CCurve());
+		m_area->m_curves.back().m_vertices.push_back(ps);
+	}
+}
+
+void AreaDxfRead::OnReadLine(const double* s, const double* e)
+{
+	StartCurveIfNecessary(s);
+	m_area->m_curves.back().m_vertices.push_back(Point(e));
+}
+
+void AreaDxfRead::OnReadArc(const double* s, const double* e, const double* c, bool dir)
+{
+	StartCurveIfNecessary(s);
+	m_area->m_curves.back().m_vertices.push_back(CVertex(dir?1:0, Point(e), Point(c)));
+}

src/Mod/Path/libarea/AreaDxf.h

@@ -0,0 +1,23 @@
+// AreaDxf.h
+// Copyright (c) 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#pragma once
+
+#include "dxf.h"
+
+class CSketch;
+class CArea;
+class CCurve;
+
+class AreaDxfRead : public CDxfRead{
+	void StartCurveIfNecessary(const double* s);
+
+public:
+	CArea* m_area;
+	AreaDxfRead(CArea* area, const char* filepath);
+
+	// AreaDxfRead's virtual functions
+	void OnReadLine(const double* s, const double* e);
+	void OnReadArc(const double* s, const double* e, const double* c, bool dir);
+};

src/Mod/Path/libarea/AreaOrderer.cpp

@@ -0,0 +1,153 @@
+// AreaOrderer.cpp
+// Copyright (c) 2010, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#include "AreaOrderer.h"
+#include "Area.h"
+
+CAreaOrderer* CInnerCurves::area_orderer = NULL;
+
+CInnerCurves::CInnerCurves(CInnerCurves* pOuter, const CCurve* curve)
+{
+	m_pOuter = pOuter;
+	m_curve = curve;
+	m_unite_area = NULL;
+}
+
+CInnerCurves::~CInnerCurves()
+{
+	delete m_unite_area;
+}
+
+void CInnerCurves::Insert(const CCurve* pcurve)
+{
+	std::list<CInnerCurves*> outside_of_these;
+	std::list<CInnerCurves*> crossing_these;
+
+	// check all inner curves
+	for(std::set<CInnerCurves*>::iterator It = m_inner_curves.begin(); It != m_inner_curves.end(); It++)
+	{
+		CInnerCurves* c = *It;
+
+		switch(GetOverlapType(*pcurve, *(c->m_curve)))
+		{
+		case eOutside:
+			outside_of_these.push_back(c);
+			break;
+
+		case eInside:
+			// insert in this inner curve
+			c->Insert(pcurve);
+			return;
+
+		case eSiblings:
+			break;
+
+		case eCrossing:
+			crossing_these.push_back(c);
+			break;
+		}
+	}
+
+	// add as a new inner
+	CInnerCurves* new_item = new CInnerCurves(this, pcurve);
+	this->m_inner_curves.insert(new_item);
+
+	for(std::list<CInnerCurves*>::iterator It = outside_of_these.begin(); It != outside_of_these.end(); It++)
+	{
+		// move items
+		CInnerCurves* c = *It;
+		c->m_pOuter = new_item;
+		new_item->m_inner_curves.insert(c);
+		this->m_inner_curves.erase(c);
+	}
+
+	for(std::list<CInnerCurves*>::iterator It = crossing_these.begin(); It != crossing_these.end(); It++)
+	{
+		// unite these
+		CInnerCurves* c = *It;
+		new_item->Unite(c);
+		this->m_inner_curves.erase(c);
+	}
+}
+
+void CInnerCurves::GetArea(CArea &area, bool outside, bool use_curve)const
+{
+	if(use_curve && m_curve)
+	{
+		area.m_curves.push_back(*m_curve);
+		outside = !outside;
+	}
+
+	std::list<const CInnerCurves*> do_after;
+
+	for(std::set<CInnerCurves*>::const_iterator It = m_inner_curves.begin(); It != m_inner_curves.end(); It++)
+	{
+		const CInnerCurves* c = *It;
+		area.m_curves.push_back(*c->m_curve);
+		if(!outside)area.m_curves.back().Reverse();
+
+		if(outside)c->GetArea(area, !outside, false);
+		else do_after.push_back(c);
+	}
+
+	for(std::list<const CInnerCurves*>::iterator It = do_after.begin(); It != do_after.end(); It++)
+	{
+		const CInnerCurves* c = *It;
+		c->GetArea(area, !outside, false);
+	}
+}
+
+void CInnerCurves::Unite(const CInnerCurves* c)
+{
+	// unite all the curves in c, with this one
+	CArea* new_area = new CArea();
+	new_area->m_curves.push_back(*m_curve);
+	delete m_unite_area;
+	m_unite_area = new_area;
+
+	CArea a2;
+	c->GetArea(a2);
+
+	m_unite_area->Union(a2);
+	m_unite_area->Reorder();
+	for(std::list<CCurve>::iterator It = m_unite_area->m_curves.begin(); It != m_unite_area->m_curves.end(); It++)
+	{
+		CCurve &curve = *It;
+		if(It == m_unite_area->m_curves.begin())
+			m_curve = &curve;
+		else
+		{
+			if(curve.IsClockwise())curve.Reverse();
+			Insert(&curve);
+		}
+	}
+}
+
+CAreaOrderer::CAreaOrderer()
+{
+	m_top_level = new CInnerCurves(NULL, NULL);
+}
+
+void CAreaOrderer::Insert(CCurve* pcurve)
+{
+	CInnerCurves::area_orderer = this;
+
+	// make them all anti-clockwise as they come in
+	if(pcurve->IsClockwise())pcurve->Reverse();
+
+	m_top_level->Insert(pcurve);
+}
+
+CArea CAreaOrderer::ResultArea()const
+{
+	CArea a;
+
+	if(m_top_level)
+	{
+		m_top_level->GetArea(a);
+	}
+
+	return a;
+}
+

src/Mod/Path/libarea/AreaOrderer.h

@@ -0,0 +1,40 @@
+// AreaOrderer.h
+// Copyright (c) 2010, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#pragma once
+#include <list>
+#include <set>
+
+class CArea;
+class CCurve;
+
+class CAreaOrderer;
+
+class CInnerCurves
+{
+	CInnerCurves* m_pOuter;
+	const CCurve* m_curve; // always empty if top level
+	std::set<CInnerCurves*> m_inner_curves;
+	CArea *m_unite_area; // new curves made by uniting are stored here
+
+public:
+	static CAreaOrderer* area_orderer;
+	CInnerCurves(CInnerCurves* pOuter, const CCurve* curve);
+	~CInnerCurves();
+
+	void Insert(const CCurve* pcurve);
+	void GetArea(CArea &area, bool outside = true, bool use_curve = true)const;
+	void Unite(const CInnerCurves* c);
+};
+
+class CAreaOrderer
+{
+public:
+	CInnerCurves* m_top_level;
+
+	CAreaOrderer();
+
+	void Insert(CCurve* pcurve);
+	CArea ResultArea()const;
+};

src/Mod/Path/libarea/AreaPocket.cpp

@@ -0,0 +1,565 @@
+// AreaPocket.cpp
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+// implements CArea::MakeOnePocketCurve
+
+#include "Area.h"
+
+#include <map>
+#include <set>
+
+static const CAreaPocketParams* pocket_params = NULL;
+
+class IslandAndOffset
+{
+public:
+	const CCurve* island;
+	CArea offset;
+	std::list<CCurve> island_inners;
+	std::list<IslandAndOffset*> touching_offsets;
+
+	IslandAndOffset(const CCurve* Island)
+	{
+		island = Island;
+
+		offset.m_curves.push_back(*island);
+		offset.m_curves.back().Reverse();
+
+		offset.Offset(-pocket_params->stepover);
+
+
+		if(offset.m_curves.size() > 1)
+		{
+			for(std::list<CCurve>::iterator It = offset.m_curves.begin(); It != offset.m_curves.end(); It++)
+			{
+				if(It == offset.m_curves.begin())continue;
+				island_inners.push_back(*It);
+				island_inners.back().Reverse();
+			}
+			offset.m_curves.resize(1);
+		}
+	}
+};
+
+class CurveTree
+{
+	static std::list<CurveTree*> to_do_list_for_MakeOffsets;
+	void MakeOffsets2();
+	static std::list<CurveTree*> islands_added;
+
+public:
+	Point point_on_parent;
+	CCurve curve;
+	std::list<CurveTree*> inners;
+	std::list<const IslandAndOffset*> offset_islands;
+	CurveTree(const CCurve &c)
+	{
+		curve = c;
+	}
+	~CurveTree(){}
+
+	void MakeOffsets();
+};
+std::list<CurveTree*> CurveTree::islands_added;
+
+class GetCurveItem
+{
+public:
+	CurveTree* curve_tree;
+	std::list<CVertex>::iterator EndIt;
+	static std::list<GetCurveItem> to_do_list;
+
+	GetCurveItem(CurveTree* ct, std::list<CVertex>::iterator EIt):curve_tree(ct), EndIt(EIt){}
+
+	void GetCurve(CCurve& output);
+	CVertex& back(){std::list<CVertex>::iterator It = EndIt; It--; return *It;}
+};
+
+std::list<GetCurveItem> GetCurveItem::to_do_list;
+std::list<CurveTree*> CurveTree::to_do_list_for_MakeOffsets;
+
+void GetCurveItem::GetCurve(CCurve& output)
+{
+	// walk around the curve adding spans to output until we get to an inner's point_on_parent
+	// then add a line from the inner's point_on_parent to inner's start point, then GetCurve from inner
+
+	// add start point
+	if(CArea::m_please_abort)return;
+	output.m_vertices.insert(this->EndIt, CVertex(curve_tree->curve.m_vertices.front()));
+
+	std::list<CurveTree*> inners_to_visit;
+	for(std::list<CurveTree*>::iterator It2 = curve_tree->inners.begin(); It2 != curve_tree->inners.end(); It2++)
+	{
+		inners_to_visit.push_back(*It2);
+	}
+
+	const CVertex* prev_vertex = NULL;
+
+	for(std::list<CVertex>::iterator It = curve_tree->curve.m_vertices.begin(); It != curve_tree->curve.m_vertices.end(); It++)
+	{
+		const CVertex& vertex = *It;
+		if(prev_vertex)
+		{
+			Span span(prev_vertex->m_p, vertex);
+
+			// order inners on this span
+			std::multimap<double, CurveTree*> ordered_inners;
+			for(std::list<CurveTree*>::iterator It2 = inners_to_visit.begin(); It2 != inners_to_visit.end();)
+			{
+				CurveTree *inner = *It2;
+				double t;
+				if(span.On(inner->point_on_parent, &t))
+				{
+					ordered_inners.insert(std::make_pair(t, inner));
+					It2 = inners_to_visit.erase(It2);
+				}
+				else
+				{
+					It2++;
+				}
+				if(CArea::m_please_abort)return;
+			}
+
+			if(CArea::m_please_abort)return;
+			for(std::multimap<double, CurveTree*>::iterator It2 = ordered_inners.begin(); It2 != ordered_inners.end(); It2++)
+			{
+				CurveTree& inner = *(It2->second);
+				if(inner.point_on_parent.dist(back().m_p) > 0.01/CArea::m_units)
+				{
+					output.m_vertices.insert(this->EndIt, CVertex(vertex.m_type, inner.point_on_parent, vertex.m_c));
+				}
+				if(CArea::m_please_abort)return;
+
+				// vertex add after GetCurve
+				std::list<CVertex>::iterator VIt = output.m_vertices.insert(this->EndIt, CVertex(inner.point_on_parent));
+
+				//inner.GetCurve(output);
+				GetCurveItem::to_do_list.push_back(GetCurveItem(&inner, VIt));
+			}
+
+			if(back().m_p != vertex.m_p)output.m_vertices.insert(this->EndIt, vertex);
+		}
+		prev_vertex = &vertex;
+	}
+
+	if(CArea::m_please_abort)return;
+	for(std::list<CurveTree*>::iterator It2 = inners_to_visit.begin(); It2 != inners_to_visit.end(); It2++)
+	{
+		CurveTree &inner = *(*It2);
+		if(inner.point_on_parent != back().m_p)
+		{
+			output.m_vertices.insert(this->EndIt, CVertex(inner.point_on_parent));
+		}
+		if(CArea::m_please_abort)return;
+
+		// vertex add after GetCurve
+		std::list<CVertex>::iterator VIt = output.m_vertices.insert(this->EndIt, CVertex(inner.point_on_parent));
+
+		//inner.GetCurve(output);
+		GetCurveItem::to_do_list.push_back(GetCurveItem(&inner, VIt));
+
+	}
+}
+
+class IslandAndOffsetLink
+{
+public:
+	const IslandAndOffset* island_and_offset;
+	CurveTree* add_to;
+	IslandAndOffsetLink(const IslandAndOffset* i, CurveTree* a){island_and_offset = i; add_to = a;}
+};
+
+static Point GetNearestPoint(CurveTree* curve_tree, std::list<CurveTree*> &islands_added, const CCurve &test_curve, CurveTree** best_curve_tree)
+{
+	// find nearest point to test_curve, from curve and all the islands in 
+	double best_dist;
+	Point best_point = curve_tree->curve.NearestPoint(test_curve, &best_dist);
+	*best_curve_tree = curve_tree;
+	for(std::list<CurveTree*>::iterator It = islands_added.begin(); It != islands_added.end(); It++)
+	{
+		CurveTree* island = *It;
+		double dist;
+		Point p = island->curve.NearestPoint(test_curve, &dist);
+		if(dist < best_dist)
+		{
+			*best_curve_tree = island;
+			best_point = p;
+			best_dist = dist;
+		}
+	}
+
+	return best_point;
+}
+
+void CurveTree::MakeOffsets2()
+{
+	// make offsets
+
+	if(CArea::m_please_abort)return;
+	CArea smaller;
+	smaller.m_curves.push_back(curve);
+	smaller.Offset(pocket_params->stepover);
+
+	if(CArea::m_please_abort)return;
+
+	// test islands
+	for(std::list<const IslandAndOffset*>::iterator It = offset_islands.begin(); It != offset_islands.end();)
+	{
+		const IslandAndOffset* island_and_offset = *It;
+
+		if(GetOverlapType(island_and_offset->offset, smaller) == eInside)
+			It++; // island is still inside
+		else
+		{
+			inners.push_back(new CurveTree(*island_and_offset->island));
+			islands_added.push_back(inners.back());
+			inners.back()->point_on_parent = curve.NearestPoint(*island_and_offset->island);
+			if(CArea::m_please_abort)return;
+			Point island_point = island_and_offset->island->NearestPoint(inners.back()->point_on_parent);
+			if(CArea::m_please_abort)return;
+			inners.back()->curve.ChangeStart(island_point);
+			if(CArea::m_please_abort)return;
+
+			// add the island offset's inner curves
+			for(std::list<CCurve>::const_iterator It2 = island_and_offset->island_inners.begin(); It2 != island_and_offset->island_inners.end(); It2++)
+			{
+				const CCurve& island_inner = *It2;
+				inners.back()->inners.push_back(new CurveTree(island_inner));
+				inners.back()->inners.back()->point_on_parent = inners.back()->curve.NearestPoint(island_inner);
+				if(CArea::m_please_abort)return;
+				Point island_point = island_inner.NearestPoint(inners.back()->inners.back()->point_on_parent);
+				if(CArea::m_please_abort)return;
+				inners.back()->inners.back()->curve.ChangeStart(island_point);
+				to_do_list_for_MakeOffsets.push_back(inners.back()->inners.back()); // do it later, in a while loop
+				if(CArea::m_please_abort)return;
+			}
+
+			smaller.Subtract(island_and_offset->offset);
+
+			std::set<const IslandAndOffset*> added;
+
+			std::list<IslandAndOffsetLink> touching_list;
+			for(std::list<IslandAndOffset*>::const_iterator It2 = island_and_offset->touching_offsets.begin(); It2 != island_and_offset->touching_offsets.end(); It2++)
+			{
+				const IslandAndOffset* touching = *It2;
+				touching_list.push_back(IslandAndOffsetLink(touching, inners.back()));
+				added.insert(touching);
+			}
+
+			while(touching_list.size() > 0)
+			{
+				IslandAndOffsetLink touching = touching_list.front();
+				touching_list.pop_front();
+				touching.add_to->inners.push_back(new CurveTree(*touching.island_and_offset->island));
+				islands_added.push_back(touching.add_to->inners.back());
+				touching.add_to->inners.back()->point_on_parent = touching.add_to->curve.NearestPoint(*touching.island_and_offset->island);
+				Point island_point = touching.island_and_offset->island->NearestPoint(touching.add_to->inners.back()->point_on_parent);
+				touching.add_to->inners.back()->curve.ChangeStart(island_point);
+				smaller.Subtract(touching.island_and_offset->offset);
+
+				// add the island offset's inner curves
+				for(std::list<CCurve>::const_iterator It2 = touching.island_and_offset->island_inners.begin(); It2 != touching.island_and_offset->island_inners.end(); It2++)
+				{
+					const CCurve& island_inner = *It2;
+					touching.add_to->inners.back()->inners.push_back(new CurveTree(island_inner));
+					touching.add_to->inners.back()->inners.back()->point_on_parent = touching.add_to->inners.back()->curve.NearestPoint(island_inner);
+					if(CArea::m_please_abort)return;
+					Point island_point = island_inner.NearestPoint(touching.add_to->inners.back()->inners.back()->point_on_parent);
+					if(CArea::m_please_abort)return;
+					touching.add_to->inners.back()->inners.back()->curve.ChangeStart(island_point);
+					to_do_list_for_MakeOffsets.push_back(touching.add_to->inners.back()->inners.back()); // do it later, in a while loop
+					if(CArea::m_please_abort)return;
+				}
+
+				for(std::list<IslandAndOffset*>::const_iterator It2 = touching.island_and_offset->touching_offsets.begin(); It2 != touching.island_and_offset->touching_offsets.end(); It2++)
+				{
+					if(added.find(*It2)==added.end() && ((*It2) != island_and_offset))
+					{
+						touching_list.push_back(IslandAndOffsetLink(*It2, touching.add_to->inners.back()));
+						added.insert(*It2);
+					}
+				}
+			}
+
+			if(CArea::m_please_abort)return;
+			It = offset_islands.erase(It);
+
+			for(std::set<const IslandAndOffset*>::iterator It2 = added.begin(); It2 != added.end(); It2++)
+			{
+				const IslandAndOffset* i = *It2;
+				offset_islands.remove(i);
+			}
+
+			if(offset_islands.size() == 0)break;
+			It = offset_islands.begin();
+		}
+	}
+
+	CArea::m_processing_done += CArea::m_MakeOffsets_increment;
+	if(CArea::m_processing_done > CArea::m_after_MakeOffsets_length)CArea::m_processing_done = CArea::m_after_MakeOffsets_length;
+
+	std::list<CArea> separate_areas;
+	smaller.Split(separate_areas);
+	if(CArea::m_please_abort)return;
+	for(std::list<CArea>::iterator It = separate_areas.begin(); It != separate_areas.end(); It++)
+	{
+		CArea& separate_area = *It;
+		CCurve& first_curve = separate_area.m_curves.front();
+
+		CurveTree* nearest_curve_tree = NULL;
+		Point near_point = GetNearestPoint(this, islands_added, first_curve, &nearest_curve_tree);
+
+		nearest_curve_tree->inners.push_back(new CurveTree(first_curve));
+
+		for(std::list<const IslandAndOffset*>::iterator It = offset_islands.begin(); It != offset_islands.end();It++)
+		{
+			const IslandAndOffset* island_and_offset = *It;
+			if(GetOverlapType(island_and_offset->offset, separate_area) == eInside)
+				nearest_curve_tree->inners.back()->offset_islands.push_back(island_and_offset);
+			if(CArea::m_please_abort)return;
+		}
+
+		nearest_curve_tree->inners.back()->point_on_parent = near_point;
+
+		if(CArea::m_please_abort)return;
+		Point first_curve_point = first_curve.NearestPoint(nearest_curve_tree->inners.back()->point_on_parent);
+		if(CArea::m_please_abort)return;
+		nearest_curve_tree->inners.back()->curve.ChangeStart(first_curve_point);
+		if(CArea::m_please_abort)return;
+		to_do_list_for_MakeOffsets.push_back(nearest_curve_tree->inners.back()); // do it later, in a while loop
+		if(CArea::m_please_abort)return;
+	}
+}
+
+void CurveTree::MakeOffsets()
+{
+	to_do_list_for_MakeOffsets.push_back(this);
+	islands_added.clear();
+
+	while(to_do_list_for_MakeOffsets.size() > 0)
+	{
+		CurveTree* curve_tree = to_do_list_for_MakeOffsets.front();
+		to_do_list_for_MakeOffsets.pop_front();
+		curve_tree->MakeOffsets2();
+	}
+}
+
+void recur(std::list<CArea> &arealist, const CArea& a1, const CAreaPocketParams &params, int level)
+{
+	//if(level > 3)return;
+
+    // this makes arealist by recursively offsetting a1 inwards
+    
+	if(a1.m_curves.size() == 0)
+		return;
+    
+	if(params.from_center)
+		arealist.push_front(a1);
+	else
+		arealist.push_back(a1);
+
+    CArea a_offset = a1;
+    a_offset.Offset(params.stepover);
+    
+    // split curves into new areas
+	if(CArea::HolesLinked())
+	{
+		for(std::list<CCurve>::iterator It = a_offset.m_curves.begin(); It != a_offset.m_curves.end(); It++)
+		{
+            CArea a2;
+			a2.m_curves.push_back(*It);
+            recur(arealist, a2, params, level + 1);
+		}
+	}
+    else
+	{
+        // split curves into new areas
+		a_offset.Reorder();
+        CArea* a2 = NULL;
+       
+		for(std::list<CCurve>::iterator It = a_offset.m_curves.begin(); It != a_offset.m_curves.end(); It++)
+		{
+			CCurve& curve = *It;
+			if(curve.IsClockwise())
+			{
+				if(a2 != NULL)
+					a2->m_curves.push_back(curve);
+			}
+			else
+			{
+				if(a2 != NULL)
+					recur(arealist, *a2, params, level + 1);
+				else
+					a2 = new CArea();
+                a2->m_curves.push_back(curve);
+			}
+		}
+
+		if(a2 != NULL)
+			recur(arealist, *a2, params, level + 1);
+	}
+}
+
+static CArea make_obround(const Point& p0, const Point& p1, double radius)
+{
+    Point dir = p1 - p0;
+    double d = dir.length();
+    dir.normalize();
+    Point right(dir.y, -dir.x);
+    CArea obround;
+    CCurve c;
+	if(fabs(radius) < 0.0000001)radius = (radius > 0.0) ? 0.002 : (-0.002);
+    Point vt0 = p0 + right * radius;
+    Point vt1 = p1 + right * radius;
+    Point vt2 = p1 - right * radius;
+    Point vt3 = p0 - right * radius;
+	c.append(vt0);
+    c.append(vt1);
+    c.append(CVertex(1, vt2, p1));
+    c.append(vt3);
+    c.append(CVertex(1, vt0, p0));
+    obround.append(c);
+    return obround;
+}
+    
+static bool feed_possible(const CArea &area_for_feed_possible, const Point& p0, const Point& p1, double tool_radius)
+{
+    CArea obround = make_obround(p0, p1, tool_radius);
+    obround.Subtract(area_for_feed_possible);
+	return obround.m_curves.size() == 0;
+}
+
+void MarkOverlappingOffsetIslands(std::list<IslandAndOffset> &offset_islands)
+{
+	for(std::list<IslandAndOffset>::iterator It1 = offset_islands.begin(); It1 != offset_islands.end(); It1++)
+	{
+		std::list<IslandAndOffset>::iterator It2 = It1;
+		It2++;
+		for(;It2 != offset_islands.end(); It2++)
+		{
+			IslandAndOffset &o1 = *It1;
+			IslandAndOffset &o2 = *It2;
+
+			if(GetOverlapType(o1.offset, o2.offset) == eCrossing)
+			{
+				o1.touching_offsets.push_back(&o2);
+				o2.touching_offsets.push_back(&o1);
+			}
+		}
+	}
+}
+
+void CArea::MakeOnePocketCurve(std::list<CCurve> &curve_list, const CAreaPocketParams &params)const
+{
+	if(CArea::m_please_abort)return;
+#if 0  // simple offsets with feed or rapid joins
+	CArea area_for_feed_possible = *this;
+
+	area_for_feed_possible.Offset(-params.tool_radius - 0.01);
+	CArea a_offset = *this;
+
+	std::list<CArea> arealist;
+	recur(arealist, a_offset, params, 0);
+
+	bool first = true;
+
+	for(std::list<CArea>::iterator It = arealist.begin(); It != arealist.end(); It++)
+	{
+		CArea& area = *It;
+		for(std::list<CCurve>::iterator It = area.m_curves.begin(); It != area.m_curves.end(); It++)
+		{
+			CCurve& curve = *It;
+			if(!first)
+			{
+				// try to join these curves with a feed move, if possible and not too long
+				CCurve &prev_curve = curve_list.back();
+				const Point &prev_p = prev_curve.m_vertices.back().m_p;
+				const Point &next_p = curve.m_vertices.front().m_p;
+
+				if(feed_possible(area_for_feed_possible, prev_p, next_p, params.tool_radius))
+				{
+					// join curves
+					prev_curve += curve;
+				}
+				else
+				{
+					curve_list.push_back(curve);
+				}
+			}
+			else
+			{
+				curve_list.push_back(curve);
+			}
+			first = false;
+		}
+	}
+#else
+	pocket_params = &params;
+	if(m_curves.size() == 0)
+	{
+		CArea::m_processing_done += CArea::m_single_area_processing_length;
+		return;
+	}
+	CurveTree top_level(m_curves.front());
+
+	std::list<IslandAndOffset> offset_islands;
+
+	for(std::list<CCurve>::const_iterator It = m_curves.begin(); It != m_curves.end(); It++)
+	{
+		const CCurve& c = *It;
+		if(It != m_curves.begin())
+		{
+			IslandAndOffset island_and_offset(&c);
+			offset_islands.push_back(island_and_offset);
+			top_level.offset_islands.push_back(&(offset_islands.back()));
+			if(m_please_abort)return;
+		}
+	}
+
+	MarkOverlappingOffsetIslands(offset_islands);
+
+	CArea::m_processing_done += CArea::m_single_area_processing_length * 0.1;
+
+	double MakeOffsets_processing_length = CArea::m_single_area_processing_length * 0.8;
+	CArea::m_after_MakeOffsets_length = CArea::m_processing_done + MakeOffsets_processing_length;
+	double guess_num_offsets = sqrt(GetArea(true)) * 0.5 / params.stepover;
+	CArea::m_MakeOffsets_increment = MakeOffsets_processing_length / guess_num_offsets;
+
+	top_level.MakeOffsets();
+	if(CArea::m_please_abort)return;
+	CArea::m_processing_done = CArea::m_after_MakeOffsets_length;
+
+	curve_list.push_back(CCurve());
+	CCurve& output = curve_list.back();
+
+	GetCurveItem::to_do_list.push_back(GetCurveItem(&top_level, output.m_vertices.end()));
+
+	while(GetCurveItem::to_do_list.size() > 0)
+	{
+		GetCurveItem item = GetCurveItem::to_do_list.front();
+		item.GetCurve(output);
+		GetCurveItem::to_do_list.pop_front();
+	}
+
+	// delete curve_trees non-recursively
+	std::list<CurveTree*> CurveTreeDestructList;
+	for(std::list<CurveTree*>::iterator It = top_level.inners.begin(); It != top_level.inners.end(); It++)
+	{
+		CurveTreeDestructList.push_back(*It);
+	}
+	while(CurveTreeDestructList.size() > 0)
+	{
+		CurveTree* curve_tree = CurveTreeDestructList.front();
+		CurveTreeDestructList.pop_front();
+		for(std::list<CurveTree*>::iterator It = curve_tree->inners.begin(); It != curve_tree->inners.end(); It++)
+		{
+			CurveTreeDestructList.push_back(*It);
+		}
+		delete curve_tree;
+	}
+
+	CArea::m_processing_done += CArea::m_single_area_processing_length * 0.1;
+#endif
+}
+

src/Mod/Path/libarea/Box2D.h

@@ -0,0 +1,70 @@
+// Box2D.h
+// Copyright (c) 2009, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#pragma once
+
+#include <string.h>	// for memcpy() prototype
+#include <math.h>	// for sqrt() prototype
+
+class CBox2D{
+public:
+	Point m_minxy;
+	Point m_maxxy;
+	bool m_valid;
+
+	CBox2D():m_valid(false){}
+	CBox2D(const Point& minxy, const Point& maxxy):m_minxy(minxy), m_maxxy(maxxy), m_valid(true){}
+
+	bool operator==( const CBox2D & rhs ) const
+	{
+		if(m_minxy != rhs.m_minxy)return false;
+		if(m_maxxy != rhs.m_maxxy)return false;
+		if (m_valid != rhs.m_valid) return(false);
+
+		return(true);
+	}
+
+	bool operator!=( const CBox2D & rhs ) const { return(! (*this == rhs)); }
+
+
+	void Insert(const Point &p){ // insert a point
+		if(m_valid){
+			if(p.x < m_minxy.x)m_minxy.x = p.x;
+			if(p.y < m_minxy.y)m_minxy.y = p.y;
+			if(p.x > m_maxxy.x)m_maxxy.x = p.x;
+			if(p.y > m_maxxy.y)m_maxxy.y = p.y;
+		}
+		else
+		{
+			m_valid = true;
+			m_minxy = p;
+			m_maxxy = p;
+		}
+	}
+
+	void Insert(const CBox2D& b){
+		if(b.m_valid){
+			if(m_valid){
+				if(b.m_minxy.x < m_minxy.x)m_minxy.x = b.m_minxy.x;
+				if(b.m_minxy.y < m_minxy.y)m_minxy.y = b.m_minxy.y;
+				if(b.m_maxxy.x > m_maxxy.x)m_maxxy.x = b.m_maxxy.x;
+				if(b.m_maxxy.y > m_maxxy.y)m_maxxy.y = b.m_maxxy.y;
+			}
+			else{
+				m_valid = b.m_valid;
+				m_minxy = b.m_minxy;
+				m_maxxy = b.m_maxxy;
+			}
+		}
+	}
+	Point Centre() const {return (m_minxy + m_maxxy) * 0.5;}
+	double Width() const {if(m_valid)return m_maxxy.x - m_minxy.x; else return 0.0;}
+	double Height() const {if(m_valid)return m_maxxy.y - m_minxy.y; else return 0.0;}
+	double Radius() const {return sqrt(Width() * Width() + Height() * Height()) /2;}
+	double MinX() const { return(m_minxy.x); }
+	double MaxX() const { return(m_maxxy.x); }
+	double MinY() const { return(m_minxy.y); }
+	double MaxY() const { return(m_maxxy.y); }
+};
+

src/Mod/Path/libarea/Circle.cpp

@@ -0,0 +1,78 @@
+// Circle.cpp
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#include "Circle.h"
+
+Circle::Circle(const Point& p0, const Point& p1, const Point& p2)
+{
+	// from TangentCircles in http://code.google.com/p/heekscad/source/browse/trunk/src/Geom.cpp
+
+	// set default values, in case this fails
+	m_radius = 0.0;
+	m_c = Point(0, 0);
+
+	double x1 = p0.x;
+	double y1 = p0.y;
+	double x2 = p1.x;
+	double y2 = p1.y;
+	double x3 = p2.x;
+	double y3 = p2.y;
+
+	double a = 2 * (x1 - x2);
+	double b = 2 * (y1 - y2);
+	double d = (x1 * x1 + y1 * y1) - (x2 * x2 + y2 * y2);
+
+	double A = 2 * (x1 - x3);
+	double B = 2 * (y1 - y3);
+	double D = (x1 * x1 + y1 * y1) - (x3 * x3 + y3 * y3);
+
+	double aBmbA = (a*B - b*A); // aB - bA
+
+	// x = k + Kr where
+	double k = (B*d - b*D) / aBmbA;
+
+	// y = l + Lr where
+	double l = (-A*d + a*D)/ aBmbA;
+
+	double qa = -1;
+	double qb = 0.0;
+	double qc = k*k + x1*x1 -2*k*x1 + l*l + y1*y1 - 2*l*y1;
+
+	// solve the quadratic equation, r = (-b +- sqrt(b*b - 4*a*c))/(2 * a)
+	for(int qs = 0; qs<2; qs++){
+		double bb = qb*qb;
+		double ac4 = 4*qa*qc;
+		if(ac4 <= bb){
+			double r = (-qb + ((qs == 0) ? 1 : -1) * sqrt(bb - ac4))/(2 * qa);
+			double x = k;
+			double y = l;
+
+			// set the circle
+			if(r >= 0.0){
+				m_c = Point(x, y);
+				m_radius = r;
+			}
+		}
+	}
+}
+
+bool Circle::PointIsOn(const Point& p, double accuracy)
+{
+	double rp = p.dist(m_c);
+	bool on = fabs(m_radius - rp) < accuracy;
+	return on;
+}
+
+bool Circle::LineIsOn(const Point& p0, const Point& p1, double accuracy)
+{
+	// checks the points are on the arc, to the given accuracy, and the mid point of the line.
+
+	if(!PointIsOn(p0, accuracy))return false;
+	if(!PointIsOn(p1, accuracy))return false;
+
+	Point mid = Point((p0 + p1)/2);
+	if(!PointIsOn(mid, accuracy))return false;
+
+	return true;
+}

src/Mod/Path/libarea/Circle.h

@@ -0,0 +1,19 @@
+// Circle.h
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#pragma once
+
+#include "Point.h"
+
+class Circle{
+public:
+	Point m_c;
+	double m_radius;
+
+	Circle(const Point& c, double radius):m_c(c), m_radius(radius){}
+	Circle(const Point& p0, const Point& p1, const Point& p2); // circle through three points
+	
+	bool PointIsOn(const Point& p, double accuracy);
+	bool LineIsOn(const Point& p0, const Point& p1, double accuracy);
+};

src/Mod/Path/libarea/Curve.h

@@ -0,0 +1,105 @@
+// Curve.h
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#pragma once
+
+#include <vector>
+#include <list>
+#include <math.h>
+#include "Point.h"
+#include "Box2D.h"
+
+class Line{
+public:
+	Point p0;
+	Point v;
+
+	// constructors
+	Line(const Point& P0, const Point& V);
+
+	double Dist(const Point& p)const;
+};
+
+class CArc;
+
+class CVertex
+{
+public:
+	int m_type; // 0 - line ( or start point ), 1 - anti-clockwise arc, -1 - clockwise arc
+	Point m_p; // end point
+	Point m_c; // centre point in absolute coordinates
+	int m_user_data;
+
+	CVertex():m_type(0), m_p(Point(0, 0)), m_c(Point(0,0)), m_user_data(0){}
+	CVertex(int type, const Point& p, const Point& c, int user_data = 0);
+	CVertex(const Point& p, int user_data = 0);
+};
+
+class Span
+{
+	Point NearestPointNotOnSpan(const Point& p)const;
+	double Parameter(const Point& p)const;
+	Point NearestPointToSpan(const Span& p, double &d)const;
+
+	static const Point null_point;
+	static const CVertex null_vertex;
+
+public:
+	bool m_start_span;
+	Point m_p;
+	CVertex m_v;
+	Span();
+	Span(const Point& p, const CVertex& v, bool start_span = false):m_start_span(start_span), m_p(p), m_v(v){}
+	Point NearestPoint(const Point& p)const;
+	Point NearestPoint(const Span& p, double *d = NULL)const;
+	void GetBox(CBox2D &box);
+	double IncludedAngle()const;
+	double GetArea()const;
+	bool On(const Point& p, double* t = NULL)const;
+	Point MidPerim(double d)const;
+	Point MidParam(double param)const;
+	double Length()const;
+	Point GetVector(double fraction)const;
+	void Intersect(const Span& s, std::list<Point> &pts)const; // finds all the intersection points between two spans
+};
+
+class CCurve
+{
+	// a closed curve, please make sure you add an end point, the same as the start point
+
+protected:
+	void AddArcOrLines(bool check_for_arc, std::list<CVertex> &new_vertices, std::list<const CVertex*>& might_be_an_arc, CArc &arc, bool &arc_found, bool &arc_added);
+	bool CheckForArc(const CVertex& prev_vt, std::list<const CVertex*>& might_be_an_arc, CArc &arc);
+
+public:
+	std::list<CVertex> m_vertices;
+	void append(const CVertex& vertex);
+
+	void FitArcs();
+	void UnFitArcs();
+	Point NearestPoint(const Point& p)const;
+	Point NearestPoint(const CCurve& p, double *d = NULL)const;
+	Point NearestPoint(const Span& p, double *d = NULL)const;
+	void GetBox(CBox2D &box);
+	void Reverse();
+	double GetArea()const;
+	bool IsClockwise()const{return GetArea()>0;}
+	bool IsClosed()const;
+	void ChangeStart(const Point &p);
+	void ChangeEnd(const Point &p);
+	bool Offset(double leftwards_value);
+	void OffsetForward(double forwards_value, bool refit_arcs = true); // for drag-knife compensation
+	void Break(const Point &p);
+	void ExtractSeparateCurves(const std::list<Point> &ordered_points, std::list<CCurve> &separate_curves)const;
+	double Perim()const;
+	Point PerimToPoint(double perim)const;
+	double PointToPerim(const Point& p)const;
+	void GetSpans(std::list<Span> &spans)const;
+	void RemoveTinySpans();
+	void operator+=(const CCurve& p);
+	void SpanIntersections(const Span& s, std::list<Point> &pts)const;
+	void CurveIntersections(const CCurve& c, std::list<Point> &pts)const;
+};
+
+void tangential_arc(const Point &p0, const Point &p1, const Point &v0, Point &c, int &dir);

src/Mod/Path/libarea/Point.h

@@ -0,0 +1,51 @@
+// Point.h
+
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#pragma once
+
+#include <cmath>
+#include "kurve/geometry.h"
+
+class Point{
+public:
+	// can be a position, or a vector
+	double x, y;
+
+	Point():x(0.0), y(0.0){}
+	Point(double X, double Y):x(X), y(Y){}
+	Point(const double* p):x(p[0]), y(p[1]){}
+	Point(const Point& p0, const Point& p1):x(p1.x - p0.x), y(p1.y - p0.y){} // vector from p0 to p1
+
+	static double tolerance;
+
+	const Point operator+(const Point& p)const{return Point(x + p.x, y + p.y);}
+	const Point operator-(const Point& p)const{return Point(x - p.x, y - p.y);}
+	const Point operator*(double d)const{return Point(x * d, y * d);}
+	const Point operator/(double d)const{return Point(x / d, y / d);}
+	bool operator==(const Point& p)const{return fabs(x-p.x)<tolerance && fabs(y-p.y)<tolerance;}
+	bool operator!=(const Point &p)const{ return !(*this == p);}
+	double dist(const Point &p)const{double dx = p.x - x; double dy = p.y - y; return sqrt(dx*dx + dy*dy);}
+    double length()const;
+    double normalize();
+	double operator*(const Point &p)const{return (x * p.x + y * p.y);}// dot product
+	double operator^(const Point &p)const{return (x * p.y - y * p.x);}// cross product m0.m1.sin a = v0 ^ v1
+	Point operator~(void)const{return Point(-y, x);}// perp to left
+	Point operator-(void)const{return Point(-x, -y);}// v1 = -v0;  (unary minus)
+	void Rotate(double cosa, double sina){// rotate vector by angle
+		double temp = -y * sina + x * cosa;
+		y = x * sina + cosa * y;
+		x = temp;
+	}	
+	void Rotate(double angle){if(fabs(angle) < 1.0e-09)return; Rotate(cos(angle), sin(angle));}
+	void Transform(const geoff_geometry::Matrix &m)
+	{
+		geoff_geometry::Point p(x,y);
+		p = p.Transform(m);
+		x = p.x;
+		y = p.y;
+	}
+};
+
+const Point operator*(const double &d, const Point &p);

src/Mod/Path/libarea/PythonStuff.cpp

@@ -0,0 +1,409 @@
+// PythonStuff.cpp
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+#include "PythonStuff.h"
+
+#include "Area.h"
+#include "Point.h"
+#include "AreaDxf.h"
+#include "kurve/geometry.h"
+
+#if _DEBUG
+#undef _DEBUG
+#include <Python.h>
+#define _DEBUG
+#else
+#include <Python.h>
+#endif
+
+#ifdef __GNUG__
+#pragma implementation
+#endif
+
+#include <boost/progress.hpp>
+#include <boost/timer.hpp>
+#include <boost/foreach.hpp>
+#include <boost/python.hpp>
+#include <boost/python/module.hpp>
+#include <boost/python/class.hpp>
+#include <boost/python/wrapper.hpp>
+#include <boost/python/call.hpp>
+
+#include "clipper.hpp"
+using namespace ClipperLib;
+
+
+namespace bp = boost::python;
+
+boost::python::list getVertices(const CCurve& curve) {
+	boost::python::list vlist;
+	BOOST_FOREACH(const CVertex& vertex, curve.m_vertices) {
+		vlist.append(vertex);
+    }
+	return vlist;
+}
+
+boost::python::list getCurves(const CArea& area) {
+	boost::python::list clist;
+	BOOST_FOREACH(const CCurve& curve, area.m_curves) {
+		clist.append(curve);
+    }
+	return clist;
+}
+
+boost::python::tuple transformed_point(const geoff_geometry::Matrix &matrix, double x, double y, double z)
+{
+	geoff_geometry::Point3d p(x,y,z);
+	p = p.Transform(matrix);
+
+	return bp::make_tuple(p.x,p.y,p.z);
+}
+
+static void print_curve(const CCurve& c)
+{
+	unsigned int nvertices = c.m_vertices.size();
+	printf("number of vertices = %d\n", nvertices);
+	int i = 0;
+	for(std::list<CVertex>::const_iterator It = c.m_vertices.begin(); It != c.m_vertices.end(); It++, i++)
+	{
+		const CVertex& vertex = *It;
+		printf("vertex %d type = %d, x = %g, y = %g", i+1, vertex.m_type, vertex.m_p.x / CArea::m_units, vertex.m_p.y / CArea::m_units);
+		if(vertex.m_type)printf(", xc = %g, yc = %g", vertex.m_c.x / CArea::m_units, vertex.m_c.y / CArea::m_units);
+		printf("\n");
+	}
+}
+
+static void print_area(const CArea &a)
+{
+	for(std::list<CCurve>::const_iterator It = a.m_curves.begin(); It != a.m_curves.end(); It++)
+	{
+		const CCurve& curve = *It;
+		print_curve(curve);
+	}
+}
+
+static unsigned int num_vertices(const CCurve& curve)
+{
+	return curve.m_vertices.size();
+}
+
+static CVertex FirstVertex(const CCurve& curve)
+{
+	return curve.m_vertices.front();
+}
+
+static CVertex LastVertex(const CCurve& curve)
+{
+	return curve.m_vertices.back();
+}
+
+static void set_units(double units)
+{
+	CArea::m_units = units;
+}
+
+static double get_units()
+{
+	return CArea::m_units;
+}
+
+static bool holes_linked()
+{
+	return CArea::HolesLinked();
+}
+
+static CArea AreaFromDxf(const char* filepath)
+{
+	CArea area;
+	AreaDxfRead dxf(&area, filepath);
+	dxf.DoRead();
+	return area;
+}
+
+static void append_point(CCurve& c, const Point& p)
+{
+	c.m_vertices.push_back(CVertex(p));
+}
+
+static boost::python::tuple nearest_point_to_curve(CCurve& c1, const CCurve& c2)
+{
+	double dist;
+	Point p = c1.NearestPoint(c2, &dist);
+
+	return bp::make_tuple(p, dist);
+}
+
+boost::python::list MakePocketToolpath(const CArea& a, double tool_radius, double extra_offset, double stepover, bool from_center, bool use_zig_zag, double zig_angle)
+{
+	std::list<CCurve> toolpath;
+
+	CAreaPocketParams params(tool_radius, extra_offset, stepover, from_center, use_zig_zag ? ZigZagPocketMode : SpiralPocketMode, zig_angle);
+	a.SplitAndMakePocketToolpath(toolpath, params);
+
+	boost::python::list clist;
+	BOOST_FOREACH(const CCurve& c, toolpath) {
+		clist.append(c);
+    }
+	return clist;
+}
+
+boost::python::list SplitArea(const CArea& a)
+{
+	std::list<CArea> areas;
+	a.Split(areas);
+
+	boost::python::list alist;
+	BOOST_FOREACH(const CArea& a, areas) {
+		alist.append(a);
+    }
+	return alist;
+}
+
+void dxfArea(CArea& area, const char* str)
+{
+	area = CArea();
+}
+
+boost::python::list getCurveSpans(const CCurve& c)
+{
+	boost::python::list span_list;
+	const Point *prev_p = NULL;
+
+	for(std::list<CVertex>::const_iterator VIt = c.m_vertices.begin(); VIt != c.m_vertices.end(); VIt++)
+	{
+		const CVertex& vertex = *VIt;
+
+		if(prev_p)
+		{
+			span_list.append(Span(*prev_p, vertex));
+		}
+		prev_p = &(vertex.m_p);
+	}
+
+	return span_list;
+}
+
+Span getFirstCurveSpan(const CCurve& c)
+{
+	if(c.m_vertices.size() < 2)return Span();
+
+	std::list<CVertex>::const_iterator VIt = c.m_vertices.begin();
+	const Point &p = (*VIt).m_p;
+	VIt++;
+	return Span(p, *VIt, true);
+}
+
+Span getLastCurveSpan(const CCurve& c)
+{
+	if(c.m_vertices.size() < 2)return Span();
+
+	std::list<CVertex>::const_reverse_iterator VIt = c.m_vertices.rbegin();
+	const CVertex &v = (*VIt);
+	VIt++;
+
+	return Span((*VIt).m_p, v, c.m_vertices.size() == 2);
+}
+
+bp::tuple TangentialArc(const Point &p0, const Point &p1, const Point &v0)
+{
+  Point c;
+  int dir;
+  tangential_arc(p0, p1, v0, c, dir);
+
+  return bp::make_tuple(c, dir);
+}
+
+boost::python::list spanIntersect(const Span& span1, const Span& span2) {
+	boost::python::list plist;
+	std::list<Point> pts;
+	span1.Intersect(span2, pts);
+	BOOST_FOREACH(const Point& p, pts) {
+		plist.append(p);
+    }
+	return plist;
+}
+
+//Matrix(boost::python::list &l){}
+
+boost::shared_ptr<geoff_geometry::Matrix> matrix_constructor(const boost::python::list& lst) {
+	double m[16] = {1,0,0,0,0,1,0,0, 0,0,1,0, 0,0,0,1};
+
+  boost::python::ssize_t n = boost::python::len(lst);
+  int j = 0;
+  for(boost::python::ssize_t i=0;i<n;i++) {
+    boost::python::object elem = lst[i];
+	m[j] = boost::python::extract<double>(elem.attr("__float__")());
+	j++;
+	if(j>=16)break;
+  }
+
+  return boost::shared_ptr<geoff_geometry::Matrix>( new geoff_geometry::Matrix(m) );
+}
+
+boost::python::list InsideCurves(const CArea& a, const CCurve& curve) {
+	boost::python::list plist;
+
+	std::list<CCurve> curves_inside;
+	a.InsideCurves(curve, curves_inside);
+	BOOST_FOREACH(const CCurve& c, curves_inside) {
+		plist.append(c);
+    }
+	return plist;
+}
+
+boost::python::list CurveIntersections(const CCurve& c1, const CCurve& c2) {
+	boost::python::list plist;
+
+	std::list<Point> pts;
+	c1.CurveIntersections(c2, pts);
+	BOOST_FOREACH(const Point& p, pts) {
+		plist.append(p);
+    }
+	return plist;
+}
+
+boost::python::list AreaIntersections(const CArea& a, const CCurve& c2) {
+	boost::python::list plist;
+
+	std::list<Point> pts;
+	a.CurveIntersections(c2, pts);
+	BOOST_FOREACH(const Point& p, pts) {
+		plist.append(p);
+    }
+	return plist;
+}
+
+double AreaGetArea(const CArea& a)
+{
+	return a.GetArea();
+}
+
+BOOST_PYTHON_MODULE(area) {
+	bp::class_<Point>("Point") 
+        .def(bp::init<double, double>())
+        .def(bp::init<Point>())
+        .def(bp::other<double>() * bp::self)
+        .def(bp::self * bp::other<double>())
+        .def(bp::self / bp::other<double>())
+        .def(bp::self * bp::other<Point>())
+        .def(bp::self - bp::other<Point>())
+        .def(bp::self + bp::other<Point>())
+        .def(bp::self ^ bp::other<Point>())
+        .def(bp::self == bp::other<Point>())
+        .def(bp::self != bp::other<Point>())
+        .def(-bp::self)
+        .def(~bp::self)
+        .def("dist", &Point::dist)
+        .def("length", &Point::length)
+        .def("normalize", &Point::normalize)
+		.def("Rotate", static_cast< void (Point::*)(double, double) >(&Point::Rotate))
+		.def("Rotate", static_cast< void (Point::*)(double) >(&Point::Rotate))
+        .def_readwrite("x", &Point::x)
+        .def_readwrite("y", &Point::y)
+		.def("Transform", &Point::Transform)
+    ;
+
+	bp::class_<CVertex>("Vertex") 
+        .def(bp::init<CVertex>())
+        .def(bp::init<int, Point, Point>())
+        .def(bp::init<Point>())
+        .def(bp::init<int, Point, Point, int>())
+        .def_readwrite("type", &CVertex::m_type)
+        .def_readwrite("p", &CVertex::m_p)
+        .def_readwrite("c", &CVertex::m_c)
+        .def_readwrite("user_data", &CVertex::m_user_data)
+    ;
+
+	bp::class_<Span>("Span") 
+        .def(bp::init<Span>())
+        .def(bp::init<Point, CVertex, bool>())
+		.def("NearestPoint", static_cast< Point (Span::*)(const Point& p)const >(&Span::NearestPoint))
+		.def("NearestPoint", static_cast< Point (Span::*)(const Span& p, double *d)const >(&Span::NearestPoint))
+		.def("GetBox", &Span::GetBox)
+		.def("IncludedAngle", &Span::IncludedAngle)
+		.def("GetArea", &Span::GetArea)
+		.def("On", &Span::On)
+		.def("MidPerim", &Span::MidPerim)
+		.def("MidParam", &Span::MidParam)
+		.def("Length", &Span::Length)
+		.def("GetVector", &Span::GetVector)
+		.def("Intersect", &spanIntersect)
+        .def_readwrite("p", &Span::m_p)
+		.def_readwrite("v", &Span::m_v)
+    ;
+
+	bp::class_<CCurve>("Curve") 
+        .def(bp::init<CCurve>())
+        .def("getVertices", &getVertices)
+        .def("append",&CCurve::append)
+        .def("append",&append_point)
+        .def("text", &print_curve)
+		.def("NearestPoint", static_cast< Point (CCurve::*)(const Point& p)const >(&CCurve::NearestPoint))
+		.def("NearestPoint", &nearest_point_to_curve)
+		.def("Reverse", &CCurve::Reverse)
+		.def("getNumVertices", &num_vertices)
+		.def("FirstVertex", &FirstVertex)
+		.def("LastVertex", &LastVertex)
+		.def("GetArea", &CCurve::GetArea)
+		.def("IsClockwise", &CCurve::IsClockwise)
+		.def("IsClosed", &CCurve::IsClosed)
+        .def("ChangeStart",&CCurve::ChangeStart)
+        .def("ChangeEnd",&CCurve::ChangeEnd)
+        .def("Offset",&CCurve::Offset)
+        .def("OffsetForward",&CCurve::OffsetForward)
+        .def("GetSpans",&getCurveSpans)
+        .def("GetFirstSpan",&getFirstCurveSpan)
+        .def("GetLastSpan",&getLastCurveSpan)
+        .def("Break",&CCurve::Break)
+        .def("Perim",&CCurve::Perim)
+        .def("PerimToPoint",&CCurve::PerimToPoint)
+        .def("PointToPerim",&CCurve::PointToPerim)
+		.def("FitArcs",&CCurve::FitArcs)
+        .def("UnFitArcs",&CCurve::UnFitArcs)
+        .def("Intersections",&CurveIntersections)
+    ;
+
+	bp::class_<CBox2D>("Box") 
+        .def(bp::init<CBox2D>())
+		.def("MinX", &CBox2D::MinX)
+		.def("MaxX", &CBox2D::MaxX)
+		.def("MinY", &CBox2D::MinY)
+		.def("MaxY", &CBox2D::MaxY)
+    ;
+
+	bp::class_<CArea>("Area") 
+        .def(bp::init<CArea>())
+        .def("getCurves", &getCurves)
+        .def("append",&CArea::append)
+        .def("Subtract",&CArea::Subtract)
+        .def("Intersect",&CArea::Intersect)
+        .def("Union",&CArea::Union)
+        .def("Offset",&CArea::Offset)
+        .def("FitArcs",&CArea::FitArcs)
+        .def("text", &print_area)
+		.def("num_curves", &CArea::num_curves)
+		.def("NearestPoint", &CArea::NearestPoint)
+		.def("GetBox", &CArea::GetBox)
+		.def("Reorder", &CArea::Reorder)
+		.def("MakePocketToolpath", &MakePocketToolpath)
+		.def("Split", &SplitArea)
+		.def("InsideCurves", &InsideCurves)
+		.def("Thicken", &CArea::Thicken)
+        .def("Intersections",&AreaIntersections)
+        .def("GetArea",&AreaGetArea)
+    ;
+
+	bp::class_<geoff_geometry::Matrix, boost::shared_ptr<geoff_geometry::Matrix> > ("Matrix")
+        .def(bp::init<geoff_geometry::Matrix>())
+	    .def("__init__", bp::make_constructor(&matrix_constructor))
+	    .def("TransformedPoint", &transformed_point)
+		.def("Multiply", &geoff_geometry::Matrix::Multiply)
+	;
+
+    bp::def("set_units", set_units);
+    bp::def("get_units", get_units);
+    bp::def("holes_linked", holes_linked);
+    bp::def("AreaFromDxf", AreaFromDxf);
+    bp::def("TangentialArc", TangentialArc);
+}

src/Mod/Path/libarea/PythonStuff.h

@@ -0,0 +1,10 @@
+// PythonStuff.h
+
+// Copyright 2011, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+
+extern void Message(const char*);
+
+void PythonInit();
+void PythonFinish();
+

src/Mod/Path/libarea/dxf.h

@@ -0,0 +1,156 @@
+// dxf.h
+// Copyright (c) 2009, Dan Heeks
+// This program is released under the BSD license. See the file COPYING for details.
+#pragma once
+
+#include <algorithm>
+#include <list>
+#include <vector>
+#include <map>
+#include <set>
+#include <fstream>
+#include <sstream>
+#include <iostream>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+//Following is required to be defined on Ubuntu with OCC 6.3.1
+#ifndef HAVE_IOSTREAM
+#define HAVE_IOSTREAM
+#endif
+
+typedef int Aci_t; // AutoCAD color index
+
+typedef enum
+{
+	eUnspecified = 0,	// Unspecified (No units)
+	eInches,
+	eFeet,
+	eMiles,
+	eMillimeters,
+	eCentimeters,
+	eMeters,
+	eKilometers,
+	eMicroinches,
+	eMils,
+	eYards,
+	eAngstroms,
+	eNanometers,
+	eMicrons,
+	eDecimeters,
+	eDekameters,
+	eHectometers,
+	eGigameters,
+	eAstronomicalUnits,
+	eLightYears,
+	eParsecs
+} eDxfUnits_t;
+
+
+struct SplineData
+{
+	double norm[3];
+	int degree;
+	int knots;
+	int control_points;
+	int fit_points;
+	int flag;
+	std::list<double> starttanx;
+	std::list<double> starttany;
+	std::list<double> starttanz;
+	std::list<double> endtanx;
+	std::list<double> endtany;
+	std::list<double> endtanz;
+	std::list<double> knot;
+	std::list<double> weight;
+	std::list<double> controlx;
+	std::list<double> controly;
+	std::list<double> controlz;
+	std::list<double> fitx;
+	std::list<double> fity;
+	std::list<double> fitz;
+};
+
+class CDxfWrite{
+private:
+	std::ofstream* m_ofs;
+	bool m_fail;
+
+public:
+	CDxfWrite(const char* filepath);
+	~CDxfWrite();
+
+	bool Failed(){return m_fail;}
+
+	void WriteLine(const double* s, const double* e, const char* layer_name );
+	void WritePoint(const double*, const char*);
+	void WriteArc(const double* s, const double* e, const double* c, bool dir, const char* layer_name );
+    void WriteEllipse(const double* c, double major_radius, double minor_radius, double rotation, double start_angle, double end_angle, bool dir, const char* layer_name );
+	void WriteCircle(const double* c, double radius, const char* layer_name );
+};
+
+// derive a class from this and implement it's virtual functions
+class CDxfRead{
+private:
+	std::ifstream* m_ifs;
+
+	bool m_fail;
+	char m_str[1024];
+	char m_unused_line[1024];
+	eDxfUnits_t m_eUnits;
+	char m_layer_name[1024];
+	char m_section_name[1024];
+	char m_block_name[1024];
+	bool m_ignore_errors;
+
+
+	typedef std::map< std::string,Aci_t > LayerAciMap_t;
+	LayerAciMap_t m_layer_aci;  // layer names -> layer color aci map
+
+	bool ReadUnits();
+	bool ReadLayer();
+	bool ReadLine();
+	bool ReadText();
+	bool ReadArc();
+	bool ReadCircle();
+	bool ReadEllipse();
+	bool ReadPoint();
+	bool ReadSpline();
+	bool ReadLwPolyLine();
+	bool ReadPolyLine();
+	bool ReadVertex(double *pVertex, bool *bulge_found, double *bulge);
+	void OnReadArc(double start_angle, double end_angle, double radius, const double* c);
+	void OnReadCircle(const double* c, double radius);
+    void OnReadEllipse(const double* c, const double* m, double ratio, double start_angle, double end_angle);
+
+	void get_line();
+	void put_line(const char *value);
+	void DerefACI();
+
+protected:
+	Aci_t m_aci; // manifest color name or 256 for layer color
+
+public:
+	CDxfRead(const char* filepath); // this opens the file
+	~CDxfRead(); // this closes the file
+
+	bool Failed(){return m_fail;}
+	void DoRead(const bool ignore_errors = false); // this reads the file and calls the following functions
+
+	double mm( const double & value ) const;
+
+	bool IgnoreErrors() const { return(m_ignore_errors); }
+
+	virtual void OnReadLine(const double* s, const double* e){}
+	virtual void OnReadPoint(const double* s){}
+	virtual void OnReadText(const double* point, const double height, const char* text){}
+	virtual void OnReadArc(const double* s, const double* e, const double* c, bool dir){}
+	virtual void OnReadCircle(const double* s, const double* c, bool dir){}
+	virtual void OnReadEllipse(const double* c, double major_radius, double minor_radius, double rotation, double start_angle, double end_angle, bool dir){}
+	virtual void OnReadSpline(struct SplineData& sd){}
+	virtual void AddGraphics() const { }
+
+    std::string LayerName() const;
+
+};

src/Mod/Path/libarea/kurve/Construction.cpp

@@ -0,0 +1,844 @@
+// ***************************************************************************************************************************************
+//                    Point, CLine & Circle classes part of geometry.lib
+//                    g.j.hawkesford August 2006 for Camtek Gmbh
+//
+// This program is released under the BSD license. See the file COPYING for details.
+//
+// ***************************************************************************************************************************************
+
+#include "geometry.h"
+using namespace geoff_geometry;
+
+namespace geoff_geometry {
+	int   UNITS = MM;
+	double TOLERANCE = 1.0e-06;
+	double TOLERANCE_SQ = TOLERANCE * TOLERANCE;
+	double TIGHT_TOLERANCE = 1.0e-09;
+	double UNIT_VECTOR_TOLERANCE = 1.0e-10;
+	double RESOLUTION = 1.0e-06;
+
+	// dummy functions
+	const wchar_t* getMessage(const wchar_t* original, int messageGroup, int stringID){return original;}
+	const wchar_t* getMessage(const wchar_t* original){return original;}
+	void FAILURE(const wchar_t* str){throw(str);}
+	void FAILURE(const std::wstring& str){throw(str);}
+
+	void set_Tolerances(int mode) {
+		UNIT_VECTOR_TOLERANCE = 1.0e-10;
+		switch (UNITS = mode)
+		{
+		case MM:
+			geoff_geometry::TOLERANCE = 1.0e-03;					// Peps
+			RESOLUTION = 1.0e-03;
+			TIGHT_TOLERANCE = 1.0e-06;
+			break;
+		case INCHES:
+			TOLERANCE = 1.0e-04;					// Peps
+			RESOLUTION = 1.0e-04;
+			TIGHT_TOLERANCE = 1.0e-7;
+			break;
+		case METRES:
+			TOLERANCE = 1.0e-06;					// p4c...SW
+			RESOLUTION = 1.0e-06;
+			TIGHT_TOLERANCE = 1.0e-09;
+			break;
+		default:
+			FAILURE(L"INVALID UNITS");
+		}
+		TOLERANCE_SQ = TOLERANCE * TOLERANCE;
+	}
+
+	double mm(double value) {
+		switch(UNITS) {
+			default:
+				return value;
+			case METRES:
+				return value * .001;
+			case INCHES:
+				return value / 25.4;
+		}
+	}
+
+	// ostream operators  = non-member overload
+	// *********************************************************************************************************
+	wostream& operator << (wostream& op, Point& p){
+		// for debug - print point to file
+		if(p.ok == false)
+			op << L" ok=\"false\"";
+		else
+			op << L" x=\"" << p.x << L"\" y=\"" << p.y << L"\"";
+		return op;
+	}
+
+	wostream& operator <<(wostream& op, CLine& cl){
+		// for debug - print cline to file
+		if(cl.ok == false)
+			op << L"(CLine UNSET)";
+		else
+			op << L"sp=" << cl.p << L" v=" << cl.v;
+		return op;
+	}
+
+	wostream& operator <<(wostream& op, Plane& pl){
+		// for debug - print plane to file stream
+		if(pl.ok == false)
+			op << L"(Plane UNSET)";
+		else
+			op << L"d=" << pl.d << L" normal=" << pl.normal;
+		return op;
+	}
+
+	ostream& operator << (ostream& op, Point3d& p){
+		// for debug - print point to file
+//		if(p.ok == false)
+//			op << "ok=\"false\"";
+//		else
+			op << "x=\"" << p.x << "\" y=\"" << p.y << "\" z=" << p.z << "\"";
+		return op;
+
+	}
+
+	wostream& operator <<(wostream& op, Vector2d& v){
+		// for debug - print vector to file
+		op << L"(" << v.getx() << L", " << v.gety() << L")";
+		return op;
+	}
+
+	wostream& operator <<(wostream& op, Vector3d& v){
+		// for debug - print vector to file
+		op << L"(" << v.getx() << L", " << v.gety() << L"," << v.getz() << L")";
+		return op;
+	}
+
+	wostream& operator <<(wostream& op, Circle& c){
+		// for debug - print circle to file
+		if(c.ok == false)
+			op << L"ok=\"false\"";
+		else
+			op << L" x=\"" << c.pc.x << L"\" y=\"" << c.pc.y << L"\" radius=\"" << c.radius << L"\"";
+		return op;
+	}
+
+	wostream& operator <<(wostream& op, Span& sp){
+		// for debug - print span to file stream		
+		op << L"p0 = " << sp.p0 << L" p1=" << sp.p1;
+		if(sp.dir) {
+			op << L" pc=" << sp.pc << L" dir=" << ((sp.dir == CW)?L"CW" : L"ACW") << L" radius=" << sp.radius;
+		}
+		return op;
+	}
+
+
+
+	// ***************************************************************************************************************************************
+	// point classes
+	// ***************************************************************************************************************************************
+	Point::Point( const Point3d& p ) {												// copy constructor  Point p1(p2);
+			x = p.x;
+			y = p.y;
+//			ok = p.ok;
+			ok = true;
+	}
+
+	Point::Point(const Vector2d& v)
+	{
+		x = v.getx(); y = v.gety();
+	}
+
+	Point3d::Point3d(const Vector3d& v) {
+			x = v.getx(); y = v.gety();  z = v.getz();// ok = true;
+	}
+
+	bool Point3d::operator==(const Point3d &p)const{
+		// p1 == p2 (uses TOLERANCE)
+		if(FNE(this->x, p.x, TOLERANCE) || FNE(this->y, p.y, TOLERANCE) || FNE(this->z, p.z, TOLERANCE)) return false;
+		return true;
+	}
+
+	Point Point::Transform(const Matrix& m) {
+		// transform Point
+		Point ret;
+		m.Transform2d(&x, &ret.x);
+		ret.ok = true;
+		return ret;
+	}
+	Point3d Point3d::Transform(const Matrix& m) {
+		// transform Point
+		Point3d ret;
+		m.Transform(&x, &ret.x);
+//		ret.ok = true;
+		return ret;
+	}
+
+	Point Point::operator+(const Vector2d &v)const{
+		return Point(x + v.getx(), y + v.gety());
+	}
+
+	Point3d Point3d::operator+(const Vector3d &v)const{
+		return Point3d(x + v.getx(), y + v.gety(), z + v.getz());
+	}
+
+	bool Point::operator==(const Point &p) const{
+		// p1 == p2 (uses TOLERANCE)
+		if(FNE(this->x, p.x, TOLERANCE) || FNE(this->y, p.y, TOLERANCE)) return false;
+		return true;
+	}
+
+
+
+	double	Point::Dist(const Point& p)const{													// distance between 2 points
+		return Vector2d(*this, p).magnitude();
+	}
+
+	double	Point::DistSq(const Point& p)const{													// distance squared between 2 points
+		return Vector2d(*this, p).magnitudesqd();
+	}
+
+	double Point3d::Dist(const Point3d& p)const {												// distance between 2 points
+		return Vector3d(*this, p).magnitude();
+	}
+
+	double Point3d::DistSq(const Point3d& p)const {			// distance squared
+		return (this->x - p.x) * (this->x - p.x) + (this->y - p.y) * (this->y - p.y) + (this->z - p.z) * (this->z - p.z);
+	}
+
+	Point Point::Mid(const Point& p1, double factor)const{
+		// Mid
+		return geoff_geometry::Mid(*this, p1, factor);
+	}
+
+	Point3d Point3d::Mid(const Point3d& p, double factor)const{
+		// Mid
+		return Vector3d(*this, p) * factor + *this;
+	}
+
+	Point Mid(const Point& p0, const Point& p1, double factor){
+		// mid or partway between 2 points
+		return Vector2d(p0, p1) * factor + p0;
+	}
+	Point Rel(const Point& p, double x0, double y0) {
+		// Relative point
+		return (p.ok)?Point(p.x + x0, p.y + y0) : INVALID_POINT;
+	}
+
+	Point Polar(const Point& p, double angle, double r) {
+		// polar from this point
+		angle *= DegreesToRadians;
+		return (p.ok)?Point(p.x + r * cos(angle), p.y + r * sin(angle)) : INVALID_POINT;
+	}
+
+	// ***************************************************************************************************************************************
+	// clines
+	// ***************************************************************************************************************************************
+	//const CLine horiz(Point(0, 0), 1, 0);											// define global horizontal line
+
+	double CLine::c() {
+		// returns c for ax + by + c = 0 format (peps format where needed)
+		return (v.getx() * p.y - v.gety() * p.x);
+	}
+	void CLine::Normalise() {
+		// normalise the cline vector
+		ok = v.normalise() >= TOLERANCE;
+	}
+
+	CLine::CLine(const Span& sp){
+		p = sp.p0;
+		v = sp.vs; 
+		ok = sp.returnSpanProperties && !sp.NullSpan;
+	}
+
+	CLine Normal(const CLine& s) {
+		// returns normal to this line
+		return CLine(s.p, ~s.v, false);
+	}
+	const CLine CLine::operator ~(void){
+		return CLine(this->p, ~v, false);
+	}
+	CLine Normal(const CLine& s, const Point& p) {
+		// returns normal to this line thro' p
+		return CLine(p, ~s.v, false);
+	}
+
+	CLine CLine::Transform(Matrix& m) {
+
+		Point p0 = this->p;
+		Point p1(p0.x + v.getx(), p0.y + v.gety());
+		return CLine(p0.Transform(m), p1.Transform(m));
+	}
+
+
+	double CLine::Dist(const Point& p0)const {
+		// distance between cline & point  >0 cw about point   <0 acw about point
+		return this->v ^ Vector2d(p0, this->p);
+	}
+
+	double Point::Dist(const CLine& cl)const {
+		// distance between cline & point  >0 cw about point   <0 acw about point
+		return cl.v ^ Vector2d(*this, cl.p);
+	}
+
+	Point CLine::Intof(const CLine& s)	{
+		// Intof 2 Clines
+		return geoff_geometry::Intof(*this, s);
+	}
+
+	Point CLine::Intof(int NF, const Circle& c)	{
+		// Intof Cline & Circleconst 
+		return geoff_geometry::Intof(NF, *this, c);
+	}
+	Point CLine::Intof(int NF, const Circle& c, Point& otherInters)	{
+		// Intof Cline & Circle & other intersection
+		return geoff_geometry::Intof(NF, *this, c, otherInters);
+	}
+
+	Point Intof(const CLine& s0, const CLine& s1)	{
+		// inters of 2 clines  (parameterise lines x = x0 + t * dx)
+		double cp = s1.v ^ s0.v;
+		if(fabs (cp) > 1.0e-6) {
+			double t = (s1.v ^ Vector2d(s0.p, s1.p)) / cp;
+			return s0.v * t + s0.p;
+		}
+		return INVALID_POINT;
+	}
+	Point XonCLine(CLine& s, double xval) {
+		// return point given X on a line
+		return Intof(s, CLine(Point(xval,0),0,1,false));
+	}
+	Point YonCLine(CLine& s, double yval) {
+		// return point given Y on a line
+		return Intof(s, CLine(Point(0,yval),1,0,false));
+	}
+	Point Along(const CLine& s, double d) {
+		// distance along line
+		return Point(s.p.x + d * s.v.getx(), s.p.y + d * s.v.gety(), s.ok);
+	}
+
+	Point Along(const CLine& s, double d, Point& p) {
+		// distance along line from point
+		return Point(p.x + d * s.v.getx(), p.y + d * s.v.gety(), p.ok);
+	}
+	Point Around(const Circle& c, double d, const Point& p) {
+		// distance around circle from point
+		CLine radial(c.pc, p);
+		if(radial.ok) {
+			if(fabs(c.radius) > TOLERANCE ) {
+				double a = sin(- d / c.radius);
+				double b = cos(- d / c.radius);
+				return Point(c.pc.x - c.radius * (radial.v.gety() * a - radial.v.getx() * b), c.pc.y + c.radius * (radial.v.gety() * b + radial.v.getx() * a));
+			}
+		}
+		return INVALID_POINT;
+	}
+	CLine AtAngle(double angle, const Point& p0, const CLine& s) {
+		// cline at angle [to a cline] thro' a point
+		angle *= DegreesToRadians;
+		Vector2d v(cos(angle), sin(angle));
+		return CLine(p0, v.getx() * s.v.getx() - v.gety() * s.v.gety(), v.gety() * s.v.getx() + v.getx() * s.v.gety());
+	}
+	CLine Parallel(int side, const CLine& s0, double distance) {
+		// parallel to line by distance
+		Vector2d v = ~s0.v;
+		return CLine(v * ((double)side * distance) + s0.p, s0.v.getx(), s0.v.gety());
+	}
+
+	CLine Parallel(const CLine& s0, Point& p) {
+		// parallel to line through point
+		return CLine(p, s0.v.getx(), s0.v.gety());
+	}
+
+	CLine CLine::Bisector(const CLine& s) {
+		//  bisector of 2 clines
+		return CLine (this->Intof(s), this->v.getx() + s.v.getx(), this->v.gety() + s.v.gety());
+	}
+
+
+
+
+	// ***************************************************************************************************************************************
+	// circle methods
+	// ***************************************************************************************************************************************
+
+	Circle::Circle(const Point& p, double rad, bool okay){
+		// Circle
+		pc = p;
+		radius = rad;
+		ok = pc.ok;
+	}
+
+	Circle::Circle( const Point& p, const Point& pc0){
+		if((ok = (p.ok && pc0.ok))) {
+			pc = pc0;
+			radius = p.Dist(pc0);
+		}
+	}
+
+	Circle::Circle( const Span& sp){
+		pc = sp.pc;
+		radius = sp.radius;
+		ok = sp.returnSpanProperties;
+	}
+
+	bool Circle::operator==(const Circle &c)const{
+		// c1 == c2 (uses TOLERANCE)
+		return FEQ(this->radius, c.radius, TOLERANCE) && (this->pc == c.pc);
+	}
+
+	Circle Circle::Transform(Matrix& m) { // transform
+		Point p0 = this->pc;
+		double scale;
+		if(m.GetScale(scale) == false) FAILURE(getMessage(L"Differential Scale not allowed for this method", GEOMETRY_ERROR_MESSAGES, MES_DIFFSCALE));
+		return Circle(p0.Transform(m), radius * scale);
+	}
+
+	Point	Circle::Intof(int LR, const Circle& c1) {
+		// intof 2 circles
+		return geoff_geometry::Intof(LR, *this, c1);
+	}
+	Point	Circle::Intof(int LR, const Circle& c1, Point& otherInters) {
+		// intof 2 circles, (returns the other intersection)
+		return geoff_geometry::Intof(LR, *this, c1, otherInters);
+	}
+	int	Circle::Intof(const Circle& c1, Point& leftInters, Point& rightInters) {
+		// intof 2 circles, (returns the other intersection)
+		return geoff_geometry::Intof(*this, c1, leftInters, rightInters);
+	}
+
+	CLine	Circle::Tanto(int AT,  double angle, const CLine& s0) const{
+		// cline tanto circle at angle to optional cline
+		return geoff_geometry::Tanto(AT, *this, angle, s0);
+	}
+
+	CLine Tanto(int AT, const Circle& c, const Point& p) {
+		// CLine tangent to a circle through a point
+		Vector2d v(p, c.pc);
+		double d = v.magnitude();
+		CLine s(p, ~v, false);								// initialise cline
+
+		if ( d <  TOLERANCE || d <  fabs(c.radius) - TOLERANCE)	// point inside circle ?
+			return INVALID_CLINE;
+		else {
+			if(d > fabs(c.radius) + TOLERANCE) {				// point outside circle
+				v.Rotate(sqrt((d - c.radius) * (d + c.radius)), - AT * c.radius);
+				s.v = v;
+			}
+		}
+		s.Normalise();
+		return s;
+	}
+
+	CLine Tanto(int AT0, const Circle& c0, int AT1, const Circle& c) {
+		// cline tanto 2 circles
+		CLine s;
+		Circle c1 = c;
+		c1.radius -= (double) (AT0 * AT1) * c0.radius;
+		s = Tanto(AT1, c1, c0.pc);
+		s.p.x += (double) AT0 * c0.radius * s.v.gety();
+		s.p.y -= (double) AT0 * c0.radius * s.v.getx();
+		return s;
+	}
+
+	CLine Tanto(int AT, const Circle& c, double angle, const CLine& s0) {
+		// cline at an angle [to a cline] tanto a circle 
+		CLine s = AtAngle(angle, c.pc, s0);
+		s.p.x += (double) AT * c.radius * s.v.gety();
+		s.p.y -= (double) AT * c.radius * s.v.getx();
+		//	s.p += ~s.v * (AT * c.radius);
+		s.ok = true;
+		return s;
+	}
+	Point AtAngle(const Circle& c, double angle) {
+		// Point at an angle on circle
+		angle *= DegreesToRadians;
+		return Point(c.pc.x + c.radius * cos(angle), c.pc.y + c.radius * sin(angle));
+	}
+
+	Point On(const CLine& s, const Point& p) {
+		// returns point that is nearest to s from p
+		double t = s.v * Vector2d(s.p, p);
+		return s.v * t + s.p;
+	}
+
+	Point On(const Circle& c, const Point& p) {
+		// returns point that is nearest to c from p
+		double r = p.Dist(c.pc);
+		if(r < TOLERANCE) FAILURE(getMessage(L",Point on Circle centre - On(Circle& c, Point& p)", GEOMETRY_ERROR_MESSAGES, MES_POINTONCENTRE));
+		return(Mid(p, c.pc, (r - c.radius) / r));
+	}
+
+
+	Point Intof( int NF, const CLine& s, const Circle& c) {
+		// inters of cline & circle  eg.     p1 = Intof(NEARINT, s1, c1);
+		Point otherInters;
+		return Intof(NF, s, c, otherInters);
+	}
+
+	Point Intof( int NF, const CLine& s, const Circle& c, Point& otherInters) {
+		// inters of cline & circle  eg.     p1 = Intof(NEARINT, s1, c1);
+		// otherInters returns the other intersection
+#if 1
+		// solving	x = x0 + dx * t			x = y0 + dy * t
+		//			x = xc + R * cos(a)		y = yc + R * sin(a)		for t
+		// gives :-  t<> (dx<64> + dy<64>) + 2t(dx*dx0 + dy*dy0) + (x0-xc)<29> + (y0-yc)<29> - R<> = 0
+		int nRoots;
+		double t, tFar, tNear, tOther;
+		Vector2d v0(c.pc, s.p);
+		if((nRoots = quadratic(1, 2 * (v0 * s.v), v0.magnitudesqd() - c.radius * c.radius, tFar, tNear)) != 0) {
+			if(nRoots == 2 && NF == NEARINT) {
+				t = tNear;
+				tOther = tFar;
+			} else {
+				t = tFar;
+				tOther = (nRoots == 2)?tNear : tFar;
+			}
+			otherInters =  s.v * tOther + s.p;
+			return s.v * t + s.p;
+		}
+		return INVALID_POINT;
+	}
+#else
+		// geometric solution - this is similar to the peps method, and it may offer better tolerancing than above??
+		Point intof;
+		CLine normal = Normal(s, c.pc);
+		intof = s.Intof(normal);
+		double d = intof.Dist(c.pc);
+
+		if(fabs(d - c.radius) < TOLERANCE) return intof;						// tangent (near enough for non-large radius I suppose?)
+
+		if(d > c.radius + TOLERANCE) return INVALID_POINT;					// no intersection
+
+		double q = (c.radius - d) * (c.radius + d);
+		if(q < 0) return intof;												// line inside tolerance
+
+		return Along(s, -(double)NF * sqrt(q), intof);						// 2 intersections (return near/far case)
+	}
+#endif
+	Point Intof( int intMode, const Circle& c0, const Circle& c1)	{
+		// inters of 2 circles		 eg.     p1 = Intof(LEFTINT, c1, c2)
+		Point otherInters;
+		return Intof(intMode, c0, c1, otherInters);
+	}
+
+	Point Intof( int intMode, const Circle& c0, const Circle& c1, Point& otherInters)	{
+		// inters of 2 circles		 eg.     p1 = Intof(LEFTINT, c1, c2);u
+		Point pLeft, pRight;
+		switch(Intof(c0, c1, pLeft, pRight)) {
+	default:
+		return INVALID_POINT;
+	case 1:
+		otherInters = pLeft;
+		return pLeft;
+	case 2:
+		if(intMode == LEFTINT) {
+			otherInters = pRight;
+			return pLeft;
+		}else {
+			otherInters = pLeft;
+			return pRight;
+		}
+		}
+	}
+
+	int Intof(const Circle& c0, const Circle& c1, Point& pLeft, Point& pRight)	{
+		// inters of 2 circles
+		// returns the number of intersctions
+		Vector2d v(c0.pc, c1.pc);
+		double d = 	v.normalise();
+		if(d < TOLERANCE)	return 0;									// co-incident circles
+
+		double sum = fabs(c0.radius) + fabs(c1.radius);
+		double diff = fabs(fabs(c0.radius) - fabs(c1.radius));
+		if(d > sum + TOLERANCE || d < diff - TOLERANCE) return 0;
+
+		// dist from centre of this circle to mid intersection
+		double d0 = 0.5 * (d + (c0.radius + c1.radius) * (c0.radius - c1.radius) / d);
+		if(d0 - c0.radius > TOLERANCE) return 0;						// circles don't intersect
+
+		double h = (c0.radius - d0) * (c0.radius + d0);				// half distance between intersects squared
+		if(h < 0) d0 = c0.radius;									// tangent
+		pLeft = v * d0 + c0.pc;										// mid-point of intersects
+		if(h < TOLERANCE_SQ) return 1;						// tangent
+		h = sqrt(h);
+
+		v = ~v;														// calculate 2 intersects
+		pRight = v * h + pLeft;
+		v = -v;
+		pLeft = v * h + pLeft; 
+		return 2;
+	}
+
+	Circle	Tanto(int NF, CLine& s0, Point& p, double rad) {
+		// circle tanto a CLine thro' a point
+		double d = s0.Dist(p);
+		if(fabs(d) > rad + TOLERANCE) return INVALID_CIRCLE;	// point too far from line
+		CLine s0offset = Parallel(RIGHTINT, s0, rad);
+
+		return Circle(Intof(NF, s0offset, Circle(p, rad)), rad);
+	}
+
+	Circle	Tanto(int AT1, CLine& s1, int AT2, CLine& s2, double rad) {
+		// circle tanto 2 clines with radius
+		CLine Offs1	= Parallel(AT1, s1, rad);
+		CLine Offs2	= Parallel(AT2, s2, rad);
+		Point pc = Intof(Offs1, Offs2);
+		return (pc.ok)? Circle(pc, rad) : INVALID_CIRCLE;
+	}
+	Circle Tanto(int AT1, CLine s1, int AT2, CLine s2, int AT3, CLine s3) {
+		// circle tanto 3 CLines
+		double s1c = s1.c(), s2c = s2.c(), s3c = s3.c();
+		double d =	  s1.v.gety() * (AT2 * s3.v.getx() - AT3 * s2.v.getx())
+			+ s2.v.gety() * (AT3 * s1.v.getx() - AT1 * s3.v.getx())
+			+ s3.v.gety() * (AT1 * s2.v.getx() - AT2 * s1.v.getx());
+		if(fabs(d) < UNIT_VECTOR_TOLERANCE) return INVALID_CIRCLE;
+		double radius =  (s1.v.gety() * (s2.v.getx() * s3c - s3.v.getx() * s2c)
+			+ s2.v.gety() * (s3.v.getx() * s1c - s1.v.getx() * s3c)
+			+ s3.v.gety() * (s1.v.getx() * s2c - s2.v.getx() * s1c)) / d ;
+		if(radius < TOLERANCE) return INVALID_CIRCLE;
+
+		CLine Offs1	= Parallel(AT1, s1, radius);
+		CLine Offs2	= Parallel(AT2, s2, radius);
+
+		Point p = Intof(Offs1, Offs2);
+		if(!p.ok) {
+			CLine Offs3	= Parallel(AT3, s3, radius);				// s1 & s2 parallel
+			p = Intof(Offs1, Offs3);
+			if(!p.ok) return INVALID_CIRCLE;						// 3 parallel lines
+		}
+		return Circle(p, radius);
+	}
+	Circle	Thro(int LR, const Point& p0, const Point& p1, double rad) {
+		// circle thro' 2 points, given radius and side
+		CLine thro(p0, p1);
+		if(thro.ok) {
+			double d = 0.5 * p0.Dist(p1);
+			Point pm = Mid(p0, p1);
+
+			if(d > rad + TOLERANCE) return INVALID_CIRCLE;
+			else if(d > rad - TOLERANCE) {
+				// within tolerance of centre of 2 points
+				return Circle(pm, d);
+			}
+			else {
+				// 2 solutions
+				return Circle(Along(Normal(thro, pm), (double)LR * sqrt((rad + d) * (rad - d)), pm), rad);
+			}
+		}
+		return INVALID_CIRCLE;
+	}
+
+	Circle	Thro(const Point& p0, const Point& p1) {
+		// circle thro 2 points (diametric)
+		return Circle(p0.Mid(p1), .5*p0.Dist(p1));
+	}
+	Circle	Thro(const Point& p0, const Point& p1, const Point& p2) {
+		// circle thro 3 points
+		CLine s0(p0, p1);
+		if(!s0.ok) return Thro(p1,p2);		// p0 & p1 coincident
+
+		CLine s1(p0, p2);
+		if(!s1.ok) return Thro(p0, p1);		// p0 & p2 coincident
+
+		CLine s2(p2, p1);
+		if(!s2.ok) return Thro(p0, p2);		// p1 & p2 coincident
+
+		Point p = Intof(Normal(s0, Mid(p0, p1)),  Normal(s1, Mid(p0, p2)));
+		return (p.ok)? Circle(p, p0.Dist(p), true) : INVALID_CIRCLE;
+	}
+	Circle	Tanto(int NF, int AT0, const CLine& s0, int AT1, const Circle &c1, double rad) {
+		// circle tanto cline & circle with radius
+		CLine Offs0	= Parallel(AT0, s0, rad);
+		Circle c2 = c1;
+		c2.radius += AT1 * rad;
+		Point pc = Intof(NF, Offs0, c2);
+		return (pc.ok)? Circle(pc, rad) : INVALID_CIRCLE;
+	}
+
+	Circle	Tanto( int LR, int AT0, const Circle& c0, const Point& p, double rad) {
+		// circle tanto circle & thro' a point
+		Circle c2 = c0;
+		c2.radius += AT0 * rad;
+		Circle c1(p, rad);
+		Point pc = Intof(LR, c2, c1);
+		return (pc.ok)? Circle(pc, rad) : INVALID_CIRCLE;
+	}
+	Circle	Tanto(int LR, int AT0, const Circle& c0, int AT1, const Circle& c1, double rad) {
+		// circle tanto 2 circles
+		Circle c2 = c0;
+		Circle c3 = c1;
+		c2.radius += AT0 * rad;
+		c3.radius += AT1 * rad;
+		Point pc = Intof(LR, c2, c3);
+		return (pc.ok)? Circle(pc, rad) : INVALID_CIRCLE;
+	}
+
+	Circle Parallel(int side, const Circle& c0, double distance) {
+		// parallel to circle by distance
+		return Circle(c0.pc, c0.radius + (double) side * distance);
+	}
+
+	// distance
+	double atn360(double dy, double dx) {
+		// angle 0 to 2pi
+		double ang = atan2(dy, dx);
+		return ((ang < 0)? 2 * PI + ang : ang);
+	}
+
+	double Dist(const Point& p0, const Circle& c, const Point& p1) {
+		// clockwise distance around c from p0 to p1 
+		double a0 = atn360(p0.y - c.pc.y, p0.x - c.pc.x);
+		double a1 = atn360(p1.y - c.pc.y ,p1.x - c.pc.x);
+		if ( a1 > a0 ) a1 -= 2 * PI ;
+		return (a0 - a1) * c.radius;
+	}
+	double Dist(const CLine& s, const Circle& c) {
+		// distance between line and circle
+		return fabs(s.Dist(c.pc)) - c.radius;
+	}
+	double Dist(const Circle& c0, const Circle& c1) {
+		// distance between 2 circles
+		return c0.pc.Dist(c1.pc) - c0.radius - c1.radius;
+	}
+	double Dist(const Circle& c, const Point& p) {
+		// distance between circle and point
+		return p.Dist(On(c, p));
+	}
+
+	double IncludedAngle(const Vector2d& v0, const Vector2d& 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. - UNIT_VECTOR_TOLERANCE) return 0;
+		if(inc_ang < -1. + UNIT_VECTOR_TOLERANCE)
+			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 IncludedAngle(const Vector3d& v0, const Vector3d& v1, const Vector3d& normal, int dir) {
+		// returns the absolute included angle between 2 vectors in the direction of dir ( 1=acw  -1=cw) about normal
+		double inc_ang = v0 * v1;
+
+		if(inc_ang >= -NEARLY_ONE) {									// dot product,   v1 . v2  =  cos ang
+			inc_ang = acos(inc_ang);									// 0 to pi radians
+
+			if(dir * (normal * (v0 ^ v1)) < 0) inc_ang = 2 * PI - inc_ang ;		// cp
+		}
+		else
+			inc_ang = PI;	// semi-cicle
+
+		return dir * inc_ang;
+	}
+
+	int corner(const Vector2d& v0, const Vector2d& v1, double cpTol) {
+		// returns corner
+		//						0 (TANGENT) = tangent
+		//						1 (LEFT)    = left turn
+		//					   -1 (RIGHT)   = right turn
+		double cp = v0 ^ v1;
+		if(fabs(cp) < cpTol) return TANGENT;
+
+		return (cp > 0)?GEOFF_LEFT : GEOFF_RIGHT;
+	}
+
+	int quadratic(double a, double b, double c, double& x0, double& x1) {
+		// solves quadratic equation ax² + bx + c = 0
+		// returns number of real roots
+//		double epsilon = 1.0e-6;
+		double epsilon = (geoff_geometry::UNITS == METRES)?1.0e-09 : 1.0e-06;
+		double epsilonsq = epsilon * epsilon;
+		if(fabs(a) < epsilon) {
+			if(fabs(b) < epsilon) return 0;		// invalid
+			x0 = - c / b;
+			return 1;
+		}
+		b /= a;
+		c /= a;
+		double s = b * b - 4 * c;
+		if(s < -epsilon) return 0;				// imaginary roots
+		x0 = - 0.5 * b;
+		if(s > epsilonsq) {
+			s = 0.5 * sqrt(s);
+			x1 = x0 - s;
+			x0 += s;
+			return 2;
+		}
+		return 1;
+	}
+
+	Plane::Plane(const Point3d& p0, const Point3d& p1, const Point3d& p2) {
+		// constructor plane from 3 points
+		normal = Vector3d(p0, p1) ^ Vector3d(p0, p2);
+		normal.normalise();
+		ok = (normal != NULL_VECTOR);
+		d = -(normal * Vector3d(p0));
+	}
+
+	Plane::Plane(const Point3d& p0, const Vector3d& v, bool normalise) {
+		// constructor plane from point & vector
+		normal = v;
+		if(normalise == true) normal.normalise();
+		d = -(normal * Vector3d(p0));
+	}
+
+	Plane::Plane(double dist, const Vector3d& n) {
+		normal = n;
+		double mag = normal.normalise();
+		if((ok = (normal != NULL_VECTOR))) d = dist / mag;
+	}
+
+	double Plane::Dist(const Point3d& p)const{
+		// returns signed distance to plane from point p
+	return (normal * Vector3d(p)) + d;
+}
+
+	Point3d Plane::Near(const Point3d& p)const {
+		// returns near point to p on the plane
+		return - normal * Dist(p) + p;
+	}
+
+	bool Plane::Intof(const Line& l, Point3d& intof, double& t) const{
+		// intersection between plane and line
+		// input this plane, line
+		// output intof
+		// method returns true for valid intersection
+		double den = l.v * this->normal;
+		if(fabs(den) < UNIT_VECTOR_TOLERANCE)	return false; // line is parallel to the plane, return false, even if the line lies on the plane
+
+		t = -(normal * Vector3d(l.p0) + d) / den;
+		intof = l.v * t + l.p0;
+		return true;
+	}
+
+	bool Plane::Intof(const Plane& pl, Line& intof)const {
+		// intersection of 2 planes
+		Vector3d d = this->normal ^ pl.normal;
+		d.normalise();
+		intof.ok = false;
+		if(d == NULL_VECTOR) return false;		// parallel planes
+
+		intof.v = d;
+		intof.length = 1;
+		
+		double dot = this->normal * pl.normal;
+
+		double den = dot * dot - 1.;
+		double a = (this->d - pl.d * dot) / den;
+		double b = (pl.d - this->d * dot) / den;
+		intof.p0 = a * this->normal + b * pl.normal;
+		intof.ok = true;
+		return true;
+	}
+
+	bool Plane::Intof(const Plane& pl0, const Plane& pl1, Point3d& intof) const{
+		// intersection of 3 planes
+		Line tmp;
+		if(Intof(pl0, tmp)) {
+			double t;
+			return pl1.Intof(tmp, intof, t);
+		}
+		return false;
+	}
+}

src/Mod/Path/libarea/kurve/Finite.cpp

@@ -0,0 +1,654 @@
+// written by g.j.hawkesford 2006 for Camtek Gmbh
+//
+// This program is released under the BSD license. See the file COPYING for details.
+//
+
+#include "geometry.h"
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+//            finite intersections
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+#ifndef WIN32
+#define __min(a,b) ((a<b)?a:b)
+#define __max(a,b) ((a>b)?a:b)
+#endif
+
+namespace geoff_geometry {
+	int Intof(const Span& sp0, const Span& sp1, Point& p0, Point& p1, double t[4])
+	{
+		// returns the number of intersects (lying within spans sp0, sp1)
+		if(sp0.box.outside(sp1.box) == true) return 0;
+		if(!sp0.dir) {
+			if(!sp1.dir) {
+				// line line
+				return LineLineIntof(sp0, sp1, p0, t);
+			}
+			else {
+				// line arc
+				return LineArcIntof(sp0, sp1, p0, p1, t);
+			}
+		}
+		else {
+			if(!sp1.dir) {
+				// arc line
+				return LineArcIntof(sp1, sp0, p0, p1, t);
+			}
+			else {
+				// arc arc
+				return ArcArcIntof(sp0, sp1, p0, p1, t);
+			}
+		}
+	}
+
+	int LineLineIntof(const Span& sp0 , const Span& sp1, Point& p, double t[2]) {
+		// intersection between 2 Line2d
+		// returns 0 for no intersection in range of either span
+		// returns 1 for intersction in range of both spans
+		// t[0] is parameter on sp0,
+		// t[1] is parameter on sp1
+		Vector2d v0(sp0.p0, sp0.p1);
+		Vector2d v1(sp1.p0, sp1.p1);
+		Vector2d v2(sp0.p0, sp1.p0);
+
+		double cp = v1 ^ v0;
+
+		if(fabs(cp) < UNIT_VECTOR_TOLERANCE ) {
+			p = INVALID_POINT;
+			return 0;					// parallel or degenerate lines
+		}
+
+		t[0] = (v1 ^ v2) / cp;
+		p = v0 * t[0] + sp0.p0;
+		p.ok = true;
+		double toler = geoff_geometry::TOLERANCE / sp0.length;				// calc a parametric tolerance
+
+		t[1] = (v0 ^ v2) / cp;
+		if(t[0] < -toler || t[0] > 1 + toler) return 0;						// intersection on first?
+		toler = geoff_geometry::TOLERANCE / sp1.length;						// calc a parametric tolerance
+		if(t[1] < -toler || t[1] > 1 + toler) return 0;						// intersection on second?
+		return 1;
+	}
+
+	int LineArcIntof(const Span& line, const Span& arc, Point& p0, Point& p1, double t[4]) {
+		// inters of line arc
+		// solving	x = x0 + dx * t			x = y0 + dy * t
+		//			x = xc + R * cos(a)		y = yc + R * sin(a)		for t
+		// gives :-  t² (dx² + dy²) + 2t(dx*dx0 + dy*dy0) + (x0-xc)² + (y0-yc)² - R² = 0
+		int nRoots;
+		Vector2d v0(arc.pc, line.p0);
+		Vector2d v1(line.p0, line.p1);
+		double s = v1.magnitudesqd();
+
+		p0.ok = p1.ok = false;
+		if((nRoots = quadratic(s, 2 * (v0 * v1), v0.magnitudesqd() - arc.radius * arc.radius, t[0], t[1])) != 0) {
+			double toler = geoff_geometry::TOLERANCE / sqrt(s);							// calc a parametric tolerance
+			if(t[0] > -toler && t[0] < 1 + toler) {
+				p0 = v1 * t[0] + line.p0;
+				p0.ok = arc.OnSpan(p0, &t[2]);
+			}
+			if(nRoots == 2) {
+				if(t[1] > -toler && t[1] < 1 + toler) {
+					p1 = v1 * t[1] + line.p0;
+					p1.ok = arc.OnSpan(p1, &t[3]);
+				}
+			}
+			if(!p0.ok && p1.ok) {
+				p0 = p1;
+				p1.ok = false;
+			}
+			nRoots = (int)p0.ok + (int)p1.ok;
+		}
+		return nRoots;
+	}
+
+	int ArcArcIntof(const Span& arc0, const Span& arc1, Point& pLeft, Point& pRight, double t[4]) {
+		// Intof 2 arcs
+		int numInts = Intof(Circle(arc0.pc, arc0.radius), Circle(arc1.pc, arc1.radius), pLeft, pRight);
+
+		if(numInts == 0) {
+			pLeft = arc0.p1;
+			pLeft.ok = false;
+			return 0;
+		}
+		int nLeft  = arc0.OnSpan(pLeft) && arc1.OnSpan(pLeft);
+		int nRight = (numInts == 2)?arc0.OnSpan(pRight) && arc1.OnSpan(pRight) : 0;
+		if(nLeft == 0 && nRight) pLeft = pRight;
+		return nLeft + nRight;
+	}
+
+	bool Span::OnSpan(const Point& p)const {
+		double t;
+		return OnSpan(p, &t);
+	}
+
+	bool Span::OnSpan(const Point& p, double* t)const {
+		// FAST OnSpan test - assumes that p lies ON the unbounded span
+#if _DEBUG
+		if(this->returnSpanProperties == false) {
+			FAILURE(L"OnSpan - properties no set, incorrect calling code");
+		}
+#endif
+#if 0
+		if(NullSpan) {
+			*t = 0.0;
+			return (p == p0);
+		}
+
+		if(p == p0) {
+			*t = 0.0;
+			return true;
+		}
+
+		if(p == p1) {
+			*t = 1.0;
+			return true;
+		}
+#endif
+		bool ret;
+//		if(p == this->p0 || p == this->p1) return true;
+
+		if(dir == LINEAR) {
+#if 1
+#if _DEBUG
+			// check p is on line
+			CLine cl(*this);
+			double d = fabs(cl.Dist(p));
+			if( d > geoff_geometry::TOLERANCE) {
+				FAILURE(L"OnSpan - point not on linear span, incorrect calling code");
+			}
+#endif
+#endif
+			Vector2d v0(p0, p);
+			*t = vs * v0;
+//			ret = (*t > - geoff_geometry::TOLERANCE && *t < length + geoff_geometry::TOLERANCE);
+
+			*t = *t / length;
+			ret = (*t >= 0 && *t <= 1.0 );
+
+		}
+		else {
+			// true if p lies on arc span sp (p must be on circle of span)
+#if 1
+#if _DEBUG
+			// check that p lies on the arc
+			double d = p.Dist(pc);
+			if(FNE(d, radius, geoff_geometry::TOLERANCE)) {
+				FAILURE(L"OnSpan - point not on circular span, incorrect calling code");
+			}
+
+#endif
+#endif
+#if 0	// alt method (faster, but doesn't provide t)
+			Vector2d v0(p0, p);
+			Vector2d v1(p0, p1);
+
+			// check angle to point from start
+			double cp;
+			ret = ((cp = (dir * (v0 ^ v1))) > 0);
+			*t = 0.0;// incorrect !!!
+#else
+			Vector2d v = ~Vector2d(pc, p);
+			v.normalise();
+			if(dir == CW) v = -v;
+
+			double ang = IncludedAngle(vs, v, dir);
+			*t = ang / angle;
+			ret = (*t >= 0 && *t <= 1.0);
+#endif
+		}
+
+		return ret;
+	}
+
+	Line::Line(const Point3d& p, const Vector3d& v0, bool boxed){
+		// constructor from point & vector
+		p0 = p;
+		v = v0;
+		length = v.magnitude();
+		if(boxed) minmax();
+		ok = (length > geoff_geometry::TOLERANCE);
+	}
+
+	Line::Line(const Point3d& p, const Point3d& p1){
+		// constructor from 2 points
+		p0 = p;
+		v = Vector3d(p, p1);
+		length = v.magnitude();
+		minmax();
+		ok = (length > geoff_geometry::TOLERANCE);
+	}
+
+	Line::Line(const Span& sp){
+		// contructor from linear span
+		p0 = sp.p0;
+		v = sp.vs * sp.length;
+		length = sp.length;
+		//	box = sp.box;
+		box.min = Point3d(sp.box.min);
+		box.max = Point3d(sp.box.max);
+		ok = !sp.NullSpan;
+	}
+
+	void Line::minmax() {
+		MinMax(this->p0, box.min, box.max);
+		MinMax(this->v + this->p0, box.min, box.max);
+	}
+
+	bool Line::atZ(double z, Point3d& p)const {
+		// returns p at z on line
+		if(FEQZ(this->v.getz())) return false;
+		double t = (z - this->p0.z) / this->v.getz();
+		p = Point3d(this->p0.x + t * this->v.getx(), this->p0.y + t * this->v.gety(), z);
+		return true;
+	}
+
+
+	bool Line::Shortest(const Line& l2, Line& lshort, double& t1, double& t2)const {
+	/*
+	Calculate the line segment PaPb that is the shortest route between
+	two lines P1P2 and P3P4. Calculate also the values of mua and mub where
+	Pa = P1 + t1 (P2 - P1)
+	Pb = P3 + t2 (P4 - P3)
+	Return FALSE if no solution exists.       P Bourke method.
+		Input this 1st line
+		Input l2   2nd line
+		Output lshort shortest line between lines (if lshort.ok == false, the line intersect at a point lshort.p0)
+		Output t1 parameter at intersection on 1st Line
+		Output t2 parameter at intersection on 2nd Line
+
+	*/
+		Vector3d v13(l2.p0, this->p0);
+		if(this->ok == false || l2.ok == false) return false;
+
+		double d1343 = v13 * l2.v;		// dot products
+		double d4321 = l2.v * this->v;
+		double d1321 = v13 * this->v;
+		double d4343 = l2.v * l2.v;
+		double d2121 = this->v * this->v;
+
+		double denom = d2121 * d4343 - d4321 * d4321;
+		if(fabs(denom) < 1.0e-09) return false;
+		double numer = d1343 * d4321 - d1321 * d4343;
+
+		t1 = numer / denom;
+		t2 = (d1343 + d4321 * t1) / d4343;
+
+		lshort = Line(t1* this->v + this->p0, t2 * l2.v + l2.p0);
+		t1 *= this->length;
+		t2 *= l2.length;		// parameter in line length for tolerance checking
+		return true;
+	}
+
+	int Intof(const Line& l0, const Line& l1, Point3d& intof)
+	{
+		/* intersection of 2 vectors
+		returns 0 for  intercept but not within either vector
+		returns 1 for intercept on both vectors
+
+		note that this routine always returns 0 for parallel vectors
+		method:
+		x = x0 + dx0 * t0	for l0
+		...
+		...
+		x = x1 + dx1 * t1	for l1
+		...
+		...
+
+		x0 + dx0 * t0 = x1 + dx1 * t1
+		dx0 * t0 - dx1 * t1 + x0 - x1 = 0
+
+		setup 3 x 3 determinent for 
+		a0 t0 + b0 t1 + c0 = 0
+		a1 t0 + b1 t1 + c1 = 0
+		a2 t0 + b2 t1 + c2 = 0
+
+		from above a = l0.v
+		b = -l1.v
+		c = Vector3d(l1, l0)
+		*/
+		//	Vector3d a = l0.v;
+		if(l0.box.outside(l1.box) == true) return 0;
+		Vector3d b = -l1.v;
+		Vector3d c = Vector3d(l1.p0, l0.p0);
+		Vector3d det = l0.v ^ b;
+		Vector3d t = b ^ c;
+
+		// choose largest determinant & corresponding parameter for accuracy
+		double t0 = t.getx();
+		double d  = det.getx();
+
+		if(fabs(det.getz()) > fabs(det.gety())) {
+			if(fabs(det.getz()) > fabs(det.getx())) {
+				t0 = t.getz();
+				d = det.getz();
+			}
+		}
+		else {
+			if(fabs(det.gety()) > fabs(det.getx())) {
+				t0 = t.gety();
+				d = det.gety();
+			}
+		}
+
+		if(fabs(d) < 1.0e-06) return 0;
+
+		t0 /= d;
+		intof = l0.v * t0 + l0.p0;
+
+		Point3d other;
+		double t1;
+		if(Dist(l1, intof, other, t1) > geoff_geometry::TOLERANCE) return 0;
+
+		t0 *= l0.length;
+		if( t0 < -geoff_geometry::TOLERANCE || t0 > l0.length + geoff_geometry::TOLERANCE || t1 < -geoff_geometry::TOLERANCE || t1 > l1.length + geoff_geometry::TOLERANCE ) return 0;
+		return 1;
+	}
+
+
+	double Dist(const Line& l, const Point3d& p, Point3d& pnear, double& t){
+		// returns the distance of a point from a line and the near point on the extended line and the parameter of the near point (0-length) in range
+		pnear = Near(l, p, t );
+		return p.Dist(pnear);
+	}
+
+	Point3d Near(const Line& l, const Point3d& p, double& t){
+		// returns the near point from a line on the extended line and the parameter of the near point (0-length) in range
+		t = (Vector3d(l.p0, p) * l.v) / l.length;		// t parametised 0 - line length
+		return l.v * (t / l.length) + l.p0;
+	}
+
+	Point3d Line::Near(const Point3d& p, double& t)const{
+		// returns the near point from a line on the extended line and the parameter of the near point (0-length) in range
+		t = (Vector3d(this->p0, p) * this->v) / this->length;		// t parametised 0 - line length
+		return this->v * (t / this->length) + this->p0;
+	}
+
+	double DistSq(const Point3d *p, const Vector3d *vl, const Point3d *pf) {
+		/// returns the distance squared of pf from the line given by p,vl
+		/// vl must be normalised
+		Vector3d v(*p, *pf);
+		Vector3d vcp = *vl ^ v;
+		double d = vcp.magnitudeSq(); // l * sina
+		return d;
+	}
+
+	double Dist(const Point3d *p, const Vector3d *vl, const Point3d *pf) {
+		/// returns the distance of pf from the line given by p,vl
+		/// vl must be normalised
+		Vector3d v(*p, *pf);
+		Vector3d vcp = *vl ^ v;
+		double d = vcp.magnitude(); // l * sina
+		return d;
+#if 0
+		// slower method requires 2 sqrts
+		Vector3d v(*p, *pf);
+		double magv = v.normalise();
+		Vector3d cp = *vl ^ v;
+		double d = magv * cp.magnitude();
+		return d;  // l * sina
+#endif
+	}
+
+	double Dist(const Span& sp, const Point& p , Point& pnear ) {
+		// returns distance of p from span, pnear is the nearpoint on the span (or endpoint)
+		if(!sp.dir) {
+			double d, t;
+			Point3d unused_pnear;
+			d = Dist(Line(sp), Point3d(p), unused_pnear, t);
+			if(t < -geoff_geometry::TOLERANCE) {
+				pnear = sp.p0;						// nearpoint
+				d = pnear.Dist(p);
+			}
+			else if(t > sp.length + geoff_geometry::TOLERANCE) {
+				pnear = sp.p1;
+				d = pnear.Dist(p);
+			}
+			return d;
+		}
+		else {
+			// put pnear on the circle
+			double radiusp;
+			Vector2d v(sp.pc, p);
+			if((radiusp = v.magnitude()) < geoff_geometry::TOLERANCE)	{
+				// point specified on circle centre - use first point as near point
+				pnear = sp.p0;						// nearpoint
+				return sp.radius;
+			}
+			else	{
+				pnear = v * (sp.radius / radiusp) + sp.pc;
+
+				// check if projected point is on the arc
+				if(sp.OnSpan(pnear)) return fabs(radiusp - sp.radius);
+				// double      h1 = pnear.x - sp.p0.x ;
+				// double      v1 = pnear.y - sp.p0.y ;
+				// double      h2 = sp.p1.x - pnear.x ;
+				// double      v2 = sp.p1.y - pnear.y ;
+				//       if ( sp.dir * ( h1 * v2 - h2 * v1 ) >= 0 )return fabs(radiusp - sp.radius);
+
+				// point not on arc so calc nearest end-point
+				double ndist = p.Dist(sp.p0);
+				double dist = p.Dist(sp.p1);
+				if(ndist >=  dist) {
+					// sp.p1 is near point
+					pnear = sp.p1;
+					return dist;
+				}
+
+				// sp.p0 is near point
+				pnear = sp.p0;						// nearpoint
+				return ndist ;
+			}
+		}
+	}
+
+	bool	OnSpan(const Span& sp, const Point& p) {
+		Point nullPoint;
+		return OnSpan(sp, p, false, nullPoint, nullPoint);
+	}
+
+	bool	OnSpan(const Span& sp, const Point& p,  bool nearPoints, Point& pNear, Point& pOnSpan) {
+		// function returns true if pNear == pOnSpan
+		//			returns pNear & pOnSpan if nearPoints true
+		//			pNear (nearest on unbound span)
+		//			pOnSpan (nearest on finite span)
+		if(sp.dir) {
+			// arc
+			if(fabs(p.Dist(sp.pc) - sp.radius) > geoff_geometry::TOLERANCE) {
+				if(!nearPoints) return false;
+			}
+
+			pNear = On(Circle(sp.pc, sp.radius), p);
+
+			if(sp.OnSpan(pNear)) {
+				if(nearPoints) pOnSpan = pNear;
+				return true; // near point is on arc - already calculated
+			}
+
+			// point not on arc return the nearest end-point
+			if(nearPoints) pOnSpan = (p.Dist(sp.p0) >= p.Dist(sp.p1)) ?sp.p1 : sp.p0;
+			return false;
+		}
+		else {
+			// straight
+			if(fabs(CLine(sp.p0, sp.vs).Dist(p)) > geoff_geometry::TOLERANCE) {
+				if(!nearPoints) return false;
+			}
+			Vector2d v(sp.p0, p);
+			double t = v * sp.vs;
+			if(nearPoints) pNear = sp.vs * t + sp.p0;
+			bool onSpan = (t > - geoff_geometry::TOLERANCE && t < sp.length + geoff_geometry::TOLERANCE);
+			if(! onSpan) {
+				if(nearPoints) pOnSpan = (p.Dist(sp.p0) >= p.Dist(sp.p1))?sp.p1 : sp.p0;
+			}
+			else {
+				if(nearPoints) pOnSpan = pNear;
+			}
+			return onSpan;
+		}
+	}
+
+	// Triangle3d Constructors
+	Triangle3d::Triangle3d(const Point3d& p1, const Point3d& p2, const Point3d& p3) {
+		vert1 = p1;
+		vert2 = p2;
+		vert3 = p3;
+		v0 = Vector3d(vert1, vert2);
+		v1 = Vector3d(vert1, vert3);
+		ok = true;
+
+		// set box
+		box.min.x = __min(__min(vert1.x, vert2.x), vert3.x);
+		box.min.y = __min(__min(vert1.y, vert2.y), vert3.y);
+		box.min.z = __min(__min(vert1.z, vert2.z), vert3.z);
+
+		box.max.x = __max(__max(vert1.x, vert2.x), vert3.x);
+		box.max.y = __max(__max(vert1.y, vert2.y), vert3.y);
+		box.max.z = __max(__max(vert1.z, vert2.z), vert3.z);
+	}
+
+	// Triangle3d methods
+	bool    Triangle3d::Intof(const Line& l, Point3d& intof)const {
+	 // returns intersection triangle to line in intof
+	// funtion returns true for intersection, false for no intersection
+	// method based on Möller & Trumbore(1997) (Barycentric coordinates)
+	// based on incorrect Pseudo code from "Geometric Tools for Computer Graphics" p.487
+		if(box.outside(l.box) == true) return false;
+
+		Vector3d line(l.v);
+		line.normalise();
+
+		Vector3d p = line ^ v1;				// cross product
+		double tmp = p * v0;				// dot product
+
+		if(FEQZ(tmp)) return false;
+
+		tmp = 1 / tmp;
+		Vector3d s(vert1, l.p0);
+
+		double u = tmp * (s * p);			// barycentric coordinate
+		if(u < 0 || u > 1) return false;	// not inside triangle
+
+		Vector3d q = s ^ v0;
+		double v = tmp * (line * q);		// barycentric coordinate
+		if(v < 0 || v > 1) return false;	// not inside triangle
+
+		if( u + v > 1) return false;		// not inside triangle
+
+		double t = tmp * (v1 * q);
+		intof = line * t + l.p0;
+		return true;
+	}
+
+
+	// box class
+	bool Box::outside(const Box& b)const {
+		// returns true if this box is outside b
+		if(b.ok == false || this->ok == false) return false;	// no box set
+		if(this->max.x < b.min.x) return true;
+		if(this->max.y < b.min.y) return true;
+		if(this->min.x > b.max.x) return true;
+		if(this->min.y > b.max.y) return true;
+		return false;
+	}
+
+	void Box::combine(const Box& b) {
+		if(b.max.x > this->max.x) this->max.x = b.max.x;
+		if(b.max.y > this->max.y) this->max.y = b.max.y;
+		if(b.min.x < this->min.x) this->min.x = b.min.x;
+		if(b.min.y < this->min.y) this->min.y = b.min.y;
+	}
+
+	void Box3d::combine(const Box3d& b) {
+		if(b.max.x > this->max.x) this->max.x = b.max.x;
+		if(b.max.y > this->max.y) this->max.y = b.max.y;
+		if(b.max.z > this->max.z) this->max.z = b.max.z;
+		if(b.min.x < this->min.x) this->min.x = b.min.x;
+		if(b.min.y < this->min.y) this->min.y = b.min.y;
+		if(b.min.z < this->min.z) this->min.z = b.min.z;
+	}
+
+	bool Box3d::outside(const Box3d& b) const{
+		// returns true if this box is outside b
+		if(b.ok == false || this->ok == false) return false;	// no box set
+		if(this->max.x < b.min.x) return true;
+		if(this->max.y < b.min.y) return true;
+		if(this->max.z < b.min.z) return true;
+		if(this->min.x > b.max.x) return true;
+		if(this->min.y > b.max.y) return true;
+		if(this->min.z > b.max.z) return true;
+		return false;
+	}
+#if 0
+	Span3d IsPtsSpan3d(const double* a, int n, double tolerance, double* deviation) {
+		// returns a span3d if all points are within tolerance
+		int np = n / 3;					// number of points
+		if(np < 2) return Span3d();		// Invalid span3d
+		Point3d sp = Point3d(&a[0]);
+		Point3d ep = Point3d(&a[n-3]);
+		Line line = IsPtsLine(a, n, tolerance, deviation);
+		if(line.ok) return Span3d(sp, ep);	// it's a line
+
+		*deviation = 0;					// cumulative deviation
+		Point3d mp = Point3d(&a[np / 2 * 3]);	// mid point
+		Plane plane(sp, mp, ep);
+		if(plane.ok) {
+			// plane of the arc is ok
+			// calculate centre point
+			Vector3d vs(mp, sp);
+			vs.normalise();
+			Vector3d ve(mp, ep);
+			ve.normalise();
+			Vector3d rs = vs ^ plane.normal;
+			Vector3d re = ve ^ plane.normal;
+
+			Line rsl(sp.Mid(mp), rs, false);
+			Line rel(ep.Mid(mp), re, false);
+
+			Point3d pc;
+			Intof(rsl, rel, pc);
+			double radius = pc.Dist(sp);
+
+			// check other points on circle
+			for(int i = 2; i < np - 1; i++) {
+				Point3d p(&a[i*3]);
+				double dp = fabs(plane.Dist(p));
+				double dr = fabs(p.Dist(pc) - radius);
+double tolerance = 10.0 * 1.0e-6;
+				if(dp > tolerance || dr > tolerance) {
+					return Span3d();
+				}
+
+			}
+			return Span3d(CW, plane.normal, sp, ep, pc);
+
+		}
+		return Span3d();
+	}
+#endif
+
+	Line IsPtsLine(const double* a, int n, double tolerance, double* deviation) {
+		// returns a Line if all points are within tolerance
+		// deviation is returned as the sum of all deviations of interior points to line(sp,ep)
+		int np = n / 3;					// number of points
+		*deviation = 0;					// cumulative deviation
+		if(np < 2) return Line();		// Invalid line
+
+		Point3d sp(&a[0]);
+		Point3d ep(&a[n-3]);
+		Line line(sp, ep);				// line start - end
+
+		if(line.ok) {
+			for(int j = 1; j < np - 1; j++) {
+				Point3d mp(&a[j * 3]);
+				double t, d=0;
+				if((d = mp.Dist(line.Near(mp, t))) > tolerance) {
+					line.ok = false;
+					return line;
+				}
+				*deviation = *deviation + d;
+			}
+		}
+		return line;
+	}
+}

src/Mod/Path/libarea/kurve/License.txt

@@ -0,0 +1,24 @@
+Copyright (c) <2015>, <Dan Heeks>
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in the
+      documentation and/or other materials provided with the distribution.
+    * Neither the name of the <organization> nor the
+      names of its contributors may be used to endorse or promote products
+      derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

src/Mod/Path/libarea/kurve/Matrix.cpp

@@ -0,0 +1,624 @@
+////////////////////////////////////////////////////////////////////////////////////////////////
+//                    3d geometry classes - implements some 3d stuff
+//
+//                    g.j.hawkesford August 2003
+//
+// This program is released under the BSD license. See the file COPYING for details.
+//
+////////////////////////////////////////////////////////////////////////////////////////////////
+
+#include "geometry.h"
+using namespace geoff_geometry;
+
+#ifdef PEPSDLL
+	#include "vdm.h"
+	#include "pepsdll.h"
+	#include "realds.h"
+#endif
+////////////////////////////////////////////////////////////////////////////////////////////////
+// matrix
+////////////////////////////////////////////////////////////////////////////////////////////////
+namespace geoff_geometry {
+
+	Matrix::Matrix(){
+		Unit();
+	}
+	Matrix::Matrix(double m[16]) {
+		memcpy(e, m, sizeof(e));
+		this->IsUnit();
+		this->IsMirrored();
+	}
+
+	Matrix::Matrix( const Matrix& m)
+	{
+		*this = m;
+	}
+
+	bool Matrix::operator==(const Matrix &m)const{
+		// m1 == m2
+		if(this->m_unit != m.m_unit || this->m_mirrored != m.m_mirrored) return false;
+		for(int i = 0; i < 16; i++)
+			if(FEQ(this->e[i], m.e[i], TIGHT_TOLERANCE) == false) return false;
+		return true;
+	}
+
+#if 0
+	const Matrix& Matrix::operator=( Matrix &m) {
+		for(int i = 0; i < 16; i++) e[i] = m.e[i];
+		m_unit = m.m_unit;
+		m_mirrored = m.m_mirrored;
+		return *this;
+	}
+#endif
+	void	Matrix::Unit()
+	{
+		// homogenous matrix - set as unit matrix
+		memset(e, 0, sizeof(e));
+		e[0] = e[5] = e[10] = e[15] = 1;
+		m_unit = true;
+		m_mirrored = false;
+	}
+
+	void	Matrix::Get(double* p) const
+	{
+		// copy the matrix
+		memcpy(p, e, sizeof(e));
+	}
+	void	Matrix::Put(double* p)
+	{
+		// assign the matrix
+		memcpy(e, p, sizeof(e));
+		m_unit = false;		// don't know
+		m_mirrored = -1;	// don't know
+
+	}
+	void	Matrix::Translate(double x, double y, double z)
+	{
+		// translation
+		e[3]  += x;
+		e[7]  += y;
+		e[11] += z;
+		m_unit = false;
+	}
+
+	void	Matrix::Rotate(double angle, Vector3d *rotAxis) {
+		/// Rotation about rotAxis with angle
+		Rotate(sin(angle), cos(angle), rotAxis);
+	}
+
+	void	Matrix::Rotate(double sinang, double cosang, Vector3d *rotAxis) {
+		/// Rotation about rotAxis with cp & dp
+		Matrix rotate;
+		double oneminusc = 1.0 - cosang;
+
+		rotate.e[0] = rotAxis->getx() * rotAxis->getx() * oneminusc + cosang;
+		rotate.e[1] = rotAxis->getx() * rotAxis->gety() * oneminusc - rotAxis->getz() * sinang;
+		rotate.e[2] = rotAxis->getx() * rotAxis->getz() * oneminusc + rotAxis->gety() * sinang;
+
+		rotate.e[4] = rotAxis->getx() * rotAxis->gety() * oneminusc + rotAxis->getz() * sinang;
+		rotate.e[5] = rotAxis->gety() * rotAxis->gety() * oneminusc + cosang;
+		rotate.e[6] = rotAxis->gety() * rotAxis->getz() * oneminusc - rotAxis->getx() * sinang;
+
+		rotate.e[8] = rotAxis->getx() * rotAxis->getz() * oneminusc - rotAxis->gety() * sinang;
+		rotate.e[9] = rotAxis->gety() * rotAxis->getz() * oneminusc + rotAxis->getx() * sinang;
+		rotate.e[10] = rotAxis->getz() * rotAxis->getz() * oneminusc  + cosang;
+		Multiply(rotate); // concatinate rotation with this matrix
+		m_unit = false;
+		m_mirrored = -1;	// don't know
+	}
+
+
+	void	Matrix::Rotate(double angle, int Axis)
+	{	// Rotation (Axis 1 = x , 2 = y , 3 = z 
+		Rotate(sin(angle), cos(angle), Axis);
+	}
+
+	void	Matrix::Rotate(double sinang, double cosang, int Axis)
+	{	// Rotation (Axis 1 = x , 2 = y , 3 = z 
+		Matrix rotate;
+		rotate.Unit();
+
+		switch(Axis)
+		{
+		case 1:
+			// about x axis
+			rotate.e[5] = rotate.e[10] = cosang;
+			rotate.e[6] = -sinang;
+			rotate.e[9] = sinang;
+			break;
+		case 2:
+			// about y axis
+			rotate.e[0] = rotate.e[10] = cosang;
+			rotate.e[2] = sinang;
+			rotate.e[8] = -sinang;
+			break;
+		case 3:
+			// about z axis
+			rotate.e[0] = rotate.e[5] = cosang;
+			rotate.e[1] = -sinang;
+			rotate.e[4] = sinang;
+			break;
+		}
+		Multiply(rotate); // concatinate rotation with this matrix
+		m_unit = false;
+		m_mirrored = -1;	// don't know
+	}
+
+	void	Matrix::Scale(double scale)
+	{
+		// add a scale
+		Scale(scale, scale, scale);
+	}
+
+	void	Matrix::Scale(double scalex, double scaley, double scalez)
+	{
+		// add a scale
+		Matrix temp;
+		temp.Unit();
+
+		temp.e[0]  = scalex;
+		temp.e[5]  = scaley;
+		temp.e[10] = scalez;
+		Multiply(temp);
+		m_unit = false;
+		m_mirrored = -1;	// don't know
+	}
+	void	Matrix::Multiply(Matrix& m)
+	{
+		// multiply this by give matrix - concatinate
+		int i, k, l;
+		Matrix ret;
+
+		for (i = 0; i < 16; i++)
+		{
+			l = i - (k = (i % 4));
+			ret.e[i] =  m.e[l] * e[k] + m.e[l+1] * e[k+4] + m.e[l+2] * e[k+8] + m.e[l+3] * e[k+12];
+		}
+		
+		*this = ret;
+		this->IsUnit();
+	}
+
+	void	Matrix::Transform(double p0[3], double p1[3]) const
+	{
+		// transform p0 thro' this matrix
+		if(m_unit)
+			memcpy(p1, p0, 3 * sizeof(double));
+		else {
+			p1[0] = p0[0] * e[0] + p0[1] * e[1] + p0[2] * e[2]  + e[3];
+			p1[1] = p0[0] * e[4] + p0[1] * e[5] + p0[2] * e[6]  + e[7];
+			p1[2] = p0[0] * e[8] + p0[1] * e[9] + p0[2] * e[10] + e[11];
+		}
+	}
+	void	Matrix::Transform2d(double p0[2], double p1[2]) const
+	{
+		// transform p0 thro' this matrix (2d only)
+		if(m_unit)
+			memcpy(p1, p0, 2 * sizeof(double));
+		else {
+			p1[0] = p0[0] * e[0] + p0[1] * e[1] + e[3];
+			p1[1] = p0[0] * e[4] + p0[1] * e[5] + e[7];
+		}
+	}
+
+	void	Matrix::Transform(double p0[3]) const
+	{
+		double p1[3];
+		if(!m_unit) {
+			Transform(p0, p1);
+			memcpy(p0, p1, 3 * sizeof(double));
+		}
+	}
+
+	int	Matrix::IsMirrored()
+	{
+		// returns true if matrix has a mirror
+		if(m_unit)
+			m_mirrored = false;
+		else if(m_mirrored == -1) {
+
+			m_mirrored =  ((e[0] * (e[5] * e[10] - e[6] * e[9])
+				- e[1] * (e[4] * e[10] - e[6] * e[8])
+				+ e[2] * (e[4] * e[9]  - e[5] * e[8])) < 0);
+		}
+		return m_mirrored;
+	}
+	int Matrix::IsUnit() {
+		// returns true if unit matrix
+		for(int i = 0; i < 16; i++) {
+			if(i == 0 || i == 5 || i == 10 || i == 15) {
+				if(e[i] != 1) return m_unit = false;
+			}
+			else {
+				if(e[i] != 0) return m_unit = false;
+			}
+		}
+		m_mirrored = false;
+		return m_unit = true;
+	}
+
+	void	Matrix::GetTranslate(double& x, double& y, double& z) const
+	{
+		// return translation
+		x = e[3];
+		y = e[7];
+		z = e[11];
+	}
+	void	Matrix::GetScale(double& sx, double& sy, double& sz) const
+	{
+		// return the scale
+		if(m_unit) {
+			sx = sy = sz = 1;
+		}
+		else {
+			sx = sqrt(e[0] * e[0] + e[1] * e[1] + e[2] * e[2]);
+			sy = sqrt(e[4] * e[4] + e[5] * e[5] + e[6] * e[6]);
+			sz = sqrt(e[8] * e[8] + e[9] * e[9] + e[10] * e[10]);
+		}
+	}
+	bool	Matrix::GetScale(double& sx) const
+	{
+		// return a uniform scale  (false if differential)
+		double sy, sz;
+		if(m_unit) {
+			sx = 1;
+			return true;
+		}
+		GetScale(sx, sy, sz);
+		return (fabs(fabs(sx) - fabs(sy)) < 0.000001)?true : false;
+	}
+	void	Matrix::GetRotation(double& ax, double& ay, double& az) const
+	{
+		// return the rotations
+		if(m_unit) {
+			ax = ay = az = 0;
+			return;
+		}
+		double    a; /* cos(bx) */
+		double    b; /* sin(bx) */
+		double    c; /* cos(by) */
+		double    d; /* sin(by) */
+		double    ee; /* cos(bz) */
+		double    f; /* sin(bz) */
+		double sx, sy, sz;
+		GetScale(sx, sy, sz);
+		if(this->m_mirrored == -1) FAILURE(L"Don't know mirror - use IsMirrored method on object");
+		if(this->m_mirrored) sx = -sx;
+
+		// solve for d and decide case and solve for a, b, c, e and f
+		d = - e[8] / sz;
+		if((c = (1 - d) * (1 + d)) > 0.001)
+		{
+			// case 1
+			c = sqrt( c );
+			a = e[10] / sz / c;
+			b = e[9]  / sz / c;
+			ee = e[0]  / sx / c;
+			f = e[4]  / sy / c;
+		}
+		else
+		{
+			// case 2
+			double   coef;
+			double   p, q;
+
+			d = ( d < 0 ) ? -1 : 1 ;
+			c = 0 ;
+			p = d * e[5] / sy - e[2] / sx;
+			q = d * e[6] / sy + e[1] / sx;
+			if((coef = sqrt( p * p + q * q )) > 0.001) {
+				a = q / coef;
+				b = p / coef;
+				ee = b;
+				f = -d * b;
+			}
+			else
+			{
+				/* dependent pairs */
+				a =  e[5] / sy;
+				b = -e[6] / sy;
+				ee = 1 ;
+				f = 0 ;
+			}
+		}
+
+		// solve and return ax, ay and az
+		ax = atan2( b, a );
+		ay = atan2( d, c );
+		az = atan2( f, ee );
+	}
+
+	Matrix	Matrix::Inverse()
+	{
+		// matrix inversion routine
+
+		// a is input matrix destroyed & replaced by inverse
+		// method used is gauss-jordan (ref ibm applications)
+
+		double  hold , biga ;
+		int i , j , k , nk , kk , ij , iz ;
+		int ki , ji , jp , jk , kj , jq , jr , ik;
+
+		int n = 4;			// 4 x 4 matrix only
+		Matrix a = *this;
+		int l[4], m[4];
+
+		if(a.m_unit) return a;	// unit matrix
+
+		// search for largest element
+		nk =  - n ;
+		for ( k = 0 ; k < n ; k++ ) {
+			nk += n ;
+			l [ k ]  = m [ k ] = k ;
+			kk = nk + k ;
+			biga = a.e[ kk ]  ;
+
+			for ( j = k ; j < n ; j++ ) {
+				iz = n * j  ;
+				for ( i = k ; i < n ; i++ ) {
+					ij = iz + i ;
+					if ( fabs ( biga )  < fabs ( a.e[ ij ]  )  ) {
+						biga = a.e[ ij ]  ;
+						l[ k ]  = i ;
+						m[ k ] = j ;
+					}
+				}
+			}
+
+
+			// interchange rows
+			j = l[ k ]  ;
+			if ( j > k ) {
+				ki = k - n ;
+
+				for ( i = 0 ; i < n ; i++ ) {
+					ki += n ;
+					hold =  - a.e[ ki ]  ;
+					ji = ki - k + j ;
+					a.e[ ki ]  = a.e[ ji ]  ;
+					a.e[ ji ]  = hold ;
+				}
+			}
+
+			// interchange columns
+			i = m[ k ]  ;
+			if ( i > k ) {
+				jp = n * i ;
+				for ( j = 0 ; j < n ; j++ ) {
+					jk = nk + j ;
+					ji = jp + j ;
+					hold =  - a.e[ jk ]  ;
+					a.e[ jk ]  = a.e[ ji ]  ;
+					a.e[ ji ]  = hold ;
+				}
+			}
+
+			// divide columns by minus pivot (value of pivot element is contained in biga)
+			if ( fabs ( biga )  < 1.0e-10 )FAILURE(getMessage(L"Singular Matrix - Inversion failure",GEOMETRY_ERROR_MESSAGES, -1));	// singular matrix
+
+			for ( i = 0 ; i < n ; i++ ) {
+				if ( i != k ) {
+					ik = nk + i ;
+					a.e[ ik ]  =  - a.e[ ik ] /biga ;
+				}
+			}
+
+			// reduce matrix
+			for ( i = 0  ; i < n ; i++ ) {
+				ik = nk + i ;
+				hold = a.e[ ik ]  ;
+				ij = i - n ;
+
+				for ( j = 0 ; j < n ; j++ ) {
+					ij = ij + n ;
+					if ( i != k && j != k ) {
+						kj = ij - i + k ;
+						a.e[ ij ] = hold * a.e[ kj ]  + a.e[ ij ]  ;
+					}
+				}
+			}
+
+			// divide row by pivot
+			kj = k - n ;
+			for ( j = 0 ; j < n ; j++ ) {
+				kj = kj + n ;
+				if ( j != k )  a.e[ kj]  = a.e[ kj ] /biga ;
+			}
+
+			// replace pivot by reciprocal
+			a.e[ kk ] = 1 / biga ;
+		}
+
+		// final row and column interchange
+		k = n - 1 ;
+
+		while ( k > 0 ) {
+			i = l[ --k ]  ;
+			if ( i > k ) {
+				jq = n * k ;
+				jr = n * i ;
+
+				for ( j = 0 ; j < n ; j++ ) {
+					jk = jq + j ;
+					hold = a.e[jk]  ;
+					ji = jr + j ;
+					a.e[jk]  =  - a.e[ji]  ;
+					a.e[ji]  = hold ;
+				}
+			}
+
+			j = m[ k ]  ;
+			if ( j > k ) {
+				ki = k - n ;
+
+				for ( i = 1 ; i <= n ; i ++ ) {
+					ki = ki + n ;
+					hold = a.e[ ki ]  ;
+					ji = ki - k + j ;
+					a.e[ ki ] =  - a.e[ ji ]  ;
+					a.e[ ji ]  = hold ;
+				}
+			}
+		}
+
+		return a;
+	}
+
+#ifdef PEPSDLL
+	void Matrix::ToPeps(int id)
+	{
+		int set = PepsVdmMake(id, VDM_MATRIX_TYPE , VDM_LOCAL);
+		if(set < 0) FAILURE(L"Failed to create Matrix VDM");
+		struct kgm_header pepsm;
+
+		Get(pepsm.matrix);
+		pepsm.off_rad = 0;
+		pepsm.off_dir = pepsm.origin_id = 0;
+
+		PepsVdmWriteTmx(set , &pepsm );
+
+		PepsVdmClose(set);
+
+	}
+
+	void Matrix::FromPeps(int id)
+	{
+		//	if(id) {
+		int set = PepsVdmOpen(id, VDM_MATRIX_TYPE , VDM_READ_ONLY | VDM_LOCAL);
+		if(set < 0) FAILURE(L"Failed to open Matrix VDM");
+
+		struct kgm_header pepsm;
+		PepsVdmReadTmx(set , &pepsm);
+		memcpy(e, pepsm.matrix, sizeof(pepsm.matrix));
+		m_unit = true;
+		for(int i = 0; i < 16; i++) {
+			// copy over matrix and check for unit matrix
+			if(i == 0 || i == 5 || i == 10 || i == 15) {
+				if((e[i] = pepsm.matrix[i]) != 1) m_unit = false;
+			}
+			else {
+				if((e[i] = pepsm.matrix[i]) != 0) m_unit = false;
+			}
+		}
+		PepsVdmClose(set);
+		m_mirrored = IsMirrored();
+		//	}
+	}
+#endif
+
+	Matrix UnitMatrix;					// a global unit matrix
+
+
+	// vector
+	 Vector2d::Vector2d(const Vector3d &v){
+		if(FEQZ(v.getz())) FAILURE(L"Converting Vector3d to Vector2d illegal");
+		dx = v.getx();
+		dy = v.gety();
+	}
+
+	 bool Vector2d::operator==(const Vector2d &v)const {
+		return FEQ(dx, v.getx(), 1.0e-06) && FEQ(dy, v.gety(), 1.0e-06);
+	 }
+
+	 void Vector2d::Transform(const Matrix& m) {
+		 // transform vector
+		 if(m.m_unit == false) {
+			double dxt = dx * m.e[0] + dy * m.e[1];
+			double dyt = dx * m.e[4] + dy * m.e[5];
+			dx = dxt;
+			dy = dyt;
+		}
+		 this->normalise();
+	 }
+
+	 void Vector3d::Transform(const Matrix& m) {
+		 // transform vector
+		 if(m.m_unit == false) {
+			double dxt = dx * m.e[0] + dy * m.e[1] + dz * m.e[2];
+			double dyt = dx * m.e[4] + dy * m.e[5] + dz * m.e[6];
+			double dzt = dx * m.e[8] + dy * m.e[9] + dz * m.e[10];
+			dx = dxt;
+			dy = dyt;
+			dz = dzt;
+		}
+		 this->normalise();
+	 }
+
+	 void Vector3d::arbitrary_axes(Vector3d& x, Vector3d& y){
+		// arbitrary axis algorithm - acad method of generating an arbitrary but
+		// consistant set of axes from a single normal ( z )
+		// arbitrary x & y axes
+
+		if ( ( fabs ( this->getx() ) < 1.0/64.0 ) && (fabs(this->gety()) < 1.0/64.0))
+			x = Y_VECTOR ^ *this;
+		else
+			x = Z_VECTOR ^ *this;
+
+		y = *this ^ x;
+	}
+
+	 int Vector3d::setCartesianAxes(Vector3d& b, Vector3d& c) {
+#define a *this
+	// computes a RH triad of Axes (Cartesian) starting from a (normalised)
+	// if a & b are perpendicular then c = a ^ b
+	// if a & c are perpendicular then b = c ^ a
+	// if neither are perpendicular to a, then return arbitrary axes from a
+
+	// calling sequence for RH cartesian
+	//	x y z
+	//  y z x
+	//  z x y
+	if(a == NULL_VECTOR) FAILURE(L"SetAxes given a NULL Vector");
+	double epsilon = 1.0e-09;
+	bool bNull = (b == NULL_VECTOR);
+	bool cNull = (c == NULL_VECTOR);
+	bool abPerp = !bNull;
+	if(abPerp) abPerp = (fabs(a * b) < epsilon);
+
+	bool acPerp = !cNull;
+	if(acPerp) acPerp = (fabs(a * c) < epsilon);
+
+	if(abPerp) {
+		c = a ^ b;
+		return 1;
+	}
+
+	if(acPerp) {
+		b = c ^ a;
+		return 1;
+	}
+
+	arbitrary_axes(b, c);
+	b.normalise();
+	c.normalise();
+	return 2;
+}
+
+
+	 void Plane::Mirrored(Matrix* tmMirrored) {
+		 // calculates a mirror transformation that mirrors 2d about plane
+	 
+		Point3d p1 = this->Near(Point3d(0.,0.,0.));
+		if(tmMirrored->m_unit == false) tmMirrored->Unit();
+
+		double nx = this->normal.getx();
+		double ny = this->normal.gety();
+		double nz = this->normal.getz();
+	   
+		// the translation
+		tmMirrored->e[ 3] = -2. * nx * this->d;
+		tmMirrored->e[ 7] = -2. * ny * this->d;
+		tmMirrored->e[11] = -2. * nz * this->d;
+	 
+		// the rest
+		tmMirrored->e[ 0] = 1. - 2. * nx * nx;
+		tmMirrored->e[ 5] = 1. - 2. * ny * ny;
+		tmMirrored->e[10] = 1. - 2. * nz * nz;
+		tmMirrored->e[ 1] = tmMirrored->e[ 4] = -2. * nx * ny;
+		tmMirrored->e[ 2] = tmMirrored->e[ 8] = -2. * nz * nx;
+		tmMirrored->e[ 6] = tmMirrored->e[ 9] = -2. * ny * nz;
+	 
+		tmMirrored->m_unit = false;
+		tmMirrored->m_mirrored = true;
+	}
+}

src/Mod/Path/libarea/kurve/offset.cpp

@@ -0,0 +1,376 @@
+////////////////////////////////////////////////////////////////////////////////////////////////
+//                    2d geometry classes - implements 2d kurve offset for use in dll
+//
+//                    g.j.hawkesford August 2003
+//
+// This program is released under the BSD license. See the file COPYING for details.
+//
+////////////////////////////////////////////////////////////////////////////////////////////////
+#include "geometry.h"
+using namespace geoff_geometry;
+
+namespace geoff_geometry {
+	static Kurve eliminateLoops(const Kurve& k , const Kurve& originalk, double offset, int& ret);
+	static bool DoesIntersInterfere(const Point& pInt, const Kurve& k, double offset);
+
+	int Kurve::Offset(vector<Kurve*>&OffsetKurves, double offset, int direction, int method, int& ret)const {
+
+		switch(method) {
+	case NO_ELIMINATION:
+	case BASIC_OFFSET:
+		{
+			Kurve* ko = new Kurve;
+			int n = OffsetMethod1(*ko, offset, direction, method, ret);
+			OffsetKurves.push_back(ko);
+			return n;
+		}
+
+	default:
+		FAILURE(L"Requested Offsetting Method not available");
+		}
+		return 0;
+	}
+
+	int Kurve::OffsetMethod1(Kurve& kOffset, double off, int direction,  int method, int& ret)const
+	{
+		// offset kurve with simple span elimination
+		// direction 1 = left,  -1 = right
+
+		// ret  = 0		- kurve offset ok
+		//		= 1		- kurve has differential scale (not allowed)
+		//		= 2		- offset failed
+		//      = 3		- offset too large
+		if(this == &kOffset) FAILURE(L"Illegal Call - 'this' must not be kOffset");
+		double offset = (direction == GEOFF_LEFT)?off : -off;
+
+		if(fabs(offset) < geoff_geometry::TOLERANCE || m_nVertices < 2) {
+			kOffset = *this;
+			ret = 0;
+			return 1;
+		}
+
+		Span curSpan, curSpanOff;	// current & offset spans
+		Span prevSpanOff;			// previous offset span
+		Point p0, p1;				// Offset span intersections
+
+		// offset Kurve
+		kOffset = Matrix(*this);
+
+		if(m_mirrored) offset = -offset;
+		int RollDir = ( off < 0 ) ? direction : - direction;				// Roll arc direction
+
+		double scalex;
+		if(!GetScale(scalex)) {
+			ret = 1;
+			return 0;	// differential scale
+		}
+		offset /= scalex;
+
+		bool bClosed = Closed();
+		int nspans = nSpans();
+		if(bClosed) {
+			Get(nspans, curSpan, true);						// assign previus span for closed
+
+			prevSpanOff = curSpan.Offset(offset);
+			nspans++; // read first again
+		}
+
+		for(int spannumber = 1; spannumber <= nspans; spannumber++) {
+			if(spannumber > nSpans())
+				Get(1, curSpan, true);						// closed kurve - read first span again
+			else
+				Get(spannumber, curSpan, true);
+
+			if(!curSpan.NullSpan) {
+				int numint = 0;
+				curSpanOff = curSpan.Offset(offset);
+				curSpanOff.ID = 0;
+				if(!kOffset.m_started) {
+					kOffset.Start(curSpanOff.p0);
+					kOffset.AddSpanID(0);
+				}
+
+				if(spannumber > 1) {
+					// see if tangent
+					double d = curSpanOff.p0.Dist(prevSpanOff.p1);
+					if((d > geoff_geometry::TOLERANCE) && (curSpanOff.NullSpan == false && prevSpanOff.NullSpan == false)) {
+						// see if offset spans intersect
+
+						double cp = prevSpanOff.ve ^ curSpanOff.vs;
+						bool inters = (cp > 0 && direction == GEOFF_LEFT) || (cp < 0 && direction == GEOFF_RIGHT);
+
+						if(inters) {
+							double t[4];
+							numint = prevSpanOff.Intof(curSpanOff, p0, p1, t);
+						}
+
+						if(numint == 1) {
+							// intersection - modify previous endpoint
+							kOffset.Replace(kOffset.m_nVertices-1, prevSpanOff.dir, p0, prevSpanOff.pc, prevSpanOff.ID);
+						}
+						else {
+							// 0 or 2 intersections, add roll around (remove -ve loops in elimination function)
+							if(kOffset.Add(RollDir, curSpanOff.p0, curSpan.p0, false))	kOffset.AddSpanID(ROLL_AROUND);
+						}
+					}
+				}
+
+				// add span
+				if(spannumber < m_nVertices) {
+					curSpanOff.ID = spannumber;
+					kOffset.Add(curSpanOff, false);
+				}
+				else if(numint == 1)		// or replace the closed first span
+					kOffset.Replace(0, 0, p0, Point(0, 0), 0);
+
+			}
+			if(!curSpanOff.NullSpan)prevSpanOff = curSpanOff;
+		}		// end of main pre-offsetting loop
+
+
+#ifdef _DEBUG
+//testDraw->AddKurve("", &kOffset, 0, GREEN);
+//		outXML oxml(L"c:\\temp\\eliminateLoops.xml");
+//		oxml.startElement(L"eliminateLoops");
+//		oxml.Write(kOffset, L"kOffset");
+//		oxml.endElement();
+#endif
+		// eliminate loops
+		if(method == NO_ELIMINATION) {
+			ret = 0;
+			return 1;
+		}
+		kOffset = eliminateLoops(kOffset, *this, offset, ret);
+
+		if(ret == 0 && bClosed) {
+			// check for inverted offsets of closed kurves
+			if(kOffset.Closed()) {
+				double a = Area();
+				int dir = (a < 0);
+				double ao = kOffset.Area();
+				int dirOffset = ao < 0;
+
+				if(dir != dirOffset)
+					ret = 3;
+				else {
+					// check area change compatible with offset direction - catastrophic failure
+					bool bigger = (a > 0 && offset > 0) || (a < 0 && offset < 0);
+					if(bigger && fabs(ao) < fabs(a)) ret = 2;
+				}
+			}
+			else
+				ret = 2;			// started closed but now open??
+		}
+		return (ret == 0)?1 : 0;
+	}
+
+
+	static Kurve eliminateLoops(const Kurve& k , const Kurve& originalk, double offset, int& ret) {
+		// a simple loop elimination routine based on first offset ideas in Peps
+		// this needs extensive work for future
+		// start point musn't disappear & only one valid offset is determined
+		//
+		// ret = 0 for ok
+		// ret = 2 for impossible geometry
+		
+		Span sp0, sp1;
+		Point pInt, pIntOther;
+
+		Kurve ko;											// eliminated output
+		ko = Matrix(k);
+		int kinVertex = 0;
+
+		while(kinVertex <= k.nSpans()) {
+			bool clipped = false ;                                       // not in a clipped section (assumption with this simple method)
+
+			sp0.dir = k.Get(kinVertex, sp0.p0, sp0.pc);
+			sp0.ID = k.GetSpanID(kinVertex++);
+			if (kinVertex == 1)	{
+				ko.Start(sp0.p0);							// start point mustn't dissappear for this simple method
+				ko.AddSpanID(sp0.ID);
+			}
+			if (kinVertex <= k.nSpans()) {   // any more?
+				int ksaveVertex = kinVertex ;
+				sp0.dir = k.Get(kinVertex, sp0.p1, sp0.pc);	// first span
+				sp0.ID = k.GetSpanID(kinVertex++);
+
+				sp0.SetProperties(true);
+
+				int ksaveVertex1 = kinVertex;									// mark position AA		
+				if (kinVertex <= k.nSpans()) {	// get the next but one span			
+					sp1.dir = k.Get(kinVertex, sp1.p0, sp1.pc);
+					sp1.ID = k.GetSpanID(kinVertex++);
+					int ksaveVertex2 = kinVertex;								// mark position BB
+
+					int fwdCount = 0;
+					while(kinVertex <= k.nSpans()) {					
+						sp1.dir = k.Get(kinVertex, sp1.p1, sp1.pc);			// check span
+						sp1.ID = k.GetSpanID(kinVertex++);
+						sp1.SetProperties(true);
+			
+						double t[4];
+						int numint = sp0.Intof(sp1, pInt, pIntOther, t);			// find span intersections
+						if(numint && sp0.p0.Dist(pInt) < geoff_geometry::TOLERANCE ) numint=0;	// check that intersection is not at the start of the check span					
+						if(numint ) {
+
+							if(numint == 2) {
+								// choose first intercept on sp0
+								Span spd = sp0;
+								spd.p1 = pInt;
+								spd.SetProperties(true);
+								double dd = spd.length;
+
+								spd.p1 = pIntOther;
+								spd.SetProperties(true);
+								if(dd > spd.length) pInt = pIntOther;
+								numint = 1;
+
+							}
+							ksaveVertex = ksaveVertex1 ;
+
+							clipped = true ;			// in a clipped section		
+							if(DoesIntersInterfere(pInt, originalk, offset) == false) {
+								sp0.p1 = pInt;			// ok so truncate this span to the intersection
+								clipped = false;		// end of clipped section
+								break;
+							}
+							// no valid intersection found so carry on
+						}
+						sp1.p0 = sp1.p1 ;		// next
+						ksaveVertex1 = ksaveVertex2 ;							// pos AA = BB
+						ksaveVertex2 = kinVertex;								// mark 
+
+						if((kinVertex > k.nSpans() || fwdCount++ > 25) && clipped == false) break;
+					}
+				}
+
+				if(clipped) {
+					ret = 2;	// still in a clipped section - error
+
+					return ko;
+				}
+
+				ko.Add(sp0, false);
+
+				kinVertex = ksaveVertex;
+			}
+		}
+		ret = 0;
+
+		return ko; // no more spans - seems ok
+	}
+
+
+	static bool DoesIntersInterfere(const Point& pInt, const Kurve& k, double offset)  {
+		// check that intersections don't interfere with the original kurve 
+		Span sp;
+		Point dummy;
+		int kCheckVertex = 0;
+		k.Get(kCheckVertex++, sp.p0, sp.pc);
+
+		offset = fabs(offset) - geoff_geometry::TOLERANCE;
+		while(kCheckVertex <= k.nSpans()) {
+			sp.dir = k.Get(kCheckVertex++, sp.p1, sp.pc);
+			sp.SetProperties(true);
+			// check for interference 
+			if(Dist(sp, pInt, dummy) < offset) return true;
+			sp.p0 = sp.p1;
+		}
+		return false;	// intersection is ok
+	}
+}
+
+
+static struct iso {
+		 Span sp;
+		 Span off;
+	} isodata;
+static void isoRadius(Span& before, Span& blend, Span& after, double radius);
+
+int Kurve::OffsetISOMethod(Kurve& kOut, double off, int direction, bool BlendAll)const {
+		// produces a special offset Kurve - observing so-called ISO radii
+		// eg line/arc/line tangent - keep arc radius constant
+		// this method also considers arc/arc/arc etc.
+		// interior radius must be smallest of triplet for above.
+
+		// parameters:-
+		// Output	kOut		resulting kurve
+		// Input	off			offset amount
+		// Input	direction	offset direction (LEFT or RIGHT)
+		// Input	BlendAall	if false only consider ISO radius for LINE/ARC/LINE
+		//						if true consider all blended radii (ARC/ARC/ARC etc.)		
+		double offset = (direction == GEOFF_LEFT)?off : -off;
+		if(FEQZ(off) || nSpans() < 1) {
+			kOut = *this;
+			return 1;
+		}
+		double cptol = 1.0e-05;
+		std::vector<iso> spans;
+		for(int i = 0; i < nSpans(); i++) {	// store all spans and offsets
+			Get(i+1, isodata.sp, true, true);
+			isodata.off = isodata.sp.Offset(offset);
+			spans.push_back(isodata);
+		}
+
+		for(int i = 0; i < nSpans() - 1; i++)		// calculate intersections for none tangent spans
+			if(fabs(spans[i].off.ve ^ spans[i+1].off.vs) > cptol)	spans[i].off.JoinSeparateSpans(spans[i+1].off);
+
+		for(int i = 1; i < nSpans() - 1; i++) {			// deal with isoradii
+			if(spans[i].off.dir) {
+				if(BlendAll) {						// interior radius should be smaller than neighbours
+					if(spans[i-1].sp.dir)
+						if(spans[i-1].sp.radius < spans[i].sp.radius) continue;
+					if(spans[i+1].sp.dir)
+						if(spans[i+1].sp.radius < spans[i].sp.radius) continue;
+				}
+				else {
+					if((spans[i-1].off.dir || spans[i+1].off.dir)) continue;		// linear neighbours only
+				}
+
+				if((fabs(spans[i-1].sp.ve ^ spans[i].sp.vs) < cptol) && (fabs(spans[i].sp.ve ^ spans[i+1].sp.vs) < cptol)) {
+					// isoradius - calculate the new offset radius and modify neighbouring spans
+					isoRadius(spans[i-1].off, spans[i].off, spans[i+1].off, spans[i].sp.radius);
+				}
+			}
+		}
+
+		kOut.Start(spans[0].off.p0);										// start point
+		for(int i = 0; i < nSpans(); i++)
+			kOut.Add(spans[i].off.dir, spans[i].off.p1, spans[i].off.pc);	// output all spans
+		return 1;
+	}
+
+static void isoRadius(Span& before, Span& blend, Span& after, double radius) {
+	// calculate the new offset radius and modify neighbouring spans
+	int direction = ((before.ve ^ after.vs) > 0)? 1 : -1;				// offset direction
+	Span beforeOff = before.Offset(direction * radius);
+	Span afterOff = after.Offset(direction * radius);
+	int turnLeft = ((before.ve ^ after.vs) > 0)? 1 : -1;
+	if(before.dir == LINEAR) {
+		CLine b(beforeOff);
+		if(after.dir == LINEAR) {
+			CLine a(afterOff);
+			blend.pc = b.Intof(a);
+		}
+		else {
+			Circle a(afterOff);
+			b.Intof(turnLeft * after.dir, a, blend.pc);
+		}
+	}
+	else {
+		Circle b(beforeOff);
+
+		if(after.dir == LINEAR) {
+			CLine a(afterOff);
+			a.Intof(-turnLeft * before.dir, b, blend.pc);
+		}
+		else {
+			// arc arc
+			Circle a(afterOff);
+			int leftright = ((Vector2d(b.pc, blend.pc) ^ Vector2d(b.pc, a.pc)) < 0)? 1 : -1;
+			b.Intof(leftright, a, blend.pc);
+		}
+	}
+	before.p1 = blend.p0 = before.Near(blend.pc);
+	after.p0 = blend.p1 = after.Near(blend.pc);
+}

src/Mod/Path/libarea/kurve/test.py

@@ -0,0 +1,9 @@
+import area
+
+p = area.Point(0,0)
+
+m = area.Matrix([1,0,0,12, 0,1,0,0, 0,0,1,0, 0,0,0,1])
+
+p.Transform(m)
+
+print p.x,  p.y