You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

490 lines
15 KiB

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