bezierTools.py revision 32c10eecffb4923e0721c395e4b80fb732543f18
191bca4244286fb519c93fe92329da96b0e6f32eejvr"""fontTools.misc.bezierTools.py -- tools for working with bezier path segments.
291bca4244286fb519c93fe92329da96b0e6f32eejvr"""
305b4b4a27160e90307372f85dd99be69a9d972ffjvr
432c10eecffb4923e0721c395e4b80fb732543f18Behdad Esfahbodfrom __future__ import print_function, division
530e691edd056ba22fa8970280e986747817bec3dBehdad Esfahbodfrom fontTools.misc.py23 import *
605b4b4a27160e90307372f85dd99be69a9d972ffjvr
710de259aec75d3ac0c34b444b2f0423fa86a7709jvr__all__ = [
891bca4244286fb519c93fe92329da96b0e6f32eejvr    "calcQuadraticBounds",
991bca4244286fb519c93fe92329da96b0e6f32eejvr    "calcCubicBounds",
1091bca4244286fb519c93fe92329da96b0e6f32eejvr    "splitLine",
1191bca4244286fb519c93fe92329da96b0e6f32eejvr    "splitQuadratic",
1291bca4244286fb519c93fe92329da96b0e6f32eejvr    "splitCubic",
1391bca4244286fb519c93fe92329da96b0e6f32eejvr    "splitQuadraticAtT",
1491bca4244286fb519c93fe92329da96b0e6f32eejvr    "splitCubicAtT",
1591bca4244286fb519c93fe92329da96b0e6f32eejvr    "solveQuadratic",
1691bca4244286fb519c93fe92329da96b0e6f32eejvr    "solveCubic",
1710de259aec75d3ac0c34b444b2f0423fa86a7709jvr]
1805b4b4a27160e90307372f85dd99be69a9d972ffjvr
1905b4b4a27160e90307372f85dd99be69a9d972ffjvrfrom fontTools.misc.arrayTools import calcBounds
2005b4b4a27160e90307372f85dd99be69a9d972ffjvr
21c53569efef1db3270a39ed2ecea89c54d05fc1b7jvrepsilon = 1e-12
22c53569efef1db3270a39ed2ecea89c54d05fc1b7jvr
2305b4b4a27160e90307372f85dd99be69a9d972ffjvr
2405b4b4a27160e90307372f85dd99be69a9d972ffjvrdef calcQuadraticBounds(pt1, pt2, pt3):
2591bca4244286fb519c93fe92329da96b0e6f32eejvr    """Return the bounding rectangle for a qudratic bezier segment.
2691bca4244286fb519c93fe92329da96b0e6f32eejvr    pt1 and pt3 are the "anchor" points, pt2 is the "handle".
2791bca4244286fb519c93fe92329da96b0e6f32eejvr
2891bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> calcQuadraticBounds((0, 0), (50, 100), (100, 0))
2991bca4244286fb519c93fe92329da96b0e6f32eejvr        (0, 0, 100, 50.0)
3091bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> calcQuadraticBounds((0, 0), (100, 0), (100, 100))
3191bca4244286fb519c93fe92329da96b0e6f32eejvr        (0.0, 0.0, 100, 100)
3291bca4244286fb519c93fe92329da96b0e6f32eejvr    """
3391bca4244286fb519c93fe92329da96b0e6f32eejvr    (ax, ay), (bx, by), (cx, cy) = calcQuadraticParameters(pt1, pt2, pt3)
3491bca4244286fb519c93fe92329da96b0e6f32eejvr    ax2 = ax*2.0
3591bca4244286fb519c93fe92329da96b0e6f32eejvr    ay2 = ay*2.0
3691bca4244286fb519c93fe92329da96b0e6f32eejvr    roots = []
3791bca4244286fb519c93fe92329da96b0e6f32eejvr    if ax2 != 0:
3891bca4244286fb519c93fe92329da96b0e6f32eejvr        roots.append(-bx/ax2)
3991bca4244286fb519c93fe92329da96b0e6f32eejvr    if ay2 != 0:
4091bca4244286fb519c93fe92329da96b0e6f32eejvr        roots.append(-by/ay2)
4191bca4244286fb519c93fe92329da96b0e6f32eejvr    points = [(ax*t*t + bx*t + cx, ay*t*t + by*t + cy) for t in roots if 0 <= t < 1] + [pt1, pt3]
4291bca4244286fb519c93fe92329da96b0e6f32eejvr    return calcBounds(points)
4305b4b4a27160e90307372f85dd99be69a9d972ffjvr
4405b4b4a27160e90307372f85dd99be69a9d972ffjvr
4505b4b4a27160e90307372f85dd99be69a9d972ffjvrdef calcCubicBounds(pt1, pt2, pt3, pt4):
4691bca4244286fb519c93fe92329da96b0e6f32eejvr    """Return the bounding rectangle for a cubic bezier segment.
4791bca4244286fb519c93fe92329da96b0e6f32eejvr    pt1 and pt4 are the "anchor" points, pt2 and pt3 are the "handles".
4891bca4244286fb519c93fe92329da96b0e6f32eejvr
4991bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> calcCubicBounds((0, 0), (25, 100), (75, 100), (100, 0))
5091bca4244286fb519c93fe92329da96b0e6f32eejvr        (0, 0, 100, 75.0)
5191bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> calcCubicBounds((0, 0), (50, 0), (100, 50), (100, 100))
5291bca4244286fb519c93fe92329da96b0e6f32eejvr        (0.0, 0.0, 100, 100)
5391bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> print "%f %f %f %f" % calcCubicBounds((50, 0), (0, 100), (100, 100), (50, 0))
5491bca4244286fb519c93fe92329da96b0e6f32eejvr        35.566243 0.000000 64.433757 75.000000
5591bca4244286fb519c93fe92329da96b0e6f32eejvr    """
5691bca4244286fb519c93fe92329da96b0e6f32eejvr    (ax, ay), (bx, by), (cx, cy), (dx, dy) = calcCubicParameters(pt1, pt2, pt3, pt4)
5791bca4244286fb519c93fe92329da96b0e6f32eejvr    # calc first derivative
5891bca4244286fb519c93fe92329da96b0e6f32eejvr    ax3 = ax * 3.0
5991bca4244286fb519c93fe92329da96b0e6f32eejvr    ay3 = ay * 3.0
6091bca4244286fb519c93fe92329da96b0e6f32eejvr    bx2 = bx * 2.0
6191bca4244286fb519c93fe92329da96b0e6f32eejvr    by2 = by * 2.0
6291bca4244286fb519c93fe92329da96b0e6f32eejvr    xRoots = [t for t in solveQuadratic(ax3, bx2, cx) if 0 <= t < 1]
6391bca4244286fb519c93fe92329da96b0e6f32eejvr    yRoots = [t for t in solveQuadratic(ay3, by2, cy) if 0 <= t < 1]
6491bca4244286fb519c93fe92329da96b0e6f32eejvr    roots = xRoots + yRoots
6591bca4244286fb519c93fe92329da96b0e6f32eejvr
6691bca4244286fb519c93fe92329da96b0e6f32eejvr    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]
6791bca4244286fb519c93fe92329da96b0e6f32eejvr    return calcBounds(points)
6805b4b4a27160e90307372f85dd99be69a9d972ffjvr
6905b4b4a27160e90307372f85dd99be69a9d972ffjvr
7005b4b4a27160e90307372f85dd99be69a9d972ffjvrdef splitLine(pt1, pt2, where, isHorizontal):
7191bca4244286fb519c93fe92329da96b0e6f32eejvr    """Split the line between pt1 and pt2 at position 'where', which
7291bca4244286fb519c93fe92329da96b0e6f32eejvr    is an x coordinate if isHorizontal is False, a y coordinate if
7391bca4244286fb519c93fe92329da96b0e6f32eejvr    isHorizontal is True. Return a list of two line segments if the
7491bca4244286fb519c93fe92329da96b0e6f32eejvr    line was successfully split, or a list containing the original
7591bca4244286fb519c93fe92329da96b0e6f32eejvr    line.
7691bca4244286fb519c93fe92329da96b0e6f32eejvr
7791bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitLine((0, 0), (100, 100), 50, True))
7891bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0, 0), (50.0, 50.0))
7991bca4244286fb519c93fe92329da96b0e6f32eejvr        ((50.0, 50.0), (100, 100))
8091bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitLine((0, 0), (100, 100), 100, True))
8191bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0, 0), (100, 100))
8291bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitLine((0, 0), (100, 100), 0, True))
8391bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0, 0), (0.0, 0.0))
8491bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (100, 100))
8591bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitLine((0, 0), (100, 100), 0, False))
8691bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0, 0), (0.0, 0.0))
8791bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (100, 100))
8891bca4244286fb519c93fe92329da96b0e6f32eejvr    """
8991bca4244286fb519c93fe92329da96b0e6f32eejvr    pt1x, pt1y = pt1
9091bca4244286fb519c93fe92329da96b0e6f32eejvr    pt2x, pt2y = pt2
9191bca4244286fb519c93fe92329da96b0e6f32eejvr
9291bca4244286fb519c93fe92329da96b0e6f32eejvr    ax = (pt2x - pt1x)
9391bca4244286fb519c93fe92329da96b0e6f32eejvr    ay = (pt2y - pt1y)
9491bca4244286fb519c93fe92329da96b0e6f32eejvr
9591bca4244286fb519c93fe92329da96b0e6f32eejvr    bx = pt1x
9691bca4244286fb519c93fe92329da96b0e6f32eejvr    by = pt1y
9791bca4244286fb519c93fe92329da96b0e6f32eejvr
9891bca4244286fb519c93fe92329da96b0e6f32eejvr    ax1 = (ax, ay)[isHorizontal]
9991bca4244286fb519c93fe92329da96b0e6f32eejvr
10091bca4244286fb519c93fe92329da96b0e6f32eejvr    if ax == 0:
10191bca4244286fb519c93fe92329da96b0e6f32eejvr        return [(pt1, pt2)]
10291bca4244286fb519c93fe92329da96b0e6f32eejvr
10332c10eecffb4923e0721c395e4b80fb732543f18Behdad Esfahbod    t = (where - (bx, by)[isHorizontal]) / ax
10491bca4244286fb519c93fe92329da96b0e6f32eejvr    if 0 <= t < 1:
10591bca4244286fb519c93fe92329da96b0e6f32eejvr        midPt = ax * t + bx, ay * t + by
10691bca4244286fb519c93fe92329da96b0e6f32eejvr        return [(pt1, midPt), (midPt, pt2)]
10791bca4244286fb519c93fe92329da96b0e6f32eejvr    else:
10891bca4244286fb519c93fe92329da96b0e6f32eejvr        return [(pt1, pt2)]
10905b4b4a27160e90307372f85dd99be69a9d972ffjvr
11005b4b4a27160e90307372f85dd99be69a9d972ffjvr
11105b4b4a27160e90307372f85dd99be69a9d972ffjvrdef splitQuadratic(pt1, pt2, pt3, where, isHorizontal):
11291bca4244286fb519c93fe92329da96b0e6f32eejvr    """Split the quadratic curve between pt1, pt2 and pt3 at position 'where',
11391bca4244286fb519c93fe92329da96b0e6f32eejvr    which is an x coordinate if isHorizontal is False, a y coordinate if
11491bca4244286fb519c93fe92329da96b0e6f32eejvr    isHorizontal is True. Return a list of curve segments.
11591bca4244286fb519c93fe92329da96b0e6f32eejvr
11691bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 150, False))
11791bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0, 0), (50, 100), (100, 0))
11891bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 50, False))
11991bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0))
12091bca4244286fb519c93fe92329da96b0e6f32eejvr        ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0))
12191bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 25, False))
12291bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (12.5, 25.0), (25.0, 37.5))
12391bca4244286fb519c93fe92329da96b0e6f32eejvr        ((25.0, 37.5), (62.5, 75.0), (100.0, 0.0))
12491bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 25, True))
12591bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (7.32233047034, 14.6446609407), (14.6446609407, 25.0))
12691bca4244286fb519c93fe92329da96b0e6f32eejvr        ((14.6446609407, 25.0), (50.0, 75.0), (85.3553390593, 25.0))
12791bca4244286fb519c93fe92329da96b0e6f32eejvr        ((85.3553390593, 25.0), (92.6776695297, 14.6446609407), (100.0, -7.1054273576e-15))
12891bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> # XXX I'm not at all sure if the following behavior is desirable:
12991bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 50, True))
13091bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0))
13191bca4244286fb519c93fe92329da96b0e6f32eejvr        ((50.0, 50.0), (50.0, 50.0), (50.0, 50.0))
13291bca4244286fb519c93fe92329da96b0e6f32eejvr        ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0))
13391bca4244286fb519c93fe92329da96b0e6f32eejvr    """
13491bca4244286fb519c93fe92329da96b0e6f32eejvr    a, b, c = calcQuadraticParameters(pt1, pt2, pt3)
13591bca4244286fb519c93fe92329da96b0e6f32eejvr    solutions = solveQuadratic(a[isHorizontal], b[isHorizontal],
13691bca4244286fb519c93fe92329da96b0e6f32eejvr        c[isHorizontal] - where)
137ac1b4359467ca3deab03186a15eae1d55eb35567Behdad Esfahbod    solutions = sorted([t for t in solutions if 0 <= t < 1])
13891bca4244286fb519c93fe92329da96b0e6f32eejvr    if not solutions:
13991bca4244286fb519c93fe92329da96b0e6f32eejvr        return [(pt1, pt2, pt3)]
14091bca4244286fb519c93fe92329da96b0e6f32eejvr    return _splitQuadraticAtT(a, b, c, *solutions)
14105b4b4a27160e90307372f85dd99be69a9d972ffjvr
14205b4b4a27160e90307372f85dd99be69a9d972ffjvr
14305b4b4a27160e90307372f85dd99be69a9d972ffjvrdef splitCubic(pt1, pt2, pt3, pt4, where, isHorizontal):
14491bca4244286fb519c93fe92329da96b0e6f32eejvr    """Split the cubic curve between pt1, pt2, pt3 and pt4 at position 'where',
14591bca4244286fb519c93fe92329da96b0e6f32eejvr    which is an x coordinate if isHorizontal is False, a y coordinate if
14691bca4244286fb519c93fe92329da96b0e6f32eejvr    isHorizontal is True. Return a list of curve segments.
14791bca4244286fb519c93fe92329da96b0e6f32eejvr
14891bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 150, False))
14991bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0, 0), (25, 100), (75, 100), (100, 0))
15091bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 50, False))
15191bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0))
15291bca4244286fb519c93fe92329da96b0e6f32eejvr        ((50.0, 75.0), (68.75, 75.0), (87.5, 50.0), (100.0, 0.0))
15391bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 25, True))
15491bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (2.2937927384, 9.17517095361), (4.79804488188, 17.5085042869), (7.47413641001, 25.0))
15591bca4244286fb519c93fe92329da96b0e6f32eejvr        ((7.47413641001, 25.0), (31.2886200204, 91.6666666667), (68.7113799796, 91.6666666667), (92.52586359, 25.0))
15691bca4244286fb519c93fe92329da96b0e6f32eejvr        ((92.52586359, 25.0), (95.2019551181, 17.5085042869), (97.7062072616, 9.17517095361), (100.0, 1.7763568394e-15))
15791bca4244286fb519c93fe92329da96b0e6f32eejvr    """
15891bca4244286fb519c93fe92329da96b0e6f32eejvr    a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)
15991bca4244286fb519c93fe92329da96b0e6f32eejvr    solutions = solveCubic(a[isHorizontal], b[isHorizontal], c[isHorizontal],
16091bca4244286fb519c93fe92329da96b0e6f32eejvr        d[isHorizontal] - where)
161ac1b4359467ca3deab03186a15eae1d55eb35567Behdad Esfahbod    solutions = sorted([t for t in solutions if 0 <= t < 1])
16291bca4244286fb519c93fe92329da96b0e6f32eejvr    if not solutions:
16391bca4244286fb519c93fe92329da96b0e6f32eejvr        return [(pt1, pt2, pt3, pt4)]
16491bca4244286fb519c93fe92329da96b0e6f32eejvr    return _splitCubicAtT(a, b, c, d, *solutions)
16510de259aec75d3ac0c34b444b2f0423fa86a7709jvr
16610de259aec75d3ac0c34b444b2f0423fa86a7709jvr
16710de259aec75d3ac0c34b444b2f0423fa86a7709jvrdef splitQuadraticAtT(pt1, pt2, pt3, *ts):
16891bca4244286fb519c93fe92329da96b0e6f32eejvr    """Split the quadratic curve between pt1, pt2 and pt3 at one or more
16991bca4244286fb519c93fe92329da96b0e6f32eejvr    values of t. Return a list of curve segments.
17086c07d2b9ad2570823119ce453ad07275a09d94cjvr
17191bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitQuadraticAtT((0, 0), (50, 100), (100, 0), 0.5))
17291bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0))
17391bca4244286fb519c93fe92329da96b0e6f32eejvr        ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0))
17491bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitQuadraticAtT((0, 0), (50, 100), (100, 0), 0.5, 0.75))
17591bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0))
17691bca4244286fb519c93fe92329da96b0e6f32eejvr        ((50.0, 50.0), (62.5, 50.0), (75.0, 37.5))
17791bca4244286fb519c93fe92329da96b0e6f32eejvr        ((75.0, 37.5), (87.5, 25.0), (100.0, 0.0))
17891bca4244286fb519c93fe92329da96b0e6f32eejvr    """
17991bca4244286fb519c93fe92329da96b0e6f32eejvr    a, b, c = calcQuadraticParameters(pt1, pt2, pt3)
18091bca4244286fb519c93fe92329da96b0e6f32eejvr    return _splitQuadraticAtT(a, b, c, *ts)
18110de259aec75d3ac0c34b444b2f0423fa86a7709jvr
18210de259aec75d3ac0c34b444b2f0423fa86a7709jvr
18310de259aec75d3ac0c34b444b2f0423fa86a7709jvrdef splitCubicAtT(pt1, pt2, pt3, pt4, *ts):
18491bca4244286fb519c93fe92329da96b0e6f32eejvr    """Split the cubic curve between pt1, pt2, pt3 and pt4 at one or more
18591bca4244286fb519c93fe92329da96b0e6f32eejvr    values of t. Return a list of curve segments.
18686c07d2b9ad2570823119ce453ad07275a09d94cjvr
18791bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5))
18891bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0))
18991bca4244286fb519c93fe92329da96b0e6f32eejvr        ((50.0, 75.0), (68.75, 75.0), (87.5, 50.0), (100.0, 0.0))
19091bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5, 0.75))
19191bca4244286fb519c93fe92329da96b0e6f32eejvr        ((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0))
19291bca4244286fb519c93fe92329da96b0e6f32eejvr        ((50.0, 75.0), (59.375, 75.0), (68.75, 68.75), (77.34375, 56.25))
19391bca4244286fb519c93fe92329da96b0e6f32eejvr        ((77.34375, 56.25), (85.9375, 43.75), (93.75, 25.0), (100.0, 0.0))
19491bca4244286fb519c93fe92329da96b0e6f32eejvr    """
19591bca4244286fb519c93fe92329da96b0e6f32eejvr    a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)
19691bca4244286fb519c93fe92329da96b0e6f32eejvr    return _splitCubicAtT(a, b, c, d, *ts)
19710de259aec75d3ac0c34b444b2f0423fa86a7709jvr
19810de259aec75d3ac0c34b444b2f0423fa86a7709jvr
19910de259aec75d3ac0c34b444b2f0423fa86a7709jvrdef _splitQuadraticAtT(a, b, c, *ts):
20091bca4244286fb519c93fe92329da96b0e6f32eejvr    ts = list(ts)
20191bca4244286fb519c93fe92329da96b0e6f32eejvr    segments = []
20291bca4244286fb519c93fe92329da96b0e6f32eejvr    ts.insert(0, 0.0)
20391bca4244286fb519c93fe92329da96b0e6f32eejvr    ts.append(1.0)
20491bca4244286fb519c93fe92329da96b0e6f32eejvr    ax, ay = a
20591bca4244286fb519c93fe92329da96b0e6f32eejvr    bx, by = b
20691bca4244286fb519c93fe92329da96b0e6f32eejvr    cx, cy = c
20791bca4244286fb519c93fe92329da96b0e6f32eejvr    for i in range(len(ts) - 1):
20891bca4244286fb519c93fe92329da96b0e6f32eejvr        t1 = ts[i]
20991bca4244286fb519c93fe92329da96b0e6f32eejvr        t2 = ts[i+1]
21091bca4244286fb519c93fe92329da96b0e6f32eejvr        delta = (t2 - t1)
21191bca4244286fb519c93fe92329da96b0e6f32eejvr        # calc new a, b and c
21291bca4244286fb519c93fe92329da96b0e6f32eejvr        a1x = ax * delta**2
21391bca4244286fb519c93fe92329da96b0e6f32eejvr        a1y = ay * delta**2
21491bca4244286fb519c93fe92329da96b0e6f32eejvr        b1x = (2*ax*t1 + bx) * delta
21591bca4244286fb519c93fe92329da96b0e6f32eejvr        b1y = (2*ay*t1 + by) * delta
21691bca4244286fb519c93fe92329da96b0e6f32eejvr        c1x = ax*t1**2 + bx*t1 + cx
21791bca4244286fb519c93fe92329da96b0e6f32eejvr        c1y = ay*t1**2 + by*t1 + cy
21891bca4244286fb519c93fe92329da96b0e6f32eejvr
21991bca4244286fb519c93fe92329da96b0e6f32eejvr        pt1, pt2, pt3 = calcQuadraticPoints((a1x, a1y), (b1x, b1y), (c1x, c1y))
22091bca4244286fb519c93fe92329da96b0e6f32eejvr        segments.append((pt1, pt2, pt3))
22191bca4244286fb519c93fe92329da96b0e6f32eejvr    return segments
22210de259aec75d3ac0c34b444b2f0423fa86a7709jvr
22310de259aec75d3ac0c34b444b2f0423fa86a7709jvr
22410de259aec75d3ac0c34b444b2f0423fa86a7709jvrdef _splitCubicAtT(a, b, c, d, *ts):
22591bca4244286fb519c93fe92329da96b0e6f32eejvr    ts = list(ts)
22691bca4244286fb519c93fe92329da96b0e6f32eejvr    ts.insert(0, 0.0)
22791bca4244286fb519c93fe92329da96b0e6f32eejvr    ts.append(1.0)
22891bca4244286fb519c93fe92329da96b0e6f32eejvr    segments = []
22991bca4244286fb519c93fe92329da96b0e6f32eejvr    ax, ay = a
23091bca4244286fb519c93fe92329da96b0e6f32eejvr    bx, by = b
23191bca4244286fb519c93fe92329da96b0e6f32eejvr    cx, cy = c
23291bca4244286fb519c93fe92329da96b0e6f32eejvr    dx, dy = d
23391bca4244286fb519c93fe92329da96b0e6f32eejvr    for i in range(len(ts) - 1):
23491bca4244286fb519c93fe92329da96b0e6f32eejvr        t1 = ts[i]
23591bca4244286fb519c93fe92329da96b0e6f32eejvr        t2 = ts[i+1]
23691bca4244286fb519c93fe92329da96b0e6f32eejvr        delta = (t2 - t1)
23791bca4244286fb519c93fe92329da96b0e6f32eejvr        # calc new a, b, c and d
23891bca4244286fb519c93fe92329da96b0e6f32eejvr        a1x = ax * delta**3
23991bca4244286fb519c93fe92329da96b0e6f32eejvr        a1y = ay * delta**3
24091bca4244286fb519c93fe92329da96b0e6f32eejvr        b1x = (3*ax*t1 + bx) * delta**2
24191bca4244286fb519c93fe92329da96b0e6f32eejvr        b1y = (3*ay*t1 + by) * delta**2
24291bca4244286fb519c93fe92329da96b0e6f32eejvr        c1x = (2*bx*t1 + cx + 3*ax*t1**2) * delta
24391bca4244286fb519c93fe92329da96b0e6f32eejvr        c1y = (2*by*t1 + cy + 3*ay*t1**2) * delta
24491bca4244286fb519c93fe92329da96b0e6f32eejvr        d1x = ax*t1**3 + bx*t1**2 + cx*t1 + dx
24591bca4244286fb519c93fe92329da96b0e6f32eejvr        d1y = ay*t1**3 + by*t1**2 + cy*t1 + dy
24691bca4244286fb519c93fe92329da96b0e6f32eejvr        pt1, pt2, pt3, pt4 = calcCubicPoints((a1x, a1y), (b1x, b1y), (c1x, c1y), (d1x, d1y))
24791bca4244286fb519c93fe92329da96b0e6f32eejvr        segments.append((pt1, pt2, pt3, pt4))
24891bca4244286fb519c93fe92329da96b0e6f32eejvr    return segments
24905b4b4a27160e90307372f85dd99be69a9d972ffjvr
25005b4b4a27160e90307372f85dd99be69a9d972ffjvr
25105b4b4a27160e90307372f85dd99be69a9d972ffjvr#
25205b4b4a27160e90307372f85dd99be69a9d972ffjvr# Equation solvers.
25305b4b4a27160e90307372f85dd99be69a9d972ffjvr#
25405b4b4a27160e90307372f85dd99be69a9d972ffjvr
25505b4b4a27160e90307372f85dd99be69a9d972ffjvrfrom math import sqrt, acos, cos, pi
25605b4b4a27160e90307372f85dd99be69a9d972ffjvr
25705b4b4a27160e90307372f85dd99be69a9d972ffjvr
25805b4b4a27160e90307372f85dd99be69a9d972ffjvrdef solveQuadratic(a, b, c,
25991bca4244286fb519c93fe92329da96b0e6f32eejvr        sqrt=sqrt):
26091bca4244286fb519c93fe92329da96b0e6f32eejvr    """Solve a quadratic equation where a, b and c are real.
26191bca4244286fb519c93fe92329da96b0e6f32eejvr        a*x*x + b*x + c = 0
26291bca4244286fb519c93fe92329da96b0e6f32eejvr    This function returns a list of roots. Note that the returned list
26391bca4244286fb519c93fe92329da96b0e6f32eejvr    is neither guaranteed to be sorted nor to contain unique values!
26491bca4244286fb519c93fe92329da96b0e6f32eejvr    """
26591bca4244286fb519c93fe92329da96b0e6f32eejvr    if abs(a) < epsilon:
26691bca4244286fb519c93fe92329da96b0e6f32eejvr        if abs(b) < epsilon:
26791bca4244286fb519c93fe92329da96b0e6f32eejvr            # We have a non-equation; therefore, we have no valid solution
26891bca4244286fb519c93fe92329da96b0e6f32eejvr            roots = []
26991bca4244286fb519c93fe92329da96b0e6f32eejvr        else:
27091bca4244286fb519c93fe92329da96b0e6f32eejvr            # We have a linear equation with 1 root.
27191bca4244286fb519c93fe92329da96b0e6f32eejvr            roots = [-c/b]
27291bca4244286fb519c93fe92329da96b0e6f32eejvr    else:
27391bca4244286fb519c93fe92329da96b0e6f32eejvr        # We have a true quadratic equation.  Apply the quadratic formula to find two roots.
27491bca4244286fb519c93fe92329da96b0e6f32eejvr        DD = b*b - 4.0*a*c
27591bca4244286fb519c93fe92329da96b0e6f32eejvr        if DD >= 0.0:
27691bca4244286fb519c93fe92329da96b0e6f32eejvr            rDD = sqrt(DD)
27791bca4244286fb519c93fe92329da96b0e6f32eejvr            roots = [(-b+rDD)/2.0/a, (-b-rDD)/2.0/a]
27891bca4244286fb519c93fe92329da96b0e6f32eejvr        else:
27991bca4244286fb519c93fe92329da96b0e6f32eejvr            # complex roots, ignore
28091bca4244286fb519c93fe92329da96b0e6f32eejvr            roots = []
28191bca4244286fb519c93fe92329da96b0e6f32eejvr    return roots
28205b4b4a27160e90307372f85dd99be69a9d972ffjvr
28305b4b4a27160e90307372f85dd99be69a9d972ffjvr
28405b4b4a27160e90307372f85dd99be69a9d972ffjvrdef solveCubic(a, b, c, d,
28591bca4244286fb519c93fe92329da96b0e6f32eejvr        abs=abs, pow=pow, sqrt=sqrt, cos=cos, acos=acos, pi=pi):
28691bca4244286fb519c93fe92329da96b0e6f32eejvr    """Solve a cubic equation where a, b, c and d are real.
28791bca4244286fb519c93fe92329da96b0e6f32eejvr        a*x*x*x + b*x*x + c*x + d = 0
28891bca4244286fb519c93fe92329da96b0e6f32eejvr    This function returns a list of roots. Note that the returned list
28991bca4244286fb519c93fe92329da96b0e6f32eejvr    is neither guaranteed to be sorted nor to contain unique values!
29091bca4244286fb519c93fe92329da96b0e6f32eejvr    """
29191bca4244286fb519c93fe92329da96b0e6f32eejvr    #
29291bca4244286fb519c93fe92329da96b0e6f32eejvr    # adapted from:
29391bca4244286fb519c93fe92329da96b0e6f32eejvr    #   CUBIC.C - Solve a cubic polynomial
29491bca4244286fb519c93fe92329da96b0e6f32eejvr    #   public domain by Ross Cottrell
29591bca4244286fb519c93fe92329da96b0e6f32eejvr    # found at: http://www.strangecreations.com/library/snippets/Cubic.C
29691bca4244286fb519c93fe92329da96b0e6f32eejvr    #
29791bca4244286fb519c93fe92329da96b0e6f32eejvr    if abs(a) < epsilon:
29891bca4244286fb519c93fe92329da96b0e6f32eejvr        # don't just test for zero; for very small values of 'a' solveCubic()
29991bca4244286fb519c93fe92329da96b0e6f32eejvr        # returns unreliable results, so we fall back to quad.
30091bca4244286fb519c93fe92329da96b0e6f32eejvr        return solveQuadratic(b, c, d)
30191bca4244286fb519c93fe92329da96b0e6f32eejvr    a = float(a)
30291bca4244286fb519c93fe92329da96b0e6f32eejvr    a1 = b/a
30391bca4244286fb519c93fe92329da96b0e6f32eejvr    a2 = c/a
30491bca4244286fb519c93fe92329da96b0e6f32eejvr    a3 = d/a
30591bca4244286fb519c93fe92329da96b0e6f32eejvr
30691bca4244286fb519c93fe92329da96b0e6f32eejvr    Q = (a1*a1 - 3.0*a2)/9.0
30791bca4244286fb519c93fe92329da96b0e6f32eejvr    R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0
30891bca4244286fb519c93fe92329da96b0e6f32eejvr    R2_Q3 = R*R - Q*Q*Q
30991bca4244286fb519c93fe92329da96b0e6f32eejvr
31091bca4244286fb519c93fe92329da96b0e6f32eejvr    if R2_Q3 < 0:
31191bca4244286fb519c93fe92329da96b0e6f32eejvr        theta = acos(R/sqrt(Q*Q*Q))
31291bca4244286fb519c93fe92329da96b0e6f32eejvr        rQ2 = -2.0*sqrt(Q)
31391bca4244286fb519c93fe92329da96b0e6f32eejvr        x0 = rQ2*cos(theta/3.0) - a1/3.0
31491bca4244286fb519c93fe92329da96b0e6f32eejvr        x1 = rQ2*cos((theta+2.0*pi)/3.0) - a1/3.0
31591bca4244286fb519c93fe92329da96b0e6f32eejvr        x2 = rQ2*cos((theta+4.0*pi)/3.0) - a1/3.0
31691bca4244286fb519c93fe92329da96b0e6f32eejvr        return [x0, x1, x2]
31791bca4244286fb519c93fe92329da96b0e6f32eejvr    else:
31891bca4244286fb519c93fe92329da96b0e6f32eejvr        if Q == 0 and R == 0:
31991bca4244286fb519c93fe92329da96b0e6f32eejvr            x = 0
32091bca4244286fb519c93fe92329da96b0e6f32eejvr        else:
32191bca4244286fb519c93fe92329da96b0e6f32eejvr            x = pow(sqrt(R2_Q3)+abs(R), 1/3.0)
32291bca4244286fb519c93fe92329da96b0e6f32eejvr            x = x + Q/x
32391bca4244286fb519c93fe92329da96b0e6f32eejvr        if R >= 0.0:
32491bca4244286fb519c93fe92329da96b0e6f32eejvr            x = -x
32591bca4244286fb519c93fe92329da96b0e6f32eejvr        x = x - a1/3.0
32691bca4244286fb519c93fe92329da96b0e6f32eejvr        return [x]
3278ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr
3288ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr
32982fdf153acf0c82926ad5f49b473b9dba71ba3d3jvr#
33082fdf153acf0c82926ad5f49b473b9dba71ba3d3jvr# Conversion routines for points to parameters and vice versa
33182fdf153acf0c82926ad5f49b473b9dba71ba3d3jvr#
33282fdf153acf0c82926ad5f49b473b9dba71ba3d3jvr
3338ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrdef calcQuadraticParameters(pt1, pt2, pt3):
33491bca4244286fb519c93fe92329da96b0e6f32eejvr    x2, y2 = pt2
33591bca4244286fb519c93fe92329da96b0e6f32eejvr    x3, y3 = pt3
33691bca4244286fb519c93fe92329da96b0e6f32eejvr    cx, cy = pt1
33791bca4244286fb519c93fe92329da96b0e6f32eejvr    bx = (x2 - cx) * 2.0
33891bca4244286fb519c93fe92329da96b0e6f32eejvr    by = (y2 - cy) * 2.0
33991bca4244286fb519c93fe92329da96b0e6f32eejvr    ax = x3 - cx - bx
34091bca4244286fb519c93fe92329da96b0e6f32eejvr    ay = y3 - cy - by
34191bca4244286fb519c93fe92329da96b0e6f32eejvr    return (ax, ay), (bx, by), (cx, cy)
3428ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr
3438ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr
3448ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrdef calcCubicParameters(pt1, pt2, pt3, pt4):
34591bca4244286fb519c93fe92329da96b0e6f32eejvr    x2, y2 = pt2
34691bca4244286fb519c93fe92329da96b0e6f32eejvr    x3, y3 = pt3
34791bca4244286fb519c93fe92329da96b0e6f32eejvr    x4, y4 = pt4
34891bca4244286fb519c93fe92329da96b0e6f32eejvr    dx, dy = pt1
34991bca4244286fb519c93fe92329da96b0e6f32eejvr    cx = (x2 -dx) * 3.0
35091bca4244286fb519c93fe92329da96b0e6f32eejvr    cy = (y2 -dy) * 3.0
35191bca4244286fb519c93fe92329da96b0e6f32eejvr    bx = (x3 - x2) * 3.0 - cx
35291bca4244286fb519c93fe92329da96b0e6f32eejvr    by = (y3 - y2) * 3.0 - cy
35391bca4244286fb519c93fe92329da96b0e6f32eejvr    ax = x4 - dx - cx - bx
35491bca4244286fb519c93fe92329da96b0e6f32eejvr    ay = y4 - dy - cy - by
35591bca4244286fb519c93fe92329da96b0e6f32eejvr    return (ax, ay), (bx, by), (cx, cy), (dx, dy)
3568ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr
3578ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr
358befd485af59eb1d553beab340a32b02c9cb717afjvrdef calcQuadraticPoints(a, b, c):
35991bca4244286fb519c93fe92329da96b0e6f32eejvr    ax, ay = a
36091bca4244286fb519c93fe92329da96b0e6f32eejvr    bx, by = b
36191bca4244286fb519c93fe92329da96b0e6f32eejvr    cx, cy = c
36291bca4244286fb519c93fe92329da96b0e6f32eejvr    x1 = cx
36391bca4244286fb519c93fe92329da96b0e6f32eejvr    y1 = cy
36491bca4244286fb519c93fe92329da96b0e6f32eejvr    x2 = (bx * 0.5) + cx
36591bca4244286fb519c93fe92329da96b0e6f32eejvr    y2 = (by * 0.5) + cy
36691bca4244286fb519c93fe92329da96b0e6f32eejvr    x3 = ax + bx + cx
36791bca4244286fb519c93fe92329da96b0e6f32eejvr    y3 = ay + by + cy
36891bca4244286fb519c93fe92329da96b0e6f32eejvr    return (x1, y1), (x2, y2), (x3, y3)
369befd485af59eb1d553beab340a32b02c9cb717afjvr
370befd485af59eb1d553beab340a32b02c9cb717afjvr
371befd485af59eb1d553beab340a32b02c9cb717afjvrdef calcCubicPoints(a, b, c, d):
37291bca4244286fb519c93fe92329da96b0e6f32eejvr    ax, ay = a
37391bca4244286fb519c93fe92329da96b0e6f32eejvr    bx, by = b
37491bca4244286fb519c93fe92329da96b0e6f32eejvr    cx, cy = c
37591bca4244286fb519c93fe92329da96b0e6f32eejvr    dx, dy = d
37691bca4244286fb519c93fe92329da96b0e6f32eejvr    x1 = dx
37791bca4244286fb519c93fe92329da96b0e6f32eejvr    y1 = dy
37891bca4244286fb519c93fe92329da96b0e6f32eejvr    x2 = (cx / 3.0) + dx
37991bca4244286fb519c93fe92329da96b0e6f32eejvr    y2 = (cy / 3.0) + dy
38091bca4244286fb519c93fe92329da96b0e6f32eejvr    x3 = (bx + cx) / 3.0 + x2
38191bca4244286fb519c93fe92329da96b0e6f32eejvr    y3 = (by + cy) / 3.0 + y2
38291bca4244286fb519c93fe92329da96b0e6f32eejvr    x4 = ax + dx + cx + bx
38391bca4244286fb519c93fe92329da96b0e6f32eejvr    y4 = ay + dy + cy + by
38491bca4244286fb519c93fe92329da96b0e6f32eejvr    return (x1, y1), (x2, y2), (x3, y3), (x4, y4)
385befd485af59eb1d553beab340a32b02c9cb717afjvr
386befd485af59eb1d553beab340a32b02c9cb717afjvr
387441fdd1e9f50e3897e740d55b3c5ffafda0538b6jvrdef _segmentrepr(obj):
38891bca4244286fb519c93fe92329da96b0e6f32eejvr    """
38991bca4244286fb519c93fe92329da96b0e6f32eejvr        >>> _segmentrepr([1, [2, 3], [], [[2, [3, 4], [0.1, 2.2]]]])
39091bca4244286fb519c93fe92329da96b0e6f32eejvr        '(1, (2, 3), (), ((2, (3, 4), (0.1, 2.2))))'
39191bca4244286fb519c93fe92329da96b0e6f32eejvr    """
39291bca4244286fb519c93fe92329da96b0e6f32eejvr    try:
39391bca4244286fb519c93fe92329da96b0e6f32eejvr        it = iter(obj)
39491bca4244286fb519c93fe92329da96b0e6f32eejvr    except TypeError:
39591bca4244286fb519c93fe92329da96b0e6f32eejvr        return str(obj)
39691bca4244286fb519c93fe92329da96b0e6f32eejvr    else:
39791bca4244286fb519c93fe92329da96b0e6f32eejvr        return "(%s)" % ", ".join([_segmentrepr(x) for x in it])
398441fdd1e9f50e3897e740d55b3c5ffafda0538b6jvr
399441fdd1e9f50e3897e740d55b3c5ffafda0538b6jvr
400441fdd1e9f50e3897e740d55b3c5ffafda0538b6jvrdef printSegments(segments):
40191bca4244286fb519c93fe92329da96b0e6f32eejvr    """Helper for the doctests, displaying each segment in a list of
40291bca4244286fb519c93fe92329da96b0e6f32eejvr    segments on a single line as a tuple.
40391bca4244286fb519c93fe92329da96b0e6f32eejvr    """
40491bca4244286fb519c93fe92329da96b0e6f32eejvr    for segment in segments:
4053ec6a258238b6068e4eef3fe579f1f5c0a06bbbaBehdad Esfahbod        print(_segmentrepr(segment))
4068ee2561bc1b45e1e0ed328c392c31137878dc0d8jvr
4078ee2561bc1b45e1e0ed328c392c31137878dc0d8jvrif __name__ == "__main__":
40891bca4244286fb519c93fe92329da96b0e6f32eejvr    import doctest
40991bca4244286fb519c93fe92329da96b0e6f32eejvr    doctest.testmod()
410