To the Math Masters: Looking for inverse function of SplitCubicAtT + Pythagoreic Triangle

yanone's picture

Dear Bezier and Math Masers,

for a piece of software I'm writing for some time now I need to calculate points on a Bezier curve. These points need to be of exactly identical absolute distance from each other.
At the moment I'm using the common SplitCubicAtT function (T of which is a float between 0 and 1) to part the curve into two with T being 0.5 at first, then calculate the distance between P1 and P1,4 using the Pythagoreic Triangle. Then I repeat the calculation with increasing or decreasing T until the difference in distance is below a certain threshold (= binary search). But this method is 1) slow and 2) not precise.

Step 1: SplitCubicAtT, result is (among other points) the point P1,4. The difference of position in x and y direction between P1,4 and P1 is used in
Step 2: to calculate the absolute distance between P1 and P1,4 using the Pythagoreic Triangle.

From what I understand there should be a possibility to invert the combination of these two funtions so one could hand over the absolute distance and receive the positions of the P1,4.

From what I can see from here this classifies as rocket science, but I thought that someone here has already stepped across this problem earlier and might have found a solution.

As SplitCubicAtT I'm using a simple interpolation function

def SplitCubicAtT(p1, p2, p3, p4, t):
	Split cubic Beziers curve at relative value t, return the two resulting segments.
	from ynlib.maths import Interpolate
	p12 = (Interpolate(p1[0], p2[0], t), Interpolate(p1[1], p2[1], t))
	p23 = (Interpolate(p2[0], p3[0], t), Interpolate(p2[1], p3[1], t))
	p34 = (Interpolate(p3[0], p4[0], t), Interpolate(p3[1], p4[1], t))

	p123 = (Interpolate(p12[0], p23[0], t), Interpolate(p12[1], p23[1], t))
	p234 = (Interpolate(p23[0], p34[0], t), Interpolate(p23[1], p34[1], t))

	p1234 = (Interpolate(p123[0], p234[0], t), Interpolate(p123[1], p234[1], t))

	return (p1, p12, p123, p1234), (p1234, p234, p34, p4)

def SameLengthSegments(segment, distance, precision, firstpoint = None):
	Finds points on a curve segment with equal distance (approximated through binary search, with given precision).

	If firstpoint is given, that would in most cases be the second last calculated point of the previous segment
	(to avoid gaps between smooth connection segments), this point is used as the starting point instead of p1.
	The distance from firstpoint to p1 should then be less than 'distance'.
	Returns a list with calculated points and the position of the last calculated point.
	from ynlib.maths import Distance
	from ynlib.beziers import SplitCubicAtT

	points = []
	p1, p2, p3, p4 = segment
	l = distance
	t = None

	segments = SplitCubicAtT(p1, p2, p3, p4, .5)

	# Use firstpoint	
	firstrun = True
	if firstrun and firstpoint != None:
		d = Distance(firstpoint, segments[0][3])
		d = Distance(p1, segments[0][3])

	count = 0

	while Distance(segments[1][0], p4) > l:

		min = 0
		max = 1
#		if t != None:
#			min = t
		t = min + (max - min) / 2.0

		segments = SplitCubicAtT(p1, p2, p3, p4, t)
		if firstrun and firstpoint != None:
			d = Distance(firstpoint, segments[0][3])
			d = Distance(p1, segments[0][3])

		# Binary search
		while (d - l) > precision or (d - l) < (precision * -1):
			if (d-l) > 0:
				max = t
			elif (d-l) < 0:
				min = t
			t = min + (max - min) / 2.0

			segments = SplitCubicAtT(p1, p2, p3, p4, t)

			# Use last point of previous curve as first point
			if firstrun and firstpoint != None:
				d = Distance(firstpoint, segments[0][3])
				d = Distance(segments[0][0], segments[0][3])
			count += 1
		p1 = segments[1][0]
		p2 = segments[1][1]
		p3 = segments[1][2]

		firstrun = False

	# List of points excluding, last point
	return points, segment[3], count

def Distance(p1, p2):
	Return distance between two points definded as (x, y).
	import math
	return math.sqrt( (p2[0] - p1[0]) ** 2 + (p2[1] - p1[1]) ** 2 )

instead of the probably more correct (but much slower) equation found in fontTools (results are identical):

def splitCubicAtT(pt1, pt2, pt3, pt4, *ts):
	"""Split the cubic curve between pt1, pt2, pt3 and pt4 at one or more
	values of t. Return a list of curve segments.

		>>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5))
		((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))
		>>> printSegments(splitCubicAtT((0, 0), (25, 100), (75, 100), (100, 0), 0.5, 0.75))
		((0.0, 0.0), (12.5, 50.0), (31.25, 75.0), (50.0, 75.0))
		((50.0, 75.0), (59.375, 75.0), (68.75, 68.75), (77.34375, 56.25))
		((77.34375, 56.25), (85.9375, 43.75), (93.75, 25.0), (100.0, 0.0))
	a, b, c, d = calcCubicParameters(pt1, pt2, pt3, pt4)
	return _splitCubicAtT(a, b, c, d, *ts)

def _splitCubicAtT(a, b, c, d, *ts):
	ts = list(ts)
	ts.insert(0, 0.0)
	segments = []
	for i in range(len(ts) - 1):
		t1 = ts[i]
		t2 = ts[i+1]
		delta = (t2 - t1)
		# calc new a, b, c and d
		a1 = a * delta**3
		b1 = (3*a*t1 + b) * delta**2
		c1 = (2*b*t1 + c + 3*a*t1**2) * delta
		d1 = a*t1**3 + b*t1**2 + c*t1 + d
		pt1, pt2, pt3, pt4 = calcCubicPoints(a1, b1, c1, d1)
		segments.append((pt1, pt2, pt3, pt4))
	return segments

def calcQuadraticParameters(pt1, pt2, pt3):
	pt1, pt2, pt3 = numpy.array((pt1, pt2, pt3))
	c = pt1
	b = (pt2 - c) * 2.0
	a = pt3 - c - b
	return a, b, c

def calcCubicPoints(a, b, c, d):
	pt1 = d
	pt2 = (c / 3.0) + d
	pt3 = (b + c) / 3.0 + pt2
	pt4 = a + d + c + b
	return pt1, pt2, pt3, pt4

curve.png5.76 KB
hrant's picture

Just ask Raph Levien.


miha's picture

You can ask the question on math forum: Basically you have to solve the equation "absolute distance(P1, P1.4) = constant" and express P1.4 point with Bezier equation with one variable t.
In case of unusual bezier segments there may be more (maybe even up to 4?) resulting points with equal distance to the P1.

butterick's picture

From what I understand there should be a possibility to invert the combination of these two funtions so one could hand over the absolute distance and receive the positions of the P1,4.

The set of all points of distance K from point P1 is a circle of radius K. If you compute the intersection of that circle with bezier P1234, that's your point.

sgh's picture

The cubic Bezier can be written as a vector-valued function (x(t),y(t)) parameterized by the single variable t for 0<=t<=1. You want to find the value of t such that the distance from (x(t),y(t)) to P1 is a fixed value. Thus, (as Miha said) you wish to solve the equation sqrt((x(t)-P1.x)^2+(y(t)-P1.y)^2)=constant (*). Moving the constant to the left side, you want to find when the function f(t)=sqrt((x(t)-P1.x)^2+(y(t)-P1.y)^2)-constant is equal to 0 for 0<=t<=1. This can be efficiently solved using Newton's Method; this will be significantly faster than the binary search approach you described above since it uses information from the derivative. Newton's Method is iterative but converges quickly, so you can quickly compute the answer to the level of precision that you need.

One trick to simplify applying Newton's Method is to square both sides of equation (*) first. This will give the same t value, but then the function f(t) becomes a polynomial in t. This is easier to differentiate. The resulting polynomial has degree 6, and so there is no closed form solution (in general) for the roots. However, Newton's Method does converge quadratically.

I'm curious as to what you are using this for. This is not one of the typical computations done with Bezier curves.

Best wishes,

raph's picture

@hrant: heh, thanks.

@sgh: You are absolutely correct, the convergence of Newton's method is much faster than bisection. However, Newton's method can fail to converge when the relationship is not approximately linear, which can definitely happen with cubic Bezier segments when the control points are "kinky." Bisection is much more robust in that case, as it's pretty much guaranteed to converge at _something_ in range of the solution.

One approach which might be more robust than Newton but converge faster than bisection is the secant method. Actually this has the added advantage that you don't need to analytically compute the derivatives, as well. I've used it quite a bit in my spline foo.

There are other approaches to the problem which could work, but they could get hairy. If you can decompose your beziers into circular arcs, then the problem reduces to computing the intersection between two circles, which is quite simple. Here's a reference on that:

And, what @sgh said, it's hard for me to imagine why you really need this.

Hope this helps.

hrant's picture

QED. ;-)

But it's great to see we have at least a few more
math experts around here. I myself have a minor in
Numerical Analysis, but there's now more rust than metal...


yanone's picture

Hi people,

thanks four your responses. Understanding them is already hairy for me.
So let me continue with answering the question to the why.
I'm working on a type design software plugin that illustrates curve speed on top of bezier outlines, live while editing outlines.
In order to calculate the amount of curvature I'm calculating on-curve points in equal distance, then calculating the angle between them. Just cutting a curve segment into equal parts using SplitCubicAtT doesn't work because distance and hence angle will change when a curve becomes tighter, so angle calculation is useless. They need to be in equal absolute distance from each other.

But, thinking about it, i just had another idea. When distance between on-curve points (returned by SplitCubicAtT) decreases with increasing curvature, there's already the numerical basis for my illustration, right?
Here's another problem, though: This might work only per segment. How to sync illustrations between segments, when an equal amount of splits per segment will always return different distances, since segment lengths differ?

hrant's picture

Cool idea. Wouldn't the first derivative of
the equation give you this relatively easily?

BTW, related and perhaps helpful:
Do I remember correctly that there's a certain type
of spline equation that ensures constant (or maybe
controllable) velocity of change?


yanone's picture

Dear Bezier masters,

I've tried to approach my problem (calculating curve speed) in a more general way, using the general cubic Bezier equation.
I've constructed the first (and to be sure, also the second) derivative, as suggested. Please correct me, if I made mistakes there. School was a looong time ago.
The function solveCubic return as the first argument p1234 (identical to the middle on-curve point calculated by the known splitCubicAtT function, I just skipped the other points for now), then the result of the first derivative, then the second.
All are vectors (or coordinates).

The segment used for the results posted below (with t as 0.1 steps) is a quarter of a circle, or a quarter of as close a circle as one can construct with Beziers. The first derivative should then have all equal (more or less) values, right, in case of a circle? In other words, curve speed should be steady at all positions in a circle.
But the results differ by 8%, even more in the second derivative. This doesn't look right. I've double-checked for a correct circle pie piece.

But maybe I have a misunderstanding of what to do with the results of the equations that are actually x/y-coordinates. So far I calculate the absolute value of the vectors (rightmost in the results, respectively).

import numpy, math

def solveCubic(p0, p1, p2, p3, t):
	p0, p1, p2, p3 = numpy.array((p0, p1, p2, p3))

	a = -p0 + 3.0 * p1 - 3.0 * p2 + p3
	b = 3.0 * p0 - 6.0 * p1 + 3.0 * p2
	c = -3.0 * p0 + 3.0 * p1
	d = p0
	# Cubic Bezier
	f = a*t**3 + b*t**2 + c*t + d
	# First derivative
	f1 = 3*a*t**2 + 2*b*t + c
	# Second derivative
	f2 = 6*a*t + 2*b
	return f, f1, f2

g = CurrentGlyph()
p1 = (g[0][0].points[0].x, g[0][0].points[0].y)
p2 = (g[0][1].points[0].x, g[0][1].points[0].y)
p3 = (g[0][1].points[1].x, g[0][1].points[1].y)
p4 = (g[0][1].points[2].x, g[0][1].points[2].y)

for t in range(11):

    f, f1, f2 = solveCubic(p1, p2, p3, p4, t / 10.0)

    print "t=%s" % (t/10.0), " -  f1=%s=abs %s" % (str(f1).ljust(17), int(math.sqrt( f1[0]**2 + f1[1]**2 ))), " -  f2=%s=abs %s" % (str(f2).ljust(13), int(math.sqrt( f2[0]**2 + f2[1]**2 )))

The results are:

t=0.0  -  f1=[-417.    0.]    =abs 417  -  f2=[ 162. -672.]=abs 691
t=0.1  -  f1=[-398.25  -64.65]=abs 403  -  f2=[ 213. -621.]=abs 656
t=0.2  -  f1=[-374.4 -124.2]  =abs 394  -  f2=[ 264. -570.]=abs 628
t=0.3  -  f1=[-345.45 -178.65]=abs 388  -  f2=[ 315. -519.]=abs 607
t=0.4  -  f1=[-311.4 -228. ]  =abs 385  -  f2=[ 366. -468.]=abs 594
t=0.5  -  f1=[-272.25 -272.25]=abs 385  -  f2=[ 417. -417.]=abs 589
t=0.6  -  f1=[-228.  -311.4]  =abs 385  -  f2=[ 468. -366.]=abs 594
t=0.7  -  f1=[-178.65 -345.45]=abs 388  -  f2=[ 519. -315.]=abs 607
t=0.8  -  f1=[-124.2 -374.4]  =abs 394  -  f2=[ 570. -264.]=abs 628
t=0.9  -  f1=[ -64.65 -398.25]=abs 403  -  f2=[ 621. -213.]=abs 656
t=1.0  -  f1=[   0. -417.]    =abs 417  -  f2=[ 672. -162.]=abs 691
miha's picture

I think the first (or second) derivative is conceptually wrong answer and what you are looking for is “curvature”, which you already compute through its “natural” definition … but you said it’s slow. There are also other types of formulas for computing curvature, but others might write more detailed answers.
When I saw your image I remembered this fine illustration by Tim Ahrens: What constitutes a "bad curve"?

yanone's picture

I found it. It's an equation that involves the first and second derivative of the general Bezier equation.
See it in action:

Syndicate content Syndicate content