Bezier curves and their acceleration

This tutorial illustrates how Bézier curves are generalized to manifolds and how to minimize their acceleration, i.e. how to get a curve that is as straight or as geodesic while fulfilling constraints

This example also illustrates how to apply the minimization on the corresponding PowerManifold manifold using a gradient_descent with ArmijoLinesearch.

Table of contents

• Setup
• de Casteljau algorithm on manifolds
• Composite Bézire curves
• Minimizing the acceleration of a Bézier curve
• Literature

Setup

We first initialize the necessary packages

using Manopt, Manifolds

and we define some colors from Paul Tol

using Colors
black = RGBA{Float64}(colorant"#000000")
TolVibrantBlue = RGBA{Float64}(colorant"#0077BB")
TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
TolVibrantMagenta = RGBA{Float64}(colorant"#EE3377")
TolVibrantCyan = RGBA{Float64}(colorant"#33BBEE")
TolVibrantTeal = RGBA{Float64}(colorant"#009988")

Then we load our data, see artificial_S2_composite_bezier_curve, a composite Bezier curve consisting of 3 segments on the Sphere

B = artificial_S2_composite_bezier_curve();

De Casteljau algorithm on manifolds

This curve can be evaluated using de Casteljau's algorithm^{[Casteljau1959][Casteljau1963]} named after Paul de Casteljau(*1930). To simplify the idea and understand this algorithm, we will first only look at the points of the first segment

M = Sphere(2)
b = B[2].pts

4-element Array{Array{Float64,1},1}:
 [0.0, -1.0, 0.0]
 [-0.5, -0.7071067811865476, -0.5]
 [-0.7071067811865476, -0.5, 0.5]
 [-1.0, 0.0, 0.0]

On Euclidean spaces Bézier curves of these $n=4$ so called control points like this segment yield polynomials of degree $3$. The resulting curve $\gamma: [0,1] \to \mathbb R^m$ is called Bezier curve or Bézier spline and is named after Piérre Bezier (1910–1999). They can be evaluated by the de Casteljau algorithm by evaluating line segments between points. While it is not easy to evaluate polynomials on a manifold, evaluating line segments generalizes to the evaluation of shortest_geodesics

We will illustrate this using these points $b=(b_1,b_2,b_3,b_4)$ on the Sphere $\mathbb S^2$. Let's evaliuate this at the point $t=\frac{1}{4}\in [0,1]$. We first compute

t = 0.66
pts1 = shortest_geodesic.(Ref(M), b[1:3], b[2:4], Ref(t))

3-element Array{Array{Float64,1},1}:
 [-0.35034218424621255, -0.8686315144381913, -0.35034218424621255]
 [-0.7309387480271131, -0.6615839039713686, 0.16743740513564936]
 [-0.9645574184577981, -0.1865864230028918, 0.1865864230028918]

We obtain 3 points on the geodesics connecting the control points. Repeating this again twice

pts2 = shortest_geodesic.(Ref(M), pts1[1:2], pts1[2:3], Ref(t))
p = shortest_geodesic(M, pts2[1], pts2[2], t)

3-element Array{Float64,1}:
 -0.8451739139029242
 -0.5204002573064785
  0.1219205784655057

we obtain the point on the Bézier curve $c(t)$. This procedure is illustrated in the following image:

From the control points (blue) and their geodesics, ont evaluation per geodesic yields three interims points (cyan), their two successive geodeics another two points (teal) and at its geodesic at $t=0.66$ we obtain the point on the curve.

In Manopt.jl, to evaluate a Bézier curve knowing its BezierSegment, use de_casteljau.

There are a few nice observations to make, that hold also for these Bézier curves on manifolds:

• The curve starts in the first controlpoint $b_0$ and ends in the last controlpoint $b_3$
• The tangent vector to the curve at the start $\dot c(0)$ is equal to $\log_{b_0}b_1 = \dot\gamma_{b_0,b_0}(0)$, where $\gamma_{a,b}$ denotes the shortest geodesic.
• The tangent vector to the curve at the end $\dot c(1)$ is equal to $-\log_{b_3}b_2 = -\dot\gamma_{b_3,b_2}(0) = \dot\gamma_{b_2,b_3}(1)$.
• the curve is differentiable.

For more details on these properties, see for example ^{[PopielNoakes2007]}.

Composite Bézier curves

With the properties of a single Bézier curve, also called Bézier segment, we can “stitch” curves together. Let $a_0,\ldots,a_n$ and $b_0,\ldots,b_m$ be two sets of controlpoints for the Bézier segments $c(t)$ and $d(t)$, respectively. We define the composite Bézier curve by $B(t) = \begin{cases} c(t) & \text{ if } 0\leq t < 1, \\ d(t-1) & \text{ if } 1\leq t \leq 2,\end{cases}$ where $t\in[0,2]$. This can of course be generalised straight forward to more than two cases. With the properties from the previous section we can now state that

• the curve $B(t)$ is continuous if $c(1)=d(0)$ or in other words $a_n=b_0$
• the curve $B(t)$ is differentiable if additionally $\dot c(1)=\dot d(0)$ or in other words $-\log_{a_n}a_{n-1} = \log_{b_0}b_1$. This is equivalent to $a_n=b_0 = \gamma_{a_{n-1}b_1}(\tfrac{1}{2})$.

One nice interpretation of the last characterization is, that the tangents $\log_{a_n}a_{n-1}$ and $\log_{b_0}b_1$ point into opposite directions. For a continuous curve, the first point of every segment (except for the first segment) can be ommitted, for a differentiable curve the first two points (except for the first segment) can be ommitted. You can reduce storage by calling get_bezier_points, though for econstruciton with get_bezier_segments you also need get_bezier_degrees. The reduced storage is represented as an array of points, i.e. an element of the corresponding PowerManifold.

For the three segment example from the beginning this looks as follows^[1]

Minimizing the acceleration of a composite Bézier curve

The motivation to minimize the acceleration of the composite Bézier curve is, that the curve should get “straighter” or more geodesic like. If we discretize the curve $B(t)$ with its control points denoted by $b_{i,j}$ for the $j$th note in the $i$th segment, the discretized model for equispaced $t_i$, $i=0,\ldots,N$ in the domain of $B$ reads^{[BergmannGousenbourger2018]}

where $\mathrm{d}_2$ denotes the second order finite difference using the mid point approach, see costTV2^{[BacakBergmannSteidlWeinmann2016]},

Another model is based on logarithmic maps, see ^{[BoumalAbsil2011]}, but that is not considered here. An advantage of the model considered here is, that it only consist of the evaluation of geodesics. This yields a gradient of $A(b)$ with respect to $b$ adjoint_Jacobi_fields. The following image shows the negative gradient (scaled)

In the following we consider two cases: Interpolation, which fixes the junction and end points of $B(t)$ and approximation, where a weight and a dataterm are additionally introduced.

Interpolation

For interpolation, the junction points are fixed and their gradient entries are hence set to zero. After transferring to the already mentioned PowerManifold, we can then perform a gradient_descent as follows

curve_samples = collect(range(0.0, 3.0; length=151)) #exactness of approximating d^2
pB = get_bezier_points(M, B, :differentiable)
N = PowerManifold(M, NestedPowerRepresentation(), length(pB))
F(pB) = cost_acceleration_bezier(M, pB, get_bezier_degrees(M, B), curve_samples)
∇F(pB) = ∇acceleration_bezier(M, pB, get_bezier_degrees(M, B), curve_samples)
x0 = pB
pB_opt_ip = gradient_descent(
    N,
    F,
    ∇F,
    x0;
    stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.5, 0.0001),
    stopping_criterion=StopWhenChangeLess(5 * 10.0^(-7)),
    debug=[
        :Iteration,
        " | ",
        :Cost,
        " | ",
        DebugGradientNorm(),
        " | ",
        DebugStepsize(),
        " | ",
        :Change,
        "\n",
        :Stop,
        10,
    ],
)

8-element Array{Array{Float64,1},1}:
 [0.0, 0.0, 1.0]
 [0.27844897623800974, -0.5873350022837311, 0.7599366833653981]
 [0.3880456421102258, -0.9052491185072907, 0.1730451186283273]
 [0.0, -1.0, 0.0]
 [-0.905249118507162, -0.388045642110108, 0.17304511862909577]
 [-1.0, 0.0, 0.0]
 [-0.5873350022830296, 0.27844897623774945, -0.7599366833660367]
 [0.0, 0.0, -1.0]

and the result looks like

Approximation

Similarly if we introduce the junction points as data fixed given $d_i$ and set (for simplicity) $p_i=b_{i,0}$ and $p_{n+1}=b_{n,4}$ and set $λ=3$ in

then $λ$ models how important closeness to the data $d_i$ is.

λ = 3.0
d = get_bezier_junctions(M, B)
F(pB) = cost_L2_acceleration_bezier(M, pB, get_bezier_degrees(M, B), curve_samples, λ, d)
∇F(pB) = ∇L2_acceleration_bezier(M, pB, get_bezier_degrees(M, B), curve_samples, λ, d)
x0 = pB
pB_opt_appr = gradient_descent(
    N,
    F,
    ∇F,
    x0;
    stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.5, 0.001),
    stopping_criterion=StopWhenChangeLess(10.0^(-5)),
    debug=[
        :Iteration,
        " | ",
        :Cost,
        " | ",
        DebugGradientNorm(),
        " | ",
        DebugStepsize(),
        " | ",
        :Change,
        "\n",
        :Stop,
        50,
    ],
)

8-element Array{Array{Float64,1},1}:
 [0.1837487271477433, -0.17275124768139632, 0.9676742280830204]
 [0.0642331672712227, -0.5426089221386995, 0.8375259147261502]
 [-0.0729969522949283, -0.8182992502159336, 0.5701383885090388]
 [-0.3363029443304467, -0.8761601552109021, 0.34531682851485007]
 [-0.8259917630675604, -0.5583393918895667, -0.07742564697145854]
 [-0.8762218316778143, -0.33656200927537233, -0.344907691424798]
 [-0.5428162372167487, 0.06395493934501625, -0.8374128601515666]
 [-0.17289300210425082, 0.1835835567412971, -0.9676802609940883]

and the result looks like

The role of $\lambda$ can be interpreted as follows: for large values of $\lambda$, the minimizer, i.e. the resulting curve, is closer to the original Bézier junction points. For small $\lambda$ the resting curve is closer to a geodesic and the control points are closer to the curve. For $\lambda=0$ any (not necessarily shortest) geodesic is a solution and the problem is ill-posed. To illustrate the effect of $\lambda$, the following image contains 1000 runs for $\lambda=10$ in dark currant to $\lambda=0.01$ in bright yellow.

The effect of the data term can also be seen in the following video, which starts a little slow and takes about 40 seconds.

Literature