bezierTools.py revision befd485af59eb1d553beab340a32b02c9cb717af
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 659369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr(splitLine((0, 0), (100, 100), 50, True)) 669369c20087caa36a4b49ec47946e1af897c05037jvr '(((0, 0), (50.0, 50.0)), ((50.0, 50.0), (100, 100)))' 679369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr(splitLine((0, 0), (100, 100), 100, True)) 689369c20087caa36a4b49ec47946e1af897c05037jvr '(((0, 0), (100, 100)))' 699369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr(splitLine((0, 0), (100, 100), 0, True)) 709369c20087caa36a4b49ec47946e1af897c05037jvr '(((0, 0), (0.0, 0.0)), ((0.0, 0.0), (100, 100)))' 719369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr(splitLine((0, 0), (100, 100), 0, False)) 729369c20087caa36a4b49ec47946e1af897c05037jvr '(((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))] 959369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr(splitQuadratic((0, 0), (50, 100), (100, 0), 50, False)) 969369c20087caa36a4b49ec47946e1af897c05037jvr '(((0.0, 0.0), (25.0, 50.0), (50.0, 50.0)), ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0)))' 979369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr(splitQuadratic((0, 0), (50, 100), (100, 0), 25, False)) 989369c20087caa36a4b49ec47946e1af897c05037jvr '(((0.0, 0.0), (12.5, 25.0), (25.0, 37.5)), ((25.0, 37.5), (62.5, 75.0), (100.0, 0.0)))' 999369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr(splitQuadratic((0, 0), (50, 100), (100, 0), 25, True)) 1009369c20087caa36a4b49ec47946e1af897c05037jvr '(((0.0, 0.0), (7.32233047034, 14.6446609407), (14.6446609407, 25.0)), ((14.6446609407, 25.0), (50.0, 75.0), (85.3553390593, 25.0)), ((85.3553390593, 25.0), (92.6776695297, 14.6446609407), (100.0, -7.1054273576e-15)))' 1018ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr >>> # XXX I'm not at all sure it the following behavior is desirable: 1029369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr(splitQuadratic((0, 0), (50, 100), (100, 0), 50, True)) 1039369c20087caa36a4b49ec47946e1af897c05037jvr '(((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 124befd485af59eb1d553beab340a32b02c9cb717afjvr pt1, pt2, pt3 = calcQuadraticPoints(a1, b1, c1) 12505b4b4a27160e90307372f85dd99be69a9d972ffjvr segments.append((pt1, pt2, pt3)) 12605b4b4a27160e90307372f85dd99be69a9d972ffjvr return segments 12705b4b4a27160e90307372f85dd99be69a9d972ffjvr 12805b4b4a27160e90307372f85dd99be69a9d972ffjvr 12905b4b4a27160e90307372f85dd99be69a9d972ffjvrdef splitCubic(pt1, pt2, pt3, pt4, where, isHorizontal): 13005b4b4a27160e90307372f85dd99be69a9d972ffjvr """Split the cubic curve between pt1, pt2, pt3 and pt4 at position 'where', 13105b4b4a27160e90307372f85dd99be69a9d972ffjvr which is an x coordinate if isHorizontal is False, a y coordinate if 1329369c20087caa36a4b49ec47946e1af897c05037jvr isHorizontal is True. Return a list of curve segments. 1339369c20087caa36a4b49ec47946e1af897c05037jvr 1349369c20087caa36a4b49ec47946e1af897c05037jvr >>> splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 150, False) 1359369c20087caa36a4b49ec47946e1af897c05037jvr [((0, 0), (25, 100), (75, 100), (100, 0))] 1369369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 50, False)) 1379369c20087caa36a4b49ec47946e1af897c05037jvr '(((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0)), ((50.0, 75.0), (68.75, 75.0), (87.5, 50.0), (100.0, 0.0)))' 1389369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 25, True)) 1399369c20087caa36a4b49ec47946e1af897c05037jvr '(((0.0, 0.0), (2.2937927384, 9.17517095361), (4.79804488188, 17.5085042869), (7.47413641001, 25.0)), ((7.47413641001, 25.0), (31.2886200204, 91.6666666667), (68.7113799796, 91.6666666667), (92.52586359, 25.0)), ((92.52586359, 25.0), (95.2019551181, 17.5085042869), (97.7062072616, 9.17517095361), (100.0, 1.7763568394e-15)))' 1409369c20087caa36a4b49ec47946e1af897c05037jvr """ 1418ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4) 14205b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions = solveCubic(a[isHorizontal], b[isHorizontal], c[isHorizontal], 14305b4b4a27160e90307372f85dd99be69a9d972ffjvr d[isHorizontal] - where) 14405b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions = [t for t in solutions if 0 <= t < 1] 14505b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions.sort() 14605b4b4a27160e90307372f85dd99be69a9d972ffjvr if not solutions: 14705b4b4a27160e90307372f85dd99be69a9d972ffjvr return [(pt1, pt2, pt3, pt4)] 14805b4b4a27160e90307372f85dd99be69a9d972ffjvr 14905b4b4a27160e90307372f85dd99be69a9d972ffjvr segments = [] 15005b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions.insert(0, 0.0) 15105b4b4a27160e90307372f85dd99be69a9d972ffjvr solutions.append(1.0) 15205b4b4a27160e90307372f85dd99be69a9d972ffjvr for i in range(len(solutions) - 1): 15305b4b4a27160e90307372f85dd99be69a9d972ffjvr t1 = solutions[i] 15405b4b4a27160e90307372f85dd99be69a9d972ffjvr t2 = solutions[i+1] 15505b4b4a27160e90307372f85dd99be69a9d972ffjvr delta = (t2 - t1) 15605b4b4a27160e90307372f85dd99be69a9d972ffjvr # calc new a, b, c and d 15705b4b4a27160e90307372f85dd99be69a9d972ffjvr a1 = a * delta**3 15805b4b4a27160e90307372f85dd99be69a9d972ffjvr b1 = (3*a*t1 + b) * delta**2 15905b4b4a27160e90307372f85dd99be69a9d972ffjvr c1 = (2*b*t1 + c + 3*a*t1**2) * delta 16005b4b4a27160e90307372f85dd99be69a9d972ffjvr d1 = a*t1**3 + b*t1**2 + c*t1 + d 161befd485af59eb1d553beab340a32b02c9cb717afjvr pt1, pt2, pt3, pt4 = calcCubicPoints(a1, b1, c1, d1) 16205b4b4a27160e90307372f85dd99be69a9d972ffjvr segments.append((pt1, pt2, pt3, pt4)) 16305b4b4a27160e90307372f85dd99be69a9d972ffjvr return segments 16405b4b4a27160e90307372f85dd99be69a9d972ffjvr 16505b4b4a27160e90307372f85dd99be69a9d972ffjvr 16605b4b4a27160e90307372f85dd99be69a9d972ffjvr# 16705b4b4a27160e90307372f85dd99be69a9d972ffjvr# Equation solvers. 16805b4b4a27160e90307372f85dd99be69a9d972ffjvr# 16905b4b4a27160e90307372f85dd99be69a9d972ffjvr 17005b4b4a27160e90307372f85dd99be69a9d972ffjvrfrom math import sqrt, acos, cos, pi 17105b4b4a27160e90307372f85dd99be69a9d972ffjvr 17205b4b4a27160e90307372f85dd99be69a9d972ffjvr 17305b4b4a27160e90307372f85dd99be69a9d972ffjvrdef solveQuadratic(a, b, c, 17405b4b4a27160e90307372f85dd99be69a9d972ffjvr sqrt=sqrt): 17505b4b4a27160e90307372f85dd99be69a9d972ffjvr """Solve a quadratic equation where a, b and c are real. 17605b4b4a27160e90307372f85dd99be69a9d972ffjvr a*x*x + b*x + c = 0 177bfadfe33dbfd14d05dac7e649e715702276cc24cjvr This function returns a list of roots. Note that the returned list 178bfadfe33dbfd14d05dac7e649e715702276cc24cjvr is neither guaranteed to be sorted nor to contain unique values! 17905b4b4a27160e90307372f85dd99be69a9d972ffjvr """ 18005b4b4a27160e90307372f85dd99be69a9d972ffjvr if a == 0.0: 18105b4b4a27160e90307372f85dd99be69a9d972ffjvr if b == 0.0: 18205b4b4a27160e90307372f85dd99be69a9d972ffjvr # We have a non-equation; therefore, we have no valid solution 18305b4b4a27160e90307372f85dd99be69a9d972ffjvr roots = [] 18405b4b4a27160e90307372f85dd99be69a9d972ffjvr else: 18505b4b4a27160e90307372f85dd99be69a9d972ffjvr # We have a linear equation with 1 root. 18605b4b4a27160e90307372f85dd99be69a9d972ffjvr roots = [-c/b] 18705b4b4a27160e90307372f85dd99be69a9d972ffjvr else: 18805b4b4a27160e90307372f85dd99be69a9d972ffjvr # We have a true quadratic equation. Apply the quadratic formula to find two roots. 18905b4b4a27160e90307372f85dd99be69a9d972ffjvr DD = b*b - 4.0*a*c 19005b4b4a27160e90307372f85dd99be69a9d972ffjvr if DD >= 0.0: 191efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr rDD = sqrt(DD) 192efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr roots = [(-b+rDD)/2.0/a, (-b-rDD)/2.0/a] 19305b4b4a27160e90307372f85dd99be69a9d972ffjvr else: 19405b4b4a27160e90307372f85dd99be69a9d972ffjvr # complex roots, ignore 19505b4b4a27160e90307372f85dd99be69a9d972ffjvr roots = [] 19605b4b4a27160e90307372f85dd99be69a9d972ffjvr return roots 19705b4b4a27160e90307372f85dd99be69a9d972ffjvr 19805b4b4a27160e90307372f85dd99be69a9d972ffjvr 19905b4b4a27160e90307372f85dd99be69a9d972ffjvrdef solveCubic(a, b, c, d, 20005b4b4a27160e90307372f85dd99be69a9d972ffjvr abs=abs, pow=pow, sqrt=sqrt, cos=cos, acos=acos, pi=pi): 20105b4b4a27160e90307372f85dd99be69a9d972ffjvr """Solve a cubic equation where a, b, c and d are real. 20205b4b4a27160e90307372f85dd99be69a9d972ffjvr a*x*x*x + b*x*x + c*x + d = 0 203bfadfe33dbfd14d05dac7e649e715702276cc24cjvr This function returns a list of roots. Note that the returned list 204bfadfe33dbfd14d05dac7e649e715702276cc24cjvr is neither guaranteed to be sorted nor to contain unique values! 20505b4b4a27160e90307372f85dd99be69a9d972ffjvr """ 20605b4b4a27160e90307372f85dd99be69a9d972ffjvr # 20705b4b4a27160e90307372f85dd99be69a9d972ffjvr # adapted from: 20805b4b4a27160e90307372f85dd99be69a9d972ffjvr # CUBIC.C - Solve a cubic polynomial 20905b4b4a27160e90307372f85dd99be69a9d972ffjvr # public domain by Ross Cottrell 21005b4b4a27160e90307372f85dd99be69a9d972ffjvr # found at: http://www.strangecreations.com/library/snippets/Cubic.C 21105b4b4a27160e90307372f85dd99be69a9d972ffjvr # 21205b4b4a27160e90307372f85dd99be69a9d972ffjvr if abs(a) < 1e-6: 21305b4b4a27160e90307372f85dd99be69a9d972ffjvr # don't just test for zero; for very small values of 'a' solveCubic() 21405b4b4a27160e90307372f85dd99be69a9d972ffjvr # returns unreliable results, so we fall back to quad. 21505b4b4a27160e90307372f85dd99be69a9d972ffjvr return solveQuadratic(b, c, d) 216efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr a = float(a) 21705b4b4a27160e90307372f85dd99be69a9d972ffjvr a1 = b/a 21805b4b4a27160e90307372f85dd99be69a9d972ffjvr a2 = c/a 21905b4b4a27160e90307372f85dd99be69a9d972ffjvr a3 = d/a 22005b4b4a27160e90307372f85dd99be69a9d972ffjvr 22105b4b4a27160e90307372f85dd99be69a9d972ffjvr Q = (a1*a1 - 3.0*a2)/9.0 22205b4b4a27160e90307372f85dd99be69a9d972ffjvr R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0 22305b4b4a27160e90307372f85dd99be69a9d972ffjvr R2_Q3 = R*R - Q*Q*Q 22405b4b4a27160e90307372f85dd99be69a9d972ffjvr 225bfadfe33dbfd14d05dac7e649e715702276cc24cjvr if R2_Q3 < 0: 22605b4b4a27160e90307372f85dd99be69a9d972ffjvr theta = acos(R/sqrt(Q*Q*Q)) 227efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr rQ2 = -2.0*sqrt(Q) 228efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr x0 = rQ2*cos(theta/3.0) - a1/3.0 229efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr x1 = rQ2*cos((theta+2.0*pi)/3.0) - a1/3.0 230efaae5245c39e47d0e0e7ec46d6acca8b766151ajvr x2 = rQ2*cos((theta+4.0*pi)/3.0) - a1/3.0 23105b4b4a27160e90307372f85dd99be69a9d972ffjvr return [x0, x1, x2] 23205b4b4a27160e90307372f85dd99be69a9d972ffjvr else: 233bfadfe33dbfd14d05dac7e649e715702276cc24cjvr if Q == 0 and R == 0: 234bfadfe33dbfd14d05dac7e649e715702276cc24cjvr x = 0 235bfadfe33dbfd14d05dac7e649e715702276cc24cjvr else: 236bfadfe33dbfd14d05dac7e649e715702276cc24cjvr x = pow(sqrt(R2_Q3)+abs(R), 1/3.0) 237bfadfe33dbfd14d05dac7e649e715702276cc24cjvr x = x + Q/x 23805b4b4a27160e90307372f85dd99be69a9d972ffjvr if R >= 0.0: 23905b4b4a27160e90307372f85dd99be69a9d972ffjvr x = -x 24005b4b4a27160e90307372f85dd99be69a9d972ffjvr x = x - a1/3.0 24105b4b4a27160e90307372f85dd99be69a9d972ffjvr return [x] 2428ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2438ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2448ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrdef calcQuadraticParameters(pt1, pt2, pt3): 2458ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr pt1, pt2, pt3 = Numeric.array((pt1, pt2, pt3)) 2468ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr c = pt1 2478ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr b = (pt2 - c) * 2.0 2488ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr a = pt3 - c - b 2498ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr return a, b, c 2508ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2518ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2528ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrdef calcCubicParameters(pt1, pt2, pt3, pt4): 2538ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr pt1, pt2, pt3, pt4 = Numeric.array((pt1, pt2, pt3, pt4)) 2548ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr d = pt1 2558ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr c = (pt2 - d) * 3.0 2568ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr b = (pt3 - pt2) * 3.0 - c 2578ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr a = pt4 - d - c - b 2588ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr return a, b, c, d 2598ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2608ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 261befd485af59eb1d553beab340a32b02c9cb717afjvrdef calcQuadraticPoints(a, b, c): 262befd485af59eb1d553beab340a32b02c9cb717afjvr pt1 = c 263befd485af59eb1d553beab340a32b02c9cb717afjvr pt2 = (b * 0.5) + c 264befd485af59eb1d553beab340a32b02c9cb717afjvr pt3 = a + b + c 265befd485af59eb1d553beab340a32b02c9cb717afjvr return pt1, pt2, pt3 266befd485af59eb1d553beab340a32b02c9cb717afjvr 267befd485af59eb1d553beab340a32b02c9cb717afjvr 268befd485af59eb1d553beab340a32b02c9cb717afjvrdef calcCubicPoints(a, b, c, d): 269befd485af59eb1d553beab340a32b02c9cb717afjvr pt1 = d 270befd485af59eb1d553beab340a32b02c9cb717afjvr pt2 = (c / 3.0) + d 271befd485af59eb1d553beab340a32b02c9cb717afjvr pt3 = (b + c) / 3.0 + pt2 272befd485af59eb1d553beab340a32b02c9cb717afjvr pt4 = a + d + c + b 273befd485af59eb1d553beab340a32b02c9cb717afjvr return pt1, pt2, pt3, pt4 274befd485af59eb1d553beab340a32b02c9cb717afjvr 275befd485af59eb1d553beab340a32b02c9cb717afjvr 2769369c20087caa36a4b49ec47946e1af897c05037jvrdef _testrepr(obj): 2778ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr """ 2789369c20087caa36a4b49ec47946e1af897c05037jvr >>> _testrepr([1, [2, 3], [], [[2, [3, 4]]]]) 2799369c20087caa36a4b49ec47946e1af897c05037jvr '(1, (2, 3), (), ((2, (3, 4))))' 2808ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr """ 2818ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr try: 2828ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr it = iter(obj) 2838ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr except TypeError: 2849369c20087caa36a4b49ec47946e1af897c05037jvr return str(obj) 2858ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr else: 2869369c20087caa36a4b49ec47946e1af897c05037jvr return "(%s)" % ", ".join([_testrepr(x) for x in it]) 2878ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2888ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 2898ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrif __name__ == "__main__": 2908ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr import doctest 2918ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr doctest.testmod() 292