bezierTools.py revision 8ee2561bc1b45e1e0ed328c392c31137878dc0d8
105b4b4a27160e90307372f85dd99be69a9d972ffjvr"""fontTools.misc.bezierTools.py -- tools for working with bezier path segments.""" 205b4b4a27160e90307372f85dd99be69a9d972ffjvr 305b4b4a27160e90307372f85dd99be69a9d972ffjvr 405b4b4a27160e90307372f85dd99be69a9d972ffjvr__all__ = ["calcQuadraticBounds", "calcCubicBounds", "splitLine", "splitQuadratic", 505b4b4a27160e90307372f85dd99be69a9d972ffjvr "splitCubic", "solveQuadratic", "solveCubic"] 605b4b4a27160e90307372f85dd99be69a9d972ffjvr 705b4b4a27160e90307372f85dd99be69a9d972ffjvr 805b4b4a27160e90307372f85dd99be69a9d972ffjvrfrom fontTools.misc.arrayTools import calcBounds 905b4b4a27160e90307372f85dd99be69a9d972ffjvrimport Numeric 1005b4b4a27160e90307372f85dd99be69a9d972ffjvr 1105b4b4a27160e90307372f85dd99be69a9d972ffjvr 1205b4b4a27160e90307372f85dd99be69a9d972ffjvrdef calcQuadraticBounds(pt1, pt2, pt3): 1305b4b4a27160e90307372f85dd99be69a9d972ffjvr """Return the bounding rectangle for a qudratic bezier segment. 148ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr pt1 and pt3 are the "anchor" points, pt2 is the "handle". 1505b4b4a27160e90307372f85dd99be69a9d972ffjvr 168ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> calcQuadraticBounds((0, 0), (50, 100), (100, 0)) 178ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (0.0, 0.0, 100.0, 50.0) 188ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> calcQuadraticBounds((0, 0), (100, 0), (100, 100)) 198ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (0.0, 0.0, 100.0, 100.0) 208ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr """ 218ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr a, b, c = calcQuadraticParameters(pt1, pt2, pt3) 2205b4b4a27160e90307372f85dd99be69a9d972ffjvr # calc first derivative 2305b4b4a27160e90307372f85dd99be69a9d972ffjvr ax, ay = a * 2 2405b4b4a27160e90307372f85dd99be69a9d972ffjvr bx, by = b 2505b4b4a27160e90307372f85dd99be69a9d972ffjvr roots = [] 2605b4b4a27160e90307372f85dd99be69a9d972ffjvr if ax != 0: 2705b4b4a27160e90307372f85dd99be69a9d972ffjvr roots.append(-bx/ax) 2805b4b4a27160e90307372f85dd99be69a9d972ffjvr if ay != 0: 2905b4b4a27160e90307372f85dd99be69a9d972ffjvr roots.append(-by/ay) 3005b4b4a27160e90307372f85dd99be69a9d972ffjvr points = [a*t*t + b*t + c for t in roots if 0 <= t < 1] + [pt1, pt3] 3105b4b4a27160e90307372f85dd99be69a9d972ffjvr return calcBounds(points) 3205b4b4a27160e90307372f85dd99be69a9d972ffjvr 3305b4b4a27160e90307372f85dd99be69a9d972ffjvr 3405b4b4a27160e90307372f85dd99be69a9d972ffjvrdef calcCubicBounds(pt1, pt2, pt3, pt4): 3505b4b4a27160e90307372f85dd99be69a9d972ffjvr """Return the bounding rectangle for a cubic bezier segment. 368ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr pt1 and pt4 are the "anchor" points, pt2 and pt3 are the "handles". 378ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 388ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> calcCubicBounds((0, 0), (25, 100), (75, 100), (100, 0)) 398ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (0.0, 0.0, 100.0, 75.0) 408ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> calcCubicBounds((0, 0), (50, 0), (100, 50), (100, 100)) 418ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (0.0, 0.0, 100.0, 100.0) 428ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> calcCubicBounds((50, 0), (0, 100), (100, 100), (50, 0)) 438ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (35.566243270259356, 0.0, 64.433756729740679, 75.0) 448ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr """ 458ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4) 4605b4b4a27160e90307372f85dd99be69a9d972ffjvr # calc first derivative 4705b4b4a27160e90307372f85dd99be69a9d972ffjvr ax, ay = a * 3.0 4805b4b4a27160e90307372f85dd99be69a9d972ffjvr bx, by = b * 2.0 4905b4b4a27160e90307372f85dd99be69a9d972ffjvr cx, cy = c 5005b4b4a27160e90307372f85dd99be69a9d972ffjvr xRoots = [t for t in solveQuadratic(ax, bx, cx) if 0 <= t < 1] 5105b4b4a27160e90307372f85dd99be69a9d972ffjvr yRoots = [t for t in solveQuadratic(ay, by, cy) if 0 <= t < 1] 5205b4b4a27160e90307372f85dd99be69a9d972ffjvr roots = xRoots + yRoots 5305b4b4a27160e90307372f85dd99be69a9d972ffjvr 5405b4b4a27160e90307372f85dd99be69a9d972ffjvr points = [(a*t*t*t + b*t*t + c * t + d) for t in roots] + [pt1, pt4] 5505b4b4a27160e90307372f85dd99be69a9d972ffjvr return calcBounds(points) 5605b4b4a27160e90307372f85dd99be69a9d972ffjvr 5705b4b4a27160e90307372f85dd99be69a9d972ffjvr 5805b4b4a27160e90307372f85dd99be69a9d972ffjvrdef splitLine(pt1, pt2, where, isHorizontal): 5905b4b4a27160e90307372f85dd99be69a9d972ffjvr """Split the line between pt1 and pt2 at position 'where', which 6005b4b4a27160e90307372f85dd99be69a9d972ffjvr is an x coordinate if isHorizontal is False, a y coordinate if 6105b4b4a27160e90307372f85dd99be69a9d972ffjvr isHorizontal is True. Return a list of two line segments if the 6205b4b4a27160e90307372f85dd99be69a9d972ffjvr line was successfully split, or a list containing the original 638ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr line. 648ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 658ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> _tuplify(splitLine((0, 0), (100, 100), 50, True)) 668ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (((0, 0), (50.0, 50.0)), ((50.0, 50.0), (100, 100))) 678ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> _tuplify(splitLine((0, 0), (100, 100), 100, True)) 688ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (((0, 0), (100, 100)),) 698ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> _tuplify(splitLine((0, 0), (100, 100), 0, True)) 708ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (((0, 0), (0.0, 0.0)), ((0.0, 0.0), (100, 100))) 718ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> _tuplify(splitLine((0, 0), (100, 100), 0, False)) 728ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (((0, 0), (0.0, 0.0)), ((0.0, 0.0), (100, 100))) 738ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr """ 7405b4b4a27160e90307372f85dd99be69a9d972ffjvr pt1, pt2 = Numeric.array((pt1, pt2)) 7505b4b4a27160e90307372f85dd99be69a9d972ffjvr a = (pt2 - pt1) 7605b4b4a27160e90307372f85dd99be69a9d972ffjvr b = pt1 7705b4b4a27160e90307372f85dd99be69a9d972ffjvr ax = a[isHorizontal] 7805b4b4a27160e90307372f85dd99be69a9d972ffjvr if ax == 0: 7905b4b4a27160e90307372f85dd99be69a9d972ffjvr return [(pt1, pt2)] 8005b4b4a27160e90307372f85dd99be69a9d972ffjvr t = float(where - b[isHorizontal]) / ax 819524c7bdd383d6c51e0c061e0158b3d2f95ff8eajvr if 0 <= t < 1: 829524c7bdd383d6c51e0c061e0158b3d2f95ff8eajvr midPt = a * t + b 839524c7bdd383d6c51e0c061e0158b3d2f95ff8eajvr return [(pt1, midPt), (midPt, pt2)] 849524c7bdd383d6c51e0c061e0158b3d2f95ff8eajvr else: 859524c7bdd383d6c51e0c061e0158b3d2f95ff8eajvr return [(pt1, pt2)] 8605b4b4a27160e90307372f85dd99be69a9d972ffjvr 8705b4b4a27160e90307372f85dd99be69a9d972ffjvr 8805b4b4a27160e90307372f85dd99be69a9d972ffjvrdef splitQuadratic(pt1, pt2, pt3, where, isHorizontal): 8905b4b4a27160e90307372f85dd99be69a9d972ffjvr """Split the quadratic curve between pt1, pt2 and pt3 at position 'where', 9005b4b4a27160e90307372f85dd99be69a9d972ffjvr which is an x coordinate if isHorizontal is False, a y coordinate if 918ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr isHorizontal is True. Return a list of curve segments. 928ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 938ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> splitQuadratic((0, 0), (50, 100), (100, 0), 150, False) 948ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr [((0, 0), (50, 100), (100, 0))] 958ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> _tuplify(splitQuadratic((0, 0), (50, 100), (100, 0), 50, False)) 968ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (((0.0, 0.0), (25.0, 50.0), (50.0, 50.0)), ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0))) 978ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> _tuplify(splitQuadratic((0, 0), (50, 100), (100, 0), 25, False)) 988ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (((0.0, 0.0), (12.5, 25.0), (25.0, 37.5)), ((25.0, 37.5), (62.5, 75.0), (100.0, 0.0))) 998ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> _tuplify(splitQuadratic((0, 0), (50, 100), (100, 0), 25, True)) 1008ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (((0.0, 0.0), (7.3223304703363103, 14.644660940672621), (14.644660940672621, 24.999999999999996)), ((14.644660940672621, 24.999999999999996), (49.999999999999993, 75.0), (85.355339059327363, 25.000000000000025)), ((85.355339059327378, 25.0), (92.677669529663689, 14.644660940672621), (100.0, -7.1054273576010019e-15))) 1018ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> # XXX I'm not at all sure it the following behavior is desirable: 1028ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> _tuplify(splitQuadratic((0, 0), (50, 100), (100, 0), 50, True)) 1038ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (((0.0, 0.0), (25.0, 50.0), (50.0, 50.0)), ((50.0, 50.0), (50.0, 50.0), (50.0, 50.0)), ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0))) 1048ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr """ 1058ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr a, b, c = calcQuadraticParameters(pt1, pt2, pt3) 10605b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions = solveQuadratic(a[isHorizontal], b[isHorizontal], 10705b4b4a27160e90307372f85dd99be69a9d972ffjvr c[isHorizontal] - where) 10805b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions = [t for t in solutions if 0 <= t < 1] 10905b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions.sort() 11005b4b4a27160e90307372f85dd99be69a9d972ffjvr if not solutions: 11105b4b4a27160e90307372f85dd99be69a9d972ffjvr return [(pt1, pt2, pt3)] 11205b4b4a27160e90307372f85dd99be69a9d972ffjvr 11305b4b4a27160e90307372f85dd99be69a9d972ffjvr segments = [] 11405b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions.insert(0, 0.0) 11505b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions.append(1.0) 11605b4b4a27160e90307372f85dd99be69a9d972ffjvr for i in range(len(solutions) - 1): 11705b4b4a27160e90307372f85dd99be69a9d972ffjvr t1 = solutions[i] 11805b4b4a27160e90307372f85dd99be69a9d972ffjvr t2 = solutions[i+1] 11905b4b4a27160e90307372f85dd99be69a9d972ffjvr delta = (t2 - t1) 12005b4b4a27160e90307372f85dd99be69a9d972ffjvr # calc new a, b and c 12105b4b4a27160e90307372f85dd99be69a9d972ffjvr a1 = a * delta**2 12205b4b4a27160e90307372f85dd99be69a9d972ffjvr b1 = (2*a*t1 + b) * delta 12305b4b4a27160e90307372f85dd99be69a9d972ffjvr c1 = a*t1**2 + b*t1 + c 12405b4b4a27160e90307372f85dd99be69a9d972ffjvr # calc new points 12505b4b4a27160e90307372f85dd99be69a9d972ffjvr pt1 = c1 12605b4b4a27160e90307372f85dd99be69a9d972ffjvr pt2 = (b1 * 0.5) + c1 12705b4b4a27160e90307372f85dd99be69a9d972ffjvr pt3 = a1 + b1 + c1 12805b4b4a27160e90307372f85dd99be69a9d972ffjvr segments.append((pt1, pt2, pt3)) 12905b4b4a27160e90307372f85dd99be69a9d972ffjvr return segments 13005b4b4a27160e90307372f85dd99be69a9d972ffjvr 13105b4b4a27160e90307372f85dd99be69a9d972ffjvr 13205b4b4a27160e90307372f85dd99be69a9d972ffjvrdef splitCubic(pt1, pt2, pt3, pt4, where, isHorizontal): 13305b4b4a27160e90307372f85dd99be69a9d972ffjvr """Split the cubic curve between pt1, pt2, pt3 and pt4 at position 'where', 13405b4b4a27160e90307372f85dd99be69a9d972ffjvr which is an x coordinate if isHorizontal is False, a y coordinate if 13505b4b4a27160e90307372f85dd99be69a9d972ffjvr isHorizontal is True. Return a list of curve segments.""" 1368ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4) 13705b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions = solveCubic(a[isHorizontal], b[isHorizontal], c[isHorizontal], 13805b4b4a27160e90307372f85dd99be69a9d972ffjvr d[isHorizontal] - where) 13905b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions = [t for t in solutions if 0 <= t < 1] 14005b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions.sort() 14105b4b4a27160e90307372f85dd99be69a9d972ffjvr if not solutions: 14205b4b4a27160e90307372f85dd99be69a9d972ffjvr return [(pt1, pt2, pt3, pt4)] 14305b4b4a27160e90307372f85dd99be69a9d972ffjvr 14405b4b4a27160e90307372f85dd99be69a9d972ffjvr segments = [] 14505b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions.insert(0, 0.0) 14605b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions.append(1.0) 14705b4b4a27160e90307372f85dd99be69a9d972ffjvr for i in range(len(solutions) - 1): 14805b4b4a27160e90307372f85dd99be69a9d972ffjvr t1 = solutions[i] 14905b4b4a27160e90307372f85dd99be69a9d972ffjvr t2 = solutions[i+1] 15005b4b4a27160e90307372f85dd99be69a9d972ffjvr delta = (t2 - t1) 15105b4b4a27160e90307372f85dd99be69a9d972ffjvr # calc new a, b, c and d 15205b4b4a27160e90307372f85dd99be69a9d972ffjvr a1 = a * delta**3 15305b4b4a27160e90307372f85dd99be69a9d972ffjvr b1 = (3*a*t1 + b) * delta**2 15405b4b4a27160e90307372f85dd99be69a9d972ffjvr c1 = (2*b*t1 + c + 3*a*t1**2) * delta 15505b4b4a27160e90307372f85dd99be69a9d972ffjvr d1 = a*t1**3 + b*t1**2 + c*t1 + d 15605b4b4a27160e90307372f85dd99be69a9d972ffjvr # calc new points 15705b4b4a27160e90307372f85dd99be69a9d972ffjvr pt1 = d1 15805b4b4a27160e90307372f85dd99be69a9d972ffjvr pt2 = (c1 / 3.0) + d1 15905b4b4a27160e90307372f85dd99be69a9d972ffjvr pt3 = (b1 + c1) / 3.0 + pt2 16005b4b4a27160e90307372f85dd99be69a9d972ffjvr pt4 = a1 + d1 + c1 + b1 16105b4b4a27160e90307372f85dd99be69a9d972ffjvr segments.append((pt1, pt2, pt3, pt4)) 16205b4b4a27160e90307372f85dd99be69a9d972ffjvr return segments 16305b4b4a27160e90307372f85dd99be69a9d972ffjvr 16405b4b4a27160e90307372f85dd99be69a9d972ffjvr 16505b4b4a27160e90307372f85dd99be69a9d972ffjvr# 16605b4b4a27160e90307372f85dd99be69a9d972ffjvr# Equation solvers. 16705b4b4a27160e90307372f85dd99be69a9d972ffjvr# 16805b4b4a27160e90307372f85dd99be69a9d972ffjvr 16905b4b4a27160e90307372f85dd99be69a9d972ffjvrfrom math import sqrt, acos, cos, pi 17005b4b4a27160e90307372f85dd99be69a9d972ffjvr 17105b4b4a27160e90307372f85dd99be69a9d972ffjvr 17205b4b4a27160e90307372f85dd99be69a9d972ffjvrdef solveQuadratic(a, b, c, 17305b4b4a27160e90307372f85dd99be69a9d972ffjvr sqrt=sqrt): 17405b4b4a27160e90307372f85dd99be69a9d972ffjvr """Solve a quadratic equation where a, b and c are real. 17505b4b4a27160e90307372f85dd99be69a9d972ffjvr a*x*x + b*x + c = 0 176bfadfe33dbfd14d05dac7e649e715702276cc24cjvr This function returns a list of roots. Note that the returned list 177bfadfe33dbfd14d05dac7e649e715702276cc24cjvr is neither guaranteed to be sorted nor to contain unique values! 17805b4b4a27160e90307372f85dd99be69a9d972ffjvr """ 17905b4b4a27160e90307372f85dd99be69a9d972ffjvr if a == 0.0: 18005b4b4a27160e90307372f85dd99be69a9d972ffjvr if b == 0.0: 18105b4b4a27160e90307372f85dd99be69a9d972ffjvr # We have a non-equation; therefore, we have no valid solution 18205b4b4a27160e90307372f85dd99be69a9d972ffjvr roots = [] 18305b4b4a27160e90307372f85dd99be69a9d972ffjvr else: 18405b4b4a27160e90307372f85dd99be69a9d972ffjvr # We have a linear equation with 1 root. 18505b4b4a27160e90307372f85dd99be69a9d972ffjvr roots = [-c/b] 18605b4b4a27160e90307372f85dd99be69a9d972ffjvr else: 18705b4b4a27160e90307372f85dd99be69a9d972ffjvr # We have a true quadratic equation. Apply the quadratic formula to find two roots. 18805b4b4a27160e90307372f85dd99be69a9d972ffjvr DD = b*b - 4.0*a*c 18905b4b4a27160e90307372f85dd99be69a9d972ffjvr if DD >= 0.0: 190efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr rDD = sqrt(DD) 191efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr roots = [(-b+rDD)/2.0/a, (-b-rDD)/2.0/a] 19205b4b4a27160e90307372f85dd99be69a9d972ffjvr else: 19305b4b4a27160e90307372f85dd99be69a9d972ffjvr # complex roots, ignore 19405b4b4a27160e90307372f85dd99be69a9d972ffjvr roots = [] 19505b4b4a27160e90307372f85dd99be69a9d972ffjvr return roots 19605b4b4a27160e90307372f85dd99be69a9d972ffjvr 19705b4b4a27160e90307372f85dd99be69a9d972ffjvr 19805b4b4a27160e90307372f85dd99be69a9d972ffjvrdef solveCubic(a, b, c, d, 19905b4b4a27160e90307372f85dd99be69a9d972ffjvr abs=abs, pow=pow, sqrt=sqrt, cos=cos, acos=acos, pi=pi): 20005b4b4a27160e90307372f85dd99be69a9d972ffjvr """Solve a cubic equation where a, b, c and d are real. 20105b4b4a27160e90307372f85dd99be69a9d972ffjvr a*x*x*x + b*x*x + c*x + d = 0 202bfadfe33dbfd14d05dac7e649e715702276cc24cjvr This function returns a list of roots. Note that the returned list 203bfadfe33dbfd14d05dac7e649e715702276cc24cjvr is neither guaranteed to be sorted nor to contain unique values! 20405b4b4a27160e90307372f85dd99be69a9d972ffjvr """ 20505b4b4a27160e90307372f85dd99be69a9d972ffjvr # 20605b4b4a27160e90307372f85dd99be69a9d972ffjvr # adapted from: 20705b4b4a27160e90307372f85dd99be69a9d972ffjvr # CUBIC.C - Solve a cubic polynomial 20805b4b4a27160e90307372f85dd99be69a9d972ffjvr # public domain by Ross Cottrell 20905b4b4a27160e90307372f85dd99be69a9d972ffjvr # found at: http://www.strangecreations.com/library/snippets/Cubic.C 21005b4b4a27160e90307372f85dd99be69a9d972ffjvr # 21105b4b4a27160e90307372f85dd99be69a9d972ffjvr if abs(a) < 1e-6: 21205b4b4a27160e90307372f85dd99be69a9d972ffjvr # don't just test for zero; for very small values of 'a' solveCubic() 21305b4b4a27160e90307372f85dd99be69a9d972ffjvr # returns unreliable results, so we fall back to quad. 21405b4b4a27160e90307372f85dd99be69a9d972ffjvr return solveQuadratic(b, c, d) 215efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr a = float(a) 21605b4b4a27160e90307372f85dd99be69a9d972ffjvr a1 = b/a 21705b4b4a27160e90307372f85dd99be69a9d972ffjvr a2 = c/a 21805b4b4a27160e90307372f85dd99be69a9d972ffjvr a3 = d/a 21905b4b4a27160e90307372f85dd99be69a9d972ffjvr 22005b4b4a27160e90307372f85dd99be69a9d972ffjvr Q = (a1*a1 - 3.0*a2)/9.0 22105b4b4a27160e90307372f85dd99be69a9d972ffjvr R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0 22205b4b4a27160e90307372f85dd99be69a9d972ffjvr R2_Q3 = R*R - Q*Q*Q 22305b4b4a27160e90307372f85dd99be69a9d972ffjvr 224bfadfe33dbfd14d05dac7e649e715702276cc24cjvr if R2_Q3 < 0: 22505b4b4a27160e90307372f85dd99be69a9d972ffjvr theta = acos(R/sqrt(Q*Q*Q)) 226efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr rQ2 = -2.0*sqrt(Q) 227efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr x0 = rQ2*cos(theta/3.0) - a1/3.0 228efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr x1 = rQ2*cos((theta+2.0*pi)/3.0) - a1/3.0 229efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr x2 = rQ2*cos((theta+4.0*pi)/3.0) - a1/3.0 23005b4b4a27160e90307372f85dd99be69a9d972ffjvr return [x0, x1, x2] 23105b4b4a27160e90307372f85dd99be69a9d972ffjvr else: 232bfadfe33dbfd14d05dac7e649e715702276cc24cjvr if Q == 0 and R == 0: 233bfadfe33dbfd14d05dac7e649e715702276cc24cjvr x = 0 234bfadfe33dbfd14d05dac7e649e715702276cc24cjvr else: 235bfadfe33dbfd14d05dac7e649e715702276cc24cjvr x = pow(sqrt(R2_Q3)+abs(R), 1/3.0) 236bfadfe33dbfd14d05dac7e649e715702276cc24cjvr x = x + Q/x 23705b4b4a27160e90307372f85dd99be69a9d972ffjvr if R >= 0.0: 23805b4b4a27160e90307372f85dd99be69a9d972ffjvr x = -x 23905b4b4a27160e90307372f85dd99be69a9d972ffjvr x = x - a1/3.0 24005b4b4a27160e90307372f85dd99be69a9d972ffjvr return [x] 2418ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2428ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2438ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrdef calcQuadraticParameters(pt1, pt2, pt3): 2448ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr pt1, pt2, pt3 = Numeric.array((pt1, pt2, pt3)) 2458ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr c = pt1 2468ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr b = (pt2 - c) * 2.0 2478ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr a = pt3 - c - b 2488ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr return a, b, c 2498ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2508ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2518ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrdef calcCubicParameters(pt1, pt2, pt3, pt4): 2528ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr pt1, pt2, pt3, pt4 = Numeric.array((pt1, pt2, pt3, pt4)) 2538ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr d = pt1 2548ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr c = (pt2 - d) * 3.0 2558ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr b = (pt3 - pt2) * 3.0 - c 2568ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr a = pt4 - d - c - b 2578ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr return a, b, c, d 2588ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2598ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2608ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrdef _tuplify(obj): 2618ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr """ 2628ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> _tuplify([1, [2, 3], [], [[2, [3, 4]]]]) 2638ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr (1, (2, 3), (), ((2, (3, 4)),)) 2648ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr """ 2658ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr try: 2668ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr it = iter(obj) 2678ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr except TypeError: 2688ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr return obj 2698ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr else: 2708ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr return tuple([_tuplify(x) for x in it]) 2718ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2728ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2738ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrif __name__ == "__main__": 2748ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr import doctest 2758ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr doctest.testmod() 276