From 02c55b36fc7069d6bc02a8f1bbf2e1c58428e21a Mon Sep 17 00:00:00 2001 From: joncruz Date: Fri, 4 Jul 2008 06:38:30 +0000 Subject: [PATCH] EOL fixup --- src/2geom/ellipse.cpp | 412 +++--- src/2geom/ellipse.h | 266 ++-- src/2geom/nearest-point.cpp | 556 +++---- src/2geom/nearest-point.h | 266 ++-- src/2geom/numeric/fitting-model.h | 846 +++++------ src/2geom/numeric/fitting-tool.h | 1064 +++++++------- src/2geom/numeric/linear_system.h | 276 ++-- src/2geom/numeric/matrix.h | 1110 +++++++------- src/2geom/numeric/vector.h | 1130 +++++++-------- src/2geom/svg-elliptical-arc.cpp | 2248 ++++++++++++++--------------- src/2geom/svg-elliptical-arc.h | 870 +++++------ 11 files changed, 4522 insertions(+), 4522 deletions(-) diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp index c9c5b9ec4..20ac0bb2b 100644 --- a/src/2geom/ellipse.cpp +++ b/src/2geom/ellipse.cpp @@ -1,206 +1,206 @@ -/* - * Ellipse Curve - * - * Authors: - * Marco Cecchetti - * - * 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 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 fitter(model, sz); - - for (size_t i = 0; i < sz; ++i) - { - fitter.append(points[i]); - } - fitter.update(); - - NL::Vector z(sz, 0.0); - model.instance(*this, fitter.result(z)); -} - - -SVGEllipticalArc -Ellipse::arc(Point const& initial, Point const& inner, Point const& final, - bool _svg_compliant) -{ - Point sp_cp = initial - center(); - Point ep_cp = final - center(); - Point ip_cp = inner - center(); - - double angle1 = angle_between(sp_cp, ep_cp); - double angle2 = angle_between(sp_cp, ip_cp); - double angle3 = angle_between(ip_cp, ep_cp); - - bool large_arc_flag = true; - bool sweep_flag = true; - - if ( angle1 > 0 ) - { - if ( angle2 > 0 && angle3 > 0 ) - { - large_arc_flag = false; - sweep_flag = true; - } - else - { - large_arc_flag = true; - sweep_flag = false; - } - } - else - { - if ( angle2 < 0 && angle3 < 0 ) - { - large_arc_flag = false; - sweep_flag = false; - } - else - { - large_arc_flag = true; - sweep_flag = true; - } - } - - SVGEllipticalArc ea( initial, ray(X), ray(Y), rot_angle(), - large_arc_flag, sweep_flag, final, _svg_compliant); - return ea; -} - - -} // end namespace Geom - - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : - - +/* + * Ellipse Curve + * + * Authors: + * Marco Cecchetti + * + * 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 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 fitter(model, sz); + + for (size_t i = 0; i < sz; ++i) + { + fitter.append(points[i]); + } + fitter.update(); + + NL::Vector z(sz, 0.0); + model.instance(*this, fitter.result(z)); +} + + +SVGEllipticalArc +Ellipse::arc(Point const& initial, Point const& inner, Point const& final, + bool _svg_compliant) +{ + Point sp_cp = initial - center(); + Point ep_cp = final - center(); + Point ip_cp = inner - center(); + + double angle1 = angle_between(sp_cp, ep_cp); + double angle2 = angle_between(sp_cp, ip_cp); + double angle3 = angle_between(ip_cp, ep_cp); + + bool large_arc_flag = true; + bool sweep_flag = true; + + if ( angle1 > 0 ) + { + if ( angle2 > 0 && angle3 > 0 ) + { + large_arc_flag = false; + sweep_flag = true; + } + else + { + large_arc_flag = true; + sweep_flag = false; + } + } + else + { + if ( angle2 < 0 && angle3 < 0 ) + { + large_arc_flag = false; + sweep_flag = false; + } + else + { + large_arc_flag = true; + sweep_flag = true; + } + } + + SVGEllipticalArc ea( initial, ray(X), ray(Y), rot_angle(), + large_arc_flag, sweep_flag, final, _svg_compliant); + return ea; +} + + +} // end namespace Geom + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : + + diff --git a/src/2geom/ellipse.h b/src/2geom/ellipse.h index af8b01e78..a10f9ced0 100644 --- a/src/2geom/ellipse.h +++ b/src/2geom/ellipse.h @@ -1,133 +1,133 @@ -/* - * Ellipse Curve - * - * Authors: - * Marco Cecchetti - * - * 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 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 const& points); - - SVGEllipticalArc - arc(Point const& initial, Point const& inner, Point const& final, - bool _svg_compliant = true); - - Point center() const - { - return m_centre; - } - - Coord center(Dim2 d) const - { - return m_centre[d]; - } - - Coord ray(Dim2 d) const - { - return m_ray[d]; - } - - Coord rot_angle() const - { - return m_angle; - } - - private: - Point m_centre, m_ray; - double m_angle; -}; - - -} // end namespace Geom - - - -#endif // _2GEOM_ELLIPSE_H_ - - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : +/* + * Ellipse Curve + * + * Authors: + * Marco Cecchetti + * + * 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 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 const& points); + + SVGEllipticalArc + arc(Point const& initial, Point const& inner, Point const& final, + bool _svg_compliant = true); + + Point center() const + { + return m_centre; + } + + Coord center(Dim2 d) const + { + return m_centre[d]; + } + + Coord ray(Dim2 d) const + { + return m_ray[d]; + } + + Coord rot_angle() const + { + return m_angle; + } + + private: + Point m_centre, m_ray; + double m_angle; +}; + + +} // end namespace Geom + + + +#endif // _2GEOM_ELLIPSE_H_ + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/nearest-point.cpp b/src/2geom/nearest-point.cpp index 59d55ba32..7fef8c120 100644 --- a/src/2geom/nearest-point.cpp +++ b/src/2geom/nearest-point.cpp @@ -1,278 +1,278 @@ -/* - * nearest point routines for D2 and Piecewise> - * - * Authors: - * - * Marco Cecchetti - * - * 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 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 const& c, - D2 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 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 -all_nearest_points( Point const& p, - D2 const& c, - D2 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 result; - SBasis dd = dot(c - p, Geom::derivative(c)); - - std::vector zeros = Geom::roots(dd); - std::vector candidates; - candidates.push_back(from); - candidates.insert(candidates.end(), zeros.begin(), zeros.end()); - candidates.push_back(to); - std::vector 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 > versions - - -double nearest_point( Point const& p, - Piecewise< D2 > 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 -all_nearest_points( Point const& p, - Piecewise< D2 > 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 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 all_t; - std::vector< std::vector > all_np; - all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) ); - std::vector 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 all_nearest; - for ( unsigned int i = 0; i < all_np.size(); ++i ) - { - for ( unsigned int j = 0; j < all_np[i].size(); ++j ) - { - all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) ); - } - } - return all_nearest; -} - -} // end namespace Geom - - +/* + * nearest point routines for D2 and Piecewise> + * + * Authors: + * + * Marco Cecchetti + * + * 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 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 const& c, + D2 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 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 +all_nearest_points( Point const& p, + D2 const& c, + D2 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 result; + SBasis dd = dot(c - p, Geom::derivative(c)); + + std::vector zeros = Geom::roots(dd); + std::vector candidates; + candidates.push_back(from); + candidates.insert(candidates.end(), zeros.begin(), zeros.end()); + candidates.push_back(to); + std::vector 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 > versions + + +double nearest_point( Point const& p, + Piecewise< D2 > 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 +all_nearest_points( Point const& p, + Piecewise< D2 > 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 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 all_t; + std::vector< std::vector > all_np; + all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) ); + std::vector 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 all_nearest; + for ( unsigned int i = 0; i < all_np.size(); ++i ) + { + for ( unsigned int j = 0; j < all_np[i].size(); ++j ) + { + all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) ); + } + } + return all_nearest; +} + +} // end namespace Geom + + diff --git a/src/2geom/nearest-point.h b/src/2geom/nearest-point.h index 0b43ce67b..c8588214b 100644 --- a/src/2geom/nearest-point.h +++ b/src/2geom/nearest-point.h @@ -1,133 +1,133 @@ -/* - * nearest point routines for D2 and Piecewise> - * - * - * Authors: - * - * Marco Cecchetti - * - * 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 - -#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 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 const& c, D2 const& dc, - double from = 0, double to = 1 ); - -inline -double nearest_point( Point const& p, - D2 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 -all_nearest_points( Point const& p, - D2 const& c, D2 const& dc, - double from = 0, double to = 1 ); - -inline -std::vector -all_nearest_points( Point const& p, - D2 const& c, - double from = 0, double to = 1 ) -{ - return all_nearest_points(p, c, Geom::derivative(c), from, to); -} - - -//////////////////////////////////////////////////////////////////////////////// -// Piecewise< D2 > versions - -double nearest_point( Point const& p, - Piecewise< D2 > const& c, - double from, double to ); - -inline -double nearest_point( Point const& p, Piecewise< D2 > const& c ) -{ - return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]); -} - - -std::vector -all_nearest_points( Point const& p, - Piecewise< D2 > const& c, - double from, double to ); - -inline -std::vector -all_nearest_points( Point const& p, Piecewise< D2 > const& c ) -{ - return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]); -} - -} // end namespace Geom - - - -#endif /*_NEAREST_POINT_H_*/ +/* + * nearest point routines for D2 and Piecewise> + * + * + * Authors: + * + * Marco Cecchetti + * + * 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 + +#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 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 const& c, D2 const& dc, + double from = 0, double to = 1 ); + +inline +double nearest_point( Point const& p, + D2 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 +all_nearest_points( Point const& p, + D2 const& c, D2 const& dc, + double from = 0, double to = 1 ); + +inline +std::vector +all_nearest_points( Point const& p, + D2 const& c, + double from = 0, double to = 1 ) +{ + return all_nearest_points(p, c, Geom::derivative(c), from, to); +} + + +//////////////////////////////////////////////////////////////////////////////// +// Piecewise< D2 > versions + +double nearest_point( Point const& p, + Piecewise< D2 > const& c, + double from, double to ); + +inline +double nearest_point( Point const& p, Piecewise< D2 > const& c ) +{ + return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]); +} + + +std::vector +all_nearest_points( Point const& p, + Piecewise< D2 > const& c, + double from, double to ); + +inline +std::vector +all_nearest_points( Point const& p, Piecewise< D2 > const& c ) +{ + return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]); +} + +} // end namespace Geom + + + +#endif /*_NEAREST_POINT_H_*/ diff --git a/src/2geom/numeric/fitting-model.h b/src/2geom/numeric/fitting-model.h index 145be40e4..cc3113372 100644 --- a/src/2geom/numeric/fitting-model.h +++ b/src/2geom/numeric/fitting-model.h @@ -1,423 +1,423 @@ -/* - * Fitting Models for Geom Types - * - * Authors: - * Marco Cecchetti - * - * 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 -{ - 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 -{ - public: - LFMPoly(size_t degree) - : LFMPowerBasis(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 -{ - 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 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 -{ - 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 -{ - 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 -{ - 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 objects -class LFMD2SBasis - : public LinearFittingModel< double, Point, D2 > -{ - 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 & 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 -{ - 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 m_bc; -}; - - -// this model generates Bezier curves -template< unsigned int N > -class LFMBezierCurve - : public LinearFittingModel< double, Point, BezierCurve > -{ - 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 & 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(bx, by); - } - - private: - LFMBezier mob; -}; - -} // end namespace NL -} // end namespace Geom - - -#endif // _NL_FITTING_MODEL_H_ - - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : +/* + * Fitting Models for Geom Types + * + * Authors: + * Marco Cecchetti + * + * 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 +{ + 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 +{ + public: + LFMPoly(size_t degree) + : LFMPowerBasis(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 +{ + 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 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 +{ + 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 +{ + 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 +{ + 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 objects +class LFMD2SBasis + : public LinearFittingModel< double, Point, D2 > +{ + 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 & 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 +{ + 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 m_bc; +}; + + +// this model generates Bezier curves +template< unsigned int N > +class LFMBezierCurve + : public LinearFittingModel< double, Point, BezierCurve > +{ + 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 & 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(bx, by); + } + + private: + LFMBezier mob; +}; + +} // end namespace NL +} // end namespace Geom + + +#endif // _NL_FITTING_MODEL_H_ + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/numeric/fitting-tool.h b/src/2geom/numeric/fitting-tool.h index edacc663a..d589d86e3 100644 --- a/src/2geom/numeric/fitting-tool.h +++ b/src/2geom/numeric/fitting-tool.h @@ -1,532 +1,532 @@ -/* - * Fitting Tools - * - * Authors: - * Marco Cecchetti - * - * 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 - - -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 - : public lsf_base -{ -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 base_type; - - using base_type::m_model; - using base_type::m_psdinv_matrix; - using base_type::total_samples; - -public: - lsf_solution( 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 - - -// 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 - : public lsf_base -{ -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 base_type; - - using base_type::m_model; - using base_type::m_psdinv_matrix; - using base_type::total_samples; - -public: - lsf_solution( model_type const& _model, - size_t forecasted_samples ) - : base_type(_model, forecasted_samples), - m_solution(_model.size(), 2) - { - } - - solution_type& result(std::vector 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 const& old_sample_values, - std::vector 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 - - - - -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 - : public lsf_solution -{ - public: - typedef ModelT model_type; - typedef typename model_type::parameter_type parameter_type; - typedef typename model_type::value_type value_type; - typedef lsf_solution 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( 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 - - -// fitting tool for partially known models -template< typename ModelT> -class lsf_with_fixed_terms - : public lsf_solution -{ - public: - typedef ModelT model_type; - typedef typename model_type::parameter_type parameter_type; - typedef typename model_type::value_type value_type; - typedef lsf_solution 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( 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() - { - if (m_vector_view != NULL) - { - delete m_vector_view; - } - } - - protected: - Vector m_vector; - VectorView* m_vector_view; - -}; // end class lsf_with_fixed_terms - - -} // 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 - : public detail::lsf_with_fixed_terms -{ - public: - typedef ModelT model_type; - typedef detail::lsf_with_fixed_terms 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( model_type const& _model, - size_t forecasted_samples ) - : base_type(_model, forecasted_samples) - { - } -}; // end class least_squeares_fitter - - -template< typename ModelT> -class least_squeares_fitter - : public detail::lsf_with_fixed_terms -{ - public: - typedef ModelT model_type; - typedef detail::lsf_with_fixed_terms 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( 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 - - -template< typename ModelT> -class least_squeares_fitter - : public detail::lsf_with_fixed_terms -{ - public: - typedef ModelT model_type; - typedef detail::lsf_with_fixed_terms 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( model_type const& _model, - size_t forecasted_samples ) - : base_type(_model, forecasted_samples) - { - } - - solution_type& result(std::vector 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 - - -} // end namespace NL -} // end namespace Geom - - - -#endif // _NL_FITTING_TOOL_H_ - - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : +/* + * Fitting Tools + * + * Authors: + * Marco Cecchetti + * + * 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 + + +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 + : public lsf_base +{ +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 base_type; + + using base_type::m_model; + using base_type::m_psdinv_matrix; + using base_type::total_samples; + +public: + lsf_solution( 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 + + +// 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 + : public lsf_base +{ +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 base_type; + + using base_type::m_model; + using base_type::m_psdinv_matrix; + using base_type::total_samples; + +public: + lsf_solution( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples), + m_solution(_model.size(), 2) + { + } + + solution_type& result(std::vector 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 const& old_sample_values, + std::vector 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 + + + + +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 + : public lsf_solution +{ + public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + typedef lsf_solution 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( 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 + + +// fitting tool for partially known models +template< typename ModelT> +class lsf_with_fixed_terms + : public lsf_solution +{ + public: + typedef ModelT model_type; + typedef typename model_type::parameter_type parameter_type; + typedef typename model_type::value_type value_type; + typedef lsf_solution 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( 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() + { + if (m_vector_view != NULL) + { + delete m_vector_view; + } + } + + protected: + Vector m_vector; + VectorView* m_vector_view; + +}; // end class lsf_with_fixed_terms + + +} // 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 + : public detail::lsf_with_fixed_terms +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms 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( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples) + { + } +}; // end class least_squeares_fitter + + +template< typename ModelT> +class least_squeares_fitter + : public detail::lsf_with_fixed_terms +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms 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( 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 + + +template< typename ModelT> +class least_squeares_fitter + : public detail::lsf_with_fixed_terms +{ + public: + typedef ModelT model_type; + typedef detail::lsf_with_fixed_terms 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( model_type const& _model, + size_t forecasted_samples ) + : base_type(_model, forecasted_samples) + { + } + + solution_type& result(std::vector 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 + + +} // end namespace NL +} // end namespace Geom + + + +#endif // _NL_FITTING_TOOL_H_ + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/numeric/linear_system.h b/src/2geom/numeric/linear_system.h index 5b516c9e6..dc2a1d7e0 100644 --- a/src/2geom/numeric/linear_system.h +++ b/src/2geom/numeric/linear_system.h @@ -1,138 +1,138 @@ -/* - * LinearSystem class wraps some gsl routines for solving linear systems - * - * Authors: - * Marco Cecchetti - * - * 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 - -#include - -#include <2geom/numeric/matrix.h> -#include <2geom/numeric/vector.h> - - -namespace Geom { namespace NL { - - -class LinearSystem -{ -public: - LinearSystem(MatrixView & _matrix, VectorView & _vector) - : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) - { - } - - LinearSystem(Matrix & _matrix, Vector & _vector) - : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) - { - } - - const Vector & LU_solve() - { - assert( matrix().rows() == matrix().columns() - && matrix().rows() == vector().size() ); - int s; - gsl_permutation * p = gsl_permutation_alloc(matrix().rows()); - gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s); - gsl_linalg_LU_solve( matrix().get_gsl_matrix(), - p, - vector().get_gsl_vector(), - m_solution.get_gsl_vector() - ); - gsl_permutation_free(p); - return solution(); - } - - const Vector & SV_solve() - { - assert( matrix().rows() >= matrix().columns() - && matrix().rows() == vector().size() ); - - gsl_matrix* U = matrix().get_gsl_matrix(); - gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns()); - gsl_vector* S = gsl_vector_alloc(matrix().columns()); - gsl_vector* work = gsl_vector_alloc(matrix().columns()); - - gsl_linalg_SV_decomp( U, V, S, work ); - - gsl_vector* b = vector().get_gsl_vector(); - gsl_vector* x = m_solution.get_gsl_vector(); - - gsl_linalg_SV_solve( U, V, S, b, x); - - gsl_matrix_free(V); - gsl_vector_free(S); - gsl_vector_free(work); - - return solution(); - } - - MatrixView & matrix() - { - return m_matrix; - } - - VectorView & vector() - { - return m_vector; - } - - const Vector & solution() const - { - return m_solution; - } - -private: - MatrixView m_matrix; - VectorView m_vector; - Vector m_solution; -}; - - -} } // end namespaces - - -#endif /*_NL_LINEAR_SYSTEM_H_*/ - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : +/* + * LinearSystem class wraps some gsl routines for solving linear systems + * + * Authors: + * Marco Cecchetti + * + * 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 + +#include + +#include <2geom/numeric/matrix.h> +#include <2geom/numeric/vector.h> + + +namespace Geom { namespace NL { + + +class LinearSystem +{ +public: + LinearSystem(MatrixView & _matrix, VectorView & _vector) + : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) + { + } + + LinearSystem(Matrix & _matrix, Vector & _vector) + : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns()) + { + } + + const Vector & LU_solve() + { + assert( matrix().rows() == matrix().columns() + && matrix().rows() == vector().size() ); + int s; + gsl_permutation * p = gsl_permutation_alloc(matrix().rows()); + gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s); + gsl_linalg_LU_solve( matrix().get_gsl_matrix(), + p, + vector().get_gsl_vector(), + m_solution.get_gsl_vector() + ); + gsl_permutation_free(p); + return solution(); + } + + const Vector & SV_solve() + { + assert( matrix().rows() >= matrix().columns() + && matrix().rows() == vector().size() ); + + gsl_matrix* U = matrix().get_gsl_matrix(); + gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns()); + gsl_vector* S = gsl_vector_alloc(matrix().columns()); + gsl_vector* work = gsl_vector_alloc(matrix().columns()); + + gsl_linalg_SV_decomp( U, V, S, work ); + + gsl_vector* b = vector().get_gsl_vector(); + gsl_vector* x = m_solution.get_gsl_vector(); + + gsl_linalg_SV_solve( U, V, S, b, x); + + gsl_matrix_free(V); + gsl_vector_free(S); + gsl_vector_free(work); + + return solution(); + } + + MatrixView & matrix() + { + return m_matrix; + } + + VectorView & vector() + { + return m_vector; + } + + const Vector & solution() const + { + return m_solution; + } + +private: + MatrixView m_matrix; + VectorView m_vector; + Vector m_solution; +}; + + +} } // end namespaces + + +#endif /*_NL_LINEAR_SYSTEM_H_*/ + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h index 64557a6f1..156b6e9a2 100644 --- a/src/2geom/numeric/matrix.h +++ b/src/2geom/numeric/matrix.h @@ -1,555 +1,555 @@ -/* - * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines; - * "views" mimic the semantic of C++ references: any operation performed - * on a "view" is actually performed on the "viewed object" - * - * Authors: - * Marco Cecchetti - * - * 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 -#include // for std::pair -#include // for std::swap -#include -#include - -#include -#include - - -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 - max_index() const - { - std::pair indices; - gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); - return indices; - } - - std::pair - min_index() const - { - std::pair 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 & -operator<< (std::basic_ostream & 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( base_type::transpose() ); - } - - Matrix & scale(double x) - { - return static_cast( base_type::scale(x) ); - } - - Matrix & translate(double x) - { - return static_cast( base_type::translate(x) ); - } - - Matrix & operator+=(base_type::base_type const& _matrix) - { - return static_cast( base_type::operator+=(_matrix) ); - } - - Matrix & operator-=(base_type::base_type const& _matrix) - { - return static_cast( 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( &(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( &(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( &(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( base_type::transpose() ); - } - - MatrixView & scale(double x) - { - return static_cast( base_type::scale(x) ); - } - - MatrixView & translate(double x) - { - return static_cast( base_type::translate(x) ); - } - - MatrixView & operator+=(base_type::base_type const& _matrix) - { - return static_cast( base_type::operator+=(_matrix) ); - } - - MatrixView & operator-=(base_type::base_type const& _matrix) - { - return static_cast( base_type::operator-=(_matrix) ); - } - - friend - void swap_view(MatrixView & m1, MatrixView & m2); - - private: - gsl_matrix_view m_matrix_view; - -}; // end class MatrixView - - -inline -void swap_view(MatrixView & m1, MatrixView & m2) -{ - assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() ); - std::swap(m1.m_matrix_view, m2.m_matrix_view); -} - -Vector operator*( detail::BaseMatrixImpl const& A, - detail::BaseVectorImpl const& v ); - -Matrix operator*( detail::BaseMatrixImpl const& A, - detail::BaseMatrixImpl const& B ); - -Matrix pseudo_inverse(detail::BaseMatrixImpl const& A); - -} } // end namespaces - -#endif /*_NL_MATRIX_H_*/ - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : +/* + * Matrix, MatrixView, ConstMatrixView classes wrap the gsl matrix routines; + * "views" mimic the semantic of C++ references: any operation performed + * on a "view" is actually performed on the "viewed object" + * + * Authors: + * Marco Cecchetti + * + * 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 +#include // for std::pair +#include // for std::swap +#include +#include + +#include +#include + + +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 + max_index() const + { + std::pair indices; + gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second)); + return indices; + } + + std::pair + min_index() const + { + std::pair 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 & +operator<< (std::basic_ostream & 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( base_type::transpose() ); + } + + Matrix & scale(double x) + { + return static_cast( base_type::scale(x) ); + } + + Matrix & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + Matrix & operator+=(base_type::base_type const& _matrix) + { + return static_cast( base_type::operator+=(_matrix) ); + } + + Matrix & operator-=(base_type::base_type const& _matrix) + { + return static_cast( 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( &(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( &(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( &(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( base_type::transpose() ); + } + + MatrixView & scale(double x) + { + return static_cast( base_type::scale(x) ); + } + + MatrixView & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + MatrixView & operator+=(base_type::base_type const& _matrix) + { + return static_cast( base_type::operator+=(_matrix) ); + } + + MatrixView & operator-=(base_type::base_type const& _matrix) + { + return static_cast( base_type::operator-=(_matrix) ); + } + + friend + void swap_view(MatrixView & m1, MatrixView & m2); + + private: + gsl_matrix_view m_matrix_view; + +}; // end class MatrixView + + +inline +void swap_view(MatrixView & m1, MatrixView & m2) +{ + assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() ); + std::swap(m1.m_matrix_view, m2.m_matrix_view); +} + +Vector operator*( detail::BaseMatrixImpl const& A, + detail::BaseVectorImpl const& v ); + +Matrix operator*( detail::BaseMatrixImpl const& A, + detail::BaseMatrixImpl const& B ); + +Matrix pseudo_inverse(detail::BaseMatrixImpl const& A); + +} } // end namespaces + +#endif /*_NL_MATRIX_H_*/ + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/numeric/vector.h b/src/2geom/numeric/vector.h index 43a39a1ac..3e53405f4 100644 --- a/src/2geom/numeric/vector.h +++ b/src/2geom/numeric/vector.h @@ -1,565 +1,565 @@ -/* - * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines; - * "views" mimic the semantic of C++ references: any operation performed - * on a "view" is actually performed on the "viewed object" - * - * Authors: - * Marco Cecchetti - * - * 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 -#include // for std::swap -#include -#include -#include - - -#include -#include - - -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 & -operator<< (std::basic_ostream & 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( base_type::scale(x) ); - } - - Vector & translate(double x) - { - return static_cast( base_type::translate(x) ); - } - - Vector & operator+=(base_type::base_type const& _vector) - { - return static_cast( base_type::operator+=(_vector) ); - } - - Vector & operator-=(base_type::base_type const& _vector) - { - return static_cast( 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( &(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( &(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( &(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( &(m_vector_view.vector) ); - } - - explicit - ConstVectorView(gsl_vector_const_view _gsl_vector_view) - : m_vector_view(_gsl_vector_view) - { - m_vector = const_cast( &(m_vector_view.vector) ); - m_size = m_vector->size; - } - - explicit - ConstVectorView(const std::vector& _vector) - : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) ) - { - m_vector = const_cast( &(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( &(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( &(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 & _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( base_type::scale(x) ); - } - - VectorView & translate(double x) - { - return static_cast( base_type::translate(x) ); - } - - VectorView & operator+=(base_type::base_type const& _vector) - { - return static_cast( base_type::operator+=(_vector) ); - } - - VectorView & operator-=(base_type::base_type const& _vector) - { - return static_cast( base_type::operator-=(_vector) ); - } - - friend - void swap_view(VectorView & v1, VectorView & v2); - - private: - gsl_vector_view m_vector_view; - -}; // end class VectorView - - -inline -void swap_view(VectorView & v1, VectorView & v2) -{ - assert( v1.size() == v2.size() ); - std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too -} - - -} } // end namespaces - - -#endif /*_NL_VECTOR_H_*/ - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : +/* + * Vector, VectorView, ConstVectorView classes wrap the gsl vector routines; + * "views" mimic the semantic of C++ references: any operation performed + * on a "view" is actually performed on the "viewed object" + * + * Authors: + * Marco Cecchetti + * + * 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 +#include // for std::swap +#include +#include +#include + + +#include +#include + + +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 & +operator<< (std::basic_ostream & 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( base_type::scale(x) ); + } + + Vector & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + Vector & operator+=(base_type::base_type const& _vector) + { + return static_cast( base_type::operator+=(_vector) ); + } + + Vector & operator-=(base_type::base_type const& _vector) + { + return static_cast( 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( &(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( &(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( &(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( &(m_vector_view.vector) ); + } + + explicit + ConstVectorView(gsl_vector_const_view _gsl_vector_view) + : m_vector_view(_gsl_vector_view) + { + m_vector = const_cast( &(m_vector_view.vector) ); + m_size = m_vector->size; + } + + explicit + ConstVectorView(const std::vector& _vector) + : m_vector_view( gsl_vector_const_view_array(&(_vector[0]), _vector.size()) ) + { + m_vector = const_cast( &(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( &(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( &(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 & _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( base_type::scale(x) ); + } + + VectorView & translate(double x) + { + return static_cast( base_type::translate(x) ); + } + + VectorView & operator+=(base_type::base_type const& _vector) + { + return static_cast( base_type::operator+=(_vector) ); + } + + VectorView & operator-=(base_type::base_type const& _vector) + { + return static_cast( base_type::operator-=(_vector) ); + } + + friend + void swap_view(VectorView & v1, VectorView & v2); + + private: + gsl_vector_view m_vector_view; + +}; // end class VectorView + + +inline +void swap_view(VectorView & v1, VectorView & v2) +{ + assert( v1.size() == v2.size() ); + std::swap(v1.m_vector_view, v2.m_vector_view); // not swap m_vector too +} + + +} } // end namespaces + + +#endif /*_NL_VECTOR_H_*/ + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : diff --git a/src/2geom/svg-elliptical-arc.cpp b/src/2geom/svg-elliptical-arc.cpp index 86b919fd2..047ae1431 100644 --- a/src/2geom/svg-elliptical-arc.cpp +++ b/src/2geom/svg-elliptical-arc.cpp @@ -1,1124 +1,1124 @@ -/* - * SVG Elliptical Arc Class - * - * Copyright 2008 Marco Cecchetti - * - * 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 -#include - -#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 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::vectorarc_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 -SVGEllipticalArc::roots(double v, Dim2 d) const -{ - if ( d > Y ) - { - THROW_RANGEERROR("dimention out of range"); - } - - std::vector 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 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 -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 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 SVGEllipticalArc:: -allNearestPoints( Point const& p, double from, double to ) const -{ - std::vector 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 solX = roots(np[X],X); -// std::vector 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 == 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 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 -> 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::max(); - double mindistsq2 = std::numeric_limits::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 SVGEllipticalArc::toSBasis() const -{ - - if (isDegenerate() && is_svg_compliant()) - return chord().toSBasis(); - - D2 arc; - // the interval of parametrization has to be [0,1] - Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() ); - Linear param(start_angle(), et); - Coord cos_rot_angle = std::cos(rotation_angle()); - Coord sin_rot_angle = std::sin(rotation_angle()); - // order = 4 seems to be enough to get a perfect looking elliptical arc - // should it be choosen in function of the arc length anyway ? - // or maybe a user settable parameter: toSBasis(unsigned int order) ? - SBasis arc_x = ray(X) * cos(param,4); - SBasis arc_y = ray(Y) * sin(param,4); - arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X)); - arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y)); - return arc; -} - - -Coord SVGEllipticalArc::map_to_02PI(Coord t) const -{ - if ( sweep_flag() ) - { - Coord angle = start_angle() + sweep_angle() * t; - if ( !(angle < 2*M_PI) ) - angle -= 2*M_PI; - return angle; - } - else - { - Coord angle = start_angle() - sweep_angle() * t; - if ( angle < 0 ) angle += 2*M_PI; - return angle; - } -} - -Coord SVGEllipticalArc::map_to_01(Coord angle) const -{ - return map_circular_arc_on_unit_interval(angle, start_angle(), - end_angle(), sweep_flag()); -} - - -namespace detail -{ - -struct ellipse_equation -{ - ellipse_equation(double a, double b, double c, double d, double e, double f) - : A(a), B(b), C(c), D(d), E(e), F(f) - { - } - - double operator()(double x, double y) const - { - // A * x * x + B * x * y + C * y * y + D * x + E * y + F; - return (A * x + B * y + D) * x + (C * y + E) * y + F; - } - - double operator()(Point const& p) const - { - return (*this)(p[X], p[Y]); - } - - Point normal(double x, double y) const - { - Point n( 2 * A * x + B * y + D, 2 * C * y + B * x + E ); - return unit_vector(n); - } - - Point normal(Point const& p) const - { - return normal(p[X], p[Y]); - } - - double A, B, C, D, E, F; -}; - -} - -make_elliptical_arc:: -make_elliptical_arc( SVGEllipticalArc& _ea, - curve_type const& _curve, - unsigned int _total_samples, - double _tolerance ) - : ea(_ea), curve(_curve), - dcurve( unitVector(derivative(curve)) ), - model(), fitter(model, _total_samples), - tolerance(_tolerance), tol_at_extr(tolerance/2), - tol_at_center(0.1), angle_tol(0.1), - initial_point(curve.at0()), final_point(curve.at1()), - N(_total_samples), last(N-1), partitions(N-1), p(N), - svg_compliant(true) -{ -} - -bool -make_elliptical_arc:: -bound_exceeded( unsigned int k, detail::ellipse_equation const & ee, - double e1x, double e1y, double e2 ) -{ - dist_err = std::fabs( ee(p[k]) ); - dist_bound = std::fabs( e1x * p[k][X] + e1y * p[k][Y] + e2 ); - angle_err = std::fabs( dot( dcurve(k/partitions), ee.normal(p[k]) ) ); - //angle_err *= angle_err; - return ( dist_err > dist_bound || angle_err > angle_tol ); -} - -bool -make_elliptical_arc:: -check_bound(double A, double B, double C, double D, double E, double F) -{ - // check error magnitude - detail::ellipse_equation ee(A, B, C, D, E, F); - - double e1x = (2*A + B) * tol_at_extr; - double e1y = (B + 2*C) * tol_at_extr; - double e2 = ((D + E) + (A + B + C) * tol_at_extr) * tol_at_extr; - if ( bound_exceeded(0, ee, e1x, e1y, e2) ) - { - print_bound_error(0); - return false; - } - if ( bound_exceeded(0, ee, e1x, e1y, e2) ) - { - print_bound_error(last); - return false; - } - - e1x = (2*A + B) * tolerance; - e1y = (B + 2*C) * tolerance; - e2 = ((D + E) + (A + B + C) * tolerance) * tolerance; -// std::cerr << "e1x = " << e1x << std::endl; -// std::cerr << "e1y = " << e1y << std::endl; -// std::cerr << "e2 = " << e2 << std::endl; - - for ( unsigned int k = 1; k < last; ++k ) - { - if ( bound_exceeded(k, ee, e1x, e1y, e2) ) - { - print_bound_error(k); - return false; - } - } - - return true; -} - -void make_elliptical_arc::fit() -{ - for (unsigned int k = 0; k < N; ++k) - { - p[k] = curve( k / partitions ); - fitter.append(p[k]); - } - fitter.update(); - - NL::Vector z(N, 0.0); - fitter.result(z); -} - -bool make_elliptical_arc::make_elliptiarc() -{ - const NL::Vector & coeff = fitter.result(); - Ellipse e; - try - { - e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]); - } - catch(LogicalError exc) - { - return false; - } - - Point inner_point = curve(0.5); - - if (svg_compliant_flag()) - { - ea = e.arc(initial_point, inner_point, final_point); - } - else - { - try - { - ea = e.arc(initial_point, inner_point, final_point, false); - } - catch(RangeError exc) - { - return false; - } - } - - - if ( !are_near( e.center(), - ea.center(), - tol_at_center * std::min(e.ray(X),e.ray(Y)) - ) - ) - { - return false; - } - return true; -} - - - -} // end namespace Geom - - - - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : - +/* + * SVG Elliptical Arc Class + * + * Copyright 2008 Marco Cecchetti + * + * 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 +#include + +#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 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::vectorarc_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 +SVGEllipticalArc::roots(double v, Dim2 d) const +{ + if ( d > Y ) + { + THROW_RANGEERROR("dimention out of range"); + } + + std::vector 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 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 +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 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 SVGEllipticalArc:: +allNearestPoints( Point const& p, double from, double to ) const +{ + std::vector 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 solX = roots(np[X],X); +// std::vector 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 == 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 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 -> 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::max(); + double mindistsq2 = std::numeric_limits::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 SVGEllipticalArc::toSBasis() const +{ + + if (isDegenerate() && is_svg_compliant()) + return chord().toSBasis(); + + D2 arc; + // the interval of parametrization has to be [0,1] + Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() ); + Linear param(start_angle(), et); + Coord cos_rot_angle = std::cos(rotation_angle()); + Coord sin_rot_angle = std::sin(rotation_angle()); + // order = 4 seems to be enough to get a perfect looking elliptical arc + // should it be choosen in function of the arc length anyway ? + // or maybe a user settable parameter: toSBasis(unsigned int order) ? + SBasis arc_x = ray(X) * cos(param,4); + SBasis arc_y = ray(Y) * sin(param,4); + arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X)); + arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y)); + return arc; +} + + +Coord SVGEllipticalArc::map_to_02PI(Coord t) const +{ + if ( sweep_flag() ) + { + Coord angle = start_angle() + sweep_angle() * t; + if ( !(angle < 2*M_PI) ) + angle -= 2*M_PI; + return angle; + } + else + { + Coord angle = start_angle() - sweep_angle() * t; + if ( angle < 0 ) angle += 2*M_PI; + return angle; + } +} + +Coord SVGEllipticalArc::map_to_01(Coord angle) const +{ + return map_circular_arc_on_unit_interval(angle, start_angle(), + end_angle(), sweep_flag()); +} + + +namespace detail +{ + +struct ellipse_equation +{ + ellipse_equation(double a, double b, double c, double d, double e, double f) + : A(a), B(b), C(c), D(d), E(e), F(f) + { + } + + double operator()(double x, double y) const + { + // A * x * x + B * x * y + C * y * y + D * x + E * y + F; + return (A * x + B * y + D) * x + (C * y + E) * y + F; + } + + double operator()(Point const& p) const + { + return (*this)(p[X], p[Y]); + } + + Point normal(double x, double y) const + { + Point n( 2 * A * x + B * y + D, 2 * C * y + B * x + E ); + return unit_vector(n); + } + + Point normal(Point const& p) const + { + return normal(p[X], p[Y]); + } + + double A, B, C, D, E, F; +}; + +} + +make_elliptical_arc:: +make_elliptical_arc( SVGEllipticalArc& _ea, + curve_type const& _curve, + unsigned int _total_samples, + double _tolerance ) + : ea(_ea), curve(_curve), + dcurve( unitVector(derivative(curve)) ), + model(), fitter(model, _total_samples), + tolerance(_tolerance), tol_at_extr(tolerance/2), + tol_at_center(0.1), angle_tol(0.1), + initial_point(curve.at0()), final_point(curve.at1()), + N(_total_samples), last(N-1), partitions(N-1), p(N), + svg_compliant(true) +{ +} + +bool +make_elliptical_arc:: +bound_exceeded( unsigned int k, detail::ellipse_equation const & ee, + double e1x, double e1y, double e2 ) +{ + dist_err = std::fabs( ee(p[k]) ); + dist_bound = std::fabs( e1x * p[k][X] + e1y * p[k][Y] + e2 ); + angle_err = std::fabs( dot( dcurve(k/partitions), ee.normal(p[k]) ) ); + //angle_err *= angle_err; + return ( dist_err > dist_bound || angle_err > angle_tol ); +} + +bool +make_elliptical_arc:: +check_bound(double A, double B, double C, double D, double E, double F) +{ + // check error magnitude + detail::ellipse_equation ee(A, B, C, D, E, F); + + double e1x = (2*A + B) * tol_at_extr; + double e1y = (B + 2*C) * tol_at_extr; + double e2 = ((D + E) + (A + B + C) * tol_at_extr) * tol_at_extr; + if ( bound_exceeded(0, ee, e1x, e1y, e2) ) + { + print_bound_error(0); + return false; + } + if ( bound_exceeded(0, ee, e1x, e1y, e2) ) + { + print_bound_error(last); + return false; + } + + e1x = (2*A + B) * tolerance; + e1y = (B + 2*C) * tolerance; + e2 = ((D + E) + (A + B + C) * tolerance) * tolerance; +// std::cerr << "e1x = " << e1x << std::endl; +// std::cerr << "e1y = " << e1y << std::endl; +// std::cerr << "e2 = " << e2 << std::endl; + + for ( unsigned int k = 1; k < last; ++k ) + { + if ( bound_exceeded(k, ee, e1x, e1y, e2) ) + { + print_bound_error(k); + return false; + } + } + + return true; +} + +void make_elliptical_arc::fit() +{ + for (unsigned int k = 0; k < N; ++k) + { + p[k] = curve( k / partitions ); + fitter.append(p[k]); + } + fitter.update(); + + NL::Vector z(N, 0.0); + fitter.result(z); +} + +bool make_elliptical_arc::make_elliptiarc() +{ + const NL::Vector & coeff = fitter.result(); + Ellipse e; + try + { + e.set(1, coeff[0], coeff[1], coeff[2], coeff[3], coeff[4]); + } + catch(LogicalError exc) + { + return false; + } + + Point inner_point = curve(0.5); + + if (svg_compliant_flag()) + { + ea = e.arc(initial_point, inner_point, final_point); + } + else + { + try + { + ea = e.arc(initial_point, inner_point, final_point, false); + } + catch(RangeError exc) + { + return false; + } + } + + + if ( !are_near( e.center(), + ea.center(), + tol_at_center * std::min(e.ray(X),e.ray(Y)) + ) + ) + { + return false; + } + return true; +} + + + +} // end namespace Geom + + + + +/* + Local Variables: + mode:c++ + c-file-style:"stroustrup" + c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) + indent-tabs-mode:nil + fill-column:99 + End: +*/ +// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : + diff --git a/src/2geom/svg-elliptical-arc.h b/src/2geom/svg-elliptical-arc.h index ae6d2e254..f129e5a65 100644 --- a/src/2geom/svg-elliptical-arc.h +++ b/src/2geom/svg-elliptical-arc.h @@ -1,435 +1,435 @@ -/* - * Elliptical Arc - implementation of the svg elliptical arc path element - * - * Authors: - * MenTaLguY - * Marco Cecchetti - * - * 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 - - - -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 roots(double v, Dim2 d) const; - - std::vector - 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 pointAndDerivatives(Coord t, unsigned int n) const; - - D2 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 - subdivide(Coord t) const - { - SVGEllipticalArc* arc1 = static_cast(portion(0, t)); - SVGEllipticalArc* arc2 = static_cast(portion(t, 1)); - assert( arc1 != NULL && arc2 != NULL); - std::pair 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 & -operator<< (std::basic_ostream & 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 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 > dcurve; - NL::LFMEllipse model; - NL::least_squeares_fitter 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 p; // sample points - double dist_err, dist_bound, angle_err; - bool svg_compliant; -}; - - -} // end namespace Geom - - - - -#endif /* _2GEOM_SVG_ELLIPTICAL_ARC_H_ */ - -/* - Local Variables: - mode:c++ - c-file-style:"stroustrup" - c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +)) - indent-tabs-mode:nil - fill-column:99 - End: -*/ -// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 : - +/* + * Elliptical Arc - implementation of the svg elliptical arc path element + * + * Authors: + * MenTaLguY + * Marco Cecchetti + * + * 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 + + + +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 roots(double v, Dim2 d) const; + + std::vector + 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 pointAndDerivatives(Coord t, unsigned int n) const; + + D2 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 + subdivide(Coord t) const + { + SVGEllipticalArc* arc1 = static_cast(portion(0, t)); + SVGEllipticalArc* arc2 = static_cast(portion(t, 1)); + assert( arc1 != NULL && arc2 != NULL); + std::pair 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 & +operator<< (std::basic_ostream & 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 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 > dcurve; + NL::LFMEllipse model; + NL::least_squeares_fitter 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 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 : + -- 2.30.2