1# -*- coding: utf-8 -*- 2"""fontTools.misc.bezierTools.py -- tools for working with bezier path segments. 3""" 4 5from __future__ import print_function, division, absolute_import 6from fontTools.misc.arrayTools import calcBounds 7from fontTools.misc.py23 import * 8import math 9 10 11__all__ = [ 12 "approximateCubicArcLength", 13 "approximateCubicArcLengthC", 14 "approximateQuadraticArcLength", 15 "approximateQuadraticArcLengthC", 16 "calcCubicArcLength", 17 "calcCubicArcLengthC", 18 "calcQuadraticArcLength", 19 "calcQuadraticArcLengthC", 20 "calcCubicBounds", 21 "calcQuadraticBounds", 22 "splitLine", 23 "splitQuadratic", 24 "splitCubic", 25 "splitQuadraticAtT", 26 "splitCubicAtT", 27 "solveQuadratic", 28 "solveCubic", 29] 30 31 32def calcCubicArcLength(pt1, pt2, pt3, pt4, tolerance=0.005): 33 """Return the arc length for a cubic bezier segment.""" 34 return calcCubicArcLengthC(complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4), tolerance) 35 36 37def _split_cubic_into_two(p0, p1, p2, p3): 38 mid = (p0 + 3 * (p1 + p2) + p3) * .125 39 deriv3 = (p3 + p2 - p1 - p0) * .125 40 return ((p0, (p0 + p1) * .5, mid - deriv3, mid), 41 (mid, mid + deriv3, (p2 + p3) * .5, p3)) 42 43def _calcCubicArcLengthCRecurse(mult, p0, p1, p2, p3): 44 arch = abs(p0-p3) 45 box = abs(p0-p1) + abs(p1-p2) + abs(p2-p3) 46 if arch * mult >= box: 47 return (arch + box) * .5 48 else: 49 one,two = _split_cubic_into_two(p0,p1,p2,p3) 50 return _calcCubicArcLengthCRecurse(mult, *one) + _calcCubicArcLengthCRecurse(mult, *two) 51 52def calcCubicArcLengthC(pt1, pt2, pt3, pt4, tolerance=0.005): 53 """Return the arc length for a cubic bezier segment using complex points.""" 54 mult = 1. + 1.5 * tolerance # The 1.5 is a empirical hack; no math 55 return _calcCubicArcLengthCRecurse(mult, pt1, pt2, pt3, pt4) 56 57 58epsilonDigits = 6 59epsilon = 1e-10 60 61 62def _dot(v1, v2): 63 return (v1 * v2.conjugate()).real 64 65 66def _intSecAtan(x): 67 # In : sympy.integrate(sp.sec(sp.atan(x))) 68 # Out: x*sqrt(x**2 + 1)/2 + asinh(x)/2 69 return x * math.sqrt(x**2 + 1)/2 + math.asinh(x)/2 70 71 72def calcQuadraticArcLength(pt1, pt2, pt3): 73 """Return the arc length for a qudratic bezier segment. 74 pt1 and pt3 are the "anchor" points, pt2 is the "handle". 75 76 >>> calcQuadraticArcLength((0, 0), (0, 0), (0, 0)) # empty segment 77 0.0 78 >>> calcQuadraticArcLength((0, 0), (50, 0), (80, 0)) # collinear points 79 80.0 80 >>> calcQuadraticArcLength((0, 0), (0, 50), (0, 80)) # collinear points vertical 81 80.0 82 >>> calcQuadraticArcLength((0, 0), (50, 20), (100, 40)) # collinear points 83 107.70329614269008 84 >>> calcQuadraticArcLength((0, 0), (0, 100), (100, 0)) 85 154.02976155645263 86 >>> calcQuadraticArcLength((0, 0), (0, 50), (100, 0)) 87 120.21581243984076 88 >>> calcQuadraticArcLength((0, 0), (50, -10), (80, 50)) 89 102.53273816445825 90 >>> calcQuadraticArcLength((0, 0), (40, 0), (-40, 0)) # collinear points, control point outside 91 66.66666666666667 92 >>> calcQuadraticArcLength((0, 0), (40, 0), (0, 0)) # collinear points, looping back 93 40.0 94 """ 95 return calcQuadraticArcLengthC(complex(*pt1), complex(*pt2), complex(*pt3)) 96 97 98def calcQuadraticArcLengthC(pt1, pt2, pt3): 99 """Return the arc length for a qudratic bezier segment using complex points. 100 pt1 and pt3 are the "anchor" points, pt2 is the "handle".""" 101 102 # Analytical solution to the length of a quadratic bezier. 103 # I'll explain how I arrived at this later. 104 d0 = pt2 - pt1 105 d1 = pt3 - pt2 106 d = d1 - d0 107 n = d * 1j 108 scale = abs(n) 109 if scale == 0.: 110 return abs(pt3-pt1) 111 origDist = _dot(n,d0) 112 if abs(origDist) < epsilon: 113 if _dot(d0,d1) >= 0: 114 return abs(pt3-pt1) 115 a, b = abs(d0), abs(d1) 116 return (a*a + b*b) / (a+b) 117 x0 = _dot(d,d0) / origDist 118 x1 = _dot(d,d1) / origDist 119 Len = abs(2 * (_intSecAtan(x1) - _intSecAtan(x0)) * origDist / (scale * (x1 - x0))) 120 return Len 121 122 123def approximateQuadraticArcLength(pt1, pt2, pt3): 124 # Approximate length of quadratic Bezier curve using Gauss-Legendre quadrature 125 # with n=3 points. 126 return approximateQuadraticArcLengthC(complex(*pt1), complex(*pt2), complex(*pt3)) 127 128 129def approximateQuadraticArcLengthC(pt1, pt2, pt3): 130 # Approximate length of quadratic Bezier curve using Gauss-Legendre quadrature 131 # with n=3 points for complex points. 132 # 133 # This, essentially, approximates the length-of-derivative function 134 # to be integrated with the best-matching fifth-degree polynomial 135 # approximation of it. 136 # 137 #https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Legendre_quadrature 138 139 # abs(BezierCurveC[2].diff(t).subs({t:T})) for T in sorted(.5, .5±sqrt(3/5)/2), 140 # weighted 5/18, 8/18, 5/18 respectively. 141 v0 = abs(-0.492943519233745*pt1 + 0.430331482911935*pt2 + 0.0626120363218102*pt3) 142 v1 = abs(pt3-pt1)*0.4444444444444444 143 v2 = abs(-0.0626120363218102*pt1 - 0.430331482911935*pt2 + 0.492943519233745*pt3) 144 145 return v0 + v1 + v2 146 147 148def calcQuadraticBounds(pt1, pt2, pt3): 149 """Return the bounding rectangle for a qudratic bezier segment. 150 pt1 and pt3 are the "anchor" points, pt2 is the "handle". 151 152 >>> calcQuadraticBounds((0, 0), (50, 100), (100, 0)) 153 (0, 0, 100, 50.0) 154 >>> calcQuadraticBounds((0, 0), (100, 0), (100, 100)) 155 (0.0, 0.0, 100, 100) 156 """ 157 (ax, ay), (bx, by), (cx, cy) = calcQuadraticParameters(pt1, pt2, pt3) 158 ax2 = ax*2.0 159 ay2 = ay*2.0 160 roots = [] 161 if ax2 != 0: 162 roots.append(-bx/ax2) 163 if ay2 != 0: 164 roots.append(-by/ay2) 165 points = [(ax*t*t + bx*t + cx, ay*t*t + by*t + cy) for t in roots if 0 <= t < 1] + [pt1, pt3] 166 return calcBounds(points) 167 168 169def approximateCubicArcLength(pt1, pt2, pt3, pt4): 170 """Return the approximate arc length for a cubic bezier segment. 171 pt1 and pt4 are the "anchor" points, pt2 and pt3 are the "handles". 172 173 >>> approximateCubicArcLength((0, 0), (25, 100), (75, 100), (100, 0)) 174 190.04332968932817 175 >>> approximateCubicArcLength((0, 0), (50, 0), (100, 50), (100, 100)) 176 154.8852074945903 177 >>> approximateCubicArcLength((0, 0), (50, 0), (100, 0), (150, 0)) # line; exact result should be 150. 178 149.99999999999991 179 >>> approximateCubicArcLength((0, 0), (50, 0), (100, 0), (-50, 0)) # cusp; exact result should be 150. 180 136.9267662156362 181 >>> approximateCubicArcLength((0, 0), (50, 0), (100, -50), (-50, 0)) # cusp 182 154.80848416537057 183 """ 184 # Approximate length of cubic Bezier curve using Gauss-Lobatto quadrature 185 # with n=5 points. 186 return approximateCubicArcLengthC(complex(*pt1), complex(*pt2), complex(*pt3), complex(*pt4)) 187 188 189def approximateCubicArcLengthC(pt1, pt2, pt3, pt4): 190 """Return the approximate arc length for a cubic bezier segment of complex points. 191 pt1 and pt4 are the "anchor" points, pt2 and pt3 are the "handles".""" 192 193 # Approximate length of cubic Bezier curve using Gauss-Lobatto quadrature 194 # with n=5 points for complex points. 195 # 196 # This, essentially, approximates the length-of-derivative function 197 # to be integrated with the best-matching seventh-degree polynomial 198 # approximation of it. 199 # 200 # https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Lobatto_rules 201 202 # abs(BezierCurveC[3].diff(t).subs({t:T})) for T in sorted(0, .5±(3/7)**.5/2, .5, 1), 203 # weighted 1/20, 49/180, 32/90, 49/180, 1/20 respectively. 204 v0 = abs(pt2-pt1)*.15 205 v1 = abs(-0.558983582205757*pt1 + 0.325650248872424*pt2 + 0.208983582205757*pt3 + 0.024349751127576*pt4) 206 v2 = abs(pt4-pt1+pt3-pt2)*0.26666666666666666 207 v3 = abs(-0.024349751127576*pt1 - 0.208983582205757*pt2 - 0.325650248872424*pt3 + 0.558983582205757*pt4) 208 v4 = abs(pt4-pt3)*.15 209 210 return v0 + v1 + v2 + v3 + v4 211 212 213def calcCubicBounds(pt1, pt2, pt3, pt4): 214 """Return the bounding rectangle for a cubic bezier segment. 215 pt1 and pt4 are the "anchor" points, pt2 and pt3 are the "handles". 216 217 >>> calcCubicBounds((0, 0), (25, 100), (75, 100), (100, 0)) 218 (0, 0, 100, 75.0) 219 >>> calcCubicBounds((0, 0), (50, 0), (100, 50), (100, 100)) 220 (0.0, 0.0, 100, 100) 221 >>> print("%f %f %f %f" % calcCubicBounds((50, 0), (0, 100), (100, 100), (50, 0))) 222 35.566243 0.000000 64.433757 75.000000 223 """ 224 (ax, ay), (bx, by), (cx, cy), (dx, dy) = calcCubicParameters(pt1, pt2, pt3, pt4) 225 # calc first derivative 226 ax3 = ax * 3.0 227 ay3 = ay * 3.0 228 bx2 = bx * 2.0 229 by2 = by * 2.0 230 xRoots = [t for t in solveQuadratic(ax3, bx2, cx) if 0 <= t < 1] 231 yRoots = [t for t in solveQuadratic(ay3, by2, cy) if 0 <= t < 1] 232 roots = xRoots + yRoots 233 234 points = [(ax*t*t*t + bx*t*t + cx * t + dx, ay*t*t*t + by*t*t + cy * t + dy) for t in roots] + [pt1, pt4] 235 return calcBounds(points) 236 237 238def splitLine(pt1, pt2, where, isHorizontal): 239 """Split the line between pt1 and pt2 at position 'where', which 240 is an x coordinate if isHorizontal is False, a y coordinate if 241 isHorizontal is True. Return a list of two line segments if the 242 line was successfully split, or a list containing the original 243 line. 244 245 >>> printSegments(splitLine((0, 0), (100, 100), 50, True)) 246 ((0, 0), (50, 50)) 247 ((50, 50), (100, 100)) 248 >>> printSegments(splitLine((0, 0), (100, 100), 100, True)) 249 ((0, 0), (100, 100)) 250 >>> printSegments(splitLine((0, 0), (100, 100), 0, True)) 251 ((0, 0), (0, 0)) 252 ((0, 0), (100, 100)) 253 >>> printSegments(splitLine((0, 0), (100, 100), 0, False)) 254 ((0, 0), (0, 0)) 255 ((0, 0), (100, 100)) 256 >>> printSegments(splitLine((100, 0), (0, 0), 50, False)) 257 ((100, 0), (50, 0)) 258 ((50, 0), (0, 0)) 259 >>> printSegments(splitLine((0, 100), (0, 0), 50, True)) 260 ((0, 100), (0, 50)) 261 ((0, 50), (0, 0)) 262 """ 263 pt1x, pt1y = pt1 264 pt2x, pt2y = pt2 265 266 ax = (pt2x - pt1x) 267 ay = (pt2y - pt1y) 268 269 bx = pt1x 270 by = pt1y 271 272 a = (ax, ay)[isHorizontal] 273 274 if a == 0: 275 return [(pt1, pt2)] 276 t = (where - (bx, by)[isHorizontal]) / a 277 if 0 <= t < 1: 278 midPt = ax * t + bx, ay * t + by 279 return [(pt1, midPt), (midPt, pt2)] 280 else: 281 return [(pt1, pt2)] 282 283 284def splitQuadratic(pt1, pt2, pt3, where, isHorizontal): 285 """Split the quadratic curve between pt1, pt2 and pt3 at position 'where', 286 which is an x coordinate if isHorizontal is False, a y coordinate if 287 isHorizontal is True. Return a list of curve segments. 288 289 >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 150, False)) 290 ((0, 0), (50, 100), (100, 0)) 291 >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 50, False)) 292 ((0, 0), (25, 50), (50, 50)) 293 ((50, 50), (75, 50), (100, 0)) 294 >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 25, False)) 295 ((0, 0), (12.5, 25), (25, 37.5)) 296 ((25, 37.5), (62.5, 75), (100, 0)) 297 >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 25, True)) 298 ((0, 0), (7.32233, 14.6447), (14.6447, 25)) 299 ((14.6447, 25), (50, 75), (85.3553, 25)) 300 ((85.3553, 25), (92.6777, 14.6447), (100, -7.10543e-15)) 301 >>> # XXX I'm not at all sure if the following behavior is desirable: 302 >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 50, True)) 303 ((0, 0), (25, 50), (50, 50)) 304 ((50, 50), (50, 50), (50, 50)) 305 ((50, 50), (75, 50), (100, 0)) 306 """ 307 a, b, c = calcQuadraticParameters(pt1, pt2, pt3) 308 solutions = solveQuadratic(a[isHorizontal], b[isHorizontal], 309 c[isHorizontal] - where) 310 solutions = sorted([t for t in solutions if 0 <= t < 1]) 311 if not solutions: 312 return [(pt1, pt2, pt3)] 313 return _splitQuadraticAtT(a, b, c, *solutions) 314 315 316def splitCubic(pt1, pt2, pt3, pt4, where, isHorizontal): 317 """Split the cubic curve between pt1, pt2, pt3 and pt4 at position 'where', 318 which is an x coordinate if isHorizontal is False, a y coordinate if 319 isHorizontal is True. Return a list of curve segments. 320 321 >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 150, False)) 322 ((0, 0), (25, 100), (75, 100), (100, 0)) 323 >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 50, False)) 324 ((0, 0), (12.5, 50), (31.25, 75), (50, 75)) 325 ((50, 75), (68.75, 75), (87.5, 50), (100, 0)) 326 >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 25, True)) 327 ((0, 0), (2.29379, 9.17517), (4.79804, 17.5085), (7.47414, 25)) 328 ((7.47414, 25), (31.2886, 91.6667), (68.7114, 91.6667), (92.5259, 25)) 329 ((92.5259, 25), (95.202, 17.5085), (97.7062, 9.17517), (100, 1.77636e-15)) 330 """ 331 a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4) 332 solutions = solveCubic(a[isHorizontal], b[isHorizontal], c[isHorizontal], 333 d[isHorizontal] - where) 334 solutions = sorted([t for t in solutions if 0 <= t < 1]) 335 if not solutions: 336 return [(pt1, pt2, pt3, pt4)] 337 return _splitCubicAtT(a, b, c, d, *solutions) 338 339 340def splitQuadraticAtT(pt1, pt2, pt3, *ts): 341 """Split the quadratic curve between pt1, pt2 and pt3 at one or more 342 values of t. Return a list of curve segments. 343 344 >>> printSegments(splitQuadraticAtT((0, 0), (50, 100), (100, 0), 0.5)) 345 ((0, 0), (25, 50), (50, 50)) 346 ((50, 50), (75, 50), (100, 0)) 347 >>> printSegments(splitQuadraticAtT((0, 0), (50, 100), (100, 0), 0.5, 0.75)) 348 ((0, 0), (25, 50), (50, 50)) 349 ((50, 50), (62.5, 50), (75, 37.5)) 350 ((75, 37.5), (87.5, 25), (100, 0)) 351 """ 352 a, b, c = calcQuadraticParameters(pt1, pt2, pt3) 353 return _splitQuadraticAtT(a, b, c, *ts) 354 355 356def splitCubicAtT(pt1, pt2, pt3, pt4, *ts): 357 """Split the cubic curve between pt1, pt2, pt3 and pt4 at one or more 358 values of t. Return a list of curve segments. 359 360 >>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5)) 361 ((0, 0), (12.5, 50), (31.25, 75), (50, 75)) 362 ((50, 75), (68.75, 75), (87.5, 50), (100, 0)) 363 >>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5, 0.75)) 364 ((0, 0), (12.5, 50), (31.25, 75), (50, 75)) 365 ((50, 75), (59.375, 75), (68.75, 68.75), (77.3438, 56.25)) 366 ((77.3438, 56.25), (85.9375, 43.75), (93.75, 25), (100, 0)) 367 """ 368 a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4) 369 return _splitCubicAtT(a, b, c, d, *ts) 370 371 372def _splitQuadraticAtT(a, b, c, *ts): 373 ts = list(ts) 374 segments = [] 375 ts.insert(0, 0.0) 376 ts.append(1.0) 377 ax, ay = a 378 bx, by = b 379 cx, cy = c 380 for i in range(len(ts) - 1): 381 t1 = ts[i] 382 t2 = ts[i+1] 383 delta = (t2 - t1) 384 # calc new a, b and c 385 delta_2 = delta*delta 386 a1x = ax * delta_2 387 a1y = ay * delta_2 388 b1x = (2*ax*t1 + bx) * delta 389 b1y = (2*ay*t1 + by) * delta 390 t1_2 = t1*t1 391 c1x = ax*t1_2 + bx*t1 + cx 392 c1y = ay*t1_2 + by*t1 + cy 393 394 pt1, pt2, pt3 = calcQuadraticPoints((a1x, a1y), (b1x, b1y), (c1x, c1y)) 395 segments.append((pt1, pt2, pt3)) 396 return segments 397 398 399def _splitCubicAtT(a, b, c, d, *ts): 400 ts = list(ts) 401 ts.insert(0, 0.0) 402 ts.append(1.0) 403 segments = [] 404 ax, ay = a 405 bx, by = b 406 cx, cy = c 407 dx, dy = d 408 for i in range(len(ts) - 1): 409 t1 = ts[i] 410 t2 = ts[i+1] 411 delta = (t2 - t1) 412 413 delta_2 = delta*delta 414 delta_3 = delta*delta_2 415 t1_2 = t1*t1 416 t1_3 = t1*t1_2 417 418 # calc new a, b, c and d 419 a1x = ax * delta_3 420 a1y = ay * delta_3 421 b1x = (3*ax*t1 + bx) * delta_2 422 b1y = (3*ay*t1 + by) * delta_2 423 c1x = (2*bx*t1 + cx + 3*ax*t1_2) * delta 424 c1y = (2*by*t1 + cy + 3*ay*t1_2) * delta 425 d1x = ax*t1_3 + bx*t1_2 + cx*t1 + dx 426 d1y = ay*t1_3 + by*t1_2 + cy*t1 + dy 427 pt1, pt2, pt3, pt4 = calcCubicPoints((a1x, a1y), (b1x, b1y), (c1x, c1y), (d1x, d1y)) 428 segments.append((pt1, pt2, pt3, pt4)) 429 return segments 430 431 432# 433# Equation solvers. 434# 435 436from math import sqrt, acos, cos, pi 437 438 439def solveQuadratic(a, b, c, 440 sqrt=sqrt): 441 """Solve a quadratic equation where a, b and c are real. 442 a*x*x + b*x + c = 0 443 This function returns a list of roots. Note that the returned list 444 is neither guaranteed to be sorted nor to contain unique values! 445 """ 446 if abs(a) < epsilon: 447 if abs(b) < epsilon: 448 # We have a non-equation; therefore, we have no valid solution 449 roots = [] 450 else: 451 # We have a linear equation with 1 root. 452 roots = [-c/b] 453 else: 454 # We have a true quadratic equation. Apply the quadratic formula to find two roots. 455 DD = b*b - 4.0*a*c 456 if DD >= 0.0: 457 rDD = sqrt(DD) 458 roots = [(-b+rDD)/2.0/a, (-b-rDD)/2.0/a] 459 else: 460 # complex roots, ignore 461 roots = [] 462 return roots 463 464 465def solveCubic(a, b, c, d): 466 """Solve a cubic equation where a, b, c and d are real. 467 a*x*x*x + b*x*x + c*x + d = 0 468 This function returns a list of roots. Note that the returned list 469 is neither guaranteed to be sorted nor to contain unique values! 470 471 >>> solveCubic(1, 1, -6, 0) 472 [-3.0, -0.0, 2.0] 473 >>> solveCubic(-10.0, -9.0, 48.0, -29.0) 474 [-2.9, 1.0, 1.0] 475 >>> solveCubic(-9.875, -9.0, 47.625, -28.75) 476 [-2.911392, 1.0, 1.0] 477 >>> solveCubic(1.0, -4.5, 6.75, -3.375) 478 [1.5, 1.5, 1.5] 479 >>> solveCubic(-12.0, 18.0, -9.0, 1.50023651123) 480 [0.5, 0.5, 0.5] 481 >>> solveCubic( 482 ... 9.0, 0.0, 0.0, -7.62939453125e-05 483 ... ) == [-0.0, -0.0, -0.0] 484 True 485 """ 486 # 487 # adapted from: 488 # CUBIC.C - Solve a cubic polynomial 489 # public domain by Ross Cottrell 490 # found at: http://www.strangecreations.com/library/snippets/Cubic.C 491 # 492 if abs(a) < epsilon: 493 # don't just test for zero; for very small values of 'a' solveCubic() 494 # returns unreliable results, so we fall back to quad. 495 return solveQuadratic(b, c, d) 496 a = float(a) 497 a1 = b/a 498 a2 = c/a 499 a3 = d/a 500 501 Q = (a1*a1 - 3.0*a2)/9.0 502 R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0 503 504 R2 = R*R 505 Q3 = Q*Q*Q 506 R2 = 0 if R2 < epsilon else R2 507 Q3 = 0 if abs(Q3) < epsilon else Q3 508 509 R2_Q3 = R2 - Q3 510 511 if R2 == 0. and Q3 == 0.: 512 x = round(-a1/3.0, epsilonDigits) 513 return [x, x, x] 514 elif R2_Q3 <= epsilon * .5: 515 # The epsilon * .5 above ensures that Q3 is not zero. 516 theta = acos(max(min(R/sqrt(Q3), 1.0), -1.0)) 517 rQ2 = -2.0*sqrt(Q) 518 a1_3 = a1/3.0 519 x0 = rQ2*cos(theta/3.0) - a1_3 520 x1 = rQ2*cos((theta+2.0*pi)/3.0) - a1_3 521 x2 = rQ2*cos((theta+4.0*pi)/3.0) - a1_3 522 x0, x1, x2 = sorted([x0, x1, x2]) 523 # Merge roots that are close-enough 524 if x1 - x0 < epsilon and x2 - x1 < epsilon: 525 x0 = x1 = x2 = round((x0 + x1 + x2) / 3., epsilonDigits) 526 elif x1 - x0 < epsilon: 527 x0 = x1 = round((x0 + x1) / 2., epsilonDigits) 528 x2 = round(x2, epsilonDigits) 529 elif x2 - x1 < epsilon: 530 x0 = round(x0, epsilonDigits) 531 x1 = x2 = round((x1 + x2) / 2., epsilonDigits) 532 else: 533 x0 = round(x0, epsilonDigits) 534 x1 = round(x1, epsilonDigits) 535 x2 = round(x2, epsilonDigits) 536 return [x0, x1, x2] 537 else: 538 x = pow(sqrt(R2_Q3)+abs(R), 1/3.0) 539 x = x + Q/x 540 if R >= 0.0: 541 x = -x 542 x = round(x - a1/3.0, epsilonDigits) 543 return [x] 544 545 546# 547# Conversion routines for points to parameters and vice versa 548# 549 550def calcQuadraticParameters(pt1, pt2, pt3): 551 x2, y2 = pt2 552 x3, y3 = pt3 553 cx, cy = pt1 554 bx = (x2 - cx) * 2.0 555 by = (y2 - cy) * 2.0 556 ax = x3 - cx - bx 557 ay = y3 - cy - by 558 return (ax, ay), (bx, by), (cx, cy) 559 560 561def calcCubicParameters(pt1, pt2, pt3, pt4): 562 x2, y2 = pt2 563 x3, y3 = pt3 564 x4, y4 = pt4 565 dx, dy = pt1 566 cx = (x2 -dx) * 3.0 567 cy = (y2 -dy) * 3.0 568 bx = (x3 - x2) * 3.0 - cx 569 by = (y3 - y2) * 3.0 - cy 570 ax = x4 - dx - cx - bx 571 ay = y4 - dy - cy - by 572 return (ax, ay), (bx, by), (cx, cy), (dx, dy) 573 574 575def calcQuadraticPoints(a, b, c): 576 ax, ay = a 577 bx, by = b 578 cx, cy = c 579 x1 = cx 580 y1 = cy 581 x2 = (bx * 0.5) + cx 582 y2 = (by * 0.5) + cy 583 x3 = ax + bx + cx 584 y3 = ay + by + cy 585 return (x1, y1), (x2, y2), (x3, y3) 586 587 588def calcCubicPoints(a, b, c, d): 589 ax, ay = a 590 bx, by = b 591 cx, cy = c 592 dx, dy = d 593 x1 = dx 594 y1 = dy 595 x2 = (cx / 3.0) + dx 596 y2 = (cy / 3.0) + dy 597 x3 = (bx + cx) / 3.0 + x2 598 y3 = (by + cy) / 3.0 + y2 599 x4 = ax + dx + cx + bx 600 y4 = ay + dy + cy + by 601 return (x1, y1), (x2, y2), (x3, y3), (x4, y4) 602 603 604def _segmentrepr(obj): 605 """ 606 >>> _segmentrepr([1, [2, 3], [], [[2, [3, 4], [0.1, 2.2]]]]) 607 '(1, (2, 3), (), ((2, (3, 4), (0.1, 2.2))))' 608 """ 609 try: 610 it = iter(obj) 611 except TypeError: 612 return "%g" % obj 613 else: 614 return "(%s)" % ", ".join([_segmentrepr(x) for x in it]) 615 616 617def printSegments(segments): 618 """Helper for the doctests, displaying each segment in a list of 619 segments on a single line as a tuple. 620 """ 621 for segment in segments: 622 print(_segmentrepr(segment)) 623 624if __name__ == "__main__": 625 import sys 626 import doctest 627 sys.exit(doctest.testmod().failed) 628