From d02f89424319afe28fe2cc6696c28108448de3ac Mon Sep 17 00:00:00 2001 From: mcecchetti Date: Tue, 20 May 2008 22:29:23 +0000 Subject: [PATCH] synchronization with 2geom library --- src/2geom/CMakeLists.txt | 2 +- src/2geom/Makefile_insert | 4 +- src/2geom/choose.h | 9 + ...-elliptical-arc.cpp => elliptical-arc.cpp} | 1850 ++++++++--------- src/2geom/exception.h | 7 + src/2geom/forward.h | 2 +- src/2geom/nearest-point.cpp | 556 ++--- src/2geom/nearest-point.h | 266 +-- src/2geom/numeric/linear_system.h | 178 +- src/2geom/numeric/matrix.h | 414 ++-- src/2geom/numeric/vector.h | 369 ++-- src/2geom/path.h | 80 +- src/2geom/poly.cpp | 1 + src/2geom/sbasis.cpp | 92 +- src/2geom/sbasis.h | 34 +- src/2geom/solve-bezier-one-d.cpp | 170 +- src/2geom/svg-path.cpp | 4 +- src/2geom/svg-path.h | 2 +- 18 files changed, 2068 insertions(+), 1972 deletions(-) rename src/2geom/{svg-elliptical-arc.cpp => elliptical-arc.cpp} (91%) diff --git a/src/2geom/CMakeLists.txt b/src/2geom/CMakeLists.txt index 309d4cdd9..ac019419d 100644 --- a/src/2geom/CMakeLists.txt +++ b/src/2geom/CMakeLists.txt @@ -20,6 +20,7 @@ crossing.h d2.h d2-sbasis.cpp d2-sbasis.h +elliptical-arc.cpp exception.h forward.h geom.cpp @@ -75,7 +76,6 @@ solve-bezier-one-d.cpp solve-bezier-parametric.cpp solver.h sturm.h -svg-elliptical-arc.cpp svg-path.cpp svg-path.h svg-path-parser.cpp diff --git a/src/2geom/Makefile_insert b/src/2geom/Makefile_insert index 6f2a3e642..e2c0d6aa3 100644 --- a/src/2geom/Makefile_insert +++ b/src/2geom/Makefile_insert @@ -27,6 +27,7 @@ 2geom/d2-sbasis.h \ 2geom/d2.h \ 2geom/exception.h \ + 2geom/elliptical-arc.cpp \ 2geom/forward.h \ 2geom/geom.cpp \ 2geom/geom.h \ @@ -82,8 +83,7 @@ 2geom/solve-bezier-one-d.cpp \ 2geom/solve-bezier-parametric.cpp \ 2geom/solver.h \ - 2geom/sturm.h \ - 2geom/svg-elliptical-arc.cpp \ + 2geom/sturm.h \ 2geom/svg-path.cpp \ 2geom/svg-path.h \ 2geom/svg-path-parser.cpp \ diff --git a/src/2geom/choose.h b/src/2geom/choose.h index 46c5b2fb4..f3a8b0f4f 100644 --- a/src/2geom/choose.h +++ b/src/2geom/choose.h @@ -64,6 +64,15 @@ T choose(unsigned n, unsigned k) { return pascals_triangle[row+k]; } +// Is it faster to store them or compute them on demand? +/*template +T choose(unsigned n, unsigned k) { + T r = 1; + for(unsigned i = 1; i <= k; i++) + r = (r*(n-k+i))/i; + return r; + }*/ + #endif /* diff --git a/src/2geom/svg-elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp similarity index 91% rename from src/2geom/svg-elliptical-arc.cpp rename to src/2geom/elliptical-arc.cpp index b3420fba6..b18991c88 100644 --- a/src/2geom/svg-elliptical-arc.cpp +++ b/src/2geom/elliptical-arc.cpp @@ -1,931 +1,919 @@ -/* - * 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 "path.h" -#include "angle.h" - -#include - -#include - - - - -namespace Geom -{ - - -Rect SVGEllipticalArc::boundsExact() const -{ - 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]) ); - -} - - -std::vector -SVGEllipticalArc::roots(double v, Dim2 d) const -{ - if ( d > Y ) - { - THROW_RANGEERROR("dimention out of range"); - } - - std::vector sol; - 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_EXCEPTION("infinite solutions"); - } - 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; - - -// return SBasisCurve(toSBasis()).roots(v, d); -} - -// D(E(t,C),t) = E(t+PI/2,O) -Curve* SVGEllipticalArc::derivative() const -{ - 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 -{ - std::vector result; - result.reserve(n); - 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(n, 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 = n / 4; - for ( unsigned int i = 1; i < m; ++i ) - { - for ( unsigned int j = 0; j < 4; ++j ) - result.push_back( result[j] ); - } - m = n - 4 * m; - for ( unsigned int i = 0; i < m; ++i ) - { - result.push_back( result[i] ); - } - if ( !result.empty() ) // n != 0 - result[0] = pointAtAngle(angle); - return result; -} - -D2 SVGEllipticalArc::toSBasis() const -{ - // 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); - D2 arc; - 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; -} - - -bool SVGEllipticalArc::containsAngle(Coord angle) const -{ - if ( sweep_flag() ) - if ( start_angle() < end_angle() ) - if ( !( angle < start_angle() || angle > end_angle() ) ) - return true; - else - return false; - else - if ( !( angle < start_angle() && angle > end_angle() ) ) - return true; - else - return false; - else - if ( start_angle() > end_angle() ) - if ( !( angle > start_angle() || angle < end_angle() ) ) - return true; - else - return false; - else - if ( !( angle > start_angle() && angle < end_angle() ) ) - return true; - else - return false; -} - - -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"); -} - - -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; -} - -// NOTE: doesn't work with 360 deg arcs -void SVGEllipticalArc::calculate_center_and_extreme_angles() -{ - 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()); - - Point sp = sweep_flag() ? initialPoint() : finalPoint(); - Point ep = sweep_flag() ? finalPoint() : initialPoint(); - - Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle, - -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle, - 0, 0 ); - Matrix im = m.inverse(); - Point sol = (ep - sp) * im; - double half_sum_angle = std::atan2(-sol[X], sol[Y]); - double half_diff_angle; - if ( are_near(std::fabs(half_sum_angle), M_PI/2) ) - { - double anti_sgn_hsa = (half_sum_angle > 0) ? -1 : 1; - double arg = anti_sgn_hsa * sol[X] / 2; - // if |arg| is a little bit > 1 acos returns nan - if ( are_near(arg, 1) ) - half_diff_angle = 0; - else if ( are_near(arg, -1) ) - half_diff_angle = M_PI; - else - { - if ( !(-1 < arg && arg < 1) ) - THROW_RANGEERROR( - "there is no ellipse that satisfies the given constraints" - ); - // assert( -1 < arg && arg < 1 ); - // if it fails - // => there is no ellipse that satisfies the given constraints - half_diff_angle = std::acos( arg ); - } - - half_diff_angle = M_PI/2 - half_diff_angle; - } - else - { - double arg = sol[Y] / ( 2 * std::cos(half_sum_angle) ); - // if |arg| is a little bit > 1 asin returns nan - if ( are_near(arg, 1) ) - half_diff_angle = M_PI/2; - else if ( are_near(arg, -1) ) - half_diff_angle = -M_PI/2; - else - { - if ( !(-1 < arg && arg < 1) ) - THROW_RANGEERROR( - "there is no ellipse that satisfies the given constraints" - ); - // assert( -1 < arg && arg < 1 ); - // if it fails - // => there is no ellipse that satisfies the given constraints - half_diff_angle = std::asin( arg ); - } - } - - if ( ( m_large_arc && half_diff_angle > 0 ) - || (!m_large_arc && half_diff_angle < 0 ) ) - { - half_diff_angle = -half_diff_angle; - } - if ( half_sum_angle < 0 ) half_sum_angle += 2*M_PI; - if ( half_diff_angle < 0 ) half_diff_angle += M_PI; - - m_start_angle = half_sum_angle - half_diff_angle; - m_end_angle = half_sum_angle + half_diff_angle; - // 0 <= m_start_angle, m_end_angle < 2PI - if ( m_start_angle < 0 ) m_start_angle += 2*M_PI; - if( !(m_end_angle < 2*M_PI) ) m_end_angle -= 2*M_PI; - sol[0] = std::cos(m_start_angle); - sol[1] = std::sin(m_start_angle); - m_center = sp - sol * m; - if ( !sweep_flag() ) - { - double angle = m_start_angle; - m_start_angle = m_end_angle; - m_end_angle = angle; - } -} - -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()); -} - - -std::vector SVGEllipticalArc:: -allNearestPoints( Point const& p, double from, double to ) const -{ - if ( from > to ) std::swap(from, to); - if ( from < 0 || to > 1 ) - { - THROW_RANGEERROR("[from,to] interval out of range"); - } - std::vector result; - 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_EXCEPTION("infinite nearest points"); - } - // 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); - double coeff[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 - { - double sol[8]; - gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc(5); - gsl_poly_complex_solve(coeff, 5, w, sol ); - gsl_poly_complex_workspace_free(w); - - for ( unsigned int i = 0; i < 4; ++i ) - { - if ( sol[2*i+1] == 0 ) real_sol.push_back(sol[2*i]); - } - } - - 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; -} - - -} // 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 "path.h" +#include "angle.h" + +#include + +#include + + + + +namespace Geom +{ + + +Rect EllipticalArc::boundsExact() const +{ + 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]) ); + +} + + +std::vector +EllipticalArc::roots(double v, Dim2 d) const +{ + if ( d > Y ) + { + THROW_RANGEERROR("dimention out of range"); + } + + std::vector sol; + 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_EXCEPTION("infinite solutions"); + } + 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; + + +// return SBasisCurve(toSBasis()).roots(v, d); +} + +// D(E(t,C),t) = E(t+PI/2,O) +Curve* EllipticalArc::derivative() const +{ + EllipticalArc* result = new EllipticalArc(*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 +EllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const +{ + std::vector result; + result.reserve(n); + double angle = map_unit_interval_on_circular_arc(t, start_angle(), + end_angle(), sweep_flag()); + EllipticalArc ea(*this); + ea.m_center = Point(0,0); + unsigned int m = std::min(n, 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 = n / 4; + for ( unsigned int i = 1; i < m; ++i ) + { + for ( unsigned int j = 0; j < 4; ++j ) + result.push_back( result[j] ); + } + m = n - 4 * m; + for ( unsigned int i = 0; i < m; ++i ) + { + result.push_back( result[i] ); + } + if ( !result.empty() ) // n != 0 + result[0] = pointAtAngle(angle); + return result; +} + +D2 EllipticalArc::toSBasis() const +{ + // 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); + D2 arc; + 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; +} + + +bool EllipticalArc::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() ) ); +} + + +double EllipticalArc::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"); +} + + +Curve* EllipticalArc::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) ) + { + EllipticalArc* arc = new EllipticalArc(); + 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; + } + EllipticalArc* arc = new EllipticalArc( *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; +} + +// NOTE: doesn't work with 360 deg arcs +void EllipticalArc::calculate_center_and_extreme_angles() +{ + 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()); + + Point sp = sweep_flag() ? initialPoint() : finalPoint(); + Point ep = sweep_flag() ? finalPoint() : initialPoint(); + + Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle, + -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle, + 0, 0 ); + Matrix im = m.inverse(); + Point sol = (ep - sp) * im; + double half_sum_angle = std::atan2(-sol[X], sol[Y]); + double half_diff_angle; + if ( are_near(std::fabs(half_sum_angle), M_PI/2) ) + { + double anti_sgn_hsa = (half_sum_angle > 0) ? -1 : 1; + double arg = anti_sgn_hsa * sol[X] / 2; + // if |arg| is a little bit > 1 acos returns nan + if ( are_near(arg, 1) ) + half_diff_angle = 0; + else if ( are_near(arg, -1) ) + half_diff_angle = M_PI; + else + { + if ( !(-1 < arg && arg < 1) ) + THROW_RANGEERROR( + "there is no ellipse that satisfies the given constraints" + ); + // assert( -1 < arg && arg < 1 ); + // if it fails + // => there is no ellipse that satisfies the given constraints + half_diff_angle = std::acos( arg ); + } + + half_diff_angle = M_PI/2 - half_diff_angle; + } + else + { + double arg = sol[Y] / ( 2 * std::cos(half_sum_angle) ); + // if |arg| is a little bit > 1 asin returns nan + if ( are_near(arg, 1) ) + half_diff_angle = M_PI/2; + else if ( are_near(arg, -1) ) + half_diff_angle = -M_PI/2; + else + { + if ( !(-1 < arg && arg < 1) ) + THROW_RANGEERROR( + "there is no ellipse that satisfies the given constraints" + ); + // assert( -1 < arg && arg < 1 ); + // if it fails + // => there is no ellipse that satisfies the given constraints + half_diff_angle = std::asin( arg ); + } + } + + if ( ( m_large_arc && half_diff_angle > 0 ) + || (!m_large_arc && half_diff_angle < 0 ) ) + { + half_diff_angle = -half_diff_angle; + } + if ( half_sum_angle < 0 ) half_sum_angle += 2*M_PI; + if ( half_diff_angle < 0 ) half_diff_angle += M_PI; + + m_start_angle = half_sum_angle - half_diff_angle; + m_end_angle = half_sum_angle + half_diff_angle; + // 0 <= m_start_angle, m_end_angle < 2PI + if ( m_start_angle < 0 ) m_start_angle += 2*M_PI; + if( !(m_end_angle < 2*M_PI) ) m_end_angle -= 2*M_PI; + sol[0] = std::cos(m_start_angle); + sol[1] = std::sin(m_start_angle); + m_center = sp - sol * m; + if ( !sweep_flag() ) + { + double angle = m_start_angle; + m_start_angle = m_end_angle; + m_end_angle = angle; + } +} + +Coord EllipticalArc::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 EllipticalArc::map_to_01(Coord angle) const +{ + return map_circular_arc_on_unit_interval(angle, start_angle(), + end_angle(), sweep_flag()); +} + + +std::vector EllipticalArc:: +allNearestPoints( Point const& p, double from, double to ) const +{ + if ( from > to ) std::swap(from, to); + if ( from < 0 || to > 1 ) + { + THROW_RANGEERROR("[from,to] interval out of range"); + } + std::vector result; + 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_EXCEPTION("infinite nearest points"); + } + // 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); + double coeff[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 + { + double sol[8]; + gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc(5); + gsl_poly_complex_solve(coeff, 5, w, sol ); + gsl_poly_complex_workspace_free(w); + + for ( unsigned int i = 0; i < 4; ++i ) + { + if ( sol[2*i+1] == 0 ) real_sol.push_back(sol[2*i]); + } + } + + 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; +} + + +} // 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/exception.h b/src/2geom/exception.h index 2749ec63a..6151f7a57 100644 --- a/src/2geom/exception.h +++ b/src/2geom/exception.h @@ -99,6 +99,13 @@ public: }; #define THROW_NOTINVERTIBLE(i) throw(NotInvertible(__FILE__, __LINE__)) +class InfiniteSolutions : public RangeError { +public: + InfiniteSolutions(const char *file, const int line) + : RangeError("There are infinite solutions", file, line) {} +}; +#define THROW_INFINITESOLUTIONS(i) throw(InfiniteSolutions(__FILE__, __LINE__)) + class ContinuityError : public RangeError { public: ContinuityError(const char *file, const int line) diff --git a/src/2geom/forward.h b/src/2geom/forward.h index a34384023..454f9bf15 100644 --- a/src/2geom/forward.h +++ b/src/2geom/forward.h @@ -78,7 +78,7 @@ typedef D2 Rect; class Shape; class Region; -class SVGEllipticalArc; +class EllipticalArc; class SVGPathSink; template class SVGPathGenerator; diff --git a/src/2geom/nearest-point.cpp b/src/2geom/nearest-point.cpp index 074de1cb3..de880713f 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 "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 "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 73ac0c3ce..8c1a6f7d7 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 "d2.h" -#include "piecewise.h" -#include "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 "d2.h" +#include "piecewise.h" +#include "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/linear_system.h b/src/2geom/numeric/linear_system.h index 2b4963998..9cb521eb2 100644 --- a/src/2geom/numeric/linear_system.h +++ b/src/2geom/numeric/linear_system.h @@ -1,89 +1,89 @@ -#ifndef _NL_LINEAR_SYSTEM_H_ -#define _NL_LINEAR_SYSTEM_H_ - - -#include - -#include - -#include "numeric/matrix.h" -#include "numeric/vector.h" - - -namespace Geom { namespace NL { - - -class LinearSystem -{ -public: - 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(); - } - - Matrix & matrix() - { - return m_matrix; - } - - Vector & vector() - { - return m_vector; - } - - const Vector & solution() const - { - return m_solution; - } - -private: - Matrix & m_matrix; - Vector & m_vector; - Vector m_solution; -}; - - -} } // end namespaces - - -#endif /*_NL_LINEAR_SYSTEM_H_*/ +#ifndef _NL_LINEAR_SYSTEM_H_ +#define _NL_LINEAR_SYSTEM_H_ + + +#include + +#include + +#include "numeric/matrix.h" +#include "numeric/vector.h" + + +namespace Geom { namespace NL { + + +class LinearSystem +{ +public: + 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(); + } + + Matrix & matrix() + { + return m_matrix; + } + + Vector & vector() + { + return m_vector; + } + + const Vector & solution() const + { + return m_solution; + } + +private: + Matrix & m_matrix; + Vector & m_vector; + Vector m_solution; +}; + + +} } // end namespaces + + +#endif /*_NL_LINEAR_SYSTEM_H_*/ diff --git a/src/2geom/numeric/matrix.h b/src/2geom/numeric/matrix.h index bdfab75ef..bdd647e53 100644 --- a/src/2geom/numeric/matrix.h +++ b/src/2geom/numeric/matrix.h @@ -1,207 +1,207 @@ -#ifndef _NL_MATRIX_H_ -#define _NL_MATRIX_H_ - - -#include -#include -#include - -#include - - - -namespace Geom { namespace NL { - -class Matrix; -void swap(Matrix & m1, Matrix & m2); - - -class Matrix -{ -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 ) - : m_rows(_matrix.rows()), m_columns(_matrix.columns()) - { - m_matrix = gsl_matrix_alloc(rows(), columns()); - gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); - } - -// explicit -// Matrix( gsl_matrix* m, size_t n1, size_t n2) -// : m_rows(n1), m_columns(n2), m_matrix(m) -// { -// } - - Matrix & operator=(Matrix const& _matrix) - { - assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); - gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); - - return *this; - } - - virtual ~Matrix() - { - gsl_matrix_free(m_matrix); - } - - void set_all( double x = 0) - { - gsl_matrix_set_all(m_matrix, x); - } - - void set_identity() - { - gsl_matrix_set_identity(m_matrix); - } - - const double & operator() (size_t i, size_t j) const - { - return *gsl_matrix_const_ptr(m_matrix, i, j); - } - - double & operator() (size_t i, size_t j) - { - return *gsl_matrix_ptr(m_matrix, i, j); - } - - gsl_matrix* get_gsl_matrix() - { - return m_matrix; - } - - 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); - } - - Matrix & transpose() - { - gsl_matrix_transpose(m_matrix); - return (*this); - } - - Matrix & scale(double x) - { - gsl_matrix_scale(m_matrix, x); - return (*this); - } - - Matrix & translate(double x) - { - gsl_matrix_add_constant(m_matrix, x); - return (*this); - } - - Matrix & operator+=(Matrix const& _matrix) - { - gsl_matrix_add(m_matrix, _matrix.m_matrix); - return (*this); - } - - Matrix & operator-=(Matrix const& _matrix) - { - gsl_matrix_sub(m_matrix, _matrix.m_matrix); - return (*this); - } - - 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; - } - - friend - void swap(Matrix & m1, Matrix & m2); - - size_t rows() const - { - return m_rows; - } - - size_t columns() const - { - return m_columns; - } - -private: - size_t m_rows, m_columns; - gsl_matrix* m_matrix; -}; - -void swap(Matrix & m1, Matrix & m2) -{ - assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() ); - std::swap(m1.m_matrix, m2.m_matrix); -} - - -} } // end namespaces - -#endif /*_NL_MATRIX_H_*/ +#ifndef _NL_MATRIX_H_ +#define _NL_MATRIX_H_ + + +#include +#include +#include + +#include + + + +namespace Geom { namespace NL { + +class Matrix; +void swap(Matrix & m1, Matrix & m2); + + +class Matrix +{ +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 ) + : m_rows(_matrix.rows()), m_columns(_matrix.columns()) + { + m_matrix = gsl_matrix_alloc(rows(), columns()); + gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); + } + +// explicit +// Matrix( gsl_matrix* m, size_t n1, size_t n2) +// : m_rows(n1), m_columns(n2), m_matrix(m) +// { +// } + + Matrix & operator=(Matrix const& _matrix) + { + assert( rows() == _matrix.rows() && columns() == _matrix.columns() ); + gsl_matrix_memcpy(m_matrix, _matrix.m_matrix); + + return *this; + } + + virtual ~Matrix() + { + gsl_matrix_free(m_matrix); + } + + void set_all( double x = 0) + { + gsl_matrix_set_all(m_matrix, x); + } + + void set_identity() + { + gsl_matrix_set_identity(m_matrix); + } + + const double & operator() (size_t i, size_t j) const + { + return *gsl_matrix_const_ptr(m_matrix, i, j); + } + + double & operator() (size_t i, size_t j) + { + return *gsl_matrix_ptr(m_matrix, i, j); + } + + gsl_matrix* get_gsl_matrix() + { + return m_matrix; + } + + 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); + } + + Matrix & transpose() + { + gsl_matrix_transpose(m_matrix); + return (*this); + } + + Matrix & scale(double x) + { + gsl_matrix_scale(m_matrix, x); + return (*this); + } + + Matrix & translate(double x) + { + gsl_matrix_add_constant(m_matrix, x); + return (*this); + } + + Matrix & operator+=(Matrix const& _matrix) + { + gsl_matrix_add(m_matrix, _matrix.m_matrix); + return (*this); + } + + Matrix & operator-=(Matrix const& _matrix) + { + gsl_matrix_sub(m_matrix, _matrix.m_matrix); + return (*this); + } + + 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; + } + + friend + void swap(Matrix & m1, Matrix & m2); + + size_t rows() const + { + return m_rows; + } + + size_t columns() const + { + return m_columns; + } + +private: + size_t m_rows, m_columns; + gsl_matrix* m_matrix; +}; + +void swap(Matrix & m1, Matrix & m2) +{ + assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() ); + std::swap(m1.m_matrix, m2.m_matrix); +} + + +} } // end namespaces + +#endif /*_NL_MATRIX_H_*/ diff --git a/src/2geom/numeric/vector.h b/src/2geom/numeric/vector.h index 136b2cc3f..b25861e22 100644 --- a/src/2geom/numeric/vector.h +++ b/src/2geom/numeric/vector.h @@ -1,184 +1,185 @@ -#ifndef _NL_VECTOR_H_ -#define _NL_VECTOR_H_ - -#include -#include - -#include - - -namespace Geom { namespace NL { - -class Vector; -void swap(Vector & v1, Vector & v2); - - -class Vector -{ -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) - : m_size(_vector.size()) - { - m_vector = gsl_vector_alloc(size()); - gsl_vector_memcpy(m_vector, _vector.m_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); - } - - void set_all(double x = 0) - { - gsl_vector_set_all(m_vector, x); - } - - void set_basis(size_t i) - { - gsl_vector_set_basis(m_vector, i); - } - - double const& operator[](size_t i) const - { - return *gsl_vector_const_ptr(m_vector, i); - } - - double & operator[](size_t i) - { - return *gsl_vector_ptr(m_vector, i); - } - - 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); - } - - Vector & scale(double x) - { - gsl_vector_scale(m_vector, x); - return (*this); - } - - Vector & translate(double x) - { - gsl_vector_add_constant(m_vector, x); - return (*this); - } - - Vector & operator+=(Vector const& _vector) - { - gsl_vector_add(m_vector, _vector.m_vector); - return (*this); - } - - Vector & operator-=(Vector const& _vector) - { - gsl_vector_sub(m_vector, _vector.m_vector); - return (*this); - } - - 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); - } - - friend - void swap(Vector & v1, Vector & v2); - - size_t size() const - { - return m_size; - } - -private: - size_t m_size; - gsl_vector* m_vector; -}; - -void swap(Vector & v1, Vector & v2) -{ - assert( v1.size() == v2.size() ); - std::swap(v1.m_vector, v2.m_vector); -} - -} } // end namespaces - - -#endif /*_NL_VECTOR_H_*/ +#ifndef _NL_VECTOR_H_ +#define _NL_VECTOR_H_ + +#include +#include + +#include + + +namespace Geom { namespace NL { + +class Vector; +void swap(Vector & v1, Vector & v2); + + +class Vector +{ +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) + : m_size(_vector.size()) + { + m_vector = gsl_vector_alloc(size()); + gsl_vector_memcpy(m_vector, _vector.m_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); + } + + void set_all(double x = 0) + { + gsl_vector_set_all(m_vector, x); + } + + void set_basis(size_t i) + { + gsl_vector_set_basis(m_vector, i); + } + + double const& operator[](size_t i) const + { + return *gsl_vector_const_ptr(m_vector, i); + } + + double & operator[](size_t i) + { + return *gsl_vector_ptr(m_vector, i); + } + + 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); + } + + Vector & scale(double x) + { + gsl_vector_scale(m_vector, x); + return (*this); + } + + Vector & translate(double x) + { + gsl_vector_add_constant(m_vector, x); + return (*this); + } + + Vector & operator+=(Vector const& _vector) + { + gsl_vector_add(m_vector, _vector.m_vector); + return (*this); + } + + Vector & operator-=(Vector const& _vector) + { + gsl_vector_sub(m_vector, _vector.m_vector); + return (*this); + } + + 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); + } + + friend + void swap(Vector & v1, Vector & v2); + + size_t size() const + { + return m_size; + } + +private: + size_t m_size; + gsl_vector* m_vector; +}; + +void swap(Vector & v1, Vector & v2) +{ + assert( v1.size() == v2.size() ); + std::swap(v1.m_vector, v2.m_vector); +} + +} } // end namespaces + + +#endif /*_NL_VECTOR_H_*/ diff --git a/src/2geom/path.h b/src/2geom/path.h index 0493b4a59..abed8604a 100644 --- a/src/2geom/path.h +++ b/src/2geom/path.h @@ -208,7 +208,7 @@ public: bool isDegenerate() const { return inner.isConstant(); } void setInitial(Point v) { setPoint(0, v); } - void setFinal(Point v) { setPoint(1, v); } + void setFinal(Point v) { setPoint(order, v); } void setPoint(unsigned ix, Point v) { inner[X].setPoint(ix, v[X]); inner[Y].setPoint(ix, v[Y]); } Point const operator[](unsigned ix) const { return Point(inner[X][ix], inner[Y][ix]); } @@ -308,10 +308,10 @@ double LineSegment::nearestPoint(Point const& p, double from, double to) const; -class SVGEllipticalArc : public Curve +class EllipticalArc : public Curve { public: - SVGEllipticalArc() + EllipticalArc() : 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) @@ -320,10 +320,10 @@ class SVGEllipticalArc : public Curve m_center = Point(0,0); } - SVGEllipticalArc( Point _initial_point, double _rx, double _ry, - double _rot_angle, bool _large_arc, bool _sweep, - Point _final_point - ) + EllipticalArc( 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) @@ -348,7 +348,7 @@ class SVGEllipticalArc : public Curve Curve* duplicate() const { - return new SVGEllipticalArc(*this); + return new EllipticalArc(*this); } double center(unsigned int i) const @@ -495,13 +495,13 @@ class SVGEllipticalArc : public Curve return pointAtAngle(tt); } - std::pair + std::pair subdivide(Coord t) const { - SVGEllipticalArc* arc1 = static_cast(portion(0, t)); - SVGEllipticalArc* arc2 = static_cast(portion(t, 1)); + EllipticalArc* arc1 = static_cast(portion(0, t)); + EllipticalArc* arc2 = static_cast(portion(t, 1)); assert( arc1 != NULL && arc2 != NULL); - std::pair arc_pair(*arc1, *arc2); + std::pair arc_pair(*arc1, *arc2); delete arc1; delete arc2; return arc_pair; @@ -512,7 +512,7 @@ class SVGEllipticalArc : public Curve // the arc is the same but traversed in the opposite direction Curve* reverse() const { - SVGEllipticalArc* rarc = new SVGEllipticalArc( *this ); + EllipticalArc* rarc = new EllipticalArc( *this ); rarc->m_sweep = !m_sweep; rarc->m_initial_point = m_final_point; rarc->m_final_point = m_initial_point; @@ -541,7 +541,7 @@ class SVGEllipticalArc : public Curve double m_start_angle, m_end_angle; Point m_center; -}; // end class SVGEllipticalArc +}; // end class EllipticalArc @@ -577,6 +577,17 @@ public: return old; } + BaseIterator &operator--() { + --impl_; + return *this; + } + + BaseIterator operator--(int) { + BaseIterator old=*this; + --(*this); + return old; + } + private: BaseIterator(IteratorImpl const &pos) : impl_(pos) {} @@ -668,14 +679,13 @@ public: Curve const &operator[](unsigned i) const { return *curves_[i]; } - iterator begin() { return curves_.begin(); } - iterator end() { return curves_.end()-1; } - Curve const &front() const { return *curves_[0]; } Curve const &back() const { return *curves_[curves_.size()-2]; } const_iterator begin() const { return curves_.begin(); } const_iterator end() const { return curves_.end()-1; } + iterator begin() { return curves_.begin(); } + iterator end() { return curves_.end()-1; } const_iterator end_open() const { return curves_.end()-1; } const_iterator end_closed() const { return curves_.end(); } @@ -910,6 +920,42 @@ public: Point initialPoint() const { return (*final_)[1]; } Point finalPoint() const { return (*final_)[0]; } + void setInitial(Point const& p) + { + if ( empty() ) return; + Curve* head = front().duplicate(); + head->setInitial(p); + Sequence::iterator replaced = curves_.begin(); + Sequence source(1, head); + try + { + do_update(replaced, replaced + 1, source.begin(), source.end()); + } + catch (...) + { + delete_range(source.begin(), source.end()); + throw; + } + } + + void setFinal(Point const& p) + { + if ( empty() ) return; + Curve* tail = back().duplicate(); + tail->setFinal(p); + Sequence::iterator replaced = curves_.end() - 2; + Sequence source(1, tail); + try + { + do_update(replaced, replaced + 1, source.begin(), source.end()); + } + catch (...) + { + delete_range(source.begin(), source.end()); + throw; + } + } + void append(Curve const &curve); void append(D2 const &curve); void append(Path const &other); diff --git a/src/2geom/poly.cpp b/src/2geom/poly.cpp index 17f286494..4f5ba6c55 100644 --- a/src/2geom/poly.cpp +++ b/src/2geom/poly.cpp @@ -11,6 +11,7 @@ Poly Poly::operator*(const Poly& p) const { } return result; } +#define HAVE_GSL #ifdef HAVE_GSL #include #endif diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp index c60132043..920fd37fe 100644 --- a/src/2geom/sbasis.cpp +++ b/src/2geom/sbasis.cpp @@ -69,22 +69,15 @@ bool SBasis::isFinite() const { } std::vector SBasis::valueAndDerivatives(double t, unsigned n) const { - std::vector ret; + std::vector ret(n); if(n==1) { ret.push_back(valueAt(t)); return ret; } - if(n==2) { - double der; - ret.push_back(valueAndDerivative(t, der)); - ret.push_back(der); - return ret; - } SBasis tmp = *this; - while(n > 0) { - ret.push_back(tmp.valueAt(t)); - tmp = derivative(tmp); - n--; + for(unsigned i = 0; i < n; i++) { + ret[i] = tmp.valueAt(t); + tmp.derive(); } return ret; } @@ -191,6 +184,7 @@ SBasis shift(Linear const &a, int sh) { return c; } +#if 0 SBasis multiply(SBasis const &a, SBasis const &b) { // c = {a0*b0 - shift(1, a.Tri*b.Tri), a1*b1 - shift(1, a.Tri*b.Tri)} @@ -199,7 +193,6 @@ SBasis multiply(SBasis const &a, SBasis const &b) { if(a.isZero() || b.isZero()) return c; c.resize(a.size() + b.size(), Linear(0,0)); - c[0] = Linear(0,0); for(unsigned j = 0; j < b.size(); j++) { for(unsigned i = j; i < a.size()+j; i++) { double tri = Tri(b[j])*Tri(a[i-j]); @@ -216,7 +209,36 @@ SBasis multiply(SBasis const &a, SBasis const &b) { //assert(!(0 == c.back()[0] && 0 == c.back()[1])); return c; } +#else +SBasis multiply_add(SBasis const &a, SBasis const &b, SBasis c) { + if(a.isZero() || b.isZero()) + return c; + c.resize(a.size() + b.size(), Linear(0,0)); + for(unsigned j = 0; j < b.size(); j++) { + for(unsigned i = j; i < a.size()+j; i++) { + double tri = Tri(b[j])*Tri(a[i-j]); + c[i+1/*shift*/] += Linear(Hat(-tri)); + } + } + for(unsigned j = 0; j < b.size(); j++) { + for(unsigned i = j; i < a.size()+j; i++) { + for(unsigned dim = 0; dim < 2; dim++) + c[i][dim] += b[j][dim]*a[i-j][dim]; + } + } + c.normalize(); + //assert(!(0 == c.back()[0] && 0 == c.back()[1])); + return c; +} + +SBasis multiply(SBasis const &a, SBasis const &b) { + SBasis c; + if(a.isZero() || b.isZero()) + return c; + return multiply_add(a, b, c); +} +#endif SBasis integral(SBasis const &c) { SBasis a; a.resize(c.size() + 1, Linear(0,0)); @@ -240,23 +262,41 @@ SBasis derivative(SBasis const &a) { SBasis c; c.resize(a.size(), Linear(0,0)); - for(unsigned k = 0; k < a.size(); k++) { - double d = (2*k+1)*Tri(a[k]); - - for(unsigned dim = 0; dim < 2; dim++) { - c[k][dim] = d; - if(k+1 < a.size()) { - if(dim) - c[k][dim] = d - (k+1)*a[k+1][dim]; - else - c[k][dim] = d + (k+1)*a[k+1][dim]; - } - } + for(unsigned k = 0; k < a.size()-1; k++) { + double d = (2*k+1)*(a[k][1] - a[k][0]); + + c[k][0] = d + (k+1)*a[k+1][0]; + c[k][1] = d - (k+1)*a[k+1][1]; + } + int k = a.size()-1; + double d = (2*k+1)*(a[k][1] - a[k][0]); + if(d == 0) + c.pop_back(); + else { + c[k][0] = d; + c[k][1] = d; } return c; } +void SBasis::derive() { // in place version + for(unsigned k = 0; k < size()-1; k++) { + double d = (2*k+1)*((*this)[k][1] - (*this)[k][0]); + + (*this)[k][0] = d + (k+1)*(*this)[k+1][0]; + (*this)[k][1] = d - (k+1)*(*this)[k+1][1]; + } + int k = size()-1; + double d = (2*k+1)*((*this)[k][1] - (*this)[k][0]); + if(d == 0) + pop_back(); + else { + (*this)[k][0] = d; + (*this)[k][1] = d; + } +} + //TODO: convert int k to unsigned k, and remove cast SBasis sqrt(SBasis const &a, int k) { SBasis c; @@ -321,7 +361,7 @@ SBasis compose(SBasis const &a, SBasis const &b) { SBasis r; for(int i = a.size()-1; i >= 0; i--) { - r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s); + r = multiply_add(r, s, SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1]); } return r; } @@ -333,7 +373,7 @@ SBasis compose(SBasis const &a, SBasis const &b, unsigned k) { SBasis r; for(int i = a.size()-1; i >= 0; i--) { - r = SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1] + multiply(r,s); + r = multiply_add(r, s, SBasis(Linear(Hat(a[i][0]))) - b*a[i][0] + b*a[i][1]); } r.truncate(k); return r; diff --git a/src/2geom/sbasis.h b/src/2geom/sbasis.h index 7f0d88f8c..0a48ff1cf 100644 --- a/src/2geom/sbasis.h +++ b/src/2geom/sbasis.h @@ -57,6 +57,9 @@ public: SBasis(Linear const & bo) { push_back(bo); } + SBasis(Linear* bo) { + push_back(*bo); + } //IMPL: FragmentConcept typedef double output_type; @@ -86,29 +89,15 @@ public: double valueAt(double t) const { double s = t*(1-t); double p0 = 0, p1 = 0; - double sk = 1; -//TODO: rewrite as horner - for(unsigned k = 0; k < size(); k++) { - p0 += sk*(*this)[k][0]; - p1 += sk*(*this)[k][1]; - sk *= s; + for(unsigned k = size(); k > 0; k--) { + const Linear &lin = (*this)[k-1]; + p0 = p0*s + lin[0]; + p1 = p1*s + lin[1]; } return (1-t)*p0 + t*p1; } - double valueAndDerivative(double t, double &der) const { - double s = t*(1-t); - double p0 = 0, p1 = 0; - double sk = 1; -//TODO: rewrite as horner - for(unsigned k = 0; k < size(); k++) { - p0 += sk*(*this)[k][0]; - p1 += sk*(*this)[k][1]; - sk *= s; - } - // p0 and p1 at this point form a linear approximation at t - der = p1 - p0; - return (1-t)*p0 + t*p1; - } + //double valueAndDerivative(double t, double &der) const { + //} double operator()(double t) const { return valueAt(t); } @@ -135,7 +124,10 @@ public: while(!empty() && 0 == back()[0] && 0 == back()[1]) pop_back(); } + void truncate(unsigned k) { if(k < size()) resize(k); } +private: + void derive(); // in place version }; //TODO: figure out how to stick this in linear, while not adding an sbasis dep @@ -244,6 +236,8 @@ inline SBasis truncate(SBasis const &a, unsigned terms) { } SBasis multiply(SBasis const &a, SBasis const &b); +// This performs a multiply and accumulate operation in about the same time as multiply. return a*b + c +SBasis multiply_add(SBasis const &a, SBasis const &b, SBasis c); SBasis integral(SBasis const &c); SBasis derivative(SBasis const &a); diff --git a/src/2geom/solve-bezier-one-d.cpp b/src/2geom/solve-bezier-one-d.cpp index 1338faa7c..cd7c36cc3 100644 --- a/src/2geom/solve-bezier-one-d.cpp +++ b/src/2geom/solve-bezier-one-d.cpp @@ -12,24 +12,37 @@ namespace Geom{ template static int SGN(t x) { return (x > 0 ? 1 : (x < 0 ? -1 : 0)); } -/* - * Forward declarations - */ -static void -Bernstein(double const *V, - unsigned degree, - double t, - double *Left, - double *Right); - -static unsigned -control_poly_flat_enough(double const *V, unsigned degree, - double left_t, double right_t); - -const unsigned MAXDEPTH = 64; /* Maximum depth for recursion */ +const unsigned MAXDEPTH = 23; // Maximum depth for recursion. Using floats means 23 bits precision max const double BEPSILON = ldexp(1.0,-MAXDEPTH-1); /*Flatness control value */ +/** + * This function is called _a lot_. We have included various manual memory management stuff to reduce the amount of mallocing that goes on. In the future it is possible that this will hurt performance. + **/ +class Bernsteins{ +public: + double *Vtemp; + unsigned N,degree; + std::vector &solutions; + Bernsteins(int degr, std::vector &so) : N(degr+1), degree(degr),solutions(so) { + Vtemp = new double[N*2]; + } + ~Bernsteins() { + delete[] Vtemp; + } + void subdivide(double const *V, + double t, + double *Left, + double *Right); + + unsigned + control_poly_flat_enough(double const *V); + + void + find_bernstein_roots(double const *w, /* The control points */ + unsigned depth, /* The depth of the recursion */ + double left_t, double right_t); +}; /* * find_bernstein_roots : Given an equation in Bernstein-Bernstein form, find all * of the roots in the open interval (0, 1). Return the number of roots found. @@ -41,10 +54,20 @@ find_bernstein_roots(double const *w, /* The control points */ unsigned depth, /* The depth of the recursion */ double left_t, double right_t) { + Bernsteins B(degree, solutions); + B.find_bernstein_roots(w, depth, left_t, right_t); +} + +void +Bernsteins::find_bernstein_roots(double const *w, /* The control points */ + unsigned depth, /* The depth of the recursion */ + double left_t, double right_t) +{ + unsigned n_crossings = 0; /* Number of zero-crossings */ int old_sign = SGN(w[0]); - for (unsigned i = 1; i <= degree; i++) { + for (unsigned i = 1; i < N; i++) { int sign = SGN(w[i]); if (sign) { if (sign != old_sign && old_sign) { @@ -54,46 +77,71 @@ find_bernstein_roots(double const *w, /* The control points */ } } - switch (n_crossings) { - case 0: /* No solutions here */ + if (n_crossings == 0) // no solutions here return; - case 1: + if (n_crossings == 1) { /* Unique solution */ /* Stop recursion when the tree is deep enough */ /* if deep enough, return 1 solution at midpoint */ if (depth >= MAXDEPTH) { + //printf("bottom out %d\n", depth); solutions.push_back((left_t + right_t) / 2.0); return; } // I thought secant method would be faster here, but it'aint. -- njh + + // The original code went to some effort to try and terminate early when the curve is sufficiently flat. However, it seems that in practice it almost always bottoms out, so leaving this code out actually speeds things up + if(0) if (control_poly_flat_enough(w)) { + //printf("flatten out %d\n", depth); + const double Ax = right_t - left_t; + const double Ay = w[degree] - w[0]; + + solutions.push_back(left_t - Ax*w[0] / Ay); + return; + } + } - if (control_poly_flat_enough(w, degree, left_t, right_t)) { - const double Ax = right_t - left_t; - const double Ay = w[degree] - w[0]; + /* Otherwise, solve recursively after subdividing control polygon */ + double Left[N], /* New left and right */ + Right[N]; /* control polygons */ + const double t = 0.5; - solutions.push_back(left_t - Ax*w[0] / Ay); - return; + +/* + * Bernstein : + * Evaluate a Bernstein function at a particular parameter value + * Fill in control points for resulting sub-curves. + * + */ + for (unsigned i = 0; i < N; i++) + Vtemp[i] = w[i]; + + /* Triangle computation */ + const double omt = (1-t); + Left[0] = Vtemp[0]; + Right[degree] = Vtemp[degree]; + double *prev_row = Vtemp; + double *row = Vtemp + N; + for (unsigned i = 1; i < N; i++) { + for (unsigned j = 0; j < N - i; j++) { + row[j] = omt*prev_row[j] + t*prev_row[j+1]; } - break; + Left[i] = row[0]; + Right[degree-i] = row[degree-i]; + std::swap(prev_row, row); } - - /* Otherwise, solve recursively after subdividing control polygon */ - double Left[degree+1], /* New left and right */ - Right[degree+1]; /* control polygons */ - const double split = 0.5; - Bernstein(w, degree, split, Left, Right); - double mid_t = left_t*(1-split) + right_t*split; + double mid_t = left_t*(1-t) + right_t*t; - find_bernstein_roots(Left, degree, solutions, depth+1, left_t, mid_t); + find_bernstein_roots(Left, depth+1, left_t, mid_t); /* Solution is exactly on the subdivision point. */ if (Right[0] == 0) solutions.push_back(mid_t); - find_bernstein_roots(Right, degree, solutions, depth+1, mid_t, right_t); + find_bernstein_roots(Right, depth+1, mid_t, right_t); } /* @@ -102,10 +150,8 @@ find_bernstein_roots(double const *w, /* The control points */ * for recursive subdivision to bottom out. * */ -static unsigned -control_poly_flat_enough(double const *V, /* Control points */ - unsigned degree, - double left_t, double right_t) /* Degree of polynomial */ +unsigned +Bernsteins::control_poly_flat_enough(double const *V) { /* Find the perpendicular distance from each interior control point to line connecting V[0] and * V[degree] */ @@ -113,8 +159,6 @@ control_poly_flat_enough(double const *V, /* Control points */ /* Derive the implicit equation for line connecting first */ /* and last control points */ const double a = V[0] - V[degree]; - const double b = right_t - left_t; - const double c = left_t * V[degree] - right_t * V[0] + a * left_t; double max_distance_above = 0.0; double max_distance_below = 0.0; @@ -122,7 +166,7 @@ control_poly_flat_enough(double const *V, /* Control points */ for (unsigned i = 1; i < degree; i++) { ii += dii; /* Compute distance from each of the points to that line */ - const double d = (a + V[i]) * ii*b + c; + const double d = (a + V[i]) * ii - a; double dist = d*d; // Find the largest distance if (d < 0.0) @@ -131,56 +175,22 @@ control_poly_flat_enough(double const *V, /* Control points */ max_distance_above = std::max(max_distance_above, dist); } - const double abSquared = (a * a) + (b * b); + const double abSquared = 1./((a * a) + 1); - const double intercept_1 = -(c + max_distance_above / abSquared); - const double intercept_2 = -(c + max_distance_below / abSquared); + const double intercept_1 = (a - max_distance_above * abSquared); + const double intercept_2 = (a - max_distance_below * abSquared); /* Compute bounding interval*/ const double left_intercept = std::min(intercept_1, intercept_2); const double right_intercept = std::max(intercept_1, intercept_2); const double error = 0.5 * (right_intercept - left_intercept); - - if (error < BEPSILON * a) - return 1; - - return 0; + //printf("error %g %g %g\n", error, a, BEPSILON * a); + return error < BEPSILON * a; } -/* - * Bernstein : - * Evaluate a Bernstein function at a particular parameter value - * Fill in control points for resulting sub-curves. - * - */ -static void -Bernstein(double const *V, /* Control pts */ - unsigned degree, /* Degree of bernstein curve */ - double t, /* Parameter value */ - double *Left, /* RETURN left half ctl pts */ - double *Right) /* RETURN right half ctl pts */ -{ - double Vtemp[degree+1][degree+1]; - - /* Copy control points */ - std::copy(V, V+degree+1, Vtemp[0]); - - /* Triangle computation */ - const double omt = (1-t); - Left[0] = Vtemp[0][0]; - Right[degree] = Vtemp[0][degree]; - for (unsigned i = 1; i <= degree; i++) { - for (unsigned j = 0; j <= degree - i; j++) { - Vtemp[i][j] = omt*Vtemp[i-1][j] + t*Vtemp[i-1][j+1]; - } - Left[i] = Vtemp[i][0]; - Right[degree-i] = Vtemp[i][degree-i]; - } -} - }; /* diff --git a/src/2geom/svg-path.cpp b/src/2geom/svg-path.cpp index c09e5c929..2c275cb4c 100644 --- a/src/2geom/svg-path.cpp +++ b/src/2geom/svg-path.cpp @@ -51,7 +51,7 @@ void output(QuadraticBezier const &curve, SVGPathSink &sink) { sink.quadTo(curve[1], curve[2]); } -void output(SVGEllipticalArc const &/*curve*/, SVGPathSink &/*sink*/) { +void output(EllipticalArc const &/*curve*/, SVGPathSink &/*sink*/) { // FIXME THROW_NOTIMPLEMENTED(); } @@ -75,7 +75,7 @@ void output_svg_path(Path &path, SVGPathSink &sink) { output_as(*iter, sink) || output_as(*iter, sink) || output_as(*iter, sink) || - output_as(*iter, sink) || + output_as(*iter, sink) || output_as(*iter, sink); } diff --git a/src/2geom/svg-path.h b/src/2geom/svg-path.h index 7dd2f31ac..870285b7f 100644 --- a/src/2geom/svg-path.h +++ b/src/2geom/svg-path.h @@ -78,7 +78,7 @@ public: void arcTo(double rx, double ry, double angle, bool large_arc, bool sweep, Point p) { - _path.template appendNew(rx, ry, angle, + _path.template appendNew(rx, ry, angle, large_arc, sweep, p); } -- 2.30.2