Movement along the curve with constant speed

Alexey Karamyshev
8 min readAug 3, 2021

--

While developing a game in quite varied genre, there may be a need to move some game object along a smooth curve with a constant or controlled speed, whether it is a truck traveling from city A to city B, a rocket fired along a tricky trajectory, or an enemy plane performing a maneuver.

Probably everyone related to the topic knows or, at least, has heard about Bezier curves, B-splines, Hermite splines and other interpolation and smoothing splines and would quite correctly suggest using one of them in the described situation, but not everything is as simple as we would like it to be.

In our case, a spline is a function that maps a set of input parameters (control points) and the value of the argument t (usually taking values from 0 to 1) to a point on a plane or in space. The resulting curve is the set of values of the spline function at t[0,1].

As an example, we can consider the cubic Bezier curve, which is given by the following equation:

Cubic Bezier curve

The figure shows two cubic Bezier curves defined by four points (the curve passes through two of them, the remaining two set the curvature).

Animation of the parameter to the curve point mapping

In order to build a more complex and variable route from the sections of curves, they can be connected along a chain, leaving one common point as a link between them:

Piecewise spline

With the parameterization of the route sorted out, now we turn to the main question. In order for our conditional plane to move along the route at a constant speed, we need to be able to calculate a point on the curve at any time depending on the distance traveled along this curve s(len), while only having the ability to calculate the position of a point on the curve by the value of the parameter t (the y(t) function). It is at this stage that the difficulties begin.

The first attempt is to make a linear mapping len[0, Length] t[0,1] and calculate y(t) from the resulting value. It’s easy, computationally cheap, I’d say, exactly what is needed.

The problem of this method is obvious— in fact, the traveled distance s depends on t non-linearly, and it is enough to place points evenly distributed along the route by the value of t to make sure of it:

“Evenly” distributed points along the path

The plane following this path will slow down in some areas and accelerate in others, which makes this method of parameterizing the curve completely inapplicable for solving the described problem (the plane was chosen solely as a visual example and the goal was not to describe its movement in a physically correct way):

Visualization of incorrect parameterization of the curve

Using the search engines, stackoverflow and youtube, I found a second way to calculate s(len) = g(t): via the curve representation as a piecewise linear function (by calculating a set of points equidistant from each other along the curve):

Representation of a curve as a piecewise linear spline

This procedure is iterative: taking a small step of the parameter t, we move along the curve, summing the distance traveled as the length of a piecewise linear spline until the specified distance is passed, after which the point is placed, and the process resumes.

The code example:

public Vector2[] CalculateEvenlySpacedPoints(float spacing, float resolution = 1)
{
List<Vector2> evenlySpacedPoints = new List<Vector2>();
evenlySpacedPoints.Add(points[0]);
Vector2 previousPoint = points[0];
float dstSinceLastEvenPoint = 0;
for (int segIdx = 0; segIdx < NumSegments; segIdx++)
{
Vector2[] p = GetPointsInSegment(segIdx);
float controlNetLength =
Vector2.Distance(p[0], p[1])
+ Vector2.Distance(p[1], p[2])
+ Vector2.Distance(p[2], p[3]);
float estimatedCurveLength =
Vector2.Distance(p[0], p[3]) + controlNetLength / 2f;
int divisions = Mathf.CeilToInt(estimatedCurveLength * resolution * 10);
float t = 0;
while (t <= 1)
{
t += 1f / divisions;
Vector2 pointOnCurve = Bezier.EvaluateCubic(p[0], p[1], p[2], p[3], t);
dstSinceLastEvenPoint += Vector2.Distance(previousPoint, pointOnCurve);
while (dstSinceLastEvenPoint >= spacing)
{
float overshootDst = dstSinceLastEvenPoint - spacing;
Vector2 newEvenlySpacedPoint = pointOnCurve + (previousPoint - pointOnCurve).normalized * overshootDst;
evenlySpacedPoints.Add(newEvenlySpacedPoint);
dstSinceLastEvenPoint = overshootDst;
previousPoint = newEvenlySpacedPoint;
}

previousPoint = pointOnCurve;
}
}

return evenlySpacedPoints.ToArray();
}

It seems like the problem is solved, because you can divide the route into a lot of segments and rejoice at how smoothly and steadily the object moves, since calculating a point on a piecewise linear function is a simple and fast. But this approach also has quite obvious disadvantages that bothered me — a rather expensive procedure for changing the division step or the curve geometry and the need to find a balance between accuracy (coming with more memory consumption) and memory savings (the route fractures become more noticeable with less divisions).

As a result, I continued my search and came across an excellent article “Moving Along a Curve with Specified Speed”, on which further arguments are based.

The velocity value can be calculated as σ(t) = |V(t)|=|dX/dt|, where X(t) is the spline function. To solve the problem, we need to find the function Y(t)=X(s), where s is the length of the arc from the starting point to the desired one, and this expression formalizes the relationship between t and s. We may apply the chain rule from Calculus to obtain

Given that the speed is a non-negative value, we get

since |dX/ds| = 1 by the condition of the velocity vector modulus immutability (generally speaking, |dX/ds| is not equal to one but to a constant, but for simplicity of calculations we assume this constant to be equal to one).

Now we get the relationship between t and s explicitly:

from where the total length of the curve L is obviously calculated as g(1). Having the resulting formula and the value of the argument t it is possible to calculate the arc length to the point at which value t is mapped.

But we are interested in the inverse operation, that is, in calculating the value of the argument t for a given arc length s:

Since there is no general analytical way to find the inverse function, we will have to look for a numerical solution. Define F(t)=g(t)-s. Given a value s, the problem is now to find a value t so that F(t) = 0. This is a root-finding problem that we might try solving using Newton’s method. Newton’s method produces a sequence

where

To calculate F(t) it is required to calculate g(t), the calculation of which, in turn, requires numerical integration.

The choice of t0=s/L as an initial approximation now looks reasonable (recall the first approach to solving the problem).

Next, we perform iterations using the Newton’s method until the error of the solution becomes acceptable, or the number of iterations performed becomes unacceptably large.

There is a potential problem when using only Newton’s method. The F(t) is a nondecreasing function because its derivative F’(t)=|dY/dt| is nonnegative. If the second derivative F”(t)>0, then the function is said to be convex and the Newton iterates are guaranteed to converge to the root. However, in our case, F”(t) may turn out to be negative, which is why the Newton iterates may converge outside the domain t[0,1]. The author of the article suggests using a modification of the Newton method, which is a hybrid of Newton’s method and bisection. If the next result of the iteration by the Newton’s method is not inside the current root-bounding interval, the midpoint of the interval is used instead. Regardless of whether a Newton iterate or bisection is used, the root-bounding interval is updated, so the interval length strictly decreases over time.

We also need to choose the numerical integration method — I used the quadratures of the Gauss method, since it provides a sufficiently accurate results at a low cost.

The code of the function that calculates t(s):

public float ArcLength(float t) => Integrate(x => TangentMagnitude(x), 0, t);

private float Parameter(float length)
{
float t = 0 + length / ArcLength(1);
float lowerBound = 0f;
float upperBound = 1f;

for (int i = 0; i < 100; ++i)
{
float f = ArcLength(t) - length;

if (Mathf.Abs(f) < 0.01f)
break;

float derivative = TangentMagnitude(t);
float candidateT = t - f / derivative;

if (f > 0)
{
upperBound = t;
if (candidateT <= 0)
t = (upperBound + lowerBound) / 2;
else
t = candidateT;
}
else
{
lowerBound = t;
if (candidateT >= 1)
t = (upperBound + lowerBound) / 2;
else
t = candidateT;
}
}
return t;
}

The numerical integration function:

private static readonly (float, float)[] CubicQuadrature =
{(-0.7745966F, 0.5555556F), (0, 0.8888889F), (0.7745966F, 0.5555556F)};

public static float Integrate(Func<float, float> f, in float lowerBound, in float uppedBound)
{
var sum = 0f;
foreach (var (arg, weight) in CubicQuadrature)
{
var t = Mathf.Lerp(lowerBound, uppedBound, Mathf.InverseLerp(-1, 1, arg));
sum += weight * f(t);
}

return sum * (upperBound - lowerBound) / 2;
}

Further we can adjust the error of the Newton’s method or choose a more accurate method of numerical integration but the problem is in fact solved — an algorithm for calculating the argument t of the spline for a given arc length was found. To combine several sections of curves into single path and write a generalized procedure for calculating s(t) you can store the lengths of all sections and pre-calculate the index of the section where you want to perform an iterative procedure using the modified Newton’s method.

Points evenly distributed along the path
Now the plane is moving evenly

Thus, we have considered several ways to parametrize the spline by the distance traveled. In the source article the author offers to solve the differential equation numerically as an option, but, according to his own remark, he prefers the modified Newton method:

I prefer the Newton’s-method approach because generally it can be computed faster than with differential equation solvers.

As a conclusion, this is a time table for calculating the position on the curve shown in the screenshots in one thread of the i5–9600K processor:

+------------------------+------------------+
| Number of calculations | Average time, ms |
+------------------------+------------------+
| 100 | 0.62 |
| 1000 | 6.24 |
| 10000 | 69.01 |
| 100000 | 672.81 |
+------------------------+------------------+

I think that such a speed of calculations makes it possible to use them in various games and simulations. I will be glad to see corrections and criticism in the comments.

--

--