Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/matplotlib/bezier.py: 16%
Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2A module providing some utility functions regarding Bézier path manipulation.
3"""
5from functools import lru_cache
6import math
7import warnings
9import numpy as np
11from matplotlib import _api
14# same algorithm as 3.8's math.comb
15@np.vectorize
16@lru_cache(maxsize=128)
17def _comb(n, k):
18 if k > n:
19 return 0
20 k = min(k, n - k)
21 i = np.arange(1, k + 1)
22 return np.prod((n + 1 - i)/i).astype(int)
25class NonIntersectingPathException(ValueError):
26 pass
29# some functions
32def get_intersection(cx1, cy1, cos_t1, sin_t1,
33 cx2, cy2, cos_t2, sin_t2):
34 """
35 Return the intersection between the line through (*cx1*, *cy1*) at angle
36 *t1* and the line through (*cx2*, *cy2*) at angle *t2*.
37 """
39 # line1 => sin_t1 * (x - cx1) - cos_t1 * (y - cy1) = 0.
40 # line1 => sin_t1 * x + cos_t1 * y = sin_t1*cx1 - cos_t1*cy1
42 line1_rhs = sin_t1 * cx1 - cos_t1 * cy1
43 line2_rhs = sin_t2 * cx2 - cos_t2 * cy2
45 # rhs matrix
46 a, b = sin_t1, -cos_t1
47 c, d = sin_t2, -cos_t2
49 ad_bc = a * d - b * c
50 if abs(ad_bc) < 1e-12:
51 raise ValueError("Given lines do not intersect. Please verify that "
52 "the angles are not equal or differ by 180 degrees.")
54 # rhs_inverse
55 a_, b_ = d, -b
56 c_, d_ = -c, a
57 a_, b_, c_, d_ = [k / ad_bc for k in [a_, b_, c_, d_]]
59 x = a_ * line1_rhs + b_ * line2_rhs
60 y = c_ * line1_rhs + d_ * line2_rhs
62 return x, y
65def get_normal_points(cx, cy, cos_t, sin_t, length):
66 """
67 For a line passing through (*cx*, *cy*) and having an angle *t*, return
68 locations of the two points located along its perpendicular line at the
69 distance of *length*.
70 """
72 if length == 0.:
73 return cx, cy, cx, cy
75 cos_t1, sin_t1 = sin_t, -cos_t
76 cos_t2, sin_t2 = -sin_t, cos_t
78 x1, y1 = length * cos_t1 + cx, length * sin_t1 + cy
79 x2, y2 = length * cos_t2 + cx, length * sin_t2 + cy
81 return x1, y1, x2, y2
84# BEZIER routines
86# subdividing bezier curve
87# http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/Bezier/bezier-sub.html
90def _de_casteljau1(beta, t):
91 next_beta = beta[:-1] * (1 - t) + beta[1:] * t
92 return next_beta
95def split_de_casteljau(beta, t):
96 """
97 Split a Bézier segment defined by its control points *beta* into two
98 separate segments divided at *t* and return their control points.
99 """
100 beta = np.asarray(beta)
101 beta_list = [beta]
102 while True:
103 beta = _de_casteljau1(beta, t)
104 beta_list.append(beta)
105 if len(beta) == 1:
106 break
107 left_beta = [beta[0] for beta in beta_list]
108 right_beta = [beta[-1] for beta in reversed(beta_list)]
110 return left_beta, right_beta
113def find_bezier_t_intersecting_with_closedpath(
114 bezier_point_at_t, inside_closedpath, t0=0., t1=1., tolerance=0.01):
115 """
116 Find the intersection of the Bézier curve with a closed path.
118 The intersection point *t* is approximated by two parameters *t0*, *t1*
119 such that *t0* <= *t* <= *t1*.
121 Search starts from *t0* and *t1* and uses a simple bisecting algorithm
122 therefore one of the end points must be inside the path while the other
123 doesn't. The search stops when the distance of the points parametrized by
124 *t0* and *t1* gets smaller than the given *tolerance*.
126 Parameters
127 ----------
128 bezier_point_at_t : callable
129 A function returning x, y coordinates of the Bézier at parameter *t*.
130 It must have the signature::
132 bezier_point_at_t(t: float) -> tuple[float, float]
134 inside_closedpath : callable
135 A function returning True if a given point (x, y) is inside the
136 closed path. It must have the signature::
138 inside_closedpath(point: tuple[float, float]) -> bool
140 t0, t1 : float
141 Start parameters for the search.
143 tolerance : float
144 Maximal allowed distance between the final points.
146 Returns
147 -------
148 t0, t1 : float
149 The Bézier path parameters.
150 """
151 start = bezier_point_at_t(t0)
152 end = bezier_point_at_t(t1)
154 start_inside = inside_closedpath(start)
155 end_inside = inside_closedpath(end)
157 if start_inside == end_inside and start != end:
158 raise NonIntersectingPathException(
159 "Both points are on the same side of the closed path")
161 while True:
163 # return if the distance is smaller than the tolerance
164 if np.hypot(start[0] - end[0], start[1] - end[1]) < tolerance:
165 return t0, t1
167 # calculate the middle point
168 middle_t = 0.5 * (t0 + t1)
169 middle = bezier_point_at_t(middle_t)
170 middle_inside = inside_closedpath(middle)
172 if start_inside ^ middle_inside:
173 t1 = middle_t
174 if end == middle:
175 # Edge case where infinite loop is possible
176 # Caused by large numbers relative to tolerance
177 return t0, t1
178 end = middle
179 else:
180 t0 = middle_t
181 if start == middle:
182 # Edge case where infinite loop is possible
183 # Caused by large numbers relative to tolerance
184 return t0, t1
185 start = middle
186 start_inside = middle_inside
189class BezierSegment:
190 """
191 A d-dimensional Bézier segment.
193 Parameters
194 ----------
195 control_points : (N, d) array
196 Location of the *N* control points.
197 """
199 def __init__(self, control_points):
200 self._cpoints = np.asarray(control_points)
201 self._N, self._d = self._cpoints.shape
202 self._orders = np.arange(self._N)
203 coeff = [math.factorial(self._N - 1)
204 // (math.factorial(i) * math.factorial(self._N - 1 - i))
205 for i in range(self._N)]
206 self._px = (self._cpoints.T * coeff).T
208 def __call__(self, t):
209 """
210 Evaluate the Bézier curve at point(s) *t* in [0, 1].
212 Parameters
213 ----------
214 t : (k,) array-like
215 Points at which to evaluate the curve.
217 Returns
218 -------
219 (k, d) array
220 Value of the curve for each point in *t*.
221 """
222 t = np.asarray(t)
223 return (np.power.outer(1 - t, self._orders[::-1])
224 * np.power.outer(t, self._orders)) @ self._px
226 def point_at_t(self, t):
227 """
228 Evaluate the curve at a single point, returning a tuple of *d* floats.
229 """
230 return tuple(self(t))
232 @property
233 def control_points(self):
234 """The control points of the curve."""
235 return self._cpoints
237 @property
238 def dimension(self):
239 """The dimension of the curve."""
240 return self._d
242 @property
243 def degree(self):
244 """Degree of the polynomial. One less the number of control points."""
245 return self._N - 1
247 @property
248 def polynomial_coefficients(self):
249 r"""
250 The polynomial coefficients of the Bézier curve.
252 .. warning:: Follows opposite convention from `numpy.polyval`.
254 Returns
255 -------
256 (n+1, d) array
257 Coefficients after expanding in polynomial basis, where :math:`n`
258 is the degree of the Bézier curve and :math:`d` its dimension.
259 These are the numbers (:math:`C_j`) such that the curve can be
260 written :math:`\sum_{j=0}^n C_j t^j`.
262 Notes
263 -----
264 The coefficients are calculated as
266 .. math::
268 {n \choose j} \sum_{i=0}^j (-1)^{i+j} {j \choose i} P_i
270 where :math:`P_i` are the control points of the curve.
271 """
272 n = self.degree
273 # matplotlib uses n <= 4. overflow plausible starting around n = 15.
274 if n > 10:
275 warnings.warn("Polynomial coefficients formula unstable for high "
276 "order Bezier curves!", RuntimeWarning)
277 P = self.control_points
278 j = np.arange(n+1)[:, None]
279 i = np.arange(n+1)[None, :] # _comb is non-zero for i <= j
280 prefactor = (-1)**(i + j) * _comb(j, i) # j on axis 0, i on axis 1
281 return _comb(n, j) * prefactor @ P # j on axis 0, self.dimension on 1
283 def axis_aligned_extrema(self):
284 """
285 Return the dimension and location of the curve's interior extrema.
287 The extrema are the points along the curve where one of its partial
288 derivatives is zero.
290 Returns
291 -------
292 dims : array of int
293 Index :math:`i` of the partial derivative which is zero at each
294 interior extrema.
295 dzeros : array of float
296 Of same size as dims. The :math:`t` such that :math:`d/dx_i B(t) =
297 0`
298 """
299 n = self.degree
300 if n <= 1:
301 return np.array([]), np.array([])
302 Cj = self.polynomial_coefficients
303 dCj = np.arange(1, n+1)[:, None] * Cj[1:]
304 dims = []
305 roots = []
306 for i, pi in enumerate(dCj.T):
307 r = np.roots(pi[::-1])
308 roots.append(r)
309 dims.append(np.full_like(r, i))
310 roots = np.concatenate(roots)
311 dims = np.concatenate(dims)
312 in_range = np.isreal(roots) & (roots >= 0) & (roots <= 1)
313 return dims[in_range], np.real(roots)[in_range]
316def split_bezier_intersecting_with_closedpath(
317 bezier, inside_closedpath, tolerance=0.01):
318 """
319 Split a Bézier curve into two at the intersection with a closed path.
321 Parameters
322 ----------
323 bezier : (N, 2) array-like
324 Control points of the Bézier segment. See `.BezierSegment`.
325 inside_closedpath : callable
326 A function returning True if a given point (x, y) is inside the
327 closed path. See also `.find_bezier_t_intersecting_with_closedpath`.
328 tolerance : float
329 The tolerance for the intersection. See also
330 `.find_bezier_t_intersecting_with_closedpath`.
332 Returns
333 -------
334 left, right
335 Lists of control points for the two Bézier segments.
336 """
338 bz = BezierSegment(bezier)
339 bezier_point_at_t = bz.point_at_t
341 t0, t1 = find_bezier_t_intersecting_with_closedpath(
342 bezier_point_at_t, inside_closedpath, tolerance=tolerance)
344 _left, _right = split_de_casteljau(bezier, (t0 + t1) / 2.)
345 return _left, _right
348# matplotlib specific
351def split_path_inout(path, inside, tolerance=0.01, reorder_inout=False):
352 """
353 Divide a path into two segments at the point where ``inside(x, y)`` becomes
354 False.
355 """
356 from .path import Path
357 path_iter = path.iter_segments()
359 ctl_points, command = next(path_iter)
360 begin_inside = inside(ctl_points[-2:]) # true if begin point is inside
362 ctl_points_old = ctl_points
364 iold = 0
365 i = 1
367 for ctl_points, command in path_iter:
368 iold = i
369 i += len(ctl_points) // 2
370 if inside(ctl_points[-2:]) != begin_inside:
371 bezier_path = np.concatenate([ctl_points_old[-2:], ctl_points])
372 break
373 ctl_points_old = ctl_points
374 else:
375 raise ValueError("The path does not intersect with the patch")
377 bp = bezier_path.reshape((-1, 2))
378 left, right = split_bezier_intersecting_with_closedpath(
379 bp, inside, tolerance)
380 if len(left) == 2:
381 codes_left = [Path.LINETO]
382 codes_right = [Path.MOVETO, Path.LINETO]
383 elif len(left) == 3:
384 codes_left = [Path.CURVE3, Path.CURVE3]
385 codes_right = [Path.MOVETO, Path.CURVE3, Path.CURVE3]
386 elif len(left) == 4:
387 codes_left = [Path.CURVE4, Path.CURVE4, Path.CURVE4]
388 codes_right = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]
389 else:
390 raise AssertionError("This should never be reached")
392 verts_left = left[1:]
393 verts_right = right[:]
395 if path.codes is None:
396 path_in = Path(np.concatenate([path.vertices[:i], verts_left]))
397 path_out = Path(np.concatenate([verts_right, path.vertices[i:]]))
399 else:
400 path_in = Path(np.concatenate([path.vertices[:iold], verts_left]),
401 np.concatenate([path.codes[:iold], codes_left]))
403 path_out = Path(np.concatenate([verts_right, path.vertices[i:]]),
404 np.concatenate([codes_right, path.codes[i:]]))
406 if reorder_inout and not begin_inside:
407 path_in, path_out = path_out, path_in
409 return path_in, path_out
412def inside_circle(cx, cy, r):
413 """
414 Return a function that checks whether a point is in a circle with center
415 (*cx*, *cy*) and radius *r*.
417 The returned function has the signature::
419 f(xy: tuple[float, float]) -> bool
420 """
421 r2 = r ** 2
423 def _f(xy):
424 x, y = xy
425 return (x - cx) ** 2 + (y - cy) ** 2 < r2
426 return _f
429# quadratic Bezier lines
431def get_cos_sin(x0, y0, x1, y1):
432 dx, dy = x1 - x0, y1 - y0
433 d = (dx * dx + dy * dy) ** .5
434 # Account for divide by zero
435 if d == 0:
436 return 0.0, 0.0
437 return dx / d, dy / d
440def check_if_parallel(dx1, dy1, dx2, dy2, tolerance=1.e-5):
441 """
442 Check if two lines are parallel.
444 Parameters
445 ----------
446 dx1, dy1, dx2, dy2 : float
447 The gradients *dy*/*dx* of the two lines.
448 tolerance : float
449 The angular tolerance in radians up to which the lines are considered
450 parallel.
452 Returns
453 -------
454 is_parallel
455 - 1 if two lines are parallel in same direction.
456 - -1 if two lines are parallel in opposite direction.
457 - False otherwise.
458 """
459 theta1 = np.arctan2(dx1, dy1)
460 theta2 = np.arctan2(dx2, dy2)
461 dtheta = abs(theta1 - theta2)
462 if dtheta < tolerance:
463 return 1
464 elif abs(dtheta - np.pi) < tolerance:
465 return -1
466 else:
467 return False
470def get_parallels(bezier2, width):
471 """
472 Given the quadratic Bézier control points *bezier2*, returns
473 control points of quadratic Bézier lines roughly parallel to given
474 one separated by *width*.
475 """
477 # The parallel Bezier lines are constructed by following ways.
478 # c1 and c2 are control points representing the start and end of the
479 # Bezier line.
480 # cm is the middle point
482 c1x, c1y = bezier2[0]
483 cmx, cmy = bezier2[1]
484 c2x, c2y = bezier2[2]
486 parallel_test = check_if_parallel(c1x - cmx, c1y - cmy,
487 cmx - c2x, cmy - c2y)
489 if parallel_test == -1:
490 _api.warn_external(
491 "Lines do not intersect. A straight line is used instead.")
492 cos_t1, sin_t1 = get_cos_sin(c1x, c1y, c2x, c2y)
493 cos_t2, sin_t2 = cos_t1, sin_t1
494 else:
495 # t1 and t2 is the angle between c1 and cm, cm, c2. They are
496 # also an angle of the tangential line of the path at c1 and c2
497 cos_t1, sin_t1 = get_cos_sin(c1x, c1y, cmx, cmy)
498 cos_t2, sin_t2 = get_cos_sin(cmx, cmy, c2x, c2y)
500 # find c1_left, c1_right which are located along the lines
501 # through c1 and perpendicular to the tangential lines of the
502 # Bezier path at a distance of width. Same thing for c2_left and
503 # c2_right with respect to c2.
504 c1x_left, c1y_left, c1x_right, c1y_right = (
505 get_normal_points(c1x, c1y, cos_t1, sin_t1, width)
506 )
507 c2x_left, c2y_left, c2x_right, c2y_right = (
508 get_normal_points(c2x, c2y, cos_t2, sin_t2, width)
509 )
511 # find cm_left which is the intersecting point of a line through
512 # c1_left with angle t1 and a line through c2_left with angle
513 # t2. Same with cm_right.
514 try:
515 cmx_left, cmy_left = get_intersection(c1x_left, c1y_left, cos_t1,
516 sin_t1, c2x_left, c2y_left,
517 cos_t2, sin_t2)
518 cmx_right, cmy_right = get_intersection(c1x_right, c1y_right, cos_t1,
519 sin_t1, c2x_right, c2y_right,
520 cos_t2, sin_t2)
521 except ValueError:
522 # Special case straight lines, i.e., angle between two lines is
523 # less than the threshold used by get_intersection (we don't use
524 # check_if_parallel as the threshold is not the same).
525 cmx_left, cmy_left = (
526 0.5 * (c1x_left + c2x_left), 0.5 * (c1y_left + c2y_left)
527 )
528 cmx_right, cmy_right = (
529 0.5 * (c1x_right + c2x_right), 0.5 * (c1y_right + c2y_right)
530 )
532 # the parallel Bezier lines are created with control points of
533 # [c1_left, cm_left, c2_left] and [c1_right, cm_right, c2_right]
534 path_left = [(c1x_left, c1y_left),
535 (cmx_left, cmy_left),
536 (c2x_left, c2y_left)]
537 path_right = [(c1x_right, c1y_right),
538 (cmx_right, cmy_right),
539 (c2x_right, c2y_right)]
541 return path_left, path_right
544def find_control_points(c1x, c1y, mmx, mmy, c2x, c2y):
545 """
546 Find control points of the Bézier curve passing through (*c1x*, *c1y*),
547 (*mmx*, *mmy*), and (*c2x*, *c2y*), at parametric values 0, 0.5, and 1.
548 """
549 cmx = .5 * (4 * mmx - (c1x + c2x))
550 cmy = .5 * (4 * mmy - (c1y + c2y))
551 return [(c1x, c1y), (cmx, cmy), (c2x, c2y)]
554def make_wedged_bezier2(bezier2, width, w1=1., wm=0.5, w2=0.):
555 """
556 Being similar to `get_parallels`, returns control points of two quadratic
557 Bézier lines having a width roughly parallel to given one separated by
558 *width*.
559 """
561 # c1, cm, c2
562 c1x, c1y = bezier2[0]
563 cmx, cmy = bezier2[1]
564 c3x, c3y = bezier2[2]
566 # t1 and t2 is the angle between c1 and cm, cm, c3.
567 # They are also an angle of the tangential line of the path at c1 and c3
568 cos_t1, sin_t1 = get_cos_sin(c1x, c1y, cmx, cmy)
569 cos_t2, sin_t2 = get_cos_sin(cmx, cmy, c3x, c3y)
571 # find c1_left, c1_right which are located along the lines
572 # through c1 and perpendicular to the tangential lines of the
573 # Bezier path at a distance of width. Same thing for c3_left and
574 # c3_right with respect to c3.
575 c1x_left, c1y_left, c1x_right, c1y_right = (
576 get_normal_points(c1x, c1y, cos_t1, sin_t1, width * w1)
577 )
578 c3x_left, c3y_left, c3x_right, c3y_right = (
579 get_normal_points(c3x, c3y, cos_t2, sin_t2, width * w2)
580 )
582 # find c12, c23 and c123 which are middle points of c1-cm, cm-c3 and
583 # c12-c23
584 c12x, c12y = (c1x + cmx) * .5, (c1y + cmy) * .5
585 c23x, c23y = (cmx + c3x) * .5, (cmy + c3y) * .5
586 c123x, c123y = (c12x + c23x) * .5, (c12y + c23y) * .5
588 # tangential angle of c123 (angle between c12 and c23)
589 cos_t123, sin_t123 = get_cos_sin(c12x, c12y, c23x, c23y)
591 c123x_left, c123y_left, c123x_right, c123y_right = (
592 get_normal_points(c123x, c123y, cos_t123, sin_t123, width * wm)
593 )
595 path_left = find_control_points(c1x_left, c1y_left,
596 c123x_left, c123y_left,
597 c3x_left, c3y_left)
598 path_right = find_control_points(c1x_right, c1y_right,
599 c123x_right, c123y_right,
600 c3x_right, c3y_right)
602 return path_left, path_right