GCC Code Coverage Report

Directory: ./
File: include/2geom/bezier.h
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 92 124 74.2%
Functions: 27 39 69.2%
Branches: 46 104 44.2%
Line Branch Exec Source
1 /**
2 * @file
3 * @brief Bernstein-Bezier polynomial
4 *//*
5 * Authors:
6 * MenTaLguY <mental@rydia.net>
7 * Michael Sloan <mgsloan@gmail.com>
8 * Nathan Hurst <njh@njhurst.com>
9 * Krzysztof Kosiński <tweenk.pl@gmail.com>
10 *
11 * Copyright 2007-2015 Authors
12 *
13 * This library is free software; you can redistribute it and/or
14 * modify it either under the terms of the GNU Lesser General Public
15 * License version 2.1 as published by the Free Software Foundation
16 * (the "LGPL") or, at your option, under the terms of the Mozilla
17 * Public License Version 1.1 (the "MPL"). If you do not alter this
18 * notice, a recipient may use your version of this file under either
19 * the MPL or the LGPL.
20 *
21 * You should have received a copy of the LGPL along with this library
22 * in the file COPYING-LGPL-2.1; if not, write to the Free Software
23 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
24 * You should have received a copy of the MPL along with this library
25 * in the file COPYING-MPL-1.1
26 *
27 * The contents of this file are subject to the Mozilla Public License
28 * Version 1.1 (the "License"); you may not use this file except in
29 * compliance with the License. You may obtain a copy of the License at
30 * http://www.mozilla.org/MPL/
31 *
32 * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
33 * OF ANY KIND, either express or implied. See the LGPL or the MPL for
34 * the specific language governing rights and limitations.
35 *
36 */
37
38 #ifndef LIB2GEOM_SEEN_BEZIER_H
39 #define LIB2GEOM_SEEN_BEZIER_H
40
41 #include <algorithm>
42 #include <valarray>
43 #include <2geom/coord.h>
44 #include <2geom/d2.h>
45 #include <2geom/math-utils.h>
46
47 namespace Geom {
48
49 /** @brief Compute the value of a Bernstein-Bezier polynomial.
50 * This method uses a Horner-like fast evaluation scheme.
51 * @param t Time value
52 * @param c_ Pointer to coefficients
53 * @param n Degree of the polynomial (number of coefficients minus one) */
54 template <typename T>
55 784396 inline T bernstein_value_at(double t, T const *c_, unsigned n) {
56 784396 double u = 1.0 - t;
57 784396 double bc = 1;
58 784396 double tn = 1;
59 784396 T tmp = c_[0]*u;
60
2/2
✓ Branch 0 taken 2769141 times.
✓ Branch 1 taken 784396 times.
 3553537 for(unsigned i=1; i<n; i++){
61 2769141 tn = tn*t;
62 2769141 bc = bc*(n-i+1)/i;
63 2769141 tmp = (tmp + tn*bc*c_[i])*u;
64 }
65 784396 return (tmp + tn*t*c_[n]);
66 }
67
68 /** @brief Perform Casteljau subdivision of a Bezier polynomial.
69 * Given an array of coefficients and a time value, computes two new Bernstein-Bezier basis
70 * polynomials corresponding to the \f$[0, t]\f$ and \f$[t, 1]\f$ intervals of the original one.
71 * @param t Time value
72 * @param v Array of input coordinates
73 * @param left Output polynomial corresponding to \f$[0, t]\f$
74 * @param right Output polynomial corresponding to \f$[t, 1]\f$
75 * @param order Order of the input polynomial, equal to one less the number of coefficients
76 * @return Value of the polynomial at @a t */
77 template <typename T>
78 22289 inline T casteljau_subdivision(double t, T const *v, T *left, T *right, unsigned order) {
79 // The Horner-like scheme gives very slightly different results, but we need
80 // the result of subdivision to match exactly with Bezier's valueAt function.
81 22289 T val = bernstein_value_at(t, v, order);
82
83
3/4
✓ Branch 0 taken 66 times.
✓ Branch 1 taken 22223 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 66 times.
 22289 if (!left && !right) {
84 ✗ return val;
85 }
86
87
2/2
✓ Branch 0 taken 74 times.
✓ Branch 1 taken 22215 times.
 22289 if (!right) {
88
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 74 times.
 74 if (left != v) {
89 ✗ std::copy(v, v + order + 1, left);
90 }
91
2/2
✓ Branch 0 taken 222 times.
✓ Branch 1 taken 74 times.
 296 for (std::size_t i = order; i > 0; --i) {
92
2/2
✓ Branch 0 taken 444 times.
✓ Branch 1 taken 222 times.
 666 for (std::size_t j = i; j <= order; ++j) {
93
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 444 left[j] = lerp(t, left[j-1], left[j]);
94 }
95 }
96 74 left[order] = val;
97 74 return left[order];
98 }
99
100
2/2
✓ Branch 0 taken 22149 times.
✓ Branch 1 taken 66 times.
 22215 if (right != v) {
101
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 22149 std::copy(v, v + order + 1, right);
102 }
103
2/2
✓ Branch 0 taken 174958 times.
✓ Branch 1 taken 22215 times.
 197173 for (std::size_t i = 1; i <= order; ++i) {
104
2/2
✓ Branch 0 taken 174760 times.
✓ Branch 1 taken 198 times.
 174958 if (left) {
105 174760 left[i-1] = right[0];
106 }
107
2/2
✓ Branch 0 taken 831550 times.
✓ Branch 1 taken 174958 times.
 1006508 for (std::size_t j = i; j > 0; --j) {
108
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 831550 right[j-1] = lerp(t, right[j-1], right[j]);
109 }
110 }
111 22215 right[0] = val;
112
2/2
✓ Branch 0 taken 22149 times.
✓ Branch 1 taken 66 times.
 22215 if (left) {
113 22149 left[order] = right[0];
114 }
115 22215 return right[0];
116 }
117
118 /**
119 * @brief Polynomial in Bernstein-Bezier basis
120 * @ingroup Fragments
121 */
122 class Bezier
123 : boost::arithmetic< Bezier, double
124 , boost::additive< Bezier
125 > >
126 {
127 private:
128 std::valarray<Coord> c_;
129
130 friend Bezier portion(const Bezier & a, Coord from, Coord to);
131 friend OptInterval bounds_fast(Bezier const & b);
132 friend Bezier derivative(const Bezier & a);
133 friend class Bernstein;
134
135 void
136 find_bezier_roots(std::vector<double> & solutions,
137 double l, double r) const;
138
139 protected:
140 Bezier(Coord const c[], unsigned ord)
141 : c_(c, ord+1)
142 {}
143
144 public:
145 10305879 unsigned order() const { return c_.size()-1;}
146 67457 unsigned degree() const { return order(); }
147 4129591 unsigned size() const { return c_.size();}
148
149 3631989 Bezier() {}
150 231816 Bezier(const Bezier& b) :c_(b.c_) {}
151 4880026 Bezier &operator=(Bezier const &other) {
152
2/2
✓ Branch 2 taken 2548894 times.
✓ Branch 3 taken 2331132 times.
 4880026 if ( c_.size() != other.c_.size() ) {
153 2548894 c_.resize(other.c_.size());
154 }
155 4880026 c_ = other.c_;
156 4880026 return *this;
157 }
158
159 bool operator==(Bezier const &other) const
160 {
161 if (degree() != other.degree()) {
162 return false;
163 }
164
165 for (size_t i = 0; i < c_.size(); i++) {
166 if (c_[i] != other.c_[i]) {
167 return false;
168 }
169 }
170 return true;
171 }
172
173 bool operator!=(Bezier const &other) const
174 {
175 return !(*this == other);
176 }
177
178 struct Order {
179 unsigned order;
180 490836 explicit Order(Bezier const &b) : order(b.order()) {}
181 62590 explicit Order(unsigned o) : order(o) {}
182 operator unsigned() const { return order; }
183 };
184
185 //Construct an arbitrary order bezier
186
1/2
✓ Branch 2 taken 553211 times.
✗ Branch 3 not taken.
 553211 Bezier(Order ord) : c_(0., ord.order+1) {
187
2/4
✓ Branch 1 taken 553211 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 553211 times.
 553211 assert(ord.order == order());
188 553211 }
189
190 /// @name Construct Bezier polynomials from their control points
191 /// @{
192
1/2
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 10 explicit Bezier(Coord c0) : c_(0., 1) {
193 10 c_[0] = c0;
194 10 }
195
1/2
✓ Branch 2 taken 1095026 times.
✗ Branch 3 not taken.
 1095026 Bezier(Coord c0, Coord c1) : c_(0., 2) {
196 1095026 c_[0] = c0; c_[1] = c1;
197 1095026 }
198 ✗ Bezier(Coord c0, Coord c1, Coord c2) : c_(0., 3) {
199 ✗ c_[0] = c0; c_[1] = c1; c_[2] = c2;
200 ✗ }
201
1/2
✓ Branch 2 taken 834900 times.
✗ Branch 3 not taken.
 834900 Bezier(Coord c0, Coord c1, Coord c2, Coord c3) : c_(0., 4) {
202 834900 c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3;
203 834900 }
204 Bezier(Coord c0, Coord c1, Coord c2, Coord c3, Coord c4) : c_(0., 5) {
205 c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3; c_[4] = c4;
206 }
207 Bezier(Coord c0, Coord c1, Coord c2, Coord c3, Coord c4,
208 Coord c5) : c_(0., 6) {
209 c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3; c_[4] = c4;
210 c_[5] = c5;
211 }
212 Bezier(Coord c0, Coord c1, Coord c2, Coord c3, Coord c4,
213 Coord c5, Coord c6) : c_(0., 7) {
214 c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3; c_[4] = c4;
215 c_[5] = c5; c_[6] = c6;
216 }
217 Bezier(Coord c0, Coord c1, Coord c2, Coord c3, Coord c4,
218 Coord c5, Coord c6, Coord c7) : c_(0., 8) {
219 c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3; c_[4] = c4;
220 c_[5] = c5; c_[6] = c6; c_[7] = c7;
221 }
222 Bezier(Coord c0, Coord c1, Coord c2, Coord c3, Coord c4,
223 Coord c5, Coord c6, Coord c7, Coord c8) : c_(0., 9) {
224 c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3; c_[4] = c4;
225 c_[5] = c5; c_[6] = c6; c_[7] = c7; c_[8] = c8;
226 }
227 Bezier(Coord c0, Coord c1, Coord c2, Coord c3, Coord c4,
228 Coord c5, Coord c6, Coord c7, Coord c8, Coord c9) : c_(0., 10) {
229 c_[0] = c0; c_[1] = c1; c_[2] = c2; c_[3] = c3; c_[4] = c4;
230 c_[5] = c5; c_[6] = c6; c_[7] = c7; c_[8] = c8; c_[9] = c9;
231 }
232
233 template <typename Iter>
234 ✗ Bezier(Iter first, Iter last) {
235 ✗ c_.resize(std::distance(first, last));
236 ✗ for (std::size_t i = 0; first != last; ++first, ++i) {
237 ✗ c_[i] = *first;
238 }
239 ✗ }
240 Bezier(std::vector<Coord> const &vec)
241 : c_(&vec[0], vec.size())
242 {}
243 /// @}
244
245 67361 void resize (unsigned int n, Coord v = 0) {
246 67361 c_.resize (n, v);
247 67361 }
248 67361 void clear() {
249 67361 c_.resize(0);
250 67361 }
251
252 //IMPL: FragmentConcept
253 typedef Coord output_type;
254 ✗ bool isZero(double eps=EPSILON) const {
255 ✗ for(unsigned i = 0; i <= order(); i++) {
256 ✗ if( ! are_near(c_[i], 0., eps) ) return false;
257 }
258 ✗ return true;
259 }
260 67929 bool isConstant(double eps=EPSILON) const {
261
2/2
✓ Branch 1 taken 69640 times.
✓ Branch 2 taken 241 times.
 69881 for(unsigned i = 1; i <= order(); i++) {
262
2/2
✓ Branch 3 taken 67688 times.
✓ Branch 4 taken 1952 times.
 69640 if( ! are_near(c_[i], c_[0], eps) ) return false;
263 }
264 241 return true;
265 }
266 ✗ bool isFinite() const {
267 ✗ for(unsigned i = 0; i <= order(); i++) {
268 ✗ if(!std::isfinite(c_[i])) return false;
269 }
270 ✗ return true;
271 }
272 2057324 Coord at0() const { return c_[0]; }
273 ✗ Coord &at0() { return c_[0]; }
274 1292427 Coord at1() const { return c_[order()]; }
275 ✗ Coord &at1() { return c_[order()]; }
276
277 533796 Coord valueAt(double t) const {
278 533796 return bernstein_value_at(t, &c_[0], order());
279 }
280 ✗ Coord operator()(double t) const { return valueAt(t); }
281
282 SBasis toSBasis() const;
283
284 15919432 Coord &operator[](unsigned ix) { return c_[ix]; }
285 7908246 Coord const &operator[](unsigned ix) const { return c_[ix]; }
286
287 void setCoeff(unsigned ix, double val) { c_[ix] = val; }
288
289 // The size of the returned vector equals n_derivs+1.
290 std::vector<Coord> valueAndDerivatives(Coord t, unsigned n_derivs) const;
291
292 void subdivide(Coord t, Bezier *left, Bezier *right) const;
293 std::pair<Bezier, Bezier> subdivide(Coord t) const;
294
295 std::vector<Coord> roots() const;
296 std::vector<Coord> roots(Interval const &ivl) const;
297
298 Bezier forward_difference(unsigned k) const;
299 Bezier elevate_degree() const;
300 Bezier reduce_degree() const;
301 Bezier elevate_to_degree(unsigned newDegree) const;
302 Bezier deflate() const;
303
304 // basic arithmetic operators
305 28 Bezier &operator+=(double v) {
306 28 c_ += v;
307 28 return *this;
308 }
309 ✗ Bezier &operator-=(double v) {
310 ✗ c_ -= v;
311 ✗ return *this;
312 }
313 68 Bezier &operator*=(double v) {
314 68 c_ *= v;
315 68 return *this;
316 }
317 Bezier &operator/=(double v) {
318 c_ /= v;
319 return *this;
320 }
321 Bezier &operator+=(Bezier const &other);
322 Bezier &operator-=(Bezier const &other);
323
324 /// Unary minus
325 Bezier operator-() const
326 {
327 Bezier result;
328 result.c_ = -c_;
329 return result;
330 }
331 };
332
333
334 void bezier_to_sbasis(SBasis &sb, Bezier const &bz);
335
336 Bezier operator*(Bezier const &f, Bezier const &g);
337 46 inline Bezier multiply(Bezier const &f, Bezier const &g) {
338 46 Bezier result = f * g;
339 46 return result;
340 }
341
342 437024 inline Bezier reverse(const Bezier & a) {
343
2/4
✓ Branch 2 taken 437024 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 437024 times.
✗ Branch 6 not taken.
 437024 Bezier result = Bezier(Bezier::Order(a));
344
3/4
✓ Branch 1 taken 2516260 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2079236 times.
✓ Branch 4 taken 437024 times.
 2516260 for(unsigned i = 0; i <= a.order(); i++)
345
1/2
✓ Branch 1 taken 2079236 times.
✗ Branch 2 not taken.
 2079236 result[i] = a[a.order() - i];
346 437024 return result;
347 ✗ }
348
349 Bezier portion(const Bezier & a, double from, double to);
350
351 // XXX Todo: how to handle differing orders
352 117 inline std::vector<Point> bezier_points(const D2<Bezier > & a) {
353 117 std::vector<Point> result;
354
3/4
✓ Branch 2 taken 593 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 476 times.
✓ Branch 5 taken 117 times.
 593 for(unsigned i = 0; i <= a[0].order(); i++) {
355 476 Point p;
356
2/2
✓ Branch 3 taken 952 times.
✓ Branch 4 taken 476 times.
 1428 for(unsigned d = 0; d < 2; d++) p[d] = a[d][i];
357
1/2
✓ Branch 1 taken 476 times.
✗ Branch 2 not taken.
 476 result.push_back(p);
358 }
359 117 return result;
360 ✗ }
361
362 Bezier derivative(Bezier const &a);
363 Bezier integral(Bezier const &a);
364 OptInterval bounds_fast(Bezier const &b);
365 OptInterval bounds_exact(Bezier const &b);
366 OptInterval bounds_local(Bezier const &b, OptInterval const &i);
367
368 /// Expand an interval to the image of a Bézier-Bernstein polynomial, assuming it already contains the initial point x0.
369 void bezier_expand_to_image(Interval &range, Coord x0, Coord x1, Coord x2);
370 void bezier_expand_to_image(Interval &range, Coord x0, Coord x1, Coord x2, Coord x3);
371
372 ✗ inline std::ostream &operator<< (std::ostream &os, const Bezier & b) {
373 ✗ os << "Bezier(";
374 ✗ for(unsigned i = 0; i < b.order(); i++) {
375 ✗ os << format_coord_nice(b[i]) << ", ";
376 }
377 ✗ os << format_coord_nice(b[b.order()]) << ")";
378 ✗ return os;
379 }
380
381 } // namespace Geom
382
383 #endif // LIB2GEOM_SEEN_BEZIER_H
384
385 /*
386 Local Variables:
387 mode:c++
388 c-file-style:"stroustrup"
389 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
390 indent-tabs-mode:nil
391 fill-column:99
392 End:
393 */
394 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
395