update to 2geom rev. 1507

(bzr r6416)

25 files changed, 576 insertions, 434 deletions

diff --git a/src/2geom/convex-cover.cpp b/src/2geom/convex-cover.cpp
index 7127a7c09..7dca7f0eb 100644
--- a/src/2geom/convex-cover.cpp
+++ b/src/2geom/convex-cover.cpp
@@ -433,6 +433,66 @@ ConvexHull graham_merge(ConvexHull a, ConvexHull b) {
     }
 }*/
 
+double ConvexHull::centroid_and_area(Geom::Point& centroid) const {
+    const unsigned n = boundary.size();
+    if (n < 2)
+        return 0;
+    if(n < 3) {
+        centroid = (boundary[0] + boundary[1])/2;
+        return 0;
+    }
+    Geom::Point centroid_tmp(0,0);
+    double atmp = 0;
+    for (unsigned i = n-1, j = 0; j < n; i = j, j++) {
+        const double ai = -cross(boundary[j], boundary[i]);
+        atmp += ai;
+        centroid_tmp += (boundary[j] + boundary[i])*ai; // first moment.
+    }
+    if (atmp != 0) {
+        centroid = centroid_tmp / (3 * atmp);
+    }
+    return atmp / 2;
+}
+
+// TODO: This can be made lg(n) using golden section/fibonacci search three starting points, say 0,
+// n/2, n-1 construct a new point, say (n/2 + n)/2 throw away the furthest boundary point iterate
+// until interval is a single value
+Point const * ConvexHull::furthest(Point direction) const {
+    Point const * p = &boundary[0];
+    double d = dot(*p, direction);
+    for(unsigned i = 1; i < boundary.size(); i++) {
+        double dd = dot(boundary[i], direction);
+        if(d < dd) {
+            p = &boundary[i];
+            d = dd;
+        }
+    }
+    return p;
+}
+
+
+// returns (a, (b,c)), three points which define the narrowest diameter of the hull as the pair of
+// lines going through b,c, and through a, parallel to b,c TODO: This can be made linear time by
+// moving point tc incrementally from the previous value (it can only move in one direction).  It
+// is currently n*O(furthest)
+double ConvexHull::narrowest_diameter(Point &a, Point &b, Point &c) {
+    Point tb = boundary.back();
+    double d = INFINITY;
+    for(unsigned i = 0; i < boundary.size(); i++) {
+        Point tc = boundary[i];
+        Point n = -rot90(tb-tc);
+        Point ta = *furthest(n);
+        double td = dot(n, ta-tb)/dot(n,n);
+        if(td < d) {
+            a = ta;
+            b = tb;
+            c = tc;
+            d = td;
+        }
+        tb = tc;
+    }
+    return d;
+}
 
 };
 
diff --git a/src/2geom/convex-cover.h b/src/2geom/convex-cover.h
index 1fe86e40d..71a47bb98 100644
--- a/src/2geom/convex-cover.h
+++ b/src/2geom/convex-cover.h
@@ -106,13 +106,18 @@ public:
     bool is_degenerate() const;
     
     // area of the convex hull
-    double area() const;
+    double centroid_and_area(Geom::Point& centroid) const;
+    double area() const {
+        Point tmp;
+        return centroid_and_area(tmp);
+    }
     
     // furthest point in a direction (lg time) 
     Point const * furthest(Point direction) const;
 
     bool is_left(Point p, int n);
     int find_left(Point p);
+    double narrowest_diameter(Point &a, Point &b, Point &c);
 };
 
 // do two convex hulls intersect?
diff --git a/src/2geom/coord.h b/src/2geom/coord.h
index b2c2a4766..99e249e29 100644
--- a/src/2geom/coord.h
+++ b/src/2geom/coord.h
@@ -32,6 +32,7 @@
 #define SEEN_Geom_COORD_H
 
 #include <cmath>
+#include <limits>
 
 namespace Geom {
 
@@ -45,7 +46,8 @@ namespace Geom {
 typedef double Coord;
 
 const Coord EPSILON = 1e-5; //1e-18;
-const Coord COORD_HUGE = 1e18;
+
+inline Coord infinity() {  return std::numeric_limits<Coord>::infinity();  }
 
 //IMPL: NearConcept
 inline bool are_near(Coord a, Coord b, double eps=EPSILON) { return fabs(a-b) <= eps; }
diff --git a/src/2geom/hvlinesegment.h b/src/2geom/hvlinesegment.h
index a34c5a962..66de31d70 100644
--- a/src/2geom/hvlinesegment.h
+++ b/src/2geom/hvlinesegment.h
@@ -41,208 +41,212 @@ namespace Geom
 class HLineSegment : public Curve
 {
   public:
-	HLineSegment()
-	{}
-	
-	HLineSegment(Coord _x0, Coord _x1, Coord _y)
-		: m_line_seg(Point(_x0, _y), Point(_x1, _y))
-	{
-	}
-	
-	HLineSegment(Point const& _p, double _length)
-		: m_line_seg(_p, Point(_p[X] + _length, _p[Y]))
-	{
-	}
-	
-	HLineSegment(Point const& _p0, Point const& _p1)
-		: m_line_seg(_p0, _p1)
-	{
-		if ( _p0[Y] != _p1[Y] )
-		{
-			THROW_RANGEERROR("HLineSegment::HLineSegment passed points should "
-					         "have the same Y value");
-		}
-	}
-	
-	Curve* duplicate() const
-	{
-		return new HLineSegment(*this);
-	}
-	
-	bool isDegenerate() const
-	{
-		return m_line_seg.isDegenerate();
-	}
-	
-	Point initialPoint() const
-	{
-		return m_line_seg.initialPoint();
-	}
-	
-	Point finalPoint() const 
-	{
-		return m_line_seg.finalPoint();
-	}
-	
-	Coord getY()
-	{
-		return initialPoint()[Y];
-	}
-	
-	void setInitial(Point _p) 
-	{ 
-		m_line_seg.setInitial( Point(_p[X], initialPoint()[Y]) ); 
-	}
-	
-	void setFinal(Point _p) 
-	{ 
-		m_line_seg.setFinal( Point(_p[X], finalPoint()[Y]) ); 
-	}
-	
-	void setX0(Coord _x)
-	{
-		m_line_seg.setInitial( Point(_x, initialPoint()[Y]) );
-	}
-	
-	void setX1(Coord _x)
-	{
-		m_line_seg.setFinal( Point(_x, finalPoint()[Y]) );
-	}
-	
-	void setY(Coord _y)
-	{
-		m_line_seg.setInitial( Point(initialPoint()[X], _y) );
-		m_line_seg.setFinal( Point(finalPoint()[X], _y) );
-	}
-
-	Rect boundsFast() const
-	{
-		return boundsExact();
-	}
-	
-	Rect boundsExact() const
-	{
-		return Rect( initialPoint(), finalPoint() );
-	}
-	
-	Rect boundsLocal(Interval i, unsigned deg) const
-	{
-		return m_line_seg.boundsLocal(i, deg);
-	}
-	
-	int winding(Point p) const
-	{
-		return m_line_seg.winding(p);
-	}
-	
-	std::vector<double>
-	roots(double v, Dim2 d) const
-	{
-		if (d < 0 || d > 1)
-		{
-			THROW_RANGEERROR("dimension argument out of range");
-		}
-		std::vector<double> result;
-		if (d == X)
-		{
-			if ( v >= initialPoint()[X] && v <= finalPoint()[X] )
-			{
-				double t = 0;
-				if (!isDegenerate())
-					t = (v - initialPoint()[X]) / (finalPoint()[X] - initialPoint()[X]);
-				result.push_back(t);
-			}
-		}
-		else
-		{
-			if (v == initialPoint()[Y])
-			{
-				if (!isDegenerate())
-					THROW_INFINITESOLUTIONS(0);
-				result.push_back(0);
-			}
-		}
-		return result;
-	}
-	
-	double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
-	{
-		if ( from > to ) std::swap(from, to);
-		double xfrom = pointAt(from)[X];
-		double xto = pointAt(to)[X];
-		if ( xfrom > xto )
-		{
-		    std::swap(xfrom, xto);
-		    std::swap(from, to);
-		}
-		if ( p[X] > xfrom && p[X] < xto )
-		{
-			return (p[X] - initialPoint()[X]) / (finalPoint()[X] - initialPoint()[X]);
-		}
-		else if ( p[X] <= xfrom )
-			return from;
-		else
-			return to;
-	}
-	
-	std::pair<HLineSegment, HLineSegment> subdivide(Coord t) const
-	{
-		std::pair<HLineSegment, HLineSegment> result;
-		Point p = pointAt(t);
-		result.first.setInitial(initialPoint());
-		result.first.setFinal(p);
-		result.second.setInitial(p);
-		result.second.setFinal(finalPoint());
-		return result;
-	}
-	
-	Curve* portion(double f, double t) const
-	{
-		Point ip = pointAt(f);
-		Point ep = pointAt(t);
-		return new HLineSegment(ip[X], ep[X], ip[Y]);
-	}
-	
-	Curve* reverse() const
-	{
-		return 
-		new HLineSegment(finalPoint()[X], initialPoint()[X], initialPoint()[Y]);
-	}
-	
-	Curve* transformed(Matrix const & m) const
-	{
-		Point ip = initialPoint() * m;
-		Point ep = finalPoint() * m;
-		return new LineSegment(ip, ep);
-	}
-	
-	Curve* derivative() const
-	{
-		double x = finalPoint()[X] - initialPoint()[X];
-		return new HLineSegment(x, x, 0);
-	}
-	
-	Point pointAt(double t) const
-	{
-		if ( t < 0 || t > 1 )
-			THROW_RANGEERROR("domain parameter out of range");
-		double x = initialPoint()[X] + t * (finalPoint()[X] - initialPoint()[X]);
-		return Point(x, initialPoint()[Y]);
-	}
-	
-	double valueAt(double t, Dim2 d) const
-	{
-		if (d < 0 || d > 1)
-		{
-			THROW_RANGEERROR("dimension argument out of range");
-		}
-		if ( t < 0 || t > 1 )
-			THROW_RANGEERROR("domain parameter out of range");
-		
-		if (d == Y) return initialPoint()[Y];
-		
-		return initialPoint()[X] + t * (finalPoint()[X] - initialPoint()[X]);
-	}
+    HLineSegment()
+    {}
+
+    HLineSegment(Coord _x0, Coord _x1, Coord _y)
+        : m_line_seg(Point(_x0, _y), Point(_x1, _y))
+    {
+    }
+
+    HLineSegment(Point const& _p, double _length)
+        : m_line_seg(_p, Point(_p[X] + _length, _p[Y]))
+    {
+    }
+
+    HLineSegment(Point const& _p0, Point const& _p1)
+        : m_line_seg(_p0, _p1)
+    {
+        if ( _p0[Y] != _p1[Y] )
+        {
+            THROW_RANGEERROR("HLineSegment::HLineSegment passed points should "
+                             "have the same Y value");
+        }
+    }
+
+    Curve* duplicate() const
+    {
+        return new HLineSegment(*this);
+    }
+
+    bool isDegenerate() const
+    {
+        return m_line_seg.isDegenerate();
+    }
+
+    Point initialPoint() const
+    {
+        return m_line_seg.initialPoint();
+    }
+
+    Point finalPoint() const
+    {
+        return m_line_seg.finalPoint();
+    }
+
+    Coord getY()
+    {
+        return initialPoint()[Y];
+    }
+
+    void setInitial(Point _p)
+    {
+        m_line_seg.setInitial( Point(_p[X], initialPoint()[Y]) );
+    }
+
+    void setFinal(Point _p)
+    {
+        m_line_seg.setFinal( Point(_p[X], finalPoint()[Y]) );
+    }
+
+    void setX0(Coord _x)
+    {
+        m_line_seg.setInitial( Point(_x, initialPoint()[Y]) );
+    }
+
+    void setX1(Coord _x)
+    {
+        m_line_seg.setFinal( Point(_x, finalPoint()[Y]) );
+    }
+
+    void setY(Coord _y)
+    {
+        m_line_seg.setInitial( Point(initialPoint()[X], _y) );
+        m_line_seg.setFinal( Point(finalPoint()[X], _y) );
+    }
+
+    Rect boundsFast() const
+    {
+        return boundsExact();
+    }
+
+    Rect boundsExact() const
+    {
+        return Rect( initialPoint(), finalPoint() );
+    }
+
+    Rect boundsLocal(Interval i, unsigned deg) const
+    {
+        return m_line_seg.boundsLocal(i, deg);
+    }
+
+    int winding(Point p) const
+    {
+        return m_line_seg.winding(p);
+    }
+
+    std::vector<double>
+    roots(double v, Dim2 d) const
+    {
+        if (d < 0 || d > 1)
+        {
+            THROW_RANGEERROR("dimension argument out of range");
+        }
+        std::vector<double> result;
+        if (d == X)
+        {
+            if ( v >= initialPoint()[X] && v <= finalPoint()[X] )
+            {
+                double t = 0;
+                if (!isDegenerate())
+                    t = (v - initialPoint()[X]) / (finalPoint()[X] - initialPoint()[X]);
+                result.push_back(t);
+            }
+        }
+        else
+        {
+            if (v == initialPoint()[Y])
+            {
+                if (!isDegenerate())
+                    THROW_INFINITESOLUTIONS(0);
+                result.push_back(0);
+            }
+        }
+        return result;
+    }
+
+    double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
+    {
+        if ( from > to ) std::swap(from, to);
+        double xfrom = pointAt(from)[X];
+        double xto = pointAt(to)[X];
+        if ( xfrom > xto )
+        {
+            std::swap(xfrom, xto);
+            std::swap(from, to);
+        }
+        if ( p[X] > xfrom && p[X] < xto )
+        {
+            return (p[X] - initialPoint()[X]) / (finalPoint()[X] - initialPoint()[X]);
+        }
+        else if ( p[X] <= xfrom )
+            return from;
+        else
+            return to;
+    }
+
+    std::pair<HLineSegment, HLineSegment> subdivide(Coord t) const
+    {
+        std::pair<HLineSegment, HLineSegment> result;
+        Point p = pointAt(t);
+        result.first.setInitial(initialPoint());
+        result.first.setFinal(p);
+        result.second.setInitial(p);
+        result.second.setFinal(finalPoint());
+        return result;
+    }
+
+    Curve* portion(double f, double t) const
+    {
+        Point ip = pointAt(f);
+        Point ep = pointAt(t);
+        return new HLineSegment(ip[X], ep[X], ip[Y]);
+    }
+
+    Curve* reverse() const
+    {
+        return
+        new HLineSegment(finalPoint()[X], initialPoint()[X], initialPoint()[Y]);
+    }
+
+    Curve* transformed(Matrix const & m) const
+    {
+        Point ip = initialPoint() * m;
+        Point ep = finalPoint() * m;
+        if (m.onlyScaleAndTranslation()) {
+            return new HLineSegment(ip[X], ep[X], ip[Y]);
+        } else {
+            return new LineSegment(ip, ep);
+        }
+    }
+
+    Curve* derivative() const
+    {
+        double x = finalPoint()[X] - initialPoint()[X];
+        return new HLineSegment(x, x, 0);
+    }
+
+    Point pointAt(double t) const
+    {
+        if ( t < 0 || t > 1 )
+            THROW_RANGEERROR("domain parameter out of range");
+        double x = initialPoint()[X] + t * (finalPoint()[X] - initialPoint()[X]);
+        return Point(x, initialPoint()[Y]);
+    }
+
+    double valueAt(double t, Dim2 d) const
+    {
+        if (d < 0 || d > 1)
+        {
+            THROW_RANGEERROR("dimension argument out of range");
+        }
+        if ( t < 0 || t > 1 )
+            THROW_RANGEERROR("domain parameter out of range");
+
+        if (d == Y) return initialPoint()[Y];
+
+        return initialPoint()[X] + t * (finalPoint()[X] - initialPoint()[X]);
+    }
 
     std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const
     {
@@ -262,222 +266,226 @@ class HLineSegment : public Curve
         return result;
     }
 
-	D2<SBasis> toSBasis() const
-	{
-		return m_line_seg.toSBasis();
-	}
-	
+    D2<SBasis> toSBasis() const
+    {
+        return m_line_seg.toSBasis();
+    }
+
   private:
-	LineSegment m_line_seg;
-	
+    LineSegment m_line_seg;
+
 };  // end class HLineSegment
 
 
 class VLineSegment : public Curve
 {
   public:
-	VLineSegment()
-	{}
-	
-	VLineSegment(Coord _x, Coord _y0, Coord _y1)
-		: m_line_seg(Point(_x, _y0), Point(_x, _y1))
-	{
-	}
-	
-	VLineSegment(Point const& _p, double _length)
-		: m_line_seg(_p, Point(_p[X], _p[Y] + _length))
-	{
-	}
-	
-	VLineSegment(Point const& _p0, Point const& _p1)
-		: m_line_seg(_p0, _p1)
-	{
-		if ( _p0[X] != _p1[X] )
-		{
-			THROW_RANGEERROR("VLineSegment::VLineSegment passed points should "
-					         "have the same X value");
-		}
-	}
-	
-	Curve* duplicate() const
-	{
-		return new VLineSegment(*this);
-	}
-	
-	bool isDegenerate() const
-	{
-		return m_line_seg.isDegenerate();
-	}
-	
-	Point initialPoint() const
-	{
-		return m_line_seg.initialPoint();
-	}
-	
-	Point finalPoint() const 
-	{
-		return m_line_seg.finalPoint();
-	}
-	
-	Coord getX()
-	{
-		return initialPoint()[X];
-	}
-	
-	void setInitial(Point _p) 
-	{ 
-		m_line_seg.setInitial( Point(initialPoint()[X], _p[Y]) ); 
-	}
-	
-	void setFinal(Point _p) 
-	{ 
-		m_line_seg.setFinal( Point(finalPoint()[X], _p[Y]) ); 
-	}
-	
-	void setY0(Coord _y)
-	{
-		m_line_seg.setInitial( Point(initialPoint()[X], _y) );
-	}
-	
-	void setY1(Coord _y)
-	{
-		m_line_seg.setFinal( Point(finalPoint()[Y], _y) );
-	}
-	
-	void setX(Coord _x)
-	{
-		m_line_seg.setInitial( Point(_x, initialPoint()[Y]) );
-		m_line_seg.setFinal( Point(_x, finalPoint()[Y]) );
-	}
-
-	Rect boundsFast() const
-	{
-		return boundsExact();
-	}
-	
-	Rect boundsExact() const
-	{
-		return Rect( initialPoint(), finalPoint() );
-	}
-	
-	Rect boundsLocal(Interval i, unsigned deg) const
-	{
-		return m_line_seg.boundsLocal(i, deg);
-	}
-	
-	int winding(Point p) const
-	{
-		return m_line_seg.winding(p);
-	}
-	
-	std::vector<double>
-	roots(double v, Dim2 d) const
-	{
-		if (d < 0 || d > 1)
-		{
-			THROW_RANGEERROR("dimension argument out of range");
-		}
-		std::vector<double> result;
-		if (d == Y)
-		{
-			if ( v >= initialPoint()[Y] && v <= finalPoint()[Y] )
-			{
-				double t = 0;
-				if (!isDegenerate())
-					t = (v - initialPoint()[Y]) / (finalPoint()[Y] - initialPoint()[Y]);
-				result.push_back(t);
-			}
-		}
-		else
-		{
-			if (v == initialPoint()[X])
-			{
-				if (!isDegenerate())
-					THROW_INFINITESOLUTIONS(0);
-				result.push_back(0);
-			}
-		}
-		return result;
-	}
-	
-	double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
-	{
-		if ( from > to ) std::swap(from, to);
-		double yfrom = pointAt(from)[Y];
-		double yto = pointAt(to)[Y];
-		if (yfrom > yto)
-		{
-		    std::swap(yfrom, yto);
-		    std::swap(from, to);
-		}
-		if ( p[Y] > yfrom && p[Y] < yto )
-		{
-			return (p[Y] - initialPoint()[Y]) / (finalPoint()[Y] - initialPoint()[Y]);
-		}
-		else if ( p[Y] <= yfrom )
-			return from;
-		else
-			return to;
-	}
-	
-	std::pair<VLineSegment, VLineSegment> subdivide(Coord t) const
-	{
-		std::pair<VLineSegment, VLineSegment> result;
-		Point p = pointAt(t);
-		result.first.setInitial(initialPoint());
-		result.first.setFinal(p);
-		result.second.setInitial(p);
-		result.second.setFinal(finalPoint());
-		return result;
-	}
-	
-	Curve* portion(double f, double t) const
-	{
-		Point ip = pointAt(f);
-		Point ep = pointAt(t);
-		return new VLineSegment(ip[X], ip[Y], ep[Y]);
-	}
-	
-	Curve* reverse() const
-	{
-		return 
-		new VLineSegment(initialPoint()[X], finalPoint()[Y], initialPoint()[Y]);
-	}
-	
-	Curve* transformed(Matrix const & m) const
-	{
-		Point ip = initialPoint() * m;
-		Point ep = finalPoint() * m;
-		return new LineSegment(ip, ep);
-	}
-	
-	Curve* derivative() const
-	{
-		double y = finalPoint()[Y] - initialPoint()[Y];
-		return new VLineSegment(0, y, y);
-	}
-	
-	Point pointAt(double t) const
-	{
-		if ( t < 0 || t > 1 )
-			THROW_RANGEERROR("domain parameter out of range");
-		double y = initialPoint()[Y] + t * (finalPoint()[Y] - initialPoint()[Y]);
-		return Point(initialPoint()[X], y);
-	}
-	
-	double valueAt(double t, Dim2 d) const
-	{
-		if (d < 0 || d > 1)
-		{
-			THROW_RANGEERROR("dimension argument out of range");
-		}
-		if ( t < 0 || t > 1 )
-			THROW_RANGEERROR("domain parameter out of range");
-		
-		if (d == X) return initialPoint()[X];
-		
-		return initialPoint()[Y] + t * (finalPoint()[Y] - initialPoint()[Y]);
-	}
+    VLineSegment()
+    {}
+
+    VLineSegment(Coord _x, Coord _y0, Coord _y1)
+        : m_line_seg(Point(_x, _y0), Point(_x, _y1))
+    {
+    }
+
+    VLineSegment(Point const& _p, double _length)
+        : m_line_seg(_p, Point(_p[X], _p[Y] + _length))
+    {
+    }
+
+    VLineSegment(Point const& _p0, Point const& _p1)
+        : m_line_seg(_p0, _p1)
+    {
+        if ( _p0[X] != _p1[X] )
+        {
+            THROW_RANGEERROR("VLineSegment::VLineSegment passed points should "
+                             "have the same X value");
+        }
+    }
+
+    Curve* duplicate() const
+    {
+        return new VLineSegment(*this);
+    }
+
+    bool isDegenerate() const
+    {
+        return m_line_seg.isDegenerate();
+    }
+
+    Point initialPoint() const
+    {
+        return m_line_seg.initialPoint();
+    }
+
+    Point finalPoint() const
+    {
+        return m_line_seg.finalPoint();
+    }
+
+    Coord getX()
+    {
+        return initialPoint()[X];
+    }
+
+    void setInitial(Point _p)
+    {
+        m_line_seg.setInitial( Point(initialPoint()[X], _p[Y]) );
+    }
+
+    void setFinal(Point _p)
+    {
+        m_line_seg.setFinal( Point(finalPoint()[X], _p[Y]) );
+    }
+
+    void setY0(Coord _y)
+    {
+        m_line_seg.setInitial( Point(initialPoint()[X], _y) );
+    }
+
+    void setY1(Coord _y)
+    {
+        m_line_seg.setFinal( Point(finalPoint()[Y], _y) );
+    }
+
+    void setX(Coord _x)
+    {
+        m_line_seg.setInitial( Point(_x, initialPoint()[Y]) );
+        m_line_seg.setFinal( Point(_x, finalPoint()[Y]) );
+    }
+
+    Rect boundsFast() const
+    {
+        return boundsExact();
+    }
+
+    Rect boundsExact() const
+    {
+        return Rect( initialPoint(), finalPoint() );
+    }
+
+    Rect boundsLocal(Interval i, unsigned deg) const
+    {
+        return m_line_seg.boundsLocal(i, deg);
+    }
+
+    int winding(Point p) const
+    {
+        return m_line_seg.winding(p);
+    }
+
+    std::vector<double>
+    roots(double v, Dim2 d) const
+    {
+        if (d < 0 || d > 1)
+        {
+            THROW_RANGEERROR("dimension argument out of range");
+        }
+        std::vector<double> result;
+        if (d == Y)
+        {
+            if ( v >= initialPoint()[Y] && v <= finalPoint()[Y] )
+            {
+                double t = 0;
+                if (!isDegenerate())
+                    t = (v - initialPoint()[Y]) / (finalPoint()[Y] - initialPoint()[Y]);
+                result.push_back(t);
+            }
+        }
+        else
+        {
+            if (v == initialPoint()[X])
+            {
+                if (!isDegenerate())
+                    THROW_INFINITESOLUTIONS(0);
+                result.push_back(0);
+            }
+        }
+        return result;
+    }
+
+    double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
+    {
+        if ( from > to ) std::swap(from, to);
+        double yfrom = pointAt(from)[Y];
+        double yto = pointAt(to)[Y];
+        if (yfrom > yto)
+        {
+            std::swap(yfrom, yto);
+            std::swap(from, to);
+        }
+        if ( p[Y] > yfrom && p[Y] < yto )
+        {
+            return (p[Y] - initialPoint()[Y]) / (finalPoint()[Y] - initialPoint()[Y]);
+        }
+        else if ( p[Y] <= yfrom )
+            return from;
+        else
+            return to;
+    }
+
+    std::pair<VLineSegment, VLineSegment> subdivide(Coord t) const
+    {
+        std::pair<VLineSegment, VLineSegment> result;
+        Point p = pointAt(t);
+        result.first.setInitial(initialPoint());
+        result.first.setFinal(p);
+        result.second.setInitial(p);
+        result.second.setFinal(finalPoint());
+        return result;
+    }
+
+    Curve* portion(double f, double t) const
+    {
+        Point ip = pointAt(f);
+        Point ep = pointAt(t);
+        return new VLineSegment(ip[X], ip[Y], ep[Y]);
+    }
+
+    Curve* reverse() const
+    {
+        return
+        new VLineSegment(initialPoint()[X], finalPoint()[Y], initialPoint()[Y]);
+    }
+
+    Curve* transformed(Matrix const & m) const
+    {
+        Point ip = initialPoint() * m;
+        Point ep = finalPoint() * m;
+        if (m.onlyScaleAndTranslation()) {
+            return new VLineSegment(ip[X], ip[Y], ep[Y]);
+        } else {
+            return new LineSegment(ip, ep);
+        }
+    }
+
+    Curve* derivative() const
+    {
+        double y = finalPoint()[Y] - initialPoint()[Y];
+        return new VLineSegment(0, y, y);
+    }
+
+    Point pointAt(double t) const
+    {
+        if ( t < 0 || t > 1 )
+            THROW_RANGEERROR("domain parameter out of range");
+        double y = initialPoint()[Y] + t * (finalPoint()[Y] - initialPoint()[Y]);
+        return Point(initialPoint()[X], y);
+    }
+
+    double valueAt(double t, Dim2 d) const
+    {
+        if (d < 0 || d > 1)
+        {
+            THROW_RANGEERROR("dimension argument out of range");
+        }
+        if ( t < 0 || t > 1 )
+            THROW_RANGEERROR("domain parameter out of range");
+
+        if (d == X) return initialPoint()[X];
+
+        return initialPoint()[Y] + t * (finalPoint()[Y] - initialPoint()[Y]);
+    }
 
     std::vector<Point> pointAndDerivatives(Coord t, unsigned n) const
     {
@@ -497,14 +505,14 @@ class VLineSegment : public Curve
         return result;
     }
 
-	D2<SBasis> toSBasis() const
-	{
-		return m_line_seg.toSBasis();
-	}
-	
+    D2<SBasis> toSBasis() const
+    {
+        return m_line_seg.toSBasis();
+    }
+
   private:
-	LineSegment m_line_seg;
-	
+    LineSegment m_line_seg;
+
 }; // end class VLineSegment
 
 
@@ -512,7 +520,7 @@ class VLineSegment : public Curve
 }  // end namespace Geom
 
 
-#endif // _2GEOM_HVLINESEGMENT_H_ 
+#endif // _2GEOM_HVLINESEGMENT_H_
 
 
 /*
diff --git a/src/2geom/interval.h b/src/2geom/interval.h
index 9887c044f..74e03aec9 100644
--- a/src/2geom/interval.h
+++ b/src/2geom/interval.h
@@ -44,16 +44,16 @@
 namespace Geom {
 
 /* Although an Interval where _b[0] > _b[1] is considered empty, for proper functioning of other methods,
- * a proper empty Interval is [+COORD_HUGE, -COORD_HUGE]. Then, expandTo(p) will set the interval to [p,p].
+ * a proper empty Interval is [+infinity, -infinity]. Then, expandTo(p) will set the interval to [p,p].
  */
 class Interval {
 private:
     Coord _b[2];
 
 public:
-    // The default constructor creates an empty interval, that ranges from +COORD_HUGE to -COORD_HUGE.
+    // The default constructor creates an empty interval, that ranges from +infinity to -infinity.
     // Doing an expandTo(p) on this empty interval will correctly set the whole interval to [p,p].
-    explicit Interval() { _b[0] = +COORD_HUGE;  _b[1] = -COORD_HUGE; }
+    explicit Interval() { _b[0] = +infinity();  _b[1] = -infinity(); }
     explicit Interval(Coord u) { _b[0] = _b[1] = u; }
     /* When creating an Interval using the constructor specifying the exact range, the created interval
      * will be [u,v] when u<=v ; and will be [v,u] when v < u !!!
diff --git a/src/2geom/matrix.cpp b/src/2geom/matrix.cpp
index d25cc8f7e..6286de86e 100644
--- a/src/2geom/matrix.cpp
+++ b/src/2geom/matrix.cpp
@@ -128,7 +128,7 @@ bool Matrix::isTranslation(Coord const eps) const {
  \return A bool representing yes/no.
  */
 bool Matrix::isScale(Coord const eps) const {
-    return !are_near(_c[0], 1.0, eps) || !are_near(_c[3], 1.0, eps) &&  //NOTE: these are the diags, and the next line opposite diags
+    return (!are_near(_c[0], 1.0, eps) || !are_near(_c[3], 1.0, eps)) &&  //NOTE: these are the diags, and the next line opposite diags
            are_near(_c[1], 0.0, eps) && are_near(_c[2], 0.0, eps) && 
            are_near(_c[4], 0.0, eps) && are_near(_c[5], 0.0, eps);
 }
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index d39b3ca92..5f321d517 100644
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
@@ -194,7 +194,7 @@ Path::allNearestPoints(Point const& _point, double from, double to) const
 	return all_nearest;	
 }
 
-double Path::nearestPoint(Point const& _point, double from, double to) const
+double Path::nearestPoint(Point const &_point, double from, double to, double *distance_squared) const
 {
 	if ( from > to ) std::swap(from, to);
 	const Path& _path = *this;
@@ -220,10 +220,11 @@ double Path::nearestPoint(Point const& _point, double from, double to) const
 	}
 	if ( si == ei )
 	{
-		double nearest =
-			_path[si].nearestPoint(_point, st, et);
+		double nearest = _path[si].nearestPoint(_point, st, et);
+        *distance_squared = distanceSq(_point, _path[si].pointAt(nearest));
 		return si + nearest;
 	}
+
 	double t;
 	double nearest = _path[si].nearestPoint(_point, st);
 	unsigned int ni = si;
@@ -254,8 +255,13 @@ double Path::nearestPoint(Point const& _point, double from, double to) const
 		{
 			nearest = t;
 			ni = ei;
+			mindistsq = dsq;
 		}
 	}
+
+    if (distance_squared) {
+        *distance_squared = mindistsq;
+    }
 	return ni + nearest;
 }
 
diff --git a/src/2geom/path.h b/src/2geom/path.h
index f391c9c24..ab91038a6 100644
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
@@ -39,7 +39,8 @@
 
 
 #include <boost/shared_ptr.hpp>
-#include <2geom/curves.h>
+#include <2geom/curve.h>
+#include <2geom/bezier-curve.h>
 
 #include <iterator>
 #include <algorithm>
@@ -357,13 +358,13 @@ public:
   }
   
   
-  double nearestPoint(Point const& _point, double from, double to) const;
+  double nearestPoint(Point const& _point, double from, double to, double *distance_squared = NULL) const;
   
-  double nearestPoint(Point const& _point) const
+  double nearestPoint(Point const& _point, double *distance_squared = NULL) const
   {
 	  unsigned int sz = size();
 	  if ( closed() ) ++sz;
-	  return nearestPoint(_point, 0, sz);
+	  return nearestPoint(_point, 0, sz, distance_squared);
   }
    
   void appendPortionTo(Path &p, double f, double t) const;
diff --git a/src/2geom/pathvector.cpp b/src/2geom/pathvector.cpp
index 1e0a20003..6d4e24c9f 100644
--- a/src/2geom/pathvector.cpp
+++ b/src/2geom/pathvector.cpp
@@ -88,6 +88,32 @@ Rect bounds_exact( PathVector const& pv )
     return bound;
 }
 
+/* Note: undefined for empty pathvectors or pathvectors with empty paths.
+ * */
+PathVectorPosition nearestPoint(PathVector const & path_in, Point const& _point, double *distance_squared)
+{
+    PathVectorPosition retval;
+
+    double mindsq = infinity();
+    unsigned int i = 0;
+    for (Geom::PathVector::const_iterator pit = path_in.begin(); pit != path_in.end(); ++pit) {
+        double dsq;
+        double t = pit->nearestPoint(_point, &dsq);
+        if (dsq < mindsq) {
+            mindsq = dsq;
+            retval.path_nr = i;
+            retval.t       = t;
+        }
+
+        ++i;
+    }
+
+    if (distance_squared) {
+        *distance_squared = mindsq;
+    }
+    return retval;
+}
+
 } // namespace Geom
 
 #endif // SEEN_GEOM_PATHVECTOR_CPP
diff --git a/src/2geom/pathvector.h b/src/2geom/pathvector.h
index 14d1efdfa..23843c52d 100644
--- a/src/2geom/pathvector.h
+++ b/src/2geom/pathvector.h
@@ -92,7 +92,16 @@ PathVector reverse_paths_and_order (PathVector const & path_in);
 
 Rect bounds_fast( PathVector const & pv );
 Rect bounds_exact( PathVector const & pv );
-}
+
+struct PathVectorPosition {
+    // pathvector[path_nr].pointAt(t) is the position
+    unsigned int path_nr;
+    double       t;
+};
+PathVectorPosition nearestPoint(PathVector const & path_in, Point const& _point, double *distance_squared = NULL);
+
+
+} // end namespace Geom
 
 #endif // SEEN_GEOM_PATHVECTOR_H
 
diff --git a/src/2geom/svg-path.h b/src/2geom/svg-path.h
index 666a249af..6da5afb7e 100644
--- a/src/2geom/svg-path.h
+++ b/src/2geom/svg-path.h
@@ -32,6 +32,7 @@
 #define SEEN_SVG_PATH_H
 
 #include <2geom/path.h>
+#include <2geom/curves.h>
 #include <iterator>
 
 namespace Geom {
diff --git a/src/display/inkscape-cairo.cpp b/src/display/inkscape-cairo.cpp
index f2892491a..d077abbd2 100644
--- a/src/display/inkscape-cairo.cpp
+++ b/src/display/inkscape-cairo.cpp
@@ -22,6 +22,8 @@
 #include "nr-arena.h"
 #include "sp-canvas.h"
 #include <2geom/pathvector.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/hvlinesegment.h>
 #include <2geom/matrix.h>
 #include <2geom/point.h>
 #include <2geom/path.h>
diff --git a/src/extension/internal/latex-pstricks.cpp b/src/extension/internal/latex-pstricks.cpp
index 343684c2e..9909d7436 100644
--- a/src/extension/internal/latex-pstricks.cpp
+++ b/src/extension/internal/latex-pstricks.cpp
@@ -29,6 +29,8 @@
 #include <libnr/nr-matrix-fns.h>
 #include <2geom/pathvector.h>
 #include <2geom/sbasis-to-bezier.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/hvlinesegment.h>
 
 #include "sp-item.h"
 
diff --git a/src/extension/internal/odf.cpp b/src/extension/internal/odf.cpp
index e4570c59a..2a71eabca 100644
--- a/src/extension/internal/odf.cpp
+++ b/src/extension/internal/odf.cpp
@@ -52,6 +52,8 @@
 #include <style.h>
 #include "display/curve.h"
 #include <2geom/pathvector.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/hvlinesegment.h>
 #include <2geom/transforms.h>
 #include <helper/geom.h>
 #include "extension/system.h"
diff --git a/src/extension/internal/pov-out.cpp b/src/extension/internal/pov-out.cpp
index 2ccbfd0f2..d1a8e833c 100644
--- a/src/extension/internal/pov-out.cpp
+++ b/src/extension/internal/pov-out.cpp
@@ -30,6 +30,8 @@
 #include <extension/system.h>
 #include <2geom/pathvector.h>
 #include <2geom/rect.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/hvlinesegment.h>
 #include "helper/geom.h"
 #include <io/sys.h>
 
diff --git a/src/extension/internal/ps.cpp b/src/extension/internal/ps.cpp
index aec3cf951..7ee28f87a 100644
--- a/src/extension/internal/ps.cpp
+++ b/src/extension/internal/ps.cpp
@@ -78,6 +78,8 @@
 #include <cmath>
 
 #include <2geom/sbasis-to-bezier.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/hvlinesegment.h>
 
 /*
 using std::atof;
diff --git a/src/helper/geom.cpp b/src/helper/geom.cpp
index f94881e0c..dcb7b3eb4 100644
--- a/src/helper/geom.cpp
+++ b/src/helper/geom.cpp
@@ -15,6 +15,8 @@
 #include <typeinfo>
 #include <2geom/pathvector.h>
 #include <2geom/path.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/hvlinesegment.h>
 #include <2geom/transforms.h>
 #include <2geom/rect.h>
 #include <2geom/coord.h>
diff --git a/src/livarot/PathCutting.cpp b/src/livarot/PathCutting.cpp
index 0ad6171c1..e84f85e0a 100644
--- a/src/livarot/PathCutting.cpp
+++ b/src/livarot/PathCutting.cpp
@@ -26,6 +26,7 @@
 #include <2geom/point.h>
 #include <2geom/matrix.h>
 #include <2geom/sbasis-to-bezier.h>
+#include <2geom/curves.h>
 #include "../display/canvas-bpath.h"
 
 void  Path::DashPolyline(float head,float tail,float body,int nbD,float *dashs,bool stPlain,float stOffset)
diff --git a/src/live_effects/lpe-copy_rotate.cpp b/src/live_effects/lpe-copy_rotate.cpp
index 1eedfa3a7..948da335a 100644
--- a/src/live_effects/lpe-copy_rotate.cpp
+++ b/src/live_effects/lpe-copy_rotate.cpp
@@ -19,6 +19,7 @@
 #include <2geom/path.h>
 #include <2geom/transforms.h>
 #include <2geom/d2-sbasis.h>
+#include <2geom/angle.h>
 
 namespace Inkscape {
 namespace LivePathEffect {
diff --git a/src/live_effects/lpe-offset.cpp b/src/live_effects/lpe-offset.cpp
index 68d61603d..027eda23f 100644
--- a/src/live_effects/lpe-offset.cpp
+++ b/src/live_effects/lpe-offset.cpp
@@ -19,7 +19,7 @@
 #include <2geom/path.h>
 #include <2geom/piecewise.h>
 #include <2geom/sbasis-geometric.h>
-#include <2geom/elliptical-arc.h>
+#include <2geom/svg-elliptical-arc.h>
 #include <2geom/transforms.h>
 
 namespace Inkscape {
diff --git a/src/live_effects/lpe-spiro.cpp b/src/live_effects/lpe-spiro.cpp
index 54c15cd96..276ed7b0b 100644
--- a/src/live_effects/lpe-spiro.cpp
+++ b/src/live_effects/lpe-spiro.cpp
@@ -11,6 +11,8 @@
 #include <typeinfo>
 #include <2geom/pathvector.h>
 #include <2geom/matrix.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/hvlinesegment.h>
 #include "helper/geom-nodetype.h"
 
 #include "live_effects/bezctx.h"
diff --git a/src/nodepath.cpp b/src/nodepath.cpp
index 387cf1b0a..c62fb2c0f 100644
--- a/src/nodepath.cpp
+++ b/src/nodepath.cpp
@@ -23,6 +23,8 @@
 #include <glibmm/i18n.h>
 #include <2geom/pathvector.h>
 #include <2geom/sbasis-to-bezier.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/hvlinesegment.h>
 #include "helper/units.h"
 #include "knot.h"
 #include "inkscape.h"
diff --git a/src/sp-path.cpp b/src/sp-path.cpp
index 0867cc3ad..654edcc61 100644
--- a/src/sp-path.cpp
+++ b/src/sp-path.cpp
@@ -25,6 +25,8 @@
 #include <libnr/nr-path.h>
 #include <libnr/nr-matrix-fns.h>
 #include <2geom/pathvector.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/hvlinesegment.h>
 
 #include "svg/svg.h"
 #include "xml/repr.h"
diff --git a/src/sp-polygon.cpp b/src/sp-polygon.cpp
index 04046a47d..fddc9d1db 100644
--- a/src/sp-polygon.cpp
+++ b/src/sp-polygon.cpp
@@ -19,6 +19,8 @@
 #include "display/curve.h"
 #include <glibmm/i18n.h>
 #include <2geom/pathvector.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/hvlinesegment.h>
 #include "svg/stringstream.h"
 #include "xml/repr.h"
 #include "document.h"
diff --git a/src/svg/svg-path.cpp b/src/svg/svg-path.cpp
index b280da07a..12b282577 100644
--- a/src/svg/svg-path.cpp
+++ b/src/svg/svg-path.cpp
@@ -41,11 +41,13 @@
 
 #include <2geom/pathvector.h>
 #include <2geom/path.h>
+#include <2geom/curves.h>
 #include <2geom/sbasis-to-bezier.h>
 #include <2geom/svg-path.h>
 #include <2geom/svg-path-parser.h>
 #include <2geom/exception.h>
 
+
 /* This module parses an SVG path element into an RsvgBpathDef.
 
    At present, there is no support for <marker> or any other contextual