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