Bézier 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 as possible while fulfilling constraints.

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

using Colors, PlutoUI, Manopt, Manifolds

We define some colors from Paul Tol

begin
    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")
end;

Let's also set up a few parameters

begin
    geo_pts = collect(range(0.0, 1.0; length=101))
    bezier_pts = collect(range(0.0, 3.0; length=201))
    camera_position = (-1.0, -0.7, 0.3)
    #render asy yes/no. If not, images included w/ markdown are assumed to be prerendered
    render_asy = false
    localpath = join(splitpath(@__FILE__)[1:(end - 1)], "/") # remove file to get files folder
    image_prefix = localpath * "/bezier"
    @info image_prefix
end;

We finally load our data, see artificial_S2_composite_bezier_curve, a composite Bézier curve consisting of 3 segments on the Sphere. The middle segment consists of the control points

begin
    B = artificial_S2_composite_bezier_curve()
    b = B[2].pts
end

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

Sphere(2, ℝ)

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 $γ: [0,1] → ℝ^m$ is called Bézier curve or Bézier spline and is named after Pierre Bézier (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 evaluate this at the point $t=\frac{1}{4}∈[0,1]$. We first compute

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

3-element Vector{Vector{Float64}}:
 [-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

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

3-element Vector{Float64}:
 -0.8451739139029242
 -0.5204002573064785
  0.1219205784655057

we obtain the point on the Bézier curve $c(t)$.

begin
    curves1 = [shortest_geodesic(M, b[i], b[i + 1], geo_pts) for i in 1:3]
    curves2 = [shortest_geodesic(M, pts1[i], pts1[i + 1], geo_pts) for i in 1:2]
    bezier_curve = [shortest_geodesic(M, pts2[1], pts2[2], geo_pts)]
end;

render_asy && begin
    asymptote_export_S2_signals(
        image_prefix * "/Casteljau-illustr.asy";
        points=[b, pts1, pts2, [p]],
        curves=[curves1..., curves2..., bezier_curve..., de_casteljau(M, B[2], geo_pts)],
        colors=Dict(
            :points =>
                [TolVibrantBlue, TolVibrantCyan, TolVibrantTeal, TolVibrantOrange],
            :curves => [
                TolVibrantBlue,
                TolVibrantBlue,
                TolVibrantBlue,
                TolVibrantCyan,
                TolVibrantCyan,
                TolVibrantTeal,
                black,
            ],
        ),
        dot_size=3.5,
        line_widths=[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.5],
        camera_position=camera_position,
    )
    render_asymptote(image_prefix * "/Casteljau-illustr.asy"; render=2)
end;

This image summarizes Casteljaus algorithm:

PlutoUI.LocalResource(image_prefix * "/Casteljau-illustr.png")

From the control points (blue) and their geodesics, one evaluation per geodesic yields three interim points (cyan), their two successive geodesics yield two other points (teal), and on the geodesic that connects these two 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 also hold 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 $\frac{1}{3}\log_{b_0}b_1 = \dot γ_{b_0,b_1}(0)$, where $γ_{a,b}$ denotes the shortest geodesic between $a$ and $b$.
The tangent vector to the curve at the end $\dot c(1)$ is equal to $-\frac{1}{3}\log_{b_3}b_2 = -γ_{b_3,b_2}(0) = \dot γ_{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,…,a_n$ and $b_0,…,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∈[0,2]$. This can of course be generalized straightforwardly 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 towards 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 a construction 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

render_asy && begin
    asymptote_export_S2_signals(
        image_prefix * "/Bezier-composite-curve.asy";
        curves=[de_casteljau(M, B, bezier_pts)],
        points=[get_bezier_junctions(M, B), get_bezier_inner_points(M, B)],
        tangent_vectors=[[
            Tuple(a) for a in zip(
                get_bezier_junctions(M, B, true),
                get_bezier_junction_tangent_vectors(M, B),
            )
        ]],
        colors=Dict(
            :curves => [black],
            :points => [TolVibrantBlue, TolVibrantTeal],
            :tvectors => [TolVibrantCyan],
        ),
        camera_position=camera_position,
        arrow_head_size=10.0,
        line_widths=[1.5, 1.5],
        dot_size=4.0,
    )
    render_asymptote(image_prefix * "/Bezier-composite-curve.asy"; render=2)
end;

Minimizing the Acceleration of a Composite Bézier Curve

The motivation to minimize the acceleration of a 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,…,N$ in the domain of $B$ reads[BergmannGousenbourger2018]

$$A(b) :eqq\sum_{i=1}^{N-1}\frac{\mathrm{d}^2_2 \bigl[ B(t_{i-1}), B(t_{i}), B(t_{i+1}) \bigr]}{\Delta_t^3},$$

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

$$d_2(x,y,z) := \min_{c ∈ \mathcal C_{x,z}} d_{\mathcal M}(c,y),\qquad x,y,z∈\mathcal M.$$

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)

render_asy && begin
    gradFullB = Manopt._grad_acceleration_bezier(
        M,
        get_bezier_points(M, B, :differentiable),
        [3, 3, 3],
        collect(range(0.0, 3.0; length=151)),
    )
    asymptote_export_S2_signals(
        image_prefix * "/Bezier-composite-curve-gradient.asy";
        curves=[de_casteljau(M, B, bezier_pts)],
        points=[get_bezier_junctions(M, B), get_bezier_inner_points(M, B)],
        tangent_vectors=[[
            Tuple(a) for a in zip(
                get_bezier_points(M, B, :continuous),
                -0.05 .* get_bezier_points(M, gradFullB, :continuous),
            )
        ]],
        colors=Dict(
            :curves => [black],
            :points => [TolVibrantBlue, TolVibrantTeal],
            :tvectors => [TolVibrantOrange],
        ),
        camera_position=camera_position,
        arrow_head_size=10.0,
        line_widths=[1.5, 1.5],
        dot_size=4.0,
    )
    render_asymptote(image_prefix * "/Bezier-composite-curve-gradient.asy"; render=2)
end;

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 data term are additionally introduced.

Interpolation

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

begin
    #number of sampling points refers to exactness of approximating d^2
    curve_samples = collect(range(0.0, 3.0; length=151))
    pB = get_bezier_points(M, B, :differentiable)
    N = PowerManifold(M, NestedPowerRepresentation(), length(pB))
    function F(M, pB)
        return cost_acceleration_bezier(
            M.manifold, pB, get_bezier_degrees(M.manifold, B), curve_samples
        )
    end
    function gradF(M, pB)
        return grad_acceleration_bezier(
            M.manifold, pB, get_bezier_degrees(M.manifold, B), curve_samples
        )
    end
    x0 = get_bezier_points(M, B, :differentiable)
    pB_opt_ip = gradient_descent(
        N,
        F,
        gradF,
        x0;
        stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.5, 0.0001),
        stopping_criterion=StopWhenChangeLess(5 * 10.0^(-7)),
    )
end;

and the result looks like

render_asy && begin
    B_opt_ip = get_bezier_segments(M, pB_opt_ip, get_bezier_degrees(M, B), :differentiable)
    asymptote_export_S2_signals(
        image_prefix * "/Bezier-IP-Min.asy";
        curves=[de_casteljau(M, B_opt_ip, bezier_pts), de_casteljau(M, B, bezier_pts)],
        points=[get_bezier_junctions(M, B_opt_ip), get_bezier_inner_points(M, B_opt_ip)],
        tangent_vectors=[[
            Tuple(a) for a in zip(
                get_bezier_junctions(M, B_opt_ip, true),
                get_bezier_junction_tangent_vectors(M, B_opt_ip),
            )
        ]],
        colors=Dict(
            :curves => [TolVibrantBlue, black],
            :points => [TolVibrantBlue, TolVibrantTeal],
            :tvectors => [TolVibrantCyan],
        ),
        camera_position=camera_position,
        arrow_head_size=10.0,
        line_widths=[1.5, 0.75, 1.5],
        dot_size=4.0,
    )
    render_asymptote(image_prefix * "/Bezier-IP-Min.asy"; render=2)
end;

Where the original curve is shown in black and the interpolating curve with minimized (discretized) acceleration is shown in blue including its junction points (also blue), tangent vectors (light blue) and control points (teal).

Approximation

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

$$\frac{λ}{2}\sum_{k=0}^3 d_{\mathcal M}(d_i,p_i)^2 + A(b),$$

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

Then we obtain the following code to minimize the acceleration while approximating the original curve

begin
    λ = 3.0
    d = get_bezier_junctions(M, B)
    function F2(M, pB)
        return cost_L2_acceleration_bezier(
            M.manifold, pB, get_bezier_degrees(M.manifold, B), curve_samples, λ, d
        )
    end
    function gradF2(M, pB)
        return grad_L2_acceleration_bezier(
            M.manifold, pB, get_bezier_degrees(M.manifold, B), curve_samples, λ, d
        )
    end
    x1 = get_bezier_points(M, B, :differentiable)
    pB_opt_appr = gradient_descent(
        N,
        F2,
        gradF2,
        x1;
        stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.5, 0.001),
        stopping_criterion=StopWhenChangeLess(10.0^(-5)),
    )
end

8-element Vector{Vector{Float64}}:
 [0.1839506847929026, -0.17268149347510886, 0.9676483076900453]
 [0.06449579457054468, -0.5425499741454511, 0.8375439200648189]
 [-0.07275785048325141, -0.818277457950548, 0.570200225362126]
 [-0.3361255397939874, -0.8761777330856455, 0.34544493272159427]
 [-0.8260046134858585, -0.5583398456161249, -0.07728515573798504]
 [-0.8762842719021493, -0.33653156307773063, -0.3447787433548887]
 [-0.5428643794614353, 0.06409564980553185, -0.8373708934443852]
 [-0.17289249065987983, 0.18371074551323974, -0.9676562140845296]

render_asy && begin
    B_opt_appr = get_bezier_segments(M, pB_opt_appr, get_bezier_degrees(M, B), :differentiable)
    asymptote_export_S2_signals(
        image_prefix * "/Bezier-Appr-Min.asy";
        curves=[de_casteljau(M, B_opt_appr, bezier_pts), de_casteljau(M, B, bezier_pts)],
        points=[
            d,
            get_bezier_junctions(M, B_opt_appr),
            get_bezier_inner_points(M, B_opt_appr),
        ],
        tangent_vectors=[[
            Tuple(a) for a in zip(
                get_bezier_junctions(M, B_opt_appr, true),
                get_bezier_junction_tangent_vectors(M, B_opt_appr),
            )
        ]],
        colors=Dict(
            :curves => [TolVibrantBlue, black],
            :points => [TolVibrantOrange, TolVibrantBlue, TolVibrantTeal],
            :tvectors => [TolVibrantCyan],
        ),
        camera_position=camera_position,
        arrow_head_size=10.0,
        line_widths=[1.5, 0.75, 1.5],
        dot_size=4.0,
    )
    render_asymptote(image_prefix * "/Bezier-Appr-Min.asy"; render=2)
end;

PlutoUI.LocalResource(image_prefix * "/Bezier-Appr-Min.png")

Additionally to the last image, the data points $d_i$ (junction points of the original curve) are shown in orange. The distance between these and the blue junction points is part of the cost function here.

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

Approximation min Acc

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

Video of the effect of lambda, the weight of the dataterm

Literature

BacakBergmannSteidlWeinmann2016

Bačák, M., Bergmann, R., Steidl, G. and Weinmann, A.: A second order nonsmooth variational model for restoring manifold-valued images, SIAM Journal on Scientific Computations, Volume 38, Number 1, pp. A567–597, doi: 10.1137/15M101988X, arXiv: 1506.02409

BergmannGousenbourger2018

Bergmann, R. and Gousenbourger, P.-Y.: A variational model for data fitting on manifolds by minimizing the acceleration of a Bézier curve. Frontiers in Applied Mathematics and Statistics, 2018. doi: 10.3389/fams.2018.00059, arXiv: 1807.10090

BoumalAbsil2011

Boumal, N. and Absil, P.-A.: A discrete regression method on manifolds and its application to data on SO(n). In: IFAC Proceedings Volumes (IFAC-PapersOnline). Vol. 18. Milano (2011). p. 2284–89. doi: 10.3182/20110828-6-IT-1002.00542, web: www

Casteljau1959

de Casteljau, P.: Outillage methodes calcul, Enveloppe Soleau 40.040 (1959), Institute National de la Propriété Industrielle, Paris.

Casteljau1963

de Casteljau, P.: Courbes et surfaces à pôles, Microfiche P 4147-1, André Citroën Automobile SA, Paris, (1963).

PopielNoakes2007

Popiel, T. and Noakes, L.: Bézier curves and $C^2$ interpolation in Riemannian manifolds. Journal of Approximation Theory (2007), 148(2), pp. 111–127.- doi: 10.1016/j.jat.2007.03.002.