bezierTools.py revision 3ec6a258238b6068e4eef3fe579f1f5c0a06bbba
191bca4244286fb519c93fe92329da96b0e6f32eejvr"""fontTools.misc.bezierTools.py -- tools for working with bezier path segments. 291bca4244286fb519c93fe92329da96b0e6f32eejvr""" 305b4b4a27160e90307372f85dd99be69a9d972ffjvr 405b4b4a27160e90307372f85dd99be69a9d972ffjvr 510de259aec75d3ac0c34b444b2f0423fa86a7709jvr__all__ = [ 691bca4244286fb519c93fe92329da96b0e6f32eejvr "calcQuadraticBounds", 791bca4244286fb519c93fe92329da96b0e6f32eejvr "calcCubicBounds", 891bca4244286fb519c93fe92329da96b0e6f32eejvr "splitLine", 991bca4244286fb519c93fe92329da96b0e6f32eejvr "splitQuadratic", 1091bca4244286fb519c93fe92329da96b0e6f32eejvr "splitCubic", 1191bca4244286fb519c93fe92329da96b0e6f32eejvr "splitQuadraticAtT", 1291bca4244286fb519c93fe92329da96b0e6f32eejvr "splitCubicAtT", 1391bca4244286fb519c93fe92329da96b0e6f32eejvr "solveQuadratic", 1491bca4244286fb519c93fe92329da96b0e6f32eejvr "solveCubic", 1510de259aec75d3ac0c34b444b2f0423fa86a7709jvr] 1605b4b4a27160e90307372f85dd99be69a9d972ffjvr 1705b4b4a27160e90307372f85dd99be69a9d972ffjvrfrom fontTools.misc.arrayTools import calcBounds 1805b4b4a27160e90307372f85dd99be69a9d972ffjvr 19c53569efef1db3270a39ed2ecea89c54d05fc1b7jvrepsilon = 1e-12 20c53569efef1db3270a39ed2ecea89c54d05fc1b7jvr 2105b4b4a27160e90307372f85dd99be69a9d972ffjvr 2205b4b4a27160e90307372f85dd99be69a9d972ffjvrdef calcQuadraticBounds(pt1, pt2, pt3): 2391bca4244286fb519c93fe92329da96b0e6f32eejvr """Return the bounding rectangle for a qudratic bezier segment. 2491bca4244286fb519c93fe92329da96b0e6f32eejvr pt1 and pt3 are the "anchor" points, pt2 is the "handle". 2591bca4244286fb519c93fe92329da96b0e6f32eejvr 2691bca4244286fb519c93fe92329da96b0e6f32eejvr >>> calcQuadraticBounds((0, 0), (50, 100), (100, 0)) 2791bca4244286fb519c93fe92329da96b0e6f32eejvr (0, 0, 100, 50.0) 2891bca4244286fb519c93fe92329da96b0e6f32eejvr >>> calcQuadraticBounds((0, 0), (100, 0), (100, 100)) 2991bca4244286fb519c93fe92329da96b0e6f32eejvr (0.0, 0.0, 100, 100) 3091bca4244286fb519c93fe92329da96b0e6f32eejvr """ 3191bca4244286fb519c93fe92329da96b0e6f32eejvr (ax, ay), (bx, by), (cx, cy) = calcQuadraticParameters(pt1, pt2, pt3) 3291bca4244286fb519c93fe92329da96b0e6f32eejvr ax2 = ax*2.0 3391bca4244286fb519c93fe92329da96b0e6f32eejvr ay2 = ay*2.0 3491bca4244286fb519c93fe92329da96b0e6f32eejvr roots = [] 3591bca4244286fb519c93fe92329da96b0e6f32eejvr if ax2 != 0: 3691bca4244286fb519c93fe92329da96b0e6f32eejvr roots.append(-bx/ax2) 3791bca4244286fb519c93fe92329da96b0e6f32eejvr if ay2 != 0: 3891bca4244286fb519c93fe92329da96b0e6f32eejvr roots.append(-by/ay2) 3991bca4244286fb519c93fe92329da96b0e6f32eejvr points = [(ax*t*t + bx*t + cx, ay*t*t + by*t + cy) for t in roots if 0 <= t < 1] + [pt1, pt3] 4091bca4244286fb519c93fe92329da96b0e6f32eejvr return calcBounds(points) 4105b4b4a27160e90307372f85dd99be69a9d972ffjvr 4205b4b4a27160e90307372f85dd99be69a9d972ffjvr 4305b4b4a27160e90307372f85dd99be69a9d972ffjvrdef calcCubicBounds(pt1, pt2, pt3, pt4): 4491bca4244286fb519c93fe92329da96b0e6f32eejvr """Return the bounding rectangle for a cubic bezier segment. 4591bca4244286fb519c93fe92329da96b0e6f32eejvr pt1 and pt4 are the "anchor" points, pt2 and pt3 are the "handles". 4691bca4244286fb519c93fe92329da96b0e6f32eejvr 4791bca4244286fb519c93fe92329da96b0e6f32eejvr >>> calcCubicBounds((0, 0), (25, 100), (75, 100), (100, 0)) 4891bca4244286fb519c93fe92329da96b0e6f32eejvr (0, 0, 100, 75.0) 4991bca4244286fb519c93fe92329da96b0e6f32eejvr >>> calcCubicBounds((0, 0), (50, 0), (100, 50), (100, 100)) 5091bca4244286fb519c93fe92329da96b0e6f32eejvr (0.0, 0.0, 100, 100) 5191bca4244286fb519c93fe92329da96b0e6f32eejvr >>> print "%f %f %f %f" % calcCubicBounds((50, 0), (0, 100), (100, 100), (50, 0)) 5291bca4244286fb519c93fe92329da96b0e6f32eejvr 35.566243 0.000000 64.433757 75.000000 5391bca4244286fb519c93fe92329da96b0e6f32eejvr """ 5491bca4244286fb519c93fe92329da96b0e6f32eejvr (ax, ay), (bx, by), (cx, cy), (dx, dy) = calcCubicParameters(pt1, pt2, pt3, pt4) 5591bca4244286fb519c93fe92329da96b0e6f32eejvr # calc first derivative 5691bca4244286fb519c93fe92329da96b0e6f32eejvr ax3 = ax * 3.0 5791bca4244286fb519c93fe92329da96b0e6f32eejvr ay3 = ay * 3.0 5891bca4244286fb519c93fe92329da96b0e6f32eejvr bx2 = bx * 2.0 5991bca4244286fb519c93fe92329da96b0e6f32eejvr by2 = by * 2.0 6091bca4244286fb519c93fe92329da96b0e6f32eejvr xRoots = [t for t in solveQuadratic(ax3, bx2, cx) if 0 <= t < 1] 6191bca4244286fb519c93fe92329da96b0e6f32eejvr yRoots = [t for t in solveQuadratic(ay3, by2, cy) if 0 <= t < 1] 6291bca4244286fb519c93fe92329da96b0e6f32eejvr roots = xRoots + yRoots 6391bca4244286fb519c93fe92329da96b0e6f32eejvr 6491bca4244286fb519c93fe92329da96b0e6f32eejvr 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] 6591bca4244286fb519c93fe92329da96b0e6f32eejvr return calcBounds(points) 6605b4b4a27160e90307372f85dd99be69a9d972ffjvr 6705b4b4a27160e90307372f85dd99be69a9d972ffjvr 6805b4b4a27160e90307372f85dd99be69a9d972ffjvrdef splitLine(pt1, pt2, where, isHorizontal): 6991bca4244286fb519c93fe92329da96b0e6f32eejvr """Split the line between pt1 and pt2 at position 'where', which 7091bca4244286fb519c93fe92329da96b0e6f32eejvr is an x coordinate if isHorizontal is False, a y coordinate if 7191bca4244286fb519c93fe92329da96b0e6f32eejvr isHorizontal is True. Return a list of two line segments if the 7291bca4244286fb519c93fe92329da96b0e6f32eejvr line was successfully split, or a list containing the original 7391bca4244286fb519c93fe92329da96b0e6f32eejvr line. 7491bca4244286fb519c93fe92329da96b0e6f32eejvr 7591bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitLine((0, 0), (100, 100), 50, True)) 7691bca4244286fb519c93fe92329da96b0e6f32eejvr ((0, 0), (50.0, 50.0)) 7791bca4244286fb519c93fe92329da96b0e6f32eejvr ((50.0, 50.0), (100, 100)) 7891bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitLine((0, 0), (100, 100), 100, True)) 7991bca4244286fb519c93fe92329da96b0e6f32eejvr ((0, 0), (100, 100)) 8091bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitLine((0, 0), (100, 100), 0, True)) 8191bca4244286fb519c93fe92329da96b0e6f32eejvr ((0, 0), (0.0, 0.0)) 8291bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (100, 100)) 8391bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitLine((0, 0), (100, 100), 0, False)) 8491bca4244286fb519c93fe92329da96b0e6f32eejvr ((0, 0), (0.0, 0.0)) 8591bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (100, 100)) 8691bca4244286fb519c93fe92329da96b0e6f32eejvr """ 8791bca4244286fb519c93fe92329da96b0e6f32eejvr pt1x, pt1y = pt1 8891bca4244286fb519c93fe92329da96b0e6f32eejvr pt2x, pt2y = pt2 8991bca4244286fb519c93fe92329da96b0e6f32eejvr 9091bca4244286fb519c93fe92329da96b0e6f32eejvr ax = (pt2x - pt1x) 9191bca4244286fb519c93fe92329da96b0e6f32eejvr ay = (pt2y - pt1y) 9291bca4244286fb519c93fe92329da96b0e6f32eejvr 9391bca4244286fb519c93fe92329da96b0e6f32eejvr bx = pt1x 9491bca4244286fb519c93fe92329da96b0e6f32eejvr by = pt1y 9591bca4244286fb519c93fe92329da96b0e6f32eejvr 9691bca4244286fb519c93fe92329da96b0e6f32eejvr ax1 = (ax, ay)[isHorizontal] 9791bca4244286fb519c93fe92329da96b0e6f32eejvr 9891bca4244286fb519c93fe92329da96b0e6f32eejvr if ax == 0: 9991bca4244286fb519c93fe92329da96b0e6f32eejvr return [(pt1, pt2)] 10091bca4244286fb519c93fe92329da96b0e6f32eejvr 10191bca4244286fb519c93fe92329da96b0e6f32eejvr t = float(where - (bx, by)[isHorizontal]) / ax 10291bca4244286fb519c93fe92329da96b0e6f32eejvr if 0 <= t < 1: 10391bca4244286fb519c93fe92329da96b0e6f32eejvr midPt = ax * t + bx, ay * t + by 10491bca4244286fb519c93fe92329da96b0e6f32eejvr return [(pt1, midPt), (midPt, pt2)] 10591bca4244286fb519c93fe92329da96b0e6f32eejvr else: 10691bca4244286fb519c93fe92329da96b0e6f32eejvr return [(pt1, pt2)] 10705b4b4a27160e90307372f85dd99be69a9d972ffjvr 10805b4b4a27160e90307372f85dd99be69a9d972ffjvr 10905b4b4a27160e90307372f85dd99be69a9d972ffjvrdef splitQuadratic(pt1, pt2, pt3, where, isHorizontal): 11091bca4244286fb519c93fe92329da96b0e6f32eejvr """Split the quadratic curve between pt1, pt2 and pt3 at position 'where', 11191bca4244286fb519c93fe92329da96b0e6f32eejvr which is an x coordinate if isHorizontal is False, a y coordinate if 11291bca4244286fb519c93fe92329da96b0e6f32eejvr isHorizontal is True. Return a list of curve segments. 11391bca4244286fb519c93fe92329da96b0e6f32eejvr 11491bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 150, False)) 11591bca4244286fb519c93fe92329da96b0e6f32eejvr ((0, 0), (50, 100), (100, 0)) 11691bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 50, False)) 11791bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0)) 11891bca4244286fb519c93fe92329da96b0e6f32eejvr ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0)) 11991bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 25, False)) 12091bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (12.5, 25.0), (25.0, 37.5)) 12191bca4244286fb519c93fe92329da96b0e6f32eejvr ((25.0, 37.5), (62.5, 75.0), (100.0, 0.0)) 12291bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 25, True)) 12391bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (7.32233047034, 14.6446609407), (14.6446609407, 25.0)) 12491bca4244286fb519c93fe92329da96b0e6f32eejvr ((14.6446609407, 25.0), (50.0, 75.0), (85.3553390593, 25.0)) 12591bca4244286fb519c93fe92329da96b0e6f32eejvr ((85.3553390593, 25.0), (92.6776695297, 14.6446609407), (100.0, -7.1054273576e-15)) 12691bca4244286fb519c93fe92329da96b0e6f32eejvr >>> # XXX I'm not at all sure if the following behavior is desirable: 12791bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 50, True)) 12891bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0)) 12991bca4244286fb519c93fe92329da96b0e6f32eejvr ((50.0, 50.0), (50.0, 50.0), (50.0, 50.0)) 13091bca4244286fb519c93fe92329da96b0e6f32eejvr ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0)) 13191bca4244286fb519c93fe92329da96b0e6f32eejvr """ 13291bca4244286fb519c93fe92329da96b0e6f32eejvr a, b, c = calcQuadraticParameters(pt1, pt2, pt3) 13391bca4244286fb519c93fe92329da96b0e6f32eejvr solutions = solveQuadratic(a[isHorizontal], b[isHorizontal], 13491bca4244286fb519c93fe92329da96b0e6f32eejvr c[isHorizontal] - where) 135ac1b4359467ca3deab03186a15eae1d55eb35567Behdad Esfahbod solutions = sorted([t for t in solutions if 0 <= t < 1]) 13691bca4244286fb519c93fe92329da96b0e6f32eejvr if not solutions: 13791bca4244286fb519c93fe92329da96b0e6f32eejvr return [(pt1, pt2, pt3)] 13891bca4244286fb519c93fe92329da96b0e6f32eejvr return _splitQuadraticAtT(a, b, c, *solutions) 13905b4b4a27160e90307372f85dd99be69a9d972ffjvr 14005b4b4a27160e90307372f85dd99be69a9d972ffjvr 14105b4b4a27160e90307372f85dd99be69a9d972ffjvrdef splitCubic(pt1, pt2, pt3, pt4, where, isHorizontal): 14291bca4244286fb519c93fe92329da96b0e6f32eejvr """Split the cubic curve between pt1, pt2, pt3 and pt4 at position 'where', 14391bca4244286fb519c93fe92329da96b0e6f32eejvr which is an x coordinate if isHorizontal is False, a y coordinate if 14491bca4244286fb519c93fe92329da96b0e6f32eejvr isHorizontal is True. Return a list of curve segments. 14591bca4244286fb519c93fe92329da96b0e6f32eejvr 14691bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 150, False)) 14791bca4244286fb519c93fe92329da96b0e6f32eejvr ((0, 0), (25, 100), (75, 100), (100, 0)) 14891bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 50, False)) 14991bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0)) 15091bca4244286fb519c93fe92329da96b0e6f32eejvr ((50.0, 75.0), (68.75, 75.0), (87.5, 50.0), (100.0, 0.0)) 15191bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 25, True)) 15291bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (2.2937927384, 9.17517095361), (4.79804488188, 17.5085042869), (7.47413641001, 25.0)) 15391bca4244286fb519c93fe92329da96b0e6f32eejvr ((7.47413641001, 25.0), (31.2886200204, 91.6666666667), (68.7113799796, 91.6666666667), (92.52586359, 25.0)) 15491bca4244286fb519c93fe92329da96b0e6f32eejvr ((92.52586359, 25.0), (95.2019551181, 17.5085042869), (97.7062072616, 9.17517095361), (100.0, 1.7763568394e-15)) 15591bca4244286fb519c93fe92329da96b0e6f32eejvr """ 15691bca4244286fb519c93fe92329da96b0e6f32eejvr a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4) 15791bca4244286fb519c93fe92329da96b0e6f32eejvr solutions = solveCubic(a[isHorizontal], b[isHorizontal], c[isHorizontal], 15891bca4244286fb519c93fe92329da96b0e6f32eejvr d[isHorizontal] - where) 159ac1b4359467ca3deab03186a15eae1d55eb35567Behdad Esfahbod solutions = sorted([t for t in solutions if 0 <= t < 1]) 16091bca4244286fb519c93fe92329da96b0e6f32eejvr if not solutions: 16191bca4244286fb519c93fe92329da96b0e6f32eejvr return [(pt1, pt2, pt3, pt4)] 16291bca4244286fb519c93fe92329da96b0e6f32eejvr return _splitCubicAtT(a, b, c, d, *solutions) 16310de259aec75d3ac0c34b444b2f0423fa86a7709jvr 16410de259aec75d3ac0c34b444b2f0423fa86a7709jvr 16510de259aec75d3ac0c34b444b2f0423fa86a7709jvrdef splitQuadraticAtT(pt1, pt2, pt3, *ts): 16691bca4244286fb519c93fe92329da96b0e6f32eejvr """Split the quadratic curve between pt1, pt2 and pt3 at one or more 16791bca4244286fb519c93fe92329da96b0e6f32eejvr values of t. Return a list of curve segments. 16886c07d2b9ad2570823119ce453ad07275a09d94cjvr 16991bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitQuadraticAtT((0, 0), (50, 100), (100, 0), 0.5)) 17091bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0)) 17191bca4244286fb519c93fe92329da96b0e6f32eejvr ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0)) 17291bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitQuadraticAtT((0, 0), (50, 100), (100, 0), 0.5, 0.75)) 17391bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0)) 17491bca4244286fb519c93fe92329da96b0e6f32eejvr ((50.0, 50.0), (62.5, 50.0), (75.0, 37.5)) 17591bca4244286fb519c93fe92329da96b0e6f32eejvr ((75.0, 37.5), (87.5, 25.0), (100.0, 0.0)) 17691bca4244286fb519c93fe92329da96b0e6f32eejvr """ 17791bca4244286fb519c93fe92329da96b0e6f32eejvr a, b, c = calcQuadraticParameters(pt1, pt2, pt3) 17891bca4244286fb519c93fe92329da96b0e6f32eejvr return _splitQuadraticAtT(a, b, c, *ts) 17910de259aec75d3ac0c34b444b2f0423fa86a7709jvr 18010de259aec75d3ac0c34b444b2f0423fa86a7709jvr 18110de259aec75d3ac0c34b444b2f0423fa86a7709jvrdef splitCubicAtT(pt1, pt2, pt3, pt4, *ts): 18291bca4244286fb519c93fe92329da96b0e6f32eejvr """Split the cubic curve between pt1, pt2, pt3 and pt4 at one or more 18391bca4244286fb519c93fe92329da96b0e6f32eejvr values of t. Return a list of curve segments. 18486c07d2b9ad2570823119ce453ad07275a09d94cjvr 18591bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5)) 18691bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0)) 18791bca4244286fb519c93fe92329da96b0e6f32eejvr ((50.0, 75.0), (68.75, 75.0), (87.5, 50.0), (100.0, 0.0)) 18891bca4244286fb519c93fe92329da96b0e6f32eejvr >>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5, 0.75)) 18991bca4244286fb519c93fe92329da96b0e6f32eejvr ((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0)) 19091bca4244286fb519c93fe92329da96b0e6f32eejvr ((50.0, 75.0), (59.375, 75.0), (68.75, 68.75), (77.34375, 56.25)) 19191bca4244286fb519c93fe92329da96b0e6f32eejvr ((77.34375, 56.25), (85.9375, 43.75), (93.75, 25.0), (100.0, 0.0)) 19291bca4244286fb519c93fe92329da96b0e6f32eejvr """ 19391bca4244286fb519c93fe92329da96b0e6f32eejvr a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4) 19491bca4244286fb519c93fe92329da96b0e6f32eejvr return _splitCubicAtT(a, b, c, d, *ts) 19510de259aec75d3ac0c34b444b2f0423fa86a7709jvr 19610de259aec75d3ac0c34b444b2f0423fa86a7709jvr 19710de259aec75d3ac0c34b444b2f0423fa86a7709jvrdef _splitQuadraticAtT(a, b, c, *ts): 19891bca4244286fb519c93fe92329da96b0e6f32eejvr ts = list(ts) 19991bca4244286fb519c93fe92329da96b0e6f32eejvr segments = [] 20091bca4244286fb519c93fe92329da96b0e6f32eejvr ts.insert(0, 0.0) 20191bca4244286fb519c93fe92329da96b0e6f32eejvr ts.append(1.0) 20291bca4244286fb519c93fe92329da96b0e6f32eejvr ax, ay = a 20391bca4244286fb519c93fe92329da96b0e6f32eejvr bx, by = b 20491bca4244286fb519c93fe92329da96b0e6f32eejvr cx, cy = c 20591bca4244286fb519c93fe92329da96b0e6f32eejvr for i in range(len(ts) - 1): 20691bca4244286fb519c93fe92329da96b0e6f32eejvr t1 = ts[i] 20791bca4244286fb519c93fe92329da96b0e6f32eejvr t2 = ts[i+1] 20891bca4244286fb519c93fe92329da96b0e6f32eejvr delta = (t2 - t1) 20991bca4244286fb519c93fe92329da96b0e6f32eejvr # calc new a, b and c 21091bca4244286fb519c93fe92329da96b0e6f32eejvr a1x = ax * delta**2 21191bca4244286fb519c93fe92329da96b0e6f32eejvr a1y = ay * delta**2 21291bca4244286fb519c93fe92329da96b0e6f32eejvr b1x = (2*ax*t1 + bx) * delta 21391bca4244286fb519c93fe92329da96b0e6f32eejvr b1y = (2*ay*t1 + by) * delta 21491bca4244286fb519c93fe92329da96b0e6f32eejvr c1x = ax*t1**2 + bx*t1 + cx 21591bca4244286fb519c93fe92329da96b0e6f32eejvr c1y = ay*t1**2 + by*t1 + cy 21691bca4244286fb519c93fe92329da96b0e6f32eejvr 21791bca4244286fb519c93fe92329da96b0e6f32eejvr pt1, pt2, pt3 = calcQuadraticPoints((a1x, a1y), (b1x, b1y), (c1x, c1y)) 21891bca4244286fb519c93fe92329da96b0e6f32eejvr segments.append((pt1, pt2, pt3)) 21991bca4244286fb519c93fe92329da96b0e6f32eejvr return segments 22010de259aec75d3ac0c34b444b2f0423fa86a7709jvr 22110de259aec75d3ac0c34b444b2f0423fa86a7709jvr 22210de259aec75d3ac0c34b444b2f0423fa86a7709jvrdef _splitCubicAtT(a, b, c, d, *ts): 22391bca4244286fb519c93fe92329da96b0e6f32eejvr ts = list(ts) 22491bca4244286fb519c93fe92329da96b0e6f32eejvr ts.insert(0, 0.0) 22591bca4244286fb519c93fe92329da96b0e6f32eejvr ts.append(1.0) 22691bca4244286fb519c93fe92329da96b0e6f32eejvr segments = [] 22791bca4244286fb519c93fe92329da96b0e6f32eejvr ax, ay = a 22891bca4244286fb519c93fe92329da96b0e6f32eejvr bx, by = b 22991bca4244286fb519c93fe92329da96b0e6f32eejvr cx, cy = c 23091bca4244286fb519c93fe92329da96b0e6f32eejvr dx, dy = d 23191bca4244286fb519c93fe92329da96b0e6f32eejvr for i in range(len(ts) - 1): 23291bca4244286fb519c93fe92329da96b0e6f32eejvr t1 = ts[i] 23391bca4244286fb519c93fe92329da96b0e6f32eejvr t2 = ts[i+1] 23491bca4244286fb519c93fe92329da96b0e6f32eejvr delta = (t2 - t1) 23591bca4244286fb519c93fe92329da96b0e6f32eejvr # calc new a, b, c and d 23691bca4244286fb519c93fe92329da96b0e6f32eejvr a1x = ax * delta**3 23791bca4244286fb519c93fe92329da96b0e6f32eejvr a1y = ay * delta**3 23891bca4244286fb519c93fe92329da96b0e6f32eejvr b1x = (3*ax*t1 + bx) * delta**2 23991bca4244286fb519c93fe92329da96b0e6f32eejvr b1y = (3*ay*t1 + by) * delta**2 24091bca4244286fb519c93fe92329da96b0e6f32eejvr c1x = (2*bx*t1 + cx + 3*ax*t1**2) * delta 24191bca4244286fb519c93fe92329da96b0e6f32eejvr c1y = (2*by*t1 + cy + 3*ay*t1**2) * delta 24291bca4244286fb519c93fe92329da96b0e6f32eejvr d1x = ax*t1**3 + bx*t1**2 + cx*t1 + dx 24391bca4244286fb519c93fe92329da96b0e6f32eejvr d1y = ay*t1**3 + by*t1**2 + cy*t1 + dy 24491bca4244286fb519c93fe92329da96b0e6f32eejvr pt1, pt2, pt3, pt4 = calcCubicPoints((a1x, a1y), (b1x, b1y), (c1x, c1y), (d1x, d1y)) 24591bca4244286fb519c93fe92329da96b0e6f32eejvr segments.append((pt1, pt2, pt3, pt4)) 24691bca4244286fb519c93fe92329da96b0e6f32eejvr return segments 24705b4b4a27160e90307372f85dd99be69a9d972ffjvr 24805b4b4a27160e90307372f85dd99be69a9d972ffjvr 24905b4b4a27160e90307372f85dd99be69a9d972ffjvr# 25005b4b4a27160e90307372f85dd99be69a9d972ffjvr# Equation solvers. 25105b4b4a27160e90307372f85dd99be69a9d972ffjvr# 25205b4b4a27160e90307372f85dd99be69a9d972ffjvr 25305b4b4a27160e90307372f85dd99be69a9d972ffjvrfrom math import sqrt, acos, cos, pi 25405b4b4a27160e90307372f85dd99be69a9d972ffjvr 25505b4b4a27160e90307372f85dd99be69a9d972ffjvr 25605b4b4a27160e90307372f85dd99be69a9d972ffjvrdef solveQuadratic(a, b, c, 25791bca4244286fb519c93fe92329da96b0e6f32eejvr sqrt=sqrt): 25891bca4244286fb519c93fe92329da96b0e6f32eejvr """Solve a quadratic equation where a, b and c are real. 25991bca4244286fb519c93fe92329da96b0e6f32eejvr a*x*x + b*x + c = 0 26091bca4244286fb519c93fe92329da96b0e6f32eejvr This function returns a list of roots. Note that the returned list 26191bca4244286fb519c93fe92329da96b0e6f32eejvr is neither guaranteed to be sorted nor to contain unique values! 26291bca4244286fb519c93fe92329da96b0e6f32eejvr """ 26391bca4244286fb519c93fe92329da96b0e6f32eejvr if abs(a) < epsilon: 26491bca4244286fb519c93fe92329da96b0e6f32eejvr if abs(b) < epsilon: 26591bca4244286fb519c93fe92329da96b0e6f32eejvr # We have a non-equation; therefore, we have no valid solution 26691bca4244286fb519c93fe92329da96b0e6f32eejvr roots = [] 26791bca4244286fb519c93fe92329da96b0e6f32eejvr else: 26891bca4244286fb519c93fe92329da96b0e6f32eejvr # We have a linear equation with 1 root. 26991bca4244286fb519c93fe92329da96b0e6f32eejvr roots = [-c/b] 27091bca4244286fb519c93fe92329da96b0e6f32eejvr else: 27191bca4244286fb519c93fe92329da96b0e6f32eejvr # We have a true quadratic equation. Apply the quadratic formula to find two roots. 27291bca4244286fb519c93fe92329da96b0e6f32eejvr DD = b*b - 4.0*a*c 27391bca4244286fb519c93fe92329da96b0e6f32eejvr if DD >= 0.0: 27491bca4244286fb519c93fe92329da96b0e6f32eejvr rDD = sqrt(DD) 27591bca4244286fb519c93fe92329da96b0e6f32eejvr roots = [(-b+rDD)/2.0/a, (-b-rDD)/2.0/a] 27691bca4244286fb519c93fe92329da96b0e6f32eejvr else: 27791bca4244286fb519c93fe92329da96b0e6f32eejvr # complex roots, ignore 27891bca4244286fb519c93fe92329da96b0e6f32eejvr roots = [] 27991bca4244286fb519c93fe92329da96b0e6f32eejvr return roots 28005b4b4a27160e90307372f85dd99be69a9d972ffjvr 28105b4b4a27160e90307372f85dd99be69a9d972ffjvr 28205b4b4a27160e90307372f85dd99be69a9d972ffjvrdef solveCubic(a, b, c, d, 28391bca4244286fb519c93fe92329da96b0e6f32eejvr abs=abs, pow=pow, sqrt=sqrt, cos=cos, acos=acos, pi=pi): 28491bca4244286fb519c93fe92329da96b0e6f32eejvr """Solve a cubic equation where a, b, c and d are real. 28591bca4244286fb519c93fe92329da96b0e6f32eejvr a*x*x*x + b*x*x + c*x + d = 0 28691bca4244286fb519c93fe92329da96b0e6f32eejvr This function returns a list of roots. Note that the returned list 28791bca4244286fb519c93fe92329da96b0e6f32eejvr is neither guaranteed to be sorted nor to contain unique values! 28891bca4244286fb519c93fe92329da96b0e6f32eejvr """ 28991bca4244286fb519c93fe92329da96b0e6f32eejvr # 29091bca4244286fb519c93fe92329da96b0e6f32eejvr # adapted from: 29191bca4244286fb519c93fe92329da96b0e6f32eejvr # CUBIC.C - Solve a cubic polynomial 29291bca4244286fb519c93fe92329da96b0e6f32eejvr # public domain by Ross Cottrell 29391bca4244286fb519c93fe92329da96b0e6f32eejvr # found at: http://www.strangecreations.com/library/snippets/Cubic.C 29491bca4244286fb519c93fe92329da96b0e6f32eejvr # 29591bca4244286fb519c93fe92329da96b0e6f32eejvr if abs(a) < epsilon: 29691bca4244286fb519c93fe92329da96b0e6f32eejvr # don't just test for zero; for very small values of 'a' solveCubic() 29791bca4244286fb519c93fe92329da96b0e6f32eejvr # returns unreliable results, so we fall back to quad. 29891bca4244286fb519c93fe92329da96b0e6f32eejvr return solveQuadratic(b, c, d) 29991bca4244286fb519c93fe92329da96b0e6f32eejvr a = float(a) 30091bca4244286fb519c93fe92329da96b0e6f32eejvr a1 = b/a 30191bca4244286fb519c93fe92329da96b0e6f32eejvr a2 = c/a 30291bca4244286fb519c93fe92329da96b0e6f32eejvr a3 = d/a 30391bca4244286fb519c93fe92329da96b0e6f32eejvr 30491bca4244286fb519c93fe92329da96b0e6f32eejvr Q = (a1*a1 - 3.0*a2)/9.0 30591bca4244286fb519c93fe92329da96b0e6f32eejvr R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0 30691bca4244286fb519c93fe92329da96b0e6f32eejvr R2_Q3 = R*R - Q*Q*Q 30791bca4244286fb519c93fe92329da96b0e6f32eejvr 30891bca4244286fb519c93fe92329da96b0e6f32eejvr if R2_Q3 < 0: 30991bca4244286fb519c93fe92329da96b0e6f32eejvr theta = acos(R/sqrt(Q*Q*Q)) 31091bca4244286fb519c93fe92329da96b0e6f32eejvr rQ2 = -2.0*sqrt(Q) 31191bca4244286fb519c93fe92329da96b0e6f32eejvr x0 = rQ2*cos(theta/3.0) - a1/3.0 31291bca4244286fb519c93fe92329da96b0e6f32eejvr x1 = rQ2*cos((theta+2.0*pi)/3.0) - a1/3.0 31391bca4244286fb519c93fe92329da96b0e6f32eejvr x2 = rQ2*cos((theta+4.0*pi)/3.0) - a1/3.0 31491bca4244286fb519c93fe92329da96b0e6f32eejvr return [x0, x1, x2] 31591bca4244286fb519c93fe92329da96b0e6f32eejvr else: 31691bca4244286fb519c93fe92329da96b0e6f32eejvr if Q == 0 and R == 0: 31791bca4244286fb519c93fe92329da96b0e6f32eejvr x = 0 31891bca4244286fb519c93fe92329da96b0e6f32eejvr else: 31991bca4244286fb519c93fe92329da96b0e6f32eejvr x = pow(sqrt(R2_Q3)+abs(R), 1/3.0) 32091bca4244286fb519c93fe92329da96b0e6f32eejvr x = x + Q/x 32191bca4244286fb519c93fe92329da96b0e6f32eejvr if R >= 0.0: 32291bca4244286fb519c93fe92329da96b0e6f32eejvr x = -x 32391bca4244286fb519c93fe92329da96b0e6f32eejvr x = x - a1/3.0 32491bca4244286fb519c93fe92329da96b0e6f32eejvr return [x] 3258ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 3268ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 32782fdf153acf0c82926ad5f49b473b9dba71ba3d3jvr# 32882fdf153acf0c82926ad5f49b473b9dba71ba3d3jvr# Conversion routines for points to parameters and vice versa 32982fdf153acf0c82926ad5f49b473b9dba71ba3d3jvr# 33082fdf153acf0c82926ad5f49b473b9dba71ba3d3jvr 3318ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrdef calcQuadraticParameters(pt1, pt2, pt3): 33291bca4244286fb519c93fe92329da96b0e6f32eejvr x2, y2 = pt2 33391bca4244286fb519c93fe92329da96b0e6f32eejvr x3, y3 = pt3 33491bca4244286fb519c93fe92329da96b0e6f32eejvr cx, cy = pt1 33591bca4244286fb519c93fe92329da96b0e6f32eejvr bx = (x2 - cx) * 2.0 33691bca4244286fb519c93fe92329da96b0e6f32eejvr by = (y2 - cy) * 2.0 33791bca4244286fb519c93fe92329da96b0e6f32eejvr ax = x3 - cx - bx 33891bca4244286fb519c93fe92329da96b0e6f32eejvr ay = y3 - cy - by 33991bca4244286fb519c93fe92329da96b0e6f32eejvr return (ax, ay), (bx, by), (cx, cy) 3408ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 3418ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 3428ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrdef calcCubicParameters(pt1, pt2, pt3, pt4): 34391bca4244286fb519c93fe92329da96b0e6f32eejvr x2, y2 = pt2 34491bca4244286fb519c93fe92329da96b0e6f32eejvr x3, y3 = pt3 34591bca4244286fb519c93fe92329da96b0e6f32eejvr x4, y4 = pt4 34691bca4244286fb519c93fe92329da96b0e6f32eejvr dx, dy = pt1 34791bca4244286fb519c93fe92329da96b0e6f32eejvr cx = (x2 -dx) * 3.0 34891bca4244286fb519c93fe92329da96b0e6f32eejvr cy = (y2 -dy) * 3.0 34991bca4244286fb519c93fe92329da96b0e6f32eejvr bx = (x3 - x2) * 3.0 - cx 35091bca4244286fb519c93fe92329da96b0e6f32eejvr by = (y3 - y2) * 3.0 - cy 35191bca4244286fb519c93fe92329da96b0e6f32eejvr ax = x4 - dx - cx - bx 35291bca4244286fb519c93fe92329da96b0e6f32eejvr ay = y4 - dy - cy - by 35391bca4244286fb519c93fe92329da96b0e6f32eejvr return (ax, ay), (bx, by), (cx, cy), (dx, dy) 3548ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 3558ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 356befd485af59eb1d553beab340a32b02c9cb717afjvrdef calcQuadraticPoints(a, b, c): 35791bca4244286fb519c93fe92329da96b0e6f32eejvr ax, ay = a 35891bca4244286fb519c93fe92329da96b0e6f32eejvr bx, by = b 35991bca4244286fb519c93fe92329da96b0e6f32eejvr cx, cy = c 36091bca4244286fb519c93fe92329da96b0e6f32eejvr x1 = cx 36191bca4244286fb519c93fe92329da96b0e6f32eejvr y1 = cy 36291bca4244286fb519c93fe92329da96b0e6f32eejvr x2 = (bx * 0.5) + cx 36391bca4244286fb519c93fe92329da96b0e6f32eejvr y2 = (by * 0.5) + cy 36491bca4244286fb519c93fe92329da96b0e6f32eejvr x3 = ax + bx + cx 36591bca4244286fb519c93fe92329da96b0e6f32eejvr y3 = ay + by + cy 36691bca4244286fb519c93fe92329da96b0e6f32eejvr return (x1, y1), (x2, y2), (x3, y3) 367befd485af59eb1d553beab340a32b02c9cb717afjvr 368befd485af59eb1d553beab340a32b02c9cb717afjvr 369befd485af59eb1d553beab340a32b02c9cb717afjvrdef calcCubicPoints(a, b, c, d): 37091bca4244286fb519c93fe92329da96b0e6f32eejvr ax, ay = a 37191bca4244286fb519c93fe92329da96b0e6f32eejvr bx, by = b 37291bca4244286fb519c93fe92329da96b0e6f32eejvr cx, cy = c 37391bca4244286fb519c93fe92329da96b0e6f32eejvr dx, dy = d 37491bca4244286fb519c93fe92329da96b0e6f32eejvr x1 = dx 37591bca4244286fb519c93fe92329da96b0e6f32eejvr y1 = dy 37691bca4244286fb519c93fe92329da96b0e6f32eejvr x2 = (cx / 3.0) + dx 37791bca4244286fb519c93fe92329da96b0e6f32eejvr y2 = (cy / 3.0) + dy 37891bca4244286fb519c93fe92329da96b0e6f32eejvr x3 = (bx + cx) / 3.0 + x2 37991bca4244286fb519c93fe92329da96b0e6f32eejvr y3 = (by + cy) / 3.0 + y2 38091bca4244286fb519c93fe92329da96b0e6f32eejvr x4 = ax + dx + cx + bx 38191bca4244286fb519c93fe92329da96b0e6f32eejvr y4 = ay + dy + cy + by 38291bca4244286fb519c93fe92329da96b0e6f32eejvr return (x1, y1), (x2, y2), (x3, y3), (x4, y4) 383befd485af59eb1d553beab340a32b02c9cb717afjvr 384befd485af59eb1d553beab340a32b02c9cb717afjvr 385441fdd1e9f50e3897e740d55b3c5ffafda0538b6jvrdef _segmentrepr(obj): 38691bca4244286fb519c93fe92329da96b0e6f32eejvr """ 38791bca4244286fb519c93fe92329da96b0e6f32eejvr >>> _segmentrepr([1, [2, 3], [], [[2, [3, 4], [0.1, 2.2]]]]) 38891bca4244286fb519c93fe92329da96b0e6f32eejvr '(1, (2, 3), (), ((2, (3, 4), (0.1, 2.2))))' 38991bca4244286fb519c93fe92329da96b0e6f32eejvr """ 39091bca4244286fb519c93fe92329da96b0e6f32eejvr try: 39191bca4244286fb519c93fe92329da96b0e6f32eejvr it = iter(obj) 39291bca4244286fb519c93fe92329da96b0e6f32eejvr except TypeError: 39391bca4244286fb519c93fe92329da96b0e6f32eejvr return str(obj) 39491bca4244286fb519c93fe92329da96b0e6f32eejvr else: 39591bca4244286fb519c93fe92329da96b0e6f32eejvr return "(%s)" % ", ".join([_segmentrepr(x) for x in it]) 396441fdd1e9f50e3897e740d55b3c5ffafda0538b6jvr 397441fdd1e9f50e3897e740d55b3c5ffafda0538b6jvr 398441fdd1e9f50e3897e740d55b3c5ffafda0538b6jvrdef printSegments(segments): 39991bca4244286fb519c93fe92329da96b0e6f32eejvr """Helper for the doctests, displaying each segment in a list of 40091bca4244286fb519c93fe92329da96b0e6f32eejvr segments on a single line as a tuple. 40191bca4244286fb519c93fe92329da96b0e6f32eejvr """ 40291bca4244286fb519c93fe92329da96b0e6f32eejvr for segment in segments: 4033ec6a258238b6068e4eef3fe579f1f5c0a06bbbaBehdad Esfahbod print(_segmentrepr(segment)) 4048ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr 4058ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrif __name__ == "__main__": 40691bca4244286fb519c93fe92329da96b0e6f32eejvr import doctest 40791bca4244286fb519c93fe92329da96b0e6f32eejvr doctest.testmod() 408