![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
summary | shortlog | log | commit | commitdiff | tree
raw | patch | inline | side by side (parent: 1e12c38)
raw | patch | inline | side by side (parent: 1e12c38)
| author | joncruz <joncruz@users.sourceforge.net> | |
| Fri, 4 Jul 2008 06:38:30 +0000 (06:38 +0000) | ||
| committer | joncruz <joncruz@users.sourceforge.net> | |
| Fri, 4 Jul 2008 06:38:30 +0000 (06:38 +0000) |
diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp
index c9c5b9ec4c03393dee2029cb158400b302e7831b..20ac0bb2bc358e54c745571b61a4fcea06384164 100644 (file)
--- a/src/2geom/ellipse.cpp
+++ b/src/2geom/ellipse.cpp
-/*\r
- * Ellipse Curve\r
- *\r
- * Authors:\r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- *\r
- * Copyright 2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#include <2geom/ellipse.h>\r
-#include <2geom/svg-elliptical-arc.h>\r
-#include <2geom/numeric/fitting-tool.h>\r
-#include <2geom/numeric/fitting-model.h>\r
-\r
-\r
-namespace Geom\r
-{\r
-\r
-void Ellipse::set(double A, double B, double C, double D, double E, double F)\r
-{\r
- double den = 4*A*C - B*B;\r
- if ( den == 0 )\r
- {\r
- THROW_LOGICALERROR("den == 0, while computing ellipse centre");\r
- }\r
- m_centre[X] = (B*E - 2*C*D) / den;\r
- m_centre[Y] = (B*D - 2*A*E) / den;\r
-\r
- // evaluate the a coefficient of the ellipse equation in normal form\r
- // E(x,y) = a*(x-cx)^2 + b*(x-cx)*(y-cy) + c*(y-cy)^2 = 1\r
- // where b = a*B , c = a*C, (cx,cy) == centre\r
- double num = A * sqr(m_centre[X])\r
- + B * m_centre[X] * m_centre[Y]\r
- + C * sqr(m_centre[Y])\r
- - A * F;\r
-\r
-\r
- //evaluate ellipse rotation angle\r
- double rot = std::atan2( -B, -(A - C) )/2;\r
-// std::cerr << "rot = " << rot << std::endl;\r
- bool swap_axes = false;\r
- if ( are_near(rot, 0) ) rot = 0;\r
- if ( are_near(rot, M_PI/2) || rot < 0 )\r
- {\r
- swap_axes = true;\r
- }\r
-\r
- // evaluate the length of the ellipse rays\r
- double cosrot = std::cos(rot);\r
- double sinrot = std::sin(rot);\r
- double cos2 = cosrot * cosrot;\r
- double sin2 = sinrot * sinrot;\r
- double cossin = cosrot * sinrot;\r
-\r
- den = A * cos2 + B * cossin + C * sin2;\r
- if ( den == 0 )\r
- {\r
- THROW_LOGICALERROR("den == 0, while computing 'rx' coefficient");\r
- }\r
- double rx2 = num/den;\r
- if ( rx2 < 0 )\r
- {\r
- THROW_LOGICALERROR("rx2 < 0, while computing 'rx' coefficient");\r
- }\r
- double rx = std::sqrt(rx2);\r
-\r
- den = C * cos2 - B * cossin + A * sin2;\r
- if ( den == 0 )\r
- {\r
- THROW_LOGICALERROR("den == 0, while computing 'ry' coefficient");\r
- }\r
- double ry2 = num/den;\r
- if ( ry2 < 0 )\r
- {\r
- THROW_LOGICALERROR("ry2 < 0, while computing 'rx' coefficient");\r
- }\r
- double ry = std::sqrt(ry2);\r
-\r
- // the solution is not unique so we choose always the ellipse\r
- // with a rotation angle between 0 and PI/2\r
- if ( swap_axes ) std::swap(rx, ry);\r
- if ( are_near(rot, M_PI/2)\r
- || are_near(rot, -M_PI/2)\r
- || are_near(rx, ry) )\r
- {\r
- rot = 0;\r
- }\r
- else if ( rot < 0 )\r
- {\r
- rot += M_PI/2;\r
- }\r
-\r
- m_ray[X] = rx;\r
- m_ray[Y] = ry;\r
- m_angle = rot;\r
-}\r
-\r
-\r
-void Ellipse::set(std::vector<Point> const& points)\r
-{\r
- size_t sz = points.size();\r
- if (sz < 5)\r
- {\r
- THROW_RANGEERROR("fitting error: too few points passed");\r
- }\r
- NL::LFMEllipse model;\r
- NL::least_squeares_fitter<NL::LFMEllipse> fitter(model, sz);\r
-\r
- for (size_t i = 0; i < sz; ++i)\r
- {\r
- fitter.append(points[i]);\r
- }\r
- fitter.update();\r
-\r
- NL::Vector z(sz, 0.0);\r
- model.instance(*this, fitter.result(z));\r
-}\r
-\r
-\r
-SVGEllipticalArc\r
-Ellipse::arc(Point const& initial, Point const& inner, Point const& final,\r
- bool _svg_compliant)\r
-{\r
- Point sp_cp = initial - center();\r
- Point ep_cp = final - center();\r
- Point ip_cp = inner - center();\r
-\r
- double angle1 = angle_between(sp_cp, ep_cp);\r
- double angle2 = angle_between(sp_cp, ip_cp);\r
- double angle3 = angle_between(ip_cp, ep_cp);\r
-\r
- bool large_arc_flag = true;\r
- bool sweep_flag = true;\r
-\r
- if ( angle1 > 0 )\r
- {\r
- if ( angle2 > 0 && angle3 > 0 )\r
- {\r
- large_arc_flag = false;\r
- sweep_flag = true;\r
- }\r
- else\r
- {\r
- large_arc_flag = true;\r
- sweep_flag = false;\r
- }\r
- }\r
- else\r
- {\r
- if ( angle2 < 0 && angle3 < 0 )\r
- {\r
- large_arc_flag = false;\r
- sweep_flag = false;\r
- }\r
- else\r
- {\r
- large_arc_flag = true;\r
- sweep_flag = true;\r
- }\r
- }\r
-\r
- SVGEllipticalArc ea( initial, ray(X), ray(Y), rot_angle(),\r
- large_arc_flag, sweep_flag, final, _svg_compliant);\r
- return ea;\r
-}\r
-\r
-\r
-} // end namespace Geom\r
-\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
-\r
-\r
+/*
+ * 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 af8b01e78f99bb85616650d5029bf50ff5c9e694..a10f9ced090bdb5c2a97d9e5329c44948fd9a3c1 100644 (file)
--- a/src/2geom/ellipse.h
+++ b/src/2geom/ellipse.h
-/*\r
- * Ellipse Curve\r
- *\r
- * Authors:\r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- *\r
- * Copyright 2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#ifndef _2GEOM_ELLIPSE_H_\r
-#define _2GEOM_ELLIPSE_H_\r
-\r
-\r
-#include <2geom/point.h>\r
-#include <2geom/exception.h>\r
-\r
-\r
-namespace Geom\r
-{\r
-\r
-class SVGEllipticalArc;\r
-\r
-class Ellipse\r
-{\r
- public:\r
- Ellipse()\r
- {}\r
-\r
- Ellipse(double cx, double cy, double rx, double ry, double a)\r
- : m_centre(cx, cy), m_ray(rx, ry), m_angle(a)\r
- {\r
- }\r
-\r
- // build an ellipse by its implicit equation:\r
- // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0\r
- Ellipse(double A, double B, double C, double D, double E, double F)\r
- {\r
- set(A, B, C, D, E, F);\r
- }\r
-\r
- Ellipse(std::vector<Point> const& points)\r
- {\r
- set(points);\r
- }\r
-\r
- void set(double cx, double cy, double rx, double ry, double a)\r
- {\r
- m_centre[X] = cx;\r
- m_centre[Y] = cy;\r
- m_ray[X] = rx;\r
- m_ray[Y] = ry;\r
- m_angle = a;\r
- }\r
-\r
- // build an ellipse by its implicit equation:\r
- // Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0\r
- void set(double A, double B, double C, double D, double E, double F);\r
-\r
- // biuld up the best fitting ellipse wrt the passed points\r
- // prerequisite: at least 5 points must be passed\r
- void set(std::vector<Point> const& points);\r
-\r
- SVGEllipticalArc\r
- arc(Point const& initial, Point const& inner, Point const& final,\r
- bool _svg_compliant = true);\r
-\r
- Point center() const\r
- {\r
- return m_centre;\r
- }\r
-\r
- Coord center(Dim2 d) const\r
- {\r
- return m_centre[d];\r
- }\r
-\r
- Coord ray(Dim2 d) const\r
- {\r
- return m_ray[d];\r
- }\r
-\r
- Coord rot_angle() const\r
- {\r
- return m_angle;\r
- }\r
-\r
- private:\r
- Point m_centre, m_ray;\r
- double m_angle;\r
-};\r
-\r
-\r
-} // end namespace Geom\r
-\r
-\r
-\r
-#endif // _2GEOM_ELLIPSE_H_\r
-\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+/*
+ * 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 :
index 59d55ba326383f5669f0a65b1aa323a83ebf871b..7fef8c12002fd216b51952ab6b9d7a472a39b82f 100644 (file)
-/*\r
- * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>\r
- *\r
- * Authors:\r
- * \r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- * \r
- * Copyright 2007-2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#include <2geom/nearest-point.h>\r
-\r
-namespace Geom\r
-{\r
-\r
-////////////////////////////////////////////////////////////////////////////////\r
-// D2<SBasis> versions\r
-\r
-/*\r
- * Return the parameter t of a nearest point on the portion of the curve "c", \r
- * related to the interval [from, to], to the point "p".\r
- * The needed curve derivative "dc" is passed as parameter.\r
- * The function return the first nearest point to "p" that is found.\r
- */\r
-\r
-double nearest_point( Point const& p, \r
- D2<SBasis> const& c, \r
- D2<SBasis> const& dc, \r
- double from, double to )\r
-{\r
- if ( from > to ) std::swap(from, to);\r
- if ( from < 0 || to > 1 )\r
- {\r
- THROW_RANGEERROR("[from,to] interval out of bounds");\r
- }\r
-\r
- SBasis dd = dot(c - p, dc); \r
- std::vector<double> zeros = Geom::roots(dd);\r
- \r
- double closest = from;\r
- double min_dist_sq = L2sq(c(from) - p);\r
- double distsq;\r
- for ( unsigned int i = 0; i < zeros.size(); ++i )\r
- {\r
- distsq = L2sq(c(zeros[i]) - p);\r
- if ( min_dist_sq > L2sq(c(zeros[i]) - p) )\r
- {\r
- closest = zeros[i];\r
- min_dist_sq = distsq;\r
- }\r
- }\r
- if ( min_dist_sq > L2sq( c(to) - p ) )\r
- closest = to;\r
- return closest;\r
-\r
-}\r
-\r
-/*\r
- * Return the parameters t of all the nearest points on the portion of \r
- * the curve "c", related to the interval [from, to], to the point "p".\r
- * The needed curve derivative "dc" is passed as parameter.\r
- */\r
-\r
-std::vector<double> \r
-all_nearest_points( Point const& p, \r
- D2<SBasis> const& c, \r
- D2<SBasis> const& /*dc*/, \r
- double from, double to )\r
-{\r
- std::swap(from, to);\r
- if ( from > to ) std::swap(from, to);\r
- if ( from < 0 || to > 1 )\r
- {\r
- THROW_RANGEERROR("[from,to] interval out of bounds");\r
- }\r
-\r
- std::vector<double> result;\r
- SBasis dd = dot(c - p, Geom::derivative(c));\r
- \r
- std::vector<double> zeros = Geom::roots(dd);\r
- std::vector<double> candidates;\r
- candidates.push_back(from);\r
- candidates.insert(candidates.end(), zeros.begin(), zeros.end());\r
- candidates.push_back(to);\r
- std::vector<double> distsq;\r
- distsq.reserve(candidates.size());\r
- for ( unsigned int i = 0; i < candidates.size(); ++i )\r
- {\r
- distsq.push_back( L2sq(c(candidates[i]) - p) );\r
- }\r
- unsigned int closest = 0;\r
- double dsq = distsq[0];\r
- for ( unsigned int i = 1; i < candidates.size(); ++i )\r
- {\r
- if ( dsq > distsq[i] )\r
- {\r
- closest = i;\r
- dsq = distsq[i];\r
- }\r
- }\r
- for ( unsigned int i = 0; i < candidates.size(); ++i )\r
- {\r
- if( distsq[closest] == distsq[i] )\r
- {\r
- result.push_back(candidates[i]);\r
- }\r
- }\r
- return result;\r
-}\r
-\r
-\r
-////////////////////////////////////////////////////////////////////////////////\r
-// Piecewise< D2<SBasis> > versions\r
-\r
-\r
-double nearest_point( Point const& p, \r
- Piecewise< D2<SBasis> > const& c,\r
- double from, double to )\r
-{\r
- if ( from > to ) std::swap(from, to);\r
- if ( from < c.cuts[0] || to > c.cuts[c.size()] )\r
- {\r
- THROW_RANGEERROR("[from,to] interval out of bounds");\r
- }\r
- \r
- unsigned int si = c.segN(from);\r
- unsigned int ei = c.segN(to);\r
- if ( si == ei )\r
- {\r
- double nearest=\r
- nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));\r
- return c.mapToDomain(nearest, si);\r
- }\r
- double t;\r
- double nearest = nearest_point(p, c[si], c.segT(from, si));\r
- unsigned int ni = si;\r
- double dsq;\r
- double mindistsq = distanceSq(p, c[si](nearest));\r
- Rect bb;\r
- for ( unsigned int i = si + 1; i < ei; ++i )\r
- {\r
- bb = bounds_fast(c[i]);\r
- dsq = distanceSq(p, bb);\r
- if ( mindistsq <= dsq ) continue;\r
- t = nearest_point(p, c[i]);\r
- dsq = distanceSq(p, c[i](t));\r
- if ( mindistsq > dsq )\r
- {\r
- nearest = t;\r
- ni = i;\r
- mindistsq = dsq;\r
- }\r
- }\r
- bb = bounds_fast(c[ei]);\r
- dsq = distanceSq(p, bb);\r
- if ( mindistsq > dsq )\r
- {\r
- t = nearest_point(p, c[ei], 0, c.segT(to, ei));\r
- dsq = distanceSq(p, c[ei](t));\r
- if ( mindistsq > dsq )\r
- {\r
- nearest = t;\r
- ni = ei;\r
- }\r
- }\r
- return c.mapToDomain(nearest, ni);\r
-}\r
-\r
-std::vector<double> \r
-all_nearest_points( Point const& p, \r
- Piecewise< D2<SBasis> > const& c, \r
- double from, double to )\r
-{\r
- if ( from > to ) std::swap(from, to);\r
- if ( from < c.cuts[0] || to > c.cuts[c.size()] )\r
- {\r
- THROW_RANGEERROR("[from,to] interval out of bounds");\r
- }\r
- \r
- unsigned int si = c.segN(from);\r
- unsigned int ei = c.segN(to);\r
- if ( si == ei )\r
- {\r
- std::vector<double> all_nearest = \r
- all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));\r
- for ( unsigned int i = 0; i < all_nearest.size(); ++i )\r
- {\r
- all_nearest[i] = c.mapToDomain(all_nearest[i], si);\r
- }\r
- return all_nearest;\r
- }\r
- std::vector<double> all_t;\r
- std::vector< std::vector<double> > all_np;\r
- all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );\r
- std::vector<unsigned int> ni;\r
- ni.push_back(si);\r
- double dsq;\r
- double mindistsq = distanceSq( p, c[si](all_np.front().front()) );\r
- Rect bb;\r
- for ( unsigned int i = si + 1; i < ei; ++i )\r
- {\r
- bb = bounds_fast(c[i]);\r
- dsq = distanceSq(p, bb);\r
- if ( mindistsq < dsq ) continue;\r
- all_t = all_nearest_points(p, c[i]);\r
- dsq = distanceSq( p, c[i](all_t.front()) );\r
- if ( mindistsq > dsq )\r
- {\r
- all_np.clear();\r
- all_np.push_back(all_t);\r
- ni.clear();\r
- ni.push_back(i);\r
- mindistsq = dsq;\r
- }\r
- else if ( mindistsq == dsq )\r
- {\r
- all_np.push_back(all_t);\r
- ni.push_back(i);\r
- }\r
- }\r
- bb = bounds_fast(c[ei]);\r
- dsq = distanceSq(p, bb);\r
- if ( mindistsq >= dsq )\r
- {\r
- all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));\r
- dsq = distanceSq( p, c[ei](all_t.front()) );\r
- if ( mindistsq > dsq )\r
- {\r
- for ( unsigned int i = 0; i < all_t.size(); ++i )\r
- {\r
- all_t[i] = c.mapToDomain(all_t[i], ei);\r
- }\r
- return all_t;\r
- }\r
- else if ( mindistsq == dsq )\r
- {\r
- all_np.push_back(all_t);\r
- ni.push_back(ei);\r
- }\r
- }\r
- std::vector<double> all_nearest;\r
- for ( unsigned int i = 0; i < all_np.size(); ++i )\r
- {\r
- for ( unsigned int j = 0; j < all_np[i].size(); ++j )\r
- {\r
- all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) );\r
- }\r
- }\r
- return all_nearest;\r
-}\r
-\r
-} // end namespace Geom\r
-\r
-\r
+/*
+ * 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
+
+
index 0b43ce67b51ee195cce8b5edb4dc266dc2bcbab4..c8588214bfc3d1f35e7e5c00c77d01374c858229 100644 (file)
-/*\r
- * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>\r
- *\r
- *\r
- * Authors:\r
- * \r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- * \r
- * Copyright 2007-2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#ifndef _NEAREST_POINT_H_\r
-#define _NEAREST_POINT_H_\r
-\r
-\r
-#include <vector>\r
-\r
-#include <2geom/d2.h>\r
-#include <2geom/piecewise.h>\r
-#include <2geom/exception.h>\r
-\r
-\r
-\r
-namespace Geom\r
-{\r
-\r
-/*\r
- * Given a line L specified by a point A and direction vector v,\r
- * return the point on L nearest to p. Note that the returned value\r
- * is with respect to the _normalized_ direction of v!\r
- */\r
-inline double nearest_point(Point const &p, Point const &A, Point const &v)\r
-{\r
- Point d(p - A);\r
- return d[0] * v[0] + d[1] * v[1];\r
-}\r
-\r
-////////////////////////////////////////////////////////////////////////////////\r
-// D2<SBasis> versions\r
-\r
-/*\r
- * Return the parameter t of a nearest point on the portion of the curve "c", \r
- * related to the interval [from, to], to the point "p".\r
- * The needed curve derivative "dc" is passed as parameter.\r
- * The function return the first nearest point to "p" that is found.\r
- */\r
-double nearest_point( Point const& p,\r
- D2<SBasis> const& c, D2<SBasis> const& dc, \r
- double from = 0, double to = 1 );\r
-\r
-inline\r
-double nearest_point( Point const& p, \r
- D2<SBasis> const& c, \r
- double from = 0, double to = 1 )\r
-{\r
- return nearest_point(p, c, Geom::derivative(c), from, to);\r
-}\r
-\r
-/*\r
- * Return the parameters t of all the nearest points on the portion of \r
- * the curve "c", related to the interval [from, to], to the point "p".\r
- * The needed curve derivative "dc" is passed as parameter.\r
- */\r
-std::vector<double> \r
-all_nearest_points( Point const& p, \r
- D2<SBasis> const& c, D2<SBasis> const& dc, \r
- double from = 0, double to = 1 );\r
-\r
-inline\r
-std::vector<double> \r
-all_nearest_points( Point const& p, \r
- D2<SBasis> const& c, \r
- double from = 0, double to = 1 )\r
-{\r
- return all_nearest_points(p, c, Geom::derivative(c), from, to);\r
-}\r
-\r
-\r
-////////////////////////////////////////////////////////////////////////////////\r
-// Piecewise< D2<SBasis> > versions\r
-\r
-double nearest_point( Point const& p, \r
- Piecewise< D2<SBasis> > const& c, \r
- double from, double to );\r
-\r
-inline\r
-double nearest_point( Point const& p, Piecewise< D2<SBasis> > const& c ) \r
-{\r
- return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]);\r
-}\r
-\r
-\r
-std::vector<double> \r
-all_nearest_points( Point const& p, \r
- Piecewise< D2<SBasis> > const& c, \r
- double from, double to );\r
-\r
-inline\r
-std::vector<double> \r
-all_nearest_points( Point const& p, Piecewise< D2<SBasis> > const& c ) \r
-{\r
- return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]);\r
-}\r
-\r
-} // end namespace Geom\r
-\r
-\r
-\r
-#endif /*_NEAREST_POINT_H_*/\r
+/*
+ * 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_*/
index 145be40e4292b884568c4cd0a2a2969b794e2af5..cc31133727079f7fb2eb2c140d9f7a73df88d169 100644 (file)
-/*\r
- * Fitting Models for Geom Types\r
- *\r
- * Authors:\r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- *\r
- * Copyright 2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#ifndef _NL_FITTING_MODEL_H_\r
-#define _NL_FITTING_MODEL_H_\r
-\r
-\r
-#include <2geom/d2.h>\r
-#include <2geom/sbasis.h>\r
-#include <2geom/bezier.h>\r
-#include <2geom/bezier-curve.h>\r
-#include <2geom/poly.h>\r
-#include <2geom/ellipse.h>\r
-#include <2geom/utils.h>\r
-\r
-\r
-namespace Geom { namespace NL {\r
-\r
-\r
-/*\r
- * completely unknown models must inherit from this template class;\r
- * example: the model a*x^2 + b*x + c = 0 to be solved wrt a, b, c;\r
- * example: the model A(t) = known_sample_value_at(t) to be solved wrt\r
- * the coefficients of the curve A(t) expressed in S-Basis form;\r
- * parameter type: the type of x and t variable in the examples above;\r
- * value type: the type of the known sample values (in the first example\r
- * is constant )\r
- * instance type: the type of the objects produced by using\r
- * the fitting raw data solution\r
- */\r
-template< typename ParameterType, typename ValueType, typename InstanceType >\r
-class LinearFittingModel\r
-{\r
- public:\r
- typedef ParameterType parameter_type;\r
- typedef ValueType value_type;\r
- typedef InstanceType instance_type;\r
-\r
- static const bool WITH_FIXED_TERMS = false;\r
-\r
- /*\r
- * a LinearFittingModel must implement the following methods:\r
- *\r
- * void feed( VectorView & vector,\r
- * parameter_type const& sample_parameter ) const;\r
- *\r
- * size_t size() const;\r
- *\r
- * void instance(instance_type &, raw_type const& raw_data) const;\r
- *\r
- */\r
-};\r
-\r
-\r
-/*\r
- * partially known models must inherit from this template class\r
- * example: the model a*x^2 + 2*x + c = 0 to be solved wrt a and c\r
- */\r
-template< typename ParameterType, typename ValueType, typename InstanceType >\r
-class LinearFittingModelWithFixedTerms\r
-{\r
- public:\r
- typedef ParameterType parameter_type;\r
- typedef ValueType value_type;\r
- typedef InstanceType instance_type;\r
-\r
- static const bool WITH_FIXED_TERMS = true;\r
-\r
- /*\r
- * a LinearFittingModelWithFixedTerms must implement the following methods:\r
- *\r
- * void feed( VectorView & vector,\r
- * value_type & fixed_term,\r
- * parameter_type const& sample_parameter ) const;\r
- *\r
- * size_t size() const;\r
- *\r
- * void instance(instance_type &, raw_type const& raw_data) const;\r
- *\r
- */\r
-\r
-\r
-};\r
-\r
-\r
-// incomplete model, it can be inherited to make up different kinds of\r
-// instance type; the raw data is a vector of coefficients of a polynomial\r
-// rapresented in standard power basis\r
-template< typename InstanceType >\r
-class LFMPowerBasis\r
- : public LinearFittingModel<double, double, InstanceType>\r
-{\r
- public:\r
- LFMPowerBasis(size_t degree)\r
- : m_size(degree + 1)\r
- {\r
- }\r
-\r
- void feed( VectorView & coeff, double sample_parameter ) const\r
- {\r
- coeff[0] = 1;\r
- double x_i = 1;\r
- for (size_t i = 1; i < coeff.size(); ++i)\r
- {\r
- x_i *= sample_parameter;\r
- coeff[i] = x_i;\r
- }\r
- }\r
-\r
- size_t size() const\r
- {\r
- return m_size;\r
- }\r
-\r
- private:\r
- size_t m_size;\r
-};\r
-\r
-\r
-// this model generates Geom::Poly objects\r
-class LFMPoly\r
- : public LFMPowerBasis<Poly>\r
-{\r
- public:\r
- LFMPoly(size_t degree)\r
- : LFMPowerBasis<Poly>(degree)\r
- {\r
- }\r
-\r
- void instance(Poly & poly, ConstVectorView const& raw_data) const\r
- {\r
- poly.clear();\r
- poly.resize(size());\r
- for (size_t i = 0; i < raw_data.size(); ++i)\r
- {\r
- poly[i] = raw_data[i];\r
- }\r
- }\r
-};\r
-\r
-\r
-// incomplete model, it can be inherited to make up different kinds of\r
-// instance type; the raw data is a vector of coefficients of a polynomial\r
-// rapresented in standard power basis with leading term coefficient equal to 1\r
-template< typename InstanceType >\r
-class LFMNormalizedPowerBasis\r
- : public LinearFittingModelWithFixedTerms<double, double, InstanceType>\r
-{\r
- public:\r
- LFMNormalizedPowerBasis(size_t _degree)\r
- : m_model( _degree - 1)\r
- {\r
- assert(_degree > 0);\r
- }\r
-\r
-\r
- void feed( VectorView & coeff,\r
- double & known_term,\r
- double sample_parameter ) const\r
- {\r
- m_model.feed(coeff, sample_parameter);\r
- known_term = coeff[m_model.size()-1] * sample_parameter;\r
- }\r
-\r
- size_t size() const\r
- {\r
- return m_model.size();\r
- }\r
-\r
- private:\r
- LFMPowerBasis<InstanceType> m_model;\r
-};\r
-\r
-\r
-// incomplete model, it can be inherited to make up different kinds of\r
-// instance type; the raw data is a vector of coefficients of the equation\r
-// of an ellipse curve\r
-template< typename InstanceType >\r
-class LFMEllipseEquation\r
- : public LinearFittingModelWithFixedTerms<Point, double, InstanceType>\r
-{\r
- public:\r
- void feed( VectorView & coeff, double & fixed_term, Point const& p ) const\r
- {\r
- coeff[0] = p[X] * p[Y];\r
- coeff[1] = p[Y] * p[Y];\r
- coeff[2] = p[X];\r
- coeff[3] = p[Y];\r
- coeff[4] = 1;\r
- fixed_term = p[X] * p[X];\r
- }\r
-\r
- size_t size() const\r
- {\r
- return 5;\r
- }\r
-};\r
-\r
-\r
-// this model generates Ellipse curves\r
-class LFMEllipse\r
- : public LFMEllipseEquation<Ellipse>\r
-{\r
- public:\r
- void instance(Ellipse & e, ConstVectorView const& coeff) const\r
- {\r
- e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);\r
- }\r
-};\r
-\r
-\r
-// this model generates SBasis objects\r
-class LFMSBasis\r
- : public LinearFittingModel<double, double, SBasis>\r
-{\r
- public:\r
- LFMSBasis( size_t _order )\r
- : m_size( 2*(_order+1) ),\r
- m_order(_order)\r
- {\r
- }\r
-\r
- void feed( VectorView & coeff, double t ) const\r
- {\r
- double u0 = 1-t;\r
- double u1 = t;\r
- double s = u0 * u1;\r
- coeff[0] = u0;\r
- coeff[1] = u1;\r
- for (size_t i = 2; i < size(); i+=2)\r
- {\r
- u0 *= s;\r
- u1 *= s;\r
- coeff[i] = u0;\r
- coeff[i+1] = u1;\r
- }\r
- }\r
-\r
- size_t size() const\r
- {\r
- return m_size;\r
- }\r
-\r
- void instance(SBasis & sb, ConstVectorView const& raw_data) const\r
- {\r
- sb.clear();\r
- sb.resize(m_order+1);\r
- for (unsigned int i = 0, k = 0; i < raw_data.size(); i+=2, ++k)\r
- {\r
- sb[k][0] = raw_data[i];\r
- sb[k][1] = raw_data[i+1];\r
- }\r
- }\r
-\r
- private:\r
- size_t m_size;\r
- size_t m_order;\r
-};\r
-\r
-\r
-// this model generates D2<SBasis> objects\r
-class LFMD2SBasis\r
- : public LinearFittingModel< double, Point, D2<SBasis> >\r
-{\r
- public:\r
- LFMD2SBasis( size_t _order )\r
- : mosb(_order)\r
- {\r
- }\r
-\r
- void feed( VectorView & coeff, double t ) const\r
- {\r
- mosb.feed(coeff, t);\r
- }\r
-\r
- size_t size() const\r
- {\r
- return mosb.size();\r
- }\r
-\r
- void instance(D2<SBasis> & d2sb, ConstMatrixView const& raw_data) const\r
- {\r
- mosb.instance(d2sb[X], raw_data.column_const_view(X));\r
- mosb.instance(d2sb[Y], raw_data.column_const_view(Y));\r
- }\r
-\r
- private:\r
- LFMSBasis mosb;\r
-};\r
-\r
-\r
-// this model generates Bezier objects\r
-class LFMBezier\r
- : public LinearFittingModel<double, double, Bezier>\r
-{\r
- public:\r
- LFMBezier( size_t _order )\r
- : m_size(_order + 1),\r
- m_order(_order)\r
- {\r
- binomial_coefficients(m_bc, m_order);\r
- }\r
-\r
- void feed( VectorView & coeff, double t ) const\r
- {\r
- double s = 1;\r
- for (size_t i = 0; i < size(); ++i)\r
- {\r
- coeff[i] = s * m_bc[i];\r
- s *= t;\r
- }\r
- double u = 1-t;\r
- s = 1;\r
- for (size_t i = size()-1; i > 0; --i)\r
- {\r
- coeff[i] *= s;\r
- s *= u;\r
- }\r
- coeff[0] *= s;\r
- }\r
-\r
- size_t size() const\r
- {\r
- return m_size;\r
- }\r
-\r
- void instance(Bezier & b, ConstVectorView const& raw_data) const\r
- {\r
- assert(b.size() == raw_data.size());\r
- for (unsigned int i = 0; i < raw_data.size(); ++i)\r
- {\r
- b[i] = raw_data[i];\r
- }\r
- }\r
-\r
- private:\r
- size_t m_size;\r
- size_t m_order;\r
- std::vector<size_t> m_bc;\r
-};\r
-\r
-\r
-// this model generates Bezier curves\r
-template< unsigned int N >\r
-class LFMBezierCurve\r
- : public LinearFittingModel< double, Point, BezierCurve<N> >\r
-{\r
- public:\r
- LFMBezierCurve( size_t _order )\r
- : mob(_order)\r
- {\r
- }\r
-\r
- void feed( VectorView & coeff, double t ) const\r
- {\r
- mob.feed(coeff, t);\r
- }\r
-\r
- size_t size() const\r
- {\r
- return mob.size();\r
- }\r
-\r
- void instance(BezierCurve<N> & bc, ConstMatrixView const& raw_data) const\r
- {\r
- Bezier bx(size()-1);\r
- Bezier by(size()-1);\r
- mob.instance(bx, raw_data.column_const_view(X));\r
- mob.instance(by, raw_data.column_const_view(Y));\r
- bc = BezierCurve<N>(bx, by);\r
- }\r
-\r
- private:\r
- LFMBezier mob;\r
-};\r
-\r
-} // end namespace NL\r
-} // end namespace Geom\r
-\r
-\r
-#endif // _NL_FITTING_MODEL_H_\r
-\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+/*
+ * 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 :
index edacc663a2b8f17cdac4b4109225089311fc2afb..d589d86e39bcdd8a882ee8874c6a11e8f0dbcfe1 100644 (file)
-/*\r
- * Fitting Tools\r
- *\r
- * Authors:\r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- *\r
- * Copyright 2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#ifndef _NL_FITTING_TOOL_H_\r
-#define _NL_FITTING_TOOL_H_\r
-\r
-\r
-#include <2geom/numeric/vector.h>\r
-#include <2geom/numeric/matrix.h>\r
-\r
-#include <2geom/point.h>\r
-\r
-#include <vector>\r
-\r
-\r
-namespace Geom { namespace NL {\r
-\r
-namespace detail {\r
-\r
-\r
-template< typename ModelT>\r
-class lsf_base\r
-{\r
- public:\r
- typedef ModelT model_type;\r
- typedef typename model_type::parameter_type parameter_type;\r
- typedef typename model_type::value_type value_type;\r
-\r
- lsf_base( model_type const& _model, size_t forecasted_samples )\r
- : m_model(_model),\r
- m_total_samples(0),\r
- m_matrix(forecasted_samples, m_model.size()),\r
- m_psdinv_matrix(NULL)\r
- {\r
- }\r
-\r
- // compute pseudo inverse\r
- void update()\r
- {\r
- if (total_samples() == 0) return;\r
- if (m_psdinv_matrix != NULL)\r
- {\r
- delete m_psdinv_matrix;\r
- }\r
- MatrixView mv(m_matrix, 0, 0, total_samples(), m_matrix.columns());\r
- m_psdinv_matrix = new Matrix( pseudo_inverse(mv) );\r
- assert(m_psdinv_matrix != NULL);\r
- }\r
-\r
- size_t total_samples() const\r
- {\r
- return m_total_samples;\r
- }\r
-\r
- bool is_full() const\r
- {\r
- return (total_samples() == m_matrix.rows());\r
- }\r
-\r
- void clear()\r
- {\r
- m_total_samples = 0;\r
- }\r
-\r
- virtual\r
- ~lsf_base()\r
- {\r
- if (m_psdinv_matrix != NULL)\r
- {\r
- delete m_psdinv_matrix;\r
- }\r
- }\r
-\r
- protected:\r
- const model_type & m_model;\r
- size_t m_total_samples;\r
- Matrix m_matrix;\r
- Matrix* m_psdinv_matrix;\r
-\r
-}; // end class lsf_base\r
-\r
-\r
-\r
-\r
-template< typename ModelT, typename ValueType = typename ModelT::value_type>\r
-class lsf_solution\r
-{\r
-};\r
-\r
-// a fitting process on samples with value of type double\r
-// produces a solution of type Vector\r
-template< typename ModelT>\r
-class lsf_solution<ModelT, double>\r
- : public lsf_base<ModelT>\r
-{\r
-public:\r
- typedef ModelT model_type;\r
- typedef typename model_type::parameter_type parameter_type;\r
- typedef typename model_type::value_type value_type;\r
- typedef Vector solution_type;\r
- typedef lsf_base<model_type> base_type;\r
-\r
- using base_type::m_model;\r
- using base_type::m_psdinv_matrix;\r
- using base_type::total_samples;\r
-\r
-public:\r
- lsf_solution<ModelT, double>( model_type const& _model,\r
- size_t forecasted_samples )\r
- : base_type(_model, forecasted_samples),\r
- m_solution(_model.size())\r
- {\r
- }\r
-\r
- template< typename VectorT >\r
- solution_type& result(VectorT const& sample_values)\r
- {\r
- assert(sample_values.size() == total_samples());\r
- ConstVectorView sv(sample_values);\r
- m_solution = (*m_psdinv_matrix) * sv;\r
- return m_solution;\r
- }\r
-\r
- // a comparison between old sample values and the new ones is performed\r
- // in order to minimize computation\r
- // prerequisite:\r
- // old_sample_values.size() == new_sample_values.size()\r
- // no update() call can be performed between two result invocations\r
- template< typename VectorT >\r
- solution_type& result( VectorT const& old_sample_values,\r
- VectorT const& new_sample_values )\r
- {\r
- assert(old_sample_values.size() == total_samples());\r
- assert(new_sample_values.size() == total_samples());\r
- Vector diff(total_samples());\r
- for (size_t i = 0; i < diff.size(); ++i)\r
- {\r
- diff[i] = new_sample_values[i] - old_sample_values[i];\r
- }\r
- Vector column(m_model.size());\r
- Vector delta(m_model.size(), 0.0);\r
- for (size_t i = 0; i < diff.size(); ++i)\r
- {\r
- if (diff[i] != 0)\r
- {\r
- column = m_psdinv_matrix->column_view(i);\r
- column.scale(diff[i]);\r
- delta += column;\r
- }\r
- }\r
- m_solution += delta;\r
- return m_solution;\r
- }\r
-\r
- solution_type& result()\r
- {\r
- return m_solution;\r
- }\r
-\r
-private:\r
- solution_type m_solution;\r
-\r
-}; // end class lsf_solution<ModelT, double>\r
-\r
-\r
-// a fitting process on samples with value of type Point\r
-// produces a solution of type Matrix (with 2 columns)\r
-template< typename ModelT>\r
-class lsf_solution<ModelT, Point>\r
- : public lsf_base<ModelT>\r
-{\r
-public:\r
- typedef ModelT model_type;\r
- typedef typename model_type::parameter_type parameter_type;\r
- typedef typename model_type::value_type value_type;\r
- typedef Matrix solution_type;\r
- typedef lsf_base<model_type> base_type;\r
-\r
- using base_type::m_model;\r
- using base_type::m_psdinv_matrix;\r
- using base_type::total_samples;\r
-\r
-public:\r
- lsf_solution<ModelT, Point>( model_type const& _model,\r
- size_t forecasted_samples )\r
- : base_type(_model, forecasted_samples),\r
- m_solution(_model.size(), 2)\r
- {\r
- }\r
-\r
- solution_type& result(std::vector<Point> const& sample_values)\r
- {\r
- assert(sample_values.size() == total_samples());\r
- Matrix svm(total_samples(), 2);\r
- for (size_t i = 0; i < total_samples(); ++i)\r
- {\r
- svm(i, X) = sample_values[i][X];\r
- svm(i, Y) = sample_values[i][Y];\r
- }\r
- m_solution = (*m_psdinv_matrix) * svm;\r
- return m_solution;\r
- }\r
-\r
- // a comparison between old sample values and the new ones is performed\r
- // in order to minimize computation\r
- // prerequisite:\r
- // old_sample_values.size() == new_sample_values.size()\r
- // no update() call can to be performed between two result invocations\r
- solution_type& result( std::vector<Point> const& old_sample_values,\r
- std::vector<Point> const& new_sample_values )\r
- {\r
- assert(old_sample_values.size() == total_samples());\r
- assert(new_sample_values.size() == total_samples());\r
- Matrix diff(total_samples(), 2);\r
- for (size_t i = 0; i < total_samples(); ++i)\r
- {\r
- diff(i, X) = new_sample_values[i][X] - old_sample_values[i][X];\r
- diff(i, Y) = new_sample_values[i][Y] - old_sample_values[i][Y];\r
- }\r
- Vector column(m_model.size());\r
- Matrix delta(m_model.size(), 2, 0.0);\r
- VectorView deltax = delta.column_view(X);\r
- VectorView deltay = delta.column_view(Y);\r
- for (size_t i = 0; i < total_samples(); ++i)\r
- {\r
- if (diff(i, X) != 0)\r
- {\r
- column = m_psdinv_matrix->column_view(i);\r
- column.scale(diff(i, X));\r
- deltax += column;\r
- }\r
- if (diff(i, Y) != 0)\r
- {\r
- column = m_psdinv_matrix->column_view(i);\r
- column.scale(diff(i, Y));\r
- deltay += column;\r
- }\r
- }\r
- m_solution += delta;\r
- return m_solution;\r
- }\r
-\r
- solution_type& result()\r
- {\r
- return m_solution;\r
- }\r
-\r
-private:\r
- solution_type m_solution;\r
-\r
-}; // end class lsf_solution<ModelT, Point>\r
-\r
-\r
-\r
-\r
-template< typename ModelT,\r
- bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >\r
-class lsf_with_fixed_terms\r
-{\r
-};\r
-\r
-\r
-// fitting tool for completely unknown models\r
-template< typename ModelT>\r
-class lsf_with_fixed_terms<ModelT, false>\r
- : public lsf_solution<ModelT>\r
-{\r
- public:\r
- typedef ModelT model_type;\r
- typedef typename model_type::parameter_type parameter_type;\r
- typedef typename model_type::value_type value_type;\r
- typedef lsf_solution<model_type> base_type;\r
- typedef typename base_type::solution_type solution_type;\r
-\r
- using base_type::total_samples;\r
- using base_type::is_full;\r
- using base_type::m_matrix;\r
- using base_type::m_total_samples;\r
- using base_type::m_model;\r
-\r
- public:\r
- lsf_with_fixed_terms<ModelT, false>( model_type const& _model,\r
- size_t forecasted_samples )\r
- : base_type(_model, forecasted_samples)\r
- {\r
- }\r
-\r
- void append(parameter_type const& sample_parameter)\r
- {\r
- assert(!is_full());\r
- VectorView row = m_matrix.row_view(total_samples());\r
- m_model.feed(row, sample_parameter);\r
- ++m_total_samples;\r
- }\r
-\r
- void append_copy(size_t sample_index)\r
- {\r
- assert(!is_full());\r
- assert(sample_index < total_samples());\r
- VectorView dest_row = m_matrix.row_view(total_samples());\r
- VectorView source_row = m_matrix.row_view(sample_index);\r
- dest_row = source_row;\r
- ++m_total_samples;\r
- }\r
-\r
-}; // end class lsf_with_fixed_terms<ModelT, false>\r
-\r
-\r
-// fitting tool for partially known models\r
-template< typename ModelT>\r
-class lsf_with_fixed_terms<ModelT, true>\r
- : public lsf_solution<ModelT>\r
-{\r
- public:\r
- typedef ModelT model_type;\r
- typedef typename model_type::parameter_type parameter_type;\r
- typedef typename model_type::value_type value_type;\r
- typedef lsf_solution<model_type> base_type;\r
- typedef typename base_type::solution_type solution_type;\r
-\r
- using base_type::total_samples;\r
- using base_type::is_full;\r
- using base_type::m_matrix;\r
- using base_type::m_total_samples;\r
- using base_type::m_model;\r
-\r
- public:\r
- lsf_with_fixed_terms<ModelT, true>( model_type const& _model,\r
- size_t forecasted_samples )\r
- : base_type(_model, forecasted_samples),\r
- m_vector(forecasted_samples),\r
- m_vector_view(NULL)\r
- {\r
- }\r
- void append(parameter_type const& sample_parameter)\r
- {\r
- assert(!is_full());\r
- VectorView row = m_matrix.row_view(total_samples());\r
- m_model.feed(row, m_vector[total_samples()], sample_parameter);\r
- ++m_total_samples;\r
- }\r
-\r
- void append_copy(size_t sample_index)\r
- {\r
- assert(!is_full());\r
- assert(sample_index < total_samples());\r
- VectorView dest_row = m_matrix.row_view(total_samples());\r
- VectorView source_row = m_matrix.row_view(sample_index);\r
- dest_row = source_row;\r
- m_vector[total_samples()] = m_vector[sample_index];\r
- ++m_total_samples;\r
- }\r
-\r
- void update()\r
- {\r
- base_type::update();\r
- if (total_samples() == 0) return;\r
- if (m_vector_view != NULL)\r
- {\r
- delete m_vector_view;\r
- }\r
- m_vector_view = new VectorView(m_vector, base_type::total_samples());\r
- assert(m_vector_view != NULL);\r
- }\r
-\r
- virtual\r
- ~lsf_with_fixed_terms<model_type, true>()\r
- {\r
- if (m_vector_view != NULL)\r
- {\r
- delete m_vector_view;\r
- }\r
- }\r
-\r
- protected:\r
- Vector m_vector;\r
- VectorView* m_vector_view;\r
-\r
-}; // end class lsf_with_fixed_terms<ModelT, true>\r
-\r
-\r
-} // end namespace detail\r
-\r
-\r
-\r
-\r
-template< typename ModelT,\r
- typename ValueType = typename ModelT::value_type,\r
- bool WITH_FIXED_TERMS = ModelT::WITH_FIXED_TERMS >\r
-class least_squeares_fitter\r
-{\r
-};\r
-\r
-\r
-template< typename ModelT, typename ValueType >\r
-class least_squeares_fitter<ModelT, ValueType, false>\r
- : public detail::lsf_with_fixed_terms<ModelT>\r
-{\r
- public:\r
- typedef ModelT model_type;\r
- typedef detail::lsf_with_fixed_terms<model_type> base_type;\r
- typedef typename base_type::parameter_type parameter_type;\r
- typedef typename base_type::value_type value_type;\r
- typedef typename base_type::solution_type solution_type;\r
-\r
- public:\r
- least_squeares_fitter<ModelT, ValueType, false>( model_type const& _model,\r
- size_t forecasted_samples )\r
- : base_type(_model, forecasted_samples)\r
- {\r
- }\r
-}; // end class least_squeares_fitter<ModelT, ValueType, true>\r
-\r
-\r
-template< typename ModelT>\r
-class least_squeares_fitter<ModelT, double, true>\r
- : public detail::lsf_with_fixed_terms<ModelT>\r
-{\r
- public:\r
- typedef ModelT model_type;\r
- typedef detail::lsf_with_fixed_terms<model_type> base_type;\r
- typedef typename base_type::parameter_type parameter_type;\r
- typedef typename base_type::value_type value_type;\r
- typedef typename base_type::solution_type solution_type;\r
-\r
- using base_type::m_vector_view;\r
- using base_type::result;\r
-\r
- public:\r
- least_squeares_fitter<ModelT, double, true>( model_type const& _model,\r
- size_t forecasted_samples )\r
- : base_type(_model, forecasted_samples)\r
- {\r
- }\r
-\r
- template< typename VectorT >\r
- solution_type& result(VectorT const& sample_values)\r
- {\r
- assert(sample_values.size() == m_vector_view->size());\r
- Vector sv(sample_values.size());\r
- for (size_t i = 0; i < sv.size(); ++i)\r
- sv[i] = sample_values[i] - (*m_vector_view)[i];\r
- return base_type::result(sv);\r
- }\r
-\r
-}; // end class least_squeares_fitter<ModelT, double, true>\r
-\r
-\r
-template< typename ModelT>\r
-class least_squeares_fitter<ModelT, Point, true>\r
- : public detail::lsf_with_fixed_terms<ModelT>\r
-{\r
- public:\r
- typedef ModelT model_type;\r
- typedef detail::lsf_with_fixed_terms<model_type> base_type;\r
- typedef typename base_type::parameter_type parameter_type;\r
- typedef typename base_type::value_type value_type;\r
- typedef typename base_type::solution_type solution_type;\r
-\r
- using base_type::m_vector_view;\r
- using base_type::result;\r
-\r
- public:\r
- least_squeares_fitter<ModelT, Point, true>( model_type const& _model,\r
- size_t forecasted_samples )\r
- : base_type(_model, forecasted_samples)\r
- {\r
- }\r
-\r
- solution_type& result(std::vector<Point> const& sample_values)\r
- {\r
- assert(sample_values.size() == m_vector_view->size());\r
- NL::Matrix sv(sample_values.size(), 2);\r
- for (size_t i = 0; i < sample_values.size(); ++i)\r
- {\r
- sv(i, X) = sample_values[i][X] - (*m_vector_view)[i];\r
- sv(i, Y) = sample_values[i][Y] - (*m_vector_view)[i];\r
- }\r
- return base_type::result(sv);\r
- }\r
-\r
-}; // end class least_squeares_fitter<ModelT, Point, true>\r
-\r
-\r
-} // end namespace NL\r
-} // end namespace Geom\r
-\r
-\r
-\r
-#endif // _NL_FITTING_TOOL_H_\r
-\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+/*
+ * 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 :
index 5b516c9e6ba6bf5ff0554ee3297f466e4cf80c34..dc2a1d7e01fe2b73e836a9026041cbebd754e62c 100644 (file)
-/*\r
- * LinearSystem class wraps some gsl routines for solving linear systems\r
- *\r
- * Authors:\r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- * \r
- * Copyright 2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#ifndef _NL_LINEAR_SYSTEM_H_\r
-#define _NL_LINEAR_SYSTEM_H_\r
-\r
-\r
-#include <cassert>\r
-\r
-#include <gsl/gsl_linalg.h>\r
-\r
-#include <2geom/numeric/matrix.h>\r
-#include <2geom/numeric/vector.h>\r
-\r
-\r
-namespace Geom { namespace NL {\r
-\r
-\r
-class LinearSystem\r
-{\r
-public:\r
- LinearSystem(MatrixView & _matrix, VectorView & _vector)\r
- : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())\r
- {\r
- }\r
- \r
- LinearSystem(Matrix & _matrix, Vector & _vector)\r
- : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())\r
- {\r
- }\r
- \r
- const Vector & LU_solve()\r
- {\r
- assert( matrix().rows() == matrix().columns() \r
- && matrix().rows() == vector().size() );\r
- int s;\r
- gsl_permutation * p = gsl_permutation_alloc(matrix().rows());\r
- gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);\r
- gsl_linalg_LU_solve( matrix().get_gsl_matrix(), \r
- p, \r
- vector().get_gsl_vector(), \r
- m_solution.get_gsl_vector()\r
- );\r
- gsl_permutation_free(p);\r
- return solution();\r
- }\r
- \r
- const Vector & SV_solve()\r
- {\r
- assert( matrix().rows() >= matrix().columns()\r
- && matrix().rows() == vector().size() );\r
- \r
- gsl_matrix* U = matrix().get_gsl_matrix();\r
- gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());\r
- gsl_vector* S = gsl_vector_alloc(matrix().columns());\r
- gsl_vector* work = gsl_vector_alloc(matrix().columns());\r
- \r
- gsl_linalg_SV_decomp( U, V, S, work );\r
- \r
- gsl_vector* b = vector().get_gsl_vector();\r
- gsl_vector* x = m_solution.get_gsl_vector();\r
- \r
- gsl_linalg_SV_solve( U, V, S, b, x);\r
- \r
- gsl_matrix_free(V);\r
- gsl_vector_free(S);\r
- gsl_vector_free(work);\r
- \r
- return solution(); \r
- }\r
- \r
- MatrixView & matrix()\r
- {\r
- return m_matrix;\r
- }\r
- \r
- VectorView & vector()\r
- {\r
- return m_vector;\r
- }\r
- \r
- const Vector & solution() const\r
- {\r
- return m_solution;\r
- }\r
- \r
-private:\r
- MatrixView m_matrix;\r
- VectorView m_vector;\r
- Vector m_solution;\r
-};\r
-\r
-\r
-} } // end namespaces\r
-\r
-\r
-#endif /*_NL_LINEAR_SYSTEM_H_*/\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+/*
+ * 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 :
index 64557a6f12ed64ffc28260871136adc3641501b9..156b6e9a25b24a2b0fae0f98918c332bd894a942 100644 (file)
-/*\r
- * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines;\r
- * "views" mimic the semantic of C++ references: any operation performed\r
- * on a "view" is actually performed on the "viewed object"\r
- *\r
- * Authors:\r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- *\r
- * Copyright 2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-\r
-\r
-#ifndef _NL_MATRIX_H_\r
-#define _NL_MATRIX_H_\r
-\r
-#include <2geom/numeric/vector.h>\r
-\r
-#include <cassert>\r
-#include <utility> // for std::pair\r
-#include <algorithm> // for std::swap\r
-#include <sstream>\r
-#include <string>\r
-\r
-#include <gsl/gsl_matrix.h>\r
-#include <gsl/gsl_linalg.h>\r
-\r
-\r
-namespace Geom { namespace NL {\r
-\r
-namespace detail\r
-{\r
-\r
-class BaseMatrixImpl\r
-{\r
- public:\r
- virtual ~BaseMatrixImpl()\r
- {\r
- }\r
-\r
- ConstVectorView row_const_view(size_t i) const\r
- {\r
- return ConstVectorView(gsl_matrix_const_row(m_matrix, i));\r
- }\r
-\r
- ConstVectorView column_const_view(size_t i) const\r
- {\r
- return ConstVectorView(gsl_matrix_const_column(m_matrix, i));\r
- }\r
-\r
- const double & operator() (size_t i, size_t j) const\r
- {\r
- return *gsl_matrix_const_ptr(m_matrix, i, j);\r
- }\r
-\r
- const gsl_matrix* get_gsl_matrix() const\r
- {\r
- return m_matrix;\r
- }\r
-\r
- bool is_zero() const\r
- {\r
- return gsl_matrix_isnull(m_matrix);\r
- }\r
-\r
- bool is_positive() const\r
- {\r
- return gsl_matrix_ispos(m_matrix);\r
- }\r
-\r
- bool is_negative() const\r
- {\r
- return gsl_matrix_isneg(m_matrix);\r
- }\r
-\r
- bool is_non_negative() const\r
- {\r
- for ( unsigned int i = 0; i < rows(); ++i )\r
- {\r
- for ( unsigned int j = 0; j < columns(); ++j )\r
- {\r
- if ( (*this)(i,j) < 0 ) return false;\r
- }\r
- }\r
- return true;\r
- }\r
-\r
- double max() const\r
- {\r
- return gsl_matrix_max(m_matrix);\r
- }\r
-\r
- double min() const\r
- {\r
- return gsl_matrix_min(m_matrix);\r
- }\r
-\r
- std::pair<size_t, size_t>\r
- max_index() const\r
- {\r
- std::pair<size_t, size_t> indices;\r
- gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));\r
- return indices;\r
- }\r
-\r
- std::pair<size_t, size_t>\r
- min_index() const\r
- {\r
- std::pair<size_t, size_t> indices;\r
- gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));\r
- return indices;\r
- }\r
-\r
- size_t rows() const\r
- {\r
- return m_rows;\r
- }\r
-\r
- size_t columns() const\r
- {\r
- return m_columns;\r
- }\r
-\r
- std::string str() const;\r
-\r
- protected:\r
- size_t m_rows, m_columns;\r
- gsl_matrix* m_matrix;\r
-\r
-}; // end class BaseMatrixImpl\r
-\r
-\r
-inline\r
-bool operator== (BaseMatrixImpl const& m1, BaseMatrixImpl const& m2)\r
-{\r
- if (m1.rows() != m2.rows() || m1.columns() != m2.columns()) return false;\r
-\r
- for (size_t i = 0; i < m1.rows(); ++i)\r
- for (size_t j = 0; j < m1.columns(); ++j)\r
- if (m1(i,j) != m2(i,j)) return false;\r
-\r
- return true;\r
-}\r
-\r
-template< class charT >\r
-inline\r
-std::basic_ostream<charT> &\r
-operator<< (std::basic_ostream<charT> & os, const BaseMatrixImpl & _matrix)\r
-{\r
- if (_matrix.rows() == 0 || _matrix.columns() == 0) return os;\r
-\r
- os << "[[" << _matrix(0,0);\r
- for (size_t j = 1; j < _matrix.columns(); ++j)\r
- {\r
- os << ", " << _matrix(0,j);\r
- }\r
- os << "]";\r
-\r
- for (size_t i = 1; i < _matrix.rows(); ++i)\r
- {\r
- os << ", [" << _matrix(i,0);\r
- for (size_t j = 1; j < _matrix.columns(); ++j)\r
- {\r
- os << ", " << _matrix(i,j);\r
- }\r
- os << "]";\r
- }\r
- os << "]";\r
- return os;\r
-}\r
-\r
-inline\r
-std::string BaseMatrixImpl::str() const\r
-{\r
- std::ostringstream oss;\r
- oss << (*this);\r
- return oss.str();\r
-}\r
-\r
-\r
-class MatrixImpl : public BaseMatrixImpl\r
-{\r
- public:\r
-\r
- typedef BaseMatrixImpl base_type;\r
-\r
- void set_all( double x )\r
- {\r
- gsl_matrix_set_all(m_matrix, x);\r
- }\r
-\r
- void set_identity()\r
- {\r
- gsl_matrix_set_identity(m_matrix);\r
- }\r
-\r
- using base_type::operator();\r
-\r
- double & operator() (size_t i, size_t j)\r
- {\r
- return *gsl_matrix_ptr(m_matrix, i, j);\r
- }\r
-\r
- using base_type::get_gsl_matrix;\r
-\r
- gsl_matrix* get_gsl_matrix()\r
- {\r
- return m_matrix;\r
- }\r
-\r
- VectorView row_view(size_t i)\r
- {\r
- return VectorView(gsl_matrix_row(m_matrix, i));\r
- }\r
-\r
- VectorView column_view(size_t i)\r
- {\r
- return VectorView(gsl_matrix_column(m_matrix, i));\r
- }\r
-\r
- void swap_rows(size_t i, size_t j)\r
- {\r
- gsl_matrix_swap_rows(m_matrix, i, j);\r
- }\r
-\r
- void swap_columns(size_t i, size_t j)\r
- {\r
- gsl_matrix_swap_columns(m_matrix, i, j);\r
- }\r
-\r
- MatrixImpl & transpose()\r
- {\r
- assert(columns() == rows());\r
- gsl_matrix_transpose(m_matrix);\r
- return (*this);\r
- }\r
-\r
- MatrixImpl & scale(double x)\r
- {\r
- gsl_matrix_scale(m_matrix, x);\r
- return (*this);\r
- }\r
-\r
- MatrixImpl & translate(double x)\r
- {\r
- gsl_matrix_add_constant(m_matrix, x);\r
- return (*this);\r
- }\r
-\r
- MatrixImpl & operator+=(base_type const& _matrix)\r
- {\r
- gsl_matrix_add(m_matrix, _matrix.get_gsl_matrix());\r
- return (*this);\r
- }\r
-\r
- MatrixImpl & operator-=(base_type const& _matrix)\r
- {\r
- gsl_matrix_sub(m_matrix, _matrix.get_gsl_matrix());\r
- return (*this);\r
- }\r
-\r
-}; // end class MatrixImpl\r
-\r
-} // end namespace detail\r
-\r
-\r
-using detail::operator==;\r
-using detail::operator<<;\r
-\r
-\r
-\r
-\r
-class Matrix: public detail::MatrixImpl\r
-{\r
- public:\r
- typedef detail::MatrixImpl base_type;\r
-\r
- public:\r
- // the matrix is not inizialized\r
- Matrix(size_t n1, size_t n2)\r
- {\r
- m_rows = n1;\r
- m_columns = n2;\r
- m_matrix = gsl_matrix_alloc(n1, n2);\r
- }\r
-\r
- Matrix(size_t n1, size_t n2, double x)\r
- {\r
- m_rows = n1;\r
- m_columns = n2;\r
- m_matrix = gsl_matrix_alloc(n1, n2);\r
- gsl_matrix_set_all(m_matrix, x);\r
- }\r
-\r
- Matrix(Matrix const& _matrix)\r
- : base_type()\r
- {\r
- m_rows = _matrix.rows();\r
- m_columns = _matrix.columns();\r
- m_matrix = gsl_matrix_alloc(rows(), columns());\r
- gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
- }\r
-\r
- explicit\r
- Matrix(base_type::base_type const& _matrix)\r
- {\r
- m_rows = _matrix.rows();\r
- m_columns = _matrix.columns();\r
- m_matrix = gsl_matrix_alloc(rows(), columns());\r
- gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
- }\r
-\r
- Matrix & operator=(Matrix const& _matrix)\r
- {\r
- assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
- gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
- return *this;\r
- }\r
-\r
- Matrix & operator=(base_type::base_type const& _matrix)\r
- {\r
- assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
- gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
- return *this;\r
- }\r
-\r
- virtual ~Matrix()\r
- {\r
- gsl_matrix_free(m_matrix);\r
- }\r
-\r
- Matrix & transpose()\r
- {\r
- return static_cast<Matrix &>( base_type::transpose() );\r
- }\r
-\r
- Matrix & scale(double x)\r
- {\r
- return static_cast<Matrix &>( base_type::scale(x) );\r
- }\r
-\r
- Matrix & translate(double x)\r
- {\r
- return static_cast<Matrix &>( base_type::translate(x) );\r
- }\r
-\r
- Matrix & operator+=(base_type::base_type const& _matrix)\r
- {\r
- return static_cast<Matrix &>( base_type::operator+=(_matrix) );\r
- }\r
-\r
- Matrix & operator-=(base_type::base_type const& _matrix)\r
- {\r
- return static_cast<Matrix &>( base_type::operator-=(_matrix) );\r
- }\r
-\r
- friend\r
- void swap(Matrix & m1, Matrix & m2);\r
- friend\r
- void swap_any(Matrix & m1, Matrix & m2);\r
-\r
-}; // end class Matrix\r
-\r
-\r
-// warning! this operation invalidates any view of the passed matrix objects\r
-inline\r
-void swap(Matrix & m1, Matrix & m2)\r
-{\r
- assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );\r
- std::swap(m1.m_matrix, m2.m_matrix);\r
-}\r
-\r
-inline\r
-void swap_any(Matrix & m1, Matrix & m2)\r
-{\r
- std::swap(m1.m_matrix, m2.m_matrix);\r
- std::swap(m1.m_rows, m2.m_rows);\r
- std::swap(m1.m_columns, m2.m_columns);\r
-}\r
-\r
-\r
-\r
-class ConstMatrixView : public detail::BaseMatrixImpl\r
-{\r
- public:\r
- typedef detail::BaseMatrixImpl base_type;\r
-\r
- public:\r
- ConstMatrixView(const base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)\r
- : m_matrix_view( gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2) )\r
- {\r
- m_rows = n1;\r
- m_columns = n2;\r
- m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
- }\r
-\r
- ConstMatrixView(const ConstMatrixView & _matrix)\r
- : base_type(),\r
- m_matrix_view(_matrix.m_matrix_view)\r
- {\r
- m_rows = _matrix.rows();\r
- m_columns = _matrix.columns();\r
- m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
- }\r
-\r
- ConstMatrixView(const base_type & _matrix)\r
- : m_matrix_view(gsl_matrix_const_submatrix(_matrix.get_gsl_matrix(), 0, 0, _matrix.rows(), _matrix.columns()))\r
- {\r
- m_rows = _matrix.rows();\r
- m_columns = _matrix.columns();\r
- m_matrix = const_cast<gsl_matrix*>( &(m_matrix_view.matrix) );\r
- }\r
-\r
- private:\r
- gsl_matrix_const_view m_matrix_view;\r
-\r
-}; // end class ConstMatrixView\r
-\r
-\r
-\r
-\r
-class MatrixView : public detail::MatrixImpl\r
-{\r
- public:\r
- typedef detail::MatrixImpl base_type;\r
-\r
- public:\r
- MatrixView(base_type & _matrix, size_t k1, size_t k2, size_t n1, size_t n2)\r
- {\r
- m_rows = n1;\r
- m_columns = n2;\r
- m_matrix_view\r
- = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), k1, k2, n1, n2);\r
- m_matrix = &(m_matrix_view.matrix);\r
- }\r
-\r
- MatrixView(const MatrixView & _matrix)\r
- : base_type()\r
- {\r
- m_rows = _matrix.rows();\r
- m_columns = _matrix.columns();\r
- m_matrix_view = _matrix.m_matrix_view;\r
- m_matrix = &(m_matrix_view.matrix);\r
- }\r
-\r
- MatrixView(Matrix & _matrix)\r
- {\r
- m_rows = _matrix.rows();\r
- m_columns = _matrix.columns();\r
- m_matrix_view\r
- = gsl_matrix_submatrix(_matrix.get_gsl_matrix(), 0, 0, rows(), columns());\r
- m_matrix = &(m_matrix_view.matrix);\r
- }\r
-\r
- MatrixView & operator=(MatrixView const& _matrix)\r
- {\r
- assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
- gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);\r
- return *this;\r
- }\r
-\r
- MatrixView & operator=(base_type::base_type const& _matrix)\r
- {\r
- assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
- gsl_matrix_memcpy(m_matrix, _matrix.get_gsl_matrix());\r
- return *this;\r
- }\r
-\r
- MatrixView & transpose()\r
- {\r
- return static_cast<MatrixView &>( base_type::transpose() );\r
- }\r
-\r
- MatrixView & scale(double x)\r
- {\r
- return static_cast<MatrixView &>( base_type::scale(x) );\r
- }\r
-\r
- MatrixView & translate(double x)\r
- {\r
- return static_cast<MatrixView &>( base_type::translate(x) );\r
- }\r
-\r
- MatrixView & operator+=(base_type::base_type const& _matrix)\r
- {\r
- return static_cast<MatrixView &>( base_type::operator+=(_matrix) );\r
- }\r
-\r
- MatrixView & operator-=(base_type::base_type const& _matrix)\r
- {\r
- return static_cast<MatrixView &>( base_type::operator-=(_matrix) );\r
- }\r
-\r
- friend\r
- void swap_view(MatrixView & m1, MatrixView & m2);\r
-\r
- private:\r
- gsl_matrix_view m_matrix_view;\r
-\r
-}; // end class MatrixView\r
-\r
-\r
-inline\r
-void swap_view(MatrixView & m1, MatrixView & m2)\r
-{\r
- assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );\r
- std::swap(m1.m_matrix_view, m2.m_matrix_view);\r
-}\r
-\r
-Vector operator*( detail::BaseMatrixImpl const& A,\r
- detail::BaseVectorImpl const& v );\r
-\r
-Matrix operator*( detail::BaseMatrixImpl const& A,\r
- detail::BaseMatrixImpl const& B );\r
-\r
-Matrix pseudo_inverse(detail::BaseMatrixImpl const& A);\r
-\r
-} } // end namespaces\r
-\r
-#endif /*_NL_MATRIX_H_*/\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+/*
+ * 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 :
index 43a39a1ac5684ba5eba76c751e8d7336e456ef07..3e53405f42d0dfef76180944bf8985aefeaaf7a5 100644 (file)
-/*\r
- * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines;\r
- * "views" mimic the semantic of C++ references: any operation performed\r
- * on a "view" is actually performed on the "viewed object"\r
- *\r
- * Authors:\r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- *\r
- * Copyright 2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-\r
-\r
-#ifndef _NL_VECTOR_H_\r
-#define _NL_VECTOR_H_\r
-\r
-#include <cassert>\r
-#include <algorithm> // for std::swap\r
-#include <vector>\r
-#include <sstream>\r
-#include <string>\r
-\r
-\r
-#include <gsl/gsl_vector.h>\r
-#include <gsl/gsl_blas.h>\r
-\r
-\r
-namespace Geom { namespace NL {\r
-\r
-namespace detail\r
-{\r
-\r
-class BaseVectorImpl\r
-{\r
- public:\r
- double const& operator[](size_t i) const\r
- {\r
- return *gsl_vector_const_ptr(m_vector, i);\r
- }\r
-\r
- const gsl_vector* get_gsl_vector() const\r
- {\r
- return m_vector;\r
- }\r
- bool is_zero() const\r
- {\r
- return gsl_vector_isnull(m_vector);\r
- }\r
-\r
- bool is_positive() const\r
- {\r
- return gsl_vector_ispos(m_vector);\r
- }\r
-\r
- bool is_negative() const\r
- {\r
- return gsl_vector_isneg(m_vector);\r
- }\r
-\r
- bool is_non_negative() const\r
- {\r
- for ( size_t i = 0; i < size(); ++i )\r
- {\r
- if ( (*this)[i] < 0 ) return false;\r
- }\r
- return true;\r
- }\r
-\r
- double max() const\r
- {\r
- return gsl_vector_max(m_vector);\r
- }\r
-\r
- double min() const\r
- {\r
- return gsl_vector_min(m_vector);\r
- }\r
-\r
- size_t max_index() const\r
- {\r
- return gsl_vector_max_index(m_vector);\r
- }\r
-\r
- size_t min_index() const\r
- {\r
- return gsl_vector_min_index(m_vector);\r
- }\r
-\r
- size_t size() const\r
- {\r
- return m_size;\r
- }\r
-\r
- std::string str() const;\r
-\r
- virtual ~BaseVectorImpl()\r
- {\r
- }\r
-\r
- protected:\r
- size_t m_size;\r
- gsl_vector* m_vector;\r
-\r
-}; // end class BaseVectorImpl\r
-\r
-\r
-inline\r
-bool operator== (BaseVectorImpl const& v1, BaseVectorImpl const& v2)\r
-{\r
- if (v1.size() != v2.size()) return false;\r
-\r
- for (size_t i = 0; i < v1.size(); ++i)\r
- {\r
- if (v1[i] != v2[i]) return false;\r
- }\r
- return true;\r
-}\r
-\r
-template< class charT >\r
-inline\r
-std::basic_ostream<charT> &\r
-operator<< (std::basic_ostream<charT> & os, const BaseVectorImpl & _vector)\r
-{\r
- if (_vector.size() == 0 ) return os;\r
- os << "[" << _vector[0];\r
- for (unsigned int i = 1; i < _vector.size(); ++i)\r
- {\r
- os << ", " << _vector[i];\r
- }\r
- os << "]";\r
- return os;\r
-}\r
-\r
-inline\r
-std::string BaseVectorImpl::str() const\r
-{\r
- std::ostringstream oss;\r
- oss << (*this);\r
- return oss.str();\r
-}\r
-\r
-inline\r
-double dot(BaseVectorImpl const& v1, BaseVectorImpl const& v2)\r
-{\r
- double result;\r
- gsl_blas_ddot(v1.get_gsl_vector(), v2.get_gsl_vector(), &result);\r
- return result;\r
-}\r
-\r
-\r
-class VectorImpl : public BaseVectorImpl\r
-{\r
- public:\r
- typedef BaseVectorImpl base_type;\r
-\r
- public:\r
- void set_all(double x)\r
- {\r
- gsl_vector_set_all(m_vector, x);\r
- }\r
-\r
- void set_basis(size_t i)\r
- {\r
- gsl_vector_set_basis(m_vector, i);\r
- }\r
-\r
- using base_type::operator[];\r
-\r
- double & operator[](size_t i)\r
- {\r
- return *gsl_vector_ptr(m_vector, i);\r
- }\r
-\r
- using base_type::get_gsl_vector;\r
-\r
- gsl_vector* get_gsl_vector()\r
- {\r
- return m_vector;\r
- }\r
-\r
- void swap_elements(size_t i, size_t j)\r
- {\r
- gsl_vector_swap_elements(m_vector, i, j);\r
- }\r
-\r
- void reverse()\r
- {\r
- gsl_vector_reverse(m_vector);\r
- }\r
-\r
- VectorImpl & scale(double x)\r
- {\r
- gsl_vector_scale(m_vector, x);\r
- return (*this);\r
- }\r
-\r
- VectorImpl & translate(double x)\r
- {\r
- gsl_vector_add_constant(m_vector, x);\r
- return (*this);\r
- }\r
-\r
- VectorImpl & operator+=(base_type const& _vector)\r
- {\r
- gsl_vector_add(m_vector, _vector.get_gsl_vector());\r
- return (*this);\r
- }\r
-\r
- VectorImpl & operator-=(base_type const& _vector)\r
- {\r
- gsl_vector_sub(m_vector, _vector.get_gsl_vector());\r
- return (*this);\r
- }\r
-\r
-}; // end class VectorImpl\r
-\r
-} // end namespace detail\r
-\r
-\r
-using detail::operator==;\r
-using detail::operator<<;\r
-\r
-class Vector : public detail::VectorImpl\r
-{\r
- public:\r
- typedef detail::VectorImpl base_type;\r
-\r
- public:\r
- Vector(size_t n)\r
- {\r
- m_size = n;\r
- m_vector = gsl_vector_alloc(n);\r
- }\r
-\r
- Vector(size_t n, double x)\r
- {\r
- m_size = n;\r
- m_vector = gsl_vector_alloc(n);\r
- gsl_vector_set_all(m_vector, x);\r
- }\r
-\r
- // create a vector with n elements all set to zero\r
- // but the i-th that is set to 1\r
- Vector(size_t n, size_t i)\r
- {\r
- m_size = n;\r
- m_vector = gsl_vector_alloc(n);\r
- gsl_vector_set_basis(m_vector, i);\r
- }\r
-\r
- Vector(Vector const& _vector)\r
- : base_type()\r
- {\r
- m_size = _vector.size();\r
- m_vector = gsl_vector_alloc(size());\r
- gsl_vector_memcpy(m_vector, _vector.m_vector);\r
- }\r
-\r
- explicit\r
- Vector(base_type::base_type const& _vector)\r
- {\r
- m_size = _vector.size();\r
- m_vector = gsl_vector_alloc(size());\r
- gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
- }\r
-\r
- virtual ~Vector()\r
- {\r
- gsl_vector_free(m_vector);\r
- }\r
-\r
-\r
- Vector & operator=(Vector const& _vector)\r
- {\r
- assert( size() == _vector.size() );\r
- gsl_vector_memcpy(m_vector, _vector.m_vector);\r
- return (*this);\r
- }\r
-\r
- Vector & operator=(base_type::base_type const& _vector)\r
- {\r
- assert( size() == _vector.size() );\r
- gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
- return (*this);\r
- }\r
-\r
- Vector & scale(double x)\r
- {\r
- return static_cast<Vector&>( base_type::scale(x) );\r
- }\r
-\r
- Vector & translate(double x)\r
- {\r
- return static_cast<Vector&>( base_type::translate(x) );\r
- }\r
-\r
- Vector & operator+=(base_type::base_type const& _vector)\r
- {\r
- return static_cast<Vector&>( base_type::operator+=(_vector) );\r
- }\r
-\r
- Vector & operator-=(base_type::base_type const& _vector)\r
- {\r
- return static_cast<Vector&>( base_type::operator-=(_vector) );\r
- }\r
-\r
- friend\r
- void swap(Vector & v1, Vector & v2);\r
- friend\r
- void swap_any(Vector & v1, Vector & v2);\r
-\r
-}; // end class Vector\r
-\r
-\r
-// warning! these operations invalidate any view of the passed vector objects\r
-inline\r
-void swap(Vector & v1, Vector & v2)\r
-{\r
- assert( v1.size() == v2.size() );\r
- std::swap(v1.m_vector, v2.m_vector);\r
-}\r
-\r
-inline\r
-void swap_any(Vector & v1, Vector & v2)\r
-{\r
- std::swap(v1.m_vector, v2.m_vector);\r
- std::swap(v1.m_size, v2.m_size);\r
-}\r
-\r
-\r
-class ConstVectorView : public detail::BaseVectorImpl\r
-{\r
- public:\r
- typedef detail::BaseVectorImpl base_type;\r
-\r
- public:\r
- ConstVectorView(const base_type & _vector, size_t n, size_t offset = 0)\r
- : m_vector_view( gsl_vector_const_subvector(_vector.get_gsl_vector(), offset, n) )\r
- {\r
- m_size = n;\r
- m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
- }\r
-\r
- ConstVectorView(const base_type & _vector, size_t n, size_t offset , size_t stride)\r
- : m_vector_view( gsl_vector_const_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n) )\r
- {\r
- m_size = n;\r
- m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
- }\r
-\r
- ConstVectorView(const double* _vector, size_t n, size_t offset = 0)\r
- : m_vector_view( gsl_vector_const_view_array(_vector + offset, n) )\r
- {\r
- m_size = n;\r
- m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
- }\r
-\r
- ConstVectorView(const double* _vector, size_t n, size_t offset, size_t stride)\r
- : m_vector_view( gsl_vector_const_view_array_with_stride(_vector + offset, stride, n) )\r
- {\r
- m_size = n;\r
- m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
- }\r
-\r
- explicit\r
- ConstVectorView(gsl_vector_const_view _gsl_vector_view)\r
- : m_vector_view(_gsl_vector_view)\r
- {\r
- m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
- m_size = m_vector->size;\r
- }\r
-\r
- explicit\r
- ConstVectorView(const std::vector<double>& _vector)\r
- : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) )\r
- {\r
- m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
- m_size = _vector.size();\r
- }\r
-\r
- ConstVectorView(const ConstVectorView & _vector)\r
- : base_type(),\r
- m_vector_view(_vector.m_vector_view)\r
- {\r
- m_size = _vector.size();\r
- m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
- }\r
-\r
- ConstVectorView(const base_type & _vector)\r
- : m_vector_view(gsl_vector_const_subvector(_vector.get_gsl_vector(), 0, _vector.size()))\r
- {\r
- m_size = _vector.size();\r
- m_vector = const_cast<gsl_vector*>( &(m_vector_view.vector) );\r
- }\r
-\r
- private:\r
- gsl_vector_const_view m_vector_view;\r
-\r
-}; // end class ConstVectorView\r
-\r
-\r
-\r
-\r
-class VectorView : public detail::VectorImpl\r
-{\r
- public:\r
- typedef detail::VectorImpl base_type;\r
-\r
- public:\r
- VectorView(base_type & _vector, size_t n, size_t offset = 0, size_t stride = 1)\r
- {\r
- m_size = n;\r
- if (stride == 1)\r
- {\r
- m_vector_view\r
- = gsl_vector_subvector(_vector.get_gsl_vector(), offset, n);\r
- m_vector = &(m_vector_view.vector);\r
- }\r
- else\r
- {\r
- m_vector_view\r
- = gsl_vector_subvector_with_stride(_vector.get_gsl_vector(), offset, stride, n);\r
- m_vector = &(m_vector_view.vector);\r
- }\r
- }\r
-\r
- VectorView(double* _vector, size_t n, size_t offset = 0, size_t stride = 1)\r
- {\r
- m_size = n;\r
- if (stride == 1)\r
- {\r
- m_vector_view\r
- = gsl_vector_view_array(_vector + offset, n);\r
- m_vector = &(m_vector_view.vector);\r
- }\r
- else\r
- {\r
- m_vector_view\r
- = gsl_vector_view_array_with_stride(_vector + offset, stride, n);\r
- m_vector = &(m_vector_view.vector);\r
- }\r
-\r
- }\r
-\r
- VectorView(const VectorView & _vector)\r
- : base_type()\r
- {\r
- m_size = _vector.size();\r
- m_vector_view = _vector.m_vector_view;\r
- m_vector = &(m_vector_view.vector);\r
- }\r
-\r
- VectorView(Vector & _vector)\r
- {\r
- m_size = _vector.size();\r
- m_vector_view = gsl_vector_subvector(_vector.get_gsl_vector(), 0, size());\r
- m_vector = &(m_vector_view.vector);\r
- }\r
-\r
- explicit\r
- VectorView(gsl_vector_view _gsl_vector_view)\r
- : m_vector_view(_gsl_vector_view)\r
- {\r
- m_vector = &(m_vector_view.vector);\r
- m_size = m_vector->size;\r
- }\r
-\r
- explicit\r
- VectorView(std::vector<double> & _vector)\r
- {\r
- m_size = _vector.size();\r
- m_vector_view = gsl_vector_view_array(&(_vector[0]), _vector.size());\r
- m_vector = &(m_vector_view.vector);\r
- }\r
-\r
- VectorView & operator=(VectorView const& _vector)\r
- {\r
- assert( size() == _vector.size() );\r
- gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
- return (*this);\r
- }\r
-\r
- VectorView & operator=(base_type::base_type const& _vector)\r
- {\r
- assert( size() == _vector.size() );\r
- gsl_vector_memcpy(m_vector, _vector.get_gsl_vector());\r
- return (*this);\r
- }\r
-\r
- VectorView & scale(double x)\r
- {\r
- return static_cast<VectorView&>( base_type::scale(x) );\r
- }\r
-\r
- VectorView & translate(double x)\r
- {\r
- return static_cast<VectorView&>( base_type::translate(x) );\r
- }\r
-\r
- VectorView & operator+=(base_type::base_type const& _vector)\r
- {\r
- return static_cast<VectorView&>( base_type::operator+=(_vector) );\r
- }\r
-\r
- VectorView & operator-=(base_type::base_type const& _vector)\r
- {\r
- return static_cast<VectorView&>( base_type::operator-=(_vector) );\r
- }\r
-\r
- friend\r
- void swap_view(VectorView & v1, VectorView & v2);\r
-\r
- private:\r
- gsl_vector_view m_vector_view;\r
-\r
-}; // end class VectorView\r
-\r
-\r
-inline\r
-void swap_view(VectorView & v1, VectorView & v2)\r
-{\r
- assert( v1.size() == v2.size() );\r
- std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too\r
-}\r
-\r
-\r
-} } // end namespaces\r
-\r
-\r
-#endif /*_NL_VECTOR_H_*/\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
+/*
+ * 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 :
index 86b919fd28eb08aac5c98d87bde9b4de7d04eb03..047ae1431585462dbc7c1984577dc741174c09c2 100644 (file)
-/*\r
- * SVG Elliptical Arc Class\r
- *\r
- * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#include <2geom/svg-elliptical-arc.h>\r
-#include <2geom/ellipse.h>\r
-#include <2geom/sbasis-geometric.h>\r
-#include <2geom/bezier-curve.h>\r
-#include <2geom/poly.h>\r
-\r
-#include <cfloat>\r
-#include <limits>\r
-\r
-#include <2geom/numeric/vector.h>\r
-#include <2geom/numeric/fitting-tool.h>\r
-#include <2geom/numeric/fitting-model.h>\r
-\r
-\r
-\r
-namespace Geom\r
-{\r
-\r
-\r
-Rect SVGEllipticalArc::boundsExact() const\r
-{\r
- if (isDegenerate() && is_svg_compliant())\r
- return chord().boundsExact();\r
-\r
- std::vector<double> extremes(4);\r
- double cosrot = std::cos(rotation_angle());\r
- double sinrot = std::sin(rotation_angle());\r
- extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );\r
- extremes[1] = extremes[0] + M_PI;\r
- if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;\r
- extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );\r
- extremes[3] = extremes[2] + M_PI;\r
- if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;\r
-\r
-\r
- std::vector<double>arc_extremes(4);\r
- arc_extremes[0] = initialPoint()[X];\r
- arc_extremes[1] = finalPoint()[X];\r
- if ( arc_extremes[0] < arc_extremes[1] )\r
- std::swap(arc_extremes[0], arc_extremes[1]);\r
- arc_extremes[2] = initialPoint()[Y];\r
- arc_extremes[3] = finalPoint()[Y];\r
- if ( arc_extremes[2] < arc_extremes[3] )\r
- std::swap(arc_extremes[2], arc_extremes[3]);\r
-\r
-\r
- if ( start_angle() < end_angle() )\r
- {\r
- if ( sweep_flag() )\r
- {\r
- for ( unsigned int i = 0; i < extremes.size(); ++i )\r
- {\r
- if ( start_angle() < extremes[i] && extremes[i] < end_angle() )\r
- {\r
- arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
- }\r
- }\r
- }\r
- else\r
- {\r
- for ( unsigned int i = 0; i < extremes.size(); ++i )\r
- {\r
- if ( start_angle() > extremes[i] || extremes[i] > end_angle() )\r
- {\r
- arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
- }\r
- }\r
- }\r
- }\r
- else\r
- {\r
- if ( sweep_flag() )\r
- {\r
- for ( unsigned int i = 0; i < extremes.size(); ++i )\r
- {\r
- if ( start_angle() < extremes[i] || extremes[i] < end_angle() )\r
- {\r
- arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
- }\r
- }\r
- }\r
- else\r
- {\r
- for ( unsigned int i = 0; i < extremes.size(); ++i )\r
- {\r
- if ( start_angle() > extremes[i] && extremes[i] > end_angle() )\r
- {\r
- arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
- }\r
- }\r
- }\r
- }\r
-\r
- return Rect( Point(arc_extremes[1], arc_extremes[3]) ,\r
- Point(arc_extremes[0], arc_extremes[2]) );\r
-}\r
-\r
-\r
-double SVGEllipticalArc::valueAtAngle(Coord t, Dim2 d) const\r
-{\r
- double sin_rot_angle = std::sin(rotation_angle());\r
- double cos_rot_angle = std::cos(rotation_angle());\r
- if ( d == X )\r
- {\r
- return ray(X) * cos_rot_angle * std::cos(t)\r
- - ray(Y) * sin_rot_angle * std::sin(t)\r
- + center(X);\r
- }\r
- else if ( d == Y )\r
- {\r
- return ray(X) * sin_rot_angle * std::cos(t)\r
- + ray(Y) * cos_rot_angle * std::sin(t)\r
- + center(Y);\r
- }\r
- THROW_RANGEERROR("dimension parameter out of range");\r
-}\r
-\r
-\r
-std::vector<double>\r
-SVGEllipticalArc::roots(double v, Dim2 d) const\r
-{\r
- if ( d > Y )\r
- {\r
- THROW_RANGEERROR("dimention out of range");\r
- }\r
-\r
- std::vector<double> sol;\r
-\r
- if (isDegenerate() && is_svg_compliant())\r
- {\r
- return chord().roots(v, d);\r
- }\r
- else\r
- {\r
- if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
- {\r
- if ( center(d) == v )\r
- sol.push_back(0);\r
- return sol;\r
- }\r
-\r
- const char* msg[2][2] =\r
- {\r
- { "d == X; ray(X) == 0; "\r
- "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "\r
- "s should be contained in [-1,1]",\r
- "d == X; ray(Y) == 0; "\r
- "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "\r
- "s should be contained in [-1,1]"\r
- },\r
- { "d == Y; ray(X) == 0; "\r
- "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "\r
- "s should be contained in [-1,1]",\r
- "d == Y; ray(Y) == 0; "\r
- "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "\r
- "s should be contained in [-1,1]"\r
- },\r
- };\r
-\r
- for ( unsigned int dim = 0; dim < 2; ++dim )\r
- {\r
- if ( are_near(ray(dim), 0) )\r
- {\r
- if ( initialPoint()[d] == v && finalPoint()[d] == v )\r
- {\r
- THROW_INFINITESOLUTIONS(0);\r
- }\r
- if ( (initialPoint()[d] < finalPoint()[d])\r
- && (initialPoint()[d] > v || finalPoint()[d] < v) )\r
- {\r
- return sol;\r
- }\r
- if ( (initialPoint()[d] > finalPoint()[d])\r
- && (finalPoint()[d] > v || initialPoint()[d] < v) )\r
- {\r
- return sol;\r
- }\r
- double ray_prj;\r
- switch(d)\r
- {\r
- case X:\r
- switch(dim)\r
- {\r
- case X: ray_prj = -ray(Y) * std::sin(rotation_angle());\r
- break;\r
- case Y: ray_prj = ray(X) * std::cos(rotation_angle());\r
- break;\r
- }\r
- break;\r
- case Y:\r
- switch(dim)\r
- {\r
- case X: ray_prj = ray(Y) * std::cos(rotation_angle());\r
- break;\r
- case Y: ray_prj = ray(X) * std::sin(rotation_angle());\r
- break;\r
- }\r
- break;\r
- }\r
-\r
- double s = (v - center(d)) / ray_prj;\r
- if ( s < -1 || s > 1 )\r
- {\r
- THROW_LOGICALERROR(msg[d][dim]);\r
- }\r
- switch(dim)\r
- {\r
- case X:\r
- s = std::asin(s); // return a value in [-PI/2,PI/2]\r
- if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) ) )\r
- {\r
- if ( s < 0 ) s += 2*M_PI;\r
- }\r
- else\r
- {\r
- s = M_PI - s;\r
- if (!(s < 2*M_PI) ) s -= 2*M_PI;\r
- }\r
- break;\r
- case Y:\r
- s = std::acos(s); // return a value in [0,PI]\r
- if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )\r
- {\r
- s = 2*M_PI - s;\r
- if ( !(s < 2*M_PI) ) s -= 2*M_PI;\r
- }\r
- break;\r
- }\r
-\r
- //std::cerr << "s = " << rad_to_deg(s);\r
- s = map_to_01(s);\r
- //std::cerr << " -> t: " << s << std::endl;\r
- if ( !(s < 0 || s > 1) )\r
- sol.push_back(s);\r
- return sol;\r
- }\r
- }\r
-\r
- }\r
-\r
- double rotx, roty;\r
- switch(d)\r
- {\r
- case X:\r
- rotx = std::cos(rotation_angle());\r
- roty = -std::sin(rotation_angle());\r
- break;\r
- case Y:\r
- rotx = std::sin(rotation_angle());\r
- roty = std::cos(rotation_angle());\r
- break;\r
- }\r
- double rxrotx = ray(X) * rotx;\r
- double c_v = center(d) - v;\r
-\r
- double a = -rxrotx + c_v;\r
- double b = ray(Y) * roty;\r
- double c = rxrotx + c_v;\r
- //std::cerr << "a = " << a << std::endl;\r
- //std::cerr << "b = " << b << std::endl;\r
- //std::cerr << "c = " << c << std::endl;\r
-\r
- if ( are_near(a,0) )\r
- {\r
- sol.push_back(M_PI);\r
- if ( !are_near(b,0) )\r
- {\r
- double s = 2 * std::atan(-c/(2*b));\r
- if ( s < 0 ) s += 2*M_PI;\r
- sol.push_back(s);\r
- }\r
- }\r
- else\r
- {\r
- double delta = b * b - a * c;\r
- //std::cerr << "delta = " << delta << std::endl;\r
- if ( are_near(delta, 0) )\r
- {\r
- double s = 2 * std::atan(-b/a);\r
- if ( s < 0 ) s += 2*M_PI;\r
- sol.push_back(s);\r
- }\r
- else if ( delta > 0 )\r
- {\r
- double sq = std::sqrt(delta);\r
- double s = 2 * std::atan( (-b - sq) / a );\r
- if ( s < 0 ) s += 2*M_PI;\r
- sol.push_back(s);\r
- s = 2 * std::atan( (-b + sq) / a );\r
- if ( s < 0 ) s += 2*M_PI;\r
- sol.push_back(s);\r
- }\r
- }\r
-\r
- std::vector<double> arc_sol;\r
- for (unsigned int i = 0; i < sol.size(); ++i )\r
- {\r
- //std::cerr << "s = " << rad_to_deg(sol[i]);\r
- sol[i] = map_to_01(sol[i]);\r
- //std::cerr << " -> t: " << sol[i] << std::endl;\r
- if ( !(sol[i] < 0 || sol[i] > 1) )\r
- arc_sol.push_back(sol[i]);\r
- }\r
- return arc_sol;\r
-}\r
-\r
-\r
-// D(E(t,C),t) = E(t+PI/2,O)\r
-Curve* SVGEllipticalArc::derivative() const\r
-{\r
- if (isDegenerate() && is_svg_compliant())\r
- return chord().derivative();\r
-\r
- SVGEllipticalArc* result = new SVGEllipticalArc(*this);\r
- result->m_center[X] = result->m_center[Y] = 0;\r
- result->m_start_angle += M_PI/2;\r
- if( !( result->m_start_angle < 2*M_PI ) )\r
- {\r
- result->m_start_angle -= 2*M_PI;\r
- }\r
- result->m_end_angle += M_PI/2;\r
- if( !( result->m_end_angle < 2*M_PI ) )\r
- {\r
- result->m_end_angle -= 2*M_PI;\r
- }\r
- result->m_initial_point = result->pointAtAngle( result->start_angle() );\r
- result->m_final_point = result->pointAtAngle( result->end_angle() );\r
- return result;\r
-}\r
-\r
-\r
-std::vector<Point>\r
-SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const\r
-{\r
- if (isDegenerate() && is_svg_compliant())\r
- return chord().pointAndDerivatives(t, n);\r
-\r
- unsigned int nn = n+1; // nn represents the size of the result vector.\r
- std::vector<Point> result;\r
- result.reserve(nn);\r
- double angle = map_unit_interval_on_circular_arc(t, start_angle(),\r
- end_angle(), sweep_flag());\r
- SVGEllipticalArc ea(*this);\r
- ea.m_center = Point(0,0);\r
- unsigned int m = std::min(nn, 4u);\r
- for ( unsigned int i = 0; i < m; ++i )\r
- {\r
- result.push_back( ea.pointAtAngle(angle) );\r
- angle += M_PI/2;\r
- if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;\r
- }\r
- m = nn / 4;\r
- for ( unsigned int i = 1; i < m; ++i )\r
- {\r
- for ( unsigned int j = 0; j < 4; ++j )\r
- result.push_back( result[j] );\r
- }\r
- m = nn - 4 * m;\r
- for ( unsigned int i = 0; i < m; ++i )\r
- {\r
- result.push_back( result[i] );\r
- }\r
- if ( !result.empty() ) // nn != 0\r
- result[0] = pointAtAngle(angle);\r
- return result;\r
-}\r
-\r
-bool SVGEllipticalArc::containsAngle(Coord angle) const\r
-{\r
- if ( sweep_flag() )\r
- if ( start_angle() < end_angle() )\r
- return ( !( angle < start_angle() || angle > end_angle() ) );\r
- else\r
- return ( !( angle < start_angle() && angle > end_angle() ) );\r
- else\r
- if ( start_angle() > end_angle() )\r
- return ( !( angle > start_angle() || angle < end_angle() ) );\r
- else\r
- return ( !( angle > start_angle() && angle < end_angle() ) );\r
-}\r
-\r
-Curve* SVGEllipticalArc::portion(double f, double t) const\r
-{\r
- if (f < 0) f = 0;\r
- if (f > 1) f = 1;\r
- if (t < 0) t = 0;\r
- if (t > 1) t = 1;\r
- if ( are_near(f, t) )\r
- {\r
- SVGEllipticalArc* arc = new SVGEllipticalArc();\r
- arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);\r
- arc->m_start_angle = arc->m_end_angle = m_start_angle;\r
- arc->m_rot_angle = m_rot_angle;\r
- arc->m_sweep = m_sweep;\r
- arc->m_large_arc = m_large_arc;\r
- }\r
- SVGEllipticalArc* arc = new SVGEllipticalArc( *this );\r
- arc->m_initial_point = pointAt(f);\r
- arc->m_final_point = pointAt(t);\r
- double sa = sweep_flag() ? sweep_angle() : -sweep_angle();\r
- arc->m_start_angle = m_start_angle + sa * f;\r
- if ( !(arc->m_start_angle < 2*M_PI) )\r
- arc->m_start_angle -= 2*M_PI;\r
- if ( arc->m_start_angle < 0 )\r
- arc->m_start_angle += 2*M_PI;\r
- arc->m_end_angle = m_start_angle + sa * t;\r
- if ( !(arc->m_end_angle < 2*M_PI) )\r
- arc->m_end_angle -= 2*M_PI;\r
- if ( arc->m_end_angle < 0 )\r
- arc->m_end_angle += 2*M_PI;\r
- if ( f > t ) arc->m_sweep = !sweep_flag();\r
- if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )\r
- arc->m_large_arc = false;\r
- return arc;\r
-}\r
-\r
-\r
-std::vector<double> SVGEllipticalArc::\r
-allNearestPoints( Point const& p, double from, double to ) const\r
-{\r
- std::vector<double> result;\r
- if (isDegenerate() && is_svg_compliant())\r
- {\r
- result.push_back( chord().nearestPoint(p, from, to) );\r
- return result;\r
- }\r
-\r
- if ( from > to ) std::swap(from, to);\r
- if ( from < 0 || to > 1 )\r
- {\r
- THROW_RANGEERROR("[from,to] interval out of range");\r
- }\r
-\r
- if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) )\r
- {\r
- result.push_back(from);\r
- return result;\r
- }\r
- else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )\r
- {\r
- LineSegment seg(pointAt(from), pointAt(to));\r
- Point np = seg.pointAt( seg.nearestPoint(p) );\r
- if ( are_near(ray(Y), 0) )\r
- {\r
- if ( are_near(rotation_angle(), M_PI/2)\r
- || are_near(rotation_angle(), 3*M_PI/2) )\r
- {\r
- result = roots(np[Y], Y);\r
- }\r
- else\r
- {\r
- result = roots(np[X], X);\r
- }\r
- }\r
- else\r
- {\r
- if ( are_near(rotation_angle(), M_PI/2)\r
- || are_near(rotation_angle(), 3*M_PI/2) )\r
- {\r
- result = roots(np[X], X);\r
- }\r
- else\r
- {\r
- result = roots(np[Y], Y);\r
- }\r
- }\r
- return result;\r
- }\r
- else if ( are_near(ray(X), ray(Y)) )\r
- {\r
- Point r = p - center();\r
- if ( are_near(r, Point(0,0)) )\r
- {\r
- THROW_INFINITESOLUTIONS(0);\r
- }\r
- // TODO: implement case r != 0\r
-// Point np = ray(X) * unit_vector(r);\r
-// std::vector<double> solX = roots(np[X],X);\r
-// std::vector<double> solY = roots(np[Y],Y);\r
-// double t;\r
-// if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))\r
-// {\r
-// t = solX[0];\r
-// }\r
-// else\r
-// {\r
-// t = solX[1];\r
-// }\r
-// if ( !(t < from || t > to) )\r
-// {\r
-// result.push_back(t);\r
-// }\r
-// else\r
-// {\r
-//\r
-// }\r
- }\r
-\r
- // solve the equation <D(E(t),t)|E(t)-p> == 0\r
- // that provides min and max distance points\r
- // on the ellipse E wrt the point p\r
- // after the substitutions:\r
- // cos(t) = (1 - s^2) / (1 + s^2)\r
- // sin(t) = 2t / (1 + s^2)\r
- // where s = tan(t/2)\r
- // we get a 4th degree equation in s\r
- /*\r
- * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +\r
- * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +\r
- * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +\r
- * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])\r
- */\r
-\r
- Point p_c = p - center();\r
- double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));\r
- double cosrot = std::cos( rotation_angle() );\r
- double sinrot = std::sin( rotation_angle() );\r
- double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);\r
- Poly coeff;\r
- coeff.resize(5);\r
- coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );\r
- coeff[3] = 2 * ( rx2_ry2 + expr1 );\r
- coeff[2] = 0;\r
- coeff[1] = 2 * ( -rx2_ry2 + expr1 );\r
- coeff[0] = -coeff[4];\r
-\r
-// for ( unsigned int i = 0; i < 5; ++i )\r
-// std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;\r
-\r
- std::vector<double> real_sol;\r
- // gsl_poly_complex_solve raises an error\r
- // if the leading coefficient is zero\r
- if ( are_near(coeff[4], 0) )\r
- {\r
- real_sol.push_back(0);\r
- if ( !are_near(coeff[3], 0) )\r
- {\r
- double sq = -coeff[1] / coeff[3];\r
- if ( sq > 0 )\r
- {\r
- double s = std::sqrt(sq);\r
- real_sol.push_back(s);\r
- real_sol.push_back(-s);\r
- }\r
- }\r
- }\r
- else\r
- {\r
- real_sol = solve_reals(coeff);\r
- }\r
-\r
- for ( unsigned int i = 0; i < real_sol.size(); ++i )\r
- {\r
- real_sol[i] = 2 * std::atan(real_sol[i]);\r
- if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;\r
- }\r
- // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0\r
- // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity\r
- if ( (real_sol.size() % 2) != 0 )\r
- {\r
- real_sol.push_back(M_PI);\r
- }\r
-\r
- double mindistsq1 = std::numeric_limits<double>::max();\r
- double mindistsq2 = std::numeric_limits<double>::max();\r
- double dsq;\r
- unsigned int mi1, mi2;\r
- for ( unsigned int i = 0; i < real_sol.size(); ++i )\r
- {\r
- dsq = distanceSq(p, pointAtAngle(real_sol[i]));\r
- if ( mindistsq1 > dsq )\r
- {\r
- mindistsq2 = mindistsq1;\r
- mi2 = mi1;\r
- mindistsq1 = dsq;\r
- mi1 = i;\r
- }\r
- else if ( mindistsq2 > dsq )\r
- {\r
- mindistsq2 = dsq;\r
- mi2 = i;\r
- }\r
- }\r
-\r
- double t = map_to_01( real_sol[mi1] );\r
- if ( !(t < from || t > to) )\r
- {\r
- result.push_back(t);\r
- }\r
-\r
- bool second_sol = false;\r
- t = map_to_01( real_sol[mi2] );\r
- if ( real_sol.size() == 4 && !(t < from || t > to) )\r
- {\r
- if ( result.empty() || are_near(mindistsq1, mindistsq2) )\r
- {\r
- result.push_back(t);\r
- second_sol = true;\r
- }\r
- }\r
-\r
- // we need to test extreme points too\r
- double dsq1 = distanceSq(p, pointAt(from));\r
- double dsq2 = distanceSq(p, pointAt(to));\r
- if ( second_sol )\r
- {\r
- if ( mindistsq2 > dsq1 )\r
- {\r
- result.clear();\r
- result.push_back(from);\r
- mindistsq2 = dsq1;\r
- }\r
- else if ( are_near(mindistsq2, dsq) )\r
- {\r
- result.push_back(from);\r
- }\r
- if ( mindistsq2 > dsq2 )\r
- {\r
- result.clear();\r
- result.push_back(to);\r
- }\r
- else if ( are_near(mindistsq2, dsq2) )\r
- {\r
- result.push_back(to);\r
- }\r
-\r
- }\r
- else\r
- {\r
- if ( result.empty() )\r
- {\r
- if ( are_near(dsq1, dsq2) )\r
- {\r
- result.push_back(from);\r
- result.push_back(to);\r
- }\r
- else if ( dsq2 > dsq1 )\r
- {\r
- result.push_back(from);\r
- }\r
- else\r
- {\r
- result.push_back(to);\r
- }\r
- }\r
- }\r
-\r
- return result;\r
-}\r
-\r
-\r
-/*\r
- * NOTE: this implementation follows Standard SVG 1.1 implementation guidelines\r
- * for elliptical arc curves. See Appendix F.6.\r
- */\r
-void SVGEllipticalArc::calculate_center_and_extreme_angles()\r
-{\r
- Point d = initialPoint() - finalPoint();\r
-\r
- if (is_svg_compliant())\r
- {\r
- if ( initialPoint() == finalPoint() )\r
- {\r
- m_rx = m_ry = m_rot_angle = m_start_angle = m_end_angle = 0;\r
- m_center = initialPoint();\r
- m_large_arc = m_sweep = false;\r
- return;\r
- }\r
-\r
- m_rx = std::fabs(m_rx);\r
- m_ry = std::fabs(m_ry);\r
-\r
- if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )\r
- {\r
- m_rx = L2(d) / 2;\r
- m_ry = 0;\r
- m_rot_angle = std::atan2(d[Y], d[X]);\r
- if (m_rot_angle < 0) m_rot_angle += 2*M_PI;\r
- m_start_angle = 0;\r
- m_end_angle = M_PI;\r
- m_center = middle_point(initialPoint(), finalPoint());\r
- m_large_arc = false;\r
- m_sweep = false;\r
- return;\r
- }\r
- }\r
- else\r
- {\r
- if ( are_near(initialPoint(), finalPoint()) )\r
- {\r
- if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
- {\r
- m_start_angle = m_end_angle = 0;\r
- m_center = initialPoint();\r
- return;\r
- }\r
- else\r
- {\r
- THROW_RANGEERROR("initial and final point are the same");\r
- }\r
- }\r
- if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
- { // but initialPoint != finalPoint\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints: "\r
- "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"\r
- );\r
- }\r
- if ( are_near(ray(Y), 0) )\r
- {\r
- Point v = initialPoint() - finalPoint();\r
- if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )\r
- {\r
- double angle = std::atan2(v[Y], v[X]);\r
- if (angle < 0) angle += 2*M_PI;\r
- if ( are_near( angle, rotation_angle() ) )\r
- {\r
- m_start_angle = 0;\r
- m_end_angle = M_PI;\r
- m_center = v/2 + finalPoint();\r
- return;\r
- }\r
- angle -= M_PI;\r
- if ( angle < 0 ) angle += 2*M_PI;\r
- if ( are_near( angle, rotation_angle() ) )\r
- {\r
- m_start_angle = M_PI;\r
- m_end_angle = 0;\r
- m_center = v/2 + finalPoint();\r
- return;\r
- }\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints: "\r
- "ray(Y) == 0 "\r
- "and slope(initialPoint - finalPoint) != rotation_angle "\r
- "and != rotation_angle + PI"\r
- );\r
- }\r
- if ( L2sq(v) > 4*ray(X)*ray(X) )\r
- {\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints: "\r
- "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"\r
- );\r
- }\r
- else\r
- {\r
- THROW_RANGEERROR(\r
- "there is infinite ellipses that satisfy the given constraints: "\r
- "ray(Y) == 0 and distance(initialPoint, finalPoint) < 2*ray(X)"\r
- );\r
- }\r
-\r
- }\r
-\r
- if ( are_near(ray(X), 0) )\r
- {\r
- Point v = initialPoint() - finalPoint();\r
- if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )\r
- {\r
- double angle = std::atan2(v[Y], v[X]);\r
- if (angle < 0) angle += 2*M_PI;\r
- double rot_angle = rotation_angle() + M_PI/2;\r
- if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;\r
- if ( are_near( angle, rot_angle ) )\r
- {\r
- m_start_angle = M_PI/2;\r
- m_end_angle = 3*M_PI/2;\r
- m_center = v/2 + finalPoint();\r
- return;\r
- }\r
- angle -= M_PI;\r
- if ( angle < 0 ) angle += 2*M_PI;\r
- if ( are_near( angle, rot_angle ) )\r
- {\r
- m_start_angle = 3*M_PI/2;\r
- m_end_angle = M_PI/2;\r
- m_center = v/2 + finalPoint();\r
- return;\r
- }\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints: "\r
- "ray(X) == 0 "\r
- "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "\r
- "and != rotation_angle + (3/2)*PI"\r
- );\r
- }\r
- if ( L2sq(v) > 4*ray(Y)*ray(Y) )\r
- {\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints: "\r
- "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"\r
- );\r
- }\r
- else\r
- {\r
- THROW_RANGEERROR(\r
- "there is infinite ellipses that satisfy the given constraints: "\r
- "ray(X) == 0 and distance(initialPoint, finalPoint) < 2*ray(Y)"\r
- );\r
- }\r
-\r
- }\r
-\r
- }\r
-\r
- double sin_rot_angle = std::sin(rotation_angle());\r
- double cos_rot_angle = std::cos(rotation_angle());\r
-\r
-\r
- Matrix m( cos_rot_angle, -sin_rot_angle,\r
- sin_rot_angle, cos_rot_angle,\r
- 0, 0 );\r
-\r
- Point p = (d / 2) * m;\r
- double rx2 = m_rx * m_rx;\r
- double ry2 = m_ry * m_ry;\r
- double rxpy = m_rx * p[Y];\r
- double rypx = m_ry * p[X];\r
- double rx2py2 = rxpy * rxpy;\r
- double ry2px2 = rypx * rypx;\r
- double num = rx2 * ry2;\r
- double den = rx2py2 + ry2px2;\r
- assert(den != 0);\r
- double rad = num / den;\r
- Point c(0,0);\r
- if (rad > 1)\r
- {\r
- rad -= 1;\r
- rad = std::sqrt(rad);\r
-\r
- if (m_large_arc == m_sweep) rad = -rad;\r
- c = rad * Point(rxpy / m_ry, -rypx / m_rx);\r
-\r
- m[1] = -m[1];\r
- m[2] = -m[2];\r
-\r
- m_center = c * m + middle_point(initialPoint(), finalPoint());\r
- }\r
- else if (rad == 1 || is_svg_compliant())\r
- {\r
- double lamda = std::sqrt(1 / rad);\r
- m_rx *= lamda;\r
- m_ry *= lamda;\r
- m_center = middle_point(initialPoint(), finalPoint());\r
- }\r
- else\r
- {\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints"\r
- );\r
- }\r
-\r
- Point sp((p[X] - c[X]) / m_rx, (p[Y] - c[Y]) / m_ry);\r
- Point ep((-p[X] - c[X]) / m_rx, (-p[Y] - c[Y]) / m_ry);\r
- Point v(1, 0);\r
- m_start_angle = angle_between(v, sp);\r
- double sweep_angle = angle_between(sp, ep);\r
- if (!m_sweep && sweep_angle > 0) sweep_angle -= 2*M_PI;\r
- if (m_sweep && sweep_angle < 0) sweep_angle += 2*M_PI;\r
-\r
- if (m_start_angle < 0) m_start_angle += 2*M_PI;\r
- m_end_angle = m_start_angle + sweep_angle;\r
- if (m_end_angle < 0) m_end_angle += 2*M_PI;\r
- if (m_end_angle >= 2*M_PI) m_end_angle -= 2*M_PI;\r
-}\r
-\r
-\r
-D2<SBasis> SVGEllipticalArc::toSBasis() const\r
-{\r
-\r
- if (isDegenerate() && is_svg_compliant())\r
- return chord().toSBasis();\r
-\r
- D2<SBasis> arc;\r
- // the interval of parametrization has to be [0,1]\r
- Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );\r
- Linear param(start_angle(), et);\r
- Coord cos_rot_angle = std::cos(rotation_angle());\r
- Coord sin_rot_angle = std::sin(rotation_angle());\r
- // order = 4 seems to be enough to get a perfect looking elliptical arc\r
- // should it be choosen in function of the arc length anyway ?\r
- // or maybe a user settable parameter: toSBasis(unsigned int order) ?\r
- SBasis arc_x = ray(X) * cos(param,4);\r
- SBasis arc_y = ray(Y) * sin(param,4);\r
- arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));\r
- arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));\r
- return arc;\r
-}\r
-\r
-\r
-Coord SVGEllipticalArc::map_to_02PI(Coord t) const\r
-{\r
- if ( sweep_flag() )\r
- {\r
- Coord angle = start_angle() + sweep_angle() * t;\r
- if ( !(angle < 2*M_PI) )\r
- angle -= 2*M_PI;\r
- return angle;\r
- }\r
- else\r
- {\r
- Coord angle = start_angle() - sweep_angle() * t;\r
- if ( angle < 0 ) angle += 2*M_PI;\r
- return angle;\r
- }\r
-}\r
-\r
-Coord SVGEllipticalArc::map_to_01(Coord angle) const\r
-{\r
- return map_circular_arc_on_unit_interval(angle, start_angle(),\r
- end_angle(), sweep_flag());\r
-}\r
-\r
-\r
-namespace detail\r
-{\r
-\r
-struct ellipse_equation\r
-{\r
- ellipse_equation(double a, double b, double c, double d, double e, double f)\r
- : A(a), B(b), C(c), D(d), E(e), F(f)\r
- {\r
- }\r
-\r
- double operator()(double x, double y) const\r
- {\r
- // A * x * x + B * x * y + C * y * y + D * x + E * y + F;\r
- return (A * x + B * y + D) * x + (C * y + E) * y + F;\r
- }\r
-\r
- double operator()(Point const& p) const\r
- {\r
- return (*this)(p[X], p[Y]);\r
- }\r
-\r
- Point normal(double x, double y) const\r
- {\r
- Point n( 2 * A * x + B * y + D, 2 * C * y + B * x + E );\r
- return unit_vector(n);\r
- }\r
-\r
- Point normal(Point const& p) const\r
- {\r
- return normal(p[X], p[Y]);\r
- }\r
-\r
- double A, B, C, D, E, F;\r
-};\r
-\r
-}\r
-\r
-make_elliptical_arc::\r
-make_elliptical_arc( SVGEllipticalArc& _ea,\r
- curve_type const& _curve,\r
- unsigned int _total_samples,\r
- double _tolerance )\r
- : ea(_ea), curve(_curve),\r
- dcurve( unitVector(derivative(curve)) ),\r
- model(), fitter(model, _total_samples),\r
- tolerance(_tolerance), tol_at_extr(tolerance/2),\r
- tol_at_center(0.1), angle_tol(0.1),\r
- initial_point(curve.at0()), final_point(curve.at1()),\r
- N(_total_samples), last(N-1), partitions(N-1), p(N),\r
- svg_compliant(true)\r
-{\r
-}\r
-\r
-bool\r
-make_elliptical_arc::\r
-bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,\r
- double e1x, double e1y, double e2 )\r
-{\r
- dist_err = std::fabs( ee(p[k]) );\r
- dist_bound = std::fabs( e1x * p[k][X] + e1y * p[k][Y] + e2 );\r
- angle_err = std::fabs( dot( dcurve(k/partitions), ee.normal(p[k]) ) );\r
- //angle_err *= angle_err;\r
- return ( dist_err > dist_bound || angle_err > angle_tol );\r
-}\r
-\r
-bool\r
-make_elliptical_arc::\r
-check_bound(double A, double B, double C, double D, double E, double F)\r
-{\r
- // check error magnitude\r
- detail::ellipse_equation ee(A, B, C, D, E, F);\r
-\r
- double e1x = (2*A + B) * tol_at_extr;\r
- double e1y = (B + 2*C) * tol_at_extr;\r
- double e2 = ((D + E) + (A + B + C) * tol_at_extr) * tol_at_extr;\r
- if ( bound_exceeded(0, ee, e1x, e1y, e2) )\r
- {\r
- print_bound_error(0);\r
- return false;\r
- }\r
- if ( bound_exceeded(0, ee, e1x, e1y, e2) )\r
- {\r
- print_bound_error(last);\r
- return false;\r
- }\r
-\r
- e1x = (2*A + B) * tolerance;\r
- e1y = (B + 2*C) * tolerance;\r
- e2 = ((D + E) + (A + B + C) * tolerance) * tolerance;\r
-// std::cerr << "e1x = " << e1x << std::endl;\r
-// std::cerr << "e1y = " << e1y << std::endl;\r
-// std::cerr << "e2 = " << e2 << std::endl;\r
-\r
- for ( unsigned int k = 1; k < last; ++k )\r
- {\r
- if ( bound_exceeded(k, ee, e1x, e1y, e2) )\r
- {\r
- print_bound_error(k);\r
- return false;\r
- }\r
- }\r
-\r
- return true;\r
-}\r
-\r
-void make_elliptical_arc::fit()\r
-{\r
- for (unsigned int k = 0; k < N; ++k)\r
- {\r
- p[k] = curve( k / partitions );\r
- fitter.append(p[k]);\r
- }\r
- fitter.update();\r
-\r
- NL::Vector z(N, 0.0);\r
- fitter.result(z);\r
-}\r
-\r
-bool make_elliptical_arc::make_elliptiarc()\r
-{\r
- const NL::Vector & coeff = fitter.result();\r
- Ellipse e;\r
- try\r
- {\r
- e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]);\r
- }\r
- catch(LogicalError exc)\r
- {\r
- return false;\r
- }\r
-\r
- Point inner_point = curve(0.5);\r
-\r
- if (svg_compliant_flag())\r
- {\r
- ea = e.arc(initial_point, inner_point, final_point);\r
- }\r
- else\r
- {\r
- try\r
- {\r
- ea = e.arc(initial_point, inner_point, final_point, false);\r
- }\r
- catch(RangeError exc)\r
- {\r
- return false;\r
- }\r
- }\r
-\r
-\r
- if ( !are_near( e.center(),\r
- ea.center(),\r
- tol_at_center * std::min(e.ray(X),e.ray(Y))\r
- )\r
- )\r
- {\r
- return false;\r
- }\r
- return true;\r
-}\r
-\r
-\r
-\r
-} // end namespace Geom\r
-\r
-\r
-\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
-\r
+/*
+ * 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 :
+
index ae6d2e2544ee08f4d2bc6455ad443b80b1624b10..f129e5a65c96193cd00f0fa8a91b431c621b954a 100644 (file)
-/*\r
- * Elliptical Arc - implementation of the svg elliptical arc path element\r
- *\r
- * Authors:\r
- * MenTaLguY <mental@rydia.net>\r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- *\r
- * Copyright 2007-2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#ifndef _2GEOM_SVG_ELLIPTICAL_ARC_H_\r
-#define _2GEOM_SVG_ELLIPTICAL_ARC_H_\r
-\r
-\r
-#include <2geom/curve.h>\r
-#include <2geom/angle.h>\r
-#include <2geom/utils.h>\r
-#include <2geom/bezier-curve.h>\r
-#include <2geom/sbasis-curve.h> // for non-native methods\r
-#include <2geom/numeric/vector.h>\r
-#include <2geom/numeric/fitting-tool.h>\r
-#include <2geom/numeric/fitting-model.h>\r
-\r
-\r
-#include <algorithm>\r
-\r
-\r
-\r
-namespace Geom\r
-{\r
-\r
-class SVGEllipticalArc : public Curve\r
-{\r
- public:\r
- SVGEllipticalArc(bool _svg_compliant = true)\r
- : m_initial_point(Point(0,0)), m_final_point(Point(0,0)),\r
- m_rx(0), m_ry(0), m_rot_angle(0),\r
- m_large_arc(true), m_sweep(true),\r
- m_svg_compliant(_svg_compliant)\r
- {\r
- m_start_angle = m_end_angle = 0;\r
- m_center = Point(0,0);\r
- }\r
-\r
- SVGEllipticalArc( Point _initial_point, double _rx, double _ry,\r
- double _rot_angle, bool _large_arc, bool _sweep,\r
- Point _final_point,\r
- bool _svg_compliant = true\r
- )\r
- : m_initial_point(_initial_point), m_final_point(_final_point),\r
- m_rx(_rx), m_ry(_ry), m_rot_angle(_rot_angle),\r
- m_large_arc(_large_arc), m_sweep(_sweep),\r
- m_svg_compliant(_svg_compliant)\r
- {\r
- calculate_center_and_extreme_angles();\r
- }\r
-\r
- void set( Point _initial_point, double _rx, double _ry,\r
- double _rot_angle, bool _large_arc, bool _sweep,\r
- Point _final_point\r
- )\r
- {\r
- m_initial_point = _initial_point;\r
- m_final_point = _final_point;\r
- m_rx = _rx;\r
- m_ry = _ry;\r
- m_rot_angle = _rot_angle;\r
- m_large_arc = _large_arc;\r
- m_sweep = _sweep;\r
- calculate_center_and_extreme_angles();\r
- }\r
-\r
- Curve* duplicate() const\r
- {\r
- return new SVGEllipticalArc(*this);\r
- }\r
-\r
- double center(unsigned int i) const\r
- {\r
- return m_center[i];\r
- }\r
-\r
- Point center() const\r
- {\r
- return m_center;\r
- }\r
-\r
- Point initialPoint() const\r
- {\r
- return m_initial_point;\r
- }\r
-\r
- Point finalPoint() const\r
- {\r
- return m_final_point;\r
- }\r
-\r
- double start_angle() const\r
- {\r
- return m_start_angle;\r
- }\r
-\r
- double end_angle() const\r
- {\r
- return m_end_angle;\r
- }\r
-\r
- double ray(unsigned int i) const\r
- {\r
- return (i == 0) ? m_rx : m_ry;\r
- }\r
-\r
- bool large_arc_flag() const\r
- {\r
- return m_large_arc;\r
- }\r
-\r
- bool sweep_flag() const\r
- {\r
- return m_sweep;\r
- }\r
-\r
- double rotation_angle() const\r
- {\r
- return m_rot_angle;\r
- }\r
-\r
- void setInitial( const Point _point)\r
- {\r
- m_initial_point = _point;\r
- calculate_center_and_extreme_angles();\r
- }\r
-\r
- void setFinal( const Point _point)\r
- {\r
- m_final_point = _point;\r
- calculate_center_and_extreme_angles();\r
- }\r
-\r
- void setExtremes( const Point& _initial_point, const Point& _final_point )\r
- {\r
- m_initial_point = _initial_point;\r
- m_final_point = _final_point;\r
- calculate_center_and_extreme_angles();\r
- }\r
-\r
- bool isDegenerate() const\r
- {\r
- return ( are_near(ray(X), 0) || are_near(ray(Y), 0) );\r
- }\r
-\r
- bool is_svg_compliant() const\r
- {\r
- return m_svg_compliant;\r
- }\r
-\r
- Rect boundsFast() const\r
- {\r
- return boundsExact();\r
- }\r
-\r
- Rect boundsExact() const;\r
-\r
- // TODO: native implementation of the following methods\r
- Rect boundsLocal(Interval i, unsigned int deg) const\r
- {\r
- if (isDegenerate() && is_svg_compliant())\r
- return chord().boundsLocal(i, deg);\r
- else\r
- return SBasisCurve(toSBasis()).boundsLocal(i, deg);\r
- }\r
-\r
- std::vector<double> roots(double v, Dim2 d) const;\r
-\r
- std::vector<double>\r
- allNearestPoints( Point const& p, double from = 0, double to = 1 ) const;\r
-\r
- double nearestPoint( Point const& p, double from = 0, double to = 1 ) const\r
- {\r
- if ( are_near(ray(X), ray(Y)) && are_near(center(), p) )\r
- {\r
- return from;\r
- }\r
- return allNearestPoints(p, from, to).front();\r
- }\r
-\r
- // TODO: native implementation of the following methods\r
- int winding(Point p) const\r
- {\r
- if (isDegenerate() && is_svg_compliant())\r
- return chord().winding(p);\r
- else\r
- return SBasisCurve(toSBasis()).winding(p);\r
- }\r
-\r
- Curve *derivative() const;\r
-\r
- // TODO: native implementation of the following methods\r
- Curve *transformed(Matrix const &m) const\r
- {\r
- return SBasisCurve(toSBasis()).transformed(m);\r
- }\r
-\r
- std::vector<Point> pointAndDerivatives(Coord t, unsigned int n) const;\r
-\r
- D2<SBasis> toSBasis() const;\r
-\r
- bool containsAngle(Coord angle) const;\r
-\r
- double valueAtAngle(Coord t, Dim2 d) const;\r
-\r
- Point pointAtAngle(Coord t) const\r
- {\r
- double sin_rot_angle = std::sin(rotation_angle());\r
- double cos_rot_angle = std::cos(rotation_angle());\r
- Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,\r
- -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,\r
- center(X), center(Y) );\r
- Point p( std::cos(t), std::sin(t) );\r
- return p * m;\r
- }\r
-\r
- double valueAt(Coord t, Dim2 d) const\r
- {\r
- if (isDegenerate() && is_svg_compliant())\r
- return chord().valueAt(t, d);\r
-\r
- Coord tt = map_to_02PI(t);\r
- return valueAtAngle(tt, d);\r
- }\r
-\r
- Point pointAt(Coord t) const\r
- {\r
- if (isDegenerate() && is_svg_compliant())\r
- return chord().pointAt(t);\r
-\r
- Coord tt = map_to_02PI(t);\r
- return pointAtAngle(tt);\r
- }\r
-\r
- std::pair<SVGEllipticalArc, SVGEllipticalArc>\r
- subdivide(Coord t) const\r
- {\r
- SVGEllipticalArc* arc1 = static_cast<SVGEllipticalArc*>(portion(0, t));\r
- SVGEllipticalArc* arc2 = static_cast<SVGEllipticalArc*>(portion(t, 1));\r
- assert( arc1 != NULL && arc2 != NULL);\r
- std::pair<SVGEllipticalArc, SVGEllipticalArc> arc_pair(*arc1, *arc2);\r
- delete arc1;\r
- delete arc2;\r
- return arc_pair;\r
- }\r
-\r
- Curve* portion(double f, double t) const;\r
-\r
- // the arc is the same but traversed in the opposite direction\r
- Curve* reverse() const\r
- {\r
- SVGEllipticalArc* rarc = new SVGEllipticalArc( *this );\r
- rarc->m_sweep = !m_sweep;\r
- rarc->m_initial_point = m_final_point;\r
- rarc->m_final_point = m_initial_point;\r
- rarc->m_start_angle = m_end_angle;\r
- rarc->m_end_angle = m_start_angle;\r
- return rarc;\r
- }\r
-\r
- double sweep_angle() const\r
- {\r
- Coord d = end_angle() - start_angle();\r
- if ( !sweep_flag() ) d = -d;\r
- if ( d < 0 )\r
- d += 2*M_PI;\r
- return d;\r
- }\r
-\r
- LineSegment chord() const\r
- {\r
- return LineSegment(initialPoint(), finalPoint());\r
- }\r
-\r
- private:\r
- Coord map_to_02PI(Coord t) const;\r
- Coord map_to_01(Coord angle) const;\r
- void calculate_center_and_extreme_angles();\r
-\r
- private:\r
- Point m_initial_point, m_final_point;\r
- double m_rx, m_ry, m_rot_angle;\r
- bool m_large_arc, m_sweep;\r
- double m_start_angle, m_end_angle;\r
- Point m_center;\r
- bool m_svg_compliant;\r
-\r
-}; // end class SVGEllipticalArc\r
-\r
-template< class charT >\r
-inline\r
-std::basic_ostream<charT> &\r
-operator<< (std::basic_ostream<charT> & os, const SVGEllipticalArc & ea)\r
-{\r
- os << "{ cx: " << ea.center(X) << ", cy: " << ea.center(Y)\r
- << ", rx: " << ea.ray(X) << ", ry: " << ea.ray(Y)\r
- << ", rot angle: " << decimal_round(rad_to_deg(ea.rotation_angle()),2)\r
- << ", start angle: " << decimal_round(rad_to_deg(ea.start_angle()),2)\r
- << ", end angle: " << decimal_round(rad_to_deg(ea.end_angle()),2)\r
- << " }";\r
-\r
- return os;\r
-}\r
-\r
-\r
-\r
-\r
-namespace detail\r
-{\r
- struct ellipse_equation;\r
-}\r
-\r
-\r
-class make_elliptical_arc\r
-{\r
- public:\r
- typedef D2<SBasis> curve_type;\r
-\r
- make_elliptical_arc( SVGEllipticalArc& _ea,\r
- curve_type const& _curve,\r
- unsigned int _total_samples,\r
- double _tolerance );\r
-\r
- private:\r
- bool bound_exceeded( unsigned int k, detail::ellipse_equation const & ee,\r
- double e1x, double e1y, double e2 );\r
-\r
- bool check_bound(double A, double B, double C, double D, double E, double F);\r
-\r
- void fit();\r
-\r
- bool make_elliptiarc();\r
-\r
- void print_bound_error(unsigned int k)\r
- {\r
- std::cerr\r
- << "tolerance error" << std::endl\r
- << "at point: " << k << std::endl\r
- << "error value: "<< dist_err << std::endl\r
- << "bound: " << dist_bound << std::endl\r
- << "angle error: " << angle_err\r
- << " (" << angle_tol << ")" << std::endl;\r
- }\r
-\r
- public:\r
- bool operator()()\r
- {\r
- const NL::Vector & coeff = fitter.result();\r
- fit();\r
- if ( !check_bound(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]) )\r
- return false;\r
- if ( !(make_elliptiarc()) ) return false;\r
- return true;\r
- }\r
-\r
- bool svg_compliant_flag() const\r
- {\r
- return svg_compliant;\r
- }\r
-\r
- void svg_compliant_on()\r
- {\r
- svg_compliant = true;\r
- }\r
-\r
- void svg_compliant_off()\r
- {\r
- svg_compliant = false;\r
- }\r
-\r
- private:\r
- SVGEllipticalArc& ea;\r
- const curve_type & curve;\r
- Piecewise<D2<SBasis> > dcurve;\r
- NL::LFMEllipse model;\r
- NL::least_squeares_fitter<NL::LFMEllipse> fitter;\r
- double tolerance, tol_at_extr, tol_at_center, angle_tol;\r
- Point initial_point, final_point;\r
- unsigned int N;\r
- unsigned int last; // N-1\r
- double partitions; // N-1\r
- std::vector<Point> p; // sample points\r
- double dist_err, dist_bound, angle_err;\r
- bool svg_compliant;\r
-};\r
-\r
-\r
-} // end namespace Geom\r
-\r
-\r
-\r
-\r
-#endif /* _2GEOM_SVG_ELLIPTICAL_ARC_H_ */\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
-\r
+/*
+ * 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 :
+