bezierTools.py revision befd485af59eb1d553beab340a32b02c9cb717af 
1"""fontTools.misc.bezierTools.py -- tools for working with bezier path segments."""
2
3
4__all__ = ["calcQuadraticBounds", "calcCubicBounds", "splitLine", "splitQuadratic",
5	"splitCubic", "solveQuadratic", "solveCubic"]
6
7
8from fontTools.misc.arrayTools import calcBounds
9import Numeric
10
11
12def calcQuadraticBounds(pt1, pt2, pt3):
13	"""Return the bounding rectangle for a qudratic bezier segment.
14	pt1 and pt3 are the "anchor" points, pt2 is the "handle".
15
16		>>> calcQuadraticBounds((0, 0), (50, 100), (100, 0))
17		(0.0, 0.0, 100.0, 50.0)
18		>>> calcQuadraticBounds((0, 0), (100, 0), (100, 100))
19		(0.0, 0.0, 100.0, 100.0)
20	"""
21	a, b, c = calcQuadraticParameters(pt1, pt2, pt3)
22	# calc first derivative
23	ax, ay = a * 2
24	bx, by = b
25	roots = []
26	if ax != 0:
27		roots.append(-bx/ax)
28	if ay != 0:
29		roots.append(-by/ay)
30	points = [a*t*t + b*t + c for t in roots if 0 <= t < 1] + [pt1, pt3]
31	return calcBounds(points)
32
33
34def calcCubicBounds(pt1, pt2, pt3, pt4):
35	"""Return the bounding rectangle for a cubic bezier segment.
36	pt1 and pt4 are the "anchor" points, pt2 and pt3 are the "handles".
37
38		>>> calcCubicBounds((0, 0), (25, 100), (75, 100), (100, 0))
39		(0.0, 0.0, 100.0, 75.0)
40		>>> calcCubicBounds((0, 0), (50, 0), (100, 50), (100, 100))
41		(0.0, 0.0, 100.0, 100.0)
42		>>> calcCubicBounds((50, 0), (0, 100), (100, 100), (50, 0))
43		(35.566243270259356, 0.0, 64.433756729740679, 75.0)
44	"""
45	a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)
46	# calc first derivative
47	ax, ay = a * 3.0
48	bx, by = b * 2.0
49	cx, cy = c
50	xRoots = [t for t in solveQuadratic(ax, bx, cx) if 0 <= t < 1]
51	yRoots = [t for t in solveQuadratic(ay, by, cy) if 0 <= t < 1]
52	roots = xRoots + yRoots
53
54	points = [(a*t*t*t + b*t*t + c * t + d) for t in roots] + [pt1, pt4]
55	return calcBounds(points)
56
57
58def splitLine(pt1, pt2, where, isHorizontal):
59	"""Split the line between pt1 and pt2 at position 'where', which
60	is an x coordinate if isHorizontal is False, a y coordinate if
61	isHorizontal is True. Return a list of two line segments if the
62	line was successfully split, or a list containing the original
63	line.
64
65		>>> _testrepr(splitLine((0, 0), (100, 100), 50, True))
66		'(((0, 0), (50.0, 50.0)), ((50.0, 50.0), (100, 100)))'
67		>>> _testrepr(splitLine((0, 0), (100, 100), 100, True))
68		'(((0, 0), (100, 100)))'
69		>>> _testrepr(splitLine((0, 0), (100, 100), 0, True))
70		'(((0, 0), (0.0, 0.0)), ((0.0, 0.0), (100, 100)))'
71		>>> _testrepr(splitLine((0, 0), (100, 100), 0, False))
72		'(((0, 0), (0.0, 0.0)), ((0.0, 0.0), (100, 100)))'
73	"""
74	pt1, pt2 = Numeric.array((pt1, pt2))
75	a = (pt2 - pt1)
76	b = pt1
77	ax = a[isHorizontal]
78	if ax == 0:
79		return [(pt1, pt2)]
80	t = float(where - b[isHorizontal]) / ax
81	if 0 <= t < 1:
82		midPt = a * t + b
83		return [(pt1, midPt), (midPt, pt2)]
84	else:
85		return [(pt1, pt2)]
86
87
88def splitQuadratic(pt1, pt2, pt3, where, isHorizontal):
89	"""Split the quadratic curve between pt1, pt2 and pt3 at position 'where',
90	which is an x coordinate if isHorizontal is False, a y coordinate if
91	isHorizontal is True. Return a list of curve segments.
92
93		>>> splitQuadratic((0, 0), (50, 100), (100, 0), 150, False)
94		[((0, 0), (50, 100), (100, 0))]
95		>>> _testrepr(splitQuadratic((0, 0), (50, 100), (100, 0), 50, False))
96		'(((0.0, 0.0), (25.0, 50.0), (50.0, 50.0)), ((50.0, 50.0), (75.0, 50.0), (100.0, 0.0)))'
97		>>> _testrepr(splitQuadratic((0, 0), (50, 100), (100, 0), 25, False))
98		'(((0.0, 0.0), (12.5, 25.0), (25.0, 37.5)), ((25.0, 37.5), (62.5, 75.0), (100.0, 0.0)))'
99		>>> _testrepr(splitQuadratic((0, 0), (50, 100), (100, 0), 25, True))
100		'(((0.0, 0.0), (7.32233047034, 14.6446609407), (14.6446609407, 25.0)), ((14.6446609407, 25.0), (50.0, 75.0), (85.3553390593, 25.0)), ((85.3553390593, 25.0), (92.6776695297, 14.6446609407), (100.0, -7.1054273576e-15)))'
101		>>> # XXX I'm not at all sure it the following behavior is desirable:
102		>>> _testrepr(splitQuadratic((0, 0), (50, 100), (100, 0), 50, True))
103		'(((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)))'
104	"""
105	a, b, c = calcQuadraticParameters(pt1, pt2, pt3)
106	solutions = solveQuadratic(a[isHorizontal], b[isHorizontal],
107		c[isHorizontal] - where)
108	solutions = [t for t in solutions if 0 <= t < 1]
109	solutions.sort()
110	if not solutions:
111		return [(pt1, pt2, pt3)]
112
113	segments = []
114	solutions.insert(0, 0.0)
115	solutions.append(1.0)
116	for i in range(len(solutions) - 1):
117		t1 = solutions[i]
118		t2 = solutions[i+1]
119		delta = (t2 - t1)
120		# calc new a, b and c
121		a1 = a * delta**2
122		b1 = (2*a*t1 + b) * delta
123		c1 = a*t1**2 + b*t1 + c
124		pt1, pt2, pt3 = calcQuadraticPoints(a1, b1, c1)
125		segments.append((pt1, pt2, pt3))
126	return segments
127
128
129def splitCubic(pt1, pt2, pt3, pt4, where, isHorizontal):
130	"""Split the cubic curve between pt1, pt2, pt3 and pt4 at position 'where',
131	which is an x coordinate if isHorizontal is False, a y coordinate if
132	isHorizontal is True. Return a list of curve segments.
133
134		>>> splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 150, False)
135		[((0, 0), (25, 100), (75, 100), (100, 0))]
136		>>> _testrepr(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 50, False))
137		'(((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0)), ((50.0, 75.0), (68.75, 75.0), (87.5, 50.0), (100.0, 0.0)))'
138		>>> _testrepr(splitCubic((0, 0), (25, 100), (75, 100), (100, 0), 25, True))
139		'(((0.0, 0.0), (2.2937927384, 9.17517095361), (4.79804488188, 17.5085042869), (7.47413641001, 25.0)), ((7.47413641001, 25.0), (31.2886200204, 91.6666666667), (68.7113799796, 91.6666666667), (92.52586359, 25.0)), ((92.52586359, 25.0), (95.2019551181, 17.5085042869), (97.7062072616, 9.17517095361), (100.0, 1.7763568394e-15)))'
140	"""
141	a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)
142	solutions = solveCubic(a[isHorizontal], b[isHorizontal], c[isHorizontal],
143		d[isHorizontal] - where)
144	solutions = [t for t in solutions if 0 <= t < 1]
145	solutions.sort()
146	if not solutions:
147		return [(pt1, pt2, pt3, pt4)]
148
149	segments = []
150	solutions.insert(0, 0.0)
151	solutions.append(1.0)
152	for i in range(len(solutions) - 1):
153		t1 = solutions[i]
154		t2 = solutions[i+1]
155		delta = (t2 - t1)
156		# calc new a, b, c and d
157		a1 = a * delta**3
158		b1 = (3*a*t1 + b) * delta**2
159		c1 = (2*b*t1 + c + 3*a*t1**2) * delta
160		d1 = a*t1**3 + b*t1**2 + c*t1 + d
161		pt1, pt2, pt3, pt4 = calcCubicPoints(a1, b1, c1, d1)
162		segments.append((pt1, pt2, pt3, pt4))
163	return segments
164
165
166#
167# Equation solvers.
168#
169
170from math import sqrt, acos, cos, pi
171
172
173def solveQuadratic(a, b, c,
174		sqrt=sqrt):
175	"""Solve a quadratic equation where a, b and c are real.
176	    a*x*x + b*x + c = 0
177	This function returns a list of roots. Note that the returned list
178	is neither guaranteed to be sorted nor to contain unique values!
179	"""
180	if a == 0.0:
181		if b == 0.0:
182			# We have a non-equation; therefore, we have no valid solution
183			roots = []
184		else:
185			# We have a linear equation with 1 root.
186			roots = [-c/b]
187	else:
188		# We have a true quadratic equation.  Apply the quadratic formula to find two roots.
189		DD = b*b - 4.0*a*c
190		if DD >= 0.0:
191			rDD = sqrt(DD)
192			roots = [(-b+rDD)/2.0/a, (-b-rDD)/2.0/a]
193		else:
194			# complex roots, ignore
195			roots = []
196	return roots
197
198
199def solveCubic(a, b, c, d,
200		abs=abs, pow=pow, sqrt=sqrt, cos=cos, acos=acos, pi=pi):
201	"""Solve a cubic equation where a, b, c and d are real.
202	    a*x*x*x + b*x*x + c*x + d = 0
203	This function returns a list of roots. Note that the returned list
204	is neither guaranteed to be sorted nor to contain unique values!
205	"""
206	#
207	# adapted from:
208	#   CUBIC.C - Solve a cubic polynomial
209	#   public domain by Ross Cottrell
210	# found at: http://www.strangecreations.com/library/snippets/Cubic.C
211	#
212	if abs(a) < 1e-6:
213		# don't just test for zero; for very small values of 'a' solveCubic()
214		# returns unreliable results, so we fall back to quad.
215		return solveQuadratic(b, c, d)
216	a = float(a)
217	a1 = b/a
218	a2 = c/a
219	a3 = d/a
220
221	Q = (a1*a1 - 3.0*a2)/9.0
222	R = (2.0*a1*a1*a1 - 9.0*a1*a2 + 27.0*a3)/54.0
223	R2_Q3 = R*R - Q*Q*Q
224
225	if R2_Q3 < 0:
226		theta = acos(R/sqrt(Q*Q*Q))
227		rQ2 = -2.0*sqrt(Q)
228		x0 = rQ2*cos(theta/3.0) - a1/3.0
229		x1 = rQ2*cos((theta+2.0*pi)/3.0) - a1/3.0
230		x2 = rQ2*cos((theta+4.0*pi)/3.0) - a1/3.0
231		return [x0, x1, x2]
232	else:
233		if Q == 0 and R == 0:
234			x = 0
235		else:
236			x = pow(sqrt(R2_Q3)+abs(R), 1/3.0)
237			x = x + Q/x
238		if R >= 0.0:
239			x = -x
240		x = x - a1/3.0
241		return [x]
242
243
244def calcQuadraticParameters(pt1, pt2, pt3):
245	pt1, pt2, pt3 = Numeric.array((pt1, pt2, pt3))
246	c = pt1
247	b = (pt2 - c) * 2.0
248	a = pt3 - c - b
249	return a, b, c
250
251
252def calcCubicParameters(pt1, pt2, pt3, pt4):
253	pt1, pt2, pt3, pt4 = Numeric.array((pt1, pt2, pt3, pt4))
254	d = pt1
255	c = (pt2 - d) * 3.0
256	b = (pt3 - pt2) * 3.0 - c
257	a = pt4 - d - c - b
258	return a, b, c, d
259
260
261def calcQuadraticPoints(a, b, c):
262	pt1 = c
263	pt2 = (b * 0.5) + c
264	pt3 = a + b + c
265	return pt1, pt2, pt3
266
267
268def calcCubicPoints(a, b, c, d):
269	pt1 = d
270	pt2 = (c / 3.0) + d
271	pt3 = (b + c) / 3.0 + pt2
272	pt4 = a + d + c + b
273	return pt1, pt2, pt3, pt4
274
275
276def _testrepr(obj):
277	"""
278		>>> _testrepr([1, [2, 3], [], [[2, [3, 4]]]])
279		'(1, (2, 3), (), ((2, (3, 4))))'
280	"""
281	try:
282		it = iter(obj)
283	except TypeError:
284		return str(obj)
285	else:
286		return "(%s)" % ", ".join([_testrepr(x) for x in it])
287
288
289if __name__ == "__main__":
290	import doctest
291	doctest.testmod()
292