bezier.py

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