GCC Code Coverage Report

Directory: ./
File: src/2geom/bezier-curve.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 231 354 65.3%
Functions: 26 37 70.3%
Branches: 245 596 41.1%
Line Branch Exec Source
1 /* Bezier curve implementation
2 *
3 * Authors:
4 * MenTaLguY <mental@rydia.net>
5 * Marco Cecchetti <mrcekets at gmail.com>
6 * Krzysztof Kosiński <tweenk.pl@gmail.com>
7 *
8 * Copyright 2007-2009 Authors
9 *
10 * This library is free software; you can redistribute it and/or
11 * modify it either under the terms of the GNU Lesser General Public
12 * License version 2.1 as published by the Free Software Foundation
13 * (the "LGPL") or, at your option, under the terms of the Mozilla
14 * Public License Version 1.1 (the "MPL"). If you do not alter this
15 * notice, a recipient may use your version of this file under either
16 * the MPL or the LGPL.
17 *
18 * You should have received a copy of the LGPL along with this library
19 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
20 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 * You should have received a copy of the MPL along with this library
22 * in the file COPYING-MPL-1.1
23 *
24 * The contents of this file are subject to the Mozilla Public License
25 * Version 1.1 (the "License"); you may not use this file except in
26 * compliance with the License. You may obtain a copy of the License at
27 * http://www.mozilla.org/MPL/
28 *
29 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
30 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
31 * the specific language governing rights and limitations.
32 */
33
34 #include <2geom/bezier-curve.h>
35 #include <2geom/path-sink.h>
36 #include <2geom/basic-intersection.h>
37 #include <2geom/nearest-time.h>
38 #include <2geom/polynomial.h>
39
40 namespace Geom
41 {
42
43 /**
44 * @class BezierCurve
45 * @brief Two-dimensional Bezier curve of arbitrary order.
46 *
47 * Bezier curves are an expansion of the concept of linear interpolation to n points.
48 * Linear segments in 2Geom are in fact Bezier curves of order 1.
49 *
50 * Let \f$\mathbf{B}_{\mathbf{p}_0\mathbf{p}_1\ldots\mathbf{p}_n}\f$ denote a Bezier curve
51 * of order \f$n\f$ defined by the points \f$\mathbf{p}_0, \mathbf{p}_1, \ldots, \mathbf{p}_n\f$.
52 * Bezier curve of order 1 is a linear interpolation curve between two points, defined as
53 * \f[ \mathbf{B}_{\mathbf{p}_0\mathbf{p}_1}(t) = (1-t)\mathbf{p}_0 + t\mathbf{p}_1 \f]
54 * If we now substitute points \f$\mathbf{p_0}\f$ and \f$\mathbf{p_1}\f$ in this definition
55 * by linear interpolations, we get the definition of a Bezier curve of order 2, also called
56 * a quadratic Bezier curve.
57 * \f{align*}{ \mathbf{B}_{\mathbf{p}_0\mathbf{p}_1\mathbf{p}_2}(t)
58 &= (1-t) \mathbf{B}_{\mathbf{p}_0\mathbf{p}_1}(t) + t \mathbf{B}_{\mathbf{p}_1\mathbf{p}_2}(t) \\
59 \mathbf{B}_{\mathbf{p}_0\mathbf{p}_1\mathbf{p}_2}(t)
60 &= (1-t)^2\mathbf{p}_0 + 2(1-t)t\mathbf{p}_1 + t^2\mathbf{p}_2 \f}
61 * By substituting points for quadratic Bezier curves in the original definition,
62 * we get a Bezier curve of order 3, called a cubic Bezier curve.
63 * \f{align*}{ \mathbf{B}_{\mathbf{p}_0\mathbf{p}_1\mathbf{p}_2\mathbf{p}_3}(t)
64 &= (1-t) \mathbf{B}_{\mathbf{p}_0\mathbf{p}_1\mathbf{p}_2}(t)
65 + t \mathbf{B}_{\mathbf{p}_1\mathbf{p}_2\mathbf{p}_3}(t) \\
66 \mathbf{B}_{\mathbf{p}_0\mathbf{p}_1\mathbf{p}_2\mathbf{p}_3}(t)
67 &= (1-t)^3\mathbf{p}_0+3(1-t)^2t\mathbf{p}_1+3(1-t)t^2\mathbf{p}_2+t^3\mathbf{p}_3 \f}
68 * In general, a Bezier curve or order \f$n\f$ can be recursively defined as
69 * \f[ \mathbf{B}_{\mathbf{p}_0\mathbf{p}_1\ldots\mathbf{p}_n}(t)
70 = (1-t) \mathbf{B}_{\mathbf{p}_0\mathbf{p}_1\ldots\mathbf{p}_{n-1}}(t)
71 + t \mathbf{B}_{\mathbf{p}_1\mathbf{p}_2\ldots\mathbf{p}_n}(t) \f]
72 *
73 * This substitution can be repeated an arbitrary number of times. To picture this, imagine
74 * the evaluation of a point on the curve as follows: first, all control points are joined with
75 * straight lines, and a point corresponding to the selected time value is marked on them.
76 * Then, the marked points are joined with straight lines and the point corresponding to
77 * the time value is marked. This is repeated until only one marked point remains, which is the
78 * point at the selected time value.
79 *
80 * @image html bezier-curve-evaluation.png "Evaluation of the Bezier curve"
81 *
82 * An important property of the Bezier curves is that their parameters (control points)
83 * have an intuitive geometric interpretation. Because of this, they are frequently used
84 * in vector graphics editors.
85 *
86 * Every Bezier curve is contained in its control polygon (the convex polygon composed
87 * of its control points). This fact is useful for sweepline algorithms and intersection.
88 *
89 * @par Implementation notes
90 * The order of a Bezier curve is immuable once it has been created. Normally, you should
91 * know the order at compile time and use the BezierCurveN template. If you need to determine
92 * the order at runtime, use the BezierCurve::create() function. It will create a BezierCurveN
93 * for orders 1, 2 and 3 (up to cubic Beziers), so you can later <tt>dynamic_cast</tt>
94 * to those types, and for higher orders it will create an instance of BezierCurve.
95 *
96 * @relates BezierCurveN
97 * @ingroup Curves
98 */
99
100 /**
101 * @class BezierCurveN
102 * @brief Bezier curve with compile-time specified order.
103 *
104 * @tparam degree unsigned value indicating the order of the Bezier curve
105 *
106 * @relates BezierCurve
107 * @ingroup Curves
108 */
109
110
111 ✗ BezierCurve::BezierCurve(std::vector<Point> const &pts)
112 ✗ : inner(pts)
113 {
114 ✗ if (pts.size() < 2) {
115 ✗ THROW_RANGEERROR("Bezier curve must have at least 2 control points");
116 }
117 ✗ }
118
119 290 bool BezierCurve::isDegenerate() const
120 {
121
1/2
✓ Branch 0 taken 290 times.
✗ Branch 1 not taken.
 290 for (unsigned d = 0; d < 2; ++d) {
122 290 Coord ic = inner[d][0];
123
1/2
✓ Branch 1 taken 309 times.
✗ Branch 2 not taken.
 309 for (unsigned i = 1; i < size(); ++i) {
124
2/2
✓ Branch 2 taken 290 times.
✓ Branch 3 taken 19 times.
 309 if (inner[d][i] != ic) return false;
125 }
126 }
127 ✗ return true;
128 }
129
130 /** Return false if there are at least 3 distinct control points, true otherwise. */
131 133 bool BezierCurve::isLineSegment() const
132 {
133
1/2
✓ Branch 1 taken 133 times.
✗ Branch 2 not taken.
 133 auto const last_idx = size() - 1;
134
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 133 times.
 133 if (last_idx == 1) {
135 ✗ return true;
136 }
137
1/2
✓ Branch 2 taken 133 times.
✗ Branch 3 not taken.
 133 auto const start = controlPoint(0);
138
1/2
✓ Branch 2 taken 133 times.
✗ Branch 3 not taken.
 133 auto const end = controlPoint(last_idx);
139
1/2
✓ Branch 0 taken 139 times.
✗ Branch 1 not taken.
 139 for (unsigned i = 1; i < last_idx; ++i) {
140
1/2
✓ Branch 2 taken 139 times.
✗ Branch 3 not taken.
 139 auto const pi = controlPoint(i);
141
5/6
✓ Branch 1 taken 133 times.
✓ Branch 2 taken 6 times.
✓ Branch 4 taken 133 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 133 times.
✓ Branch 7 taken 6 times.
 139 if (pi != start && pi != end) {
142 133 return false;
143 }
144 }
145 ✗ return true;
146 }
147
148 ✗ void BezierCurve::expandToTransformed(Rect &bbox, Affine const &transform) const
149 {
150 ✗ bbox |= bounds_exact(inner * transform);
151 ✗ }
152
153 30010 Coord BezierCurve::length(Coord tolerance) const
154 {
155
1/5
✗ Branch 1 not taken.
✓ Branch 2 taken 30010 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 30010 switch (order())
156 {
157 ✗ case 0:
158 ✗ return 0.0;
159 30010 case 1:
160
3/6
✓ Branch 2 taken 30010 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 30010 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 30010 times.
✗ Branch 10 not taken.
 30010 return distance(initialPoint(), finalPoint());
161 ✗ case 2:
162 {
163 ✗ std::vector<Point> pts = controlPoints();
164 ✗ return bezier_length(pts[0], pts[1], pts[2], tolerance);
165 ✗ }
166 ✗ case 3:
167 {
168 ✗ std::vector<Point> pts = controlPoints();
169 ✗ return bezier_length(pts[0], pts[1], pts[2], pts[3], tolerance);
170 ✗ }
171 ✗ default:
172 ✗ return bezier_length(controlPoints(), tolerance);
173 }
174 }
175
176 std::vector<CurveIntersection>
177 59 BezierCurve::intersect(Curve const &other, Coord eps) const
178 {
179 59 std::vector<CurveIntersection> result;
180
181 // in case we encounter an order-1 curve created from a vector
182 // or a degenerate elliptical arc
183
2/4
✓ Branch 1 taken 59 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 59 times.
 59 if (isLineSegment()) {
184 ✗ LineSegment ls(initialPoint(), finalPoint());
185 ✗ result = ls.intersect(other);
186 ✗ return result;
187 ✗ }
188
189 // here we are sure that this curve is at least a quadratic Bezier
190
1/2
✓ Branch 0 taken 59 times.
✗ Branch 1 not taken.
 59 BezierCurve const *bez = dynamic_cast<BezierCurve const *>(&other);
191
2/2
✓ Branch 0 taken 58 times.
✓ Branch 1 taken 1 times.
 59 if (bez) {
192 58 std::vector<std::pair<double, double> > xs;
193
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 58 find_intersections(xs, inner, bez->inner, eps);
194
2/2
✓ Branch 7 taken 84 times.
✓ Branch 8 taken 58 times.
 142 for (auto & i : xs) {
195
1/2
✓ Branch 2 taken 84 times.
✗ Branch 3 not taken.
 84 CurveIntersection x(*this, other, i.first, i.second);
196
1/2
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
 84 result.push_back(x);
197 }
198 return result;
199 58 }
200
201 // pass other intersection types to the other curve
202
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 result = other.intersect(*this, eps);
203 1 transpose_in_place(result);
204 1 return result;
205 ✗ }
206
207 66 bool BezierCurve::isNear(Curve const &c, Coord precision) const
208 {
209
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 57 times.
 66 if (this == &c) return true;
210
211
1/2
✓ Branch 0 taken 57 times.
✗ Branch 1 not taken.
 57 BezierCurve const *other = dynamic_cast<BezierCurve const *>(&c);
212
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 55 times.
 57 if (!other) return false;
213
214
2/4
✓ Branch 5 taken 55 times.
✗ Branch 6 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 55 times.
 55 if (!are_near(inner.at0(), other->inner.at0(), precision)) return false;
215
5/8
✓ Branch 2 taken 55 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 55 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 55 times.
✗ Branch 10 not taken.
✓ Branch 13 taken 6 times.
✓ Branch 14 taken 49 times.
 55 if (!are_near(inner.at1(), other->inner.at1(), precision)) return false;
216
217
2/2
✓ Branch 2 taken 47 times.
✓ Branch 3 taken 2 times.
 49 if (size() == other->size()) {
218
2/2
✓ Branch 1 taken 23 times.
✓ Branch 2 taken 46 times.
 69 for (unsigned i = 1; i < order(); ++i) {
219
3/4
✓ Branch 5 taken 23 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 1 times.
✓ Branch 10 taken 22 times.
 23 if (!are_near(inner.point(i), other->inner.point(i), precision)) {
220 1 return false;
221 }
222 }
223 46 return true;
224 } else {
225 // Must equalize the degrees before comparing
226
2/4
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
 2 BezierCurve elevated_this, elevated_other;
227
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
 6 for (size_t dim = 0; dim < 2; dim++) {
228
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 unsigned const our_degree = inner[dim].degree();
229
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 4 unsigned const other_degree = other->inner[dim].degree();
230
231
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if (our_degree < other_degree) {
232 // Elevate our degree
233 ✗ elevated_this.inner[dim] = inner[dim].elevate_to_degree(other_degree);
234 ✗ elevated_other.inner[dim] = other->inner[dim];
235
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 } else if (our_degree > other_degree) {
236 // Elevate the other's degree
237
1/2
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 4 elevated_this.inner[dim] = inner[dim];
238
2/4
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 4 elevated_other.inner[dim] = other->inner[dim].elevate_to_degree(our_degree);
239 } else {
240 // Equal degrees: just copy
241 ✗ elevated_this.inner[dim] = inner[dim];
242 ✗ elevated_other.inner[dim] = other->inner[dim];
243 }
244 }
245
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2 times.
 2 assert(elevated_other.size() == elevated_this.size());
246
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 return elevated_this.isNear(elevated_other, precision);
247 2 }
248 }
249
250 ✗ Curve *BezierCurve::portion(Coord f, Coord t) const
251 {
252 ✗ if (f == 0.0 && t == 1.0) {
253 ✗ return duplicate();
254 }
255 ✗ if (f == 1.0 && t == 0.0) {
256 ✗ return reverse();
257 }
258 ✗ return new BezierCurve(Geom::portion(inner, f, t));
259 }
260
261 299 bool BezierCurve::_equalTo(Curve const &c) const
262 {
263
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 259 times.
 299 if (this == &c) return true;
264
1/2
✓ Branch 0 taken 259 times.
✗ Branch 1 not taken.
 259 auto other = dynamic_cast<BezierCurve const *>(&c);
265
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 259 times.
 259 if (!other) return false;
266
267
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 259 times.
 259 if (size() != other->size()) return false;
268
269
2/2
✓ Branch 1 taken 680 times.
✓ Branch 2 taken 259 times.
 939 for (unsigned i = 0; i < size(); ++i) {
270
3/6
✓ Branch 2 taken 680 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 680 times.
✗ Branch 7 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 680 times.
 680 if (controlPoint(i) != other->controlPoint(i)) return false;
271 }
272
273 259 return true;
274 }
275
276 23 Coord BezierCurve::nearestTime(Point const &p, Coord from, Coord to) const
277 {
278 23 return nearest_time(p, inner, from, to);
279 }
280
281 ✗ void BezierCurve::feed(PathSink &sink, bool moveto_initial) const
282 {
283 ✗ if (size() > 4) {
284 ✗ Curve::feed(sink, moveto_initial);
285 ✗ return;
286 }
287
288 ✗ Point ip = controlPoint(0);
289 ✗ if (moveto_initial) {
290 ✗ sink.moveTo(ip);
291 }
292 ✗ switch (size()) {
293 ✗ case 2:
294 ✗ sink.lineTo(controlPoint(1));
295 ✗ break;
296 ✗ case 3:
297 ✗ sink.quadTo(controlPoint(1), controlPoint(2));
298 ✗ break;
299 ✗ case 4:
300 ✗ sink.curveTo(controlPoint(1), controlPoint(2), controlPoint(3));
301 ✗ break;
302 ✗ default:
303 // TODO: add a path sink method that accepts a vector of control points
304 // and converts to cubic spline by default
305 ✗ assert(false);
306 break;
307 }
308 }
309
310 ✗ BezierCurve *BezierCurve::create(std::vector<Point> const &pts)
311 {
312 ✗ switch (pts.size()) {
313 ✗ case 0:
314 case 1:
315 ✗ THROW_LOGICALERROR("BezierCurve::create: too few points in vector");
316 return NULL;
317 ✗ case 2:
318 ✗ return new LineSegment(pts[0], pts[1]);
319 ✗ case 3:
320 ✗ return new QuadraticBezier(pts[0], pts[1], pts[2]);
321 ✗ case 4:
322 ✗ return new CubicBezier(pts[0], pts[1], pts[2], pts[3]);
323 ✗ default:
324 ✗ return new BezierCurve(pts);
325 }
326 }
327
328 /**
329 * Computes the times where the radius of curvature of the bezier curve equals the given radius.
330 */
331 11 std::vector<Coord> BezierCurve::timesWithRadiusOfCurvature(double radius) const
332 {
333 /**
334 * The algorithm works as follows:
335 * Find the solutions of curvature of curve at t = curvatureValue
336 * This is equivalent to (dx*ddy-ddx*dy)/curvatureValue = (dx**2+dy**2)**(3/2)
337 * When the left side is positive, taking the square gives
338 * ((dx*ddy-ddx*dy)/curvatureValue)**2 - (dx**2+dy**2)**3 = 0
339 * This is a polyomial for BezierCurves and can be solved with root finding algos.
340 */
341
342
3/4
✓ Branch 1 taken 11 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 9 times.
 11 if (this->size() <= 2) {
343 2 return {};
344 }
345
346 9 std::vector<Coord> res;
347
348
1/2
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
 9 auto const dx = Geom::derivative(inner[X]);
349
1/2
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
 9 auto const dy = Geom::derivative(inner[Y]);
350
1/2
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
 9 auto const ddx = Geom::derivative(dx);
351
1/2
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
 9 auto const ddy = Geom::derivative(dy);
352
4/8
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 9 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 9 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 9 times.
✗ Branch 15 not taken.
 9 auto const c0 = (dx * ddy - ddx * dy) * radius;
353
3/6
✓ Branch 3 taken 9 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 9 times.
✗ Branch 11 not taken.
 9 auto const c1 = dx * dx + dy * dy;
354
4/8
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 9 times.
✗ Branch 8 not taken.
✓ Branch 11 taken 9 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 9 times.
✗ Branch 15 not taken.
 9 auto const p = c0 * c0 - c1 * c1 * c1;
355
1/2
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
 9 auto const candidates = p.roots();
356 // check which candidates have positive nominator
357 // as squaring also give negative (spurious) results
358
2/2
✓ Branch 8 taken 40 times.
✓ Branch 9 taken 9 times.
 49 for (Coord const candidate : candidates) {
359
3/4
✓ Branch 1 taken 40 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 32 times.
 40 if (c0.valueAt(candidate) > 0) {
360
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 res.push_back(candidate);
361 }
362 }
363
364 9 return res;
365 9 }
366
367 // optimized specializations for LineSegment
368
369 template <>
370 85 Curve *BezierCurveN<1>::derivative() const {
371 85 double dx = inner[X][1] - inner[X][0], dy = inner[Y][1] - inner[Y][0];
372
1/4
✓ Branch 6 taken 85 times.
✗ Branch 7 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 85 return new BezierCurveN<1>(Point(dx,dy),Point(dx,dy));
373 }
374
375 template<>
376 67 Coord BezierCurveN<1>::nearestTime(Point const& p, Coord from, Coord to) const
377 {
378 using std::swap;
379
380
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 67 times.
 67 if ( from > to ) swap(from, to);
381
1/2
✓ Branch 2 taken 67 times.
✗ Branch 3 not taken.
 67 Point ip = pointAt(from);
382
1/2
✓ Branch 2 taken 67 times.
✗ Branch 3 not taken.
 67 Point fp = pointAt(to);
383 67 Point v = fp - ip;
384 67 Coord l2v = L2sq(v);
385
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 67 times.
 67 if (l2v == 0) return 0;
386 67 Coord t = dot( p - ip, v ) / l2v;
387
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 62 times.
 67 if ( t <= 0 ) return from;
388
2/2
✓ Branch 0 taken 36 times.
✓ Branch 1 taken 26 times.
 62 else if ( t >= 1 ) return to;
389 26 else return from + t*(to-from);
390 }
391
392 /* Specialized intersection routine for line segments.
393 *
394 * NOTE: if the segments overlap in part or in full, the function returns the start and end
395 * of the overlapping subsegment as intersections. This behavior is more useful than throwing
396 * Geom::InfinitelyManySolutions.
397 */
398 template <>
399 380108 std::vector<CurveIntersection> BezierCurveN<1>::intersect(Curve const &other, Coord eps) const
400 {
401 380108 std::vector<CurveIntersection> result;
402
403 // only handle intersections with other LineSegments here
404
3/4
✓ Branch 1 taken 380108 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 380100 times.
 380108 if (!other.isLineSegment()) {
405 // pass all other types to the other curve
406
1/2
✓ Branch 2 taken 8 times.
✗ Branch 3 not taken.
 8 result = other.intersect(*this, eps);
407 8 transpose_in_place(result);
408 8 return result;
409 }
410
411
2/4
✓ Branch 3 taken 380100 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 380100 times.
✗ Branch 8 not taken.
 380100 Point const u = finalPoint() - initialPoint();
412
2/4
✓ Branch 3 taken 380100 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 380100 times.
✗ Branch 8 not taken.
 380100 Point const v = other.initialPoint() - other.finalPoint();
413
3/6
✓ Branch 1 taken 380100 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 380100 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 380100 times.
 380100 if (u.isZero() || v.isZero()) {
414 ✗ return {};
415 }
416 380100 Coord const uv = u.length() * v.length();
417 380100 Coord const u_cross_v = cross(u, v);
418 380100 bool const segments_are_parallel = std::abs(u_cross_v) < eps * uv;
419
420
2/2
✓ Branch 0 taken 190023 times.
✓ Branch 1 taken 190077 times.
 380100 if (segments_are_parallel) {
421 // We check if the segments lie on the same line.
422
1/2
✓ Branch 1 taken 190023 times.
✗ Branch 2 not taken.
 190023 Coord const distance_between_lines = std::abs(cross(u.normalized(),
423
2/4
✓ Branch 3 taken 190023 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 190023 times.
✗ Branch 8 not taken.
 380046 other.initialPoint() - initialPoint()));
424
2/2
✓ Branch 0 taken 69847 times.
✓ Branch 1 taken 120176 times.
 190023 if (distance_between_lines > eps) {
425 // Segments are parallel but aren't part of the same line => no intersections.
426 69847 return {};
427 }
428 // The segments are on the same line, so they may overlap in part or in full.
429 // Look for the times on this segment's line at which the line passes through
430 // the initial and final points of the other segment.
431 120176 Coord const ulen_inverse = 1.0 / u.length();
432 120176 auto const time_of_passage = [&](Point const &point_on_line) -> Coord {
433
2/4
✓ Branch 3 taken 240352 times.
✗ Branch 4 not taken.
✓ Branch 8 taken 240352 times.
✗ Branch 9 not taken.
 240352 return dot(u.normalized(), point_on_line - initialPoint()) * ulen_inverse;
434 120176 };
435 // Find the range of times on our segment where we travel through the other segment.
436
1/2
✓ Branch 2 taken 120176 times.
✗ Branch 3 not taken.
 360528 auto time_in_other = Interval(time_of_passage(other.initialPoint()),
437
4/8
✓ Branch 2 taken 120176 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 120176 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 120176 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 120176 times.
✗ Branch 13 not taken.
 360528 time_of_passage(other.finalPoint()));
438 120176 Coord const eps_utime = eps * ulen_inverse;
439
5/6
✓ Branch 1 taken 120176 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✓ Branch 5 taken 120167 times.
✓ Branch 6 taken 9 times.
✓ Branch 7 taken 120167 times.
 120176 if (time_in_other.min() > 1 + eps_utime || time_in_other.max() < -eps_utime) {
440 9 return {};
441 }
442
443 // Create two intersections, one at each end of the overlap interval.
444 120167 Coord last_time = infinity();
445
2/2
✓ Branch 7 taken 240334 times.
✓ Branch 8 taken 120167 times.
 360501 for (Coord t : {time_in_other.min(), time_in_other.max()}) {
446 240334 t = std::clamp(t, 0.0, 1.0);
447
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 240322 times.
 240334 if (t == last_time) {
448 12 continue;
449 }
450 240322 last_time = t;
451
1/2
✓ Branch 2 taken 240322 times.
✗ Branch 3 not taken.
 240322 auto const point = pointAt(t);
452
2/4
✓ Branch 7 taken 240322 times.
✗ Branch 8 not taken.
✓ Branch 12 taken 240322 times.
✗ Branch 13 not taken.
 240322 Coord const other_t = std::clamp(dot(v.normalized(), other.initialPoint() - point) / v.length(), 0.0, 1.0);
453
1/2
✓ Branch 2 taken 240322 times.
✗ Branch 3 not taken.
 240322 auto const other_pt = other.pointAt(other_t);
454
2/4
✓ Branch 1 taken 240322 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 240322 times.
 240322 if (distance(point, other_pt) > eps) {
455 ✗ continue;
456 }
457
2/4
✓ Branch 2 taken 240322 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 240322 times.
✗ Branch 6 not taken.
 240322 result.emplace_back(t, other_t, middle_point(point, other_pt));
458 }
459 120167 return result;
460 } else {
461 // Segments are not collinear - there is 0 or 1 intersection.
462 // This is the typical case so it's important for it to be fast.
463
2/4
✓ Branch 3 taken 190077 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 190077 times.
✗ Branch 8 not taken.
 190077 Point const w = other.initialPoint() - initialPoint();
464 190077 Coord candidate_time_this = cross(w, v) / u_cross_v;
465
466 /// Filter out candidate times outside of the interval [0, 1] fuzzed so as to accommodate
467 /// for epsilon.
468 190077 auto const is_outside_seg = [eps](Coord t, Coord len) -> bool {
469 380149 Coord const arclen_coord = t * len;
470
4/4
✓ Branch 0 taken 380147 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 54828 times.
✓ Branch 3 taken 325319 times.
 380149 return arclen_coord < -eps || arclen_coord > len + eps;
471 190077 };
472
473
2/2
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 190072 times.
 190077 if (is_outside_seg(candidate_time_this, u.length())) {
474 5 return {};
475 }
476
477 190072 Coord candidate_time_othr = cross(u, w) / u_cross_v;
478
2/2
✓ Branch 2 taken 54825 times.
✓ Branch 3 taken 135247 times.
 190072 if (is_outside_seg(candidate_time_othr, v.length())) {
479 54825 return {};
480 }
481
482 // Even if there was some fuzz in the interval, we must now restrict to [0, 1] and verify
483 // that the intersection found is still within epsilon precision (the fuzz may have been too
484 // large, depending on the angle between the intervals).
485 135247 candidate_time_this = std::clamp(candidate_time_this, 0.0, 1.0);
486 135247 candidate_time_othr = std::clamp(candidate_time_othr, 0.0, 1.0);
487
488
1/2
✓ Branch 2 taken 135247 times.
✗ Branch 3 not taken.
 135247 auto const pt_this = pointAt(candidate_time_this);
489
1/2
✓ Branch 2 taken 135247 times.
✗ Branch 3 not taken.
 135247 auto const pt_othr = other.pointAt(candidate_time_othr);
490
2/4
✓ Branch 1 taken 135247 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 135247 times.
 135247 if (distance(pt_this, pt_othr) > eps) {
491 ✗ return {};
492 }
493
2/4
✓ Branch 3 taken 135247 times.
✗ Branch 4 not taken.
✓ Branch 8 taken 135247 times.
✗ Branch 9 not taken.
 405741 return { CurveIntersection(candidate_time_this, candidate_time_othr, middle_point(pt_this, pt_othr)) };
494 }
495 380108 }
496
497 /** @brief Find intersections of a low-degree Bézier curve with a line segment.
498 *
499 * Uses algebraic solutions to low-degree polynomial equations which may be faster
500 * and more precise than iterative methods.
501 *
502 * @tparam degree The degree of the Bézier curve; must be 2 or 3.
503 * @param curve A Bézier curve of the given degree.
504 * @param line A line (but really a segment).
505 * @return Intersections between the passed curve and the fundamental segment of the line
506 * (the segment where the time parameter lies in the unit interval).
507 */
508 template <unsigned degree>
509 40048 static std::vector<CurveIntersection> bezier_line_intersections(BezierCurveN<degree> const &curve, Line const &line)
510 {
511 static_assert(degree == 2 || degree == 3, "bezier_line_intersections<degree>() error: degree must be 2 or 3.");
512
513
1/2
✓ Branch 5 taken 20024 times.
✗ Branch 6 not taken.
 40048 auto const length = distance(line.initialPoint(), line.finalPoint());
514
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 20024 times.
 40048 if (length == 0) {
515 ✗ return {};
516 }
517 40048 std::vector<CurveIntersection> result;
518
519 // Find the isometry mapping the line to the x-axis, taking the initial point to the origin
520 // and the final point to (length, 0). Apply this transformation to the Bézier curve and
521 // extract the y-coordinate polynomial.
522
1/2
✓ Branch 2 taken 20024 times.
✗ Branch 3 not taken.
 40048 auto const transform = line.rotationToZero(Y);
523
2/4
✓ Branch 4 taken 20024 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 20024 times.
✗ Branch 9 not taken.
 40048 Bezier const y = (curve.fragment() * transform)[Y];
524 40048 std::vector<double> roots;
525
526 // Find roots of the polynomial y.
527 {
528 40048 double const c2 = y[0] + y[2] - 2.0 * y[1];
529 40048 double const c1 = y[1] - y[0];
530 40048 double const c0 = y[0];
531
532 if constexpr (degree == 2) {
533
1/2
✓ Branch 2 taken 10007 times.
✗ Branch 3 not taken.
 20014 roots = solve_quadratic(c2, 2.0 * c1, c0);
534 } else if constexpr (degree == 3) {
535 20034 double const c3 = y[3] - y[0] + 3.0 * (y[1] - y[2]);
536
1/2
✓ Branch 2 taken 10017 times.
✗ Branch 3 not taken.
 20034 roots = solve_cubic(c3, 3.0 * c2, 3.0 * c1 , c0);
537 }
538 }
539
540 // Filter the roots and assemble intersections.
541
2/2
✓ Branch 8 taken 36051 times.
✓ Branch 9 taken 20024 times.
 112150 for (double root : roots) {
542
4/4
✓ Branch 0 taken 28050 times.
✓ Branch 1 taken 8001 times.
✓ Branch 2 taken 8007 times.
✓ Branch 3 taken 20043 times.
 72102 if (root < 0.0 || root > 1.0) {
543 32016 continue;
544 }
545
2/4
✓ Branch 3 taken 20043 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 20043 times.
✗ Branch 7 not taken.
 40086 Coord x = (curve.pointAt(root) * transform)[X];
546
3/4
✓ Branch 0 taken 20043 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 20035 times.
 40086 if (x < 0.0 || x > length) {
547 16 continue;
548 }
549
1/2
✓ Branch 2 taken 20035 times.
✗ Branch 3 not taken.
 40070 result.emplace_back(curve, line, root, x / length);
550 }
551 40048 return result;
552 40048 }
553
554 /* Specialized intersection routine for quadratic Bézier curves. */
555 template <>
556 10009 std::vector<CurveIntersection> BezierCurveN<2>::intersect(Curve const &other, Coord eps) const
557 {
558
2/4
✓ Branch 0 taken 10009 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10009 times.
✗ Branch 3 not taken.
 10009 if (auto other_bezier = dynamic_cast<BezierCurve const *>(&other)) {
559 10009 auto const other_degree = other_bezier->order();
560
2/2
✓ Branch 0 taken 10007 times.
✓ Branch 1 taken 2 times.
 10009 if (other_degree == 1) {
561 // Use the exact method to intersect a quadratic Bézier with a line segment.
562
2/4
✓ Branch 3 taken 10007 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 10007 times.
✗ Branch 8 not taken.
 10007 auto line = Line(other_bezier->initialPoint(), other_bezier->finalPoint());
563
1/2
✓ Branch 1 taken 10007 times.
✗ Branch 2 not taken.
 10007 return bezier_line_intersections<2>(*this, line);
564 }
565 // TODO: implement exact intersection of two quadratic Béziers using the method of resultants.
566 }
567 2 return BezierCurve::intersect(other, eps);
568 }
569
570 /* Specialized intersection routine for cubic Bézier curves. */
571 template <>
572 10070 std::vector<CurveIntersection> BezierCurveN<3>::intersect(Curve const &other, Coord eps) const
573 {
574
3/4
✓ Branch 0 taken 10070 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10069 times.
✓ Branch 3 taken 1 times.
 10070 if (auto other_bezier = dynamic_cast<BezierCurve const *>(&other)) {
575
2/2
✓ Branch 1 taken 10017 times.
✓ Branch 2 taken 52 times.
 10069 if (other_bezier->order() == 1) {
576 // Use the exact method to intersect a cubic Bézier with a line segment.
577
2/4
✓ Branch 3 taken 10017 times.
✗ Branch 4 not taken.
✓ Branch 7 taken 10017 times.
✗ Branch 8 not taken.
 10017 auto line = Line(other_bezier->initialPoint(), other_bezier->finalPoint());
578
1/2
✓ Branch 1 taken 10017 times.
✗ Branch 2 not taken.
 10017 return bezier_line_intersections<3>(*this, line);
579 }
580 }
581 53 return BezierCurve::intersect(other, eps);
582 }
583
584 template <>
585 61331 int BezierCurveN<1>::winding(Point const &p) const
586 {
587
1/2
✓ Branch 4 taken 61331 times.
✗ Branch 5 not taken.
 61331 Point ip = inner.at0(), fp = inner.at1();
588
1/2
✗ Branch 4 not taken.
✓ Branch 5 taken 61331 times.
 61331 if (p[Y] == std::max(ip[Y], fp[Y])) return 0;
589
590 61331 Point v = fp - ip;
591
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 61331 times.
 61331 assert(v[Y] != 0);
592 61331 Coord t = (p[Y] - ip[Y]) / v[Y];
593
2/4
✓ Branch 0 taken 61331 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 61331 times.
✗ Branch 3 not taken.
 61331 assert(t >= 0 && t <= 1);
594 61331 Coord xcross = lerp(t, ip[X], fp[X]);
595
2/2
✓ Branch 1 taken 30682 times.
✓ Branch 2 taken 30649 times.
 61331 if (xcross > p[X]) {
596
2/2
✓ Branch 1 taken 17503 times.
✓ Branch 2 taken 13179 times.
 30682 return v[Y] > 0 ? 1 : -1;
597 }
598 30649 return 0;
599 }
600
601 template <>
602 1696 void BezierCurveN<1>::feed(PathSink &sink, bool moveto_initial) const
603 {
604
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1696 times.
 1696 if (moveto_initial) {
605 ✗ sink.moveTo(controlPoint(0));
606 }
607
2/4
✓ Branch 2 taken 1696 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1696 times.
✗ Branch 6 not taken.
 1696 sink.lineTo(controlPoint(1));
608 1696 }
609
610 template <>
611 65 void BezierCurveN<2>::feed(PathSink &sink, bool moveto_initial) const
612 {
613
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 65 times.
 65 if (moveto_initial) {
614 ✗ sink.moveTo(controlPoint(0));
615 }
616
3/6
✓ Branch 2 taken 65 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 65 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 65 times.
✗ Branch 10 not taken.
 65 sink.quadTo(controlPoint(1), controlPoint(2));
617 65 }
618
619 template <>
620 23200 void BezierCurveN<3>::feed(PathSink &sink, bool moveto_initial) const
621 {
622
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 23200 times.
 23200 if (moveto_initial) {
623 ✗ sink.moveTo(controlPoint(0));
624 }
625
4/8
✓ Branch 2 taken 23200 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 23200 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 23200 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 23200 times.
✗ Branch 14 not taken.
 23200 sink.curveTo(controlPoint(1), controlPoint(2), controlPoint(3));
626 23200 }
627
628 50 static void bezier_expand_to_image(Rect &range, Point const &x0, Point const &x1, Point const &x2)
629 {
630
2/2
✓ Branch 4 taken 100 times.
✓ Branch 5 taken 50 times.
 150 for (auto i : { X, Y }) {
631
1/2
✓ Branch 5 taken 100 times.
✗ Branch 6 not taken.
 100 bezier_expand_to_image(range[i], x0[i], x1[i], x2[i]);
632 }
633 50 }
634
635 50 static void bezier_expand_to_image(Rect &range, Point const &x0, Point const &x1, Point const &x2, Point const &x3)
636 {
637
2/2
✓ Branch 4 taken 100 times.
✓ Branch 5 taken 50 times.
 150 for (auto i : { X, Y }) {
638
1/2
✓ Branch 6 taken 100 times.
✗ Branch 7 not taken.
 100 bezier_expand_to_image(range[i], x0[i], x1[i], x2[i], x3[i]);
639 }
640 50 }
641
642 template <>
643 50 void BezierCurveN<1>::expandToTransformed(Rect &bbox, Affine const &transform) const
644 {
645
2/4
✓ Branch 3 taken 50 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 50 times.
✗ Branch 7 not taken.
 50 bbox.expandTo(finalPoint() * transform);
646 50 }
647
648 template <>
649 50 void BezierCurveN<2>::expandToTransformed(Rect &bbox, Affine const &transform) const
650 {
651
3/6
✓ Branch 5 taken 50 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 50 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 50 times.
✗ Branch 12 not taken.
 150 bezier_expand_to_image(bbox, controlPoint(0) * transform,
652
2/4
✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 50 times.
✗ Branch 6 not taken.
 100 controlPoint(1) * transform,
653
2/4
✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 50 times.
✗ Branch 6 not taken.
 100 controlPoint(2) * transform);
654 50 }
655
656 template <>
657 50 void BezierCurveN<3>::expandToTransformed(Rect &bbox, Affine const &transform) const
658 {
659
3/6
✓ Branch 6 taken 50 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 50 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 50 times.
✗ Branch 13 not taken.
 200 bezier_expand_to_image(bbox, controlPoint(0) * transform,
660
2/4
✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 50 times.
✗ Branch 6 not taken.
 100 controlPoint(1) * transform,
661
2/4
✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 50 times.
✗ Branch 6 not taken.
 100 controlPoint(2) * transform,
662
2/4
✓ Branch 2 taken 50 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 50 times.
✗ Branch 6 not taken.
 100 controlPoint(3) * transform);
663 50 }
664
665 ✗ static Coord bezier_length_internal(std::vector<Point> &v1, Coord tolerance, int level)
666 {
667 /* The Bezier length algorithm used in 2Geom utilizes a simple fact:
668 * the Bezier curve is longer than the distance between its endpoints
669 * but shorter than the length of the polyline formed by its control
670 * points. When the difference between the two values is smaller than the
671 * error tolerance, we can be sure that the true value is no further than
672 * 0.5 * tolerance from their arithmetic mean. When it's larger, we recursively
673 * subdivide the Bezier curve into two parts and add their lengths.
674 *
675 * We cap the maximum number of subdivisions at 256, which corresponds to 8 levels.
676 */
677 ✗ Coord lower = distance(v1.front(), v1.back());
678 ✗ Coord upper = 0.0;
679 ✗ for (size_t i = 0; i < v1.size() - 1; ++i) {
680 ✗ upper += distance(v1[i], v1[i+1]);
681 }
682 ✗ if (upper - lower <= 2*tolerance || level >= 8) {
683 ✗ return (lower + upper) / 2;
684 }
685
686
687 ✗ std::vector<Point> v2 = v1;
688
689 /* Compute the right subdivision directly in v1 and the left one in v2.
690 * Explanation of the algorithm used:
691 * We have to compute the left and right edges of this triangle in which
692 * the top row are the control points of the Bezier curve, and each cell
693 * is equal to the arithmetic mean of the cells directly above it
694 * to the right and left. This corresponds to subdividing the Bezier curve
695 * at time value 0.5: the left edge has the control points of the first
696 * portion of the Bezier curve and the right edge - the second one.
697 * In the example we subdivide a curve with 5 control points (order 4).
698 *
699 * Start:
700 * 0 1 2 3 4
701 * ? ? ? ?
702 * ? ? ?
703 * ? ?
704 * ?
705 * # means we have overwritten the value, ? means we don't know
706 * the value yet. Numbers mean the value is at i-th position in the vector.
707 *
708 * After loop with i==1
709 * # 1 2 3 4
710 * 0 ? ? ? -> write 0 to v2[1]
711 * ? ? ?
712 * ? ?
713 * ?
714 *
715 * After loop with i==2
716 * # # 2 3 4
717 * # 1 ? ?
718 * 0 ? ? -> write 0 to v2[2]
719 * ? ?
720 * ?
721 *
722 * After loop with i==3
723 * # # # 3 4
724 * # # 2 ?
725 * # 1 ?
726 * 0 ? -> write 0 to v2[3]
727 * ?
728 *
729 * After loop with i==4, we have the right edge of the triangle in v1,
730 * and we write the last value needed for the left edge in v2[4].
731 */
732
733 ✗ for (size_t i = 1; i < v1.size(); ++i) {
734 ✗ for (size_t j = i; j > 0; --j) {
735 ✗ v1[j-1] = 0.5 * (v1[j-1] + v1[j]);
736 }
737 ✗ v2[i] = v1[0];
738 }
739
740 ✗ return bezier_length_internal(v1, 0.5 * tolerance, level + 1) +
741 ✗ bezier_length_internal(v2, 0.5 * tolerance, level + 1);
742 ✗ }
743
744 /** @brief Compute the length of a bezier curve given by a vector of its control points
745 * @relatesalso BezierCurve */
746 ✗ Coord bezier_length(std::vector<Point> const &points, Coord tolerance)
747 {
748 ✗ if (points.size() < 2) return 0.0;
749 ✗ std::vector<Point> v1 = points;
750 ✗ return bezier_length_internal(v1, tolerance, 0);
751 ✗ }
752
753 ✗ static Coord bezier_length_internal(Point a0, Point a1, Point a2, Coord tolerance, int level)
754 {
755 ✗ Coord lower = distance(a0, a2);
756 ✗ Coord upper = distance(a0, a1) + distance(a1, a2);
757
758 ✗ if (upper - lower <= 2*tolerance || level >= 8) {
759 ✗ return (lower + upper) / 2;
760 }
761
762 Point // Casteljau subdivision
763 // b0 = a0,
764 // c0 = a2,
765 ✗ b1 = 0.5*(a0 + a1),
766 ✗ c1 = 0.5*(a1 + a2),
767 ✗ b2 = 0.5*(b1 + c1); // == c2
768 ✗ return bezier_length_internal(a0, b1, b2, 0.5 * tolerance, level + 1) +
769 ✗ bezier_length_internal(b2, c1, a2, 0.5 * tolerance, level + 1);
770 }
771
772 /** @brief Compute the length of a quadratic bezier curve given by its control points
773 * @relatesalso QuadraticBezier */
774 ✗ Coord bezier_length(Point a0, Point a1, Point a2, Coord tolerance)
775 {
776 ✗ return bezier_length_internal(a0, a1, a2, tolerance, 0);
777 }
778
779 ✗ static Coord bezier_length_internal(Point a0, Point a1, Point a2, Point a3, Coord tolerance, int level)
780 {
781 ✗ Coord lower = distance(a0, a3);
782 ✗ Coord upper = distance(a0, a1) + distance(a1, a2) + distance(a2, a3);
783
784 ✗ if (upper - lower <= 2*tolerance || level >= 8) {
785 ✗ return (lower + upper) / 2;
786 }
787
788 Point // Casteljau subdivision
789 // b0 = a0,
790 // c0 = a3,
791 ✗ b1 = 0.5*(a0 + a1),
792 ✗ t0 = 0.5*(a1 + a2),
793 ✗ c1 = 0.5*(a2 + a3),
794 ✗ b2 = 0.5*(b1 + t0),
795 ✗ c2 = 0.5*(t0 + c1),
796 ✗ b3 = 0.5*(b2 + c2); // == c3
797 ✗ return bezier_length_internal(a0, b1, b2, b3, 0.5 * tolerance, level + 1) +
798 ✗ bezier_length_internal(b3, c2, c1, a3, 0.5 * tolerance, level + 1);
799 }
800
801 /** @brief Compute the length of a cubic bezier curve given by its control points
802 * @relatesalso CubicBezier */
803 ✗ Coord bezier_length(Point a0, Point a1, Point a2, Point a3, Coord tolerance)
804 {
805 ✗ return bezier_length_internal(a0, a1, a2, a3, tolerance, 0);
806 }
807
808 } // end namespace Geom
809
810 /*
811 Local Variables:
812 mode:c++
813 c-file-style:"stroustrup"
814 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
815 indent-tabs-mode:nil
816 fill-column:99
817 End:
818 */
819 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
820