bezier.py 19 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618
  1. """
  2. A module providing some utility functions regarding Bezier path manipulation.
  3. """
  4. from functools import lru_cache
  5. import math
  6. import warnings
  7. import numpy as np
  8. import matplotlib.cbook as cbook
  9. # same algorithm as 3.8's math.comb
  10. @np.vectorize
  11. @lru_cache(maxsize=128)
  12. def _comb(n, k):
  13. if k > n:
  14. return 0
  15. k = min(k, n - k)
  16. i = np.arange(1, k + 1)
  17. return np.prod((n + 1 - i)/i).astype(int)
  18. class NonIntersectingPathException(ValueError):
  19. pass
  20. # some functions
  21. def get_intersection(cx1, cy1, cos_t1, sin_t1,
  22. cx2, cy2, cos_t2, sin_t2):
  23. """
  24. Return the intersection between the line through (*cx1*, *cy1*) at angle
  25. *t1* and the line through (*cx2*, *cy2*) at angle *t2*.
  26. """
  27. # line1 => sin_t1 * (x - cx1) - cos_t1 * (y - cy1) = 0.
  28. # line1 => sin_t1 * x + cos_t1 * y = sin_t1*cx1 - cos_t1*cy1
  29. line1_rhs = sin_t1 * cx1 - cos_t1 * cy1
  30. line2_rhs = sin_t2 * cx2 - cos_t2 * cy2
  31. # rhs matrix
  32. a, b = sin_t1, -cos_t1
  33. c, d = sin_t2, -cos_t2
  34. ad_bc = a * d - b * c
  35. if abs(ad_bc) < 1e-12:
  36. raise ValueError("Given lines do not intersect. Please verify that "
  37. "the angles are not equal or differ by 180 degrees.")
  38. # rhs_inverse
  39. a_, b_ = d, -b
  40. c_, d_ = -c, a
  41. a_, b_, c_, d_ = [k / ad_bc for k in [a_, b_, c_, d_]]
  42. x = a_ * line1_rhs + b_ * line2_rhs
  43. y = c_ * line1_rhs + d_ * line2_rhs
  44. return x, y
  45. def get_normal_points(cx, cy, cos_t, sin_t, length):
  46. """
  47. For a line passing through (*cx*, *cy*) and having an angle *t*, return
  48. locations of the two points located along its perpendicular line at the
  49. distance of *length*.
  50. """
  51. if length == 0.:
  52. return cx, cy, cx, cy
  53. cos_t1, sin_t1 = sin_t, -cos_t
  54. cos_t2, sin_t2 = -sin_t, cos_t
  55. x1, y1 = length * cos_t1 + cx, length * sin_t1 + cy
  56. x2, y2 = length * cos_t2 + cx, length * sin_t2 + cy
  57. return x1, y1, x2, y2
  58. # BEZIER routines
  59. # subdividing bezier curve
  60. # http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/Bezier/bezier-sub.html
  61. def _de_casteljau1(beta, t):
  62. next_beta = beta[:-1] * (1 - t) + beta[1:] * t
  63. return next_beta
  64. def split_de_casteljau(beta, t):
  65. """
  66. Split a Bezier segment defined by its control points *beta* into two
  67. separate segments divided at *t* and return their control points.
  68. """
  69. beta = np.asarray(beta)
  70. beta_list = [beta]
  71. while True:
  72. beta = _de_casteljau1(beta, t)
  73. beta_list.append(beta)
  74. if len(beta) == 1:
  75. break
  76. left_beta = [beta[0] for beta in beta_list]
  77. right_beta = [beta[-1] for beta in reversed(beta_list)]
  78. return left_beta, right_beta
  79. def find_bezier_t_intersecting_with_closedpath(
  80. bezier_point_at_t, inside_closedpath, t0=0., t1=1., tolerance=0.01):
  81. """
  82. Find the intersection of the Bezier curve with a closed path.
  83. The intersection point *t* is approximated by two parameters *t0*, *t1*
  84. such that *t0* <= *t* <= *t1*.
  85. Search starts from *t0* and *t1* and uses a simple bisecting algorithm
  86. therefore one of the end points must be inside the path while the other
  87. doesn't. The search stops when the distance of the points parametrized by
  88. *t0* and *t1* gets smaller than the given *tolerance*.
  89. Parameters
  90. ----------
  91. bezier_point_at_t : callable
  92. A function returning x, y coordinates of the Bezier at parameter *t*.
  93. It must have the signature::
  94. bezier_point_at_t(t: float) -> Tuple[float, float]
  95. inside_closedpath : callable
  96. A function returning True if a given point (x, y) is inside the
  97. closed path. It must have the signature::
  98. inside_closedpath(point: Tuple[float, float]) -> bool
  99. t0, t1 : float
  100. Start parameters for the search.
  101. tolerance : float
  102. Maximal allowed distance between the final points.
  103. Returns
  104. -------
  105. t0, t1 : float
  106. The Bezier path parameters.
  107. """
  108. start = bezier_point_at_t(t0)
  109. end = bezier_point_at_t(t1)
  110. start_inside = inside_closedpath(start)
  111. end_inside = inside_closedpath(end)
  112. if start_inside == end_inside and start != end:
  113. raise NonIntersectingPathException(
  114. "Both points are on the same side of the closed path")
  115. while True:
  116. # return if the distance is smaller than the tolerance
  117. if np.hypot(start[0] - end[0], start[1] - end[1]) < tolerance:
  118. return t0, t1
  119. # calculate the middle point
  120. middle_t = 0.5 * (t0 + t1)
  121. middle = bezier_point_at_t(middle_t)
  122. middle_inside = inside_closedpath(middle)
  123. if start_inside ^ middle_inside:
  124. t1 = middle_t
  125. end = middle
  126. end_inside = middle_inside
  127. else:
  128. t0 = middle_t
  129. start = middle
  130. start_inside = middle_inside
  131. class BezierSegment:
  132. """
  133. A d-dimensional Bezier segment.
  134. Parameters
  135. ----------
  136. control_points : (N, d) array
  137. Location of the *N* control points.
  138. """
  139. def __init__(self, control_points):
  140. self._cpoints = np.asarray(control_points)
  141. self._N, self._d = self._cpoints.shape
  142. self._orders = np.arange(self._N)
  143. coeff = [math.factorial(self._N - 1)
  144. // (math.factorial(i) * math.factorial(self._N - 1 - i))
  145. for i in range(self._N)]
  146. self._px = (self._cpoints.T * coeff).T
  147. def __call__(self, t):
  148. """
  149. Evaluate the Bezier curve at point(s) t in [0, 1].
  150. Parameters
  151. ----------
  152. t : float (k,), array_like
  153. Points at which to evaluate the curve.
  154. Returns
  155. -------
  156. float (k, d), array_like
  157. Value of the curve for each point in *t*.
  158. """
  159. t = np.asarray(t)
  160. return (np.power.outer(1 - t, self._orders[::-1])
  161. * np.power.outer(t, self._orders)) @ self._px
  162. def point_at_t(self, t):
  163. """Evaluate curve at a single point *t*. Returns a Tuple[float*d]."""
  164. return tuple(self(t))
  165. @property
  166. def control_points(self):
  167. """The control points of the curve."""
  168. return self._cpoints
  169. @property
  170. def dimension(self):
  171. """The dimension of the curve."""
  172. return self._d
  173. @property
  174. def degree(self):
  175. """Degree of the polynomial. One less the number of control points."""
  176. return self._N - 1
  177. @property
  178. def polynomial_coefficients(self):
  179. r"""
  180. The polynomial coefficients of the Bezier curve.
  181. .. warning:: Follows opposite convention from `numpy.polyval`.
  182. Returns
  183. -------
  184. float, (n+1, d) array_like
  185. Coefficients after expanding in polynomial basis, where :math:`n`
  186. is the degree of the bezier curve and :math:`d` its dimension.
  187. These are the numbers (:math:`C_j`) such that the curve can be
  188. written :math:`\sum_{j=0}^n C_j t^j`.
  189. Notes
  190. -----
  191. The coefficients are calculated as
  192. .. math::
  193. {n \choose j} \sum_{i=0}^j (-1)^{i+j} {j \choose i} P_i
  194. where :math:`P_i` are the control points of the curve.
  195. """
  196. n = self.degree
  197. # matplotlib uses n <= 4. overflow plausible starting around n = 15.
  198. if n > 10:
  199. warnings.warn("Polynomial coefficients formula unstable for high "
  200. "order Bezier curves!", RuntimeWarning)
  201. P = self.control_points
  202. j = np.arange(n+1)[:, None]
  203. i = np.arange(n+1)[None, :] # _comb is non-zero for i <= j
  204. prefactor = (-1)**(i + j) * _comb(j, i) # j on axis 0, i on axis 1
  205. return _comb(n, j) * prefactor @ P # j on axis 0, self.dimension on 1
  206. def axis_aligned_extrema(self):
  207. """
  208. Return the dimension and location of the curve's interior extrema.
  209. The extrema are the points along the curve where one of its partial
  210. derivatives is zero.
  211. Returns
  212. -------
  213. dims : int, array_like
  214. Index :math:`i` of the partial derivative which is zero at each
  215. interior extrema.
  216. dzeros : float, array_like
  217. Of same size as dims. The :math:`t` such that :math:`d/dx_i B(t) =
  218. 0`
  219. """
  220. n = self.degree
  221. if n <= 1:
  222. return np.array([]), np.array([])
  223. Cj = self.polynomial_coefficients
  224. dCj = np.arange(1, n+1)[:, None] * Cj[1:]
  225. dims = []
  226. roots = []
  227. for i, pi in enumerate(dCj.T):
  228. r = np.roots(pi[::-1])
  229. roots.append(r)
  230. dims.append(np.full_like(r, i))
  231. roots = np.concatenate(roots)
  232. dims = np.concatenate(dims)
  233. in_range = np.isreal(roots) & (roots >= 0) & (roots <= 1)
  234. return dims[in_range], np.real(roots)[in_range]
  235. def split_bezier_intersecting_with_closedpath(
  236. bezier, inside_closedpath, tolerance=0.01):
  237. """
  238. Split a Bezier curve into two at the intersection with a closed path.
  239. Parameters
  240. ----------
  241. bezier : array-like(N, 2)
  242. Control points of the Bezier segment. See `.BezierSegment`.
  243. inside_closedpath : callable
  244. A function returning True if a given point (x, y) is inside the
  245. closed path. See also `.find_bezier_t_intersecting_with_closedpath`.
  246. tolerance : float
  247. The tolerance for the intersection. See also
  248. `.find_bezier_t_intersecting_with_closedpath`.
  249. Returns
  250. -------
  251. left, right
  252. Lists of control points for the two Bezier segments.
  253. """
  254. bz = BezierSegment(bezier)
  255. bezier_point_at_t = bz.point_at_t
  256. t0, t1 = find_bezier_t_intersecting_with_closedpath(
  257. bezier_point_at_t, inside_closedpath, tolerance=tolerance)
  258. _left, _right = split_de_casteljau(bezier, (t0 + t1) / 2.)
  259. return _left, _right
  260. # matplotlib specific
  261. def split_path_inout(path, inside, tolerance=0.01, reorder_inout=False):
  262. """
  263. Divide a path into two segments at the point where ``inside(x, y)`` becomes
  264. False.
  265. """
  266. from .path import Path
  267. path_iter = path.iter_segments()
  268. ctl_points, command = next(path_iter)
  269. begin_inside = inside(ctl_points[-2:]) # true if begin point is inside
  270. ctl_points_old = ctl_points
  271. iold = 0
  272. i = 1
  273. for ctl_points, command in path_iter:
  274. iold = i
  275. i += len(ctl_points) // 2
  276. if inside(ctl_points[-2:]) != begin_inside:
  277. bezier_path = np.concatenate([ctl_points_old[-2:], ctl_points])
  278. break
  279. ctl_points_old = ctl_points
  280. else:
  281. raise ValueError("The path does not intersect with the patch")
  282. bp = bezier_path.reshape((-1, 2))
  283. left, right = split_bezier_intersecting_with_closedpath(
  284. bp, inside, tolerance)
  285. if len(left) == 2:
  286. codes_left = [Path.LINETO]
  287. codes_right = [Path.MOVETO, Path.LINETO]
  288. elif len(left) == 3:
  289. codes_left = [Path.CURVE3, Path.CURVE3]
  290. codes_right = [Path.MOVETO, Path.CURVE3, Path.CURVE3]
  291. elif len(left) == 4:
  292. codes_left = [Path.CURVE4, Path.CURVE4, Path.CURVE4]
  293. codes_right = [Path.MOVETO, Path.CURVE4, Path.CURVE4, Path.CURVE4]
  294. else:
  295. raise AssertionError("This should never be reached")
  296. verts_left = left[1:]
  297. verts_right = right[:]
  298. if path.codes is None:
  299. path_in = Path(np.concatenate([path.vertices[:i], verts_left]))
  300. path_out = Path(np.concatenate([verts_right, path.vertices[i:]]))
  301. else:
  302. path_in = Path(np.concatenate([path.vertices[:iold], verts_left]),
  303. np.concatenate([path.codes[:iold], codes_left]))
  304. path_out = Path(np.concatenate([verts_right, path.vertices[i:]]),
  305. np.concatenate([codes_right, path.codes[i:]]))
  306. if reorder_inout and not begin_inside:
  307. path_in, path_out = path_out, path_in
  308. return path_in, path_out
  309. def inside_circle(cx, cy, r):
  310. """
  311. Return a function that checks whether a point is in a circle with center
  312. (*cx*, *cy*) and radius *r*.
  313. The returned function has the signature::
  314. f(xy: Tuple[float, float]) -> bool
  315. """
  316. r2 = r ** 2
  317. def _f(xy):
  318. x, y = xy
  319. return (x - cx) ** 2 + (y - cy) ** 2 < r2
  320. return _f
  321. # quadratic Bezier lines
  322. def get_cos_sin(x0, y0, x1, y1):
  323. dx, dy = x1 - x0, y1 - y0
  324. d = (dx * dx + dy * dy) ** .5
  325. # Account for divide by zero
  326. if d == 0:
  327. return 0.0, 0.0
  328. return dx / d, dy / d
  329. def check_if_parallel(dx1, dy1, dx2, dy2, tolerance=1.e-5):
  330. """
  331. Check if two lines are parallel.
  332. Parameters
  333. ----------
  334. dx1, dy1, dx2, dy2 : float
  335. The gradients *dy*/*dx* of the two lines.
  336. tolerance : float
  337. The angular tolerance in radians up to which the lines are considered
  338. parallel.
  339. Returns
  340. -------
  341. is_parallel
  342. - 1 if two lines are parallel in same direction.
  343. - -1 if two lines are parallel in opposite direction.
  344. - False otherwise.
  345. """
  346. theta1 = np.arctan2(dx1, dy1)
  347. theta2 = np.arctan2(dx2, dy2)
  348. dtheta = abs(theta1 - theta2)
  349. if dtheta < tolerance:
  350. return 1
  351. elif abs(dtheta - np.pi) < tolerance:
  352. return -1
  353. else:
  354. return False
  355. def get_parallels(bezier2, width):
  356. """
  357. Given the quadratic Bezier control points *bezier2*, returns
  358. control points of quadratic Bezier lines roughly parallel to given
  359. one separated by *width*.
  360. """
  361. # The parallel Bezier lines are constructed by following ways.
  362. # c1 and c2 are control points representing the begin and end of the
  363. # Bezier line.
  364. # cm is the middle point
  365. c1x, c1y = bezier2[0]
  366. cmx, cmy = bezier2[1]
  367. c2x, c2y = bezier2[2]
  368. parallel_test = check_if_parallel(c1x - cmx, c1y - cmy,
  369. cmx - c2x, cmy - c2y)
  370. if parallel_test == -1:
  371. cbook._warn_external(
  372. "Lines do not intersect. A straight line is used instead.")
  373. cos_t1, sin_t1 = get_cos_sin(c1x, c1y, c2x, c2y)
  374. cos_t2, sin_t2 = cos_t1, sin_t1
  375. else:
  376. # t1 and t2 is the angle between c1 and cm, cm, c2. They are
  377. # also a angle of the tangential line of the path at c1 and c2
  378. cos_t1, sin_t1 = get_cos_sin(c1x, c1y, cmx, cmy)
  379. cos_t2, sin_t2 = get_cos_sin(cmx, cmy, c2x, c2y)
  380. # find c1_left, c1_right which are located along the lines
  381. # through c1 and perpendicular to the tangential lines of the
  382. # Bezier path at a distance of width. Same thing for c2_left and
  383. # c2_right with respect to c2.
  384. c1x_left, c1y_left, c1x_right, c1y_right = (
  385. get_normal_points(c1x, c1y, cos_t1, sin_t1, width)
  386. )
  387. c2x_left, c2y_left, c2x_right, c2y_right = (
  388. get_normal_points(c2x, c2y, cos_t2, sin_t2, width)
  389. )
  390. # find cm_left which is the intersecting point of a line through
  391. # c1_left with angle t1 and a line through c2_left with angle
  392. # t2. Same with cm_right.
  393. try:
  394. cmx_left, cmy_left = get_intersection(c1x_left, c1y_left, cos_t1,
  395. sin_t1, c2x_left, c2y_left,
  396. cos_t2, sin_t2)
  397. cmx_right, cmy_right = get_intersection(c1x_right, c1y_right, cos_t1,
  398. sin_t1, c2x_right, c2y_right,
  399. cos_t2, sin_t2)
  400. except ValueError:
  401. # Special case straight lines, i.e., angle between two lines is
  402. # less than the threshold used by get_intersection (we don't use
  403. # check_if_parallel as the threshold is not the same).
  404. cmx_left, cmy_left = (
  405. 0.5 * (c1x_left + c2x_left), 0.5 * (c1y_left + c2y_left)
  406. )
  407. cmx_right, cmy_right = (
  408. 0.5 * (c1x_right + c2x_right), 0.5 * (c1y_right + c2y_right)
  409. )
  410. # the parallel Bezier lines are created with control points of
  411. # [c1_left, cm_left, c2_left] and [c1_right, cm_right, c2_right]
  412. path_left = [(c1x_left, c1y_left),
  413. (cmx_left, cmy_left),
  414. (c2x_left, c2y_left)]
  415. path_right = [(c1x_right, c1y_right),
  416. (cmx_right, cmy_right),
  417. (c2x_right, c2y_right)]
  418. return path_left, path_right
  419. def find_control_points(c1x, c1y, mmx, mmy, c2x, c2y):
  420. """
  421. Find control points of the Bezier curve passing through (*c1x*, *c1y*),
  422. (*mmx*, *mmy*), and (*c2x*, *c2y*), at parametric values 0, 0.5, and 1.
  423. """
  424. cmx = .5 * (4 * mmx - (c1x + c2x))
  425. cmy = .5 * (4 * mmy - (c1y + c2y))
  426. return [(c1x, c1y), (cmx, cmy), (c2x, c2y)]
  427. def make_wedged_bezier2(bezier2, width, w1=1., wm=0.5, w2=0.):
  428. """
  429. Being similar to get_parallels, returns control points of two quadratic
  430. Bezier lines having a width roughly parallel to given one separated by
  431. *width*.
  432. """
  433. # c1, cm, c2
  434. c1x, c1y = bezier2[0]
  435. cmx, cmy = bezier2[1]
  436. c3x, c3y = bezier2[2]
  437. # t1 and t2 is the angle between c1 and cm, cm, c3.
  438. # They are also a angle of the tangential line of the path at c1 and c3
  439. cos_t1, sin_t1 = get_cos_sin(c1x, c1y, cmx, cmy)
  440. cos_t2, sin_t2 = get_cos_sin(cmx, cmy, c3x, c3y)
  441. # find c1_left, c1_right which are located along the lines
  442. # through c1 and perpendicular to the tangential lines of the
  443. # Bezier path at a distance of width. Same thing for c3_left and
  444. # c3_right with respect to c3.
  445. c1x_left, c1y_left, c1x_right, c1y_right = (
  446. get_normal_points(c1x, c1y, cos_t1, sin_t1, width * w1)
  447. )
  448. c3x_left, c3y_left, c3x_right, c3y_right = (
  449. get_normal_points(c3x, c3y, cos_t2, sin_t2, width * w2)
  450. )
  451. # find c12, c23 and c123 which are middle points of c1-cm, cm-c3 and
  452. # c12-c23
  453. c12x, c12y = (c1x + cmx) * .5, (c1y + cmy) * .5
  454. c23x, c23y = (cmx + c3x) * .5, (cmy + c3y) * .5
  455. c123x, c123y = (c12x + c23x) * .5, (c12y + c23y) * .5
  456. # tangential angle of c123 (angle between c12 and c23)
  457. cos_t123, sin_t123 = get_cos_sin(c12x, c12y, c23x, c23y)
  458. c123x_left, c123y_left, c123x_right, c123y_right = (
  459. get_normal_points(c123x, c123y, cos_t123, sin_t123, width * wm)
  460. )
  461. path_left = find_control_points(c1x_left, c1y_left,
  462. c123x_left, c123y_left,
  463. c3x_left, c3y_left)
  464. path_right = find_control_points(c1x_right, c1y_right,
  465. c123x_right, c123y_right,
  466. c3x_right, c3y_right)
  467. return path_left, path_right
  468. @cbook.deprecated(
  469. "3.3", alternative="Path.cleaned() and remove the final STOP if needed")
  470. def make_path_regular(p):
  471. """
  472. If the ``codes`` attribute of `.Path` *p* is None, return a copy of *p*
  473. with ``codes`` set to (MOVETO, LINETO, LINETO, ..., LINETO); otherwise
  474. return *p* itself.
  475. """
  476. from .path import Path
  477. c = p.codes
  478. if c is None:
  479. c = np.full(len(p.vertices), Path.LINETO, dtype=Path.code_type)
  480. c[0] = Path.MOVETO
  481. return Path(p.vertices, c)
  482. else:
  483. return p
  484. @cbook.deprecated("3.3", alternative="Path.make_compound_path()")
  485. def concatenate_paths(paths):
  486. """Concatenate a list of paths into a single path."""
  487. from .path import Path
  488. return Path.make_compound_path(*paths)