GCC Code Coverage Report

Directory: ./
File: src/2geom/bezier-clipping.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 177 421 42.0%
Functions: 15 31 48.4%
Branches: 166 530 31.3%
Line Branch Exec Source
1 /*
2 * Implement the Bezier clipping algorithm for finding
3 * Bezier curve intersection points and collinear normals
4 *
5 * Authors:
6 * Marco Cecchetti <mrcekets at gmail.com>
7 *
8 * Copyright 2008 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
35
36
37 #include <2geom/basic-intersection.h>
38 #include <2geom/choose.h>
39 #include <2geom/point.h>
40 #include <2geom/interval.h>
41 #include <2geom/bezier.h>
42 #include <2geom/numeric/matrix.h>
43 #include <2geom/convex-hull.h>
44 #include <2geom/line.h>
45
46 #include <cassert>
47 #include <vector>
48 #include <algorithm>
49 #include <utility>
50 //#include <iomanip>
51
52 using std::swap;
53
54
55 #define VERBOSE 0
56 #define CHECK 0
57
58 namespace Geom {
59
60 namespace detail { namespace bezier_clipping {
61
62 ////////////////////////////////////////////////////////////////////////////////
63 // for debugging
64 //
65
66 ✗ void print(std::vector<Point> const& cp, const char* msg = "")
67 {
68 ✗ std::cerr << msg << std::endl;
69 ✗ for (size_t i = 0; i < cp.size(); ++i)
70 ✗ std::cerr << i << " : " << cp[i] << std::endl;
71 ✗ }
72
73 template< class charT >
74 std::basic_ostream<charT> &
75 operator<< (std::basic_ostream<charT> & os, const Interval & I)
76 {
77 os << "[" << I.min() << ", " << I.max() << "]";
78 return os;
79 }
80
81 ✗ double angle (std::vector<Point> const& A)
82 {
83 ✗ size_t n = A.size() -1;
84 ✗ double a = std::atan2(A[n][Y] - A[0][Y], A[n][X] - A[0][X]);
85 ✗ return (180 * a / M_PI);
86 }
87
88 ✗ size_t get_precision(Interval const& I)
89 {
90 ✗ double d = I.extent();
91 ✗ double e = 0.1, p = 10;
92 ✗ int n = 0;
93 ✗ while (n < 16 && d < e)
94 {
95 ✗ p *= 10;
96 ✗ e = 1/p;
97 ✗ ++n;
98 }
99 ✗ return n;
100 }
101
102 ✗ void range_assertion(int k, int m, int n, const char* msg)
103 {
104 ✗ if ( k < m || k > n)
105 {
106 std::cerr << "range assertion failed: \n"
107 ✗ << msg << std::endl
108 ✗ << "value: " << k
109 ✗ << " range: " << m << ", " << n << std::endl;
110 ✗ assert (k >= m && k <= n);
111 }
112 ✗ }
113
114
115 ////////////////////////////////////////////////////////////////////////////////
116 // numerical routines
117
118 /*
119 * Compute the determinant of the 2x2 matrix with column the point P1, P2
120 */
121 ✗ double det(Point const& P1, Point const& P2)
122 {
123 ✗ return P1[X]*P2[Y] - P1[Y]*P2[X];
124 }
125
126 /*
127 * Solve the linear system [P1,P2] * P = Q
128 * in case there isn't exactly one solution the routine returns false
129 */
130 ✗ bool solve(Point & P, Point const& P1, Point const& P2, Point const& Q)
131 {
132 ✗ double d = det(P1, P2);
133 ✗ if (d == 0) return false;
134 ✗ d = 1 / d;
135 ✗ P[X] = det(Q, P2) * d;
136 ✗ P[Y] = det(P1, Q) * d;
137 ✗ return true;
138 }
139
140 ////////////////////////////////////////////////////////////////////////////////
141 // interval routines
142
143 /*
144 * Map the sub-interval I in [0,1] into the interval J and assign it to J
145 */
146 966 void map_to(Interval & J, Interval const& I)
147 {
148 966 J.setEnds(J.valueAt(I.min()), J.valueAt(I.max()));
149 966 }
150
151 ////////////////////////////////////////////////////////////////////////////////
152 // bezier curve routines
153
154 /*
155 * Return true if all the Bezier curve control points are near,
156 * false otherwise
157 */
158 // Bezier.isConstant(precision)
159 2121 bool is_constant(std::vector<Point> const& A, double precision)
160 {
161
2/2
✓ Branch 1 taken 3009 times.
✓ Branch 2 taken 387 times.
 3396 for (unsigned int i = 1; i < A.size(); ++i)
162 {
163
2/2
✓ Branch 3 taken 1734 times.
✓ Branch 4 taken 1275 times.
 3009 if(!are_near(A[i], A[0], precision))
164 1734 return false;
165 }
166 387 return true;
167 }
168
169 /*
170 * Compute the hodograph of the bezier curve B and return it in D
171 */
172 // derivative(Bezier)
173 ✗ void derivative(std::vector<Point> & D, std::vector<Point> const& B)
174 {
175 ✗ D.clear();
176 ✗ size_t sz = B.size();
177 ✗ if (sz == 0) return;
178 ✗ if (sz == 1)
179 {
180 ✗ D.resize(1, Point(0,0));
181 ✗ return;
182 }
183 ✗ size_t n = sz-1;
184 ✗ D.reserve(n);
185 ✗ for (size_t i = 0; i < n; ++i)
186 {
187 ✗ D.push_back(n*(B[i+1] - B[i]));
188 }
189 }
190
191 /*
192 * Compute the hodograph of the Bezier curve B rotated of 90 degree
193 * and return it in D; we have N(t) orthogonal to B(t) for any t
194 */
195 // rot90(derivative(Bezier))
196 ✗ void normal(std::vector<Point> & N, std::vector<Point> const& B)
197 {
198 ✗ derivative(N,B);
199 ✗ for (auto & i : N)
200 {
201 ✗ i = rot90(i);
202 }
203 ✗ }
204
205 /*
206 * Compute the portion of the Bezier curve "B" wrt the interval [0,t]
207 */
208 // portion(Bezier, 0, t)
209 584 void left_portion(Coord t, std::vector<Point> & B)
210 {
211 584 size_t n = B.size();
212
2/2
✓ Branch 0 taken 2061 times.
✓ Branch 1 taken 584 times.
 2645 for (size_t i = 1; i < n; ++i)
213 {
214
2/2
✓ Branch 0 taken 5109 times.
✓ Branch 1 taken 2061 times.
 7170 for (size_t j = n-1; j > i-1 ; --j)
215 {
216 5109 B[j] = lerp(t, B[j-1], B[j]);
217 }
218 }
219 584 }
220
221 /*
222 * Compute the portion of the Bezier curve "B" wrt the interval [t,1]
223 */
224 // portion(Bezier, t, 1)
225 567 void right_portion(Coord t, std::vector<Point> & B)
226 {
227 567 size_t n = B.size();
228
2/2
✓ Branch 0 taken 2013 times.
✓ Branch 1 taken 567 times.
 2580 for (size_t i = 1; i < n; ++i)
229 {
230
2/2
✓ Branch 0 taken 5017 times.
✓ Branch 1 taken 2013 times.
 7030 for (size_t j = 0; j < n-i; ++j)
231 {
232 5017 B[j] = lerp(t, B[j], B[j+1]);
233 }
234 }
235 567 }
236
237 /*
238 * Compute the portion of the Bezier curve "B" wrt the interval "I"
239 */
240 // portion(Bezier, I)
241 966 void portion (std::vector<Point> & B , Interval const& I)
242 {
243
2/2
✓ Branch 1 taken 399 times.
✓ Branch 2 taken 567 times.
 966 if (I.min() == 0)
244 {
245
2/2
✓ Branch 1 taken 133 times.
✓ Branch 2 taken 266 times.
 399 if (I.max() == 1) return;
246 266 left_portion(I.max(), B);
247 266 return;
248 }
249 567 right_portion(I.min(), B);
250
2/2
✓ Branch 1 taken 249 times.
✓ Branch 2 taken 318 times.
 567 if (I.max() == 1) return;
251 318 double t = I.extent() / (1 - I.min());
252 318 left_portion(t, B);
253 }
254
255
256 ////////////////////////////////////////////////////////////////////////////////
257 // tags
258
259 struct intersection_point_tag;
260 struct collinear_normal_tag;
261 template <typename Tag>
262 OptInterval clip(std::vector<Point> const& A,
263 std::vector<Point> const& B,
264 double precision);
265 template <typename Tag>
266 void iterate(std::vector<Interval>& domsA,
267 std::vector<Interval>& domsB,
268 std::vector<Point> const& A,
269 std::vector<Point> const& B,
270 Interval const& domA,
271 Interval const& domB,
272 double precision );
273
274
275 ////////////////////////////////////////////////////////////////////////////////
276 // intersection
277
278 /*
279 * Make up an orientation line using the control points c[i] and c[j]
280 * the line is returned in the output parameter "l" in the form of a 3 element
281 * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
282 */
283 // Line(c[i], c[j])
284 ✗ void orientation_line (std::vector<double> & l,
285 std::vector<Point> const& c,
286 size_t i, size_t j)
287 {
288 ✗ l[0] = c[j][Y] - c[i][Y];
289 ✗ l[1] = c[i][X] - c[j][X];
290 ✗ l[2] = cross(c[j], c[i]);
291 ✗ double length = std::sqrt(l[0] * l[0] + l[1] * l[1]);
292 ✗ assert (length != 0);
293 ✗ l[0] /= length;
294 ✗ l[1] /= length;
295 ✗ l[2] /= length;
296 ✗ }
297
298 /*
299 * Pick up an orientation line for the Bezier curve "c" and return it in
300 * the output parameter "l"
301 */
302 697 Line pick_orientation_line (std::vector<Point> const &c, double precision)
303 {
304 697 size_t i = c.size();
305
3/6
✓ Branch 0 taken 697 times.
✗ Branch 1 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 697 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 697 times.
 697 while (--i > 0 && are_near(c[0], c[i], precision))
306 {}
307
308 // this should never happen because when a new curve portion is created
309 // we check that it is not constant;
310 // however this requires that the precision used in the is_constant
311 // routine has to be the same used here in the are_near test
312
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 697 times.
 697 assert(i != 0);
313
314 697 Line line(c[0], c[i]);
315 697 return line;
316 //std::cerr << "i = " << i << std::endl;
317 }
318
319 /*
320 * Make up an orientation line for constant bezier curve;
321 * the orientation line is made up orthogonal to the other curve base line;
322 * the line is returned in the output parameter "l" in the form of a 3 element
323 * vector : l[0] * x + l[1] * y + l[2] == 0; the line is normalized.
324 */
325 76 Line orthogonal_orientation_line (std::vector<Point> const &c,
326 Point const &p,
327 double precision)
328 {
329 // this should never happen
330
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 76 times.
 76 assert(!is_constant(c, precision));
331
332 76 Line line(p, (c.back() - c.front()).cw() + p);
333 76 return line;
334 }
335
336 /*
337 * Compute the signed distance of the point "P" from the normalized line l
338 */
339 6813 double signed_distance(Point const &p, Line const &l)
340 {
341 6813 Coord a, b, c;
342
1/2
✓ Branch 1 taken 6813 times.
✗ Branch 2 not taken.
 6813 l.coefficients(a, b, c);
343 13626 return a * p[X] + b * p[Y] + c;
344 }
345
346 /*
347 * Compute the min and max distance of the control points of the Bezier
348 * curve "c" from the normalized orientation line "l".
349 * This bounds are returned through the output Interval parameter"bound".
350 */
351 773 Interval fat_line_bounds (std::vector<Point> const &c,
352 Line const &l)
353 {
354 773 Interval bound(0, 0);
355
2/2
✓ Branch 8 taken 3448 times.
✓ Branch 9 taken 773 times.
 4221 for (auto i : c) {
356
1/2
✓ Branch 1 taken 3448 times.
✗ Branch 2 not taken.
 3448 bound.expandTo(signed_distance(i, l));
357 }
358 773 return bound;
359 }
360
361 /*
362 * return the x component of the intersection point between the line
363 * passing through points p1, p2 and the line Y = "y"
364 */
365 1756 double intersect (Point const& p1, Point const& p2, double y)
366 {
367 // we are sure that p2[Y] != p1[Y] because this routine is called
368 // only when the lower or the upper bound is crossed
369 1756 double dy = (p2[Y] - p1[Y]);
370 1756 double s = (y - p1[Y]) / dy;
371 1756 return (p2[X]-p1[X])*s + p1[X];
372 }
373
374 /*
375 * Clip the Bezier curve "B" wrt the fat line defined by the orientation
376 * line "l" and the interval range "bound", the new parameter interval for
377 * the clipped curve is returned through the output parameter "dom"
378 */
379 773 OptInterval clip_interval (std::vector<Point> const& B,
380 Line const &l,
381 Interval const &bound)
382 {
383 773 double n = B.size() - 1; // number of sub-intervals
384 773 std::vector<Point> D; // distance curve control points
385
1/2
✓ Branch 2 taken 773 times.
✗ Branch 3 not taken.
 773 D.reserve (B.size());
386
2/2
✓ Branch 1 taken 3365 times.
✓ Branch 2 taken 773 times.
 4138 for (size_t i = 0; i < B.size(); ++i)
387 {
388
1/2
✓ Branch 3 taken 3365 times.
✗ Branch 4 not taken.
 3365 const double d = signed_distance(B[i], l);
389
1/2
✓ Branch 2 taken 3365 times.
✗ Branch 3 not taken.
 3365 D.emplace_back(i/n, d);
390 }
391 //print(D);
392
393 773 ConvexHull p;
394
1/2
✓ Branch 1 taken 773 times.
✗ Branch 2 not taken.
 773 p.swap(D);
395 //print(p);
396
397 bool plower, phigher;
398 bool clower, chigher;
399 773 double t, tmin = 1, tmax = 0;
400 // std::cerr << "bound : " << bound << std::endl;
401
402 773 plower = (p[0][Y] < bound.min());
403 773 phigher = (p[0][Y] > bound.max());
404
4/4
✓ Branch 0 taken 491 times.
✓ Branch 1 taken 282 times.
✓ Branch 2 taken 232 times.
✓ Branch 3 taken 259 times.
 773 if (!(plower || phigher)) // inside the fat line
405 {
406
1/2
✓ Branch 2 taken 232 times.
✗ Branch 3 not taken.
 232 if (tmin > p[0][X]) tmin = p[0][X];
407
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 232 times.
 232 if (tmax < p[0][X]) tmax = p[0][X];
408 // std::cerr << "0 : inside " << p[0]
409 // << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
410 }
411
412
2/2
✓ Branch 1 taken 2103 times.
✓ Branch 2 taken 773 times.
 2876 for (size_t i = 1; i < p.size(); ++i)
413 {
414 2103 clower = (p[i][Y] < bound.min());
415 2103 chigher = (p[i][Y] > bound.max());
416
4/4
✓ Branch 0 taken 1133 times.
✓ Branch 1 taken 970 times.
✓ Branch 2 taken 351 times.
✓ Branch 3 taken 782 times.
 2103 if (!(clower || chigher)) // inside the fat line
417 {
418
2/2
✓ Branch 2 taken 31 times.
✓ Branch 3 taken 320 times.
 351 if (tmin > p[i][X]) tmin = p[i][X];
419
2/2
✓ Branch 2 taken 285 times.
✓ Branch 3 taken 66 times.
 351 if (tmax < p[i][X]) tmax = p[i][X];
420 // std::cerr << i << " : inside " << p[i]
421 // << " : tmin = " << tmin << ", tmax = " << tmax
422 // << std::endl;
423 }
424
2/2
✓ Branch 0 taken 628 times.
✓ Branch 1 taken 1475 times.
 2103 if (clower != plower) // cross the lower bound
425 {
426
1/2
✓ Branch 4 taken 628 times.
✗ Branch 5 not taken.
 628 t = intersect(p[i-1], p[i], bound.min());
427
2/2
✓ Branch 0 taken 359 times.
✓ Branch 1 taken 269 times.
 628 if (tmin > t) tmin = t;
428
2/2
✓ Branch 0 taken 520 times.
✓ Branch 1 taken 108 times.
 628 if (tmax < t) tmax = t;
429 628 plower = clower;
430 // std::cerr << i << " : lower " << p[i]
431 // << " : tmin = " << tmin << ", tmax = " << tmax
432 // << std::endl;
433 }
434
2/2
✓ Branch 0 taken 575 times.
✓ Branch 1 taken 1528 times.
 2103 if (chigher != phigher) // cross the upper bound
435 {
436
1/2
✓ Branch 4 taken 575 times.
✗ Branch 5 not taken.
 575 t = intersect(p[i-1], p[i], bound.max());
437
2/2
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 335 times.
 575 if (tmin > t) tmin = t;
438
2/2
✓ Branch 0 taken 163 times.
✓ Branch 1 taken 412 times.
 575 if (tmax < t) tmax = t;
439 575 phigher = chigher;
440 // std::cerr << i << " : higher " << p[i]
441 // << " : tmin = " << tmin << ", tmax = " << tmax
442 // << std::endl;
443 }
444 }
445
446 // we have to test the closing segment for intersection
447 773 size_t last = p.size() - 1;
448 773 clower = (p[0][Y] < bound.min());
449 773 chigher = (p[0][Y] > bound.max());
450
2/2
✓ Branch 0 taken 254 times.
✓ Branch 1 taken 519 times.
 773 if (clower != plower) // cross the lower bound
451 {
452
1/2
✓ Branch 4 taken 254 times.
✗ Branch 5 not taken.
 254 t = intersect(p[last], p[0], bound.min());
453
2/2
✓ Branch 0 taken 137 times.
✓ Branch 1 taken 117 times.
 254 if (tmin > t) tmin = t;
454
2/2
✓ Branch 0 taken 105 times.
✓ Branch 1 taken 149 times.
 254 if (tmax < t) tmax = t;
455 // std::cerr << "0 : lower " << p[0]
456 // << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
457 }
458
2/2
✓ Branch 0 taken 299 times.
✓ Branch 1 taken 474 times.
 773 if (chigher != phigher) // cross the upper bound
459 {
460
1/2
✓ Branch 4 taken 299 times.
✗ Branch 5 not taken.
 299 t = intersect(p[last], p[0], bound.max());
461
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 299 times.
 299 if (tmin > t) tmin = t;
462
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 299 times.
 299 if (tmax < t) tmax = t;
463 // std::cerr << "0 : higher " << p[0]
464 // << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
465 }
466
467
4/4
✓ Branch 0 taken 163 times.
✓ Branch 1 taken 610 times.
✓ Branch 2 taken 141 times.
✓ Branch 3 taken 22 times.
 773 if (tmin == 1 && tmax == 0) {
468 141 return OptInterval();
469 } else {
470
2/4
✓ Branch 2 taken 632 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 632 times.
✗ Branch 6 not taken.
 1264 return Interval(tmin, tmax);
471 }
472 773 }
473
474 /*
475 * Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating
476 * intersection points the new parameter interval for the clipped curve
477 * is returned through the output parameter "dom"
478 */
479 template <>
480 773 OptInterval clip<intersection_point_tag> (std::vector<Point> const& A,
481 std::vector<Point> const& B,
482 double precision)
483 {
484 773 Line bl;
485
3/4
✓ Branch 1 taken 773 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 76 times.
✓ Branch 4 taken 697 times.
 773 if (is_constant(A, precision)) {
486
1/2
✓ Branch 4 taken 76 times.
✗ Branch 5 not taken.
 76 Point M = middle_point(A.front(), A.back());
487
1/2
✓ Branch 1 taken 76 times.
✗ Branch 2 not taken.
 76 bl = orthogonal_orientation_line(B, M, precision);
488 } else {
489
1/2
✓ Branch 1 taken 697 times.
✗ Branch 2 not taken.
 697 bl = pick_orientation_line(A, precision);
490 }
491
1/2
✓ Branch 1 taken 773 times.
✗ Branch 2 not taken.
 773 bl.normalize();
492
1/2
✓ Branch 2 taken 773 times.
✗ Branch 3 not taken.
 773 Interval bound = fat_line_bounds(A, bl);
493
1/2
✓ Branch 1 taken 773 times.
✗ Branch 2 not taken.
 1546 return clip_interval(B, bl, bound);
494 }
495
496
497 ///////////////////////////////////////////////////////////////////////////////
498 // collinear normal
499
500 /*
501 * Compute a closed focus for the Bezier curve B and return it in F
502 * A focus is any curve through which all lines perpendicular to B(t) pass.
503 */
504 ✗ void make_focus (std::vector<Point> & F, std::vector<Point> const& B)
505 {
506 ✗ assert (B.size() > 2);
507 ✗ size_t n = B.size() - 1;
508 ✗ normal(F, B);
509 ✗ Point c(1, 1);
510 #if VERBOSE
511 if (!solve(c, F[0], -F[n-1], B[n]-B[0]))
512 {
513 std::cerr << "make_focus: unable to make up a closed focus" << std::endl;
514 }
515 #else
516 ✗ solve(c, F[0], -F[n-1], B[n]-B[0]);
517 #endif
518 // std::cerr << "c = " << c << std::endl;
519
520
521 // B(t) + c(t) * N(t)
522 ✗ double n_inv = 1 / (double)(n);
523 ✗ Point c0ni;
524 ✗ F.push_back(c[1] * F[n-1]);
525 ✗ F[n] += B[n];
526 ✗ for (size_t i = n-1; i > 0; --i)
527 {
528 ✗ F[i] *= -c[0];
529 ✗ c0ni = F[i];
530 ✗ F[i] += (c[1] * F[i-1]);
531 ✗ F[i] *= (i * n_inv);
532 ✗ F[i] -= c0ni;
533 ✗ F[i] += B[i];
534 }
535 ✗ F[0] *= c[0];
536 ✗ F[0] += B[0];
537 ✗ }
538
539 /*
540 * Compute the projection on the plane (t, d) of the control points
541 * (t, u, D(t,u)) where D(t,u) = <(B(t) - F(u)), B'(t)> with 0 <= t, u <= 1
542 * B is a Bezier curve and F is a focus of another Bezier curve.
543 * See Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping.
544 */
545 ✗ void distance_control_points (std::vector<Point> & D,
546 std::vector<Point> const& B,
547 std::vector<Point> const& F)
548 {
549 ✗ assert (B.size() > 1);
550 ✗ assert (!F.empty());
551 ✗ const size_t n = B.size() - 1;
552 ✗ const size_t m = F.size() - 1;
553 ✗ const size_t r = 2 * n - 1;
554 ✗ const double r_inv = 1 / (double)(r);
555 ✗ D.clear();
556 ✗ D.reserve (B.size() * F.size());
557
558 ✗ std::vector<Point> dB;
559 ✗ dB.reserve(n);
560 ✗ for (size_t k = 0; k < n; ++k)
561 {
562 ✗ dB.push_back (B[k+1] - B[k]);
563 }
564 ✗ NL::Matrix dBB(n,B.size());
565 ✗ for (size_t i = 0; i < n; ++i)
566 ✗ for (size_t j = 0; j < B.size(); ++j)
567 ✗ dBB(i,j) = dot (dB[i], B[j]);
568 ✗ NL::Matrix dBF(n, F.size());
569 ✗ for (size_t i = 0; i < n; ++i)
570 ✗ for (size_t j = 0; j < F.size(); ++j)
571 ✗ dBF(i,j) = dot (dB[i], F[j]);
572
573 size_t l;
574 double bc;
575 ✗ Point dij;
576 ✗ std::vector<double> d(F.size());
577 ✗ int rci = 1;
578 ✗ int b1 = 1;
579 ✗ for (size_t i = 0; i <= r; ++i)
580 {
581 ✗ for (size_t j = 0; j <= m; ++j)
582 {
583 ✗ d[j] = 0;
584 }
585 ✗ const size_t k0 = std::max(i, n) - n;
586 ✗ const size_t kn = std::min(i, n-1);
587 ✗ const double bri = (double)n / rci;
588
589 // assert(rci == binomial(r, i));
590 ✗ binomial_increment_k(rci, r, i);
591
592 ✗ int b2 = b1;
593 ✗ for (size_t k = k0; k <= kn; ++k)
594 {
595 //if (k > i || (i-k) > n) continue;
596 ✗ l = i - k;
597 #if CHECK
598 assert (l <= n);
599 #endif
600 ✗ bc = bri * b2;
601
602 // assert(b2 == binomial(n, l) * binomial(n - 1, k));
603 ✗ binomial_decrement_k(b2, n, l);
604 ✗ binomial_increment_k(b2, n - 1, k);
605
606 ✗ for (size_t j = 0; j <= m; ++j)
607 {
608 //d[j] += bc * dot(dB[k], B[l] - F[j]);
609 ✗ d[j] += bc * (dBB(k,l) - dBF(k,j));
610 }
611 }
612
613 // assert(b1 == binomial(n, i - k0) * binomial(n - 1, k0));
614 ✗ if (i < n) {
615 ✗ binomial_increment_k(b1, n, i);
616 } else {
617 ✗ binomial_increment_k(b1, n - 1, k0);
618 }
619
620 double dmin, dmax;
621 ✗ dmin = dmax = d[m];
622 ✗ for (size_t j = 0; j < m; ++j)
623 {
624 ✗ if (dmin > d[j]) dmin = d[j];
625 ✗ if (dmax < d[j]) dmax = d[j];
626 }
627 ✗ dij[0] = i * r_inv;
628 ✗ dij[1] = dmin;
629 ✗ D.push_back (dij);
630 ✗ dij[1] = dmax;
631 ✗ D.push_back (dij);
632 }
633 ✗ }
634
635 /*
636 * Clip the Bezier curve "B" wrt the focus "F"; the new parameter interval for
637 * the clipped curve is returned through the output parameter "dom"
638 */
639 ✗ OptInterval clip_interval (std::vector<Point> const& B,
640 std::vector<Point> const& F)
641 {
642 ✗ std::vector<Point> D; // distance curve control points
643 ✗ distance_control_points(D, B, F);
644 //print(D, "D");
645 // ConvexHull chD(D);
646 // std::vector<Point>& p = chD.boundary; // convex hull vertices
647
648 ✗ ConvexHull p;
649 ✗ p.swap(D);
650 //print(p, "CH(D)");
651
652 bool plower, clower;
653 ✗ double t, tmin = 1, tmax = 0;
654
655 ✗ plower = (p[0][Y] < 0);
656 ✗ if (p[0][Y] == 0) // on the x axis
657 {
658 ✗ if (tmin > p[0][X]) tmin = p[0][X];
659 ✗ if (tmax < p[0][X]) tmax = p[0][X];
660 // std::cerr << "0 : on x axis " << p[0]
661 // << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
662 }
663
664 ✗ for (size_t i = 1; i < p.size(); ++i)
665 {
666 ✗ clower = (p[i][Y] < 0);
667 ✗ if (p[i][Y] == 0) // on x axis
668 {
669 ✗ if (tmin > p[i][X]) tmin = p[i][X];
670 ✗ if (tmax < p[i][X]) tmax = p[i][X];
671 // std::cerr << i << " : on x axis " << p[i]
672 // << " : tmin = " << tmin << ", tmax = " << tmax
673 // << std::endl;
674 }
675 ✗ else if (clower != plower) // cross the x axis
676 {
677 ✗ t = intersect(p[i-1], p[i], 0);
678 ✗ if (tmin > t) tmin = t;
679 ✗ if (tmax < t) tmax = t;
680 ✗ plower = clower;
681 // std::cerr << i << " : lower " << p[i]
682 // << " : tmin = " << tmin << ", tmax = " << tmax
683 // << std::endl;
684 }
685 }
686
687 // we have to test the closing segment for intersection
688 ✗ size_t last = p.size() - 1;
689 ✗ clower = (p[0][Y] < 0);
690 ✗ if (clower != plower) // cross the x axis
691 {
692 ✗ t = intersect(p[last], p[0], 0);
693 ✗ if (tmin > t) tmin = t;
694 ✗ if (tmax < t) tmax = t;
695 // std::cerr << "0 : lower " << p[0]
696 // << " : tmin = " << tmin << ", tmax = " << tmax << std::endl;
697 }
698 ✗ if (tmin == 1 && tmax == 0) {
699 ✗ return OptInterval();
700 } else {
701 ✗ return Interval(tmin, tmax);
702 }
703 ✗ }
704
705 /*
706 * Clip the Bezier curve "B" wrt the Bezier curve "A" for individuating
707 * points which have collinear normals; the new parameter interval
708 * for the clipped curve is returned through the output parameter "dom"
709 */
710 template <>
711 ✗ OptInterval clip<collinear_normal_tag> (std::vector<Point> const& A,
712 std::vector<Point> const& B,
713 double /*precision*/)
714 {
715 ✗ std::vector<Point> F;
716 ✗ make_focus(F, A);
717 ✗ return clip_interval(B, F);
718 ✗ }
719
720
721
722 const double MAX_PRECISION = 1e-8;
723 const double MIN_CLIPPED_SIZE_THRESHOLD = 0.8;
724 const Interval UNIT_INTERVAL(0,1);
725 const OptInterval EMPTY_INTERVAL;
726 const Interval H1_INTERVAL(0, 0.5);
727 const Interval H2_INTERVAL(nextafter(0.5, 1.0), 1.0);
728
729 /*
730 * iterate
731 *
732 * input:
733 * A, B: control point sets of two bezier curves
734 * domA, domB: real parameter intervals of the two curves
735 * precision: required computational precision of the returned parameter ranges
736 * output:
737 * domsA, domsB: sets of parameter intervals
738 *
739 * The parameter intervals are computed by using a Bezier clipping algorithm,
740 * in case the clipping doesn't shrink the initial interval more than 20%,
741 * a subdivision step is performed.
742 * If during the computation both curves collapse to a single point
743 * the routine exits independently by the precision reached in the computation
744 * of the curve intervals.
745 */
746 template <>
747 392 void iterate<intersection_point_tag> (std::vector<Interval>& domsA,
748 std::vector<Interval>& domsB,
749 std::vector<Point> const& A,
750 std::vector<Point> const& B,
751 Interval const& domA,
752 Interval const& domB,
753 double precision )
754 {
755 // in order to limit recursion
756 static size_t counter = 0;
757
7/10
✓ Branch 1 taken 392 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 58 times.
✓ Branch 4 taken 334 times.
✓ Branch 6 taken 58 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 58 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 58 times.
✓ Branch 11 taken 334 times.
 392 if (domA.extent() == 1 && domB.extent() == 1) counter = 0;
758
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 392 times.
 392 if (++counter > 100) return;
759 #if VERBOSE
760 std::cerr << std::fixed << std::setprecision(16);
761 std::cerr << ">> curve subdision performed <<" << std::endl;
762 std::cerr << "dom(A) : " << domA << std::endl;
763 std::cerr << "dom(B) : " << domB << std::endl;
764 // std::cerr << "angle(A) : " << angle(A) << std::endl;
765 // std::cerr << "angle(B) : " << angle(B) << std::endl;
766 #endif
767
768
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 392 times.
 392 if (precision < MAX_PRECISION)
769 ✗ precision = MAX_PRECISION;
770
771
1/2
✓ Branch 2 taken 392 times.
✗ Branch 3 not taken.
 392 std::vector<Point> pA = A;
772
1/2
✓ Branch 2 taken 392 times.
✗ Branch 3 not taken.
 392 std::vector<Point> pB = B;
773 392 std::vector<Point>* C1 = &pA;
774 392 std::vector<Point>* C2 = &pB;
775
776 392 Interval dompA = domA;
777 392 Interval dompB = domB;
778 392 Interval* dom1 = &dompA;
779 392 Interval* dom2 = &dompB;
780
781 392 OptInterval dom;
782
783
3/10
✓ Branch 1 taken 392 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 392 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 392 times.
 392 if ( is_constant(A, precision) && is_constant(B, precision) ){
784 ✗ Point M1 = middle_point(C1->front(), C1->back());
785 ✗ Point M2 = middle_point(C2->front(), C2->back());
786 ✗ if (are_near(M1,M2)){
787 ✗ domsA.push_back(domA);
788 ✗ domsB.push_back(domB);
789 }
790 ✗ return;
791 }
792
793 392 size_t iter = 0;
794 392 while (++iter < 100
795
9/12
✓ Branch 0 taken 794 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 794 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 58 times.
✓ Branch 6 taken 736 times.
✓ Branch 8 taken 58 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 37 times.
✓ Branch 11 taken 21 times.
✓ Branch 12 taken 773 times.
✓ Branch 13 taken 21 times.
 794 && (dompA.extent() >= precision || dompB.extent() >= precision))
796 {
797 #if VERBOSE
798 std::cerr << "iter: " << iter << std::endl;
799 #endif
800
1/2
✓ Branch 1 taken 773 times.
✗ Branch 2 not taken.
 773 dom = clip<intersection_point_tag>(*C1, *C2, precision);
801
802
2/2
✓ Branch 1 taken 141 times.
✓ Branch 2 taken 632 times.
 773 if (dom.empty())
803 {
804 #if VERBOSE
805 std::cerr << "dom: empty" << std::endl;
806 #endif
807 141 return;
808 }
809 #if VERBOSE
810 std::cerr << "dom : " << dom << std::endl;
811 #endif
812 // all other cases where dom[0] > dom[1] are invalid
813
1/2
✗ Branch 4 not taken.
✓ Branch 5 taken 632 times.
 632 assert(dom->min() <= dom->max());
814
815
1/2
✓ Branch 2 taken 632 times.
✗ Branch 3 not taken.
 632 map_to(*dom2, *dom);
816
817
1/2
✓ Branch 2 taken 632 times.
✗ Branch 3 not taken.
 632 portion(*C2, *dom);
818
8/10
✓ Branch 1 taken 632 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 248 times.
✓ Branch 4 taken 384 times.
✓ Branch 6 taken 248 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 63 times.
✓ Branch 9 taken 185 times.
✓ Branch 10 taken 63 times.
✓ Branch 11 taken 569 times.
 632 if (is_constant(*C2, precision) && is_constant(*C1, precision))
819 {
820
1/2
✓ Branch 4 taken 63 times.
✗ Branch 5 not taken.
 63 Point M1 = middle_point(C1->front(), C1->back());
821
1/2
✓ Branch 4 taken 63 times.
✗ Branch 5 not taken.
 63 Point M2 = middle_point(C2->front(), C2->back());
822 #if VERBOSE
823 std::cerr << "both curves are constant: \n"
824 << "M1: " << M1 << "\n"
825 << "M2: " << M2 << std::endl;
826 print(*C2, "C2");
827 print(*C1, "C1");
828 #endif
829
2/4
✓ Branch 1 taken 63 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 63 times.
✗ Branch 4 not taken.
 63 if (are_near(M1,M2))
830 63 break; // append the new interval
831 else
832 ✗ return; // exit without appending any new interval
833 }
834
835
836 // if we have clipped less than 20% than we need to subdive the curve
837 // with the largest domain into two sub-curves
838
3/4
✓ Branch 2 taken 569 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 167 times.
✓ Branch 5 taken 402 times.
 569 if (dom->extent() > MIN_CLIPPED_SIZE_THRESHOLD)
839 {
840 #if VERBOSE
841 std::cerr << "clipped less than 20% : " << dom->extent() << std::endl;
842 std::cerr << "angle(pA) : " << angle(pA) << std::endl;
843 std::cerr << "angle(pB) : " << angle(pB) << std::endl;
844 #endif
845 167 std::vector<Point> pC1, pC2;
846 167 Interval dompC1, dompC2;
847
4/6
✓ Branch 1 taken 167 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 167 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 113 times.
✓ Branch 7 taken 54 times.
 167 if (dompA.extent() > dompB.extent())
848 {
849
2/4
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 113 times.
✗ Branch 5 not taken.
 113 pC1 = pC2 = pA;
850
1/2
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
 113 portion(pC1, H1_INTERVAL);
851
1/2
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
 113 portion(pC2, H2_INTERVAL);
852 113 dompC1 = dompC2 = dompA;
853
1/2
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
 113 map_to(dompC1, H1_INTERVAL);
854
1/2
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
 113 map_to(dompC2, H2_INTERVAL);
855
1/2
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
 113 iterate<intersection_point_tag>(domsA, domsB, pC1, pB,
856 dompC1, dompB, precision);
857
1/2
✓ Branch 1 taken 113 times.
✗ Branch 2 not taken.
 113 iterate<intersection_point_tag>(domsA, domsB, pC2, pB,
858 dompC2, dompB, precision);
859 }
860 else
861 {
862
2/4
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 54 times.
✗ Branch 5 not taken.
 54 pC1 = pC2 = pB;
863
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 portion(pC1, H1_INTERVAL);
864
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 portion(pC2, H2_INTERVAL);
865 54 dompC1 = dompC2 = dompB;
866
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 map_to(dompC1, H1_INTERVAL);
867
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 map_to(dompC2, H2_INTERVAL);
868
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 iterate<intersection_point_tag>(domsB, domsA, pC1, pA,
869 dompC1, dompA, precision);
870
1/2
✓ Branch 1 taken 54 times.
✗ Branch 2 not taken.
 54 iterate<intersection_point_tag>(domsB, domsA, pC2, pA,
871 dompC2, dompA, precision);
872 }
873 167 return;
874 167 }
875
876 402 swap(C1, C2);
877 402 swap(dom1, dom2);
878 #if VERBOSE
879 std::cerr << "dom(pA) : " << dompA << std::endl;
880 std::cerr << "dom(pB) : " << dompB << std::endl;
881 #endif
882 }
883
1/2
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
 84 domsA.push_back(dompA);
884
1/2
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
 84 domsB.push_back(dompB);
885
4/4
✓ Branch 1 taken 84 times.
✓ Branch 2 taken 308 times.
✓ Branch 4 taken 84 times.
✓ Branch 5 taken 308 times.
 700 }
886
887
888 /*
889 * iterate
890 *
891 * input:
892 * A, B: control point sets of two bezier curves
893 * domA, domB: real parameter intervals of the two curves
894 * precision: required computational precision of the returned parameter ranges
895 * output:
896 * domsA, domsB: sets of parameter intervals
897 *
898 * The parameter intervals are computed by using a Bezier clipping algorithm,
899 * in case the clipping doesn't shrink the initial interval more than 20%,
900 * a subdivision step is performed.
901 * If during the computation one of the two curve interval length becomes less
902 * than MAX_PRECISION the routine exits independently by the precision reached
903 * in the computation of the other curve interval.
904 */
905 template <>
906 ✗ void iterate<collinear_normal_tag> (std::vector<Interval>& domsA,
907 std::vector<Interval>& domsB,
908 std::vector<Point> const& A,
909 std::vector<Point> const& B,
910 Interval const& domA,
911 Interval const& domB,
912 double precision)
913 {
914 // in order to limit recursion
915 static size_t counter = 0;
916 ✗ if (domA.extent() == 1 && domB.extent() == 1) counter = 0;
917 ✗ if (++counter > 100) return;
918 #if VERBOSE
919 std::cerr << std::fixed << std::setprecision(16);
920 std::cerr << ">> curve subdision performed <<" << std::endl;
921 std::cerr << "dom(A) : " << domA << std::endl;
922 std::cerr << "dom(B) : " << domB << std::endl;
923 // std::cerr << "angle(A) : " << angle(A) << std::endl;
924 // std::cerr << "angle(B) : " << angle(B) << std::endl;
925 #endif
926
927 ✗ if (precision < MAX_PRECISION)
928 ✗ precision = MAX_PRECISION;
929
930 ✗ std::vector<Point> pA = A;
931 ✗ std::vector<Point> pB = B;
932 ✗ std::vector<Point>* C1 = &pA;
933 ✗ std::vector<Point>* C2 = &pB;
934
935 ✗ Interval dompA = domA;
936 ✗ Interval dompB = domB;
937 ✗ Interval* dom1 = &dompA;
938 ✗ Interval* dom2 = &dompB;
939
940 ✗ OptInterval dom;
941
942 ✗ size_t iter = 0;
943 ✗ while (++iter < 100
944 ✗ && (dompA.extent() >= precision || dompB.extent() >= precision))
945 {
946 #if VERBOSE
947 std::cerr << "iter: " << iter << std::endl;
948 #endif
949 ✗ dom = clip<collinear_normal_tag>(*C1, *C2, precision);
950
951 ✗ if (dom.empty()) {
952 #if VERBOSE
953 std::cerr << "dom: empty" << std::endl;
954 #endif
955 ✗ return;
956 }
957 #if VERBOSE
958 std::cerr << "dom : " << dom << std::endl;
959 #endif
960 ✗ assert(dom->min() <= dom->max());
961
962 ✗ map_to(*dom2, *dom);
963
964 // it's better to stop before losing computational precision
965 ✗ if (iter > 1 && (dom2->extent() <= MAX_PRECISION))
966 {
967 #if VERBOSE
968 std::cerr << "beyond max precision limit" << std::endl;
969 #endif
970 ✗ break;
971 }
972
973 ✗ portion(*C2, *dom);
974 ✗ if (iter > 1 && is_constant(*C2, precision))
975 {
976 #if VERBOSE
977 std::cerr << "new curve portion pC1 is constant" << std::endl;
978 #endif
979 ✗ break;
980 }
981
982
983 // if we have clipped less than 20% than we need to subdive the curve
984 // with the largest domain into two sub-curves
985 ✗ if ( dom->extent() > MIN_CLIPPED_SIZE_THRESHOLD)
986 {
987 #if VERBOSE
988 std::cerr << "clipped less than 20% : " << dom->extent() << std::endl;
989 std::cerr << "angle(pA) : " << angle(pA) << std::endl;
990 std::cerr << "angle(pB) : " << angle(pB) << std::endl;
991 #endif
992 ✗ std::vector<Point> pC1, pC2;
993 ✗ Interval dompC1, dompC2;
994 ✗ if (dompA.extent() > dompB.extent())
995 {
996 ✗ if ((dompA.extent() / 2) < MAX_PRECISION)
997 {
998 ✗ break;
999 }
1000 ✗ pC1 = pC2 = pA;
1001 ✗ portion(pC1, H1_INTERVAL);
1002 if (false && is_constant(pC1, precision))
1003 {
1004 #if VERBOSE
1005 std::cerr << "new curve portion pC1 is constant" << std::endl;
1006 #endif
1007 break;
1008 }
1009 ✗ portion(pC2, H2_INTERVAL);
1010 ✗ if (is_constant(pC2, precision))
1011 {
1012 #if VERBOSE
1013 std::cerr << "new curve portion pC2 is constant" << std::endl;
1014 #endif
1015 ✗ break;
1016 }
1017 ✗ dompC1 = dompC2 = dompA;
1018 ✗ map_to(dompC1, H1_INTERVAL);
1019 ✗ map_to(dompC2, H2_INTERVAL);
1020 ✗ iterate<collinear_normal_tag>(domsA, domsB, pC1, pB,
1021 dompC1, dompB, precision);
1022 ✗ iterate<collinear_normal_tag>(domsA, domsB, pC2, pB,
1023 dompC2, dompB, precision);
1024 }
1025 else
1026 {
1027 ✗ if ((dompB.extent() / 2) < MAX_PRECISION)
1028 {
1029 ✗ break;
1030 }
1031 ✗ pC1 = pC2 = pB;
1032 ✗ portion(pC1, H1_INTERVAL);
1033 ✗ if (is_constant(pC1, precision))
1034 {
1035 #if VERBOSE
1036 std::cerr << "new curve portion pC1 is constant" << std::endl;
1037 #endif
1038 ✗ break;
1039 }
1040 ✗ portion(pC2, H2_INTERVAL);
1041 ✗ if (is_constant(pC2, precision))
1042 {
1043 #if VERBOSE
1044 std::cerr << "new curve portion pC2 is constant" << std::endl;
1045 #endif
1046 ✗ break;
1047 }
1048 ✗ dompC1 = dompC2 = dompB;
1049 ✗ map_to(dompC1, H1_INTERVAL);
1050 ✗ map_to(dompC2, H2_INTERVAL);
1051 ✗ iterate<collinear_normal_tag>(domsB, domsA, pC1, pA,
1052 dompC1, dompA, precision);
1053 ✗ iterate<collinear_normal_tag>(domsB, domsA, pC2, pA,
1054 dompC2, dompA, precision);
1055 }
1056 ✗ return;
1057 ✗ }
1058
1059 ✗ swap(C1, C2);
1060 ✗ swap(dom1, dom2);
1061 #if VERBOSE
1062 std::cerr << "dom(pA) : " << dompA << std::endl;
1063 std::cerr << "dom(pB) : " << dompB << std::endl;
1064 #endif
1065 }
1066 ✗ domsA.push_back(dompA);
1067 ✗ domsB.push_back(dompB);
1068 ✗ }
1069
1070
1071 /*
1072 * get_solutions
1073 *
1074 * input: A, B - set of control points of two Bezier curve
1075 * input: precision - required precision of computation
1076 * input: clip - the routine used for clipping
1077 * output: xs - set of pairs of parameter values
1078 * at which the clipping algorithm converges
1079 *
1080 * This routine is based on the Bezier Clipping Algorithm,
1081 * see: Sederberg - Computer Aided Geometric Design
1082 */
1083 template <typename Tag>
1084 116 void get_solutions (std::vector< std::pair<double, double> >& xs,
1085 std::vector<Point> const& A,
1086 std::vector<Point> const& B,
1087 double precision)
1088 {
1089 116 std::pair<double, double> ci;
1090 116 std::vector<Interval> domsA, domsB;
1091
1/2
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
 116 iterate<Tag> (domsA, domsB, A, B, UNIT_INTERVAL, UNIT_INTERVAL, precision);
1092
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 58 times.
 116 if (domsA.size() != domsB.size())
1093 {
1094 ✗ assert (domsA.size() == domsB.size());
1095 }
1096 116 xs.clear();
1097
1/2
✓ Branch 2 taken 58 times.
✗ Branch 3 not taken.
 116 xs.reserve(domsA.size());
1098
2/2
✓ Branch 1 taken 84 times.
✓ Branch 2 taken 58 times.
 284 for (size_t i = 0; i < domsA.size(); ++i)
1099 {
1100 #if VERBOSE
1101 std::cerr << i << " : domA : " << domsA[i] << std::endl;
1102 std::cerr << "extent A: " << domsA[i].extent() << " ";
1103 std::cerr << "precision A: " << get_precision(domsA[i]) << std::endl;
1104 std::cerr << i << " : domB : " << domsB[i] << std::endl;
1105 std::cerr << "extent B: " << domsB[i].extent() << " ";
1106 std::cerr << "precision B: " << get_precision(domsB[i]) << std::endl;
1107 #endif
1108 168 ci.first = domsA[i].middle();
1109 168 ci.second = domsB[i].middle();
1110
1/2
✓ Branch 1 taken 84 times.
✗ Branch 2 not taken.
 168 xs.push_back(ci);
1111 }
1112 232 }
1113
1114 } /* end namespace bezier_clipping */ } /* end namespace detail */
1115
1116
1117 /*
1118 * find_collinear_normal
1119 *
1120 * input: A, B - set of control points of two Bezier curve
1121 * input: precision - required precision of computation
1122 * output: xs - set of pairs of parameter values
1123 * at which there are collinear normals
1124 *
1125 * This routine is based on the Bezier Clipping Algorithm,
1126 * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
1127 */
1128 ✗ void find_collinear_normal (std::vector< std::pair<double, double> >& xs,
1129 std::vector<Point> const& A,
1130 std::vector<Point> const& B,
1131 double precision)
1132 {
1133 using detail::bezier_clipping::get_solutions;
1134 using detail::bezier_clipping::collinear_normal_tag;
1135 ✗ get_solutions<collinear_normal_tag>(xs, A, B, precision);
1136 ✗ }
1137
1138
1139 /*
1140 * find_intersections_bezier_clipping
1141 *
1142 * input: A, B - set of control points of two Bezier curve
1143 * input: precision - required precision of computation
1144 * output: xs - set of pairs of parameter values
1145 * at which crossing happens
1146 *
1147 * This routine is based on the Bezier Clipping Algorithm,
1148 * see: Sederberg, Nishita, 1990 - Curve intersection using Bezier clipping
1149 */
1150 58 void find_intersections_bezier_clipping (std::vector< std::pair<double, double> >& xs,
1151 std::vector<Point> const& A,
1152 std::vector<Point> const& B,
1153 double precision)
1154 {
1155 using detail::bezier_clipping::get_solutions;
1156 using detail::bezier_clipping::intersection_point_tag;
1157 58 get_solutions<intersection_point_tag>(xs, A, B, precision);
1158 58 }
1159
1160 } // end namespace Geom
1161
1162
1163
1164
1165 /*
1166 Local Variables:
1167 mode:c++
1168 c-file-style:"stroustrup"
1169 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
1170 indent-tabs-mode:nil
1171 fill-column:99
1172 End:
1173 */
1174 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
1175