![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
summary | shortlog | log | commit | commitdiff | tree
raw | patch | inline | side by side (parent: b07935a)
raw | patch | inline | side by side (parent: b07935a)
| author | mcecchetti <mcecchetti@users.sourceforge.net> | |
| Tue, 20 May 2008 22:29:23 +0000 (22:29 +0000) | ||
| committer | mcecchetti <mcecchetti@users.sourceforge.net> | |
| Tue, 20 May 2008 22:29:23 +0000 (22:29 +0000) |
19 files changed:
| src/2geom/CMakeLists.txt | patch | blob | history | |
| src/2geom/Makefile_insert | patch | blob | history | |
| src/2geom/choose.h | patch | blob | history | |
| src/2geom/elliptical-arc.cpp | [new file with mode: 0644] | patch | blob |
| src/2geom/exception.h | patch | blob | history | |
| src/2geom/forward.h | patch | blob | history | |
| src/2geom/nearest-point.cpp | patch | blob | history | |
| src/2geom/nearest-point.h | patch | blob | history | |
| src/2geom/numeric/linear_system.h | patch | blob | history | |
| src/2geom/numeric/matrix.h | patch | blob | history | |
| src/2geom/numeric/vector.h | patch | blob | history | |
| src/2geom/path.h | patch | blob | history | |
| src/2geom/poly.cpp | patch | blob | history | |
| src/2geom/sbasis.cpp | patch | blob | history | |
| src/2geom/sbasis.h | patch | blob | history | |
| src/2geom/solve-bezier-one-d.cpp | patch | blob | history | |
| src/2geom/svg-elliptical-arc.cpp | [deleted file] | patch | blob | history |
| src/2geom/svg-path.cpp | patch | blob | history | |
| src/2geom/svg-path.h | patch | blob | history |
index 309d4cdd99e0a72b7456d0be12a5da80a6ab1edc..ac019419ddd93c152d4dc6a09ccd64e9fe261677 100644 (file)
--- a/src/2geom/CMakeLists.txt
+++ b/src/2geom/CMakeLists.txt
d2.h
d2-sbasis.cpp
d2-sbasis.h
+elliptical-arc.cpp
exception.h
forward.h
geom.cpp
solve-bezier-parametric.cpp
solver.h
sturm.h
-svg-elliptical-arc.cpp
svg-path.cpp
svg-path.h
svg-path-parser.cpp
index 6f2a3e642cbfcf16687dd797fda297575a80ee58..e2c0d6aa37fd4a54bb4d6dbfab40bcb2a4977137 100644 (file)
2geom/d2-sbasis.h \\r
2geom/d2.h \\r
2geom/exception.h \\r
+ 2geom/elliptical-arc.cpp \\r
2geom/forward.h \\r
2geom/geom.cpp \\r
2geom/geom.h \\r
2geom/solve-bezier-one-d.cpp \\r
2geom/solve-bezier-parametric.cpp \\r
2geom/solver.h \\r
- 2geom/sturm.h \\r
- 2geom/svg-elliptical-arc.cpp \
+ 2geom/sturm.h \
2geom/svg-path.cpp \\r
2geom/svg-path.h \\r
2geom/svg-path-parser.cpp \\r
diff --git a/src/2geom/choose.h b/src/2geom/choose.h
index 46c5b2fb47cc81855107e788d353f2ba45a14d05..f3a8b0f4fcda87fd796c06080bd179cf62be3031 100644 (file)
--- a/src/2geom/choose.h
+++ b/src/2geom/choose.h
return pascals_triangle[row+k];
}
+// Is it faster to store them or compute them on demand?
+/*template <typename T>
+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/elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp
--- /dev/null
@@ -0,0 +1,919 @@
+/*
+ * SVG Elliptical Arc Class
+ *
+ * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include "path.h"
+#include "angle.h"
+
+#include <gsl/gsl_poly.h>
+
+#include <cfloat>
+
+
+
+
+namespace Geom
+{
+
+
+Rect EllipticalArc::boundsExact() const
+{
+ std::vector<double> extremes(4);
+ double cosrot = std::cos(rotation_angle());
+ double sinrot = std::sin(rotation_angle());
+ extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );
+ extremes[1] = extremes[0] + M_PI;
+ if ( extremes[0] < 0 ) extremes[0] += 2*M_PI;
+ extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );
+ extremes[3] = extremes[2] + M_PI;
+ if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;
+
+
+ std::vector<double>arc_extremes(4);
+ arc_extremes[0] = initialPoint()[X];
+ arc_extremes[1] = finalPoint()[X];
+ if ( arc_extremes[0] < arc_extremes[1] )
+ std::swap(arc_extremes[0], arc_extremes[1]);
+ arc_extremes[2] = initialPoint()[Y];
+ arc_extremes[3] = finalPoint()[Y];
+ if ( arc_extremes[2] < arc_extremes[3] )
+ std::swap(arc_extremes[2], arc_extremes[3]);
+
+
+ if ( start_angle() < end_angle() )
+ {
+ if ( sweep_flag() )
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() < extremes[i] && extremes[i] < end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ else
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() > extremes[i] || extremes[i] > end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ }
+ else
+ {
+ if ( sweep_flag() )
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() < extremes[i] || extremes[i] < end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ else
+ {
+ for ( unsigned int i = 0; i < extremes.size(); ++i )
+ {
+ if ( start_angle() > extremes[i] && extremes[i] > end_angle() )
+ {
+ arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];
+ }
+ }
+ }
+ }
+
+ return Rect( Point(arc_extremes[1], arc_extremes[3]) ,
+ Point(arc_extremes[0], arc_extremes[2]) );
+
+}
+
+
+std::vector<double>
+EllipticalArc::roots(double v, Dim2 d) const
+{
+ if ( d > Y )
+ {
+ THROW_RANGEERROR("dimention out of range");
+ }
+
+ std::vector<double> 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<double> arc_sol;
+ for (unsigned int i = 0; i < sol.size(); ++i )
+ {
+ //std::cerr << "s = " << rad_to_deg(sol[i]);
+ sol[i] = map_to_01(sol[i]);
+ //std::cerr << " -> t: " << sol[i] << std::endl;
+ if ( !(sol[i] < 0 || sol[i] > 1) )
+ arc_sol.push_back(sol[i]);
+ }
+ return arc_sol;
+
+
+// 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<Point>
+EllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const
+{
+ std::vector<Point> 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<SBasis> 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<SBasis> 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<double> 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<double> 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<double> solX = roots(np[X],X);
+// std::vector<double> solY = roots(np[Y],Y);
+// double t;
+// if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))
+// {
+// t = solX[0];
+// }
+// else
+// {
+// t = solX[1];
+// }
+// if ( !(t < from || t > to) )
+// {
+// result.push_back(t);
+// }
+// else
+// {
+//
+// }
+ }
+
+ // solve the equation <D(E(t),t)|E(t)-p> == 0
+ // that provides min and max distance points
+ // on the ellipse E wrt the point p
+ // after the substitutions:
+ // cos(t) = (1 - s^2) / (1 + s^2)
+ // sin(t) = 2t / (1 + s^2)
+ // where s = tan(t/2)
+ // we get a 4th degree equation in s
+ /*
+ * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) +
+ * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) +
+ * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) +
+ * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])
+ */
+
+ Point p_c = p - center();
+ double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));
+ double cosrot = std::cos( rotation_angle() );
+ double sinrot = std::sin( rotation_angle() );
+ double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);
+ 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<double> real_sol;
+ // gsl_poly_complex_solve raises an error
+ // if the leading coefficient is zero
+ if ( are_near(coeff[4], 0) )
+ {
+ real_sol.push_back(0);
+ if ( !are_near(coeff[3], 0) )
+ {
+ double sq = -coeff[1] / coeff[3];
+ if ( sq > 0 )
+ {
+ double s = std::sqrt(sq);
+ real_sol.push_back(s);
+ real_sol.push_back(-s);
+ }
+ }
+ }
+ else
+ {
+ 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 <D(E)|E-p> -> 0 iff coeff[4] == 0
+ // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity
+ if ( (real_sol.size() % 2) != 0 )
+ {
+ real_sol.push_back(M_PI);
+ }
+
+ double mindistsq1 = std::numeric_limits<double>::max();
+ double mindistsq2 = std::numeric_limits<double>::max();
+ double dsq;
+ unsigned int mi1, mi2;
+ for ( unsigned int i = 0; i < real_sol.size(); ++i )
+ {
+ dsq = distanceSq(p, pointAtAngle(real_sol[i]));
+ if ( mindistsq1 > dsq )
+ {
+ mindistsq2 = mindistsq1;
+ mi2 = mi1;
+ mindistsq1 = dsq;
+ mi1 = i;
+ }
+ else if ( mindistsq2 > dsq )
+ {
+ mindistsq2 = dsq;
+ mi2 = i;
+ }
+ }
+
+ double t = map_to_01( real_sol[mi1] );
+ if ( !(t < from || t > to) )
+ {
+ result.push_back(t);
+ }
+
+ bool second_sol = false;
+ t = map_to_01( real_sol[mi2] );
+ if ( real_sol.size() == 4 && !(t < from || t > to) )
+ {
+ if ( result.empty() || are_near(mindistsq1, mindistsq2) )
+ {
+ result.push_back(t);
+ second_sol = true;
+ }
+ }
+
+ // we need to test extreme points too
+ double dsq1 = distanceSq(p, pointAt(from));
+ double dsq2 = distanceSq(p, pointAt(to));
+ if ( second_sol )
+ {
+ if ( mindistsq2 > dsq1 )
+ {
+ result.clear();
+ result.push_back(from);
+ mindistsq2 = dsq1;
+ }
+ else if ( are_near(mindistsq2, dsq) )
+ {
+ result.push_back(from);
+ }
+ if ( mindistsq2 > dsq2 )
+ {
+ result.clear();
+ result.push_back(to);
+ }
+ else if ( are_near(mindistsq2, dsq2) )
+ {
+ result.push_back(to);
+ }
+
+ }
+ else
+ {
+ if ( result.empty() )
+ {
+ if ( are_near(dsq1, dsq2) )
+ {
+ result.push_back(from);
+ result.push_back(to);
+ }
+ else if ( dsq2 > dsq1 )
+ {
+ result.push_back(from);
+ }
+ else
+ {
+ result.push_back(to);
+ }
+ }
+ }
+
+ return result;
+}
+
+
+} // 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 2749ec63a4515645bc39b4bc4a8c0644d4313f58..6151f7a57355be8b7b942c6e25e89e29d71c1ccf 100644 (file)
--- a/src/2geom/exception.h
+++ b/src/2geom/exception.h
};
#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 a34384023c33e8c826a8ef8838645259fa52698f..454f9bf157164342820816deed66a5c15d1bd1ee 100644 (file)
--- a/src/2geom/forward.h
+++ b/src/2geom/forward.h
class Shape;
class Region;
-class SVGEllipticalArc;
+class EllipticalArc;
class SVGPathSink;
template <typename> class SVGPathGenerator;
index 074de1cb3364c54780731ffb0588c23e33eb9be3..de880713f651bcc8feb710c1bc3d9f364c5f4c83 100644 (file)
-/*\r
- * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>\r
- *\r
- * Authors:\r
- * \r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- * \r
- * Copyright 2007-2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#include "nearest-point.h"\r
-\r
-namespace Geom\r
-{\r
-\r
-////////////////////////////////////////////////////////////////////////////////\r
-// D2<SBasis> versions\r
-\r
-/*\r
- * Return the parameter t of a nearest point on the portion of the curve "c", \r
- * related to the interval [from, to], to the point "p".\r
- * The needed curve derivative "dc" is passed as parameter.\r
- * The function return the first nearest point to "p" that is found.\r
- */\r
-\r
-double nearest_point( Point const& p, \r
- D2<SBasis> const& c, \r
- D2<SBasis> const& dc, \r
- double from, double to )\r
-{\r
- if ( from > to ) std::swap(from, to);\r
- if ( from < 0 || to > 1 )\r
- {\r
- THROW_RANGEERROR("[from,to] interval out of bounds");\r
- }\r
-\r
- SBasis dd = dot(c - p, dc); \r
- std::vector<double> zeros = Geom::roots(dd);\r
- \r
- double closest = from;\r
- double min_dist_sq = L2sq(c(from) - p);\r
- double distsq;\r
- for ( unsigned int i = 0; i < zeros.size(); ++i )\r
- {\r
- distsq = L2sq(c(zeros[i]) - p);\r
- if ( min_dist_sq > L2sq(c(zeros[i]) - p) )\r
- {\r
- closest = zeros[i];\r
- min_dist_sq = distsq;\r
- }\r
- }\r
- if ( min_dist_sq > L2sq( c(to) - p ) )\r
- closest = to;\r
- return closest;\r
-\r
-}\r
-\r
-/*\r
- * Return the parameters t of all the nearest points on the portion of \r
- * the curve "c", related to the interval [from, to], to the point "p".\r
- * The needed curve derivative "dc" is passed as parameter.\r
- */\r
-\r
-std::vector<double> \r
-all_nearest_points( Point const& p, \r
- D2<SBasis> const& c, \r
- D2<SBasis> const& dc, \r
- double from, double to )\r
-{\r
- std::swap(from, to);\r
- if ( from > to ) std::swap(from, to);\r
- if ( from < 0 || to > 1 )\r
- {\r
- THROW_RANGEERROR("[from,to] interval out of bounds");\r
- }\r
-\r
- std::vector<double> result;\r
- SBasis dd = dot(c - p, Geom::derivative(c));\r
- \r
- std::vector<double> zeros = Geom::roots(dd);\r
- std::vector<double> candidates;\r
- candidates.push_back(from);\r
- candidates.insert(candidates.end(), zeros.begin(), zeros.end());\r
- candidates.push_back(to);\r
- std::vector<double> distsq;\r
- distsq.reserve(candidates.size());\r
- for ( unsigned int i = 0; i < candidates.size(); ++i )\r
- {\r
- distsq.push_back( L2sq(c(candidates[i]) - p) );\r
- }\r
- unsigned int closest = 0;\r
- double dsq = distsq[0];\r
- for ( unsigned int i = 1; i < candidates.size(); ++i )\r
- {\r
- if ( dsq > distsq[i] )\r
- {\r
- closest = i;\r
- dsq = distsq[i];\r
- }\r
- }\r
- for ( unsigned int i = 0; i < candidates.size(); ++i )\r
- {\r
- if( distsq[closest] == distsq[i] )\r
- {\r
- result.push_back(candidates[i]);\r
- }\r
- }\r
- return result;\r
-}\r
-\r
-\r
-////////////////////////////////////////////////////////////////////////////////\r
-// Piecewise< D2<SBasis> > versions\r
-\r
-\r
-double nearest_point( Point const& p, \r
- Piecewise< D2<SBasis> > const& c,\r
- double from, double to )\r
-{\r
- if ( from > to ) std::swap(from, to);\r
- if ( from < c.cuts[0] || to > c.cuts[c.size()] )\r
- {\r
- THROW_RANGEERROR("[from,to] interval out of bounds");\r
- }\r
- \r
- unsigned int si = c.segN(from);\r
- unsigned int ei = c.segN(to);\r
- if ( si == ei )\r
- {\r
- double nearest=\r
- nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));\r
- return c.mapToDomain(nearest, si);\r
- }\r
- double t;\r
- double nearest = nearest_point(p, c[si], c.segT(from, si));\r
- unsigned int ni = si;\r
- double dsq;\r
- double mindistsq = distanceSq(p, c[si](nearest));\r
- Rect bb;\r
- for ( unsigned int i = si + 1; i < ei; ++i )\r
- {\r
- bb = bounds_fast(c[i]);\r
- dsq = distanceSq(p, bb);\r
- if ( mindistsq <= dsq ) continue;\r
- t = nearest_point(p, c[i]);\r
- dsq = distanceSq(p, c[i](t));\r
- if ( mindistsq > dsq )\r
- {\r
- nearest = t;\r
- ni = i;\r
- mindistsq = dsq;\r
- }\r
- }\r
- bb = bounds_fast(c[ei]);\r
- dsq = distanceSq(p, bb);\r
- if ( mindistsq > dsq )\r
- {\r
- t = nearest_point(p, c[ei], 0, c.segT(to, ei));\r
- dsq = distanceSq(p, c[ei](t));\r
- if ( mindistsq > dsq )\r
- {\r
- nearest = t;\r
- ni = ei;\r
- }\r
- }\r
- return c.mapToDomain(nearest, ni);\r
-}\r
-\r
-std::vector<double> \r
-all_nearest_points( Point const& p, \r
- Piecewise< D2<SBasis> > const& c, \r
- double from, double to )\r
-{\r
- if ( from > to ) std::swap(from, to);\r
- if ( from < c.cuts[0] || to > c.cuts[c.size()] )\r
- {\r
- THROW_RANGEERROR("[from,to] interval out of bounds");\r
- }\r
- \r
- unsigned int si = c.segN(from);\r
- unsigned int ei = c.segN(to);\r
- if ( si == ei )\r
- {\r
- std::vector<double> all_nearest = \r
- all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));\r
- for ( unsigned int i = 0; i < all_nearest.size(); ++i )\r
- {\r
- all_nearest[i] = c.mapToDomain(all_nearest[i], si);\r
- }\r
- return all_nearest;\r
- }\r
- std::vector<double> all_t;\r
- std::vector< std::vector<double> > all_np;\r
- all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );\r
- std::vector<unsigned int> ni;\r
- ni.push_back(si);\r
- double dsq;\r
- double mindistsq = distanceSq( p, c[si](all_np.front().front()) );\r
- Rect bb;\r
- for ( unsigned int i = si + 1; i < ei; ++i )\r
- {\r
- bb = bounds_fast(c[i]);\r
- dsq = distanceSq(p, bb);\r
- if ( mindistsq < dsq ) continue;\r
- all_t = all_nearest_points(p, c[i]);\r
- dsq = distanceSq( p, c[i](all_t.front()) );\r
- if ( mindistsq > dsq )\r
- {\r
- all_np.clear();\r
- all_np.push_back(all_t);\r
- ni.clear();\r
- ni.push_back(i);\r
- mindistsq = dsq;\r
- }\r
- else if ( mindistsq == dsq )\r
- {\r
- all_np.push_back(all_t);\r
- ni.push_back(i);\r
- }\r
- }\r
- bb = bounds_fast(c[ei]);\r
- dsq = distanceSq(p, bb);\r
- if ( mindistsq >= dsq )\r
- {\r
- all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));\r
- dsq = distanceSq( p, c[ei](all_t.front()) );\r
- if ( mindistsq > dsq )\r
- {\r
- for ( unsigned int i = 0; i < all_t.size(); ++i )\r
- {\r
- all_t[i] = c.mapToDomain(all_t[i], ei);\r
- }\r
- return all_t;\r
- }\r
- else if ( mindistsq == dsq )\r
- {\r
- all_np.push_back(all_t);\r
- ni.push_back(ei);\r
- }\r
- }\r
- std::vector<double> all_nearest;\r
- for ( unsigned int i = 0; i < all_np.size(); ++i )\r
- {\r
- for ( unsigned int j = 0; j < all_np[i].size(); ++j )\r
- {\r
- all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) );\r
- }\r
- }\r
- return all_nearest;\r
-}\r
-\r
-} // end namespace Geom\r
-\r
-\r
+/*
+ * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>
+ *
+ * Authors:
+ *
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2007-2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#include "nearest-point.h"
+
+namespace Geom
+{
+
+////////////////////////////////////////////////////////////////////////////////
+// D2<SBasis> versions
+
+/*
+ * Return the parameter t of a nearest point on the portion of the curve "c",
+ * related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ * The function return the first nearest point to "p" that is found.
+ */
+
+double nearest_point( Point const& p,
+ D2<SBasis> const& c,
+ D2<SBasis> const& dc,
+ double from, double to )
+{
+ if ( from > to ) std::swap(from, to);
+ if ( from < 0 || to > 1 )
+ {
+ THROW_RANGEERROR("[from,to] interval out of bounds");
+ }
+
+ SBasis dd = dot(c - p, dc);
+ std::vector<double> zeros = Geom::roots(dd);
+
+ double closest = from;
+ double min_dist_sq = L2sq(c(from) - p);
+ double distsq;
+ for ( unsigned int i = 0; i < zeros.size(); ++i )
+ {
+ distsq = L2sq(c(zeros[i]) - p);
+ if ( min_dist_sq > L2sq(c(zeros[i]) - p) )
+ {
+ closest = zeros[i];
+ min_dist_sq = distsq;
+ }
+ }
+ if ( min_dist_sq > L2sq( c(to) - p ) )
+ closest = to;
+ return closest;
+
+}
+
+/*
+ * Return the parameters t of all the nearest points on the portion of
+ * the curve "c", related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ */
+
+std::vector<double>
+all_nearest_points( Point const& p,
+ D2<SBasis> const& c,
+ D2<SBasis> const& dc,
+ double from, double to )
+{
+ std::swap(from, to);
+ if ( from > to ) std::swap(from, to);
+ if ( from < 0 || to > 1 )
+ {
+ THROW_RANGEERROR("[from,to] interval out of bounds");
+ }
+
+ std::vector<double> result;
+ SBasis dd = dot(c - p, Geom::derivative(c));
+
+ std::vector<double> zeros = Geom::roots(dd);
+ std::vector<double> candidates;
+ candidates.push_back(from);
+ candidates.insert(candidates.end(), zeros.begin(), zeros.end());
+ candidates.push_back(to);
+ std::vector<double> distsq;
+ distsq.reserve(candidates.size());
+ for ( unsigned int i = 0; i < candidates.size(); ++i )
+ {
+ distsq.push_back( L2sq(c(candidates[i]) - p) );
+ }
+ unsigned int closest = 0;
+ double dsq = distsq[0];
+ for ( unsigned int i = 1; i < candidates.size(); ++i )
+ {
+ if ( dsq > distsq[i] )
+ {
+ closest = i;
+ dsq = distsq[i];
+ }
+ }
+ for ( unsigned int i = 0; i < candidates.size(); ++i )
+ {
+ if( distsq[closest] == distsq[i] )
+ {
+ result.push_back(candidates[i]);
+ }
+ }
+ return result;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Piecewise< D2<SBasis> > versions
+
+
+double nearest_point( Point const& p,
+ Piecewise< D2<SBasis> > const& c,
+ double from, double to )
+{
+ if ( from > to ) std::swap(from, to);
+ if ( from < c.cuts[0] || to > c.cuts[c.size()] )
+ {
+ THROW_RANGEERROR("[from,to] interval out of bounds");
+ }
+
+ unsigned int si = c.segN(from);
+ unsigned int ei = c.segN(to);
+ if ( si == ei )
+ {
+ double nearest=
+ nearest_point(p, c[si], c.segT(from, si), c.segT(to, si));
+ return c.mapToDomain(nearest, si);
+ }
+ double t;
+ double nearest = nearest_point(p, c[si], c.segT(from, si));
+ unsigned int ni = si;
+ double dsq;
+ double mindistsq = distanceSq(p, c[si](nearest));
+ Rect bb;
+ for ( unsigned int i = si + 1; i < ei; ++i )
+ {
+ bb = bounds_fast(c[i]);
+ dsq = distanceSq(p, bb);
+ if ( mindistsq <= dsq ) continue;
+ t = nearest_point(p, c[i]);
+ dsq = distanceSq(p, c[i](t));
+ if ( mindistsq > dsq )
+ {
+ nearest = t;
+ ni = i;
+ mindistsq = dsq;
+ }
+ }
+ bb = bounds_fast(c[ei]);
+ dsq = distanceSq(p, bb);
+ if ( mindistsq > dsq )
+ {
+ t = nearest_point(p, c[ei], 0, c.segT(to, ei));
+ dsq = distanceSq(p, c[ei](t));
+ if ( mindistsq > dsq )
+ {
+ nearest = t;
+ ni = ei;
+ }
+ }
+ return c.mapToDomain(nearest, ni);
+}
+
+std::vector<double>
+all_nearest_points( Point const& p,
+ Piecewise< D2<SBasis> > const& c,
+ double from, double to )
+{
+ if ( from > to ) std::swap(from, to);
+ if ( from < c.cuts[0] || to > c.cuts[c.size()] )
+ {
+ THROW_RANGEERROR("[from,to] interval out of bounds");
+ }
+
+ unsigned int si = c.segN(from);
+ unsigned int ei = c.segN(to);
+ if ( si == ei )
+ {
+ std::vector<double> all_nearest =
+ all_nearest_points(p, c[si], c.segT(from, si), c.segT(to, si));
+ for ( unsigned int i = 0; i < all_nearest.size(); ++i )
+ {
+ all_nearest[i] = c.mapToDomain(all_nearest[i], si);
+ }
+ return all_nearest;
+ }
+ std::vector<double> all_t;
+ std::vector< std::vector<double> > all_np;
+ all_np.push_back( all_nearest_points(p, c[si], c.segT(from, si)) );
+ std::vector<unsigned int> ni;
+ ni.push_back(si);
+ double dsq;
+ double mindistsq = distanceSq( p, c[si](all_np.front().front()) );
+ Rect bb;
+ for ( unsigned int i = si + 1; i < ei; ++i )
+ {
+ bb = bounds_fast(c[i]);
+ dsq = distanceSq(p, bb);
+ if ( mindistsq < dsq ) continue;
+ all_t = all_nearest_points(p, c[i]);
+ dsq = distanceSq( p, c[i](all_t.front()) );
+ if ( mindistsq > dsq )
+ {
+ all_np.clear();
+ all_np.push_back(all_t);
+ ni.clear();
+ ni.push_back(i);
+ mindistsq = dsq;
+ }
+ else if ( mindistsq == dsq )
+ {
+ all_np.push_back(all_t);
+ ni.push_back(i);
+ }
+ }
+ bb = bounds_fast(c[ei]);
+ dsq = distanceSq(p, bb);
+ if ( mindistsq >= dsq )
+ {
+ all_t = all_nearest_points(p, c[ei], 0, c.segT(to, ei));
+ dsq = distanceSq( p, c[ei](all_t.front()) );
+ if ( mindistsq > dsq )
+ {
+ for ( unsigned int i = 0; i < all_t.size(); ++i )
+ {
+ all_t[i] = c.mapToDomain(all_t[i], ei);
+ }
+ return all_t;
+ }
+ else if ( mindistsq == dsq )
+ {
+ all_np.push_back(all_t);
+ ni.push_back(ei);
+ }
+ }
+ std::vector<double> all_nearest;
+ for ( unsigned int i = 0; i < all_np.size(); ++i )
+ {
+ for ( unsigned int j = 0; j < all_np[i].size(); ++j )
+ {
+ all_nearest.push_back( c.mapToDomain(all_np[i][j], ni[i]) );
+ }
+ }
+ return all_nearest;
+}
+
+} // end namespace Geom
+
+
index 73ac0c3ce73f9785a68d19d3afb9f83bdf310f79..8c1a6f7d7be90dde4374ea5e63b2cb89f0ca2efd 100644 (file)
-/*\r
- * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>\r
- *\r
- *\r
- * Authors:\r
- * \r
- * Marco Cecchetti <mrcekets at gmail.com>\r
- * \r
- * Copyright 2007-2008 authors\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#ifndef _NEAREST_POINT_H_\r
-#define _NEAREST_POINT_H_\r
-\r
-\r
-#include <vector>\r
-\r
-#include "d2.h"\r
-#include "piecewise.h"\r
-#include "exception.h"\r
-\r
-\r
-\r
-namespace Geom\r
-{\r
-\r
-/*\r
- * Given a line L specified by a point A and direction vector v,\r
- * return the point on L nearest to p. Note that the returned value\r
- * is with respect to the _normalized_ direction of v!\r
- */\r
-inline double nearest_point(Point const &p, Point const &A, Point const &v)\r
-{\r
- Point d(p - A);\r
- return d[0] * v[0] + d[1] * v[1];\r
-}\r
-\r
-////////////////////////////////////////////////////////////////////////////////\r
-// D2<SBasis> versions\r
-\r
-/*\r
- * Return the parameter t of a nearest point on the portion of the curve "c", \r
- * related to the interval [from, to], to the point "p".\r
- * The needed curve derivative "dc" is passed as parameter.\r
- * The function return the first nearest point to "p" that is found.\r
- */\r
-double nearest_point( Point const& p,\r
- D2<SBasis> const& c, D2<SBasis> const& dc, \r
- double from = 0, double to = 1 );\r
-\r
-inline\r
-double nearest_point( Point const& p, \r
- D2<SBasis> const& c, \r
- double from = 0, double to = 1 )\r
-{\r
- return nearest_point(p, c, Geom::derivative(c), from, to);\r
-}\r
-\r
-/*\r
- * Return the parameters t of all the nearest points on the portion of \r
- * the curve "c", related to the interval [from, to], to the point "p".\r
- * The needed curve derivative "dc" is passed as parameter.\r
- */\r
-std::vector<double> \r
-all_nearest_points( Point const& p, \r
- D2<SBasis> const& c, D2<SBasis> const& dc, \r
- double from = 0, double to = 1 );\r
-\r
-inline\r
-std::vector<double> \r
-all_nearest_points( Point const& p, \r
- D2<SBasis> const& c, \r
- double from = 0, double to = 1 )\r
-{\r
- return all_nearest_points(p, c, Geom::derivative(c), from, to);\r
-}\r
-\r
-\r
-////////////////////////////////////////////////////////////////////////////////\r
-// Piecewise< D2<SBasis> > versions\r
-\r
-double nearest_point( Point const& p, \r
- Piecewise< D2<SBasis> > const& c, \r
- double from, double to );\r
-\r
-inline\r
-double nearest_point( Point const& p, Piecewise< D2<SBasis> > const& c ) \r
-{\r
- return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]);\r
-}\r
-\r
-\r
-std::vector<double> \r
-all_nearest_points( Point const& p, \r
- Piecewise< D2<SBasis> > const& c, \r
- double from, double to );\r
-\r
-inline\r
-std::vector<double> \r
-all_nearest_points( Point const& p, Piecewise< D2<SBasis> > const& c ) \r
-{\r
- return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]);\r
-}\r
-\r
-} // end namespace Geom\r
-\r
-\r
-\r
-#endif /*_NEAREST_POINT_H_*/\r
+/*
+ * nearest point routines for D2<SBasis> and Piecewise<D2<SBasis>>
+ *
+ *
+ * Authors:
+ *
+ * Marco Cecchetti <mrcekets at gmail.com>
+ *
+ * Copyright 2007-2008 authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+
+#ifndef _NEAREST_POINT_H_
+#define _NEAREST_POINT_H_
+
+
+#include <vector>
+
+#include "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<SBasis> versions
+
+/*
+ * Return the parameter t of a nearest point on the portion of the curve "c",
+ * related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ * The function return the first nearest point to "p" that is found.
+ */
+double nearest_point( Point const& p,
+ D2<SBasis> const& c, D2<SBasis> const& dc,
+ double from = 0, double to = 1 );
+
+inline
+double nearest_point( Point const& p,
+ D2<SBasis> const& c,
+ double from = 0, double to = 1 )
+{
+ return nearest_point(p, c, Geom::derivative(c), from, to);
+}
+
+/*
+ * Return the parameters t of all the nearest points on the portion of
+ * the curve "c", related to the interval [from, to], to the point "p".
+ * The needed curve derivative "dc" is passed as parameter.
+ */
+std::vector<double>
+all_nearest_points( Point const& p,
+ D2<SBasis> const& c, D2<SBasis> const& dc,
+ double from = 0, double to = 1 );
+
+inline
+std::vector<double>
+all_nearest_points( Point const& p,
+ D2<SBasis> const& c,
+ double from = 0, double to = 1 )
+{
+ return all_nearest_points(p, c, Geom::derivative(c), from, to);
+}
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Piecewise< D2<SBasis> > versions
+
+double nearest_point( Point const& p,
+ Piecewise< D2<SBasis> > const& c,
+ double from, double to );
+
+inline
+double nearest_point( Point const& p, Piecewise< D2<SBasis> > const& c )
+{
+ return nearest_point(p, c, c.cuts[0], c.cuts[c.size()]);
+}
+
+
+std::vector<double>
+all_nearest_points( Point const& p,
+ Piecewise< D2<SBasis> > const& c,
+ double from, double to );
+
+inline
+std::vector<double>
+all_nearest_points( Point const& p, Piecewise< D2<SBasis> > const& c )
+{
+ return all_nearest_points(p, c, c.cuts[0], c.cuts[c.size()]);
+}
+
+} // end namespace Geom
+
+
+
+#endif /*_NEAREST_POINT_H_*/
index 2b49639983fde5f446e85e399d72f1e06ab24abe..9cb521eb2a2f919a8d7982d46117109362c5643f 100644 (file)
-#ifndef _NL_LINEAR_SYSTEM_H_\r
-#define _NL_LINEAR_SYSTEM_H_\r
-\r
-\r
-#include <cassert>\r
-\r
-#include <gsl/gsl_linalg.h>\r
-\r
-#include "numeric/matrix.h"\r
-#include "numeric/vector.h"\r
-\r
-\r
-namespace Geom { namespace NL {\r
-\r
-\r
-class LinearSystem\r
-{\r
-public:\r
- LinearSystem(Matrix & _matrix, Vector & _vector)\r
- : m_matrix(_matrix), m_vector(_vector), m_solution(_matrix.columns())\r
- {\r
- }\r
- \r
- const Vector & LU_solve()\r
- {\r
- assert( matrix().rows() == matrix().columns() \r
- && matrix().rows() == vector().size() );\r
- int s;\r
- gsl_permutation * p = gsl_permutation_alloc(matrix().rows());\r
- gsl_linalg_LU_decomp (matrix().get_gsl_matrix(), p, &s);\r
- gsl_linalg_LU_solve( matrix().get_gsl_matrix(), \r
- p, \r
- vector().get_gsl_vector(), \r
- m_solution.get_gsl_vector()\r
- );\r
- gsl_permutation_free(p);\r
- return solution();\r
- }\r
- \r
- const Vector & SV_solve()\r
- {\r
- assert( matrix().rows() >= matrix().columns()\r
- && matrix().rows() == vector().size() );\r
- \r
- gsl_matrix* U = matrix().get_gsl_matrix();\r
- gsl_matrix* V = gsl_matrix_alloc(matrix().columns(), matrix().columns());\r
- gsl_vector* S = gsl_vector_alloc(matrix().columns());\r
- gsl_vector* work = gsl_vector_alloc(matrix().columns());\r
- \r
- gsl_linalg_SV_decomp( U, V, S, work );\r
- \r
- gsl_vector* b = vector().get_gsl_vector();\r
- gsl_vector* x = m_solution.get_gsl_vector();\r
- \r
- gsl_linalg_SV_solve( U, V, S, b, x);\r
- \r
- gsl_matrix_free(V);\r
- gsl_vector_free(S);\r
- gsl_vector_free(work);\r
- \r
- return solution(); \r
- }\r
- \r
- Matrix & matrix()\r
- {\r
- return m_matrix;\r
- }\r
- \r
- Vector & vector()\r
- {\r
- return m_vector;\r
- }\r
- \r
- const Vector & solution() const\r
- {\r
- return m_solution;\r
- }\r
- \r
-private:\r
- Matrix & m_matrix;\r
- Vector & m_vector;\r
- Vector m_solution;\r
-};\r
-\r
-\r
-} } // end namespaces\r
-\r
-\r
-#endif /*_NL_LINEAR_SYSTEM_H_*/\r
+#ifndef _NL_LINEAR_SYSTEM_H_
+#define _NL_LINEAR_SYSTEM_H_
+
+
+#include <cassert>
+
+#include <gsl/gsl_linalg.h>
+
+#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_*/
index bdfab75ef7d45c34d2ba92f87c041e40dc193596..bdd647e53167fad3d34d143122d9497ebfdcd45c 100644 (file)
-#ifndef _NL_MATRIX_H_\r
-#define _NL_MATRIX_H_\r
-\r
-\r
-#include <cassert>\r
-#include <utility>\r
-#include <algorithm>\r
-\r
-#include <gsl/gsl_matrix.h>\r
-\r
-\r
-\r
-namespace Geom { namespace NL {\r
-\r
-class Matrix;\r
-void swap(Matrix & m1, Matrix & m2);\r
-\r
-\r
-class Matrix\r
-{\r
-public:\r
- // the matrix is not inizialized\r
- Matrix(size_t n1, size_t n2)\r
- : m_rows(n1), m_columns(n2)\r
- {\r
- m_matrix = gsl_matrix_alloc(n1, n2);\r
- }\r
- \r
- Matrix(size_t n1, size_t n2, double x)\r
- :m_rows(n1), m_columns(n2)\r
- {\r
- m_matrix = gsl_matrix_alloc(n1, n2);\r
- gsl_matrix_set_all(m_matrix, x);\r
- }\r
- \r
- Matrix( Matrix const& _matrix )\r
- : m_rows(_matrix.rows()), m_columns(_matrix.columns())\r
- {\r
- m_matrix = gsl_matrix_alloc(rows(), columns());\r
- gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);\r
- }\r
- \r
-// explicit\r
-// Matrix( gsl_matrix* m, size_t n1, size_t n2)\r
-// : m_rows(n1), m_columns(n2), m_matrix(m)\r
-// { \r
-// }\r
- \r
- Matrix & operator=(Matrix const& _matrix)\r
- {\r
- assert( rows() == _matrix.rows() && columns() == _matrix.columns() );\r
- gsl_matrix_memcpy(m_matrix, _matrix.m_matrix);\r
- \r
- return *this;\r
- }\r
- \r
- virtual ~Matrix()\r
- {\r
- gsl_matrix_free(m_matrix);\r
- }\r
- \r
- void set_all( double x = 0)\r
- {\r
- gsl_matrix_set_all(m_matrix, x);\r
- }\r
- \r
- void set_identity()\r
- {\r
- gsl_matrix_set_identity(m_matrix);\r
- }\r
- \r
- const double & operator() (size_t i, size_t j) const\r
- {\r
- return *gsl_matrix_const_ptr(m_matrix, i, j);\r
- }\r
- \r
- double & operator() (size_t i, size_t j)\r
- {\r
- return *gsl_matrix_ptr(m_matrix, i, j);\r
- }\r
- \r
- gsl_matrix* get_gsl_matrix()\r
- {\r
- return m_matrix;\r
- }\r
- \r
- void swap_rows(size_t i, size_t j)\r
- {\r
- gsl_matrix_swap_rows(m_matrix, i, j);\r
- }\r
- \r
- void swap_columns(size_t i, size_t j)\r
- {\r
- gsl_matrix_swap_columns(m_matrix, i, j);\r
- }\r
- \r
- Matrix & transpose()\r
- {\r
- gsl_matrix_transpose(m_matrix);\r
- return (*this);\r
- }\r
- \r
- Matrix & scale(double x)\r
- {\r
- gsl_matrix_scale(m_matrix, x);\r
- return (*this);\r
- }\r
- \r
- Matrix & translate(double x)\r
- {\r
- gsl_matrix_add_constant(m_matrix, x);\r
- return (*this);\r
- }\r
- \r
- Matrix & operator+=(Matrix const& _matrix)\r
- {\r
- gsl_matrix_add(m_matrix, _matrix.m_matrix);\r
- return (*this);\r
- }\r
- \r
- Matrix & operator-=(Matrix const& _matrix)\r
- {\r
- gsl_matrix_sub(m_matrix, _matrix.m_matrix);\r
- return (*this);\r
- }\r
- \r
- bool is_zero() const\r
- {\r
- return gsl_matrix_isnull(m_matrix);\r
- }\r
- \r
- bool is_positive() const\r
- {\r
- return gsl_matrix_ispos(m_matrix);\r
- }\r
- \r
- bool is_negative() const\r
- {\r
- return gsl_matrix_isneg(m_matrix);\r
- }\r
- \r
- bool is_non_negative() const\r
- {\r
- for ( unsigned int i = 0; i < rows(); ++i )\r
- {\r
- for ( unsigned int j = 0; j < columns(); ++j )\r
- {\r
- if ( (*this)(i,j) < 0 ) return false;\r
- }\r
- }\r
- return true;\r
- }\r
- \r
- double max() const\r
- {\r
- return gsl_matrix_max(m_matrix);\r
- }\r
- \r
- double min() const\r
- {\r
- return gsl_matrix_min(m_matrix);\r
- }\r
- \r
- std::pair<size_t, size_t>\r
- max_index() const\r
- {\r
- std::pair<size_t, size_t> indices;\r
- gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));\r
- return indices;\r
- }\r
- \r
- std::pair<size_t, size_t>\r
- min_index() const\r
- {\r
- std::pair<size_t, size_t> indices;\r
- gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));\r
- return indices;\r
- }\r
- \r
- friend \r
- void swap(Matrix & m1, Matrix & m2);\r
- \r
- size_t rows() const\r
- {\r
- return m_rows;\r
- }\r
- \r
- size_t columns() const\r
- {\r
- return m_columns;\r
- }\r
- \r
-private:\r
- size_t m_rows, m_columns;\r
- gsl_matrix* m_matrix;\r
-};\r
-\r
-void swap(Matrix & m1, Matrix & m2)\r
-{\r
- assert( m1.rows() == m2.rows() && m1.columns() == m2.columns() );\r
- std::swap(m1.m_matrix, m2.m_matrix);\r
-}\r
-\r
-\r
-} } // end namespaces\r
-\r
-#endif /*_NL_MATRIX_H_*/\r
+#ifndef _NL_MATRIX_H_
+#define _NL_MATRIX_H_
+
+
+#include <cassert>
+#include <utility>
+#include <algorithm>
+
+#include <gsl/gsl_matrix.h>
+
+
+
+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<size_t, size_t>
+ max_index() const
+ {
+ std::pair<size_t, size_t> indices;
+ gsl_matrix_max_index(m_matrix, &(indices.first), &(indices.second));
+ return indices;
+ }
+
+ std::pair<size_t, size_t>
+ min_index() const
+ {
+ std::pair<size_t, size_t> indices;
+ gsl_matrix_min_index(m_matrix, &(indices.first), &(indices.second));
+ return indices;
+ }
+
+ 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_*/
index 136b2cc3fdd47db1b46bff03b59812b41ec2bee8..b25861e22c7d353fac06412d462dd8e1869856ec 100644 (file)
-#ifndef _NL_VECTOR_H_\r
-#define _NL_VECTOR_H_\r
-\r
-#include <cassert>\r
-#include <utility>\r
-\r
-#include <gsl/gsl_vector.h>\r
-\r
-\r
-namespace Geom { namespace NL {\r
-\r
-class Vector;\r
-void swap(Vector & v1, Vector & v2);\r
-\r
-\r
-class Vector\r
-{\r
-public:\r
- Vector(size_t n)\r
- : m_size(n)\r
- {\r
- m_vector = gsl_vector_alloc(n);\r
- }\r
- \r
- Vector(size_t n, double x)\r
- : m_size(n)\r
- {\r
- m_vector = gsl_vector_alloc(n);\r
- gsl_vector_set_all(m_vector, x);\r
- }\r
- \r
- // create a vector with n elements all set to zero \r
- // but the i-th that is set to 1\r
- Vector(size_t n, size_t i)\r
- : m_size(n)\r
- {\r
- m_vector = gsl_vector_alloc(n);\r
- gsl_vector_set_basis(m_vector, i);\r
- }\r
- \r
- Vector(Vector const& _vector)\r
- : m_size(_vector.size())\r
- {\r
- m_vector = gsl_vector_alloc(size());\r
- gsl_vector_memcpy(m_vector, _vector.m_vector);\r
- }\r
- \r
- virtual ~Vector()\r
- {\r
- gsl_vector_free(m_vector);\r
- }\r
- \r
- Vector & operator=(Vector const& _vector)\r
- {\r
- assert( size() == _vector.size() );\r
- gsl_vector_memcpy(m_vector, _vector.m_vector);\r
- }\r
- \r
- void set_all(double x = 0)\r
- {\r
- gsl_vector_set_all(m_vector, x);\r
- }\r
-\r
- void set_basis(size_t i)\r
- {\r
- gsl_vector_set_basis(m_vector, i);\r
- }\r
- \r
- double const& operator[](size_t i) const\r
- {\r
- return *gsl_vector_const_ptr(m_vector, i);\r
- }\r
- \r
- double & operator[](size_t i)\r
- {\r
- return *gsl_vector_ptr(m_vector, i);\r
- }\r
- \r
- gsl_vector* get_gsl_vector()\r
- {\r
- return m_vector;\r
- }\r
- \r
- void swap_elements(size_t i, size_t j)\r
- {\r
- gsl_vector_swap_elements(m_vector, i, j);\r
- }\r
- \r
- void reverse()\r
- {\r
- gsl_vector_reverse(m_vector);\r
- }\r
- \r
- Vector & scale(double x)\r
- {\r
- gsl_vector_scale(m_vector, x);\r
- return (*this);\r
- }\r
- \r
- Vector & translate(double x)\r
- {\r
- gsl_vector_add_constant(m_vector, x);\r
- return (*this);\r
- }\r
- \r
- Vector & operator+=(Vector const& _vector)\r
- {\r
- gsl_vector_add(m_vector, _vector.m_vector);\r
- return (*this);\r
- }\r
- \r
- Vector & operator-=(Vector const& _vector)\r
- {\r
- gsl_vector_sub(m_vector, _vector.m_vector);\r
- return (*this);\r
- }\r
- \r
- bool is_zero() const\r
- {\r
- return gsl_vector_isnull(m_vector);\r
- }\r
- \r
- bool is_positive() const\r
- {\r
- return gsl_vector_ispos(m_vector);\r
- }\r
- \r
- bool is_negative() const\r
- {\r
- return gsl_vector_isneg(m_vector);\r
- }\r
- \r
- bool is_non_negative() const\r
- {\r
- for ( size_t i = 0; i < size(); ++i )\r
- {\r
- if ( (*this)[i] < 0 ) return false;\r
- }\r
- return true;\r
- }\r
- \r
- double max() const\r
- {\r
- return gsl_vector_max(m_vector);\r
- }\r
- \r
- double min() const\r
- {\r
- return gsl_vector_min(m_vector);\r
- }\r
- \r
- size_t max_index() const\r
- {\r
- return gsl_vector_max_index(m_vector);\r
- }\r
- \r
- size_t min_index() const\r
- {\r
- return gsl_vector_min_index(m_vector);\r
- }\r
- \r
- friend\r
- void swap(Vector & v1, Vector & v2);\r
- \r
- size_t size() const\r
- {\r
- return m_size;\r
- }\r
- \r
-private:\r
- size_t m_size;\r
- gsl_vector* m_vector;\r
-};\r
-\r
-void swap(Vector & v1, Vector & v2)\r
-{\r
- assert( v1.size() == v2.size() );\r
- std::swap(v1.m_vector, v2.m_vector);\r
-}\r
-\r
-} } // end namespaces\r
-\r
-\r
-#endif /*_NL_VECTOR_H_*/\r
+#ifndef _NL_VECTOR_H_
+#define _NL_VECTOR_H_
+
+#include <cassert>
+#include <utility>
+
+#include <gsl/gsl_vector.h>
+
+
+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 0493b4a593555d0d7641776d6a41b793ecd496b7..abed8604ab963243492b2df0bb3f773a2832ab89 100644 (file)
--- a/src/2geom/path.h
+++ b/src/2geom/path.h
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)
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)
Curve* duplicate() const
{
- return new SVGEllipticalArc(*this);
+ return new EllipticalArc(*this);
}
double center(unsigned int i) const
return pointAtAngle(tt);
}
- std::pair<SVGEllipticalArc, SVGEllipticalArc>
+ std::pair<EllipticalArc, EllipticalArc>
subdivide(Coord t) const
{
- SVGEllipticalArc* arc1 = static_cast<SVGEllipticalArc*>(portion(0, t));
- SVGEllipticalArc* arc2 = static_cast<SVGEllipticalArc*>(portion(t, 1));
+ EllipticalArc* arc1 = static_cast<EllipticalArc*>(portion(0, t));
+ EllipticalArc* arc2 = static_cast<EllipticalArc*>(portion(t, 1));
assert( arc1 != NULL && arc2 != NULL);
- std::pair<SVGEllipticalArc, SVGEllipticalArc> arc_pair(*arc1, *arc2);
+ std::pair<EllipticalArc, EllipticalArc> arc_pair(*arc1, *arc2);
delete arc1;
delete arc2;
return arc_pair;
// 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;
double m_start_angle, m_end_angle;
Point m_center;
-}; // end class SVGEllipticalArc
+}; // end class EllipticalArc
return old;
}
+ BaseIterator &operator--() {
+ --impl_;
+ return *this;
+ }
+
+ BaseIterator operator--(int) {
+ BaseIterator old=*this;
+ --(*this);
+ return old;
+ }
+
private:
BaseIterator(IteratorImpl const &pos) : impl_(pos) {}
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(); }
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<SBasis> const &curve);
void append(Path const &other);
diff --git a/src/2geom/poly.cpp b/src/2geom/poly.cpp
index 17f28649401d43bce4939de5b45321371666b333..4f5ba6c557ea25227c6672c4dd9327ad2a200447 100644 (file)
--- a/src/2geom/poly.cpp
+++ b/src/2geom/poly.cpp
}
return result;
}
+#define HAVE_GSL
#ifdef HAVE_GSL
#include <gsl/gsl_poly.h>
#endif
diff --git a/src/2geom/sbasis.cpp b/src/2geom/sbasis.cpp
index c6013204374d19a2b69bb506d4e535fab6aeaf28..920fd37fe7ed29b66c58d5cfa904b18a806c1326 100644 (file)
--- a/src/2geom/sbasis.cpp
+++ b/src/2geom/sbasis.cpp
}
std::vector<double> SBasis::valueAndDerivatives(double t, unsigned n) const {
- std::vector<double> ret;
+ std::vector<double> 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;
}
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)}
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]);
//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));
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;
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;
}
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 7f0d88f8cd2211e7443e3bcc6f71ec935a07d9fa..0a48ff1cf78039b0c4fe22d686a6ea8ca51b1dc5 100644 (file)
--- a/src/2geom/sbasis.h
+++ b/src/2geom/sbasis.h
SBasis(Linear const & bo) {
push_back(bo);
}
+ SBasis(Linear* bo) {
+ push_back(*bo);
+ }
//IMPL: FragmentConcept
typedef double output_type;
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);
}
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
}
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);
index 1338faa7c5aa6fab05a14a46842331bd730c59db..cd7c36cc33e7a71bcff15e9e34f504ccfda2ba6b 100644 (file)
template<class t>
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<double> &solutions;
+ Bernsteins(int degr, std::vector<double> &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.
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) {
}
}
- 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);
}
/*
* 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] */
/* 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;
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)
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-elliptical-arc.cpp b/src/2geom/svg-elliptical-arc.cpp
+++ /dev/null
@@ -1,931 +0,0 @@
-/*\r
- * SVG Elliptical Arc Class\r
- *\r
- * Copyright 2008 Marco Cecchetti <mrcekets at gmail.com>\r
- *\r
- * This library is free software; you can redistribute it and/or\r
- * modify it either under the terms of the GNU Lesser General Public\r
- * License version 2.1 as published by the Free Software Foundation\r
- * (the "LGPL") or, at your option, under the terms of the Mozilla\r
- * Public License Version 1.1 (the "MPL"). If you do not alter this\r
- * notice, a recipient may use your version of this file under either\r
- * the MPL or the LGPL.\r
- *\r
- * You should have received a copy of the LGPL along with this library\r
- * in the file COPYING-LGPL-2.1; if not, write to the Free Software\r
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA\r
- * You should have received a copy of the MPL along with this library\r
- * in the file COPYING-MPL-1.1\r
- *\r
- * The contents of this file are subject to the Mozilla Public License\r
- * Version 1.1 (the "License"); you may not use this file except in\r
- * compliance with the License. You may obtain a copy of the License at\r
- * http://www.mozilla.org/MPL/\r
- *\r
- * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY\r
- * OF ANY KIND, either express or implied. See the LGPL or the MPL for\r
- * the specific language governing rights and limitations.\r
- */\r
-\r
-\r
-#include "path.h"\r
-#include "angle.h"\r
-\r
-#include <gsl/gsl_poly.h>\r
-\r
-#include <cfloat>\r
-\r
-\r
-\r
-\r
-namespace Geom\r
-{\r
-\r
-\r
-Rect SVGEllipticalArc::boundsExact() const\r
-{\r
- std::vector<double> extremes(4);\r
- double cosrot = std::cos(rotation_angle());\r
- double sinrot = std::sin(rotation_angle());\r
- extremes[0] = std::atan2( -ray(Y) * sinrot, ray(X) * cosrot );\r
- extremes[1] = extremes[0] + M_PI;\r
- if ( extremes[0] < 0 ) extremes[0] += 2*M_PI; \r
- extremes[2] = std::atan2( ray(Y) * cosrot, ray(X) * sinrot );\r
- extremes[3] = extremes[2] + M_PI;\r
- if ( extremes[2] < 0 ) extremes[2] += 2*M_PI;\r
- \r
- \r
- std::vector<double>arc_extremes(4);\r
- arc_extremes[0] = initialPoint()[X];\r
- arc_extremes[1] = finalPoint()[X];\r
- if ( arc_extremes[0] < arc_extremes[1] ) \r
- std::swap(arc_extremes[0], arc_extremes[1]);\r
- arc_extremes[2] = initialPoint()[Y];\r
- arc_extremes[3] = finalPoint()[Y];\r
- if ( arc_extremes[2] < arc_extremes[3] ) \r
- std::swap(arc_extremes[2], arc_extremes[3]);\r
- \r
- \r
- if ( start_angle() < end_angle() )\r
- {\r
- if ( sweep_flag() )\r
- {\r
- for ( unsigned int i = 0; i < extremes.size(); ++i )\r
- {\r
- if ( start_angle() < extremes[i] && extremes[i] < end_angle() )\r
- {\r
- arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
- }\r
- }\r
- }\r
- else\r
- {\r
- for ( unsigned int i = 0; i < extremes.size(); ++i )\r
- {\r
- if ( start_angle() > extremes[i] || extremes[i] > end_angle() )\r
- {\r
- arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
- }\r
- }\r
- }\r
- }\r
- else\r
- {\r
- if ( sweep_flag() )\r
- {\r
- for ( unsigned int i = 0; i < extremes.size(); ++i )\r
- {\r
- if ( start_angle() < extremes[i] || extremes[i] < end_angle() )\r
- {\r
- arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
- }\r
- }\r
- }\r
- else\r
- {\r
- for ( unsigned int i = 0; i < extremes.size(); ++i )\r
- {\r
- if ( start_angle() > extremes[i] && extremes[i] > end_angle() )\r
- {\r
- arc_extremes[i] = pointAtAngle(extremes[i])[i >> 1];\r
- }\r
- } \r
- }\r
- }\r
- \r
- return Rect( Point(arc_extremes[1], arc_extremes[3]) , \r
- Point(arc_extremes[0], arc_extremes[2]) );\r
-\r
-}\r
-\r
-\r
-std::vector<double> \r
-SVGEllipticalArc::roots(double v, Dim2 d) const\r
-{\r
- if ( d > Y )\r
- {\r
- THROW_RANGEERROR("dimention out of range");\r
- }\r
- \r
- std::vector<double> sol;\r
- if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
- {\r
- if ( center(d) == v )\r
- sol.push_back(0);\r
- return sol;\r
- }\r
- \r
- const char* msg[2][2] = \r
- {\r
- { "d == X; ray(X) == 0; "\r
- "s = (v - center(X)) / ( -ray(Y) * std::sin(rotation_angle()) ); "\r
- "s should be contained in [-1,1]",\r
- "d == X; ray(Y) == 0; "\r
- "s = (v - center(X)) / ( ray(X) * std::cos(rotation_angle()) ); "\r
- "s should be contained in [-1,1]"\r
- },\r
- { "d == Y; ray(X) == 0; "\r
- "s = (v - center(X)) / ( ray(Y) * std::cos(rotation_angle()) ); "\r
- "s should be contained in [-1,1]",\r
- "d == Y; ray(Y) == 0; "\r
- "s = (v - center(X)) / ( ray(X) * std::sin(rotation_angle()) ); "\r
- "s should be contained in [-1,1]"\r
- },\r
- }; \r
- \r
- for ( unsigned int dim = 0; dim < 2; ++dim )\r
- {\r
- if ( are_near(ray(dim), 0) )\r
- {\r
- \r
- if ( initialPoint()[d] == v && finalPoint()[d] == v )\r
- {\r
- THROW_EXCEPTION("infinite solutions");\r
- }\r
- if ( (initialPoint()[d] < finalPoint()[d])\r
- && (initialPoint()[d] > v || finalPoint()[d] < v) )\r
- {\r
- return sol;\r
- }\r
- if ( (initialPoint()[d] > finalPoint()[d])\r
- && (finalPoint()[d] > v || initialPoint()[d] < v) )\r
- {\r
- return sol;\r
- }\r
- double ray_prj;\r
- switch(d)\r
- {\r
- case X: \r
- switch(dim)\r
- { \r
- case X: ray_prj = -ray(Y) * std::sin(rotation_angle());\r
- break;\r
- case Y: ray_prj = ray(X) * std::cos(rotation_angle());\r
- break;\r
- }\r
- break;\r
- case Y:\r
- switch(dim)\r
- { \r
- case X: ray_prj = ray(Y) * std::cos(rotation_angle());\r
- break;\r
- case Y: ray_prj = ray(X) * std::sin(rotation_angle());\r
- break;\r
- }\r
- break;\r
- }\r
- \r
- double s = (v - center(d)) / ray_prj;\r
- if ( s < -1 || s > 1 )\r
- {\r
- THROW_LOGICALERROR(msg[d][dim]);\r
- }\r
- switch(dim)\r
- { \r
- case X: \r
- s = std::asin(s); // return a value in [-PI/2,PI/2]\r
- if ( logical_xor( sweep_flag(), are_near(start_angle(), M_PI/2) ) )\r
- {\r
- if ( s < 0 ) s += 2*M_PI;\r
- }\r
- else\r
- {\r
- s = M_PI - s;\r
- if (!(s < 2*M_PI) ) s -= 2*M_PI;\r
- }\r
- break;\r
- case Y: \r
- s = std::acos(s); // return a value in [0,PI]\r
- if ( logical_xor( sweep_flag(), are_near(start_angle(), 0) ) )\r
- {\r
- s = 2*M_PI - s;\r
- if ( !(s < 2*M_PI) ) s -= 2*M_PI;\r
- }\r
- break;\r
- }\r
- \r
- //std::cerr << "s = " << rad_to_deg(s);\r
- s = map_to_01(s);\r
- //std::cerr << " -> t: " << s << std::endl;\r
- if ( !(s < 0 || s > 1) )\r
- sol.push_back(s);\r
- return sol;\r
- }\r
- }\r
- \r
- double rotx, roty;\r
- switch(d)\r
- {\r
- case X: \r
- rotx = std::cos(rotation_angle());\r
- roty = -std::sin(rotation_angle());\r
- break;\r
- case Y:\r
- rotx = std::sin(rotation_angle());\r
- roty = std::cos(rotation_angle());\r
- break;\r
- }\r
- double rxrotx = ray(X) * rotx;\r
- double c_v = center(d) - v;\r
-\r
- double a = -rxrotx + c_v;\r
- double b = ray(Y) * roty;\r
- double c = rxrotx + c_v;\r
- //std::cerr << "a = " << a << std::endl;\r
- //std::cerr << "b = " << b << std::endl;\r
- //std::cerr << "c = " << c << std::endl;\r
- \r
- if ( are_near(a,0) )\r
- {\r
- sol.push_back(M_PI);\r
- if ( !are_near(b,0) )\r
- {\r
- double s = 2 * std::atan(-c/(2*b));\r
- if ( s < 0 ) s += 2*M_PI;\r
- sol.push_back(s);\r
- }\r
- }\r
- else\r
- {\r
- double delta = b * b - a * c;\r
- //std::cerr << "delta = " << delta << std::endl;\r
- if ( are_near(delta, 0) )\r
- {\r
- double s = 2 * std::atan(-b/a);\r
- if ( s < 0 ) s += 2*M_PI;\r
- sol.push_back(s);\r
- }\r
- else if ( delta > 0 )\r
- { \r
- double sq = std::sqrt(delta);\r
- double s = 2 * std::atan( (-b - sq) / a );\r
- if ( s < 0 ) s += 2*M_PI;\r
- sol.push_back(s);\r
- s = 2 * std::atan( (-b + sq) / a );\r
- if ( s < 0 ) s += 2*M_PI;\r
- sol.push_back(s);\r
- }\r
- }\r
- \r
- std::vector<double> arc_sol;\r
- for (unsigned int i = 0; i < sol.size(); ++i )\r
- {\r
- //std::cerr << "s = " << rad_to_deg(sol[i]);\r
- sol[i] = map_to_01(sol[i]);\r
- //std::cerr << " -> t: " << sol[i] << std::endl;\r
- if ( !(sol[i] < 0 || sol[i] > 1) )\r
- arc_sol.push_back(sol[i]);\r
- }\r
- return arc_sol;\r
- \r
- \r
-// return SBasisCurve(toSBasis()).roots(v, d);\r
-}\r
-\r
-// D(E(t,C),t) = E(t+PI/2,O)\r
-Curve* SVGEllipticalArc::derivative() const\r
-{\r
- SVGEllipticalArc* result = new SVGEllipticalArc(*this);\r
- result->m_center[X] = result->m_center[Y] = 0;\r
- result->m_start_angle += M_PI/2;\r
- if( !( result->m_start_angle < 2*M_PI ) )\r
- {\r
- result->m_start_angle -= 2*M_PI;\r
- }\r
- result->m_end_angle += M_PI/2;\r
- if( !( result->m_end_angle < 2*M_PI ) )\r
- {\r
- result->m_end_angle -= 2*M_PI;\r
- }\r
- result->m_initial_point = result->pointAtAngle( result->start_angle() );\r
- result->m_final_point = result->pointAtAngle( result->end_angle() );\r
- return result;\r
- \r
-}\r
-\r
-std::vector<Point> \r
-SVGEllipticalArc::pointAndDerivatives(Coord t, unsigned int n) const\r
-{\r
- std::vector<Point> result;\r
- result.reserve(n);\r
- double angle = map_unit_interval_on_circular_arc(t, start_angle(), \r
- end_angle(), sweep_flag());\r
- SVGEllipticalArc ea(*this);\r
- ea.m_center = Point(0,0);\r
- unsigned int m = std::min(n, 4u);\r
- for ( unsigned int i = 0; i < m; ++i )\r
- {\r
- result.push_back( ea.pointAtAngle(angle) );\r
- angle += M_PI/2;\r
- if ( !(angle < 2*M_PI) ) angle -= 2*M_PI;\r
- }\r
- m = n / 4;\r
- for ( unsigned int i = 1; i < m; ++i )\r
- {\r
- for ( unsigned int j = 0; j < 4; ++j )\r
- result.push_back( result[j] );\r
- }\r
- m = n - 4 * m;\r
- for ( unsigned int i = 0; i < m; ++i )\r
- {\r
- result.push_back( result[i] );\r
- }\r
- if ( !result.empty() ) // n != 0\r
- result[0] = pointAtAngle(angle);\r
- return result;\r
-}\r
-\r
-D2<SBasis> SVGEllipticalArc::toSBasis() const\r
-{\r
- // the interval of parametrization has to be [0,1]\r
- Coord et = start_angle() + ( sweep_flag() ? sweep_angle() : -sweep_angle() );\r
- Linear param(start_angle(), et);\r
- Coord cos_rot_angle = std::cos(rotation_angle());\r
- Coord sin_rot_angle = std::sin(rotation_angle());\r
- // order = 4 seems to be enough to get a perfect looking elliptical arc\r
- // should it be choosen in function of the arc length anyway ?\r
- // or maybe a user settable parameter: toSBasis(unsigned int order) ?\r
- SBasis arc_x = ray(X) * cos(param,4);\r
- SBasis arc_y = ray(Y) * sin(param,4);\r
- D2<SBasis> arc;\r
- arc[0] = arc_x * cos_rot_angle - arc_y * sin_rot_angle + Linear(center(X),center(X));\r
- arc[1] = arc_x * sin_rot_angle + arc_y * cos_rot_angle + Linear(center(Y),center(Y));\r
- return arc;\r
-}\r
-\r
-\r
-bool SVGEllipticalArc::containsAngle(Coord angle) const\r
-{\r
- if ( sweep_flag() )\r
- if ( start_angle() < end_angle() )\r
- if ( !( angle < start_angle() || angle > end_angle() ) )\r
- return true;\r
- else\r
- return false;\r
- else\r
- if ( !( angle < start_angle() && angle > end_angle() ) )\r
- return true;\r
- else\r
- return false;\r
- else\r
- if ( start_angle() > end_angle() )\r
- if ( !( angle > start_angle() || angle < end_angle() ) )\r
- return true;\r
- else\r
- return false;\r
- else\r
- if ( !( angle > start_angle() && angle < end_angle() ) )\r
- return true;\r
- else\r
- return false;\r
-}\r
-\r
-\r
-double SVGEllipticalArc::valueAtAngle(Coord t, Dim2 d) const\r
-{\r
- double sin_rot_angle = std::sin(rotation_angle());\r
- double cos_rot_angle = std::cos(rotation_angle());\r
- if ( d == X )\r
- {\r
- return ray(X) * cos_rot_angle * std::cos(t) \r
- - ray(Y) * sin_rot_angle * std::sin(t) \r
- + center(X);\r
- }\r
- else if ( d == Y )\r
- {\r
- return ray(X) * sin_rot_angle * std::cos(t) \r
- + ray(Y) * cos_rot_angle * std::sin(t) \r
- + center(Y);\r
- }\r
- THROW_RANGEERROR("dimension parameter out of range");\r
-}\r
-\r
-\r
-Curve* SVGEllipticalArc::portion(double f, double t) const \r
-{\r
- if (f < 0) f = 0;\r
- if (f > 1) f = 1;\r
- if (t < 0) t = 0;\r
- if (t > 1) t = 1;\r
- if ( are_near(f, t) )\r
- {\r
- SVGEllipticalArc* arc = new SVGEllipticalArc();\r
- arc->m_center = arc->m_initial_point = arc->m_final_point = pointAt(f);\r
- arc->m_start_angle = arc->m_end_angle = m_start_angle;\r
- arc->m_rot_angle = m_rot_angle;\r
- arc->m_sweep = m_sweep;\r
- arc->m_large_arc = m_large_arc;\r
- }\r
- SVGEllipticalArc* arc = new SVGEllipticalArc( *this );\r
- arc->m_initial_point = pointAt(f);\r
- arc->m_final_point = pointAt(t);\r
- double sa = sweep_flag() ? sweep_angle() : -sweep_angle();\r
- arc->m_start_angle = m_start_angle + sa * f;\r
- if ( !(arc->m_start_angle < 2*M_PI) )\r
- arc->m_start_angle -= 2*M_PI;\r
- if ( arc->m_start_angle < 0 )\r
- arc->m_start_angle += 2*M_PI;\r
- arc->m_end_angle = m_start_angle + sa * t;\r
- if ( !(arc->m_end_angle < 2*M_PI) )\r
- arc->m_end_angle -= 2*M_PI;\r
- if ( arc->m_end_angle < 0 )\r
- arc->m_end_angle += 2*M_PI;\r
- if ( f > t ) arc->m_sweep = !sweep_flag();\r
- if ( large_arc_flag() && (arc->sweep_angle() < M_PI) )\r
- arc->m_large_arc = false;\r
- return arc;\r
-}\r
-\r
-// NOTE: doesn't work with 360 deg arcs\r
-void SVGEllipticalArc::calculate_center_and_extreme_angles()\r
-{\r
- if ( are_near(initialPoint(), finalPoint()) )\r
- {\r
- if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
- {\r
- m_start_angle = m_end_angle = 0;\r
- m_center = initialPoint();\r
- return;\r
- }\r
- else\r
- {\r
- THROW_RANGEERROR("initial and final point are the same");\r
- }\r
- }\r
- if ( are_near(ray(X), 0) && are_near(ray(Y), 0) )\r
- { // but initialPoint != finalPoint\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints: "\r
- "ray(X) == 0 && ray(Y) == 0 but initialPoint != finalPoint"\r
- );\r
- }\r
- if ( are_near(ray(Y), 0) )\r
- {\r
- Point v = initialPoint() - finalPoint();\r
- if ( are_near(L2sq(v), 4*ray(X)*ray(X)) )\r
- {\r
- double angle = std::atan2(v[Y], v[X]);\r
- if (angle < 0) angle += 2*M_PI;\r
- if ( are_near( angle, rotation_angle() ) )\r
- {\r
- m_start_angle = 0;\r
- m_end_angle = M_PI;\r
- m_center = v/2 + finalPoint();\r
- return;\r
- }\r
- angle -= M_PI;\r
- if ( angle < 0 ) angle += 2*M_PI;\r
- if ( are_near( angle, rotation_angle() ) )\r
- {\r
- m_start_angle = M_PI;\r
- m_end_angle = 0;\r
- m_center = v/2 + finalPoint();\r
- return;\r
- }\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints: "\r
- "ray(Y) == 0 "\r
- "and slope(initialPoint - finalPoint) != rotation_angle "\r
- "and != rotation_angle + PI"\r
- );\r
- }\r
- if ( L2sq(v) > 4*ray(X)*ray(X) )\r
- {\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints: "\r
- "ray(Y) == 0 and distance(initialPoint, finalPoint) > 2*ray(X)"\r
- );\r
- }\r
- else\r
- {\r
- THROW_RANGEERROR(\r
- "there is infinite ellipses that satisfy the given constraints: "\r
- "ray(Y) == 0 and distance(initialPoint, finalPoint) < 2*ray(X)"\r
- );\r
- }\r
- \r
- }\r
- \r
- if ( are_near(ray(X), 0) )\r
- {\r
- Point v = initialPoint() - finalPoint();\r
- if ( are_near(L2sq(v), 4*ray(Y)*ray(Y)) )\r
- {\r
- double angle = std::atan2(v[Y], v[X]);\r
- if (angle < 0) angle += 2*M_PI;\r
- double rot_angle = rotation_angle() + M_PI/2;\r
- if ( !(rot_angle < 2*M_PI) ) rot_angle -= 2*M_PI;\r
- if ( are_near( angle, rot_angle ) )\r
- {\r
- m_start_angle = M_PI/2;\r
- m_end_angle = 3*M_PI/2;\r
- m_center = v/2 + finalPoint();\r
- return;\r
- }\r
- angle -= M_PI;\r
- if ( angle < 0 ) angle += 2*M_PI;\r
- if ( are_near( angle, rot_angle ) )\r
- {\r
- m_start_angle = 3*M_PI/2;\r
- m_end_angle = M_PI/2;\r
- m_center = v/2 + finalPoint();\r
- return;\r
- }\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints: "\r
- "ray(X) == 0 "\r
- "and slope(initialPoint - finalPoint) != rotation_angle + PI/2 "\r
- "and != rotation_angle + (3/2)*PI"\r
- );\r
- }\r
- if ( L2sq(v) > 4*ray(Y)*ray(Y) )\r
- {\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints: "\r
- "ray(X) == 0 and distance(initialPoint, finalPoint) > 2*ray(Y)"\r
- );\r
- }\r
- else\r
- {\r
- THROW_RANGEERROR(\r
- "there is infinite ellipses that satisfy the given constraints: "\r
- "ray(X) == 0 and distance(initialPoint, finalPoint) < 2*ray(Y)"\r
- );\r
- }\r
- \r
- }\r
- \r
- double sin_rot_angle = std::sin(rotation_angle());\r
- double cos_rot_angle = std::cos(rotation_angle());\r
-\r
- Point sp = sweep_flag() ? initialPoint() : finalPoint();\r
- Point ep = sweep_flag() ? finalPoint() : initialPoint();\r
-\r
- Matrix m( ray(X) * cos_rot_angle, ray(X) * sin_rot_angle,\r
- -ray(Y) * sin_rot_angle, ray(Y) * cos_rot_angle,\r
- 0, 0 );\r
- Matrix im = m.inverse();\r
- Point sol = (ep - sp) * im;\r
- double half_sum_angle = std::atan2(-sol[X], sol[Y]);\r
- double half_diff_angle;\r
- if ( are_near(std::fabs(half_sum_angle), M_PI/2) )\r
- {\r
- double anti_sgn_hsa = (half_sum_angle > 0) ? -1 : 1;\r
- double arg = anti_sgn_hsa * sol[X] / 2;\r
- // if |arg| is a little bit > 1 acos returns nan\r
- if ( are_near(arg, 1) )\r
- half_diff_angle = 0;\r
- else if ( are_near(arg, -1) )\r
- half_diff_angle = M_PI;\r
- else\r
- {\r
- if ( !(-1 < arg && arg < 1) )\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints"\r
- );\r
- // assert( -1 < arg && arg < 1 );\r
- // if it fails \r
- // => there is no ellipse that satisfies the given constraints\r
- half_diff_angle = std::acos( arg );\r
- }\r
-\r
- half_diff_angle = M_PI/2 - half_diff_angle;\r
- }\r
- else\r
- {\r
- double arg = sol[Y] / ( 2 * std::cos(half_sum_angle) );\r
- // if |arg| is a little bit > 1 asin returns nan\r
- if ( are_near(arg, 1) ) \r
- half_diff_angle = M_PI/2;\r
- else if ( are_near(arg, -1) )\r
- half_diff_angle = -M_PI/2;\r
- else\r
- {\r
- if ( !(-1 < arg && arg < 1) )\r
- THROW_RANGEERROR(\r
- "there is no ellipse that satisfies the given constraints"\r
- );\r
- // assert( -1 < arg && arg < 1 );\r
- // if it fails \r
- // => there is no ellipse that satisfies the given constraints\r
- half_diff_angle = std::asin( arg );\r
- }\r
- }\r
-\r
- if ( ( m_large_arc && half_diff_angle > 0 ) \r
- || (!m_large_arc && half_diff_angle < 0 ) )\r
- {\r
- half_diff_angle = -half_diff_angle;\r
- }\r
- if ( half_sum_angle < 0 ) half_sum_angle += 2*M_PI;\r
- if ( half_diff_angle < 0 ) half_diff_angle += M_PI;\r
- \r
- m_start_angle = half_sum_angle - half_diff_angle;\r
- m_end_angle = half_sum_angle + half_diff_angle;\r
- // 0 <= m_start_angle, m_end_angle < 2PI\r
- if ( m_start_angle < 0 ) m_start_angle += 2*M_PI;\r
- if( !(m_end_angle < 2*M_PI) ) m_end_angle -= 2*M_PI;\r
- sol[0] = std::cos(m_start_angle);\r
- sol[1] = std::sin(m_start_angle);\r
- m_center = sp - sol * m;\r
- if ( !sweep_flag() )\r
- {\r
- double angle = m_start_angle;\r
- m_start_angle = m_end_angle;\r
- m_end_angle = angle;\r
- }\r
-}\r
-\r
-Coord SVGEllipticalArc::map_to_02PI(Coord t) const\r
-{\r
- if ( sweep_flag() )\r
- {\r
- Coord angle = start_angle() + sweep_angle() * t;\r
- if ( !(angle < 2*M_PI) )\r
- angle -= 2*M_PI;\r
- return angle;\r
- }\r
- else\r
- {\r
- Coord angle = start_angle() - sweep_angle() * t;\r
- if ( angle < 0 ) angle += 2*M_PI;\r
- return angle;\r
- }\r
-}\r
-\r
-Coord SVGEllipticalArc::map_to_01(Coord angle) const \r
-{\r
- return map_circular_arc_on_unit_interval(angle, start_angle(), \r
- end_angle(), sweep_flag());\r
-}\r
-\r
-\r
-std::vector<double> SVGEllipticalArc::\r
-allNearestPoints( Point const& p, double from, double to ) const\r
-{\r
- if ( from > to ) std::swap(from, to);\r
- if ( from < 0 || to > 1 )\r
- {\r
- THROW_RANGEERROR("[from,to] interval out of range");\r
- }\r
- std::vector<double> result;\r
- if ( ( are_near(ray(X), 0) && are_near(ray(Y), 0) ) || are_near(from, to) )\r
- {\r
- result.push_back(from);\r
- return result;\r
- }\r
- else if ( are_near(ray(X), 0) || are_near(ray(Y), 0) )\r
- {\r
- LineSegment seg(pointAt(from), pointAt(to));\r
- Point np = seg.pointAt( seg.nearestPoint(p) );\r
- if ( are_near(ray(Y), 0) )\r
- {\r
- if ( are_near(rotation_angle(), M_PI/2) \r
- || are_near(rotation_angle(), 3*M_PI/2) )\r
- {\r
- result = roots(np[Y], Y);\r
- }\r
- else\r
- {\r
- result = roots(np[X], X);\r
- }\r
- }\r
- else\r
- {\r
- if ( are_near(rotation_angle(), M_PI/2) \r
- || are_near(rotation_angle(), 3*M_PI/2) )\r
- {\r
- result = roots(np[X], X);\r
- }\r
- else\r
- {\r
- result = roots(np[Y], Y);\r
- }\r
- }\r
- return result;\r
- }\r
- else if ( are_near(ray(X), ray(Y)) )\r
- {\r
- Point r = p - center();\r
- if ( are_near(r, Point(0,0)) )\r
- {\r
- THROW_EXCEPTION("infinite nearest points");\r
- }\r
- // TODO: implement case r != 0\r
-// Point np = ray(X) * unit_vector(r);\r
-// std::vector<double> solX = roots(np[X],X);\r
-// std::vector<double> solY = roots(np[Y],Y);\r
-// double t;\r
-// if ( are_near(solX[0], solY[0]) || are_near(solX[0], solY[1]))\r
-// {\r
-// t = solX[0];\r
-// }\r
-// else\r
-// {\r
-// t = solX[1];\r
-// }\r
-// if ( !(t < from || t > to) )\r
-// {\r
-// result.push_back(t);\r
-// }\r
-// else\r
-// {\r
-// \r
-// }\r
- }\r
- \r
- // solve the equation <D(E(t),t)|E(t)-p> == 0\r
- // that provides min and max distance points \r
- // on the ellipse E wrt the point p\r
- // after the substitutions: \r
- // cos(t) = (1 - s^2) / (1 + s^2)\r
- // sin(t) = 2t / (1 + s^2)\r
- // where s = tan(t/2)\r
- // we get a 4th degree equation in s\r
- /*\r
- * ry s^4 ((-cy + py) Cos[Phi] + (cx - px) Sin[Phi]) + \r
- * ry ((cy - py) Cos[Phi] + (-cx + px) Sin[Phi]) + \r
- * 2 s^3 (rx^2 - ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi]) + \r
- * 2 s (-rx^2 + ry^2 + (-cx + px) rx Cos[Phi] + (-cy + py) rx Sin[Phi])\r
- */\r
-\r
- Point p_c = p - center();\r
- double rx2_ry2 = (ray(X) - ray(Y)) * (ray(X) + ray(Y));\r
- double cosrot = std::cos( rotation_angle() );\r
- double sinrot = std::sin( rotation_angle() );\r
- double expr1 = ray(X) * (p_c[X] * cosrot + p_c[Y] * sinrot);\r
- double coeff[5];\r
- coeff[4] = ray(Y) * ( p_c[Y] * cosrot - p_c[X] * sinrot );\r
- coeff[3] = 2 * ( rx2_ry2 + expr1 );\r
- coeff[2] = 0;\r
- coeff[1] = 2 * ( -rx2_ry2 + expr1 );\r
- coeff[0] = -coeff[4];\r
- \r
-// for ( unsigned int i = 0; i < 5; ++i )\r
-// std::cerr << "c[" << i << "] = " << coeff[i] << std::endl;\r
- \r
- std::vector<double> real_sol;\r
- // gsl_poly_complex_solve raises an error \r
- // if the leading coefficient is zero\r
- if ( are_near(coeff[4], 0) ) \r
- {\r
- real_sol.push_back(0);\r
- if ( !are_near(coeff[3], 0) )\r
- {\r
- double sq = -coeff[1] / coeff[3];\r
- if ( sq > 0 )\r
- {\r
- double s = std::sqrt(sq);\r
- real_sol.push_back(s);\r
- real_sol.push_back(-s);\r
- }\r
- }\r
- }\r
- else\r
- {\r
- double sol[8];\r
- gsl_poly_complex_workspace * w = gsl_poly_complex_workspace_alloc(5);\r
- gsl_poly_complex_solve(coeff, 5, w, sol );\r
- gsl_poly_complex_workspace_free(w);\r
- \r
- for ( unsigned int i = 0; i < 4; ++i )\r
- {\r
- if ( sol[2*i+1] == 0 ) real_sol.push_back(sol[2*i]);\r
- }\r
- }\r
- \r
- for ( unsigned int i = 0; i < real_sol.size(); ++i )\r
- {\r
- real_sol[i] = 2 * std::atan(real_sol[i]);\r
- if ( real_sol[i] < 0 ) real_sol[i] += 2*M_PI;\r
- }\r
- // when s -> Infinity then <D(E)|E-p> -> 0 iff coeff[4] == 0\r
- // so we add M_PI to the solutions being lim arctan(s) = PI when s->Infinity\r
- if ( (real_sol.size() % 2) != 0 )\r
- {\r
- real_sol.push_back(M_PI);\r
- }\r
- \r
- double mindistsq1 = std::numeric_limits<double>::max();\r
- double mindistsq2 = std::numeric_limits<double>::max();\r
- double dsq;\r
- unsigned int mi1, mi2;\r
- for ( unsigned int i = 0; i < real_sol.size(); ++i )\r
- {\r
- dsq = distanceSq(p, pointAtAngle(real_sol[i]));\r
- if ( mindistsq1 > dsq )\r
- {\r
- mindistsq2 = mindistsq1;\r
- mi2 = mi1;\r
- mindistsq1 = dsq;\r
- mi1 = i;\r
- }\r
- else if ( mindistsq2 > dsq )\r
- {\r
- mindistsq2 = dsq;\r
- mi2 = i;\r
- }\r
- }\r
- \r
- double t = map_to_01( real_sol[mi1] );\r
- if ( !(t < from || t > to) )\r
- {\r
- result.push_back(t);\r
- }\r
- \r
- bool second_sol = false; \r
- t = map_to_01( real_sol[mi2] );\r
- if ( real_sol.size() == 4 && !(t < from || t > to) )\r
- {\r
- if ( result.empty() || are_near(mindistsq1, mindistsq2) )\r
- {\r
- result.push_back(t);\r
- second_sol = true;\r
- }\r
- }\r
- \r
- // we need to test extreme points too\r
- double dsq1 = distanceSq(p, pointAt(from));\r
- double dsq2 = distanceSq(p, pointAt(to));\r
- if ( second_sol )\r
- {\r
- if ( mindistsq2 > dsq1 )\r
- {\r
- result.clear();\r
- result.push_back(from);\r
- mindistsq2 = dsq1;\r
- }\r
- else if ( are_near(mindistsq2, dsq) )\r
- {\r
- result.push_back(from);\r
- }\r
- if ( mindistsq2 > dsq2 )\r
- {\r
- result.clear();\r
- result.push_back(to);\r
- }\r
- else if ( are_near(mindistsq2, dsq2) )\r
- {\r
- result.push_back(to);\r
- }\r
- \r
- }\r
- else\r
- {\r
- if ( result.empty() )\r
- {\r
- if ( are_near(dsq1, dsq2) )\r
- {\r
- result.push_back(from);\r
- result.push_back(to);\r
- }\r
- else if ( dsq2 > dsq1 )\r
- {\r
- result.push_back(from);\r
- }\r
- else\r
- {\r
- result.push_back(to);\r
- }\r
- }\r
- }\r
- \r
- return result;\r
-}\r
-\r
-\r
-} // end namespace Geom\r
-\r
-\r
-/*\r
- Local Variables:\r
- mode:c++\r
- c-file-style:"stroustrup"\r
- c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))\r
- indent-tabs-mode:nil\r
- fill-column:99\r
- End:\r
-*/\r
-// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :\r
-\r
-\r
diff --git a/src/2geom/svg-path.cpp b/src/2geom/svg-path.cpp
index c09e5c929e408874cf933d7a8b268653dbfb98d8..2c275cb4c5d247222c487ff15f25087164086b5c 100644 (file)
--- a/src/2geom/svg-path.cpp
+++ b/src/2geom/svg-path.cpp
sink.quadTo(curve[1], curve[2]);
}
-void output(SVGEllipticalArc const &/*curve*/, SVGPathSink &/*sink*/) {
+void output(EllipticalArc const &/*curve*/, SVGPathSink &/*sink*/) {
// FIXME
THROW_NOTIMPLEMENTED();
}
output_as<LineSegment>(*iter, sink) ||
output_as<CubicBezier>(*iter, sink) ||
output_as<QuadraticBezier>(*iter, sink) ||
- output_as<SVGEllipticalArc>(*iter, sink) ||
+ output_as<EllipticalArc>(*iter, sink) ||
output_as<Curve>(*iter, sink);
}
diff --git a/src/2geom/svg-path.h b/src/2geom/svg-path.h
index 7dd2f31ac5ec76a8b368d52164e539a233b7ea35..870285b7fd0a22ced0ec85edee4fd1ac8e904c17 100644 (file)
--- a/src/2geom/svg-path.h
+++ b/src/2geom/svg-path.h
void arcTo(double rx, double ry, double angle,
bool large_arc, bool sweep, Point p)
{
- _path.template appendNew<SVGEllipticalArc>(rx, ry, angle,
+ _path.template appendNew<EllipticalArc>(rx, ry, angle,
large_arc, sweep, p);
}