Coverage for /usr/lib/python3/dist-packages/matplotlib/bezier.py: 16%
224 statements
« prev ^ index » next coverage.py v7.4.4, created at 2025-06-14 15:25 +0200
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 end = middle
175 else:
176 t0 = middle_t
177 start = middle
178 start_inside = middle_inside
181class BezierSegment:
182 """
183 A d-dimensional Bézier segment.
185 Parameters
186 ----------
187 control_points : (N, d) array
188 Location of the *N* control points.
189 """
191 def __init__(self, control_points):
192 self._cpoints = np.asarray(control_points)
193 self._N, self._d = self._cpoints.shape
194 self._orders = np.arange(self._N)
195 coeff = [math.factorial(self._N - 1)
196 // (math.factorial(i) * math.factorial(self._N - 1 - i))
197 for i in range(self._N)]
198 self._px = (self._cpoints.T * coeff).T
200 def __call__(self, t):
201 """
202 Evaluate the Bézier curve at point(s) *t* in [0, 1].
204 Parameters
205 ----------
206 t : (k,) array-like
207 Points at which to evaluate the curve.
209 Returns
210 -------
211 (k, d) array
212 Value of the curve for each point in *t*.
213 """
214 t = np.asarray(t)
215 return (np.power.outer(1 - t, self._orders[::-1])
216 * np.power.outer(t, self._orders)) @ self._px
218 def point_at_t(self, t):
219 """
220 Evaluate the curve at a single point, returning a tuple of *d* floats.
221 """
222 return tuple(self(t))
224 @property
225 def control_points(self):
226 """The control points of the curve."""
227 return self._cpoints
229 @property
230 def dimension(self):
231 """The dimension of the curve."""
232 return self._d
234 @property
235 def degree(self):
236 """Degree of the polynomial. One less the number of control points."""
237 return self._N - 1
239 @property
240 def polynomial_coefficients(self):
241 r"""
242 The polynomial coefficients of the Bézier curve.
244 .. warning:: Follows opposite convention from `numpy.polyval`.
246 Returns
247 -------
248 (n+1, d) array
249 Coefficients after expanding in polynomial basis, where :math:`n`
250 is the degree of the Bézier curve and :math:`d` its dimension.
251 These are the numbers (:math:`C_j`) such that the curve can be
252 written :math:`\sum_{j=0}^n C_j t^j`.
254 Notes
255 -----
256 The coefficients are calculated as
258 .. math::
260 {n \choose j} \sum_{i=0}^j (-1)^{i+j} {j \choose i} P_i
262 where :math:`P_i` are the control points of the curve.
263 """
264 n = self.degree
265 # matplotlib uses n <= 4. overflow plausible starting around n = 15.
266 if n > 10:
267 warnings.warn("Polynomial coefficients formula unstable for high "
268 "order Bezier curves!", RuntimeWarning)
269 P = self.control_points
270 j = np.arange(n+1)[:, None]
271 i = np.arange(n+1)[None, :] # _comb is non-zero for i <= j
272 prefactor = (-1)**(i + j) * _comb(j, i) # j on axis 0, i on axis 1
273 return _comb(n, j) * prefactor @ P # j on axis 0, self.dimension on 1
275 def axis_aligned_extrema(self):
276 """
277 Return the dimension and location of the curve's interior extrema.
279 The extrema are the points along the curve where one of its partial
280 derivatives is zero.
282 Returns
283 -------
284 dims : array of int
285 Index :math:`i` of the partial derivative which is zero at each
286 interior extrema.
287 dzeros : array of float
288 Of same size as dims. The :math:`t` such that :math:`d/dx_i B(t) =
289 0`
290 """
291 n = self.degree
292 if n <= 1:
293 return np.array([]), np.array([])
294 Cj = self.polynomial_coefficients
295 dCj = np.arange(1, n+1)[:, None] * Cj[1:]
296 dims = []
297 roots = []
298 for i, pi in enumerate(dCj.T):
299 r = np.roots(pi[::-1])
300 roots.append(r)
301 dims.append(np.full_like(r, i))
302 roots = np.concatenate(roots)
303 dims = np.concatenate(dims)
304 in_range = np.isreal(roots) & (roots >= 0) & (roots <= 1)
305 return dims[in_range], np.real(roots)[in_range]
308def split_bezier_intersecting_with_closedpath(
309 bezier, inside_closedpath, tolerance=0.01):
310 """
311 Split a Bézier curve into two at the intersection with a closed path.
313 Parameters
314 ----------
315 bezier : (N, 2) array-like
316 Control points of the Bézier segment. See `.BezierSegment`.
317 inside_closedpath : callable
318 A function returning True if a given point (x, y) is inside the
319 closed path. See also `.find_bezier_t_intersecting_with_closedpath`.
320 tolerance : float
321 The tolerance for the intersection. See also
322 `.find_bezier_t_intersecting_with_closedpath`.
324 Returns
325 -------
326 left, right
327 Lists of control points for the two Bézier segments.
328 """
330 bz = BezierSegment(bezier)
331 bezier_point_at_t = bz.point_at_t
333 t0, t1 = find_bezier_t_intersecting_with_closedpath(
334 bezier_point_at_t, inside_closedpath, tolerance=tolerance)
336 _left, _right = split_de_casteljau(bezier, (t0 + t1) / 2.)
337 return _left, _right
340# matplotlib specific
343def split_path_inout(path, inside, tolerance=0.01, reorder_inout=False):
344 """
345 Divide a path into two segments at the point where ``inside(x, y)`` becomes
346 False.
347 """
348 from .path import Path
349 path_iter = path.iter_segments()
351 ctl_points, command = next(path_iter)
352 begin_inside = inside(ctl_points[-2:]) # true if begin point is inside
354 ctl_points_old = ctl_points
356 iold = 0
357 i = 1
359 for ctl_points, command in path_iter:
360 iold = i
361 i += len(ctl_points) // 2
362 if inside(ctl_points[-2:]) != begin_inside:
363 bezier_path = np.concatenate([ctl_points_old[-2:], ctl_points])
364 break
365 ctl_points_old = ctl_points
366 else:
367 raise ValueError("The path does not intersect with the patch")
369 bp = bezier_path.reshape((-1, 2))
370 left, right = split_bezier_intersecting_with_closedpath(
371 bp, inside, tolerance)
372 if len(left) == 2:
373 codes_left = [Path.LINETO]
374 codes_right = [Path.MOVETO, Path.LINETO]
375 elif len(left) == 3:
376 codes_left = [Path.CURVE3, Path.CURVE3]
377 codes_right = [Path.MOVETO, Path.CURVE3, Path.CURVE3]
378 elif len(left) == 4:
379 codes_left = [Path.CURVE4, Path.CURVE4, Path.CURVE4]
380 codes_right = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]
381 else:
382 raise AssertionError("This should never be reached")
384 verts_left = left[1:]
385 verts_right = right[:]
387 if path.codes is None:
388 path_in = Path(np.concatenate([path.vertices[:i], verts_left]))
389 path_out = Path(np.concatenate([verts_right, path.vertices[i:]]))
391 else:
392 path_in = Path(np.concatenate([path.vertices[:iold], verts_left]),
393 np.concatenate([path.codes[:iold], codes_left]))
395 path_out = Path(np.concatenate([verts_right, path.vertices[i:]]),
396 np.concatenate([codes_right, path.codes[i:]]))
398 if reorder_inout and not begin_inside:
399 path_in, path_out = path_out, path_in
401 return path_in, path_out
404def inside_circle(cx, cy, r):
405 """
406 Return a function that checks whether a point is in a circle with center
407 (*cx*, *cy*) and radius *r*.
409 The returned function has the signature::
411 f(xy: tuple[float, float]) -> bool
412 """
413 r2 = r ** 2
415 def _f(xy):
416 x, y = xy
417 return (x - cx) ** 2 + (y - cy) ** 2 < r2
418 return _f
421# quadratic Bezier lines
423def get_cos_sin(x0, y0, x1, y1):
424 dx, dy = x1 - x0, y1 - y0
425 d = (dx * dx + dy * dy) ** .5
426 # Account for divide by zero
427 if d == 0:
428 return 0.0, 0.0
429 return dx / d, dy / d
432def check_if_parallel(dx1, dy1, dx2, dy2, tolerance=1.e-5):
433 """
434 Check if two lines are parallel.
436 Parameters
437 ----------
438 dx1, dy1, dx2, dy2 : float
439 The gradients *dy*/*dx* of the two lines.
440 tolerance : float
441 The angular tolerance in radians up to which the lines are considered
442 parallel.
444 Returns
445 -------
446 is_parallel
447 - 1 if two lines are parallel in same direction.
448 - -1 if two lines are parallel in opposite direction.
449 - False otherwise.
450 """
451 theta1 = np.arctan2(dx1, dy1)
452 theta2 = np.arctan2(dx2, dy2)
453 dtheta = abs(theta1 - theta2)
454 if dtheta < tolerance:
455 return 1
456 elif abs(dtheta - np.pi) < tolerance:
457 return -1
458 else:
459 return False
462def get_parallels(bezier2, width):
463 """
464 Given the quadratic Bézier control points *bezier2*, returns
465 control points of quadratic Bézier lines roughly parallel to given
466 one separated by *width*.
467 """
469 # The parallel Bezier lines are constructed by following ways.
470 # c1 and c2 are control points representing the start and end of the
471 # Bezier line.
472 # cm is the middle point
474 c1x, c1y = bezier2[0]
475 cmx, cmy = bezier2[1]
476 c2x, c2y = bezier2[2]
478 parallel_test = check_if_parallel(c1x - cmx, c1y - cmy,
479 cmx - c2x, cmy - c2y)
481 if parallel_test == -1:
482 _api.warn_external(
483 "Lines do not intersect. A straight line is used instead.")
484 cos_t1, sin_t1 = get_cos_sin(c1x, c1y, c2x, c2y)
485 cos_t2, sin_t2 = cos_t1, sin_t1
486 else:
487 # t1 and t2 is the angle between c1 and cm, cm, c2. They are
488 # also an angle of the tangential line of the path at c1 and c2
489 cos_t1, sin_t1 = get_cos_sin(c1x, c1y, cmx, cmy)
490 cos_t2, sin_t2 = get_cos_sin(cmx, cmy, c2x, c2y)
492 # find c1_left, c1_right which are located along the lines
493 # through c1 and perpendicular to the tangential lines of the
494 # Bezier path at a distance of width. Same thing for c2_left and
495 # c2_right with respect to c2.
496 c1x_left, c1y_left, c1x_right, c1y_right = (
497 get_normal_points(c1x, c1y, cos_t1, sin_t1, width)
498 )
499 c2x_left, c2y_left, c2x_right, c2y_right = (
500 get_normal_points(c2x, c2y, cos_t2, sin_t2, width)
501 )
503 # find cm_left which is the intersecting point of a line through
504 # c1_left with angle t1 and a line through c2_left with angle
505 # t2. Same with cm_right.
506 try:
507 cmx_left, cmy_left = get_intersection(c1x_left, c1y_left, cos_t1,
508 sin_t1, c2x_left, c2y_left,
509 cos_t2, sin_t2)
510 cmx_right, cmy_right = get_intersection(c1x_right, c1y_right, cos_t1,
511 sin_t1, c2x_right, c2y_right,
512 cos_t2, sin_t2)
513 except ValueError:
514 # Special case straight lines, i.e., angle between two lines is
515 # less than the threshold used by get_intersection (we don't use
516 # check_if_parallel as the threshold is not the same).
517 cmx_left, cmy_left = (
518 0.5 * (c1x_left + c2x_left), 0.5 * (c1y_left + c2y_left)
519 )
520 cmx_right, cmy_right = (
521 0.5 * (c1x_right + c2x_right), 0.5 * (c1y_right + c2y_right)
522 )
524 # the parallel Bezier lines are created with control points of
525 # [c1_left, cm_left, c2_left] and [c1_right, cm_right, c2_right]
526 path_left = [(c1x_left, c1y_left),
527 (cmx_left, cmy_left),
528 (c2x_left, c2y_left)]
529 path_right = [(c1x_right, c1y_right),
530 (cmx_right, cmy_right),
531 (c2x_right, c2y_right)]
533 return path_left, path_right
536def find_control_points(c1x, c1y, mmx, mmy, c2x, c2y):
537 """
538 Find control points of the Bézier curve passing through (*c1x*, *c1y*),
539 (*mmx*, *mmy*), and (*c2x*, *c2y*), at parametric values 0, 0.5, and 1.
540 """
541 cmx = .5 * (4 * mmx - (c1x + c2x))
542 cmy = .5 * (4 * mmy - (c1y + c2y))
543 return [(c1x, c1y), (cmx, cmy), (c2x, c2y)]
546def make_wedged_bezier2(bezier2, width, w1=1., wm=0.5, w2=0.):
547 """
548 Being similar to `get_parallels`, returns control points of two quadratic
549 Bézier lines having a width roughly parallel to given one separated by
550 *width*.
551 """
553 # c1, cm, c2
554 c1x, c1y = bezier2[0]
555 cmx, cmy = bezier2[1]
556 c3x, c3y = bezier2[2]
558 # t1 and t2 is the angle between c1 and cm, cm, c3.
559 # They are also an angle of the tangential line of the path at c1 and c3
560 cos_t1, sin_t1 = get_cos_sin(c1x, c1y, cmx, cmy)
561 cos_t2, sin_t2 = get_cos_sin(cmx, cmy, c3x, c3y)
563 # find c1_left, c1_right which are located along the lines
564 # through c1 and perpendicular to the tangential lines of the
565 # Bezier path at a distance of width. Same thing for c3_left and
566 # c3_right with respect to c3.
567 c1x_left, c1y_left, c1x_right, c1y_right = (
568 get_normal_points(c1x, c1y, cos_t1, sin_t1, width * w1)
569 )
570 c3x_left, c3y_left, c3x_right, c3y_right = (
571 get_normal_points(c3x, c3y, cos_t2, sin_t2, width * w2)
572 )
574 # find c12, c23 and c123 which are middle points of c1-cm, cm-c3 and
575 # c12-c23
576 c12x, c12y = (c1x + cmx) * .5, (c1y + cmy) * .5
577 c23x, c23y = (cmx + c3x) * .5, (cmy + c3y) * .5
578 c123x, c123y = (c12x + c23x) * .5, (c12y + c23y) * .5
580 # tangential angle of c123 (angle between c12 and c23)
581 cos_t123, sin_t123 = get_cos_sin(c12x, c12y, c23x, c23y)
583 c123x_left, c123y_left, c123x_right, c123y_right = (
584 get_normal_points(c123x, c123y, cos_t123, sin_t123, width * wm)
585 )
587 path_left = find_control_points(c1x_left, c1y_left,
588 c123x_left, c123y_left,
589 c3x_left, c3y_left)
590 path_right = find_control_points(c1x_right, c1y_right,
591 c123x_right, c123y_right,
592 c3x_right, c3y_right)
594 return path_left, path_right