bezier.py

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