bezier.py 17 KB

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