```/* <![CDATA[ */
2"""
3
4from __future__ import print_function, division, absolute_import
5from fontTools.misc.py23 import *
6
7__all__ = [
9    "calcCubicBounds",
10    "splitLine",
12    "splitCubic",
14    "splitCubicAtT",
16    "solveCubic",
17]
18
19from fontTools.misc.arrayTools import calcBounds
20
21epsilon = 1e-12
22
23
24def calcQuadraticBounds(pt1, pt2, pt3):
25    """Return the bounding rectangle for a qudratic bezier segment.
26    pt1 and pt3 are the "anchor" points, pt2 is the "handle".
27
28        >>> calcQuadraticBounds((0, 0), (50, 100), (100, 0))
29        (0, 0, 100, 50.0)
30        >>> calcQuadraticBounds((0, 0), (100, 0), (100, 100))
31        (0.0, 0.0, 100, 100)
32    """
33    (ax, ay), (bx, by), (cx, cy) = calcQuadraticParameters(pt1, pt2, pt3)
34    ax2 = ax*2.0
35    ay2 = ay*2.0
36    roots = []
37    if ax2 != 0:
38        roots.append(-bx/ax2)
39    if ay2 != 0:
40        roots.append(-by/ay2)
41    points = [(ax*t*t + bx*t + cx, ay*t*t + by*t + cy) for t in roots if 0 <= t < 1] + [pt1, pt3]
42    return calcBounds(points)
43
44
45def calcCubicBounds(pt1, pt2, pt3, pt4):
46    """Return the bounding rectangle for a cubic bezier segment.
47    pt1 and pt4 are the "anchor" points, pt2 and pt3 are the "handles".
48
49        >>> calcCubicBounds((0, 0), (25, 100), (75, 100), (100, 0))
50        (0, 0, 100, 75.0)
51        >>> calcCubicBounds((0, 0), (50, 0), (100, 50), (100, 100))
52        (0.0, 0.0, 100, 100)
53        >>> print "%f %f %f %f" % calcCubicBounds((50, 0), (0, 100), (100, 100), (50, 0))
54        35.566243 0.000000 64.433757 75.000000
55    """
56    (ax, ay), (bx, by), (cx, cy), (dx, dy) = calcCubicParameters(pt1, pt2, pt3, pt4)
57    # calc first derivative
58    ax3 = ax * 3.0
59    ay3 = ay * 3.0
60    bx2 = bx * 2.0
61    by2 = by * 2.0
62    xRoots = [t for t in solveQuadratic(ax3, bx2, cx) if 0 <= t < 1]
63    yRoots = [t for t in solveQuadratic(ay3, by2, cy) if 0 <= t < 1]
64    roots = xRoots + yRoots
65
66    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]
67    return calcBounds(points)
68
69
70def splitLine(pt1, pt2, where, isHorizontal):
71    """Split the line between pt1 and pt2 at position 'where', which
72    is an x coordinate if isHorizontal is False, a y coordinate if
73    isHorizontal is True. Return a list of two line segments if the
74    line was successfully split, or a list containing the original
75    line.
76
77        >>> printSegments(splitLine((0, 0), (100, 100), 50, True))
78        ((0, 0), (50.0, 50.0))
79        ((50.0, 50.0), (100, 100))
80        >>> printSegments(splitLine((0, 0), (100, 100), 100, True))
81        ((0, 0), (100, 100))
82        >>> printSegments(splitLine((0, 0), (100, 100), 0, True))
83        ((0, 0), (0.0, 0.0))
84        ((0.0, 0.0), (100, 100))
85        >>> printSegments(splitLine((0, 0), (100, 100), 0, False))
86        ((0, 0), (0.0, 0.0))
87        ((0.0, 0.0), (100, 100))
88    """
89    pt1x, pt1y = pt1
90    pt2x, pt2y = pt2
91
92    ax = (pt2x - pt1x)
93    ay = (pt2y - pt1y)
94
95    bx = pt1x
96    by = pt1y
97
98    if ax == 0:
99        return [(pt1, pt2)]
100
101    t = (where - (bx, by)[isHorizontal]) / ax
102    if 0 <= t < 1:
103        midPt = ax * t + bx, ay * t + by
104        return [(pt1, midPt), (midPt, pt2)]
105    else:
106        return [(pt1, pt2)]
107
108
109def splitQuadratic(pt1, pt2, pt3, where, isHorizontal):
110    """Split the quadratic curve between pt1, pt2 and pt3 at position 'where',
111    which is an x coordinate if isHorizontal is False, a y coordinate if
112    isHorizontal is True. Return a list of curve segments.
113
114        >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 150, False))
115        ((0, 0), (50, 100), (100, 0))
116        >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 50, False))
117        ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0))
118        ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0))
119        >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 25, False))
120        ((0.0, 0.0), (12.5, 25.0), (25.0, 37.5))
121        ((25.0, 37.5), (62.5, 75.0), (100.0, 0.0))
122        >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 25, True))
123        ((0.0, 0.0), (7.32233047034, 14.6446609407), (14.6446609407, 25.0))
124        ((14.6446609407, 25.0), (50.0, 75.0), (85.3553390593, 25.0))
125        ((85.3553390593, 25.0), (92.6776695297, 14.6446609407), (100.0, -7.1054273576e-15))
126        >>> # XXX I'm not at all sure if the following behavior is desirable:
127        >>> printSegments(splitQuadratic((0, 0), (50, 100), (100, 0), 50, True))
128        ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0))
129        ((50.0, 50.0), (50.0, 50.0), (50.0, 50.0))
130        ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0))
131    """
132    a, b, c = calcQuadraticParameters(pt1, pt2, pt3)
133    solutions = solveQuadratic(a[isHorizontal], b[isHorizontal],
134        c[isHorizontal] - where)
135    solutions = sorted([t for t in solutions if 0 <= t < 1])
136    if not solutions:
137        return [(pt1, pt2, pt3)]
138    return _splitQuadraticAtT(a, b, c, *solutions)
139
140
141def splitCubic(pt1, pt2, pt3, pt4, where, isHorizontal):
142    """Split the cubic curve between pt1, pt2, pt3 and pt4 at position 'where',
143    which is an x coordinate if isHorizontal is False, a y coordinate if
144    isHorizontal is True. Return a list of curve segments.
145
146        >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 150, False))
147        ((0, 0), (25, 100), (75, 100), (100, 0))
148        >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 50, False))
149        ((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0))
150        ((50.0, 75.0), (68.75, 75.0), (87.5, 50.0), (100.0, 0.0))
151        >>> printSegments(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 25, True))
152        ((0.0, 0.0), (2.2937927384, 9.17517095361), (4.79804488188, 17.5085042869), (7.47413641001, 25.0))
153        ((7.47413641001, 25.0), (31.2886200204, 91.6666666667), (68.7113799796, 91.6666666667), (92.52586359, 25.0))
154        ((92.52586359, 25.0), (95.2019551181, 17.5085042869), (97.7062072616, 9.17517095361), (100.0, 1.7763568394e-15))
155    """
156    a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)
157    solutions = solveCubic(a[isHorizontal], b[isHorizontal], c[isHorizontal],
158        d[isHorizontal] - where)
159    solutions = sorted([t for t in solutions if 0 <= t < 1])
160    if not solutions:
161        return [(pt1, pt2, pt3, pt4)]
162    return _splitCubicAtT(a, b, c, d, *solutions)
163
164
165def splitQuadraticAtT(pt1, pt2, pt3, *ts):
166    """Split the quadratic curve between pt1, pt2 and pt3 at one or more
167    values of t. Return a list of curve segments.
168
169        >>> printSegments(splitQuadraticAtT((0, 0), (50, 100), (100, 0), 0.5))
170        ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0))
171        ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0))
172        >>> printSegments(splitQuadraticAtT((0, 0), (50, 100), (100, 0), 0.5, 0.75))
173        ((0.0, 0.0), (25.0, 50.0), (50.0, 50.0))
174        ((50.0, 50.0), (62.5, 50.0), (75.0, 37.5))
175        ((75.0, 37.5), (87.5, 25.0), (100.0, 0.0))
176    """
177    a, b, c = calcQuadraticParameters(pt1, pt2, pt3)
178    return _splitQuadraticAtT(a, b, c, *ts)
179
180
181def splitCubicAtT(pt1, pt2, pt3, pt4, *ts):
182    """Split the cubic curve between pt1, pt2, pt3 and pt4 at one or more
183    values of t. Return a list of curve segments.
184
185        >>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5))
186        ((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0))
187        ((50.0, 75.0), (68.75, 75.0), (87.5, 50.0), (100.0, 0.0))
188        >>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5, 0.75))
189        ((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0))
190        ((50.0, 75.0), (59.375, 75.0), (68.75, 68.75), (77.34375, 56.25))
191        ((77.34375, 56.25), (85.9375, 43.75), (93.75, 25.0), (100.0, 0.0))
192    """
193    a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)
194    return _splitCubicAtT(a, b, c, d, *ts)
195
196
197def _splitQuadraticAtT(a, b, c, *ts):
198    ts = list(ts)
199    segments = []
200    ts.insert(0, 0.0)
201    ts.append(1.0)
202    ax, ay = a
203    bx, by = b
204    cx, cy = c
205    for i in range(len(ts) - 1):
206        t1 = ts[i]
207        t2 = ts[i+1]
208        delta = (t2 - t1)
209        # calc new a, b and c
210        a1x = ax * delta**2
211        a1y = ay * delta**2
212        b1x = (2*ax*t1 + bx) * delta
213        b1y = (2*ay*t1 + by) * delta
214        c1x = ax*t1**2 + bx*t1 + cx
215        c1y = ay*t1**2 + by*t1 + cy
216
217        pt1, pt2, pt3 = calcQuadraticPoints((a1x, a1y), (b1x, b1y), (c1x, c1y))
218        segments.append((pt1, pt2, pt3))
219    return segments
220
221
222def _splitCubicAtT(a, b, c, d, *ts):
223    ts = list(ts)
224    ts.insert(0, 0.0)
225    ts.append(1.0)
226    segments = []
227    ax, ay = a
228    bx, by = b
229    cx, cy = c
230    dx, dy = d
231    for i in range(len(ts) - 1):
232        t1 = ts[i]
233        t2 = ts[i+1]
234        delta = (t2 - t1)
235        # calc new a, b, c and d
236        a1x = ax * delta**3
237        a1y = ay * delta**3
238        b1x = (3*ax*t1 + bx) * delta**2
239        b1y = (3*ay*t1 + by) * delta**2
240        c1x = (2*bx*t1 + cx + 3*ax*t1**2) * delta
241        c1y = (2*by*t1 + cy + 3*ay*t1**2) * delta
242        d1x = ax*t1**3 + bx*t1**2 + cx*t1 + dx
243        d1y = ay*t1**3 + by*t1**2 + cy*t1 + dy
244        pt1, pt2, pt3, pt4 = calcCubicPoints((a1x, a1y), (b1x, b1y), (c1x, c1y), (d1x, d1y))
245        segments.append((pt1, pt2, pt3, pt4))
246    return segments
247
248
249#
250# Equation solvers.
251#
252
253from math import sqrt, acos, cos, pi
254
255
256def solveQuadratic(a, b, c,
257        sqrt=sqrt):
258    """Solve a quadratic equation where a, b and c are real.
259        a*x*x + b*x + c = 0
260    This function returns a list of roots. Note that the returned list
261    is neither guaranteed to be sorted nor to contain unique values!
262    """
263    if abs(a) < epsilon:
264        if abs(b) < epsilon:
265            # We have a non-equation; therefore, we have no valid solution
266            roots = []
267        else:
268            # We have a linear equation with 1 root.
269            roots = [-c/b]
270    else:
271        # We have a true quadratic equation.  Apply the quadratic formula to find two roots.
272        DD = b*b - 4.0*a*c
273        if DD >= 0.0:
274            rDD = sqrt(DD)
275            roots = [(-b+rDD)/2.0/a, (-b-rDD)/2.0/a]
276        else:
277            # complex roots, ignore
278            roots = []
279    return roots
280
281
282def solveCubic(a, b, c, d):
283    """Solve a cubic equation where a, b, c and d are real.
284        a*x*x*x + b*x*x + c*x + d = 0
285    This function returns a list of roots. Note that the returned list
286    is neither guaranteed to be sorted nor to contain unique values!
287    """
288    #
289    # adapted from:
290    #   CUBIC.C - Solve a cubic polynomial
291    #   public domain by Ross Cottrell
292    # found at: http://www.strangecreations.com/library/snippets/Cubic.C
293    #
294    if abs(a) < epsilon:
295        # don't just test for zero; for very small values of 'a' solveCubic()
296        # returns unreliable results, so we fall back to quad.
297        return solveQuadratic(b, c, d)
298    a = float(a)
299    a1 = b/a
300    a2 = c/a
301    a3 = d/a
302
303    Q = (a1*a1 - 3.0*a2)/9.0
304    R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0
305    R2_Q3 = R*R - Q*Q*Q
306
307    if R2_Q3 < 0:
308        theta = acos(R/sqrt(Q*Q*Q))
309        rQ2 = -2.0*sqrt(Q)
310        x0 = rQ2*cos(theta/3.0) - a1/3.0
311        x1 = rQ2*cos((theta+2.0*pi)/3.0) - a1/3.0
312        x2 = rQ2*cos((theta+4.0*pi)/3.0) - a1/3.0
313        return [x0, x1, x2]
314    else:
315        if Q == 0 and R == 0:
316            x = 0
317        else:
318            x = pow(sqrt(R2_Q3)+abs(R), 1/3.0)
319            x = x + Q/x
320        if R >= 0.0:
321            x = -x
322        x = x - a1/3.0
323        return [x]
324
325
326#
327# Conversion routines for points to parameters and vice versa
328#
329
330def calcQuadraticParameters(pt1, pt2, pt3):
331    x2, y2 = pt2
332    x3, y3 = pt3
333    cx, cy = pt1
334    bx = (x2 - cx) * 2.0
335    by = (y2 - cy) * 2.0
336    ax = x3 - cx - bx
337    ay = y3 - cy - by
338    return (ax, ay), (bx, by), (cx, cy)
339
340
341def calcCubicParameters(pt1, pt2, pt3, pt4):
342    x2, y2 = pt2
343    x3, y3 = pt3
344    x4, y4 = pt4
345    dx, dy = pt1
346    cx = (x2 -dx) * 3.0
347    cy = (y2 -dy) * 3.0
348    bx = (x3 - x2) * 3.0 - cx
349    by = (y3 - y2) * 3.0 - cy
350    ax = x4 - dx - cx - bx
351    ay = y4 - dy - cy - by
352    return (ax, ay), (bx, by), (cx, cy), (dx, dy)
353
354
355def calcQuadraticPoints(a, b, c):
356    ax, ay = a
357    bx, by = b
358    cx, cy = c
359    x1 = cx
360    y1 = cy
361    x2 = (bx * 0.5) + cx
362    y2 = (by * 0.5) + cy
363    x3 = ax + bx + cx
364    y3 = ay + by + cy
365    return (x1, y1), (x2, y2), (x3, y3)
366
367
368def calcCubicPoints(a, b, c, d):
369    ax, ay = a
370    bx, by = b
371    cx, cy = c
372    dx, dy = d
373    x1 = dx
374    y1 = dy
375    x2 = (cx / 3.0) + dx
376    y2 = (cy / 3.0) + dy
377    x3 = (bx + cx) / 3.0 + x2
378    y3 = (by + cy) / 3.0 + y2
379    x4 = ax + dx + cx + bx
380    y4 = ay + dy + cy + by
381    return (x1, y1), (x2, y2), (x3, y3), (x4, y4)
382
383
384def _segmentrepr(obj):
385    """
386        >>> _segmentrepr([1, [2, 3], [], [[2, [3, 4], [0.1, 2.2]]]])
387        '(1, (2, 3), (), ((2, (3, 4), (0.1, 2.2))))'
388    """
389    try:
390        it = iter(obj)
391    except TypeError:
392        return str(obj)
393    else:
394        return "(%s)" % ", ".join([_segmentrepr(x) for x in it])
395
396
397def printSegments(segments):
398    """Helper for the doctests, displaying each segment in a list of
399    segments on a single line as a tuple.
400    """
401    for segment in segments:
402        print(_segmentrepr(segment))
403
404if __name__ == "__main__":
405    import doctest
406    doctest.testmod()
407```