EOL fixup

(bzr r6149)

11 files changed, 4522 insertions, 4522 deletions

diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp
index c9c5b9ec4..20ac0bb2b 100644
--- a/src/2geom/ellipse.cpp
+++ b/src/2geom/ellipse.cpp
@@ -1,206 +1,206 @@
-/*

- * Ellipse Curve

- *

- * Authors:

- *      Marco Cecchetti <mrcekets at gmail.com>

- *

- * Copyright 2008  authors

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-#include <2geom/ellipse.h>

-#include <2geom/svg-elliptical-arc.h>

-#include <2geom/numeric/fitting-tool.h>

-#include <2geom/numeric/fitting-model.h>

-

-

-namespace Geom

-{

-

-void Ellipse::set(double A, double B, double C, double D, double E, double F)

-{

-    double den = 4*A*C - B*B;

-    if ( den == 0 )

-    {

-        THROW_LOGICALERROR("den == 0, while computing ellipse centre");

-    }

-    m_centre[X] = (B*E - 2*C*D) / den;

-    m_centre[Y] = (B*D - 2*A*E) / den;

-

-    // evaluate the a coefficient of the ellipse equation in normal form

-    // E(x,y) = a*(x-cx)^2 + b*(x-cx)*(y-cy) + c*(y-cy)^2 = 1

-    // where b = a*B , c = a*C, (cx,cy) == centre

-    double num =   A * sqr(m_centre[X])

-                 + B * m_centre[X] * m_centre[Y]

-                 + C * sqr(m_centre[Y])

-                 - A * F;

-

-

-    //evaluate ellipse rotation angle

-    double rot = std::atan2( -B, -(A - C) )/2;

-//      std::cerr << "rot = " << rot << std::endl;

-    bool swap_axes = false;

-    if ( are_near(rot, 0) ) rot = 0;

-    if ( are_near(rot, M_PI/2)  || rot < 0 )

-    {

-        swap_axes = true;

-    }

-

-    // evaluate the length of the ellipse rays

-    double cosrot = std::cos(rot);

-    double sinrot = std::sin(rot);

-    double cos2 = cosrot * cosrot;

-    double sin2 = sinrot * sinrot;

-    double cossin = cosrot * sinrot;

-

-    den = A * cos2 + B * cossin + C * sin2;

-    if ( den == 0 )

-    {

-        THROW_LOGICALERROR("den == 0, while computing 'rx' coefficient");

-    }

-    double rx2 =  num/den;

-    if ( rx2 < 0 )

-    {

-        THROW_LOGICALERROR("rx2 < 0, while computing 'rx' coefficient");

-    }

-    double rx = std::sqrt(rx2);

-

-    den = C * cos2 - B * cossin + A * sin2;

-    if ( den == 0 )

-    {

-        THROW_LOGICALERROR("den == 0, while computing 'ry' coefficient");

-    }

-    double ry2 =  num/den;

-    if ( ry2 < 0 )

-    {

-        THROW_LOGICALERROR("ry2 < 0, while computing 'rx' coefficient");

-    }

-    double ry = std::sqrt(ry2);

-

-    // the solution is not unique so we choose always the ellipse

-    // with a rotation angle between 0 and PI/2

-    if ( swap_axes ) std::swap(rx, ry);

-    if (    are_near(rot,  M_PI/2)

-         || are_near(rot, -M_PI/2)

-         || are_near(rx, ry)       )

-    {

-        rot = 0;

-    }

-    else if ( rot < 0 )

-    {

-        rot += M_PI/2;

-    }

-

-    m_ray[X] = rx;

-    m_ray[Y] = ry;

-    m_angle = rot;

-}

-

-

-void Ellipse::set(std::vector<Point> const& points)

-{

-    size_t sz = points.size();

-    if (sz < 5)

-    {

-        THROW_RANGEERROR("fitting error: too few points passed");

-    }

-    NL::LFMEllipse model;

-    NL::least_squeares_fitter<NL::LFMEllipse> fitter(model, sz);

-

-    for (size_t i = 0; i < sz; ++i)

-    {

-        fitter.append(points[i]);

-    }

-    fitter.update();

-

-    NL::Vector z(sz, 0.0);

-    model.instance(*this, fitter.result(z));

-}

-

-

-SVGEllipticalArc

-Ellipse::arc(Point const& initial, Point const& inner, Point const& final,

-             bool _svg_compliant)

-{

-    Point sp_cp = initial - center();

-    Point ep_cp = final   - center();

-    Point ip_cp = inner   - center();

-

-    double angle1 = angle_between(sp_cp, ep_cp);

-    double angle2 = angle_between(sp_cp, ip_cp);

-    double angle3 = angle_between(ip_cp, ep_cp);

-

-    bool large_arc_flag = true;

-    bool sweep_flag = true;

-

-    if ( angle1 > 0 )

-    {

-        if ( angle2 > 0 && angle3 > 0 )

-        {

-            large_arc_flag = false;

-            sweep_flag = true;

-        }

-        else

-        {

-            large_arc_flag = true;

-            sweep_flag = false;

-        }

-    }

-    else

-    {

-        if ( angle2 < 0 && angle3 < 0 )

-        {

-            large_arc_flag = false;

-            sweep_flag = false;

-        }

-        else

-        {

-            large_arc_flag = true;

-            sweep_flag = true;

-        }

-    }

-

-    SVGEllipticalArc ea( initial, ray(X), ray(Y), rot_angle(),

-                      large_arc_flag, sweep_flag, final, _svg_compliant);

-    return ea;

-}

-

-

-}  // end namespace Geom

-

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

-

-

+/*
+ * Ellipse Curve
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include <2geom/ellipse.h>
+#include <2geom/svg-elliptical-arc.h>
+#include <2geom/numeric/fitting-tool.h>
+#include <2geom/numeric/fitting-model.h>
+
+
+namespace Geom
+{
+
+void Ellipse::set(double A, double B, double C, double D, double E, double F)
+{
+    double den = 4*A*C - B*B;
+    if ( den == 0 )
+    {
+        THROW_LOGICALERROR("den == 0, while computing ellipse centre");
+    }
+    m_centre[X] = (B*E - 2*C*D) / den;
+    m_centre[Y] = (B*D - 2*A*E) / den;
+
+    // evaluate the a coefficient of the ellipse equation in normal form
+    // E(x,y) = a*(x-cx)^2 + b*(x-cx)*(y-cy) + c*(y-cy)^2 = 1
+    // where b = a*B , c = a*C, (cx,cy) == centre
+    double num =   A * sqr(m_centre[X])
+                 + B * m_centre[X] * m_centre[Y]
+                 + C * sqr(m_centre[Y])
+                 - A * F;
+
+
+    //evaluate ellipse rotation angle
+    double rot = std::atan2( -B, -(A - C) )/2;
+//      std::cerr << "rot = " << rot << std::endl;
+    bool swap_axes = false;
+    if ( are_near(rot, 0) ) rot = 0;
+    if ( are_near(rot, M_PI/2)  || rot < 0 )
+    {
+        swap_axes = true;
+    }
+
+    // evaluate the length of the ellipse rays
+    double cosrot = std::cos(rot);
+    double sinrot = std::sin(rot);
+    double cos2 = cosrot * cosrot;
+    double sin2 = sinrot * sinrot;
+    double cossin = cosrot * sinrot;
+
+    den = A * cos2 + B * cossin + C * sin2;
+    if ( den == 0 )
+    {
+        THROW_LOGICALERROR("den == 0, while computing 'rx' coefficient");
+    }
+    double rx2 =  num/den;
+    if ( rx2 < 0 )
+    {
+        THROW_LOGICALERROR("rx2 < 0, while computing 'rx' coefficient");
+    }
+    double rx = std::sqrt(rx2);
+
+    den = C * cos2 - B * cossin + A * sin2;
+    if ( den == 0 )
+    {
+        THROW_LOGICALERROR("den == 0, while computing 'ry' coefficient");
+    }
+    double ry2 =  num/den;
+    if ( ry2 < 0 )
+    {
+        THROW_LOGICALERROR("ry2 < 0, while computing 'rx' coefficient");
+    }
+    double ry = std::sqrt(ry2);
+
+    // the solution is not unique so we choose always the ellipse
+    // with a rotation angle between 0 and PI/2
+    if ( swap_axes ) std::swap(rx, ry);
+    if (    are_near(rot,  M_PI/2)
+         || are_near(rot, -M_PI/2)
+         || are_near(rx, ry)       )
+    {
+        rot = 0;
+    }
+    else if ( rot < 0 )
+    {
+        rot += M_PI/2;
+    }
+
+    m_ray[X] = rx;
+    m_ray[Y] = ry;
+    m_angle = rot;
+}
+
+
+void Ellipse::set(std::vector<Point> const& points)
+{
+    size_t sz = points.size();
+    if (sz < 5)
+    {
+        THROW_RANGEERROR("fitting error: too few points passed");
+    }
+    NL::LFMEllipse model;
+    NL::least_squeares_fitter<NL::LFMEllipse> fitter(model, sz);
+
+    for (size_t i = 0; i < sz; ++i)
+    {
+        fitter.append(points[i]);
+    }
+    fitter.update();
+
+    NL::Vector z(sz, 0.0);
+    model.instance(*this, fitter.result(z));
+}
+
+
+SVGEllipticalArc
+Ellipse::arc(Point const& initial, Point const& inner, Point const& final,
+             bool _svg_compliant)
+{
+    Point sp_cp = initial - center();
+    Point ep_cp = final   - center();
+    Point ip_cp = inner   - center();
+
+    double angle1 = angle_between(sp_cp, ep_cp);
+    double angle2 = angle_between(sp_cp, ip_cp);
+    double angle3 = angle_between(ip_cp, ep_cp);
+
+    bool large_arc_flag = true;
+    bool sweep_flag = true;
+
+    if ( angle1 > 0 )
+    {
+        if ( angle2 > 0 && angle3 > 0 )
+        {
+            large_arc_flag = false;
+            sweep_flag = true;
+        }
+        else
+        {
+            large_arc_flag = true;
+            sweep_flag = false;
+        }
+    }
+    else
+    {
+        if ( angle2 < 0 && angle3 < 0 )
+        {
+            large_arc_flag = false;
+            sweep_flag = false;
+        }
+        else
+        {
+            large_arc_flag = true;
+            sweep_flag = true;
+        }
+    }
+
+    SVGEllipticalArc ea( initial, ray(X), ray(Y), rot_angle(),
+                      large_arc_flag, sweep_flag, final, _svg_compliant);
+    return ea;
+}
+
+
+}  // end namespace Geom
+
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+
+
diff --git a/src/2geom/ellipse.h b/src/2geom/ellipse.h
index af8b01e78..a10f9ced0 100644
--- a/src/2geom/ellipse.h
+++ b/src/2geom/ellipse.h
@@ -1,133 +1,133 @@
-/*

- * Ellipse Curve

- *

- * Authors:

- *      Marco Cecchetti <mrcekets at gmail.com>

- *

- * Copyright 2008  authors

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-#ifndef _2GEOM_ELLIPSE_H_

-#define _2GEOM_ELLIPSE_H_

-

-

-#include <2geom/point.h>

-#include <2geom/exception.h>

-

-

-namespace Geom

-{

-

-class SVGEllipticalArc;

-

-class Ellipse

-{

-  public:

-    Ellipse()

-    {}

-

-    Ellipse(double cx, double cy, double rx, double ry, double a)

-        : m_centre(cx, cy), m_ray(rx, ry), m_angle(a)

-    {

-    }

-

-    // build an ellipse by its implicit equation:

-    // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0

-    Ellipse(double A, double B, double C, double D, double E, double F)

-    {

-        set(A, B, C, D, E, F);

-    }

-

-    Ellipse(std::vector<Point> const& points)

-    {

-        set(points);

-    }

-

-    void set(double cx, double cy, double rx, double ry, double a)

-    {

-        m_centre[X] = cx;

-        m_centre[Y] = cy;

-        m_ray[X] = rx;

-        m_ray[Y] = ry;

-        m_angle = a;

-    }

-

-    // build an ellipse by its implicit equation:

-    // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0

-    void set(double A, double B, double C, double D, double E, double F);

-

-    // biuld up the best fitting ellipse wrt the passed points

-    // prerequisite: at least 5 points must be passed

-    void set(std::vector<Point> const& points);

-

-    SVGEllipticalArc

-    arc(Point const& initial, Point const& inner, Point const& final,

-        bool _svg_compliant = true);

-

-    Point center() const

-    {

-        return m_centre;

-    }

-

-    Coord center(Dim2 d) const

-    {

-        return m_centre[d];

-    }

-

-    Coord ray(Dim2 d) const

-    {

-        return m_ray[d];

-    }

-

-    Coord rot_angle() const

-    {

-        return m_angle;

-    }

-

-  private:

-    Point m_centre, m_ray;

-    double m_angle;

-};

-

-

-} // end namespace Geom

-

-

-

-#endif // _2GEOM_ELLIPSE_H_

-

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

+/*
+ * Ellipse Curve
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _2GEOM_ELLIPSE_H_
+#define _2GEOM_ELLIPSE_H_
+
+
+#include <2geom/point.h>
+#include <2geom/exception.h>
+
+
+namespace Geom
+{
+
+class SVGEllipticalArc;
+
+class Ellipse
+{
+  public:
+    Ellipse()
+    {}
+
+    Ellipse(double cx, double cy, double rx, double ry, double a)
+        : m_centre(cx, cy), m_ray(rx, ry), m_angle(a)
+    {
+    }
+
+    // build an ellipse by its implicit equation:
+    // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
+    Ellipse(double A, double B, double C, double D, double E, double F)
+    {
+        set(A, B, C, D, E, F);
+    }
+
+    Ellipse(std::vector<Point> const& points)
+    {
+        set(points);
+    }
+
+    void set(double cx, double cy, double rx, double ry, double a)
+    {
+        m_centre[X] = cx;
+        m_centre[Y] = cy;
+        m_ray[X] = rx;
+        m_ray[Y] = ry;
+        m_angle = a;
+    }
+
+    // build an ellipse by its implicit equation:
+    // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0
+    void set(double A, double B, double C, double D, double E, double F);
+
+    // biuld up the best fitting ellipse wrt the passed points
+    // prerequisite: at least 5 points must be passed
+    void set(std::vector<Point> const& points);
+
+    SVGEllipticalArc
+    arc(Point const& initial, Point const& inner, Point const& final,
+        bool _svg_compliant = true);
+
+    Point center() const
+    {
+        return m_centre;
+    }
+
+    Coord center(Dim2 d) const
+    {
+        return m_centre[d];
+    }
+
+    Coord ray(Dim2 d) const
+    {
+        return m_ray[d];
+    }
+
+    Coord rot_angle() const
+    {
+        return m_angle;
+    }
+
+  private:
+    Point m_centre, m_ray;
+    double m_angle;
+};
+
+
+} // end namespace Geom
+
+
+
+#endif // _2GEOM_ELLIPSE_H_
+
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/nearest-point.cpp b/src/2geom/nearest-point.cpp
index 59d55ba32..7fef8c120 100644
--- a/src/2geom/nearest-point.cpp
+++ b/src/2geom/nearest-point.cpp
@@ -1,278 +1,278 @@
-/*

- * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>

- *

- * Authors:

- * 		

- * 		Marco Cecchetti <mrcekets at gmail.com>

- * 

- * Copyright 2007-2008  authors

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-#include <2geom/nearest-point.h>

-

-namespace Geom

-{

-

-////////////////////////////////////////////////////////////////////////////////

-// D2<SBasis> versions

-

-/*

- * Return the parameter t of a nearest point on the portion of the curve "c", 

- * related to the interval [from, to], to the point "p".

- * The needed curve derivative "dc" is passed as parameter.

- * The function return the first nearest point to "p" that is found.

- */

-

-double nearest_point( Point const& p, 

-					  D2<SBasis> const& c, 

-					  D2<SBasis> const& dc, 

-		              double from, double to )

-{

-	if ( from > to ) std::swap(from, to);

-	if ( from < 0 || to > 1 )

-	{

-		THROW_RANGEERROR("[from,to] interval out of bounds");

-	}

-

-	SBasis dd = dot(c - p, dc);	

-	std::vector<double> zeros = Geom::roots(dd);

-	

-	double closest = from;

-	double min_dist_sq = L2sq(c(from) - p);

-	double distsq;

-	for ( unsigned int i = 0; i < zeros.size(); ++i )

-	{

-		distsq = L2sq(c(zeros[i]) - p);

-		if ( min_dist_sq > L2sq(c(zeros[i]) - p) )

-		{

-			closest = zeros[i];

-			min_dist_sq = distsq;

-		}

-	}

-	if ( min_dist_sq > L2sq( c(to) - p ) )

-		closest = to;

-	return closest;

-

-}

-

-/*

- * Return the parameters t of all the nearest points on the portion of 

- * the curve "c", related to the interval [from, to], to the point "p".

- * The needed curve derivative "dc" is passed as parameter.

- */

-

-std::vector<double> 

-all_nearest_points( Point const& p, 

-		    D2<SBasis> const& c, 

-		    D2<SBasis> const& /*dc*/, 

-		    double from, double to )

-{

-	std::swap(from, to);

-	if ( from > to ) std::swap(from, to);

-	if ( from < 0 || to > 1 )

-	{

-		THROW_RANGEERROR("[from,to] interval out of bounds");

-	}

-

-	std::vector<double> result;

-	SBasis dd = dot(c - p, Geom::derivative(c));

-	

-	std::vector<double> zeros = Geom::roots(dd);

-	std::vector<double> candidates;

-	candidates.push_back(from);

-	candidates.insert(candidates.end(), zeros.begin(), zeros.end());

-	candidates.push_back(to);

-	std::vector<double> distsq;

-	distsq.reserve(candidates.size());

-	for ( unsigned int i = 0; i < candidates.size(); ++i )

-	{

-		distsq.push_back( L2sq(c(candidates[i]) - p) );

-	}

-	unsigned int closest = 0;

-	double dsq = distsq[0];

-	for ( unsigned int i = 1; i < candidates.size(); ++i )

-	{

-		if ( dsq > distsq[i] )

-		{

-			closest = i;

-			dsq = distsq[i];

-		}

-	}

-	for ( unsigned int i = 0; i < candidates.size(); ++i )

-	{

-		if( distsq[closest] == distsq[i] )

-		{

-			result.push_back(candidates[i]);

-		}

-	}

-	return result;

-}

-

-

-////////////////////////////////////////////////////////////////////////////////

-// Piecewise< D2<SBasis> > versions

-

-

-double nearest_point( Point const& p,  

-					  Piecewise< D2<SBasis> > const& c,

-		              double from, double to )

-{

-	if ( from > to ) std::swap(from, to);

-	if ( from < c.cuts[0] || to > c.cuts[c.size()] )

-	{

-		THROW_RANGEERROR("[from,to] interval out of bounds");

-	}

-	

-	unsigned int si = c.segN(from);

-	unsigned int ei = c.segN(to);

-	if ( si == ei )

-	{

-		double nearest=

-			nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));

-		return c.mapToDomain(nearest, si);

-	}

-	double t;

-	double nearest = nearest_point(p, c[si], c.segT(from, si));

-	unsigned int ni = si;

-	double dsq;

-	double mindistsq = distanceSq(p, c[si](nearest));

-	Rect bb;

-	for ( unsigned int i = si + 1; i < ei; ++i )

-	{

-		bb = bounds_fast(c[i]);

-		dsq = distanceSq(p, bb);

-		if ( mindistsq <= dsq ) continue;

-		t = nearest_point(p, c[i]);

-		dsq = distanceSq(p, c[i](t));

-		if ( mindistsq > dsq )

-		{

-			nearest = t;

-			ni = i;

-			mindistsq = dsq;

-		}

-	}

-	bb = bounds_fast(c[ei]);

-	dsq = distanceSq(p, bb);

-	if ( mindistsq > dsq )

-	{

-		t = nearest_point(p, c[ei], 0, c.segT(to, ei));

-		dsq = distanceSq(p, c[ei](t));

-		if ( mindistsq > dsq )

-		{

-			nearest = t;

-			ni = ei;

-		}

-	}

-	return c.mapToDomain(nearest, ni);

-}

-

-std::vector<double> 

-all_nearest_points( Point const& p, 

-                    Piecewise< D2<SBasis> > const& c, 

-		            double from, double to )

-{

-	if ( from > to ) std::swap(from, to);

-	if ( from < c.cuts[0] || to > c.cuts[c.size()] )

-	{

-		THROW_RANGEERROR("[from,to] interval out of bounds");

-	}

-	

-	unsigned int si = c.segN(from);

-	unsigned int ei = c.segN(to);

-	if ( si == ei )

-	{

-		std::vector<double>	all_nearest = 

-			all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));

-		for ( unsigned int i = 0; i < all_nearest.size(); ++i )

-		{

-			all_nearest[i] = c.mapToDomain(all_nearest[i], si);

-		}

-		return all_nearest;

-	}

-	std::vector<double> all_t;

-	std::vector< std::vector<double> > all_np;

-	all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );

-	std::vector<unsigned int> ni;

-	ni.push_back(si);

-	double dsq;

-	double mindistsq = distanceSq( p, c[si](all_np.front().front()) );

-	Rect bb;

-	for ( unsigned int i = si + 1; i < ei; ++i )

-	{

-		bb = bounds_fast(c[i]);

-		dsq = distanceSq(p, bb);

-		if ( mindistsq < dsq ) continue;

-		all_t = all_nearest_points(p, c[i]);

-		dsq = distanceSq( p, c[i](all_t.front()) );

-		if ( mindistsq > dsq )

-		{

-			all_np.clear();

-			all_np.push_back(all_t);

-			ni.clear();

-			ni.push_back(i);

-			mindistsq = dsq;

-		}

-		else if ( mindistsq == dsq )

-		{

-			all_np.push_back(all_t);

-			ni.push_back(i);

-		}

-	}

-	bb = bounds_fast(c[ei]);

-	dsq = distanceSq(p, bb);

-	if ( mindistsq >= dsq )

-	{

-		all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));

-		dsq = distanceSq( p, c[ei](all_t.front()) );

-		if ( mindistsq > dsq )

-		{

-			for ( unsigned int i = 0; i < all_t.size(); ++i )

-			{

-				all_t[i] = c.mapToDomain(all_t[i], ei);

-			}

-			return all_t;

-		}

-		else if ( mindistsq == dsq )

-		{

-			all_np.push_back(all_t);

-			ni.push_back(ei);

-		}

-	}

-	std::vector<double> all_nearest;

-	for ( unsigned int i = 0; i < all_np.size(); ++i )

-	{

-		for ( unsigned int j = 0; j < all_np[i].size(); ++j )

-		{

-			all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) );

-		}

-	}

-	return all_nearest;

-}

-

-} // end namespace Geom

-

-

+/*
+ * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>
+ *
+ * Authors:
+ * 		
+ * 		Marco Cecchetti <mrcekets at gmail.com>
+ * 
+ * Copyright 2007-2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include <2geom/nearest-point.h>
+
+namespace Geom
+{
+
+////////////////////////////////////////////////////////////////////////////////
+// D2<SBasis> versions
+
+/*
+ * Return the parameter t of a nearest point on the portion of the curve "c", 
+ * related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ * The function return the first nearest point to "p" that is found.
+ */
+
+double nearest_point( Point const& p, 
+					  D2<SBasis> const& c, 
+					  D2<SBasis> const& dc, 
+		              double from, double to )
+{
+	if ( from > to ) std::swap(from, to);
+	if ( from < 0 || to > 1 )
+	{
+		THROW_RANGEERROR("[from,to] interval out of bounds");
+	}
+
+	SBasis dd = dot(c - p, dc);	
+	std::vector<double> zeros = Geom::roots(dd);
+	
+	double closest = from;
+	double min_dist_sq = L2sq(c(from) - p);
+	double distsq;
+	for ( unsigned int i = 0; i < zeros.size(); ++i )
+	{
+		distsq = L2sq(c(zeros[i]) - p);
+		if ( min_dist_sq > L2sq(c(zeros[i]) - p) )
+		{
+			closest = zeros[i];
+			min_dist_sq = distsq;
+		}
+	}
+	if ( min_dist_sq > L2sq( c(to) - p ) )
+		closest = to;
+	return closest;
+
+}
+
+/*
+ * Return the parameters t of all the nearest points on the portion of 
+ * the curve "c", related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ */
+
+std::vector<double> 
+all_nearest_points( Point const& p, 
+		    D2<SBasis> const& c, 
+		    D2<SBasis> const& /*dc*/, 
+		    double from, double to )
+{
+	std::swap(from, to);
+	if ( from > to ) std::swap(from, to);
+	if ( from < 0 || to > 1 )
+	{
+		THROW_RANGEERROR("[from,to] interval out of bounds");
+	}
+
+	std::vector<double> result;
+	SBasis dd = dot(c - p, Geom::derivative(c));
+	
+	std::vector<double> zeros = Geom::roots(dd);
+	std::vector<double> candidates;
+	candidates.push_back(from);
+	candidates.insert(candidates.end(), zeros.begin(), zeros.end());
+	candidates.push_back(to);
+	std::vector<double> distsq;
+	distsq.reserve(candidates.size());
+	for ( unsigned int i = 0; i < candidates.size(); ++i )
+	{
+		distsq.push_back( L2sq(c(candidates[i]) - p) );
+	}
+	unsigned int closest = 0;
+	double dsq = distsq[0];
+	for ( unsigned int i = 1; i < candidates.size(); ++i )
+	{
+		if ( dsq > distsq[i] )
+		{
+			closest = i;
+			dsq = distsq[i];
+		}
+	}
+	for ( unsigned int i = 0; i < candidates.size(); ++i )
+	{
+		if( distsq[closest] == distsq[i] )
+		{
+			result.push_back(candidates[i]);
+		}
+	}
+	return result;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Piecewise< D2<SBasis> > versions
+
+
+double nearest_point( Point const& p,  
+					  Piecewise< D2<SBasis> > const& c,
+		              double from, double to )
+{
+	if ( from > to ) std::swap(from, to);
+	if ( from < c.cuts[0] || to > c.cuts[c.size()] )
+	{
+		THROW_RANGEERROR("[from,to] interval out of bounds");
+	}
+	
+	unsigned int si = c.segN(from);
+	unsigned int ei = c.segN(to);
+	if ( si == ei )
+	{
+		double nearest=
+			nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));
+		return c.mapToDomain(nearest, si);
+	}
+	double t;
+	double nearest = nearest_point(p, c[si], c.segT(from, si));
+	unsigned int ni = si;
+	double dsq;
+	double mindistsq = distanceSq(p, c[si](nearest));
+	Rect bb;
+	for ( unsigned int i = si + 1; i < ei; ++i )
+	{
+		bb = bounds_fast(c[i]);
+		dsq = distanceSq(p, bb);
+		if ( mindistsq <= dsq ) continue;
+		t = nearest_point(p, c[i]);
+		dsq = distanceSq(p, c[i](t));
+		if ( mindistsq > dsq )
+		{
+			nearest = t;
+			ni = i;
+			mindistsq = dsq;
+		}
+	}
+	bb = bounds_fast(c[ei]);
+	dsq = distanceSq(p, bb);
+	if ( mindistsq > dsq )
+	{
+		t = nearest_point(p, c[ei], 0, c.segT(to, ei));
+		dsq = distanceSq(p, c[ei](t));
+		if ( mindistsq > dsq )
+		{
+			nearest = t;
+			ni = ei;
+		}
+	}
+	return c.mapToDomain(nearest, ni);
+}
+
+std::vector<double> 
+all_nearest_points( Point const& p, 
+                    Piecewise< D2<SBasis> > const& c, 
+		            double from, double to )
+{
+	if ( from > to ) std::swap(from, to);
+	if ( from < c.cuts[0] || to > c.cuts[c.size()] )
+	{
+		THROW_RANGEERROR("[from,to] interval out of bounds");
+	}
+	
+	unsigned int si = c.segN(from);
+	unsigned int ei = c.segN(to);
+	if ( si == ei )
+	{
+		std::vector<double>	all_nearest = 
+			all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));
+		for ( unsigned int i = 0; i < all_nearest.size(); ++i )
+		{
+			all_nearest[i] = c.mapToDomain(all_nearest[i], si);
+		}
+		return all_nearest;
+	}
+	std::vector<double> all_t;
+	std::vector< std::vector<double> > all_np;
+	all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );
+	std::vector<unsigned int> ni;
+	ni.push_back(si);
+	double dsq;
+	double mindistsq = distanceSq( p, c[si](all_np.front().front()) );
+	Rect bb;
+	for ( unsigned int i = si + 1; i < ei; ++i )
+	{
+		bb = bounds_fast(c[i]);
+		dsq = distanceSq(p, bb);
+		if ( mindistsq < dsq ) continue;
+		all_t = all_nearest_points(p, c[i]);
+		dsq = distanceSq( p, c[i](all_t.front()) );
+		if ( mindistsq > dsq )
+		{
+			all_np.clear();
+			all_np.push_back(all_t);
+			ni.clear();
+			ni.push_back(i);
+			mindistsq = dsq;
+		}
+		else if ( mindistsq == dsq )
+		{
+			all_np.push_back(all_t);
+			ni.push_back(i);
+		}
+	}
+	bb = bounds_fast(c[ei]);
+	dsq = distanceSq(p, bb);
+	if ( mindistsq >= dsq )
+	{
+		all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));
+		dsq = distanceSq( p, c[ei](all_t.front()) );
+		if ( mindistsq > dsq )
+		{
+			for ( unsigned int i = 0; i < all_t.size(); ++i )
+			{
+				all_t[i] = c.mapToDomain(all_t[i], ei);
+			}
+			return all_t;
+		}
+		else if ( mindistsq == dsq )
+		{
+			all_np.push_back(all_t);
+			ni.push_back(ei);
+		}
+	}
+	std::vector<double> all_nearest;
+	for ( unsigned int i = 0; i < all_np.size(); ++i )
+	{
+		for ( unsigned int j = 0; j < all_np[i].size(); ++j )
+		{
+			all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) );
+		}
+	}
+	return all_nearest;
+}
+
+} // end namespace Geom
+
+
diff --git a/src/2geom/nearest-point.h b/src/2geom/nearest-point.h
index 0b43ce67b..c8588214b 100644
--- a/src/2geom/nearest-point.h
+++ b/src/2geom/nearest-point.h
@@ -1,133 +1,133 @@
-/*

- * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>

- *

- *

- * Authors:

- * 		

- * 		Marco Cecchetti <mrcekets at gmail.com>

- * 

- * Copyright 2007-2008  authors

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-#ifndef _NEAREST_POINT_H_

-#define _NEAREST_POINT_H_

-

-

-#include <vector>

-

-#include <2geom/d2.h>

-#include <2geom/piecewise.h>

-#include <2geom/exception.h>

-

-

-

-namespace Geom

-{

-

-/*

- * Given a line L specified by a point A and direction vector v,

- * return the point on L nearest to p. Note that the returned value

- * is with respect to the _normalized_ direction of v!

- */

-inline double nearest_point(Point const &p, Point const &A, Point const &v)

-{

-    Point d(p - A);

-    return d[0] * v[0] + d[1] * v[1];

-}

-

-////////////////////////////////////////////////////////////////////////////////

-// D2<SBasis> versions

-

-/*

- * Return the parameter t of a nearest point on the portion of the curve "c", 

- * related to the interval [from, to], to the point "p".

- * The needed curve derivative "dc" is passed as parameter.

- * The function return the first nearest point to "p" that is found.

- */

-double nearest_point( Point const& p,

-		              D2<SBasis> const& c, D2<SBasis> const& dc, 

-		              double from = 0, double to = 1 );

-

-inline

-double nearest_point( Point const& p, 

-		              D2<SBasis> const& c, 

-		              double from = 0, double to = 1 )

-{

-	return nearest_point(p, c, Geom::derivative(c), from, to);

-}

-

-/*

- * Return the parameters t of all the nearest points on the portion of 

- * the curve "c", related to the interval [from, to], to the point "p".

- * The needed curve derivative "dc" is passed as parameter.

- */

-std::vector<double> 

-all_nearest_points( Point const& p, 

-		            D2<SBasis> const& c, D2<SBasis> const& dc, 

-		            double from = 0, double to = 1 );

-

-inline

-std::vector<double> 

-all_nearest_points( Point const& p, 

-		            D2<SBasis> const& c, 

-		            double from = 0, double to = 1 )

-{

-	return all_nearest_points(p, c,  Geom::derivative(c), from, to);

-}

-

-

-////////////////////////////////////////////////////////////////////////////////

-// Piecewise< D2<SBasis> > versions

-

-double nearest_point( Point const& p, 

-		              Piecewise< D2<SBasis> > const& c, 

-		              double from, double to );

-

-inline

-double nearest_point( Point const& p, Piecewise< D2<SBasis> > const& c ) 

-{

-	return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]);

-}

-

-

-std::vector<double> 

-all_nearest_points( Point const& p, 

-					Piecewise< D2<SBasis> > const& c, 

-		            double from, double to );

-

-inline

-std::vector<double> 

-all_nearest_points( Point const& p, Piecewise< D2<SBasis> > const& c ) 

-{

-	return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]);

-}

-

-} // end namespace Geom

-

-

-

-#endif /*_NEAREST_POINT_H_*/

+/*
+ * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>
+ *
+ *
+ * Authors:
+ * 		
+ * 		Marco Cecchetti <mrcekets at gmail.com>
+ * 
+ * Copyright 2007-2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NEAREST_POINT_H_
+#define _NEAREST_POINT_H_
+
+
+#include <vector>
+
+#include <2geom/d2.h>
+#include <2geom/piecewise.h>
+#include <2geom/exception.h>
+
+
+
+namespace Geom
+{
+
+/*
+ * Given a line L specified by a point A and direction vector v,
+ * return the point on L nearest to p. Note that the returned value
+ * is with respect to the _normalized_ direction of v!
+ */
+inline double nearest_point(Point const &p, Point const &A, Point const &v)
+{
+    Point d(p - A);
+    return d[0] * v[0] + d[1] * v[1];
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// D2<SBasis> versions
+
+/*
+ * Return the parameter t of a nearest point on the portion of the curve "c", 
+ * related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ * The function return the first nearest point to "p" that is found.
+ */
+double nearest_point( Point const& p,
+		              D2<SBasis> const& c, D2<SBasis> const& dc, 
+		              double from = 0, double to = 1 );
+
+inline
+double nearest_point( Point const& p, 
+		              D2<SBasis> const& c, 
+		              double from = 0, double to = 1 )
+{
+	return nearest_point(p, c, Geom::derivative(c), from, to);
+}
+
+/*
+ * Return the parameters t of all the nearest points on the portion of 
+ * the curve "c", related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ */
+std::vector<double> 
+all_nearest_points( Point const& p, 
+		            D2<SBasis> const& c, D2<SBasis> const& dc, 
+		            double from = 0, double to = 1 );
+
+inline
+std::vector<double> 
+all_nearest_points( Point const& p, 
+		            D2<SBasis> const& c, 
+		            double from = 0, double to = 1 )
+{
+	return all_nearest_points(p, c,  Geom::derivative(c), from, to);
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Piecewise< D2<SBasis> > versions
+
+double nearest_point( Point const& p, 
+		              Piecewise< D2<SBasis> > const& c, 
+		              double from, double to );
+
+inline
+double nearest_point( Point const& p, Piecewise< D2<SBasis> > const& c ) 
+{
+	return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]);
+}
+
+
+std::vector<double> 
+all_nearest_points( Point const& p, 
+					Piecewise< D2<SBasis> > const& c, 
+		            double from, double to );
+
+inline
+std::vector<double> 
+all_nearest_points( Point const& p, Piecewise< D2<SBasis> > const& c ) 
+{
+	return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]);
+}
+
+} // end namespace Geom
+
+
+
+#endif /*_NEAREST_POINT_H_*/
diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h
index 145be40e4..cc3113372 100644
--- a/src/2geom/numeric/fitting-model.h
+++ b/src/2geom/numeric/fitting-model.h
@@ -1,423 +1,423 @@
-/*

- * Fitting Models for Geom Types

- *

- * Authors:

- *      Marco Cecchetti <mrcekets at gmail.com>

- *

- * Copyright 2008  authors

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-#ifndef _NL_FITTING_MODEL_H_

-#define _NL_FITTING_MODEL_H_

-

-

-#include <2geom/d2.h>

-#include <2geom/sbasis.h>

-#include <2geom/bezier.h>

-#include <2geom/bezier-curve.h>

-#include <2geom/poly.h>

-#include <2geom/ellipse.h>

-#include <2geom/utils.h>

-

-

-namespace Geom { namespace NL {

-

-

-/*

- *   completely unknown models must inherit from this template class;

- *   example: the model a*x^2 + b*x + c = 0 to be solved wrt a, b, c;

- *   example: the model A(t) = known_sample_value_at(t) to be solved wrt

- *       the coefficients of the curve A(t) expressed in S-Basis form;

- *   parameter type: the type of x and t variable in the examples above;

- *   value type:     the type of the known sample values (in the first example

- *                   is constant )

- *   instance type:  the type of the objects produced by using

- *                   the fitting raw data solution

- */

-template< typename ParameterType, typename ValueType, typename InstanceType >

-class LinearFittingModel

-{

-  public:

-    typedef ParameterType       parameter_type;

-    typedef ValueType           value_type;

-    typedef InstanceType        instance_type;

-

-    static const bool WITH_FIXED_TERMS = false;

-

-    /*

-     * a LinearFittingModel must implement the following methods:

-     *

-     * void feed( VectorView & vector,

-     *            parameter_type const& sample_parameter ) const;

-     *

-     * size_t size() const;

-     *

-     * void instance(instance_type &, raw_type const& raw_data) const;

-     *

-     */

-};

-

-

-/*

- *   partially known models must inherit from this template class

- *   example: the model a*x^2 + 2*x + c = 0 to be solved wrt a and c

- */

-template< typename ParameterType, typename ValueType, typename InstanceType >

-class LinearFittingModelWithFixedTerms

-{

-  public:

-    typedef ParameterType       parameter_type;

-    typedef ValueType           value_type;

-    typedef InstanceType        instance_type;

-

-    static const bool WITH_FIXED_TERMS = true;

-

-    /*

-     * a LinearFittingModelWithFixedTerms must implement the following methods:

-     *

-     * void feed( VectorView & vector,

-     *            value_type & fixed_term,

-     *            parameter_type const& sample_parameter ) const;

-     *

-     * size_t size() const;

-     *

-     * void instance(instance_type &, raw_type const& raw_data) const;

-     *

-     */

-

-

-};

-

-

-// incomplete model, it can be inherited to make up different kinds of

-// instance type; the raw data is a vector of coefficients of a polynomial

-// rapresented in standard power basis

-template< typename InstanceType >

-class LFMPowerBasis

-    : public LinearFittingModel<double, double, InstanceType>

-{

-  public:

-    LFMPowerBasis(size_t degree)

-        : m_size(degree + 1)

-    {

-    }

-

-    void feed( VectorView & coeff, double sample_parameter ) const

-    {

-        coeff[0] = 1;

-        double x_i = 1;

-        for (size_t i = 1; i < coeff.size(); ++i)

-        {

-          x_i *= sample_parameter;

-          coeff[i] = x_i;

-        }

-    }

-

-    size_t size() const

-    {

-        return m_size;

-    }

-

-  private:

-    size_t m_size;

-};

-

-

-// this model generates Geom::Poly objects

-class LFMPoly

-    : public LFMPowerBasis<Poly>

-{

-  public:

-    LFMPoly(size_t degree)

-        : LFMPowerBasis<Poly>(degree)

-    {

-    }

-

-    void instance(Poly & poly, ConstVectorView const& raw_data) const

-    {

-        poly.clear();

-        poly.resize(size());

-        for (size_t i = 0; i < raw_data.size(); ++i)

-        {

-            poly[i] =  raw_data[i];

-        }

-    }

-};

-

-

-// incomplete model, it can be inherited to make up different kinds of

-// instance type; the raw data is a vector of coefficients of a polynomial

-// rapresented in standard power basis with leading term coefficient equal to 1

-template< typename InstanceType >

-class LFMNormalizedPowerBasis

-    : public LinearFittingModelWithFixedTerms<double, double, InstanceType>

-{

-  public:

-    LFMNormalizedPowerBasis(size_t _degree)

-        : m_model( _degree - 1)

-    {

-        assert(_degree > 0);

-    }

-

-

-    void feed( VectorView & coeff,

-               double & known_term,

-               double sample_parameter ) const

-    {

-        m_model.feed(coeff, sample_parameter);

-        known_term = coeff[m_model.size()-1] * sample_parameter;

-    }

-

-    size_t size() const

-    {

-        return m_model.size();

-    }

-

-  private:

-    LFMPowerBasis<InstanceType> m_model;

-};

-

-

-// incomplete model, it can be inherited to make up different kinds of

-// instance type; the raw data is a vector of coefficients of the equation

-// of an ellipse curve

-template< typename InstanceType >

-class LFMEllipseEquation

-    : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>

-{

-  public:

-    void feed( VectorView & coeff, double & fixed_term, Point const& p ) const

-    {

-        coeff[0] = p[X] * p[Y];

-        coeff[1] = p[Y] * p[Y];

-        coeff[2] = p[X];

-        coeff[3] = p[Y];

-        coeff[4] = 1;

-        fixed_term = p[X] * p[X];

-    }

-

-    size_t size() const

-    {

-        return 5;

-    }

-};

-

-

-// this model generates Ellipse curves

-class LFMEllipse

-    : public LFMEllipseEquation<Ellipse>

-{

-  public:

-    void instance(Ellipse & e, ConstVectorView const& coeff) const

-    {

-        e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);

-    }

-};

-

-

-// this model generates SBasis objects

-class LFMSBasis

-    : public LinearFittingModel<double, double, SBasis>

-{

-  public:

-    LFMSBasis( size_t _order )

-        : m_size( 2*(_order+1) ),

-          m_order(_order)

-    {

-    }

-

-    void feed( VectorView & coeff, double t ) const

-    {

-        double u0 = 1-t;

-        double u1 = t;

-        double s = u0 * u1;

-        coeff[0] = u0;

-        coeff[1] = u1;

-        for (size_t i = 2; i < size(); i+=2)

-        {

-            u0 *= s;

-            u1 *= s;

-            coeff[i] = u0;

-            coeff[i+1] = u1;

-        }

-    }

-

-    size_t size() const

-    {

-        return m_size;

-    }

-

-    void instance(SBasis & sb, ConstVectorView const& raw_data) const

-    {

-        sb.clear();

-        sb.resize(m_order+1);

-        for (unsigned int i = 0, k = 0; i < raw_data.size(); i+=2, ++k)

-        {

-            sb[k][0] = raw_data[i];

-            sb[k][1] = raw_data[i+1];

-        }

-    }

-

-  private:

-    size_t m_size;

-    size_t m_order;

-};

-

-

-// this model generates D2<SBasis> objects

-class LFMD2SBasis

-    : public LinearFittingModel< double, Point, D2<SBasis> >

-{

-  public:

-    LFMD2SBasis( size_t _order )

-        : mosb(_order)

-    {

-    }

-

-    void feed( VectorView & coeff, double t ) const

-    {

-        mosb.feed(coeff, t);

-    }

-

-    size_t size() const

-    {

-        return mosb.size();

-    }

-

-    void instance(D2<SBasis> & d2sb, ConstMatrixView const& raw_data) const

-    {

-        mosb.instance(d2sb[X], raw_data.column_const_view(X));

-        mosb.instance(d2sb[Y], raw_data.column_const_view(Y));

-    }

-

-  private:

-    LFMSBasis mosb;

-};

-

-

-// this model generates Bezier objects

-class LFMBezier

-    : public LinearFittingModel<double, double, Bezier>

-{

-  public:

-    LFMBezier( size_t _order )

-        : m_size(_order + 1),

-          m_order(_order)

-    {

-        binomial_coefficients(m_bc, m_order);

-    }

-

-    void feed( VectorView & coeff, double t ) const

-    {

-        double s = 1;

-        for (size_t i = 0; i < size(); ++i)

-        {

-            coeff[i] = s * m_bc[i];

-            s *= t;

-        }

-        double u = 1-t;

-        s = 1;

-        for (size_t i = size()-1; i > 0; --i)

-        {

-            coeff[i] *= s;

-            s *= u;

-        }

-        coeff[0] *= s;

-    }

-

-    size_t size() const

-    {

-        return m_size;

-    }

-

-    void instance(Bezier & b, ConstVectorView const& raw_data) const

-    {

-        assert(b.size() == raw_data.size());

-        for (unsigned int i = 0; i < raw_data.size(); ++i)

-        {

-            b[i] = raw_data[i];

-        }

-    }

-

-  private:

-    size_t m_size;

-    size_t m_order;

-    std::vector<size_t> m_bc;

-};

-

-

-// this model generates Bezier curves

-template< unsigned int N >

-class LFMBezierCurve

-    : public LinearFittingModel< double, Point, BezierCurve<N> >

-{

-  public:

-    LFMBezierCurve( size_t _order )

-        : mob(_order)

-    {

-    }

-

-    void feed( VectorView & coeff, double t ) const

-    {

-        mob.feed(coeff, t);

-    }

-

-    size_t size() const

-    {

-        return mob.size();

-    }

-

-    void instance(BezierCurve<N> & bc, ConstMatrixView const& raw_data) const

-    {

-        Bezier bx(size()-1);

-        Bezier by(size()-1);

-        mob.instance(bx, raw_data.column_const_view(X));

-        mob.instance(by, raw_data.column_const_view(Y));

-        bc = BezierCurve<N>(bx, by);

-    }

-

-  private:

-    LFMBezier mob;

-};

-

-}  // end namespace NL

-}  // end namespace Geom

-

-

-#endif // _NL_FITTING_MODEL_H_

-

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

+/*
+ * Fitting Models for Geom Types
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NL_FITTING_MODEL_H_
+#define _NL_FITTING_MODEL_H_
+
+
+#include <2geom/d2.h>
+#include <2geom/sbasis.h>
+#include <2geom/bezier.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/poly.h>
+#include <2geom/ellipse.h>
+#include <2geom/utils.h>
+
+
+namespace Geom { namespace NL {
+
+
+/*
+ *   completely unknown models must inherit from this template class;
+ *   example: the model a*x^2 + b*x + c = 0 to be solved wrt a, b, c;
+ *   example: the model A(t) = known_sample_value_at(t) to be solved wrt
+ *       the coefficients of the curve A(t) expressed in S-Basis form;
+ *   parameter type: the type of x and t variable in the examples above;
+ *   value type:     the type of the known sample values (in the first example
+ *                   is constant )
+ *   instance type:  the type of the objects produced by using
+ *                   the fitting raw data solution
+ */
+template< typename ParameterType, typename ValueType, typename InstanceType >
+class LinearFittingModel
+{
+  public:
+    typedef ParameterType       parameter_type;
+    typedef ValueType           value_type;
+    typedef InstanceType        instance_type;
+
+    static const bool WITH_FIXED_TERMS = false;
+
+    /*
+     * a LinearFittingModel must implement the following methods:
+     *
+     * void feed( VectorView & vector,
+     *            parameter_type const& sample_parameter ) const;
+     *
+     * size_t size() const;
+     *
+     * void instance(instance_type &, raw_type const& raw_data) const;
+     *
+     */
+};
+
+
+/*
+ *   partially known models must inherit from this template class
+ *   example: the model a*x^2 + 2*x + c = 0 to be solved wrt a and c
+ */
+template< typename ParameterType, typename ValueType, typename InstanceType >
+class LinearFittingModelWithFixedTerms
+{
+  public:
+    typedef ParameterType       parameter_type;
+    typedef ValueType           value_type;
+    typedef InstanceType        instance_type;
+
+    static const bool WITH_FIXED_TERMS = true;
+
+    /*
+     * a LinearFittingModelWithFixedTerms must implement the following methods:
+     *
+     * void feed( VectorView & vector,
+     *            value_type & fixed_term,
+     *            parameter_type const& sample_parameter ) const;
+     *
+     * size_t size() const;
+     *
+     * void instance(instance_type &, raw_type const& raw_data) const;
+     *
+     */
+
+
+};
+
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of a polynomial
+// rapresented in standard power basis
+template< typename InstanceType >
+class LFMPowerBasis
+    : public LinearFittingModel<double, double, InstanceType>
+{
+  public:
+    LFMPowerBasis(size_t degree)
+        : m_size(degree + 1)
+    {
+    }
+
+    void feed( VectorView & coeff, double sample_parameter ) const
+    {
+        coeff[0] = 1;
+        double x_i = 1;
+        for (size_t i = 1; i < coeff.size(); ++i)
+        {
+          x_i *= sample_parameter;
+          coeff[i] = x_i;
+        }
+    }
+
+    size_t size() const
+    {
+        return m_size;
+    }
+
+  private:
+    size_t m_size;
+};
+
+
+// this model generates Geom::Poly objects
+class LFMPoly
+    : public LFMPowerBasis<Poly>
+{
+  public:
+    LFMPoly(size_t degree)
+        : LFMPowerBasis<Poly>(degree)
+    {
+    }
+
+    void instance(Poly & poly, ConstVectorView const& raw_data) const
+    {
+        poly.clear();
+        poly.resize(size());
+        for (size_t i = 0; i < raw_data.size(); ++i)
+        {
+            poly[i] =  raw_data[i];
+        }
+    }
+};
+
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of a polynomial
+// rapresented in standard power basis with leading term coefficient equal to 1
+template< typename InstanceType >
+class LFMNormalizedPowerBasis
+    : public LinearFittingModelWithFixedTerms<double, double, InstanceType>
+{
+  public:
+    LFMNormalizedPowerBasis(size_t _degree)
+        : m_model( _degree - 1)
+    {
+        assert(_degree > 0);
+    }
+
+
+    void feed( VectorView & coeff,
+               double & known_term,
+               double sample_parameter ) const
+    {
+        m_model.feed(coeff, sample_parameter);
+        known_term = coeff[m_model.size()-1] * sample_parameter;
+    }
+
+    size_t size() const
+    {
+        return m_model.size();
+    }
+
+  private:
+    LFMPowerBasis<InstanceType> m_model;
+};
+
+
+// incomplete model, it can be inherited to make up different kinds of
+// instance type; the raw data is a vector of coefficients of the equation
+// of an ellipse curve
+template< typename InstanceType >
+class LFMEllipseEquation
+    : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>
+{
+  public:
+    void feed( VectorView & coeff, double & fixed_term, Point const& p ) const
+    {
+        coeff[0] = p[X] * p[Y];
+        coeff[1] = p[Y] * p[Y];
+        coeff[2] = p[X];
+        coeff[3] = p[Y];
+        coeff[4] = 1;
+        fixed_term = p[X] * p[X];
+    }
+
+    size_t size() const
+    {
+        return 5;
+    }
+};
+
+
+// this model generates Ellipse curves
+class LFMEllipse
+    : public LFMEllipseEquation<Ellipse>
+{
+  public:
+    void instance(Ellipse & e, ConstVectorView const& coeff) const
+    {
+        e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
+    }
+};
+
+
+// this model generates SBasis objects
+class LFMSBasis
+    : public LinearFittingModel<double, double, SBasis>
+{
+  public:
+    LFMSBasis( size_t _order )
+        : m_size( 2*(_order+1) ),
+          m_order(_order)
+    {
+    }
+
+    void feed( VectorView & coeff, double t ) const
+    {
+        double u0 = 1-t;
+        double u1 = t;
+        double s = u0 * u1;
+        coeff[0] = u0;
+        coeff[1] = u1;
+        for (size_t i = 2; i < size(); i+=2)
+        {
+            u0 *= s;
+            u1 *= s;
+            coeff[i] = u0;
+            coeff[i+1] = u1;
+        }
+    }
+
+    size_t size() const
+    {
+        return m_size;
+    }
+
+    void instance(SBasis & sb, ConstVectorView const& raw_data) const
+    {
+        sb.clear();
+        sb.resize(m_order+1);
+        for (unsigned int i = 0, k = 0; i < raw_data.size(); i+=2, ++k)
+        {
+            sb[k][0] = raw_data[i];
+            sb[k][1] = raw_data[i+1];
+        }
+    }
+
+  private:
+    size_t m_size;
+    size_t m_order;
+};
+
+
+// this model generates D2<SBasis> objects
+class LFMD2SBasis
+    : public LinearFittingModel< double, Point, D2<SBasis> >
+{
+  public:
+    LFMD2SBasis( size_t _order )
+        : mosb(_order)
+    {
+    }
+
+    void feed( VectorView & coeff, double t ) const
+    {
+        mosb.feed(coeff, t);
+    }
+
+    size_t size() const
+    {
+        return mosb.size();
+    }
+
+    void instance(D2<SBasis> & d2sb, ConstMatrixView const& raw_data) const
+    {
+        mosb.instance(d2sb[X], raw_data.column_const_view(X));
+        mosb.instance(d2sb[Y], raw_data.column_const_view(Y));
+    }
+
+  private:
+    LFMSBasis mosb;
+};
+
+
+// this model generates Bezier objects
+class LFMBezier
+    : public LinearFittingModel<double, double, Bezier>
+{
+  public:
+    LFMBezier( size_t _order )
+        : m_size(_order + 1),
+          m_order(_order)
+    {
+        binomial_coefficients(m_bc, m_order);
+    }
+
+    void feed( VectorView & coeff, double t ) const
+    {
+        double s = 1;
+        for (size_t i = 0; i < size(); ++i)
+        {
+            coeff[i] = s * m_bc[i];
+            s *= t;
+        }
+        double u = 1-t;
+        s = 1;
+        for (size_t i = size()-1; i > 0; --i)
+        {
+            coeff[i] *= s;
+            s *= u;
+        }
+        coeff[0] *= s;
+    }
+
+    size_t size() const
+    {
+        return m_size;
+    }
+
+    void instance(Bezier & b, ConstVectorView const& raw_data) const
+    {
+        assert(b.size() == raw_data.size());
+        for (unsigned int i = 0; i < raw_data.size(); ++i)
+        {
+            b[i] = raw_data[i];
+        }
+    }
+
+  private:
+    size_t m_size;
+    size_t m_order;
+    std::vector<size_t> m_bc;
+};
+
+
+// this model generates Bezier curves
+template< unsigned int N >
+class LFMBezierCurve
+    : public LinearFittingModel< double, Point, BezierCurve<N> >
+{
+  public:
+    LFMBezierCurve( size_t _order )
+        : mob(_order)
+    {
+    }
+
+    void feed( VectorView & coeff, double t ) const
+    {
+        mob.feed(coeff, t);
+    }
+
+    size_t size() const
+    {
+        return mob.size();
+    }
+
+    void instance(BezierCurve<N> & bc, ConstMatrixView const& raw_data) const
+    {
+        Bezier bx(size()-1);
+        Bezier by(size()-1);
+        mob.instance(bx, raw_data.column_const_view(X));
+        mob.instance(by, raw_data.column_const_view(Y));
+        bc = BezierCurve<N>(bx, by);
+    }
+
+  private:
+    LFMBezier mob;
+};
+
+}  // end namespace NL
+}  // end namespace Geom
+
+
+#endif // _NL_FITTING_MODEL_H_
+
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/fitting-tool.h b/src/2geom/numeric/fitting-tool.h
index edacc663a..d589d86e3 100644
--- a/src/2geom/numeric/fitting-tool.h
+++ b/src/2geom/numeric/fitting-tool.h
@@ -1,532 +1,532 @@
-/*

- * Fitting Tools

- *

- * Authors:

- *      Marco Cecchetti <mrcekets at gmail.com>

- *

- * Copyright 2008  authors

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-#ifndef _NL_FITTING_TOOL_H_

-#define _NL_FITTING_TOOL_H_

-

-

-#include <2geom/numeric/vector.h>

-#include <2geom/numeric/matrix.h>

-

-#include <2geom/point.h>

-

-#include <vector>

-

-

-namespace Geom { namespace NL {

-

-namespace detail {

-

-

-template< typename ModelT>

-class lsf_base

-{

-  public:

-    typedef ModelT                                   model_type;

-    typedef typename model_type::parameter_type      parameter_type;

-    typedef typename model_type::value_type          value_type;

-

-    lsf_base( model_type const& _model, size_t forecasted_samples )

-        : m_model(_model),

-          m_total_samples(0),

-          m_matrix(forecasted_samples, m_model.size()),

-          m_psdinv_matrix(NULL)

-    {

-    }

-

-    // compute pseudo inverse

-    void update()

-    {

-        if (total_samples() == 0) return;

-        if (m_psdinv_matrix != NULL)

-        {

-            delete m_psdinv_matrix;

-        }

-        MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns());

-        m_psdinv_matrix = new Matrix( pseudo_inverse(mv) );

-        assert(m_psdinv_matrix != NULL);

-    }

-

-    size_t total_samples() const

-    {

-        return m_total_samples;

-    }

-

-    bool is_full() const

-    {

-        return (total_samples() == m_matrix.rows());

-    }

-

-    void clear()

-    {

-        m_total_samples = 0;

-    }

-

-    virtual

-    ~lsf_base()

-    {

-        if (m_psdinv_matrix != NULL)

-        {

-            delete m_psdinv_matrix;

-        }

-    }

-

-  protected:

-    const model_type &      m_model;

-    size_t                  m_total_samples;

-    Matrix                  m_matrix;

-    Matrix*                 m_psdinv_matrix;

-

-};  // end class lsf_base

-

-

-

-

-template< typename ModelT, typename ValueType = typename ModelT::value_type>

-class lsf_solution

-{

-};

-

-// a fitting process on samples with value of type double

-// produces a solution of type Vector

-template< typename ModelT>

-class lsf_solution<ModelT, double>

-    : public lsf_base<ModelT>

-{

-public:

-    typedef ModelT                                  model_type;

-    typedef typename model_type::parameter_type     parameter_type;

-    typedef typename model_type::value_type         value_type;

-    typedef Vector                                  solution_type;

-    typedef lsf_base<model_type>                    base_type;

-

-    using base_type::m_model;

-    using base_type::m_psdinv_matrix;

-    using base_type::total_samples;

-

-public:

-    lsf_solution<ModelT, double>( model_type const& _model,

-                                  size_t forecasted_samples )

-        : base_type(_model, forecasted_samples),

-          m_solution(_model.size())

-    {

-    }

-

-    template< typename VectorT >

-    solution_type& result(VectorT const& sample_values)

-    {

-        assert(sample_values.size() == total_samples());

-        ConstVectorView sv(sample_values);

-        m_solution = (*m_psdinv_matrix) * sv;

-        return m_solution;

-    }

-

-    // a comparison between old sample values and the new ones is performed

-    // in order to minimize computation

-    // prerequisite:

-    //     old_sample_values.size() == new_sample_values.size()

-    //     no update() call can be performed between two result invocations

-    template< typename VectorT >

-    solution_type& result( VectorT const& old_sample_values,

-                           VectorT const& new_sample_values )

-    {

-        assert(old_sample_values.size() == total_samples());

-        assert(new_sample_values.size() == total_samples());

-        Vector diff(total_samples());

-        for (size_t i = 0; i < diff.size(); ++i)

-        {

-            diff[i] = new_sample_values[i] - old_sample_values[i];

-        }

-        Vector column(m_model.size());

-        Vector delta(m_model.size(), 0.0);

-        for (size_t i = 0; i < diff.size(); ++i)

-        {

-            if (diff[i] != 0)

-            {

-                column = m_psdinv_matrix->column_view(i);

-                column.scale(diff[i]);

-                delta += column;

-            }

-        }

-        m_solution += delta;

-        return m_solution;

-    }

-

-    solution_type& result()

-    {

-        return m_solution;

-    }

-

-private:

-    solution_type m_solution;

-

-}; // end class lsf_solution<ModelT, double>

-

-

-// a fitting process on samples with value of type Point

-// produces a solution of type Matrix (with 2 columns)

-template< typename ModelT>

-class lsf_solution<ModelT, Point>

-    : public lsf_base<ModelT>

-{

-public:

-    typedef ModelT                                  model_type;

-    typedef typename model_type::parameter_type     parameter_type;

-    typedef typename model_type::value_type         value_type;

-    typedef Matrix                                  solution_type;

-    typedef lsf_base<model_type>                    base_type;

-

-    using base_type::m_model;

-    using base_type::m_psdinv_matrix;

-    using base_type::total_samples;

-

-public:

-    lsf_solution<ModelT, Point>( model_type const& _model,

-                                 size_t forecasted_samples )

-        : base_type(_model, forecasted_samples),

-          m_solution(_model.size(), 2)

-    {

-    }

-

-    solution_type& result(std::vector<Point> const& sample_values)

-    {

-        assert(sample_values.size() == total_samples());

-        Matrix svm(total_samples(), 2);

-        for (size_t i = 0; i < total_samples(); ++i)

-        {

-            svm(i, X) = sample_values[i][X];

-            svm(i, Y) = sample_values[i][Y];

-        }

-        m_solution = (*m_psdinv_matrix) * svm;

-        return m_solution;

-    }

-

-    // a comparison between old sample values and the new ones is performed

-    // in order to minimize computation

-    // prerequisite:

-    //     old_sample_values.size() == new_sample_values.size()

-    //     no update() call can to be performed between two result invocations

-    solution_type& result( std::vector<Point> const& old_sample_values,

-                           std::vector<Point> const& new_sample_values )

-    {

-        assert(old_sample_values.size() == total_samples());

-        assert(new_sample_values.size() == total_samples());

-        Matrix diff(total_samples(), 2);

-        for (size_t i = 0; i < total_samples(); ++i)

-        {

-            diff(i, X) =  new_sample_values[i][X] - old_sample_values[i][X];

-            diff(i, Y) =  new_sample_values[i][Y] - old_sample_values[i][Y];

-        }

-        Vector column(m_model.size());

-        Matrix delta(m_model.size(), 2, 0.0);

-        VectorView deltax = delta.column_view(X);

-        VectorView deltay = delta.column_view(Y);

-        for (size_t i = 0; i < total_samples(); ++i)

-        {

-            if (diff(i, X) != 0)

-            {

-                column = m_psdinv_matrix->column_view(i);

-                column.scale(diff(i, X));

-                deltax += column;

-            }

-            if (diff(i, Y) != 0)

-            {

-                column = m_psdinv_matrix->column_view(i);

-                column.scale(diff(i, Y));

-                deltay += column;

-            }

-        }

-        m_solution += delta;

-        return m_solution;

-    }

-

-    solution_type& result()

-    {

-        return m_solution;

-    }

-

-private:

-    solution_type m_solution;

-

-}; // end class lsf_solution<ModelT, Point>

-

-

-

-

-template< typename ModelT,

-          bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >

-class lsf_with_fixed_terms

-{

-};

-

-

-// fitting tool for completely unknown models

-template< typename ModelT>

-class lsf_with_fixed_terms<ModelT, false>

-    : public lsf_solution<ModelT>

-{

-  public:

-    typedef ModelT                                      model_type;

-    typedef typename model_type::parameter_type         parameter_type;

-    typedef typename model_type::value_type             value_type;

-    typedef lsf_solution<model_type>                    base_type;

-    typedef typename base_type::solution_type           solution_type;

-

-    using base_type::total_samples;

-    using base_type::is_full;

-    using base_type::m_matrix;

-    using base_type::m_total_samples;

-    using base_type::m_model;

-

-  public:

-      lsf_with_fixed_terms<ModelT, false>( model_type const& _model,

-                                           size_t forecasted_samples )

-        : base_type(_model, forecasted_samples)

-    {

-    }

-

-    void append(parameter_type const& sample_parameter)

-    {

-        assert(!is_full());

-        VectorView row = m_matrix.row_view(total_samples());

-        m_model.feed(row, sample_parameter);

-        ++m_total_samples;

-    }

-

-    void append_copy(size_t sample_index)

-    {

-        assert(!is_full());

-        assert(sample_index < total_samples());

-        VectorView dest_row = m_matrix.row_view(total_samples());

-        VectorView source_row = m_matrix.row_view(sample_index);

-        dest_row = source_row;

-        ++m_total_samples;

-    }

-

-}; // end class lsf_with_fixed_terms<ModelT, false>

-

-

-// fitting tool for partially known models

-template< typename ModelT>

-class lsf_with_fixed_terms<ModelT, true>

-    : public lsf_solution<ModelT>

-{

-  public:

-    typedef ModelT                                      model_type;

-    typedef typename model_type::parameter_type         parameter_type;

-    typedef typename model_type::value_type             value_type;

-    typedef lsf_solution<model_type>                    base_type;

-    typedef typename base_type::solution_type           solution_type;

-

-    using base_type::total_samples;

-    using base_type::is_full;

-    using base_type::m_matrix;

-    using base_type::m_total_samples;

-    using base_type::m_model;

-

-  public:

-    lsf_with_fixed_terms<ModelT, true>( model_type const& _model,

-                                        size_t forecasted_samples )

-        : base_type(_model, forecasted_samples),

-          m_vector(forecasted_samples),

-          m_vector_view(NULL)

-    {

-    }

-    void append(parameter_type const& sample_parameter)

-    {

-        assert(!is_full());

-        VectorView row = m_matrix.row_view(total_samples());

-        m_model.feed(row, m_vector[total_samples()], sample_parameter);

-        ++m_total_samples;

-    }

-

-    void append_copy(size_t sample_index)

-    {

-        assert(!is_full());

-        assert(sample_index < total_samples());

-        VectorView dest_row = m_matrix.row_view(total_samples());

-        VectorView source_row = m_matrix.row_view(sample_index);

-        dest_row = source_row;

-        m_vector[total_samples()] = m_vector[sample_index];

-        ++m_total_samples;

-    }

-

-    void update()

-    {

-        base_type::update();

-        if (total_samples() == 0) return;

-        if (m_vector_view != NULL)

-        {

-            delete m_vector_view;

-        }

-        m_vector_view = new VectorView(m_vector, base_type::total_samples());

-        assert(m_vector_view != NULL);

-    }

-

-    virtual

-    ~lsf_with_fixed_terms<model_type, true>()

-    {

-        if (m_vector_view != NULL)

-        {

-            delete m_vector_view;

-        }

-    }

-

-  protected:

-    Vector          m_vector;

-    VectorView*     m_vector_view;

-

-}; // end class lsf_with_fixed_terms<ModelT, true>

-

-

-} // end namespace detail

-

-

-

-

-template< typename ModelT,

-          typename ValueType = typename ModelT::value_type,

-          bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >

-class least_squeares_fitter

-{

-};

-

-

-template< typename ModelT, typename ValueType >

-class least_squeares_fitter<ModelT, ValueType, false>

-    : public detail::lsf_with_fixed_terms<ModelT>

-{

-  public:

-    typedef ModelT                                      model_type;

-    typedef detail::lsf_with_fixed_terms<model_type>    base_type;

-    typedef typename base_type::parameter_type          parameter_type;

-    typedef typename base_type::value_type              value_type;

-    typedef typename base_type::solution_type           solution_type;

-

-  public:

-    least_squeares_fitter<ModelT, ValueType, false>( model_type const& _model,

-                                                     size_t forecasted_samples )

-          : base_type(_model, forecasted_samples)

-  {

-  }

-}; // end class least_squeares_fitter<ModelT, ValueType, true>

-

-

-template< typename ModelT>

-class least_squeares_fitter<ModelT, double, true>

-    : public detail::lsf_with_fixed_terms<ModelT>

-{

-  public:

-    typedef ModelT                                      model_type;

-    typedef detail::lsf_with_fixed_terms<model_type>    base_type;

-    typedef typename base_type::parameter_type          parameter_type;

-    typedef typename base_type::value_type              value_type;

-    typedef typename base_type::solution_type           solution_type;

-

-    using base_type::m_vector_view;

-    using base_type::result;

-

-  public:

-    least_squeares_fitter<ModelT, double, true>( model_type const& _model,

-                                                 size_t forecasted_samples )

-        : base_type(_model, forecasted_samples)

-    {

-    }

-

-    template< typename VectorT >

-    solution_type& result(VectorT const& sample_values)

-    {

-        assert(sample_values.size() == m_vector_view->size());

-        Vector sv(sample_values.size());

-        for (size_t i = 0; i < sv.size(); ++i)

-            sv[i] = sample_values[i] - (*m_vector_view)[i];

-        return base_type::result(sv);

-    }

-

-}; // end class least_squeares_fitter<ModelT, double, true>

-

-

-template< typename ModelT>

-class least_squeares_fitter<ModelT, Point, true>

-    : public detail::lsf_with_fixed_terms<ModelT>

-{

-  public:

-    typedef ModelT                                      model_type;

-    typedef detail::lsf_with_fixed_terms<model_type>    base_type;

-    typedef typename base_type::parameter_type          parameter_type;

-    typedef typename base_type::value_type              value_type;

-    typedef typename base_type::solution_type           solution_type;

-

-    using base_type::m_vector_view;

-    using base_type::result;

-

-  public:

-    least_squeares_fitter<ModelT, Point, true>( model_type const& _model,

-                                                size_t forecasted_samples )

-        : base_type(_model, forecasted_samples)

-    {

-    }

-

-    solution_type& result(std::vector<Point> const& sample_values)

-    {

-        assert(sample_values.size() == m_vector_view->size());

-        NL::Matrix sv(sample_values.size(), 2);

-        for (size_t i = 0; i < sample_values.size(); ++i)

-        {

-            sv(i, X) = sample_values[i][X] - (*m_vector_view)[i];

-            sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i];

-        }

-        return base_type::result(sv);

-    }

-

-};  // end class least_squeares_fitter<ModelT, Point, true>

-

-

-}  // end namespace NL

-}  // end namespace Geom

-

-

-

-#endif  // _NL_FITTING_TOOL_H_

-

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

+/*
+ * Fitting Tools
+ *
+ * Authors:
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NL_FITTING_TOOL_H_
+#define _NL_FITTING_TOOL_H_
+
+
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/matrix.h>
+
+#include <2geom/point.h>
+
+#include <vector>
+
+
+namespace Geom { namespace NL {
+
+namespace detail {
+
+
+template< typename ModelT>
+class lsf_base
+{
+  public:
+    typedef ModelT                                   model_type;
+    typedef typename model_type::parameter_type      parameter_type;
+    typedef typename model_type::value_type          value_type;
+
+    lsf_base( model_type const& _model, size_t forecasted_samples )
+        : m_model(_model),
+          m_total_samples(0),
+          m_matrix(forecasted_samples, m_model.size()),
+          m_psdinv_matrix(NULL)
+    {
+    }
+
+    // compute pseudo inverse
+    void update()
+    {
+        if (total_samples() == 0) return;
+        if (m_psdinv_matrix != NULL)
+        {
+            delete m_psdinv_matrix;
+        }
+        MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns());
+        m_psdinv_matrix = new Matrix( pseudo_inverse(mv) );
+        assert(m_psdinv_matrix != NULL);
+    }
+
+    size_t total_samples() const
+    {
+        return m_total_samples;
+    }
+
+    bool is_full() const
+    {
+        return (total_samples() == m_matrix.rows());
+    }
+
+    void clear()
+    {
+        m_total_samples = 0;
+    }
+
+    virtual
+    ~lsf_base()
+    {
+        if (m_psdinv_matrix != NULL)
+        {
+            delete m_psdinv_matrix;
+        }
+    }
+
+  protected:
+    const model_type &      m_model;
+    size_t                  m_total_samples;
+    Matrix                  m_matrix;
+    Matrix*                 m_psdinv_matrix;
+
+};  // end class lsf_base
+
+
+
+
+template< typename ModelT, typename ValueType = typename ModelT::value_type>
+class lsf_solution
+{
+};
+
+// a fitting process on samples with value of type double
+// produces a solution of type Vector
+template< typename ModelT>
+class lsf_solution<ModelT, double>
+    : public lsf_base<ModelT>
+{
+public:
+    typedef ModelT                                  model_type;
+    typedef typename model_type::parameter_type     parameter_type;
+    typedef typename model_type::value_type         value_type;
+    typedef Vector                                  solution_type;
+    typedef lsf_base<model_type>                    base_type;
+
+    using base_type::m_model;
+    using base_type::m_psdinv_matrix;
+    using base_type::total_samples;
+
+public:
+    lsf_solution<ModelT, double>( model_type const& _model,
+                                  size_t forecasted_samples )
+        : base_type(_model, forecasted_samples),
+          m_solution(_model.size())
+    {
+    }
+
+    template< typename VectorT >
+    solution_type& result(VectorT const& sample_values)
+    {
+        assert(sample_values.size() == total_samples());
+        ConstVectorView sv(sample_values);
+        m_solution = (*m_psdinv_matrix) * sv;
+        return m_solution;
+    }
+
+    // a comparison between old sample values and the new ones is performed
+    // in order to minimize computation
+    // prerequisite:
+    //     old_sample_values.size() == new_sample_values.size()
+    //     no update() call can be performed between two result invocations
+    template< typename VectorT >
+    solution_type& result( VectorT const& old_sample_values,
+                           VectorT const& new_sample_values )
+    {
+        assert(old_sample_values.size() == total_samples());
+        assert(new_sample_values.size() == total_samples());
+        Vector diff(total_samples());
+        for (size_t i = 0; i < diff.size(); ++i)
+        {
+            diff[i] = new_sample_values[i] - old_sample_values[i];
+        }
+        Vector column(m_model.size());
+        Vector delta(m_model.size(), 0.0);
+        for (size_t i = 0; i < diff.size(); ++i)
+        {
+            if (diff[i] != 0)
+            {
+                column = m_psdinv_matrix->column_view(i);
+                column.scale(diff[i]);
+                delta += column;
+            }
+        }
+        m_solution += delta;
+        return m_solution;
+    }
+
+    solution_type& result()
+    {
+        return m_solution;
+    }
+
+private:
+    solution_type m_solution;
+
+}; // end class lsf_solution<ModelT, double>
+
+
+// a fitting process on samples with value of type Point
+// produces a solution of type Matrix (with 2 columns)
+template< typename ModelT>
+class lsf_solution<ModelT, Point>
+    : public lsf_base<ModelT>
+{
+public:
+    typedef ModelT                                  model_type;
+    typedef typename model_type::parameter_type     parameter_type;
+    typedef typename model_type::value_type         value_type;
+    typedef Matrix                                  solution_type;
+    typedef lsf_base<model_type>                    base_type;
+
+    using base_type::m_model;
+    using base_type::m_psdinv_matrix;
+    using base_type::total_samples;
+
+public:
+    lsf_solution<ModelT, Point>( model_type const& _model,
+                                 size_t forecasted_samples )
+        : base_type(_model, forecasted_samples),
+          m_solution(_model.size(), 2)
+    {
+    }
+
+    solution_type& result(std::vector<Point> const& sample_values)
+    {
+        assert(sample_values.size() == total_samples());
+        Matrix svm(total_samples(), 2);
+        for (size_t i = 0; i < total_samples(); ++i)
+        {
+            svm(i, X) = sample_values[i][X];
+            svm(i, Y) = sample_values[i][Y];
+        }
+        m_solution = (*m_psdinv_matrix) * svm;
+        return m_solution;
+    }
+
+    // a comparison between old sample values and the new ones is performed
+    // in order to minimize computation
+    // prerequisite:
+    //     old_sample_values.size() == new_sample_values.size()
+    //     no update() call can to be performed between two result invocations
+    solution_type& result( std::vector<Point> const& old_sample_values,
+                           std::vector<Point> const& new_sample_values )
+    {
+        assert(old_sample_values.size() == total_samples());
+        assert(new_sample_values.size() == total_samples());
+        Matrix diff(total_samples(), 2);
+        for (size_t i = 0; i < total_samples(); ++i)
+        {
+            diff(i, X) =  new_sample_values[i][X] - old_sample_values[i][X];
+            diff(i, Y) =  new_sample_values[i][Y] - old_sample_values[i][Y];
+        }
+        Vector column(m_model.size());
+        Matrix delta(m_model.size(), 2, 0.0);
+        VectorView deltax = delta.column_view(X);
+        VectorView deltay = delta.column_view(Y);
+        for (size_t i = 0; i < total_samples(); ++i)
+        {
+            if (diff(i, X) != 0)
+            {
+                column = m_psdinv_matrix->column_view(i);
+                column.scale(diff(i, X));
+                deltax += column;
+            }
+            if (diff(i, Y) != 0)
+            {
+                column = m_psdinv_matrix->column_view(i);
+                column.scale(diff(i, Y));
+                deltay += column;
+            }
+        }
+        m_solution += delta;
+        return m_solution;
+    }
+
+    solution_type& result()
+    {
+        return m_solution;
+    }
+
+private:
+    solution_type m_solution;
+
+}; // end class lsf_solution<ModelT, Point>
+
+
+
+
+template< typename ModelT,
+          bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
+class lsf_with_fixed_terms
+{
+};
+
+
+// fitting tool for completely unknown models
+template< typename ModelT>
+class lsf_with_fixed_terms<ModelT, false>
+    : public lsf_solution<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef typename model_type::parameter_type         parameter_type;
+    typedef typename model_type::value_type             value_type;
+    typedef lsf_solution<model_type>                    base_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::total_samples;
+    using base_type::is_full;
+    using base_type::m_matrix;
+    using base_type::m_total_samples;
+    using base_type::m_model;
+
+  public:
+      lsf_with_fixed_terms<ModelT, false>( model_type const& _model,
+                                           size_t forecasted_samples )
+        : base_type(_model, forecasted_samples)
+    {
+    }
+
+    void append(parameter_type const& sample_parameter)
+    {
+        assert(!is_full());
+        VectorView row = m_matrix.row_view(total_samples());
+        m_model.feed(row, sample_parameter);
+        ++m_total_samples;
+    }
+
+    void append_copy(size_t sample_index)
+    {
+        assert(!is_full());
+        assert(sample_index < total_samples());
+        VectorView dest_row = m_matrix.row_view(total_samples());
+        VectorView source_row = m_matrix.row_view(sample_index);
+        dest_row = source_row;
+        ++m_total_samples;
+    }
+
+}; // end class lsf_with_fixed_terms<ModelT, false>
+
+
+// fitting tool for partially known models
+template< typename ModelT>
+class lsf_with_fixed_terms<ModelT, true>
+    : public lsf_solution<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef typename model_type::parameter_type         parameter_type;
+    typedef typename model_type::value_type             value_type;
+    typedef lsf_solution<model_type>                    base_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::total_samples;
+    using base_type::is_full;
+    using base_type::m_matrix;
+    using base_type::m_total_samples;
+    using base_type::m_model;
+
+  public:
+    lsf_with_fixed_terms<ModelT, true>( model_type const& _model,
+                                        size_t forecasted_samples )
+        : base_type(_model, forecasted_samples),
+          m_vector(forecasted_samples),
+          m_vector_view(NULL)
+    {
+    }
+    void append(parameter_type const& sample_parameter)
+    {
+        assert(!is_full());
+        VectorView row = m_matrix.row_view(total_samples());
+        m_model.feed(row, m_vector[total_samples()], sample_parameter);
+        ++m_total_samples;
+    }
+
+    void append_copy(size_t sample_index)
+    {
+        assert(!is_full());
+        assert(sample_index < total_samples());
+        VectorView dest_row = m_matrix.row_view(total_samples());
+        VectorView source_row = m_matrix.row_view(sample_index);
+        dest_row = source_row;
+        m_vector[total_samples()] = m_vector[sample_index];
+        ++m_total_samples;
+    }
+
+    void update()
+    {
+        base_type::update();
+        if (total_samples() == 0) return;
+        if (m_vector_view != NULL)
+        {
+            delete m_vector_view;
+        }
+        m_vector_view = new VectorView(m_vector, base_type::total_samples());
+        assert(m_vector_view != NULL);
+    }
+
+    virtual
+    ~lsf_with_fixed_terms<model_type, true>()
+    {
+        if (m_vector_view != NULL)
+        {
+            delete m_vector_view;
+        }
+    }
+
+  protected:
+    Vector          m_vector;
+    VectorView*     m_vector_view;
+
+}; // end class lsf_with_fixed_terms<ModelT, true>
+
+
+} // end namespace detail
+
+
+
+
+template< typename ModelT,
+          typename ValueType = typename ModelT::value_type,
+          bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >
+class least_squeares_fitter
+{
+};
+
+
+template< typename ModelT, typename ValueType >
+class least_squeares_fitter<ModelT, ValueType, false>
+    : public detail::lsf_with_fixed_terms<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef detail::lsf_with_fixed_terms<model_type>    base_type;
+    typedef typename base_type::parameter_type          parameter_type;
+    typedef typename base_type::value_type              value_type;
+    typedef typename base_type::solution_type           solution_type;
+
+  public:
+    least_squeares_fitter<ModelT, ValueType, false>( model_type const& _model,
+                                                     size_t forecasted_samples )
+          : base_type(_model, forecasted_samples)
+  {
+  }
+}; // end class least_squeares_fitter<ModelT, ValueType, true>
+
+
+template< typename ModelT>
+class least_squeares_fitter<ModelT, double, true>
+    : public detail::lsf_with_fixed_terms<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef detail::lsf_with_fixed_terms<model_type>    base_type;
+    typedef typename base_type::parameter_type          parameter_type;
+    typedef typename base_type::value_type              value_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::m_vector_view;
+    using base_type::result;
+
+  public:
+    least_squeares_fitter<ModelT, double, true>( model_type const& _model,
+                                                 size_t forecasted_samples )
+        : base_type(_model, forecasted_samples)
+    {
+    }
+
+    template< typename VectorT >
+    solution_type& result(VectorT const& sample_values)
+    {
+        assert(sample_values.size() == m_vector_view->size());
+        Vector sv(sample_values.size());
+        for (size_t i = 0; i < sv.size(); ++i)
+            sv[i] = sample_values[i] - (*m_vector_view)[i];
+        return base_type::result(sv);
+    }
+
+}; // end class least_squeares_fitter<ModelT, double, true>
+
+
+template< typename ModelT>
+class least_squeares_fitter<ModelT, Point, true>
+    : public detail::lsf_with_fixed_terms<ModelT>
+{
+  public:
+    typedef ModelT                                      model_type;
+    typedef detail::lsf_with_fixed_terms<model_type>    base_type;
+    typedef typename base_type::parameter_type          parameter_type;
+    typedef typename base_type::value_type              value_type;
+    typedef typename base_type::solution_type           solution_type;
+
+    using base_type::m_vector_view;
+    using base_type::result;
+
+  public:
+    least_squeares_fitter<ModelT, Point, true>( model_type const& _model,
+                                                size_t forecasted_samples )
+        : base_type(_model, forecasted_samples)
+    {
+    }
+
+    solution_type& result(std::vector<Point> const& sample_values)
+    {
+        assert(sample_values.size() == m_vector_view->size());
+        NL::Matrix sv(sample_values.size(), 2);
+        for (size_t i = 0; i < sample_values.size(); ++i)
+        {
+            sv(i, X) = sample_values[i][X] - (*m_vector_view)[i];
+            sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i];
+        }
+        return base_type::result(sv);
+    }
+
+};  // end class least_squeares_fitter<ModelT, Point, true>
+
+
+}  // end namespace NL
+}  // end namespace Geom
+
+
+
+#endif  // _NL_FITTING_TOOL_H_
+
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/linear_system.h b/src/2geom/numeric/linear_system.h
index 5b516c9e6..dc2a1d7e0 100644
--- a/src/2geom/numeric/linear_system.h
+++ b/src/2geom/numeric/linear_system.h
@@ -1,138 +1,138 @@
-/*

- * LinearSystem class wraps some gsl routines for solving linear systems

- *

- * Authors:

- * 		Marco Cecchetti <mrcekets at gmail.com>

- * 

- * Copyright 2008  authors

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-#ifndef _NL_LINEAR_SYSTEM_H_

-#define _NL_LINEAR_SYSTEM_H_

-

-

-#include <cassert>

-

-#include <gsl/gsl_linalg.h>

-

-#include <2geom/numeric/matrix.h>

-#include <2geom/numeric/vector.h>

-

-

-namespace Geom { namespace NL {

-

-

-class LinearSystem

-{

-public:

-	LinearSystem(MatrixView & _matrix, VectorView & _vector)

-		: m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())

-	{

-	}

-	

-	LinearSystem(Matrix & _matrix, Vector & _vector)

-		: m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())

-	{

-	}

-	

-	const Vector & LU_solve()

-	{

-		assert( matrix().rows() == matrix().columns() 

-				&& matrix().rows() == vector().size() );

-		int s;

-		gsl_permutation * p = gsl_permutation_alloc(matrix().rows());

-		gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);

-		gsl_linalg_LU_solve( matrix().get_gsl_matrix(), 

-							 p, 

-				             vector().get_gsl_vector(), 

-				             m_solution.get_gsl_vector()

-				           );

-		gsl_permutation_free(p);

-		return solution();

-	}

-	

-	const Vector & SV_solve()

-	{

-		assert( matrix().rows() >= matrix().columns()

-				&& matrix().rows() == vector().size() );

-		

-		gsl_matrix* U = matrix().get_gsl_matrix();

-		gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());

-		gsl_vector* S = gsl_vector_alloc(matrix().columns());

-		gsl_vector* work = gsl_vector_alloc(matrix().columns());

-		

-		gsl_linalg_SV_decomp( U, V, S, work );

-		

-		gsl_vector* b = vector().get_gsl_vector();

-		gsl_vector* x = m_solution.get_gsl_vector();

-		

-		gsl_linalg_SV_solve( U, V, S, b, x);

-		

-		gsl_matrix_free(V);

-		gsl_vector_free(S);

-		gsl_vector_free(work);

-		

-		return solution();			  

-	}

-	

-	MatrixView & matrix()

-	{

-		return m_matrix;

-	}

-	

-	VectorView & vector()

-	{

-		return m_vector;

-	}

-	

-	const Vector & solution() const

-	{

-		return m_solution;

-	}

-	

-private:

-	MatrixView m_matrix;

-	VectorView m_vector;

-	Vector m_solution;

-};

-

-

-} } // end namespaces

-

-

-#endif /*_NL_LINEAR_SYSTEM_H_*/

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

+/*
+ * LinearSystem class wraps some gsl routines for solving linear systems
+ *
+ * Authors:
+ * 		Marco Cecchetti <mrcekets at gmail.com>
+ * 
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NL_LINEAR_SYSTEM_H_
+#define _NL_LINEAR_SYSTEM_H_
+
+
+#include <cassert>
+
+#include <gsl/gsl_linalg.h>
+
+#include <2geom/numeric/matrix.h>
+#include <2geom/numeric/vector.h>
+
+
+namespace Geom { namespace NL {
+
+
+class LinearSystem
+{
+public:
+	LinearSystem(MatrixView & _matrix, VectorView & _vector)
+		: m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
+	{
+	}
+	
+	LinearSystem(Matrix & _matrix, Vector & _vector)
+		: m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())
+	{
+	}
+	
+	const Vector & LU_solve()
+	{
+		assert( matrix().rows() == matrix().columns() 
+				&& matrix().rows() == vector().size() );
+		int s;
+		gsl_permutation * p = gsl_permutation_alloc(matrix().rows());
+		gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);
+		gsl_linalg_LU_solve( matrix().get_gsl_matrix(), 
+							 p, 
+				             vector().get_gsl_vector(), 
+				             m_solution.get_gsl_vector()
+				           );
+		gsl_permutation_free(p);
+		return solution();
+	}
+	
+	const Vector & SV_solve()
+	{
+		assert( matrix().rows() >= matrix().columns()
+				&& matrix().rows() == vector().size() );
+		
+		gsl_matrix* U = matrix().get_gsl_matrix();
+		gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());
+		gsl_vector* S = gsl_vector_alloc(matrix().columns());
+		gsl_vector* work = gsl_vector_alloc(matrix().columns());
+		
+		gsl_linalg_SV_decomp( U, V, S, work );
+		
+		gsl_vector* b = vector().get_gsl_vector();
+		gsl_vector* x = m_solution.get_gsl_vector();
+		
+		gsl_linalg_SV_solve( U, V, S, b, x);
+		
+		gsl_matrix_free(V);
+		gsl_vector_free(S);
+		gsl_vector_free(work);
+		
+		return solution();			  
+	}
+	
+	MatrixView & matrix()
+	{
+		return m_matrix;
+	}
+	
+	VectorView & vector()
+	{
+		return m_vector;
+	}
+	
+	const Vector & solution() const
+	{
+		return m_solution;
+	}
+	
+private:
+	MatrixView m_matrix;
+	VectorView m_vector;
+	Vector m_solution;
+};
+
+
+} } // end namespaces
+
+
+#endif /*_NL_LINEAR_SYSTEM_H_*/
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h
index 64557a6f1..156b6e9a2 100644
--- a/src/2geom/numeric/matrix.h
+++ b/src/2geom/numeric/matrix.h
@@ -1,555 +1,555 @@
-/*

- * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;

- * "views" mimic the semantic of C++ references: any operation performed

- * on a "view" is actually performed on the "viewed object"

- *

- * Authors:

- * 		Marco Cecchetti <mrcekets at gmail.com>

- *

- * Copyright 2008  authors

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-

-

-#ifndef _NL_MATRIX_H_

-#define _NL_MATRIX_H_

-

-#include <2geom/numeric/vector.h>

-

-#include <cassert>

-#include <utility>    // for std::pair

-#include <algorithm>  // for std::swap

-#include <sstream>

-#include <string>

-

-#include <gsl/gsl_matrix.h>

-#include <gsl/gsl_linalg.h>

-

-

-namespace Geom { namespace NL {

-

-namespace detail

-{

-

-class BaseMatrixImpl

-{

-  public:

-	virtual ~BaseMatrixImpl()

-	{

-	}

-

-	ConstVectorView row_const_view(size_t i) const

-	{

-		return ConstVectorView(gsl_matrix_const_row(m_matrix, i));

-	}

-

-	ConstVectorView column_const_view(size_t i) const

-	{

-		return ConstVectorView(gsl_matrix_const_column(m_matrix, i));

-	}

-

-	const double & operator() (size_t i, size_t j) const

-	{

-		return *gsl_matrix_const_ptr(m_matrix, i, j);

-	}

-

-	const gsl_matrix* get_gsl_matrix() const

-	{

-		return m_matrix;

-	}

-

-	bool is_zero() const

-	{

-		return gsl_matrix_isnull(m_matrix);

-	}

-

-	bool is_positive() const

-	{

-		return gsl_matrix_ispos(m_matrix);

-	}

-

-	bool is_negative() const

-	{

-		return gsl_matrix_isneg(m_matrix);

-	}

-

-	bool is_non_negative() const

-	{

-		for ( unsigned int i = 0; i < rows(); ++i )

-		{

-			for ( unsigned int j = 0; j < columns(); ++j )

-			{

-				if ( (*this)(i,j) < 0 ) return false;

-			}

-		}

-		return true;

-	}

-

-	double max() const

-	{

-		return gsl_matrix_max(m_matrix);

-	}

-

-	double min() const

-	{

-		return gsl_matrix_min(m_matrix);

-	}

-

-	std::pair<size_t, size_t>

-	max_index() const

-	{

-		std::pair<size_t, size_t> indices;

-		gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));

-		return indices;

-	}

-

-	std::pair<size_t, size_t>

-	min_index() const

-	{

-		std::pair<size_t, size_t> indices;

-		gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));

-		return indices;

-	}

-

-	size_t rows() const

-	{

-		return m_rows;

-	}

-

-	size_t columns() const

-	{

-		return m_columns;

-	}

-

-	std::string str() const;

-

-  protected:

-	size_t m_rows, m_columns;

-	gsl_matrix* m_matrix;

-

-};  // end class BaseMatrixImpl

-

-

-inline

-bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)

-{

-	if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;

-

-	for (size_t i = 0; i < m1.rows(); ++i)

-		for (size_t j = 0; j < m1.columns(); ++j)

-			if (m1(i,j) != m2(i,j)) return false;

-

-	return true;

-}

-

-template< class charT >

-inline

-std::basic_ostream<charT> &

-operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix)

-{

-	if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;

-

-	os << "[[" << _matrix(0,0);

-	for (size_t j = 1; j < _matrix.columns(); ++j)

-	{

-		os << ", " << _matrix(0,j);

-	}

-	os << "]";

-

-	for (size_t i = 1; i < _matrix.rows(); ++i)

-	{

-		os << ", [" << _matrix(i,0);

-		for (size_t j = 1; j < _matrix.columns(); ++j)

-		{

-			os << ", " << _matrix(i,j);

-		}

-		os << "]";

-	}

-	os << "]";

-	return os;

-}

-

-inline

-std::string BaseMatrixImpl::str() const

-{

-	std::ostringstream oss;

-	oss << (*this);

-	return oss.str();

-}

-

-

-class MatrixImpl : public BaseMatrixImpl

-{

-  public:

-

-	typedef BaseMatrixImpl base_type;

-

-	void set_all( double x )

-	{

-		gsl_matrix_set_all(m_matrix, x);

-	}

-

-	void set_identity()

-	{

-		gsl_matrix_set_identity(m_matrix);

-	}

-

-	using base_type::operator();

-

-	double & operator() (size_t i, size_t j)

-	{

-		return *gsl_matrix_ptr(m_matrix, i, j);

-	}

-

-	using base_type::get_gsl_matrix;

-

-	gsl_matrix* get_gsl_matrix()

-	{

-		return m_matrix;

-	}

-

-	VectorView row_view(size_t i)

-	{

-		return VectorView(gsl_matrix_row(m_matrix, i));

-	}

-

-	VectorView column_view(size_t i)

-	{

-		return VectorView(gsl_matrix_column(m_matrix, i));

-	}

-

-	void swap_rows(size_t i, size_t j)

-	{

-		 gsl_matrix_swap_rows(m_matrix, i, j);

-	}

-

-	void swap_columns(size_t i, size_t j)

-	{

-		gsl_matrix_swap_columns(m_matrix, i, j);

-	}

-

-	MatrixImpl & transpose()

-	{

-	    assert(columns() == rows());

-		gsl_matrix_transpose(m_matrix);

-		return (*this);

-	}

-

-	MatrixImpl & scale(double x)

-	{

-		gsl_matrix_scale(m_matrix, x);

-		return (*this);

-	}

-

-	MatrixImpl & translate(double x)

-	{

-		 gsl_matrix_add_constant(m_matrix, x);

-		 return (*this);

-	}

-

-	MatrixImpl & operator+=(base_type const& _matrix)

-	{

-		gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());

-		return (*this);

-	}

-

-	MatrixImpl & operator-=(base_type const& _matrix)

-	{

-		gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());

-		return (*this);

-	}

-

-}; // end class MatrixImpl

-

-}  // end namespace detail

-

-

-using detail::operator==;

-using detail::operator<<;

-

-

-

-

-class Matrix: public detail::MatrixImpl

-{

-  public:

-	typedef detail::MatrixImpl base_type;

-

-  public:

-	// the matrix is not inizialized

-	Matrix(size_t n1, size_t n2)

-	{

-		m_rows = n1;

-		m_columns = n2;

-		m_matrix = gsl_matrix_alloc(n1, n2);

-	}

-

-	Matrix(size_t n1, size_t n2, double x)

-	{

-		m_rows = n1;

-		m_columns = n2;

-		m_matrix = gsl_matrix_alloc(n1, n2);

-		gsl_matrix_set_all(m_matrix, x);

-	}

-

-	Matrix(Matrix const& _matrix)

-        : base_type()

-	{

-		m_rows = _matrix.rows();

-		m_columns = _matrix.columns();

-		m_matrix = gsl_matrix_alloc(rows(), columns());

-		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());

-	}

-

-	explicit

-	Matrix(base_type::base_type const& _matrix)

-	{

-		m_rows = _matrix.rows();

-		m_columns = _matrix.columns();

-		m_matrix = gsl_matrix_alloc(rows(), columns());

-		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());

-	}

-

-	Matrix & operator=(Matrix const& _matrix)

-	{

-		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );

-		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());

-		return *this;

-	}

-

-	Matrix & operator=(base_type::base_type const& _matrix)

-	{

-		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );

-		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());

-		return *this;

-	}

-

-	virtual ~Matrix()

-	{

-		gsl_matrix_free(m_matrix);

-	}

-

-	Matrix & transpose()

-	{

-		return static_cast<Matrix &>( base_type::transpose() );

-	}

-

-	Matrix & scale(double x)

-	{

-		return static_cast<Matrix &>( base_type::scale(x) );

-	}

-

-	Matrix & translate(double x)

-	{

-		return static_cast<Matrix &>( base_type::translate(x) );

-	}

-

-	Matrix & operator+=(base_type::base_type const& _matrix)

-	{

-		return static_cast<Matrix &>( base_type::operator+=(_matrix) );

-	}

-

-	Matrix & operator-=(base_type::base_type const& _matrix)

-	{

-		return static_cast<Matrix &>( base_type::operator-=(_matrix) );

-	}

-

-	friend

-	void swap(Matrix & m1, Matrix & m2);

-	friend

-	void swap_any(Matrix & m1, Matrix & m2);

-

-};  // end class Matrix

-

-

-// warning! this operation invalidates any view of the passed matrix objects

-inline

-void swap(Matrix & m1, Matrix & m2)

-{

-	assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );

-	std::swap(m1.m_matrix, m2.m_matrix);

-}

-

-inline

-void swap_any(Matrix & m1, Matrix & m2)

-{

-    std::swap(m1.m_matrix, m2.m_matrix);

-    std::swap(m1.m_rows, m2.m_rows);

-    std::swap(m1.m_columns, m2.m_columns);

-}

-

-

-

-class ConstMatrixView : public detail::BaseMatrixImpl

-{

-  public:

-	typedef detail::BaseMatrixImpl base_type;

-

-  public:

-	ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)

-		: m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )

-	{

-		m_rows = n1;

-		m_columns = n2;

-		m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );

-	}

-

-	ConstMatrixView(const ConstMatrixView & _matrix)

-		: base_type(),

-		  m_matrix_view(_matrix.m_matrix_view)

-	{

-		m_rows = _matrix.rows();

-		m_columns = _matrix.columns();

-		m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );

-	}

-

-	ConstMatrixView(const base_type & _matrix)

-		: m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))

-	{

-		m_rows = _matrix.rows();

-		m_columns = _matrix.columns();

-		m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );

-	}

-

-  private:

-	gsl_matrix_const_view m_matrix_view;

-

-};  // end class ConstMatrixView

-

-

-

-

-class MatrixView : public detail::MatrixImpl

-{

-  public:

-	typedef detail::MatrixImpl base_type;

-

-  public:

-	MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)

-	{

-		m_rows = n1;

-		m_columns = n2;

-		m_matrix_view

-			= gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);

-		m_matrix = &(m_matrix_view.matrix);

-	}

-

-	MatrixView(const MatrixView & _matrix)

-        : base_type()

-	{

-		m_rows = _matrix.rows();

-		m_columns = _matrix.columns();

-		m_matrix_view = _matrix.m_matrix_view;

-		m_matrix = &(m_matrix_view.matrix);

-	}

-

-	MatrixView(Matrix & _matrix)

-	{

-		m_rows = _matrix.rows();

-		m_columns = _matrix.columns();

-		m_matrix_view

-			= gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());

-		m_matrix = &(m_matrix_view.matrix);

-	}

-

-	MatrixView & operator=(MatrixView const& _matrix)

-	{

-		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );

-		gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);

-		return *this;

-	}

-

-	MatrixView & operator=(base_type::base_type const& _matrix)

-	{

-		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );

-		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());

-		return *this;

-	}

-

-	MatrixView & transpose()

-	{

-		return static_cast<MatrixView &>( base_type::transpose() );

-	}

-

-	MatrixView & scale(double x)

-	{

-		return static_cast<MatrixView &>( base_type::scale(x) );

-	}

-

-	MatrixView & translate(double x)

-	{

-		return static_cast<MatrixView &>( base_type::translate(x) );

-	}

-

-	MatrixView & operator+=(base_type::base_type const& _matrix)

-	{

-		return static_cast<MatrixView &>( base_type::operator+=(_matrix) );

-	}

-

-	MatrixView & operator-=(base_type::base_type const& _matrix)

-	{

-		return static_cast<MatrixView &>( base_type::operator-=(_matrix) );

-	}

-

-	friend

-	void swap_view(MatrixView & m1, MatrixView & m2);

-

-  private:

-	gsl_matrix_view m_matrix_view;

-

-};  // end class MatrixView

-

-

-inline

-void swap_view(MatrixView & m1, MatrixView & m2)

-{

-	assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );

-	std::swap(m1.m_matrix_view, m2.m_matrix_view);

-}

-

-Vector operator*( detail::BaseMatrixImpl const& A,

-                  detail::BaseVectorImpl const& v );

-

-Matrix operator*( detail::BaseMatrixImpl const& A,

-                  detail::BaseMatrixImpl const& B );

-

-Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);

-

-} } // end namespaces

-

-#endif /*_NL_MATRIX_H_*/

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

+/*
+ * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;
+ * "views" mimic the semantic of C++ references: any operation performed
+ * on a "view" is actually performed on the "viewed object"
+ *
+ * Authors:
+ * 		Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+
+
+#ifndef _NL_MATRIX_H_
+#define _NL_MATRIX_H_
+
+#include <2geom/numeric/vector.h>
+
+#include <cassert>
+#include <utility>    // for std::pair
+#include <algorithm>  // for std::swap
+#include <sstream>
+#include <string>
+
+#include <gsl/gsl_matrix.h>
+#include <gsl/gsl_linalg.h>
+
+
+namespace Geom { namespace NL {
+
+namespace detail
+{
+
+class BaseMatrixImpl
+{
+  public:
+	virtual ~BaseMatrixImpl()
+	{
+	}
+
+	ConstVectorView row_const_view(size_t i) const
+	{
+		return ConstVectorView(gsl_matrix_const_row(m_matrix, i));
+	}
+
+	ConstVectorView column_const_view(size_t i) const
+	{
+		return ConstVectorView(gsl_matrix_const_column(m_matrix, i));
+	}
+
+	const double & operator() (size_t i, size_t j) const
+	{
+		return *gsl_matrix_const_ptr(m_matrix, i, j);
+	}
+
+	const gsl_matrix* get_gsl_matrix() const
+	{
+		return m_matrix;
+	}
+
+	bool is_zero() const
+	{
+		return gsl_matrix_isnull(m_matrix);
+	}
+
+	bool is_positive() const
+	{
+		return gsl_matrix_ispos(m_matrix);
+	}
+
+	bool is_negative() const
+	{
+		return gsl_matrix_isneg(m_matrix);
+	}
+
+	bool is_non_negative() const
+	{
+		for ( unsigned int i = 0; i < rows(); ++i )
+		{
+			for ( unsigned int j = 0; j < columns(); ++j )
+			{
+				if ( (*this)(i,j) < 0 ) return false;
+			}
+		}
+		return true;
+	}
+
+	double max() const
+	{
+		return gsl_matrix_max(m_matrix);
+	}
+
+	double min() const
+	{
+		return gsl_matrix_min(m_matrix);
+	}
+
+	std::pair<size_t, size_t>
+	max_index() const
+	{
+		std::pair<size_t, size_t> indices;
+		gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));
+		return indices;
+	}
+
+	std::pair<size_t, size_t>
+	min_index() const
+	{
+		std::pair<size_t, size_t> indices;
+		gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));
+		return indices;
+	}
+
+	size_t rows() const
+	{
+		return m_rows;
+	}
+
+	size_t columns() const
+	{
+		return m_columns;
+	}
+
+	std::string str() const;
+
+  protected:
+	size_t m_rows, m_columns;
+	gsl_matrix* m_matrix;
+
+};  // end class BaseMatrixImpl
+
+
+inline
+bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)
+{
+	if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;
+
+	for (size_t i = 0; i < m1.rows(); ++i)
+		for (size_t j = 0; j < m1.columns(); ++j)
+			if (m1(i,j) != m2(i,j)) return false;
+
+	return true;
+}
+
+template< class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix)
+{
+	if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;
+
+	os << "[[" << _matrix(0,0);
+	for (size_t j = 1; j < _matrix.columns(); ++j)
+	{
+		os << ", " << _matrix(0,j);
+	}
+	os << "]";
+
+	for (size_t i = 1; i < _matrix.rows(); ++i)
+	{
+		os << ", [" << _matrix(i,0);
+		for (size_t j = 1; j < _matrix.columns(); ++j)
+		{
+			os << ", " << _matrix(i,j);
+		}
+		os << "]";
+	}
+	os << "]";
+	return os;
+}
+
+inline
+std::string BaseMatrixImpl::str() const
+{
+	std::ostringstream oss;
+	oss << (*this);
+	return oss.str();
+}
+
+
+class MatrixImpl : public BaseMatrixImpl
+{
+  public:
+
+	typedef BaseMatrixImpl base_type;
+
+	void set_all( double x )
+	{
+		gsl_matrix_set_all(m_matrix, x);
+	}
+
+	void set_identity()
+	{
+		gsl_matrix_set_identity(m_matrix);
+	}
+
+	using base_type::operator();
+
+	double & operator() (size_t i, size_t j)
+	{
+		return *gsl_matrix_ptr(m_matrix, i, j);
+	}
+
+	using base_type::get_gsl_matrix;
+
+	gsl_matrix* get_gsl_matrix()
+	{
+		return m_matrix;
+	}
+
+	VectorView row_view(size_t i)
+	{
+		return VectorView(gsl_matrix_row(m_matrix, i));
+	}
+
+	VectorView column_view(size_t i)
+	{
+		return VectorView(gsl_matrix_column(m_matrix, i));
+	}
+
+	void swap_rows(size_t i, size_t j)
+	{
+		 gsl_matrix_swap_rows(m_matrix, i, j);
+	}
+
+	void swap_columns(size_t i, size_t j)
+	{
+		gsl_matrix_swap_columns(m_matrix, i, j);
+	}
+
+	MatrixImpl & transpose()
+	{
+	    assert(columns() == rows());
+		gsl_matrix_transpose(m_matrix);
+		return (*this);
+	}
+
+	MatrixImpl & scale(double x)
+	{
+		gsl_matrix_scale(m_matrix, x);
+		return (*this);
+	}
+
+	MatrixImpl & translate(double x)
+	{
+		 gsl_matrix_add_constant(m_matrix, x);
+		 return (*this);
+	}
+
+	MatrixImpl & operator+=(base_type const& _matrix)
+	{
+		gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());
+		return (*this);
+	}
+
+	MatrixImpl & operator-=(base_type const& _matrix)
+	{
+		gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());
+		return (*this);
+	}
+
+}; // end class MatrixImpl
+
+}  // end namespace detail
+
+
+using detail::operator==;
+using detail::operator<<;
+
+
+
+
+class Matrix: public detail::MatrixImpl
+{
+  public:
+	typedef detail::MatrixImpl base_type;
+
+  public:
+	// the matrix is not inizialized
+	Matrix(size_t n1, size_t n2)
+	{
+		m_rows = n1;
+		m_columns = n2;
+		m_matrix = gsl_matrix_alloc(n1, n2);
+	}
+
+	Matrix(size_t n1, size_t n2, double x)
+	{
+		m_rows = n1;
+		m_columns = n2;
+		m_matrix = gsl_matrix_alloc(n1, n2);
+		gsl_matrix_set_all(m_matrix, x);
+	}
+
+	Matrix(Matrix const& _matrix)
+        : base_type()
+	{
+		m_rows = _matrix.rows();
+		m_columns = _matrix.columns();
+		m_matrix = gsl_matrix_alloc(rows(), columns());
+		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+	}
+
+	explicit
+	Matrix(base_type::base_type const& _matrix)
+	{
+		m_rows = _matrix.rows();
+		m_columns = _matrix.columns();
+		m_matrix = gsl_matrix_alloc(rows(), columns());
+		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+	}
+
+	Matrix & operator=(Matrix const& _matrix)
+	{
+		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );
+		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+		return *this;
+	}
+
+	Matrix & operator=(base_type::base_type const& _matrix)
+	{
+		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );
+		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+		return *this;
+	}
+
+	virtual ~Matrix()
+	{
+		gsl_matrix_free(m_matrix);
+	}
+
+	Matrix & transpose()
+	{
+		return static_cast<Matrix &>( base_type::transpose() );
+	}
+
+	Matrix & scale(double x)
+	{
+		return static_cast<Matrix &>( base_type::scale(x) );
+	}
+
+	Matrix & translate(double x)
+	{
+		return static_cast<Matrix &>( base_type::translate(x) );
+	}
+
+	Matrix & operator+=(base_type::base_type const& _matrix)
+	{
+		return static_cast<Matrix &>( base_type::operator+=(_matrix) );
+	}
+
+	Matrix & operator-=(base_type::base_type const& _matrix)
+	{
+		return static_cast<Matrix &>( base_type::operator-=(_matrix) );
+	}
+
+	friend
+	void swap(Matrix & m1, Matrix & m2);
+	friend
+	void swap_any(Matrix & m1, Matrix & m2);
+
+};  // end class Matrix
+
+
+// warning! this operation invalidates any view of the passed matrix objects
+inline
+void swap(Matrix & m1, Matrix & m2)
+{
+	assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );
+	std::swap(m1.m_matrix, m2.m_matrix);
+}
+
+inline
+void swap_any(Matrix & m1, Matrix & m2)
+{
+    std::swap(m1.m_matrix, m2.m_matrix);
+    std::swap(m1.m_rows, m2.m_rows);
+    std::swap(m1.m_columns, m2.m_columns);
+}
+
+
+
+class ConstMatrixView : public detail::BaseMatrixImpl
+{
+  public:
+	typedef detail::BaseMatrixImpl base_type;
+
+  public:
+	ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)
+		: m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )
+	{
+		m_rows = n1;
+		m_columns = n2;
+		m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
+	}
+
+	ConstMatrixView(const ConstMatrixView & _matrix)
+		: base_type(),
+		  m_matrix_view(_matrix.m_matrix_view)
+	{
+		m_rows = _matrix.rows();
+		m_columns = _matrix.columns();
+		m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
+	}
+
+	ConstMatrixView(const base_type & _matrix)
+		: m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))
+	{
+		m_rows = _matrix.rows();
+		m_columns = _matrix.columns();
+		m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );
+	}
+
+  private:
+	gsl_matrix_const_view m_matrix_view;
+
+};  // end class ConstMatrixView
+
+
+
+
+class MatrixView : public detail::MatrixImpl
+{
+  public:
+	typedef detail::MatrixImpl base_type;
+
+  public:
+	MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)
+	{
+		m_rows = n1;
+		m_columns = n2;
+		m_matrix_view
+			= gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);
+		m_matrix = &(m_matrix_view.matrix);
+	}
+
+	MatrixView(const MatrixView & _matrix)
+        : base_type()
+	{
+		m_rows = _matrix.rows();
+		m_columns = _matrix.columns();
+		m_matrix_view = _matrix.m_matrix_view;
+		m_matrix = &(m_matrix_view.matrix);
+	}
+
+	MatrixView(Matrix & _matrix)
+	{
+		m_rows = _matrix.rows();
+		m_columns = _matrix.columns();
+		m_matrix_view
+			= gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());
+		m_matrix = &(m_matrix_view.matrix);
+	}
+
+	MatrixView & operator=(MatrixView const& _matrix)
+	{
+		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );
+		gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);
+		return *this;
+	}
+
+	MatrixView & operator=(base_type::base_type const& _matrix)
+	{
+		assert( rows() == _matrix.rows() && columns() ==  _matrix.columns() );
+		gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());
+		return *this;
+	}
+
+	MatrixView & transpose()
+	{
+		return static_cast<MatrixView &>( base_type::transpose() );
+	}
+
+	MatrixView & scale(double x)
+	{
+		return static_cast<MatrixView &>( base_type::scale(x) );
+	}
+
+	MatrixView & translate(double x)
+	{
+		return static_cast<MatrixView &>( base_type::translate(x) );
+	}
+
+	MatrixView & operator+=(base_type::base_type const& _matrix)
+	{
+		return static_cast<MatrixView &>( base_type::operator+=(_matrix) );
+	}
+
+	MatrixView & operator-=(base_type::base_type const& _matrix)
+	{
+		return static_cast<MatrixView &>( base_type::operator-=(_matrix) );
+	}
+
+	friend
+	void swap_view(MatrixView & m1, MatrixView & m2);
+
+  private:
+	gsl_matrix_view m_matrix_view;
+
+};  // end class MatrixView
+
+
+inline
+void swap_view(MatrixView & m1, MatrixView & m2)
+{
+	assert( m1.rows() == m2.rows() && m1.columns() ==  m2.columns() );
+	std::swap(m1.m_matrix_view, m2.m_matrix_view);
+}
+
+Vector operator*( detail::BaseMatrixImpl const& A,
+                  detail::BaseVectorImpl const& v );
+
+Matrix operator*( detail::BaseMatrixImpl const& A,
+                  detail::BaseMatrixImpl const& B );
+
+Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);
+
+} } // end namespaces
+
+#endif /*_NL_MATRIX_H_*/
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/numeric/vector.h b/src/2geom/numeric/vector.h
index 43a39a1ac..3e53405f4 100644
--- a/src/2geom/numeric/vector.h
+++ b/src/2geom/numeric/vector.h
@@ -1,565 +1,565 @@
-/*

- * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines;

- * "views" mimic the semantic of C++ references: any operation performed

- * on a "view" is actually performed on the "viewed object"

- *

- * Authors:

- * 		Marco Cecchetti <mrcekets at gmail.com>

- *

- * Copyright 2008  authors

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-

-

-#ifndef _NL_VECTOR_H_

-#define _NL_VECTOR_H_

-

-#include <cassert>

-#include <algorithm> // for std::swap

-#include <vector>

-#include <sstream>

-#include <string>

-

-

-#include <gsl/gsl_vector.h>

-#include <gsl/gsl_blas.h>

-

-

-namespace Geom { namespace NL {

-

-namespace detail

-{

-

-class BaseVectorImpl

-{

-  public:

-	double const& operator[](size_t i) const

-	{

-		return *gsl_vector_const_ptr(m_vector, i);

-	}

-

-	const gsl_vector* get_gsl_vector() const

-	{

-		return m_vector;

-	}

-	bool is_zero() const

-	{

-		return gsl_vector_isnull(m_vector);

-	}

-

-	bool is_positive() const

-	{

-		return gsl_vector_ispos(m_vector);

-	}

-

-	bool is_negative() const

-	{

-		return gsl_vector_isneg(m_vector);

-	}

-

-	bool is_non_negative() const

-	{

-		for ( size_t i = 0; i < size(); ++i )

-		{

-			if ( (*this)[i] < 0 ) return false;

-		}

-		return true;

-	}

-

-	double max() const

-	{

-		return gsl_vector_max(m_vector);

-	}

-

-	double min() const

-	{

-		return gsl_vector_min(m_vector);

-	}

-

-	size_t max_index() const

-	{

-		return gsl_vector_max_index(m_vector);

-	}

-

-	size_t min_index() const

-	{

-		return gsl_vector_min_index(m_vector);

-	}

-

-	size_t size() const

-	{

-		return m_size;

-	}

-

-	std::string str() const;

-

-	virtual ~BaseVectorImpl()

-	{

-	}

-

-  protected:

-	size_t m_size;

-	gsl_vector* m_vector;

-

-};  // end class BaseVectorImpl

-

-

-inline

-bool operator== (BaseVectorImpl const& v1, BaseVectorImpl const& v2)

-{

-	if (v1.size() != v2.size())	return false;

-

-	for (size_t i = 0; i < v1.size(); ++i)

-	{

-		if (v1[i] != v2[i])  return false;

-	}

-	return true;

-}

-

-template< class charT >

-inline

-std::basic_ostream<charT> &

-operator<< (std::basic_ostream<charT> & os, const BaseVectorImpl & _vector)

-{

-	if (_vector.size() == 0 ) return os;

-	os << "[" << _vector[0];

-	for (unsigned int i = 1; i < _vector.size(); ++i)

-	{

-		os << ", " << _vector[i];

-	}

-	os << "]";

-	return os;

-}

-

-inline

-std::string BaseVectorImpl::str() const

-{

-	std::ostringstream oss;

-	oss << (*this);

-	return oss.str();

-}

-

-inline

-double dot(BaseVectorImpl const& v1, BaseVectorImpl const& v2)

-{

-    double result;

-    gsl_blas_ddot(v1.get_gsl_vector(), v2.get_gsl_vector(), &result);

-    return result;

-}

-

-

-class VectorImpl : public BaseVectorImpl

-{

-  public:

-	typedef BaseVectorImpl base_type;

-

-  public:

-	void set_all(double x)

-	{

-		gsl_vector_set_all(m_vector, x);

-	}

-

-	void set_basis(size_t i)

-	{

-		gsl_vector_set_basis(m_vector, i);

-	}

-

-	using base_type::operator[];

-

-	double & operator[](size_t i)

-	{

-		return *gsl_vector_ptr(m_vector, i);

-	}

-

-	using base_type::get_gsl_vector;

-

-	gsl_vector* get_gsl_vector()

-	{

-		return m_vector;

-	}

-

-	void swap_elements(size_t i, size_t j)

-	{

-		gsl_vector_swap_elements(m_vector, i, j);

-	}

-

-	void reverse()

-	{

-		gsl_vector_reverse(m_vector);

-	}

-

-	VectorImpl & scale(double x)

-	{

-		gsl_vector_scale(m_vector, x);

-		return (*this);

-	}

-

-	VectorImpl & translate(double x)

-	{

-		gsl_vector_add_constant(m_vector, x);

-		return (*this);

-	}

-

-	VectorImpl & operator+=(base_type const& _vector)

-	{

-		gsl_vector_add(m_vector, _vector.get_gsl_vector());

-		return (*this);

-	}

-

-	VectorImpl & operator-=(base_type const& _vector)

-	{

-		gsl_vector_sub(m_vector, _vector.get_gsl_vector());

-		return (*this);

-	}

-

-};  // end class VectorImpl

-

-}  // end namespace detail

-

-

-using detail::operator==;

-using detail::operator<<;

-

-class Vector : public detail::VectorImpl

-{

-  public:

-	typedef detail::VectorImpl base_type;

-

-  public:

-	Vector(size_t n)

-	{

-		m_size = n;

-		m_vector = gsl_vector_alloc(n);

-	}

-

-	Vector(size_t n, double x)

-	{

-		m_size = n;

-		m_vector = gsl_vector_alloc(n);

-		gsl_vector_set_all(m_vector, x);

-	}

-

-	// create a vector with n elements all set to zero

-	// but the i-th that is set to 1

-	Vector(size_t n, size_t i)

-	{

-		m_size = n;

-		m_vector = gsl_vector_alloc(n);

-		gsl_vector_set_basis(m_vector, i);

-	}

-

-	Vector(Vector const& _vector)

-        : base_type()

-	{

-		m_size = _vector.size();

-		m_vector = gsl_vector_alloc(size());

-		gsl_vector_memcpy(m_vector, _vector.m_vector);

-	}

-

-	explicit

-	Vector(base_type::base_type const& _vector)

-	{

-		m_size = _vector.size();

-		m_vector = gsl_vector_alloc(size());

-		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());

-	}

-

-	virtual ~Vector()

-	{

-		gsl_vector_free(m_vector);

-	}

-

-

-	Vector & operator=(Vector const& _vector)

-	{

-		assert( size() == _vector.size() );

-		gsl_vector_memcpy(m_vector, _vector.m_vector);

-		return (*this);

-	}

-

-	Vector & operator=(base_type::base_type const& _vector)

-	{

-		assert( size() == _vector.size() );

-		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());

-		return (*this);

-	}

-

-	Vector & scale(double x)

-	{

-		return static_cast<Vector&>( base_type::scale(x) );

-	}

-

-	Vector & translate(double x)

-	{

-		return static_cast<Vector&>( base_type::translate(x) );

-	}

-

-	Vector & operator+=(base_type::base_type const& _vector)

-	{

-		return static_cast<Vector&>( base_type::operator+=(_vector) );

-	}

-

-	Vector & operator-=(base_type::base_type const& _vector)

-	{

-		return static_cast<Vector&>( base_type::operator-=(_vector) );

-	}

-

-	friend

-	void swap(Vector & v1, Vector & v2);

-	friend

-	void swap_any(Vector & v1, Vector & v2);

-

-}; // end class Vector

-

-

-// warning! these operations invalidate any view of the passed vector objects

-inline

-void swap(Vector & v1, Vector & v2)

-{

-	assert( v1.size() == v2.size() );

-	std::swap(v1.m_vector, v2.m_vector);

-}

-

-inline

-void swap_any(Vector & v1, Vector & v2)

-{

-    std::swap(v1.m_vector, v2.m_vector);

-    std::swap(v1.m_size, v2.m_size);

-}

-

-

-class ConstVectorView : public detail::BaseVectorImpl

-{

-  public:

-	typedef detail::BaseVectorImpl base_type;

-

-  public:

-	ConstVectorView(const base_type & _vector, size_t n, size_t offset = 0)

-		: m_vector_view( gsl_vector_const_subvector(_vector.get_gsl_vector(), offset, n) )

-	{

-		m_size = n;

-		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

-	}

-

-	ConstVectorView(const base_type & _vector, size_t n, size_t offset , size_t stride)

-		: m_vector_view( gsl_vector_const_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n) )

-	{

-		m_size = n;

-		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

-	}

-

-    ConstVectorView(const double* _vector, size_t n, size_t offset = 0)

-        : m_vector_view( gsl_vector_const_view_array(_vector + offset, n) )

-    {

-        m_size = n;

-        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

-    }

-

-    ConstVectorView(const double* _vector, size_t n, size_t offset, size_t stride)

-        : m_vector_view( gsl_vector_const_view_array_with_stride(_vector + offset, stride, n) )

-    {

-        m_size = n;

-        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

-    }

-

-	explicit

-	ConstVectorView(gsl_vector_const_view  _gsl_vector_view)

-		: m_vector_view(_gsl_vector_view)

-	{

-		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

-		m_size = m_vector->size;

-	}

-

-    explicit

-    ConstVectorView(const std::vector<double>&  _vector)

-        : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) )

-    {

-        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

-        m_size = _vector.size();

-    }

-

-	ConstVectorView(const ConstVectorView & _vector)

-		: base_type(),

-		  m_vector_view(_vector.m_vector_view)

-	{

-		m_size = _vector.size();

-		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

-	}

-

-	ConstVectorView(const base_type & _vector)

-		: m_vector_view(gsl_vector_const_subvector(_vector.get_gsl_vector(), 0, _vector.size()))

-	{

-		m_size = _vector.size();

-		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );

-	}

-

-  private:

-	gsl_vector_const_view m_vector_view;

-

-}; // end class  ConstVectorView

-

-

-

-

-class VectorView : public detail::VectorImpl

-{

-  public:

-	typedef detail::VectorImpl base_type;

-

-  public:

-	VectorView(base_type & _vector, size_t n, size_t offset = 0, size_t stride = 1)

-	{

-		m_size = n;

-		if (stride == 1)

-		{

-			m_vector_view

-				= gsl_vector_subvector(_vector.get_gsl_vector(), offset, n);

-			m_vector = &(m_vector_view.vector);

-		}

-		else

-		{

-			m_vector_view

-				= gsl_vector_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n);

-			m_vector = &(m_vector_view.vector);

-		}

-	}

-

-    VectorView(double* _vector, size_t n, size_t offset = 0, size_t stride = 1)

-    {

-        m_size = n;

-        if (stride == 1)

-        {

-            m_vector_view

-                = gsl_vector_view_array(_vector + offset, n);

-            m_vector = &(m_vector_view.vector);

-        }

-        else

-        {

-            m_vector_view

-                = gsl_vector_view_array_with_stride(_vector + offset, stride, n);

-            m_vector = &(m_vector_view.vector);

-        }

-

-    }

-

-	VectorView(const VectorView & _vector)

-        : base_type()

-	{

-		m_size = _vector.size();

-		m_vector_view = _vector.m_vector_view;

-		m_vector = &(m_vector_view.vector);

-	}

-

-	VectorView(Vector & _vector)

-	{

-		m_size = _vector.size();

-		m_vector_view = gsl_vector_subvector(_vector.get_gsl_vector(), 0, size());

-		m_vector = &(m_vector_view.vector);

-	}

-

-	explicit

-	VectorView(gsl_vector_view _gsl_vector_view)

-		: m_vector_view(_gsl_vector_view)

-	{

-		m_vector = &(m_vector_view.vector);

-		m_size = m_vector->size;

-	}

-

-	explicit

-	VectorView(std::vector<double> & _vector)

-	{

-	    m_size = _vector.size();

-	    m_vector_view = gsl_vector_view_array(&(_vector[0]), _vector.size());

-	    m_vector = &(m_vector_view.vector);

-	}

-

-	VectorView & operator=(VectorView const& _vector)

-	{

-		assert( size() == _vector.size() );

-		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());

-		return (*this);

-	}

-

-	VectorView & operator=(base_type::base_type const& _vector)

-	{

-		assert( size() == _vector.size() );

-		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());

-		return (*this);

-	}

-

-	VectorView & scale(double x)

-	{

-		return static_cast<VectorView&>( base_type::scale(x) );

-	}

-

-	VectorView & translate(double x)

-	{

-		return static_cast<VectorView&>( base_type::translate(x) );

-	}

-

-	VectorView & operator+=(base_type::base_type const& _vector)

-	{

-		return static_cast<VectorView&>( base_type::operator+=(_vector) );

-	}

-

-	VectorView & operator-=(base_type::base_type const& _vector)

-	{

-		return static_cast<VectorView&>( base_type::operator-=(_vector) );

-	}

-

-	friend

-	void swap_view(VectorView & v1, VectorView & v2);

-

-  private:

-	gsl_vector_view m_vector_view;

-

-}; // end class VectorView

-

-

-inline

-void swap_view(VectorView & v1, VectorView & v2)

-{

-	assert( v1.size() == v2.size() );

-	std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too

-}

-

-

-} } // end namespaces

-

-

-#endif /*_NL_VECTOR_H_*/

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

+/*
+ * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines;
+ * "views" mimic the semantic of C++ references: any operation performed
+ * on a "view" is actually performed on the "viewed object"
+ *
+ * Authors:
+ * 		Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+
+
+#ifndef _NL_VECTOR_H_
+#define _NL_VECTOR_H_
+
+#include <cassert>
+#include <algorithm> // for std::swap
+#include <vector>
+#include <sstream>
+#include <string>
+
+
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+
+
+namespace Geom { namespace NL {
+
+namespace detail
+{
+
+class BaseVectorImpl
+{
+  public:
+	double const& operator[](size_t i) const
+	{
+		return *gsl_vector_const_ptr(m_vector, i);
+	}
+
+	const gsl_vector* get_gsl_vector() const
+	{
+		return m_vector;
+	}
+	bool is_zero() const
+	{
+		return gsl_vector_isnull(m_vector);
+	}
+
+	bool is_positive() const
+	{
+		return gsl_vector_ispos(m_vector);
+	}
+
+	bool is_negative() const
+	{
+		return gsl_vector_isneg(m_vector);
+	}
+
+	bool is_non_negative() const
+	{
+		for ( size_t i = 0; i < size(); ++i )
+		{
+			if ( (*this)[i] < 0 ) return false;
+		}
+		return true;
+	}
+
+	double max() const
+	{
+		return gsl_vector_max(m_vector);
+	}
+
+	double min() const
+	{
+		return gsl_vector_min(m_vector);
+	}
+
+	size_t max_index() const
+	{
+		return gsl_vector_max_index(m_vector);
+	}
+
+	size_t min_index() const
+	{
+		return gsl_vector_min_index(m_vector);
+	}
+
+	size_t size() const
+	{
+		return m_size;
+	}
+
+	std::string str() const;
+
+	virtual ~BaseVectorImpl()
+	{
+	}
+
+  protected:
+	size_t m_size;
+	gsl_vector* m_vector;
+
+};  // end class BaseVectorImpl
+
+
+inline
+bool operator== (BaseVectorImpl const& v1, BaseVectorImpl const& v2)
+{
+	if (v1.size() != v2.size())	return false;
+
+	for (size_t i = 0; i < v1.size(); ++i)
+	{
+		if (v1[i] != v2[i])  return false;
+	}
+	return true;
+}
+
+template< class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const BaseVectorImpl & _vector)
+{
+	if (_vector.size() == 0 ) return os;
+	os << "[" << _vector[0];
+	for (unsigned int i = 1; i < _vector.size(); ++i)
+	{
+		os << ", " << _vector[i];
+	}
+	os << "]";
+	return os;
+}
+
+inline
+std::string BaseVectorImpl::str() const
+{
+	std::ostringstream oss;
+	oss << (*this);
+	return oss.str();
+}
+
+inline
+double dot(BaseVectorImpl const& v1, BaseVectorImpl const& v2)
+{
+    double result;
+    gsl_blas_ddot(v1.get_gsl_vector(), v2.get_gsl_vector(), &result);
+    return result;
+}
+
+
+class VectorImpl : public BaseVectorImpl
+{
+  public:
+	typedef BaseVectorImpl base_type;
+
+  public:
+	void set_all(double x)
+	{
+		gsl_vector_set_all(m_vector, x);
+	}
+
+	void set_basis(size_t i)
+	{
+		gsl_vector_set_basis(m_vector, i);
+	}
+
+	using base_type::operator[];
+
+	double & operator[](size_t i)
+	{
+		return *gsl_vector_ptr(m_vector, i);
+	}
+
+	using base_type::get_gsl_vector;
+
+	gsl_vector* get_gsl_vector()
+	{
+		return m_vector;
+	}
+
+	void swap_elements(size_t i, size_t j)
+	{
+		gsl_vector_swap_elements(m_vector, i, j);
+	}
+
+	void reverse()
+	{
+		gsl_vector_reverse(m_vector);
+	}
+
+	VectorImpl & scale(double x)
+	{
+		gsl_vector_scale(m_vector, x);
+		return (*this);
+	}
+
+	VectorImpl & translate(double x)
+	{
+		gsl_vector_add_constant(m_vector, x);
+		return (*this);
+	}
+
+	VectorImpl & operator+=(base_type const& _vector)
+	{
+		gsl_vector_add(m_vector, _vector.get_gsl_vector());
+		return (*this);
+	}
+
+	VectorImpl & operator-=(base_type const& _vector)
+	{
+		gsl_vector_sub(m_vector, _vector.get_gsl_vector());
+		return (*this);
+	}
+
+};  // end class VectorImpl
+
+}  // end namespace detail
+
+
+using detail::operator==;
+using detail::operator<<;
+
+class Vector : public detail::VectorImpl
+{
+  public:
+	typedef detail::VectorImpl base_type;
+
+  public:
+	Vector(size_t n)
+	{
+		m_size = n;
+		m_vector = gsl_vector_alloc(n);
+	}
+
+	Vector(size_t n, double x)
+	{
+		m_size = n;
+		m_vector = gsl_vector_alloc(n);
+		gsl_vector_set_all(m_vector, x);
+	}
+
+	// create a vector with n elements all set to zero
+	// but the i-th that is set to 1
+	Vector(size_t n, size_t i)
+	{
+		m_size = n;
+		m_vector = gsl_vector_alloc(n);
+		gsl_vector_set_basis(m_vector, i);
+	}
+
+	Vector(Vector const& _vector)
+        : base_type()
+	{
+		m_size = _vector.size();
+		m_vector = gsl_vector_alloc(size());
+		gsl_vector_memcpy(m_vector, _vector.m_vector);
+	}
+
+	explicit
+	Vector(base_type::base_type const& _vector)
+	{
+		m_size = _vector.size();
+		m_vector = gsl_vector_alloc(size());
+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+	}
+
+	virtual ~Vector()
+	{
+		gsl_vector_free(m_vector);
+	}
+
+
+	Vector & operator=(Vector const& _vector)
+	{
+		assert( size() == _vector.size() );
+		gsl_vector_memcpy(m_vector, _vector.m_vector);
+		return (*this);
+	}
+
+	Vector & operator=(base_type::base_type const& _vector)
+	{
+		assert( size() == _vector.size() );
+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+		return (*this);
+	}
+
+	Vector & scale(double x)
+	{
+		return static_cast<Vector&>( base_type::scale(x) );
+	}
+
+	Vector & translate(double x)
+	{
+		return static_cast<Vector&>( base_type::translate(x) );
+	}
+
+	Vector & operator+=(base_type::base_type const& _vector)
+	{
+		return static_cast<Vector&>( base_type::operator+=(_vector) );
+	}
+
+	Vector & operator-=(base_type::base_type const& _vector)
+	{
+		return static_cast<Vector&>( base_type::operator-=(_vector) );
+	}
+
+	friend
+	void swap(Vector & v1, Vector & v2);
+	friend
+	void swap_any(Vector & v1, Vector & v2);
+
+}; // end class Vector
+
+
+// warning! these operations invalidate any view of the passed vector objects
+inline
+void swap(Vector & v1, Vector & v2)
+{
+	assert( v1.size() == v2.size() );
+	std::swap(v1.m_vector, v2.m_vector);
+}
+
+inline
+void swap_any(Vector & v1, Vector & v2)
+{
+    std::swap(v1.m_vector, v2.m_vector);
+    std::swap(v1.m_size, v2.m_size);
+}
+
+
+class ConstVectorView : public detail::BaseVectorImpl
+{
+  public:
+	typedef detail::BaseVectorImpl base_type;
+
+  public:
+	ConstVectorView(const base_type & _vector, size_t n, size_t offset = 0)
+		: m_vector_view( gsl_vector_const_subvector(_vector.get_gsl_vector(), offset, n) )
+	{
+		m_size = n;
+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+	}
+
+	ConstVectorView(const base_type & _vector, size_t n, size_t offset , size_t stride)
+		: m_vector_view( gsl_vector_const_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n) )
+	{
+		m_size = n;
+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+	}
+
+    ConstVectorView(const double* _vector, size_t n, size_t offset = 0)
+        : m_vector_view( gsl_vector_const_view_array(_vector + offset, n) )
+    {
+        m_size = n;
+        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+    }
+
+    ConstVectorView(const double* _vector, size_t n, size_t offset, size_t stride)
+        : m_vector_view( gsl_vector_const_view_array_with_stride(_vector + offset, stride, n) )
+    {
+        m_size = n;
+        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+    }
+
+	explicit
+	ConstVectorView(gsl_vector_const_view  _gsl_vector_view)
+		: m_vector_view(_gsl_vector_view)
+	{
+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+		m_size = m_vector->size;
+	}
+
+    explicit
+    ConstVectorView(const std::vector<double>&  _vector)
+        : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) )
+    {
+        m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+        m_size = _vector.size();
+    }
+
+	ConstVectorView(const ConstVectorView & _vector)
+		: base_type(),
+		  m_vector_view(_vector.m_vector_view)
+	{
+		m_size = _vector.size();
+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+	}
+
+	ConstVectorView(const base_type & _vector)
+		: m_vector_view(gsl_vector_const_subvector(_vector.get_gsl_vector(), 0, _vector.size()))
+	{
+		m_size = _vector.size();
+		m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );
+	}
+
+  private:
+	gsl_vector_const_view m_vector_view;
+
+}; // end class  ConstVectorView
+
+
+
+
+class VectorView : public detail::VectorImpl
+{
+  public:
+	typedef detail::VectorImpl base_type;
+
+  public:
+	VectorView(base_type & _vector, size_t n, size_t offset = 0, size_t stride = 1)
+	{
+		m_size = n;
+		if (stride == 1)
+		{
+			m_vector_view
+				= gsl_vector_subvector(_vector.get_gsl_vector(), offset, n);
+			m_vector = &(m_vector_view.vector);
+		}
+		else
+		{
+			m_vector_view
+				= gsl_vector_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n);
+			m_vector = &(m_vector_view.vector);
+		}
+	}
+
+    VectorView(double* _vector, size_t n, size_t offset = 0, size_t stride = 1)
+    {
+        m_size = n;
+        if (stride == 1)
+        {
+            m_vector_view
+                = gsl_vector_view_array(_vector + offset, n);
+            m_vector = &(m_vector_view.vector);
+        }
+        else
+        {
+            m_vector_view
+                = gsl_vector_view_array_with_stride(_vector + offset, stride, n);
+            m_vector = &(m_vector_view.vector);
+        }
+
+    }
+
+	VectorView(const VectorView & _vector)
+        : base_type()
+	{
+		m_size = _vector.size();
+		m_vector_view = _vector.m_vector_view;
+		m_vector = &(m_vector_view.vector);
+	}
+
+	VectorView(Vector & _vector)
+	{
+		m_size = _vector.size();
+		m_vector_view = gsl_vector_subvector(_vector.get_gsl_vector(), 0, size());
+		m_vector = &(m_vector_view.vector);
+	}
+
+	explicit
+	VectorView(gsl_vector_view _gsl_vector_view)
+		: m_vector_view(_gsl_vector_view)
+	{
+		m_vector = &(m_vector_view.vector);
+		m_size = m_vector->size;
+	}
+
+	explicit
+	VectorView(std::vector<double> & _vector)
+	{
+	    m_size = _vector.size();
+	    m_vector_view = gsl_vector_view_array(&(_vector[0]), _vector.size());
+	    m_vector = &(m_vector_view.vector);
+	}
+
+	VectorView & operator=(VectorView const& _vector)
+	{
+		assert( size() == _vector.size() );
+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+		return (*this);
+	}
+
+	VectorView & operator=(base_type::base_type const& _vector)
+	{
+		assert( size() == _vector.size() );
+		gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());
+		return (*this);
+	}
+
+	VectorView & scale(double x)
+	{
+		return static_cast<VectorView&>( base_type::scale(x) );
+	}
+
+	VectorView & translate(double x)
+	{
+		return static_cast<VectorView&>( base_type::translate(x) );
+	}
+
+	VectorView & operator+=(base_type::base_type const& _vector)
+	{
+		return static_cast<VectorView&>( base_type::operator+=(_vector) );
+	}
+
+	VectorView & operator-=(base_type::base_type const& _vector)
+	{
+		return static_cast<VectorView&>( base_type::operator-=(_vector) );
+	}
+
+	friend
+	void swap_view(VectorView & v1, VectorView & v2);
+
+  private:
+	gsl_vector_view m_vector_view;
+
+}; // end class VectorView
+
+
+inline
+void swap_view(VectorView & v1, VectorView & v2)
+{
+	assert( v1.size() == v2.size() );
+	std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too
+}
+
+
+} } // end namespaces
+
+
+#endif /*_NL_VECTOR_H_*/
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
diff --git a/src/2geom/svg-elliptical-arc.cpp b/src/2geom/svg-elliptical-arc.cpp
index 86b919fd2..047ae1431 100644
--- a/src/2geom/svg-elliptical-arc.cpp
+++ b/src/2geom/svg-elliptical-arc.cpp
@@ -1,1124 +1,1124 @@
-/*

- * SVG Elliptical Arc Class

- *

- * Copyright 2008  Marco Cecchetti <mrcekets at gmail.com>

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-#include <2geom/svg-elliptical-arc.h>

-#include <2geom/ellipse.h>

-#include <2geom/sbasis-geometric.h>

-#include <2geom/bezier-curve.h>

-#include <2geom/poly.h>

-

-#include <cfloat>

-#include <limits>

-

-#include <2geom/numeric/vector.h>

-#include <2geom/numeric/fitting-tool.h>

-#include <2geom/numeric/fitting-model.h>

-

-

-

-namespace Geom

-{

-

-

-Rect SVGEllipticalArc::boundsExact() const

-{

-    if (isDegenerate() && is_svg_compliant())

-        return chord().boundsExact();

-

-    std::vector<double> extremes(4);

-    double cosrot = std::cos(rotation_angle());

-    double sinrot = std::sin(rotation_angle());

-    extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );

-    extremes[1] = extremes[0] + M_PI;

-    if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;

-    extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );

-    extremes[3] = extremes[2] + M_PI;

-    if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;

-

-

-    std::vector<double>arc_extremes(4);

-    arc_extremes[0] = initialPoint()[X];

-    arc_extremes[1] = finalPoint()[X];

-    if ( arc_extremes[0] < arc_extremes[1] )

-        std::swap(arc_extremes[0], arc_extremes[1]);

-    arc_extremes[2] = initialPoint()[Y];

-    arc_extremes[3] = finalPoint()[Y];

-    if ( arc_extremes[2] < arc_extremes[3] )

-        std::swap(arc_extremes[2], arc_extremes[3]);

-

-

-    if ( start_angle() < end_angle() )

-    {

-        if ( sweep_flag() )

-        {

-            for ( unsigned int i = 0; i < extremes.size(); ++i )

-            {

-                if ( start_angle() < extremes[i] && extremes[i] < end_angle() )

-                {

-                    arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];

-                }

-            }

-        }

-        else

-        {

-            for ( unsigned int i = 0; i < extremes.size(); ++i )

-            {

-                if ( start_angle() > extremes[i] || extremes[i] > end_angle() )

-                {

-                    arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];

-                }

-            }

-        }

-    }

-    else

-    {

-        if ( sweep_flag() )

-        {

-            for ( unsigned int i = 0; i < extremes.size(); ++i )

-            {

-                if ( start_angle() < extremes[i] || extremes[i] < end_angle() )

-                {

-                    arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];

-                }

-            }

-        }

-        else

-        {

-            for ( unsigned int i = 0; i < extremes.size(); ++i )

-            {

-                if ( start_angle() > extremes[i] && extremes[i] > end_angle() )

-                {

-                    arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];

-                }

-            }

-        }

-    }

-

-    return Rect( Point(arc_extremes[1], arc_extremes[3]) ,

-                 Point(arc_extremes[0], arc_extremes[2]) );

-}

-

-

-double SVGEllipticalArc::valueAtAngle(Coord t, Dim2 d) const

-{

-    double sin_rot_angle = std::sin(rotation_angle());

-    double cos_rot_angle = std::cos(rotation_angle());

-    if ( d == X )

-    {

-        return    ray(X) * cos_rot_angle * std::cos(t)

-                - ray(Y) * sin_rot_angle * std::sin(t)

-                + center(X);

-    }

-    else if ( d == Y )

-    {

-        return    ray(X) * sin_rot_angle * std::cos(t)

-                + ray(Y) * cos_rot_angle * std::sin(t)

-                + center(Y);

-    }

-    THROW_RANGEERROR("dimension parameter out of range");

-}

-

-

-std::vector<double>

-SVGEllipticalArc::roots(double v, Dim2 d) const

-{

-    if ( d > Y )

-    {

-        THROW_RANGEERROR("dimention out of range");

-    }

-

-    std::vector<double> sol;

-

-    if (isDegenerate() && is_svg_compliant())

-    {

-        return chord().roots(v, d);

-    }

-    else

-    {

-        if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )

-        {

-            if ( center(d) == v )

-                sol.push_back(0);

-            return sol;

-        }

-

-        const char* msg[2][2] =

-        {

-            { "d == X; ray(X) == 0; "

-              "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "

-              "s should be contained in [-1,1]",

-              "d == X; ray(Y) == 0; "

-              "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "

-              "s should be contained in [-1,1]"

-            },

-            { "d == Y; ray(X) == 0; "

-              "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "

-              "s should be contained in [-1,1]",

-              "d == Y; ray(Y) == 0; "

-              "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "

-              "s should be contained in [-1,1]"

-            },

-        };

-

-        for ( unsigned int dim = 0; dim < 2; ++dim )

-        {

-            if ( are_near(ray(dim), 0) )

-            {

-                if ( initialPoint()[d] == v && finalPoint()[d] == v )

-                {

-                    THROW_INFINITESOLUTIONS(0);

-                }

-                if ( (initialPoint()[d] < finalPoint()[d])

-                     && (initialPoint()[d] > v || finalPoint()[d] < v) )

-                {

-                    return sol;

-                }

-                if ( (initialPoint()[d] > finalPoint()[d])

-                     && (finalPoint()[d] > v || initialPoint()[d] < v) )

-                {

-                    return sol;

-                }

-                double ray_prj;

-                switch(d)

-                {

-                    case X:

-                        switch(dim)

-                        {

-                            case X: ray_prj = -ray(Y) * std::sin(rotation_angle());

-                                    break;

-                            case Y: ray_prj = ray(X) * std::cos(rotation_angle());

-                                    break;

-                        }

-                        break;

-                    case Y:

-                        switch(dim)

-                        {

-                            case X: ray_prj = ray(Y) * std::cos(rotation_angle());

-                                    break;

-                            case Y: ray_prj = ray(X) * std::sin(rotation_angle());

-                                    break;

-                        }

-                        break;

-                }

-

-                double s = (v - center(d)) / ray_prj;

-                if ( s < -1 || s > 1 )

-                {

-                    THROW_LOGICALERROR(msg[d][dim]);

-                }

-                switch(dim)

-                {

-                    case X:

-                        s = std::asin(s); // return a value in [-PI/2,PI/2]

-                        if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) )  )

-                        {

-                            if ( s < 0 ) s += 2*M_PI;

-                        }

-                        else

-                        {

-                            s = M_PI - s;

-                            if (!(s < 2*M_PI) ) s -= 2*M_PI;

-                        }

-                        break;

-                    case Y:

-                        s = std::acos(s); // return a value in [0,PI]

-                        if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )

-                        {

-                            s = 2*M_PI - s;

-                            if ( !(s < 2*M_PI) ) s -= 2*M_PI;

-                        }

-                        break;

-                }

-

-                //std::cerr << "s = " << rad_to_deg(s);

-                s = map_to_01(s);

-                //std::cerr << " -> t: " << s << std::endl;

-                if ( !(s < 0 || s > 1) )

-                    sol.push_back(s);

-                return sol;

-            }

-        }

-

-    }

-

-    double rotx, roty;

-    switch(d)

-    {

-        case X:

-            rotx = std::cos(rotation_angle());

-            roty = -std::sin(rotation_angle());

-            break;

-        case Y:

-            rotx = std::sin(rotation_angle());

-            roty = std::cos(rotation_angle());

-            break;

-    }

-    double rxrotx = ray(X) * rotx;

-    double c_v = center(d) - v;

-

-    double a = -rxrotx + c_v;

-    double b = ray(Y) * roty;

-    double c = rxrotx + c_v;

-    //std::cerr << "a = " << a << std::endl;

-    //std::cerr << "b = " << b << std::endl;

-    //std::cerr << "c = " << c << std::endl;

-

-    if ( are_near(a,0) )

-    {

-        sol.push_back(M_PI);

-        if ( !are_near(b,0) )

-        {

-            double s = 2 * std::atan(-c/(2*b));

-            if ( s < 0 ) s += 2*M_PI;

-            sol.push_back(s);

-        }

-    }

-    else

-    {

-        double delta = b * b - a * c;

-        //std::cerr << "delta = " << delta << std::endl;

-        if ( are_near(delta, 0) )

-        {

-            double s = 2 * std::atan(-b/a);

-            if ( s < 0 ) s += 2*M_PI;

-            sol.push_back(s);

-        }

-        else if ( delta > 0 )

-        {

-            double sq = std::sqrt(delta);

-            double s = 2 * std::atan( (-b - sq) / a );

-            if ( s < 0 ) s += 2*M_PI;

-            sol.push_back(s);

-            s = 2 * std::atan( (-b + sq) / a );

-            if ( s < 0 ) s += 2*M_PI;

-            sol.push_back(s);

-        }

-    }

-

-    std::vector<double> arc_sol;

-    for (unsigned int i = 0; i < sol.size(); ++i )

-    {

-        //std::cerr << "s = " << rad_to_deg(sol[i]);

-        sol[i] = map_to_01(sol[i]);

-        //std::cerr << " -> t: " << sol[i] << std::endl;

-        if ( !(sol[i] < 0 || sol[i] > 1) )

-            arc_sol.push_back(sol[i]);

-    }

-    return arc_sol;

-}

-

-

-// D(E(t,C),t) = E(t+PI/2,O)

-Curve* SVGEllipticalArc::derivative() const

-{

-    if (isDegenerate() && is_svg_compliant())

-            return chord().derivative();

-

-    SVGEllipticalArc* result = new SVGEllipticalArc(*this);

-    result->m_center[X] = result->m_center[Y] = 0;

-    result->m_start_angle += M_PI/2;

-    if( !( result->m_start_angle < 2*M_PI ) )

-    {

-        result->m_start_angle -= 2*M_PI;

-    }

-    result->m_end_angle += M_PI/2;

-    if( !( result->m_end_angle < 2*M_PI ) )

-    {

-        result->m_end_angle -= 2*M_PI;

-    }

-    result->m_initial_point = result->pointAtAngle( result->start_angle() );

-    result->m_final_point = result->pointAtAngle( result->end_angle() );

-    return result;

-}

-

-

-std::vector<Point>

-SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const

-{

-    if (isDegenerate() && is_svg_compliant())

-            return chord().pointAndDerivatives(t, n);

-

-    unsigned int nn = n+1; // nn represents the size of the result vector.

-    std::vector<Point> result;

-    result.reserve(nn);

-    double angle = map_unit_interval_on_circular_arc(t, start_angle(),

-                                                     end_angle(), sweep_flag());

-    SVGEllipticalArc ea(*this);

-    ea.m_center = Point(0,0);

-    unsigned int m = std::min(nn, 4u);

-    for ( unsigned int i = 0; i < m; ++i )

-    {

-        result.push_back( ea.pointAtAngle(angle) );

-        angle += M_PI/2;

-        if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;

-    }

-    m = nn / 4;

-    for ( unsigned int i = 1; i < m; ++i )

-    {

-        for ( unsigned int j = 0; j < 4; ++j )

-            result.push_back( result[j] );

-    }

-    m = nn - 4 * m;

-    for ( unsigned int i = 0; i < m; ++i )

-    {

-        result.push_back( result[i] );

-    }

-    if ( !result.empty() ) // nn != 0

-        result[0] = pointAtAngle(angle);

-    return result;

-}

-

-bool SVGEllipticalArc::containsAngle(Coord angle) const

-{

-    if ( sweep_flag() )

-        if ( start_angle() < end_angle() )

-            return ( !( angle < start_angle() || angle > end_angle() ) );

-        else

-            return ( !( angle < start_angle() && angle > end_angle() ) );

-    else

-        if ( start_angle() > end_angle() )

-            return ( !( angle > start_angle() || angle < end_angle() ) );

-        else

-            return ( !( angle > start_angle() && angle < end_angle() ) );

-}

-

-Curve* SVGEllipticalArc::portion(double f, double t) const

-{

-    if (f < 0) f = 0;

-    if (f > 1) f = 1;

-    if (t < 0) t = 0;

-    if (t > 1) t = 1;

-    if ( are_near(f, t) )

-    {

-        SVGEllipticalArc* arc = new SVGEllipticalArc();

-        arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);

-        arc->m_start_angle = arc->m_end_angle = m_start_angle;

-        arc->m_rot_angle = m_rot_angle;

-        arc->m_sweep = m_sweep;

-        arc->m_large_arc = m_large_arc;

-    }

-    SVGEllipticalArc* arc = new SVGEllipticalArc( *this );

-    arc->m_initial_point = pointAt(f);

-    arc->m_final_point = pointAt(t);

-    double sa = sweep_flag() ? sweep_angle() : -sweep_angle();

-    arc->m_start_angle = m_start_angle + sa * f;

-    if ( !(arc->m_start_angle < 2*M_PI) )

-        arc->m_start_angle -= 2*M_PI;

-    if ( arc->m_start_angle < 0 )

-        arc->m_start_angle += 2*M_PI;

-    arc->m_end_angle = m_start_angle + sa * t;

-    if ( !(arc->m_end_angle < 2*M_PI) )

-        arc->m_end_angle -= 2*M_PI;

-    if ( arc->m_end_angle < 0 )

-        arc->m_end_angle += 2*M_PI;

-    if ( f > t ) arc->m_sweep = !sweep_flag();

-    if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )

-        arc->m_large_arc = false;

-    return arc;

-}

-

-

-std::vector<double> SVGEllipticalArc::

-allNearestPoints( Point const& p, double from, double to ) const

-{

-    std::vector<double> result;

-    if (isDegenerate() && is_svg_compliant())

-    {

-        result.push_back( chord().nearestPoint(p, from, to) );

-        return result;

-    }

-

-    if ( from > to ) std::swap(from, to);

-    if ( from < 0 || to > 1 )

-    {

-        THROW_RANGEERROR("[from,to] interval out of range");

-    }

-

-    if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) )  || are_near(from, to) )

-    {

-        result.push_back(from);

-        return result;

-    }

-    else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )

-    {

-        LineSegment seg(pointAt(from), pointAt(to));

-        Point np = seg.pointAt( seg.nearestPoint(p) );

-        if ( are_near(ray(Y), 0) )

-        {

-            if ( are_near(rotation_angle(), M_PI/2)

-                 || are_near(rotation_angle(), 3*M_PI/2) )

-            {

-                result = roots(np[Y], Y);

-            }

-            else

-            {

-                result = roots(np[X], X);

-            }

-        }

-        else

-        {

-            if ( are_near(rotation_angle(), M_PI/2)

-                 || are_near(rotation_angle(), 3*M_PI/2) )

-            {

-                result = roots(np[X], X);

-            }

-            else

-            {

-                result = roots(np[Y], Y);

-            }

-        }

-        return result;

-    }

-    else if ( are_near(ray(X), ray(Y)) )

-    {

-        Point r = p - center();

-        if ( are_near(r, Point(0,0)) )

-        {

-            THROW_INFINITESOLUTIONS(0);

-        }

-        // TODO: implement case r != 0

-//      Point np = ray(X) * unit_vector(r);

-//      std::vector<double> solX = roots(np[X],X);

-//      std::vector<double> solY = roots(np[Y],Y);

-//      double t;

-//      if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))

-//      {

-//          t = solX[0];

-//      }

-//      else

-//      {

-//          t = solX[1];

-//      }

-//      if ( !(t < from || t > to) )

-//      {

-//          result.push_back(t);

-//      }

-//      else

-//      {

-//

-//      }

-    }

-

-    // solve the equation <D(E(t),t)|E(t)-p> == 0

-    // that provides min and max distance points

-    // on the ellipse E wrt the point p

-    // after the substitutions:

-    // cos(t) = (1 - s^2) / (1 + s^2)

-    // sin(t) = 2t / (1 + s^2)

-    // where s = tan(t/2)

-    // we get a 4th degree equation in s

-    /*

-     *  ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +

-     *  ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +

-     *  2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +

-     *  2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])

-     */

-

-    Point p_c = p - center();

-    double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));

-    double cosrot = std::cos( rotation_angle() );

-    double sinrot = std::sin( rotation_angle() );

-    double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);

-    Poly coeff;

-    coeff.resize(5);

-    coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );

-    coeff[3] = 2 * ( rx2_ry2 + expr1 );

-    coeff[2] = 0;

-    coeff[1] = 2 * ( -rx2_ry2 + expr1 );

-    coeff[0] = -coeff[4];

-

-//  for ( unsigned int i = 0; i < 5; ++i )

-//      std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;

-

-    std::vector<double> real_sol;

-    // gsl_poly_complex_solve raises an error

-    // if the leading coefficient is zero

-    if ( are_near(coeff[4], 0) )

-    {

-        real_sol.push_back(0);

-        if ( !are_near(coeff[3], 0) )

-        {

-            double sq = -coeff[1] / coeff[3];

-            if ( sq > 0 )

-            {

-                double s = std::sqrt(sq);

-                real_sol.push_back(s);

-                real_sol.push_back(-s);

-            }

-        }

-    }

-    else

-    {

-        real_sol = solve_reals(coeff);

-    }

-

-    for ( unsigned int i = 0; i < real_sol.size(); ++i )

-    {

-        real_sol[i] = 2 * std::atan(real_sol[i]);

-        if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;

-    }

-    // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0

-    // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity

-    if ( (real_sol.size() % 2) != 0 )

-    {

-        real_sol.push_back(M_PI);

-    }

-

-    double mindistsq1 = std::numeric_limits<double>::max();

-    double mindistsq2 = std::numeric_limits<double>::max();

-    double dsq;

-    unsigned int mi1, mi2;

-    for ( unsigned int i = 0; i < real_sol.size(); ++i )

-    {

-        dsq = distanceSq(p, pointAtAngle(real_sol[i]));

-        if ( mindistsq1 > dsq )

-        {

-            mindistsq2 = mindistsq1;

-            mi2 = mi1;

-            mindistsq1 = dsq;

-            mi1 = i;

-        }

-        else if ( mindistsq2 > dsq )

-        {

-            mindistsq2 = dsq;

-            mi2 = i;

-        }

-    }

-

-    double t = map_to_01( real_sol[mi1] );

-    if ( !(t < from || t > to) )

-    {

-        result.push_back(t);

-    }

-

-    bool second_sol = false;

-    t = map_to_01( real_sol[mi2] );

-    if ( real_sol.size() == 4 && !(t < from || t > to) )

-    {

-        if ( result.empty() || are_near(mindistsq1, mindistsq2) )

-        {

-            result.push_back(t);

-            second_sol = true;

-        }

-    }

-

-    // we need to test extreme points too

-    double dsq1 = distanceSq(p, pointAt(from));

-    double dsq2 = distanceSq(p, pointAt(to));

-    if ( second_sol )

-    {

-        if ( mindistsq2 > dsq1 )

-        {

-            result.clear();

-            result.push_back(from);

-            mindistsq2 = dsq1;

-        }

-        else if ( are_near(mindistsq2, dsq) )

-        {

-            result.push_back(from);

-        }

-        if ( mindistsq2 > dsq2 )

-        {

-            result.clear();

-            result.push_back(to);

-        }

-        else if ( are_near(mindistsq2, dsq2) )

-        {

-            result.push_back(to);

-        }

-

-    }

-    else

-    {

-        if ( result.empty() )

-        {

-            if ( are_near(dsq1, dsq2) )

-            {

-                result.push_back(from);

-                result.push_back(to);

-            }

-            else if ( dsq2 > dsq1 )

-            {

-                result.push_back(from);

-            }

-            else

-            {

-                result.push_back(to);

-            }

-        }

-    }

-

-    return result;

-}

-

-

-/*

- * NOTE: this implementation follows Standard SVG 1.1 implementation guidelines

- * for elliptical arc curves. See Appendix F.6.

- */

-void SVGEllipticalArc::calculate_center_and_extreme_angles()

-{

-    Point d = initialPoint() - finalPoint();

-

-    if (is_svg_compliant())

-    {

-        if ( initialPoint() == finalPoint() )

-        {

-            m_rx = m_ry = m_rot_angle = m_start_angle = m_end_angle = 0;

-            m_center = initialPoint();

-            m_large_arc = m_sweep = false;

-            return;

-        }

-

-        m_rx = std::fabs(m_rx);

-        m_ry = std::fabs(m_ry);

-

-        if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )

-        {

-            m_rx = L2(d) / 2;

-            m_ry = 0;

-            m_rot_angle = std::atan2(d[Y], d[X]);

-            if (m_rot_angle < 0) m_rot_angle += 2*M_PI;

-            m_start_angle = 0;

-            m_end_angle = M_PI;

-            m_center = middle_point(initialPoint(), finalPoint());

-            m_large_arc = false;

-            m_sweep = false;

-            return;

-        }

-    }

-    else

-    {

-        if ( are_near(initialPoint(), finalPoint()) )

-        {

-            if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )

-            {

-                m_start_angle = m_end_angle = 0;

-                m_center = initialPoint();

-                return;

-            }

-            else

-            {

-                THROW_RANGEERROR("initial and final point are the same");

-            }

-        }

-        if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )

-        { // but initialPoint != finalPoint

-            THROW_RANGEERROR(

-                "there is no ellipse that satisfies the given constraints: "

-                "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"

-            );

-        }

-        if ( are_near(ray(Y), 0) )

-        {

-            Point v = initialPoint() - finalPoint();

-            if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )

-            {

-                double angle = std::atan2(v[Y], v[X]);

-                if (angle < 0) angle += 2*M_PI;

-                if ( are_near( angle, rotation_angle() ) )

-                {

-                    m_start_angle = 0;

-                    m_end_angle = M_PI;

-                    m_center = v/2 + finalPoint();

-                    return;

-                }

-                angle -= M_PI;

-                if ( angle < 0 ) angle += 2*M_PI;

-                if ( are_near( angle, rotation_angle() ) )

-                {

-                    m_start_angle = M_PI;

-                    m_end_angle = 0;

-                    m_center = v/2 + finalPoint();

-                    return;

-                }

-                THROW_RANGEERROR(

-                    "there is no ellipse that satisfies the given constraints: "

-                    "ray(Y) == 0 "

-                    "and slope(initialPoint - finalPoint) != rotation_angle "

-                    "and != rotation_angle + PI"

-                );

-            }

-            if ( L2sq(v) > 4*ray(X)*ray(X) )

-            {

-                THROW_RANGEERROR(

-                    "there is no ellipse that satisfies the given constraints: "

-                    "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"

-                );

-            }

-            else

-            {

-                THROW_RANGEERROR(

-                    "there is infinite ellipses that satisfy the given constraints: "

-                    "ray(Y) == 0  and distance(initialPoint, finalPoint) < 2*ray(X)"

-                );

-            }

-

-        }

-

-        if ( are_near(ray(X), 0) )

-        {

-            Point v = initialPoint() - finalPoint();

-            if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )

-            {

-                double angle = std::atan2(v[Y], v[X]);

-                if (angle < 0) angle += 2*M_PI;

-                double rot_angle = rotation_angle() + M_PI/2;

-                if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;

-                if ( are_near( angle, rot_angle ) )

-                {

-                    m_start_angle = M_PI/2;

-                    m_end_angle = 3*M_PI/2;

-                    m_center = v/2 + finalPoint();

-                    return;

-                }

-                angle -= M_PI;

-                if ( angle < 0 ) angle += 2*M_PI;

-                if ( are_near( angle, rot_angle ) )

-                {

-                    m_start_angle = 3*M_PI/2;

-                    m_end_angle = M_PI/2;

-                    m_center = v/2 + finalPoint();

-                    return;

-                }

-                THROW_RANGEERROR(

-                    "there is no ellipse that satisfies the given constraints: "

-                    "ray(X) == 0 "

-                    "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "

-                    "and != rotation_angle + (3/2)*PI"

-                );

-            }

-            if ( L2sq(v) > 4*ray(Y)*ray(Y) )

-            {

-                THROW_RANGEERROR(

-                    "there is no ellipse that satisfies the given constraints: "

-                    "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"

-                );

-            }

-            else

-            {

-                THROW_RANGEERROR(

-                    "there is infinite ellipses that satisfy the given constraints: "

-                    "ray(X) == 0  and distance(initialPoint, finalPoint) < 2*ray(Y)"

-                );

-            }

-

-        }

-

-    }

-

-    double sin_rot_angle = std::sin(rotation_angle());

-    double cos_rot_angle = std::cos(rotation_angle());

-

-

-    Matrix m( cos_rot_angle, -sin_rot_angle,

-              sin_rot_angle, cos_rot_angle,

-              0,             0              );

-

-    Point p = (d / 2) * m;

-    double rx2 = m_rx * m_rx;

-    double ry2 = m_ry * m_ry;

-    double rxpy = m_rx * p[Y];

-    double rypx = m_ry * p[X];

-    double rx2py2 = rxpy * rxpy;

-    double ry2px2 = rypx * rypx;

-    double num = rx2 * ry2;

-    double den = rx2py2 + ry2px2;

-    assert(den != 0);

-    double rad = num / den;

-    Point c(0,0);

-    if (rad > 1)

-    {

-        rad -= 1;

-        rad = std::sqrt(rad);

-

-        if (m_large_arc == m_sweep) rad = -rad;

-        c = rad * Point(rxpy / m_ry, -rypx / m_rx);

-

-        m[1] = -m[1];

-        m[2] = -m[2];

-

-        m_center = c * m + middle_point(initialPoint(), finalPoint());

-    }

-    else if (rad == 1 || is_svg_compliant())

-    {

-        double lamda = std::sqrt(1 / rad);

-        m_rx *= lamda;

-        m_ry *= lamda;

-        m_center = middle_point(initialPoint(), finalPoint());

-    }

-    else

-    {

-        THROW_RANGEERROR(

-            "there is no ellipse that satisfies the given constraints"

-        );

-    }

-

-    Point sp((p[X] - c[X]) / m_rx, (p[Y] - c[Y]) / m_ry);

-    Point ep((-p[X] - c[X]) / m_rx, (-p[Y] - c[Y]) / m_ry);

-    Point v(1, 0);

-    m_start_angle = angle_between(v, sp);

-    double sweep_angle = angle_between(sp, ep);

-    if (!m_sweep && sweep_angle > 0) sweep_angle -= 2*M_PI;

-    if (m_sweep && sweep_angle < 0) sweep_angle += 2*M_PI;

-

-    if (m_start_angle < 0) m_start_angle += 2*M_PI;

-    m_end_angle = m_start_angle + sweep_angle;

-    if (m_end_angle < 0) m_end_angle += 2*M_PI;

-    if (m_end_angle >= 2*M_PI) m_end_angle -= 2*M_PI;

-}

-

-

-D2<SBasis> SVGEllipticalArc::toSBasis() const

-{

-

-    if (isDegenerate() && is_svg_compliant())

-        return chord().toSBasis();

-

-    D2<SBasis> arc;

-    // the interval of parametrization has to be [0,1]

-    Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );

-    Linear param(start_angle(), et);

-    Coord cos_rot_angle = std::cos(rotation_angle());

-    Coord sin_rot_angle = std::sin(rotation_angle());

-    // order = 4 seems to be enough to get a perfect looking elliptical arc

-    // should it be choosen in function of the arc length anyway ?

-    // or maybe a user settable parameter: toSBasis(unsigned int order) ?

-    SBasis arc_x = ray(X) * cos(param,4);

-    SBasis arc_y = ray(Y) * sin(param,4);

-    arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));

-    arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));

-    return arc;

-}

-

-

-Coord SVGEllipticalArc::map_to_02PI(Coord t) const

-{

-    if ( sweep_flag() )

-    {

-        Coord angle = start_angle() + sweep_angle() * t;

-        if ( !(angle < 2*M_PI) )

-            angle -= 2*M_PI;

-        return angle;

-    }

-    else

-    {

-        Coord angle = start_angle() - sweep_angle() * t;

-        if ( angle < 0 ) angle += 2*M_PI;

-        return angle;

-    }

-}

-

-Coord SVGEllipticalArc::map_to_01(Coord angle) const

-{

-    return map_circular_arc_on_unit_interval(angle, start_angle(),

-                                             end_angle(), sweep_flag());

-}

-

-

-namespace detail

-{

-

-struct ellipse_equation

-{

-    ellipse_equation(double a, double b, double c, double d, double e, double f)

-        : A(a), B(b), C(c), D(d), E(e), F(f)

-    {

-    }

-

-    double operator()(double x, double y) const

-    {

-        // A * x * x + B * x * y + C * y * y + D * x + E * y + F;

-        return (A * x + B * y + D) * x + (C * y + E) * y + F;

-    }

-

-    double operator()(Point const& p) const

-    {

-        return (*this)(p[X], p[Y]);

-    }

-

-    Point normal(double x, double y) const

-    {

-        Point n( 2 * A * x + B * y + D, 2 * C * y + B * x + E );

-        return unit_vector(n);

-    }

-

-    Point normal(Point const& p) const

-    {

-        return normal(p[X], p[Y]);

-    }

-

-    double A, B, C, D, E, F;

-};

-

-}

-

-make_elliptical_arc::

-make_elliptical_arc( SVGEllipticalArc& _ea,

-                     curve_type const& _curve,

-                     unsigned int _total_samples,

-                     double _tolerance )

-    : ea(_ea), curve(_curve),

-      dcurve( unitVector(derivative(curve)) ),

-      model(), fitter(model, _total_samples),

-      tolerance(_tolerance), tol_at_extr(tolerance/2),

-      tol_at_center(0.1), angle_tol(0.1),

-      initial_point(curve.at0()), final_point(curve.at1()),

-      N(_total_samples), last(N-1), partitions(N-1), p(N),

-      svg_compliant(true)

-{

-}

-

-bool

-make_elliptical_arc::

-bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,

-                double e1x, double e1y, double e2 )

-{

-    dist_err = std::fabs( ee(p[k]) );

-    dist_bound = std::fabs( e1x * p[k][X] + e1y * p[k][Y] + e2 );

-    angle_err = std::fabs( dot( dcurve(k/partitions), ee.normal(p[k]) ) );

-    //angle_err *= angle_err;

-    return ( dist_err  > dist_bound || angle_err > angle_tol );

-}

-

-bool

-make_elliptical_arc::

-check_bound(double A, double B, double C, double D, double E, double F)

-{

-    // check error magnitude

-    detail::ellipse_equation ee(A, B, C, D, E, F);

-

-    double e1x = (2*A + B) * tol_at_extr;

-    double e1y = (B + 2*C) * tol_at_extr;

-    double e2 = ((D + E)  + (A + B + C) * tol_at_extr) * tol_at_extr;

-    if ( bound_exceeded(0, ee, e1x, e1y, e2) )

-    {

-        print_bound_error(0);

-        return false;

-    }

-    if ( bound_exceeded(0, ee, e1x, e1y, e2) )

-    {

-        print_bound_error(last);

-        return false;

-    }

-

-    e1x = (2*A + B) * tolerance;

-    e1y = (B + 2*C) * tolerance;

-    e2 = ((D + E)  + (A + B + C) * tolerance) * tolerance;

-//  std::cerr << "e1x = " << e1x << std::endl;

-//  std::cerr << "e1y = " << e1y << std::endl;

-//  std::cerr << "e2 = " << e2 << std::endl;

-

-    for ( unsigned int k = 1; k < last; ++k )

-    {

-        if ( bound_exceeded(k, ee, e1x, e1y, e2) )

-        {

-            print_bound_error(k);

-            return false;

-        }

-    }

-

-    return true;

-}

-

-void make_elliptical_arc::fit()

-{

-    for (unsigned int k = 0; k < N; ++k)

-    {

-        p[k] = curve( k / partitions );

-        fitter.append(p[k]);

-    }

-    fitter.update();

-

-    NL::Vector z(N, 0.0);

-    fitter.result(z);

-}

-

-bool make_elliptical_arc::make_elliptiarc()

-{

-    const NL::Vector & coeff = fitter.result();

-    Ellipse e;

-    try

-    {

-        e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);

-    }

-    catch(LogicalError exc)

-    {

-        return false;

-    }

-

-    Point inner_point = curve(0.5);

-

-    if (svg_compliant_flag())

-    {

-        ea = e.arc(initial_point, inner_point, final_point);

-    }

-    else

-    {

-        try

-        {

-            ea = e.arc(initial_point, inner_point, final_point, false);

-        }

-        catch(RangeError exc)

-        {

-            return false;

-        }

-    }

-

-

-    if ( !are_near( e.center(),

-                    ea.center(),

-                    tol_at_center * std::min(e.ray(X),e.ray(Y))

-                  )

-       )

-    {

-        return false;

-    }

-    return true;

-}

-

-

-

-} // end namespace Geom

-

-

-

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

-

+/*
+ * SVG Elliptical Arc Class
+ *
+ * Copyright 2008  Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include <2geom/svg-elliptical-arc.h>
+#include <2geom/ellipse.h>
+#include <2geom/sbasis-geometric.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/poly.h>
+
+#include <cfloat>
+#include <limits>
+
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/fitting-tool.h>
+#include <2geom/numeric/fitting-model.h>
+
+
+
+namespace Geom
+{
+
+
+Rect SVGEllipticalArc::boundsExact() const
+{
+    if (isDegenerate() && is_svg_compliant())
+        return chord().boundsExact();
+
+    std::vector<double> extremes(4);
+    double cosrot = std::cos(rotation_angle());
+    double sinrot = std::sin(rotation_angle());
+    extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
+    extremes[1] = extremes[0] + M_PI;
+    if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;
+    extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
+    extremes[3] = extremes[2] + M_PI;
+    if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;
+
+
+    std::vector<double>arc_extremes(4);
+    arc_extremes[0] = initialPoint()[X];
+    arc_extremes[1] = finalPoint()[X];
+    if ( arc_extremes[0] < arc_extremes[1] )
+        std::swap(arc_extremes[0], arc_extremes[1]);
+    arc_extremes[2] = initialPoint()[Y];
+    arc_extremes[3] = finalPoint()[Y];
+    if ( arc_extremes[2] < arc_extremes[3] )
+        std::swap(arc_extremes[2], arc_extremes[3]);
+
+
+    if ( start_angle() < end_angle() )
+    {
+        if ( sweep_flag() )
+        {
+            for ( unsigned int i = 0; i < extremes.size(); ++i )
+            {
+                if ( start_angle() < extremes[i] && extremes[i] < end_angle() )
+                {
+                    arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+                }
+            }
+        }
+        else
+        {
+            for ( unsigned int i = 0; i < extremes.size(); ++i )
+            {
+                if ( start_angle() > extremes[i] || extremes[i] > end_angle() )
+                {
+                    arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+                }
+            }
+        }
+    }
+    else
+    {
+        if ( sweep_flag() )
+        {
+            for ( unsigned int i = 0; i < extremes.size(); ++i )
+            {
+                if ( start_angle() < extremes[i] || extremes[i] < end_angle() )
+                {
+                    arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+                }
+            }
+        }
+        else
+        {
+            for ( unsigned int i = 0; i < extremes.size(); ++i )
+            {
+                if ( start_angle() > extremes[i] && extremes[i] > end_angle() )
+                {
+                    arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+                }
+            }
+        }
+    }
+
+    return Rect( Point(arc_extremes[1], arc_extremes[3]) ,
+                 Point(arc_extremes[0], arc_extremes[2]) );
+}
+
+
+double SVGEllipticalArc::valueAtAngle(Coord t, Dim2 d) const
+{
+    double sin_rot_angle = std::sin(rotation_angle());
+    double cos_rot_angle = std::cos(rotation_angle());
+    if ( d == X )
+    {
+        return    ray(X) * cos_rot_angle * std::cos(t)
+                - ray(Y) * sin_rot_angle * std::sin(t)
+                + center(X);
+    }
+    else if ( d == Y )
+    {
+        return    ray(X) * sin_rot_angle * std::cos(t)
+                + ray(Y) * cos_rot_angle * std::sin(t)
+                + center(Y);
+    }
+    THROW_RANGEERROR("dimension parameter out of range");
+}
+
+
+std::vector<double>
+SVGEllipticalArc::roots(double v, Dim2 d) const
+{
+    if ( d > Y )
+    {
+        THROW_RANGEERROR("dimention out of range");
+    }
+
+    std::vector<double> sol;
+
+    if (isDegenerate() && is_svg_compliant())
+    {
+        return chord().roots(v, d);
+    }
+    else
+    {
+        if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
+        {
+            if ( center(d) == v )
+                sol.push_back(0);
+            return sol;
+        }
+
+        const char* msg[2][2] =
+        {
+            { "d == X; ray(X) == 0; "
+              "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "
+              "s should be contained in [-1,1]",
+              "d == X; ray(Y) == 0; "
+              "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "
+              "s should be contained in [-1,1]"
+            },
+            { "d == Y; ray(X) == 0; "
+              "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "
+              "s should be contained in [-1,1]",
+              "d == Y; ray(Y) == 0; "
+              "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "
+              "s should be contained in [-1,1]"
+            },
+        };
+
+        for ( unsigned int dim = 0; dim < 2; ++dim )
+        {
+            if ( are_near(ray(dim), 0) )
+            {
+                if ( initialPoint()[d] == v && finalPoint()[d] == v )
+                {
+                    THROW_INFINITESOLUTIONS(0);
+                }
+                if ( (initialPoint()[d] < finalPoint()[d])
+                     && (initialPoint()[d] > v || finalPoint()[d] < v) )
+                {
+                    return sol;
+                }
+                if ( (initialPoint()[d] > finalPoint()[d])
+                     && (finalPoint()[d] > v || initialPoint()[d] < v) )
+                {
+                    return sol;
+                }
+                double ray_prj;
+                switch(d)
+                {
+                    case X:
+                        switch(dim)
+                        {
+                            case X: ray_prj = -ray(Y) * std::sin(rotation_angle());
+                                    break;
+                            case Y: ray_prj = ray(X) * std::cos(rotation_angle());
+                                    break;
+                        }
+                        break;
+                    case Y:
+                        switch(dim)
+                        {
+                            case X: ray_prj = ray(Y) * std::cos(rotation_angle());
+                                    break;
+                            case Y: ray_prj = ray(X) * std::sin(rotation_angle());
+                                    break;
+                        }
+                        break;
+                }
+
+                double s = (v - center(d)) / ray_prj;
+                if ( s < -1 || s > 1 )
+                {
+                    THROW_LOGICALERROR(msg[d][dim]);
+                }
+                switch(dim)
+                {
+                    case X:
+                        s = std::asin(s); // return a value in [-PI/2,PI/2]
+                        if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) )  )
+                        {
+                            if ( s < 0 ) s += 2*M_PI;
+                        }
+                        else
+                        {
+                            s = M_PI - s;
+                            if (!(s < 2*M_PI) ) s -= 2*M_PI;
+                        }
+                        break;
+                    case Y:
+                        s = std::acos(s); // return a value in [0,PI]
+                        if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )
+                        {
+                            s = 2*M_PI - s;
+                            if ( !(s < 2*M_PI) ) s -= 2*M_PI;
+                        }
+                        break;
+                }
+
+                //std::cerr << "s = " << rad_to_deg(s);
+                s = map_to_01(s);
+                //std::cerr << " -> t: " << s << std::endl;
+                if ( !(s < 0 || s > 1) )
+                    sol.push_back(s);
+                return sol;
+            }
+        }
+
+    }
+
+    double rotx, roty;
+    switch(d)
+    {
+        case X:
+            rotx = std::cos(rotation_angle());
+            roty = -std::sin(rotation_angle());
+            break;
+        case Y:
+            rotx = std::sin(rotation_angle());
+            roty = std::cos(rotation_angle());
+            break;
+    }
+    double rxrotx = ray(X) * rotx;
+    double c_v = center(d) - v;
+
+    double a = -rxrotx + c_v;
+    double b = ray(Y) * roty;
+    double c = rxrotx + c_v;
+    //std::cerr << "a = " << a << std::endl;
+    //std::cerr << "b = " << b << std::endl;
+    //std::cerr << "c = " << c << std::endl;
+
+    if ( are_near(a,0) )
+    {
+        sol.push_back(M_PI);
+        if ( !are_near(b,0) )
+        {
+            double s = 2 * std::atan(-c/(2*b));
+            if ( s < 0 ) s += 2*M_PI;
+            sol.push_back(s);
+        }
+    }
+    else
+    {
+        double delta = b * b - a * c;
+        //std::cerr << "delta = " << delta << std::endl;
+        if ( are_near(delta, 0) )
+        {
+            double s = 2 * std::atan(-b/a);
+            if ( s < 0 ) s += 2*M_PI;
+            sol.push_back(s);
+        }
+        else if ( delta > 0 )
+        {
+            double sq = std::sqrt(delta);
+            double s = 2 * std::atan( (-b - sq) / a );
+            if ( s < 0 ) s += 2*M_PI;
+            sol.push_back(s);
+            s = 2 * std::atan( (-b + sq) / a );
+            if ( s < 0 ) s += 2*M_PI;
+            sol.push_back(s);
+        }
+    }
+
+    std::vector<double> arc_sol;
+    for (unsigned int i = 0; i < sol.size(); ++i )
+    {
+        //std::cerr << "s = " << rad_to_deg(sol[i]);
+        sol[i] = map_to_01(sol[i]);
+        //std::cerr << " -> t: " << sol[i] << std::endl;
+        if ( !(sol[i] < 0 || sol[i] > 1) )
+            arc_sol.push_back(sol[i]);
+    }
+    return arc_sol;
+}
+
+
+// D(E(t,C),t) = E(t+PI/2,O)
+Curve* SVGEllipticalArc::derivative() const
+{
+    if (isDegenerate() && is_svg_compliant())
+            return chord().derivative();
+
+    SVGEllipticalArc* result = new SVGEllipticalArc(*this);
+    result->m_center[X] = result->m_center[Y] = 0;
+    result->m_start_angle += M_PI/2;
+    if( !( result->m_start_angle < 2*M_PI ) )
+    {
+        result->m_start_angle -= 2*M_PI;
+    }
+    result->m_end_angle += M_PI/2;
+    if( !( result->m_end_angle < 2*M_PI ) )
+    {
+        result->m_end_angle -= 2*M_PI;
+    }
+    result->m_initial_point = result->pointAtAngle( result->start_angle() );
+    result->m_final_point = result->pointAtAngle( result->end_angle() );
+    return result;
+}
+
+
+std::vector<Point>
+SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
+{
+    if (isDegenerate() && is_svg_compliant())
+            return chord().pointAndDerivatives(t, n);
+
+    unsigned int nn = n+1; // nn represents the size of the result vector.
+    std::vector<Point> result;
+    result.reserve(nn);
+    double angle = map_unit_interval_on_circular_arc(t, start_angle(),
+                                                     end_angle(), sweep_flag());
+    SVGEllipticalArc ea(*this);
+    ea.m_center = Point(0,0);
+    unsigned int m = std::min(nn, 4u);
+    for ( unsigned int i = 0; i < m; ++i )
+    {
+        result.push_back( ea.pointAtAngle(angle) );
+        angle += M_PI/2;
+        if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;
+    }
+    m = nn / 4;
+    for ( unsigned int i = 1; i < m; ++i )
+    {
+        for ( unsigned int j = 0; j < 4; ++j )
+            result.push_back( result[j] );
+    }
+    m = nn - 4 * m;
+    for ( unsigned int i = 0; i < m; ++i )
+    {
+        result.push_back( result[i] );
+    }
+    if ( !result.empty() ) // nn != 0
+        result[0] = pointAtAngle(angle);
+    return result;
+}
+
+bool SVGEllipticalArc::containsAngle(Coord angle) const
+{
+    if ( sweep_flag() )
+        if ( start_angle() < end_angle() )
+            return ( !( angle < start_angle() || angle > end_angle() ) );
+        else
+            return ( !( angle < start_angle() && angle > end_angle() ) );
+    else
+        if ( start_angle() > end_angle() )
+            return ( !( angle > start_angle() || angle < end_angle() ) );
+        else
+            return ( !( angle > start_angle() && angle < end_angle() ) );
+}
+
+Curve* SVGEllipticalArc::portion(double f, double t) const
+{
+    if (f < 0) f = 0;
+    if (f > 1) f = 1;
+    if (t < 0) t = 0;
+    if (t > 1) t = 1;
+    if ( are_near(f, t) )
+    {
+        SVGEllipticalArc* arc = new SVGEllipticalArc();
+        arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);
+        arc->m_start_angle = arc->m_end_angle = m_start_angle;
+        arc->m_rot_angle = m_rot_angle;
+        arc->m_sweep = m_sweep;
+        arc->m_large_arc = m_large_arc;
+    }
+    SVGEllipticalArc* arc = new SVGEllipticalArc( *this );
+    arc->m_initial_point = pointAt(f);
+    arc->m_final_point = pointAt(t);
+    double sa = sweep_flag() ? sweep_angle() : -sweep_angle();
+    arc->m_start_angle = m_start_angle + sa * f;
+    if ( !(arc->m_start_angle < 2*M_PI) )
+        arc->m_start_angle -= 2*M_PI;
+    if ( arc->m_start_angle < 0 )
+        arc->m_start_angle += 2*M_PI;
+    arc->m_end_angle = m_start_angle + sa * t;
+    if ( !(arc->m_end_angle < 2*M_PI) )
+        arc->m_end_angle -= 2*M_PI;
+    if ( arc->m_end_angle < 0 )
+        arc->m_end_angle += 2*M_PI;
+    if ( f > t ) arc->m_sweep = !sweep_flag();
+    if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )
+        arc->m_large_arc = false;
+    return arc;
+}
+
+
+std::vector<double> SVGEllipticalArc::
+allNearestPoints( Point const& p, double from, double to ) const
+{
+    std::vector<double> result;
+    if (isDegenerate() && is_svg_compliant())
+    {
+        result.push_back( chord().nearestPoint(p, from, to) );
+        return result;
+    }
+
+    if ( from > to ) std::swap(from, to);
+    if ( from < 0 || to > 1 )
+    {
+        THROW_RANGEERROR("[from,to] interval out of range");
+    }
+
+    if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) )  || are_near(from, to) )
+    {
+        result.push_back(from);
+        return result;
+    }
+    else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
+    {
+        LineSegment seg(pointAt(from), pointAt(to));
+        Point np = seg.pointAt( seg.nearestPoint(p) );
+        if ( are_near(ray(Y), 0) )
+        {
+            if ( are_near(rotation_angle(), M_PI/2)
+                 || are_near(rotation_angle(), 3*M_PI/2) )
+            {
+                result = roots(np[Y], Y);
+            }
+            else
+            {
+                result = roots(np[X], X);
+            }
+        }
+        else
+        {
+            if ( are_near(rotation_angle(), M_PI/2)
+                 || are_near(rotation_angle(), 3*M_PI/2) )
+            {
+                result = roots(np[X], X);
+            }
+            else
+            {
+                result = roots(np[Y], Y);
+            }
+        }
+        return result;
+    }
+    else if ( are_near(ray(X), ray(Y)) )
+    {
+        Point r = p - center();
+        if ( are_near(r, Point(0,0)) )
+        {
+            THROW_INFINITESOLUTIONS(0);
+        }
+        // TODO: implement case r != 0
+//      Point np = ray(X) * unit_vector(r);
+//      std::vector<double> solX = roots(np[X],X);
+//      std::vector<double> solY = roots(np[Y],Y);
+//      double t;
+//      if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))
+//      {
+//          t = solX[0];
+//      }
+//      else
+//      {
+//          t = solX[1];
+//      }
+//      if ( !(t < from || t > to) )
+//      {
+//          result.push_back(t);
+//      }
+//      else
+//      {
+//
+//      }
+    }
+
+    // solve the equation <D(E(t),t)|E(t)-p> == 0
+    // that provides min and max distance points
+    // on the ellipse E wrt the point p
+    // after the substitutions:
+    // cos(t) = (1 - s^2) / (1 + s^2)
+    // sin(t) = 2t / (1 + s^2)
+    // where s = tan(t/2)
+    // we get a 4th degree equation in s
+    /*
+     *  ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +
+     *  ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +
+     *  2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +
+     *  2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])
+     */
+
+    Point p_c = p - center();
+    double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
+    double cosrot = std::cos( rotation_angle() );
+    double sinrot = std::sin( rotation_angle() );
+    double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
+    Poly coeff;
+    coeff.resize(5);
+    coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );
+    coeff[3] = 2 * ( rx2_ry2 + expr1 );
+    coeff[2] = 0;
+    coeff[1] = 2 * ( -rx2_ry2 + expr1 );
+    coeff[0] = -coeff[4];
+
+//  for ( unsigned int i = 0; i < 5; ++i )
+//      std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;
+
+    std::vector<double> real_sol;
+    // gsl_poly_complex_solve raises an error
+    // if the leading coefficient is zero
+    if ( are_near(coeff[4], 0) )
+    {
+        real_sol.push_back(0);
+        if ( !are_near(coeff[3], 0) )
+        {
+            double sq = -coeff[1] / coeff[3];
+            if ( sq > 0 )
+            {
+                double s = std::sqrt(sq);
+                real_sol.push_back(s);
+                real_sol.push_back(-s);
+            }
+        }
+    }
+    else
+    {
+        real_sol = solve_reals(coeff);
+    }
+
+    for ( unsigned int i = 0; i < real_sol.size(); ++i )
+    {
+        real_sol[i] = 2 * std::atan(real_sol[i]);
+        if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;
+    }
+    // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0
+    // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity
+    if ( (real_sol.size() % 2) != 0 )
+    {
+        real_sol.push_back(M_PI);
+    }
+
+    double mindistsq1 = std::numeric_limits<double>::max();
+    double mindistsq2 = std::numeric_limits<double>::max();
+    double dsq;
+    unsigned int mi1, mi2;
+    for ( unsigned int i = 0; i < real_sol.size(); ++i )
+    {
+        dsq = distanceSq(p, pointAtAngle(real_sol[i]));
+        if ( mindistsq1 > dsq )
+        {
+            mindistsq2 = mindistsq1;
+            mi2 = mi1;
+            mindistsq1 = dsq;
+            mi1 = i;
+        }
+        else if ( mindistsq2 > dsq )
+        {
+            mindistsq2 = dsq;
+            mi2 = i;
+        }
+    }
+
+    double t = map_to_01( real_sol[mi1] );
+    if ( !(t < from || t > to) )
+    {
+        result.push_back(t);
+    }
+
+    bool second_sol = false;
+    t = map_to_01( real_sol[mi2] );
+    if ( real_sol.size() == 4 && !(t < from || t > to) )
+    {
+        if ( result.empty() || are_near(mindistsq1, mindistsq2) )
+        {
+            result.push_back(t);
+            second_sol = true;
+        }
+    }
+
+    // we need to test extreme points too
+    double dsq1 = distanceSq(p, pointAt(from));
+    double dsq2 = distanceSq(p, pointAt(to));
+    if ( second_sol )
+    {
+        if ( mindistsq2 > dsq1 )
+        {
+            result.clear();
+            result.push_back(from);
+            mindistsq2 = dsq1;
+        }
+        else if ( are_near(mindistsq2, dsq) )
+        {
+            result.push_back(from);
+        }
+        if ( mindistsq2 > dsq2 )
+        {
+            result.clear();
+            result.push_back(to);
+        }
+        else if ( are_near(mindistsq2, dsq2) )
+        {
+            result.push_back(to);
+        }
+
+    }
+    else
+    {
+        if ( result.empty() )
+        {
+            if ( are_near(dsq1, dsq2) )
+            {
+                result.push_back(from);
+                result.push_back(to);
+            }
+            else if ( dsq2 > dsq1 )
+            {
+                result.push_back(from);
+            }
+            else
+            {
+                result.push_back(to);
+            }
+        }
+    }
+
+    return result;
+}
+
+
+/*
+ * NOTE: this implementation follows Standard SVG 1.1 implementation guidelines
+ * for elliptical arc curves. See Appendix F.6.
+ */
+void SVGEllipticalArc::calculate_center_and_extreme_angles()
+{
+    Point d = initialPoint() - finalPoint();
+
+    if (is_svg_compliant())
+    {
+        if ( initialPoint() == finalPoint() )
+        {
+            m_rx = m_ry = m_rot_angle = m_start_angle = m_end_angle = 0;
+            m_center = initialPoint();
+            m_large_arc = m_sweep = false;
+            return;
+        }
+
+        m_rx = std::fabs(m_rx);
+        m_ry = std::fabs(m_ry);
+
+        if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )
+        {
+            m_rx = L2(d) / 2;
+            m_ry = 0;
+            m_rot_angle = std::atan2(d[Y], d[X]);
+            if (m_rot_angle < 0) m_rot_angle += 2*M_PI;
+            m_start_angle = 0;
+            m_end_angle = M_PI;
+            m_center = middle_point(initialPoint(), finalPoint());
+            m_large_arc = false;
+            m_sweep = false;
+            return;
+        }
+    }
+    else
+    {
+        if ( are_near(initialPoint(), finalPoint()) )
+        {
+            if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
+            {
+                m_start_angle = m_end_angle = 0;
+                m_center = initialPoint();
+                return;
+            }
+            else
+            {
+                THROW_RANGEERROR("initial and final point are the same");
+            }
+        }
+        if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )
+        { // but initialPoint != finalPoint
+            THROW_RANGEERROR(
+                "there is no ellipse that satisfies the given constraints: "
+                "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"
+            );
+        }
+        if ( are_near(ray(Y), 0) )
+        {
+            Point v = initialPoint() - finalPoint();
+            if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )
+            {
+                double angle = std::atan2(v[Y], v[X]);
+                if (angle < 0) angle += 2*M_PI;
+                if ( are_near( angle, rotation_angle() ) )
+                {
+                    m_start_angle = 0;
+                    m_end_angle = M_PI;
+                    m_center = v/2 + finalPoint();
+                    return;
+                }
+                angle -= M_PI;
+                if ( angle < 0 ) angle += 2*M_PI;
+                if ( are_near( angle, rotation_angle() ) )
+                {
+                    m_start_angle = M_PI;
+                    m_end_angle = 0;
+                    m_center = v/2 + finalPoint();
+                    return;
+                }
+                THROW_RANGEERROR(
+                    "there is no ellipse that satisfies the given constraints: "
+                    "ray(Y) == 0 "
+                    "and slope(initialPoint - finalPoint) != rotation_angle "
+                    "and != rotation_angle + PI"
+                );
+            }
+            if ( L2sq(v) > 4*ray(X)*ray(X) )
+            {
+                THROW_RANGEERROR(
+                    "there is no ellipse that satisfies the given constraints: "
+                    "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"
+                );
+            }
+            else
+            {
+                THROW_RANGEERROR(
+                    "there is infinite ellipses that satisfy the given constraints: "
+                    "ray(Y) == 0  and distance(initialPoint, finalPoint) < 2*ray(X)"
+                );
+            }
+
+        }
+
+        if ( are_near(ray(X), 0) )
+        {
+            Point v = initialPoint() - finalPoint();
+            if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )
+            {
+                double angle = std::atan2(v[Y], v[X]);
+                if (angle < 0) angle += 2*M_PI;
+                double rot_angle = rotation_angle() + M_PI/2;
+                if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;
+                if ( are_near( angle, rot_angle ) )
+                {
+                    m_start_angle = M_PI/2;
+                    m_end_angle = 3*M_PI/2;
+                    m_center = v/2 + finalPoint();
+                    return;
+                }
+                angle -= M_PI;
+                if ( angle < 0 ) angle += 2*M_PI;
+                if ( are_near( angle, rot_angle ) )
+                {
+                    m_start_angle = 3*M_PI/2;
+                    m_end_angle = M_PI/2;
+                    m_center = v/2 + finalPoint();
+                    return;
+                }
+                THROW_RANGEERROR(
+                    "there is no ellipse that satisfies the given constraints: "
+                    "ray(X) == 0 "
+                    "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "
+                    "and != rotation_angle + (3/2)*PI"
+                );
+            }
+            if ( L2sq(v) > 4*ray(Y)*ray(Y) )
+            {
+                THROW_RANGEERROR(
+                    "there is no ellipse that satisfies the given constraints: "
+                    "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"
+                );
+            }
+            else
+            {
+                THROW_RANGEERROR(
+                    "there is infinite ellipses that satisfy the given constraints: "
+                    "ray(X) == 0  and distance(initialPoint, finalPoint) < 2*ray(Y)"
+                );
+            }
+
+        }
+
+    }
+
+    double sin_rot_angle = std::sin(rotation_angle());
+    double cos_rot_angle = std::cos(rotation_angle());
+
+
+    Matrix m( cos_rot_angle, -sin_rot_angle,
+              sin_rot_angle, cos_rot_angle,
+              0,             0              );
+
+    Point p = (d / 2) * m;
+    double rx2 = m_rx * m_rx;
+    double ry2 = m_ry * m_ry;
+    double rxpy = m_rx * p[Y];
+    double rypx = m_ry * p[X];
+    double rx2py2 = rxpy * rxpy;
+    double ry2px2 = rypx * rypx;
+    double num = rx2 * ry2;
+    double den = rx2py2 + ry2px2;
+    assert(den != 0);
+    double rad = num / den;
+    Point c(0,0);
+    if (rad > 1)
+    {
+        rad -= 1;
+        rad = std::sqrt(rad);
+
+        if (m_large_arc == m_sweep) rad = -rad;
+        c = rad * Point(rxpy / m_ry, -rypx / m_rx);
+
+        m[1] = -m[1];
+        m[2] = -m[2];
+
+        m_center = c * m + middle_point(initialPoint(), finalPoint());
+    }
+    else if (rad == 1 || is_svg_compliant())
+    {
+        double lamda = std::sqrt(1 / rad);
+        m_rx *= lamda;
+        m_ry *= lamda;
+        m_center = middle_point(initialPoint(), finalPoint());
+    }
+    else
+    {
+        THROW_RANGEERROR(
+            "there is no ellipse that satisfies the given constraints"
+        );
+    }
+
+    Point sp((p[X] - c[X]) / m_rx, (p[Y] - c[Y]) / m_ry);
+    Point ep((-p[X] - c[X]) / m_rx, (-p[Y] - c[Y]) / m_ry);
+    Point v(1, 0);
+    m_start_angle = angle_between(v, sp);
+    double sweep_angle = angle_between(sp, ep);
+    if (!m_sweep && sweep_angle > 0) sweep_angle -= 2*M_PI;
+    if (m_sweep && sweep_angle < 0) sweep_angle += 2*M_PI;
+
+    if (m_start_angle < 0) m_start_angle += 2*M_PI;
+    m_end_angle = m_start_angle + sweep_angle;
+    if (m_end_angle < 0) m_end_angle += 2*M_PI;
+    if (m_end_angle >= 2*M_PI) m_end_angle -= 2*M_PI;
+}
+
+
+D2<SBasis> SVGEllipticalArc::toSBasis() const
+{
+
+    if (isDegenerate() && is_svg_compliant())
+        return chord().toSBasis();
+
+    D2<SBasis> arc;
+    // the interval of parametrization has to be [0,1]
+    Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );
+    Linear param(start_angle(), et);
+    Coord cos_rot_angle = std::cos(rotation_angle());
+    Coord sin_rot_angle = std::sin(rotation_angle());
+    // order = 4 seems to be enough to get a perfect looking elliptical arc
+    // should it be choosen in function of the arc length anyway ?
+    // or maybe a user settable parameter: toSBasis(unsigned int order) ?
+    SBasis arc_x = ray(X) * cos(param,4);
+    SBasis arc_y = ray(Y) * sin(param,4);
+    arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));
+    arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));
+    return arc;
+}
+
+
+Coord SVGEllipticalArc::map_to_02PI(Coord t) const
+{
+    if ( sweep_flag() )
+    {
+        Coord angle = start_angle() + sweep_angle() * t;
+        if ( !(angle < 2*M_PI) )
+            angle -= 2*M_PI;
+        return angle;
+    }
+    else
+    {
+        Coord angle = start_angle() - sweep_angle() * t;
+        if ( angle < 0 ) angle += 2*M_PI;
+        return angle;
+    }
+}
+
+Coord SVGEllipticalArc::map_to_01(Coord angle) const
+{
+    return map_circular_arc_on_unit_interval(angle, start_angle(),
+                                             end_angle(), sweep_flag());
+}
+
+
+namespace detail
+{
+
+struct ellipse_equation
+{
+    ellipse_equation(double a, double b, double c, double d, double e, double f)
+        : A(a), B(b), C(c), D(d), E(e), F(f)
+    {
+    }
+
+    double operator()(double x, double y) const
+    {
+        // A * x * x + B * x * y + C * y * y + D * x + E * y + F;
+        return (A * x + B * y + D) * x + (C * y + E) * y + F;
+    }
+
+    double operator()(Point const& p) const
+    {
+        return (*this)(p[X], p[Y]);
+    }
+
+    Point normal(double x, double y) const
+    {
+        Point n( 2 * A * x + B * y + D, 2 * C * y + B * x + E );
+        return unit_vector(n);
+    }
+
+    Point normal(Point const& p) const
+    {
+        return normal(p[X], p[Y]);
+    }
+
+    double A, B, C, D, E, F;
+};
+
+}
+
+make_elliptical_arc::
+make_elliptical_arc( SVGEllipticalArc& _ea,
+                     curve_type const& _curve,
+                     unsigned int _total_samples,
+                     double _tolerance )
+    : ea(_ea), curve(_curve),
+      dcurve( unitVector(derivative(curve)) ),
+      model(), fitter(model, _total_samples),
+      tolerance(_tolerance), tol_at_extr(tolerance/2),
+      tol_at_center(0.1), angle_tol(0.1),
+      initial_point(curve.at0()), final_point(curve.at1()),
+      N(_total_samples), last(N-1), partitions(N-1), p(N),
+      svg_compliant(true)
+{
+}
+
+bool
+make_elliptical_arc::
+bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,
+                double e1x, double e1y, double e2 )
+{
+    dist_err = std::fabs( ee(p[k]) );
+    dist_bound = std::fabs( e1x * p[k][X] + e1y * p[k][Y] + e2 );
+    angle_err = std::fabs( dot( dcurve(k/partitions), ee.normal(p[k]) ) );
+    //angle_err *= angle_err;
+    return ( dist_err  > dist_bound || angle_err > angle_tol );
+}
+
+bool
+make_elliptical_arc::
+check_bound(double A, double B, double C, double D, double E, double F)
+{
+    // check error magnitude
+    detail::ellipse_equation ee(A, B, C, D, E, F);
+
+    double e1x = (2*A + B) * tol_at_extr;
+    double e1y = (B + 2*C) * tol_at_extr;
+    double e2 = ((D + E)  + (A + B + C) * tol_at_extr) * tol_at_extr;
+    if ( bound_exceeded(0, ee, e1x, e1y, e2) )
+    {
+        print_bound_error(0);
+        return false;
+    }
+    if ( bound_exceeded(0, ee, e1x, e1y, e2) )
+    {
+        print_bound_error(last);
+        return false;
+    }
+
+    e1x = (2*A + B) * tolerance;
+    e1y = (B + 2*C) * tolerance;
+    e2 = ((D + E)  + (A + B + C) * tolerance) * tolerance;
+//  std::cerr << "e1x = " << e1x << std::endl;
+//  std::cerr << "e1y = " << e1y << std::endl;
+//  std::cerr << "e2 = " << e2 << std::endl;
+
+    for ( unsigned int k = 1; k < last; ++k )
+    {
+        if ( bound_exceeded(k, ee, e1x, e1y, e2) )
+        {
+            print_bound_error(k);
+            return false;
+        }
+    }
+
+    return true;
+}
+
+void make_elliptical_arc::fit()
+{
+    for (unsigned int k = 0; k < N; ++k)
+    {
+        p[k] = curve( k / partitions );
+        fitter.append(p[k]);
+    }
+    fitter.update();
+
+    NL::Vector z(N, 0.0);
+    fitter.result(z);
+}
+
+bool make_elliptical_arc::make_elliptiarc()
+{
+    const NL::Vector & coeff = fitter.result();
+    Ellipse e;
+    try
+    {
+        e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);
+    }
+    catch(LogicalError exc)
+    {
+        return false;
+    }
+
+    Point inner_point = curve(0.5);
+
+    if (svg_compliant_flag())
+    {
+        ea = e.arc(initial_point, inner_point, final_point);
+    }
+    else
+    {
+        try
+        {
+            ea = e.arc(initial_point, inner_point, final_point, false);
+        }
+        catch(RangeError exc)
+        {
+            return false;
+        }
+    }
+
+
+    if ( !are_near( e.center(),
+                    ea.center(),
+                    tol_at_center * std::min(e.ray(X),e.ray(Y))
+                  )
+       )
+    {
+        return false;
+    }
+    return true;
+}
+
+
+
+} // end namespace Geom
+
+
+
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+
diff --git a/src/2geom/svg-elliptical-arc.h b/src/2geom/svg-elliptical-arc.h
index ae6d2e254..f129e5a65 100644
--- a/src/2geom/svg-elliptical-arc.h
+++ b/src/2geom/svg-elliptical-arc.h
@@ -1,435 +1,435 @@
-/*

- * Elliptical Arc - implementation of the svg elliptical arc path element

- *

- * Authors:

- *      MenTaLguY <mental@rydia.net>

- *      Marco Cecchetti <mrcekets at gmail.com>

- *

- * Copyright 2007-2008  authors

- *

- * This library is free software; you can redistribute it and/or

- * modify it either under the terms of the GNU Lesser General Public

- * License version 2.1 as published by the Free Software Foundation

- * (the "LGPL") or, at your option, under the terms of the Mozilla

- * Public License Version 1.1 (the "MPL"). If you do not alter this

- * notice, a recipient may use your version of this file under either

- * the MPL or the LGPL.

- *

- * You should have received a copy of the LGPL along with this library

- * in the file COPYING-LGPL-2.1; if not, write to the Free Software

- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

- * You should have received a copy of the MPL along with this library

- * in the file COPYING-MPL-1.1

- *

- * The contents of this file are subject to the Mozilla Public License

- * Version 1.1 (the "License"); you may not use this file except in

- * compliance with the License. You may obtain a copy of the License at

- * http://www.mozilla.org/MPL/

- *

- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY

- * OF ANY KIND, either express or implied. See the LGPL or the MPL for

- * the specific language governing rights and limitations.

- */

-

-

-#ifndef _2GEOM_SVG_ELLIPTICAL_ARC_H_

-#define _2GEOM_SVG_ELLIPTICAL_ARC_H_

-

-

-#include <2geom/curve.h>

-#include <2geom/angle.h>

-#include <2geom/utils.h>

-#include <2geom/bezier-curve.h>

-#include <2geom/sbasis-curve.h>  // for non-native methods

-#include <2geom/numeric/vector.h>

-#include <2geom/numeric/fitting-tool.h>

-#include <2geom/numeric/fitting-model.h>

-

-

-#include <algorithm>

-

-

-

-namespace Geom

-{

-

-class SVGEllipticalArc : public Curve

-{

-  public:

-    SVGEllipticalArc(bool _svg_compliant = true)

-        : m_initial_point(Point(0,0)), m_final_point(Point(0,0)),

-          m_rx(0), m_ry(0), m_rot_angle(0),

-          m_large_arc(true), m_sweep(true),

-          m_svg_compliant(_svg_compliant)

-    {

-        m_start_angle = m_end_angle = 0;

-        m_center = Point(0,0);

-    }

-

-    SVGEllipticalArc( Point _initial_point, double _rx, double _ry,

-                      double _rot_angle, bool _large_arc, bool _sweep,

-                      Point _final_point,

-                      bool _svg_compliant = true

-                    )

-        : m_initial_point(_initial_point), m_final_point(_final_point),

-          m_rx(_rx), m_ry(_ry), m_rot_angle(_rot_angle),

-          m_large_arc(_large_arc), m_sweep(_sweep),

-          m_svg_compliant(_svg_compliant)

-    {

-            calculate_center_and_extreme_angles();

-    }

-

-    void set( Point _initial_point, double _rx, double _ry,

-              double _rot_angle, bool _large_arc, bool _sweep,

-              Point _final_point

-             )

-    {

-        m_initial_point = _initial_point;

-        m_final_point = _final_point;

-        m_rx = _rx;

-        m_ry = _ry;

-        m_rot_angle = _rot_angle;

-        m_large_arc = _large_arc;

-        m_sweep = _sweep;

-        calculate_center_and_extreme_angles();

-    }

-

-    Curve* duplicate() const

-    {

-        return new SVGEllipticalArc(*this);

-    }

-

-    double center(unsigned int i) const

-    {

-        return m_center[i];

-    }

-

-    Point center() const

-    {

-        return m_center;

-    }

-

-    Point initialPoint() const

-    {

-        return m_initial_point;

-    }

-

-    Point finalPoint() const

-    {

-        return m_final_point;

-    }

-

-    double start_angle() const

-    {

-        return m_start_angle;

-    }

-

-    double end_angle() const

-    {

-        return m_end_angle;

-    }

-

-    double ray(unsigned int i) const

-    {

-        return (i == 0) ? m_rx : m_ry;

-    }

-

-    bool large_arc_flag() const

-    {

-        return m_large_arc;

-    }

-

-    bool sweep_flag() const

-    {

-        return m_sweep;

-    }

-

-    double rotation_angle() const

-    {

-        return m_rot_angle;

-    }

-

-    void setInitial( const Point _point)

-    {

-        m_initial_point = _point;

-        calculate_center_and_extreme_angles();

-    }

-

-    void setFinal( const Point _point)

-    {

-        m_final_point = _point;

-        calculate_center_and_extreme_angles();

-    }

-

-    void setExtremes( const Point& _initial_point, const Point& _final_point )

-    {

-        m_initial_point = _initial_point;

-        m_final_point = _final_point;

-        calculate_center_and_extreme_angles();

-    }

-

-    bool isDegenerate() const

-    {

-        return ( are_near(ray(X), 0) || are_near(ray(Y), 0) );

-    }

-

-    bool is_svg_compliant() const

-    {

-        return m_svg_compliant;

-    }

-

-    Rect boundsFast() const

-    {

-        return boundsExact();

-    }

-

-    Rect boundsExact() const;

-

-    // TODO: native implementation of the following methods

-    Rect boundsLocal(Interval i, unsigned int deg) const

-    {

-        if (isDegenerate() && is_svg_compliant())

-            return chord().boundsLocal(i, deg);

-        else

-            return SBasisCurve(toSBasis()).boundsLocal(i, deg);

-    }

-

-    std::vector<double> roots(double v, Dim2 d) const;

-

-    std::vector<double>

-    allNearestPoints( Point const& p, double from = 0, double to = 1 ) const;

-

-    double nearestPoint( Point const& p, double from = 0, double to = 1 ) const

-    {

-        if ( are_near(ray(X), ray(Y)) && are_near(center(), p) )

-        {

-            return from;

-        }

-        return allNearestPoints(p, from, to).front();

-    }

-

-    // TODO: native implementation of the following methods

-    int winding(Point p) const

-    {

-        if (isDegenerate() && is_svg_compliant())

-            return chord().winding(p);

-        else

-            return SBasisCurve(toSBasis()).winding(p);

-    }

-

-    Curve *derivative() const;

-

-    // TODO: native implementation of the following methods

-    Curve *transformed(Matrix const &m) const

-    {

-        return SBasisCurve(toSBasis()).transformed(m);

-    }

-

-    std::vector<Point> pointAndDerivatives(Coord t, unsigned int n) const;

-

-    D2<SBasis> toSBasis() const;

-

-    bool containsAngle(Coord angle) const;

-

-    double valueAtAngle(Coord t, Dim2 d) const;

-

-    Point pointAtAngle(Coord t) const

-    {

-        double sin_rot_angle = std::sin(rotation_angle());

-        double cos_rot_angle = std::cos(rotation_angle());

-        Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,

-                 -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,

-                  center(X),              center(Y) );

-        Point p( std::cos(t), std::sin(t) );

-        return p * m;

-    }

-

-    double valueAt(Coord t, Dim2 d) const

-    {

-        if (isDegenerate() && is_svg_compliant())

-            return chord().valueAt(t, d);

-

-        Coord tt = map_to_02PI(t);

-        return valueAtAngle(tt, d);

-    }

-

-    Point pointAt(Coord t) const

-    {

-        if (isDegenerate() && is_svg_compliant())

-            return chord().pointAt(t);

-

-        Coord tt = map_to_02PI(t);

-        return pointAtAngle(tt);

-    }

-

-    std::pair<SVGEllipticalArc, SVGEllipticalArc>

-    subdivide(Coord t) const

-    {

-        SVGEllipticalArc* arc1 = static_cast<SVGEllipticalArc*>(portion(0, t));

-        SVGEllipticalArc* arc2 = static_cast<SVGEllipticalArc*>(portion(t, 1));

-        assert( arc1 != NULL && arc2 != NULL);

-        std::pair<SVGEllipticalArc, SVGEllipticalArc> arc_pair(*arc1, *arc2);

-        delete arc1;

-        delete arc2;

-        return arc_pair;

-    }

-

-    Curve* portion(double f, double t) const;

-

-    // the arc is the same but traversed in the opposite direction

-    Curve* reverse() const

-    {

-        SVGEllipticalArc* rarc = new SVGEllipticalArc( *this );

-        rarc->m_sweep = !m_sweep;

-        rarc->m_initial_point = m_final_point;

-        rarc->m_final_point = m_initial_point;

-        rarc->m_start_angle = m_end_angle;

-        rarc->m_end_angle = m_start_angle;

-        return rarc;

-    }

-

-    double sweep_angle() const

-    {

-        Coord d = end_angle() - start_angle();

-        if ( !sweep_flag() ) d = -d;

-        if ( d < 0 )

-            d += 2*M_PI;

-        return d;

-    }

-

-    LineSegment chord() const

-    {

-        return LineSegment(initialPoint(), finalPoint());

-    }

-

-  private:

-    Coord map_to_02PI(Coord t) const;

-    Coord map_to_01(Coord angle) const;

-    void calculate_center_and_extreme_angles();

-

-  private:

-    Point m_initial_point, m_final_point;

-    double m_rx, m_ry, m_rot_angle;

-    bool m_large_arc, m_sweep;

-    double m_start_angle, m_end_angle;

-    Point m_center;

-    bool m_svg_compliant;

-

-}; // end class SVGEllipticalArc

-

-template< class charT >

-inline

-std::basic_ostream<charT> &

-operator<< (std::basic_ostream<charT> & os, const SVGEllipticalArc & ea)

-{

-    os << "{ cx: " << ea.center(X) << ", cy: " <<  ea.center(Y)

-       << ", rx: " << ea.ray(X) << ", ry: " << ea.ray(Y)

-       << ", rot angle: " << decimal_round(rad_to_deg(ea.rotation_angle()),2)

-       << ", start angle: " << decimal_round(rad_to_deg(ea.start_angle()),2)

-       << ", end angle: " << decimal_round(rad_to_deg(ea.end_angle()),2)

-       << " }";

-

-    return os;

-}

-

-

-

-

-namespace detail

-{

-    struct ellipse_equation;

-}

-

-

-class make_elliptical_arc

-{

-  public:

-    typedef D2<SBasis> curve_type;

-

-    make_elliptical_arc( SVGEllipticalArc& _ea,

-                         curve_type const& _curve,

-                         unsigned int _total_samples,

-                         double _tolerance );

-

-  private:

-    bool bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,

-                         double e1x, double e1y, double e2 );

-

-    bool check_bound(double A, double B, double C, double D, double E, double F);

-

-    void fit();

-

-    bool make_elliptiarc();

-

-    void print_bound_error(unsigned int k)

-    {

-        std::cerr

-            << "tolerance error" << std::endl

-            << "at point: " << k << std::endl

-            << "error value: "<< dist_err << std::endl

-            << "bound: " << dist_bound << std::endl

-            << "angle error: " << angle_err

-            << " (" << angle_tol << ")" << std::endl;

-    }

-

-  public:

-    bool operator()()

-    {

-        const NL::Vector & coeff = fitter.result();

-        fit();

-        if ( !check_bound(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]) )

-            return false;

-        if ( !(make_elliptiarc()) ) return false;

-        return true;

-    }

-

-    bool svg_compliant_flag() const

-    {

-        return svg_compliant;

-    }

-

-    void svg_compliant_on()

-    {

-        svg_compliant = true;

-    }

-

-    void svg_compliant_off()

-    {

-        svg_compliant = false;

-    }

-

-  private:

-      SVGEllipticalArc& ea;

-      const curve_type & curve;

-      Piecewise<D2<SBasis> > dcurve;

-      NL::LFMEllipse model;

-      NL::least_squeares_fitter<NL::LFMEllipse> fitter;

-      double tolerance, tol_at_extr, tol_at_center, angle_tol;

-      Point initial_point, final_point;

-      unsigned int N;

-      unsigned int last; // N-1

-      double partitions; // N-1

-      std::vector<Point> p; // sample points

-      double dist_err, dist_bound, angle_err;

-      bool svg_compliant;

-};

-

-

-} // end namespace Geom

-

-

-

-

-#endif /* _2GEOM_SVG_ELLIPTICAL_ARC_H_ */

-

-/*

-  Local Variables:

-  mode:c++

-  c-file-style:"stroustrup"

-  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))

-  indent-tabs-mode:nil

-  fill-column:99

-  End:

-*/

-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :

-

+/*
+ * Elliptical Arc - implementation of the svg elliptical arc path element
+ *
+ * Authors:
+ *      MenTaLguY <mental@rydia.net>
+ *      Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2007-2008  authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _2GEOM_SVG_ELLIPTICAL_ARC_H_
+#define _2GEOM_SVG_ELLIPTICAL_ARC_H_
+
+
+#include <2geom/curve.h>
+#include <2geom/angle.h>
+#include <2geom/utils.h>
+#include <2geom/bezier-curve.h>
+#include <2geom/sbasis-curve.h>  // for non-native methods
+#include <2geom/numeric/vector.h>
+#include <2geom/numeric/fitting-tool.h>
+#include <2geom/numeric/fitting-model.h>
+
+
+#include <algorithm>
+
+
+
+namespace Geom
+{
+
+class SVGEllipticalArc : public Curve
+{
+  public:
+    SVGEllipticalArc(bool _svg_compliant = true)
+        : m_initial_point(Point(0,0)), m_final_point(Point(0,0)),
+          m_rx(0), m_ry(0), m_rot_angle(0),
+          m_large_arc(true), m_sweep(true),
+          m_svg_compliant(_svg_compliant)
+    {
+        m_start_angle = m_end_angle = 0;
+        m_center = Point(0,0);
+    }
+
+    SVGEllipticalArc( Point _initial_point, double _rx, double _ry,
+                      double _rot_angle, bool _large_arc, bool _sweep,
+                      Point _final_point,
+                      bool _svg_compliant = true
+                    )
+        : m_initial_point(_initial_point), m_final_point(_final_point),
+          m_rx(_rx), m_ry(_ry), m_rot_angle(_rot_angle),
+          m_large_arc(_large_arc), m_sweep(_sweep),
+          m_svg_compliant(_svg_compliant)
+    {
+            calculate_center_and_extreme_angles();
+    }
+
+    void set( Point _initial_point, double _rx, double _ry,
+              double _rot_angle, bool _large_arc, bool _sweep,
+              Point _final_point
+             )
+    {
+        m_initial_point = _initial_point;
+        m_final_point = _final_point;
+        m_rx = _rx;
+        m_ry = _ry;
+        m_rot_angle = _rot_angle;
+        m_large_arc = _large_arc;
+        m_sweep = _sweep;
+        calculate_center_and_extreme_angles();
+    }
+
+    Curve* duplicate() const
+    {
+        return new SVGEllipticalArc(*this);
+    }
+
+    double center(unsigned int i) const
+    {
+        return m_center[i];
+    }
+
+    Point center() const
+    {
+        return m_center;
+    }
+
+    Point initialPoint() const
+    {
+        return m_initial_point;
+    }
+
+    Point finalPoint() const
+    {
+        return m_final_point;
+    }
+
+    double start_angle() const
+    {
+        return m_start_angle;
+    }
+
+    double end_angle() const
+    {
+        return m_end_angle;
+    }
+
+    double ray(unsigned int i) const
+    {
+        return (i == 0) ? m_rx : m_ry;
+    }
+
+    bool large_arc_flag() const
+    {
+        return m_large_arc;
+    }
+
+    bool sweep_flag() const
+    {
+        return m_sweep;
+    }
+
+    double rotation_angle() const
+    {
+        return m_rot_angle;
+    }
+
+    void setInitial( const Point _point)
+    {
+        m_initial_point = _point;
+        calculate_center_and_extreme_angles();
+    }
+
+    void setFinal( const Point _point)
+    {
+        m_final_point = _point;
+        calculate_center_and_extreme_angles();
+    }
+
+    void setExtremes( const Point& _initial_point, const Point& _final_point )
+    {
+        m_initial_point = _initial_point;
+        m_final_point = _final_point;
+        calculate_center_and_extreme_angles();
+    }
+
+    bool isDegenerate() const
+    {
+        return ( are_near(ray(X), 0) || are_near(ray(Y), 0) );
+    }
+
+    bool is_svg_compliant() const
+    {
+        return m_svg_compliant;
+    }
+
+    Rect boundsFast() const
+    {
+        return boundsExact();
+    }
+
+    Rect boundsExact() const;
+
+    // TODO: native implementation of the following methods
+    Rect boundsLocal(Interval i, unsigned int deg) const
+    {
+        if (isDegenerate() && is_svg_compliant())
+            return chord().boundsLocal(i, deg);
+        else
+            return SBasisCurve(toSBasis()).boundsLocal(i, deg);
+    }
+
+    std::vector<double> roots(double v, Dim2 d) const;
+
+    std::vector<double>
+    allNearestPoints( Point const& p, double from = 0, double to = 1 ) const;
+
+    double nearestPoint( Point const& p, double from = 0, double to = 1 ) const
+    {
+        if ( are_near(ray(X), ray(Y)) && are_near(center(), p) )
+        {
+            return from;
+        }
+        return allNearestPoints(p, from, to).front();
+    }
+
+    // TODO: native implementation of the following methods
+    int winding(Point p) const
+    {
+        if (isDegenerate() && is_svg_compliant())
+            return chord().winding(p);
+        else
+            return SBasisCurve(toSBasis()).winding(p);
+    }
+
+    Curve *derivative() const;
+
+    // TODO: native implementation of the following methods
+    Curve *transformed(Matrix const &m) const
+    {
+        return SBasisCurve(toSBasis()).transformed(m);
+    }
+
+    std::vector<Point> pointAndDerivatives(Coord t, unsigned int n) const;
+
+    D2<SBasis> toSBasis() const;
+
+    bool containsAngle(Coord angle) const;
+
+    double valueAtAngle(Coord t, Dim2 d) const;
+
+    Point pointAtAngle(Coord t) const
+    {
+        double sin_rot_angle = std::sin(rotation_angle());
+        double cos_rot_angle = std::cos(rotation_angle());
+        Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,
+                 -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,
+                  center(X),              center(Y) );
+        Point p( std::cos(t), std::sin(t) );
+        return p * m;
+    }
+
+    double valueAt(Coord t, Dim2 d) const
+    {
+        if (isDegenerate() && is_svg_compliant())
+            return chord().valueAt(t, d);
+
+        Coord tt = map_to_02PI(t);
+        return valueAtAngle(tt, d);
+    }
+
+    Point pointAt(Coord t) const
+    {
+        if (isDegenerate() && is_svg_compliant())
+            return chord().pointAt(t);
+
+        Coord tt = map_to_02PI(t);
+        return pointAtAngle(tt);
+    }
+
+    std::pair<SVGEllipticalArc, SVGEllipticalArc>
+    subdivide(Coord t) const
+    {
+        SVGEllipticalArc* arc1 = static_cast<SVGEllipticalArc*>(portion(0, t));
+        SVGEllipticalArc* arc2 = static_cast<SVGEllipticalArc*>(portion(t, 1));
+        assert( arc1 != NULL && arc2 != NULL);
+        std::pair<SVGEllipticalArc, SVGEllipticalArc> arc_pair(*arc1, *arc2);
+        delete arc1;
+        delete arc2;
+        return arc_pair;
+    }
+
+    Curve* portion(double f, double t) const;
+
+    // the arc is the same but traversed in the opposite direction
+    Curve* reverse() const
+    {
+        SVGEllipticalArc* rarc = new SVGEllipticalArc( *this );
+        rarc->m_sweep = !m_sweep;
+        rarc->m_initial_point = m_final_point;
+        rarc->m_final_point = m_initial_point;
+        rarc->m_start_angle = m_end_angle;
+        rarc->m_end_angle = m_start_angle;
+        return rarc;
+    }
+
+    double sweep_angle() const
+    {
+        Coord d = end_angle() - start_angle();
+        if ( !sweep_flag() ) d = -d;
+        if ( d < 0 )
+            d += 2*M_PI;
+        return d;
+    }
+
+    LineSegment chord() const
+    {
+        return LineSegment(initialPoint(), finalPoint());
+    }
+
+  private:
+    Coord map_to_02PI(Coord t) const;
+    Coord map_to_01(Coord angle) const;
+    void calculate_center_and_extreme_angles();
+
+  private:
+    Point m_initial_point, m_final_point;
+    double m_rx, m_ry, m_rot_angle;
+    bool m_large_arc, m_sweep;
+    double m_start_angle, m_end_angle;
+    Point m_center;
+    bool m_svg_compliant;
+
+}; // end class SVGEllipticalArc
+
+template< class charT >
+inline
+std::basic_ostream<charT> &
+operator<< (std::basic_ostream<charT> & os, const SVGEllipticalArc & ea)
+{
+    os << "{ cx: " << ea.center(X) << ", cy: " <<  ea.center(Y)
+       << ", rx: " << ea.ray(X) << ", ry: " << ea.ray(Y)
+       << ", rot angle: " << decimal_round(rad_to_deg(ea.rotation_angle()),2)
+       << ", start angle: " << decimal_round(rad_to_deg(ea.start_angle()),2)
+       << ", end angle: " << decimal_round(rad_to_deg(ea.end_angle()),2)
+       << " }";
+
+    return os;
+}
+
+
+
+
+namespace detail
+{
+    struct ellipse_equation;
+}
+
+
+class make_elliptical_arc
+{
+  public:
+    typedef D2<SBasis> curve_type;
+
+    make_elliptical_arc( SVGEllipticalArc& _ea,
+                         curve_type const& _curve,
+                         unsigned int _total_samples,
+                         double _tolerance );
+
+  private:
+    bool bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,
+                         double e1x, double e1y, double e2 );
+
+    bool check_bound(double A, double B, double C, double D, double E, double F);
+
+    void fit();
+
+    bool make_elliptiarc();
+
+    void print_bound_error(unsigned int k)
+    {
+        std::cerr
+            << "tolerance error" << std::endl
+            << "at point: " << k << std::endl
+            << "error value: "<< dist_err << std::endl
+            << "bound: " << dist_bound << std::endl
+            << "angle error: " << angle_err
+            << " (" << angle_tol << ")" << std::endl;
+    }
+
+  public:
+    bool operator()()
+    {
+        const NL::Vector & coeff = fitter.result();
+        fit();
+        if ( !check_bound(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]) )
+            return false;
+        if ( !(make_elliptiarc()) ) return false;
+        return true;
+    }
+
+    bool svg_compliant_flag() const
+    {
+        return svg_compliant;
+    }
+
+    void svg_compliant_on()
+    {
+        svg_compliant = true;
+    }
+
+    void svg_compliant_off()
+    {
+        svg_compliant = false;
+    }
+
+  private:
+      SVGEllipticalArc& ea;
+      const curve_type & curve;
+      Piecewise<D2<SBasis> > dcurve;
+      NL::LFMEllipse model;
+      NL::least_squeares_fitter<NL::LFMEllipse> fitter;
+      double tolerance, tol_at_extr, tol_at_center, angle_tol;
+      Point initial_point, final_point;
+      unsigned int N;
+      unsigned int last; // N-1
+      double partitions; // N-1
+      std::vector<Point> p; // sample points
+      double dist_err, dist_bound, angle_err;
+      bool svg_compliant;
+};
+
+
+} // end namespace Geom
+
+
+
+
+#endif /* _2GEOM_SVG_ELLIPTICAL_ARC_H_ */
+
+/*
+  Local Variables:
+  mode:c++
+  c-file-style:"stroustrup"
+  c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+  indent-tabs-mode:nil
+  fill-column:99
+  End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
+