First derivative of a Bézier curve (full development)

Development of the velocity of a Bézier curve. - 08/2023 - #Jumble

$ % This version will make mathjax ignore the width of the annotation and keep the main % formula centered. Obviously this mean the annotation is forcefully placed to the right % and may bleed out: \newcommand{Annotation}[1]{ \rlap{~~~~~~\raise 0.8em {\color{grey}\scriptsize{ #1 }}} } \newcommand{annotation}[1]{ {~~~~~~\raise 0.8em {\color{grey}\scriptsize{ #1 }}} } \newcommand{red}[1]{ { \color{red} #1 } } \newcommand{purple}[1]{ { \color{purple} #1 } } \newcommand{green}[1]{ { \color{green} #1 } } $

Here I detail the algebraic manipulation leading to the analytical expression of the first derivative (a.k.a velocity $\vec v(u)$) of a Bézier curve. I also provide 2nd derivative and general nth derivative formula for any Bézier curve.

Javascript applet

Here some ThreeJs (a Javascript / WebGL library) code demonstrating a simple cubic Bézier curve (degree 3, 4 control points). As an exercise the code could be extend to handle general Bernstein polynomial representation and compute a Frenet frame using the 1st and 2nd derivative presented below. You may play with the code directly into code pen

import * as THREE from "https://cdn.jsdelivr.net/npm/three@0.114/build/three.module.js";
import {
    OrbitControls
}
from "https://cdn.jsdelivr.net/npm/three@0.114/examples/jsm/controls/OrbitControls.js";
import {
    DragControls
}
from "https://cdn.jsdelivr.net/npm/three@0.114/examples/jsm/controls/DragControls.js";
"use strict";

/*
up axis: y
2d plane coordinates: x,z
 */

// for debugging purposes we allow switching backto perspective:
let g_orthographic = true;
var g_camera, g_scene, g_renderer;
var g_controls;

const g_frustumSize = 50;
let g_cam_params = {
    aspect_ratio() {
        return window.innerWidth / window.innerHeight;
    },
    left() {
        return g_frustumSize * this.aspect_ratio() * (-0.5);
    },
    right() {
        return g_frustumSize * this.aspect_ratio() * 0.5;
    },
    top() {
        return g_frustumSize * 0.5;
    },
    bottom() {
        return g_frustumSize * (-0.5);
    },
    near: 1,
    far: 1000,
};

// ---------------------------------------------------------------
let g_bezier = {
    ctrlPointInitList: [
        new THREE.Vector3(-5, -5, 0),
        new THREE.Vector3(-2, 5, 0),
        new THREE.Vector3(2, 5, 0),
        new THREE.Vector3(5, -5, 0),
    ],

    ctrlPointMeshes:[],
  
    ctrlPoint(i) {
        return this.ctrlPointMeshes[i].position.clone();
    },

    // Add control points to scene
    addCtrlPoints(scene, draggableList) {
        const pointGeometry = new THREE.CircleGeometry(0.3 /*radius*/, 32);
        const pointMaterial = new THREE.MeshBasicMaterial({
            color: 0xaa0000
        });
        for (const p of this.ctrlPointInitList) {            
            const point = new THREE.Mesh(pointGeometry, pointMaterial);
            this.ctrlPointMeshes.push(point);
          
            point.position.set(p.x, p.y, p.z);
            //point.rotation.x = -Math.PI * 0.5;
            
            scene.add(point);
            draggableList.push(point);
        }
    },

    // evaluate cubic bezize at t in [0.0, 1.0]
    eval(t) {
        const p1 = this.ctrlPoint(0).multiplyScalar((1 - t) ** 3);
        const p2 = this.ctrlPoint(1).multiplyScalar(3 * ((1 - t) ** 2) * t);
        const p3 = this.ctrlPoint(2).multiplyScalar(3 * (1 - t) * (t ** 2));
        const p4 = this.ctrlPoint(3).multiplyScalar(t ** 3);

        const x = p1.x + p2.x + p3.x + p4.x;
        const y = p1.y + p2.y + p3.y + p4.y;
        const z = p1.z + p2.z + p3.z + p4.z;

        return new THREE.Vector3(x, y, z);
    },

    sampleCurve(curvePoints, segments) {
        for (let t = 0; t < segments; t += 1) {
            //console.log( t / (segments-1.0) );
          
            let p = this.eval( t / (segments-1) );
            curvePoints.push( p );
        }
    },

    resolution: 25,
    curveLine: null,
  
    addCurve(scene){
        let curveGeometry = new THREE.BufferGeometry();  
        const positions = new Float32Array( g_bezier.resolution * 3 ); // 3 vertices per point
        curveGeometry.setAttribute( 'position', new THREE.BufferAttribute( positions, 3 ) );
        const curveMaterial = new THREE.LineBasicMaterial({
            color: 0x00ff02
        });

        this.curveLine = new THREE.Line(curveGeometry, curveMaterial);
        this.curveLine.position.set(0, 0, 0);
        //g_bezier.curveLine.rotation.x = Math.PI * 0.5;
        this.curveLine.geometry.attributes.position.needsUpdate = true;
        scene.add(g_bezier.curveLine);
    }
}; 
// ------------------------------------------------------------------------------


function init() {
    let container = document.createElement('div');
    document.body.appendChild(container);

    if (!g_orthographic) {
        g_camera = new THREE.PerspectiveCamera(
                60,
                g_cam_params.aspect_ratio(),
                g_cam_params.near, g_cam_params.far);
    } else {
        g_camera = new THREE.OrthographicCamera(
                g_cam_params.left(),
                g_cam_params.right(),
                g_cam_params.top(),
                g_cam_params.bottom(),
                g_cam_params.near,
                g_cam_params.far);
    }
    // OrthographicCamera( left , right , top , bottom , near , far )

    g_camera.position.set(0, 0, g_frustumSize);
    g_camera.lookAt(0, 0, 0);

    g_scene = new THREE.Scene();
    g_scene.background = new THREE.Color(0xffffff); //0xa0a0a0

    if (true) {
        let light1 = new THREE.HemisphereLight(0xffffff, 0x444444);
        light1.position.set(0, 200, 0);
        g_scene.add(light1);

        let light2 = new THREE.DirectionalLight(0xbbbbbb);
        light2.position.set(0, 200, 100);
        light2.castShadow = true;
        light2.shadow.camera.top = 180;
        light2.shadow.camera.bottom =  - 100;
        light2.shadow.camera.left =  - 120;
        light2.shadow.camera.right = 120;
        g_scene.add(light2);
    }
    //
    //scene.add(new THREE.CameraHelper(light.shadow.camera));

    // ground
    /*
    let mesh = new THREE.Mesh(new THREE.PlaneBufferGeometry(2000, 2000), new THREE.MeshPhongMaterial({ color: 0x999999, depthWrite: false }));
    mesh.rotation.x = - Math.PI / 2;
    mesh.receiveShadow = true;
    g_scene.add(mesh);
     */

    let grid = new THREE.GridHelper(
            g_frustumSize /*size*/,
            4 /*divisions*/,
            0x000000, 0x999999);
    //grid.divisions = 2;
    grid.position.set(0,0,-1);
    grid.rotation.x = -Math.PI * 0.5;
    grid.material.opacity = 1.0;
    grid.material.transparent = true;
    g_scene.add(grid);

    g_renderer = new THREE.WebGLRenderer({
        antialias: true
    });
    g_renderer.setPixelRatio(window.devicePixelRatio);
    g_renderer.setSize(window.innerWidth, window.innerHeight);
    g_renderer.shadowMap.enabled = true;
    document.body.style.margin = 0;
    document.body.style.padding = 0;
    document.body.style.overflow = 'hidden';
    document.body.style.position = 'fixed';
    container.appendChild(g_renderer.domElement);
    window.addEventListener('resize', onWindowResize, false);

    g_controls = new /*THREE.*/ OrbitControls(g_camera, g_renderer.domElement);
    g_controls.enableRotate = !g_orthographic;
    g_controls.screenSpacePanning = true;
    g_controls.target.set(0, 0, 0);
    g_controls.update();

    var draggableObjects = [];

    //draggableObjects.push(mesh);
    //g_scene.add(g_cogPointMesh);


    // Create a group to hold the points
    //const pointsGroup = new THREE.Group();
    //g_scene.add(pointsGroup);

    // Add Bézier curve and control points:
    g_bezier.addCtrlPoints(g_scene, draggableObjects);    
    g_bezier.addCurve(g_scene);
  
  

    /* Add blue box
    let boxGeometry = new THREE.BoxBufferGeometry(10, 10, 10);
    let blue = new THREE.MeshPhongMaterial({
    color: 0x3399dd
    });
    let white = new THREE.MeshLambertMaterial({
    color: 0x888888
    });

    let box = new THREE.Mesh(boxGeometry, blue);
    g_scene.add(box);
     */

    var dragControls = new /*THREE.*/ DragControls(draggableObjects, g_camera, g_renderer.domElement);

    dragControls.addEventListener('dragstart', function () {
        g_controls.enabled = false;
    });
    dragControls.addEventListener('dragend', function () {
        g_controls.enabled = true;
    });

}

// ---------------------------------------------------------------

function onWindowResize() {
    g_camera.left = g_cam_params.left();
    g_camera.right = g_cam_params.right();
    g_camera.top = g_cam_params.top();
    g_camera.bottom = g_cam_params.bottom();
    g_camera.aspect = g_cam_params.aspect_ratio();
    g_camera.updateProjectionMatrix();
    g_renderer.setSize(window.innerWidth, window.innerHeight);
}

// ------------------------------------------------------------------

function animate() {

    
    const curvePoints = [];
    g_bezier.sampleCurve(curvePoints, g_bezier.resolution); 
  
    for (let v = 0; v < g_bezier.resolution; v++) {
        let p = curvePoints[v];
        g_bezier.curveLine.geometry.attributes.position.setXYZ(v, p.x, p.y, p.z);
    }
    g_bezier.curveLine.geometry.attributes.position.needsUpdate = true; 
  
    //g_bezier.curveLine.geometry.computeBoundingBox();
    //g_bezier.curveLine.geometry.computeBoundingSphere();
 
  
    requestAnimationFrame(animate);
    g_renderer.render(g_scene, g_camera);

}

// -------------------------------------------------------------------


init();
animate();

First derivative

First, I'll give the result. Let $\vec p(u) : \mathbb R \rightarrow \mathbb R^3$ be our parametric function (a.k.a vectored valued function) of a Bézier curve:

$$ \vec p\phantom{'}(u) = \phantom{n} \sum_{i=0}^{n} B_i^{n}(u) \quad \vec p_i \\ $$

Where $\vec p_i$ are the control points of the curve, and $B_i^n(u) : \mathbb R \rightarrow \mathbb R$ are the Bernstein's polynoms of degree $n$ with a total of $n+1$ polynoms where $i \in [0, n]$ represent its index :

$$ B_i^n(u) = {n \choose i} u^i \ (1-u)^{n-i} \quad \text{ith Bernstein polynomial associated to } \vec p_i $$

Then the first derivative is expressed has:

$$ \vec v(u) = \vec p'(u) = n \ \sum_{i=0}^{n-1} B_i^{n-1}(u) \ \Delta \vec p_i \quad \text{ with } \Delta \vec p_i = \vec p_{i+1} - \vec p_i $$

Derivation

We take the derivative according to $u$ of $\vec p(u)$

$$ \begin{aligned} \vec p'(u) &= \left [ \sum_{i=0}^{n} B_i^{n}(u) \quad \vec p_i \right ]' \\ &= \sum_{i=0}^{n} \Big [ B_i^{n}(u)\Big ]' \quad \vec p_i \annotation{ \text{we can differentiate each term separatly } \\ \text{also control points are constant }} \\ &= \sum_{i=0}^{n} \left [ {n \choose i} u^i \ (1-u)^{n-i} \right ]' \quad \vec p_i \\ &= \sum_{i=0}^{n} {n \choose i} \Big [ u^i \ (1-u)^{n-i} \Big ]' \vec p_i \annotation{ \text{binomial terms are constant }} \end{aligned} $$

Let's develop $\Big [ u^i \ (1-u)^{n-i} \Big ]'$:

$$ \begin{aligned} &= \left [ u^i \ (1-u)^{n-i} \right ]' \\ &= [u^i]' \ (1-u)^{n-i} + u^i \ \left [ (1-u)^{n-i} \right ]' \annotation{\text{ product rule: }(uv)' = u'v + v'u} \\ &= \red{i.u^{i-1}} \ (1-u)^{n-i} + u^i \ \left [ (1-u)^{n-i} \right ]' \annotation{\text{ exponent rule: } x^n = n.x^{n-1} } \\ &= i.u^{i-1} \ (1-u)^{n-i} + u^i \ . \red{(-1).(n-i).(1-u)^{n-i-1} } \annotation{\text{ chain rule: } \\ \left[ f(g(x)) \right ]' = g'(x) . f'(g(x))} \\ &= i.u^{i-1} \ (1-u)^{n-i} \red{-} u^i \ (n-i).(1-u)^{n-i-1} \end{aligned} $$

Next we will show $ \left[ B^n_i(u)\right ]' = n \left ( B^{n-1}_{i-1}(u) - B^{n-1}_{i}(u) \right )$, to this end, let's consider back the binomial term and develop the expression $ \left[ B^n_i(u)\right ]' = {n \choose i} \left [ u^i \ (1-u)^{n-i} \right ]' $:

$$ \begin{aligned} & = {n \choose i} \Big ( i.u^{i-1} \ (1-u)^{n-i} - u^i \ (n-i).(1-u)^{n-i-1} \Big) \\ & = \red{ {n \choose i} } \ i.u^{i-1} \ (1-u)^{n-i} - \red{ {n \choose i} } \ u^i \ (n-i).(1-u)^{n-i-1} \\ & = \red{ \frac{n!}{i!(n-i)!} } \ i.u^{i-1} \ (1-u)^{n-i} - \red{ \frac{n!}{i!(n-i)!} } \ u^i \ (n-i).(1-u)^{n-i-1} \\ & = \frac{n!}{i!(n-i)!} \ i.u^{i-1} \ (1-u)^{n-i} -\red{ n } . \frac{ \red{ (n-1)}!}{i!(n-i)!} \ u^i \ (n-i).(1-u)^{n-i-1} \annotation{n! = n.(n-1)!} \\ & = \frac{n!}{i!(n-i)!} \ i.u^{i-1} \ (1-u)^{n-i} - n . \frac{ (n-1)!}{i!( \red{(n-1)}-i)!} \ u^i \ .(1-u)^{n-i-1} \annotation{\frac{n-i}{(n-i)!} =\frac{n-i}{(n-i) (n-i-1)!} } \\ &= \frac{n!}{i!(n-i)!} \ i.u^{i-1} \ (1-u)^{n-i} - n . \frac{ (n-1)!}{i!( (n-1)-i)!} \ u^i \ .(1-u)^{\red{(n-1)} - i} \\ &= \frac{n!}{i!(n-i)!} \ i.u^{i-1} \ (1-u)^{n-i} - n . \red{ B^{n-1}_i } \\ &= \red{ n } . \frac{ \red{ (n-1)}!}{i!(n-i)!} \ i.u^{i-1} \ (1-u)^{n-i} - n . B^{n-1}_i \\ &= n . \frac{(n-1)!}{\red{(i-1)}!(n-i)!} \ u^{i-1} \ (1-u)^{n-i} - n . B^{n-1}_i \annotation{\frac{i}{i!} =\frac{i}{(i)(i-1)!} } \\ &= n . \frac{(n-1)!}{(i-1)!(\red{(n-1)-(i-1)})!} \ u^{i-1} \ (1-u)^{n-i} - n . B^{n-1}_i \annotation{ \phantom{=} (n-1)-(i-1) \\= n-1-i+1 \\= n-i } \\ &= n . \frac{(n-1)!}{(i-1)!((n-1)-(i-1))!} \ u^{i-1} \ (1-u)^{\red{(n-1)-(i-1)}} - n . B^{n-1}_i \annotation{ \phantom{=} (n-1)-(i-1) \\= n-1-i+1 \\= n-i } \\ ~ \\ &= n . B^{n-1}_{i-1} - n . B^{n-1}_i \end{aligned} $$

So far we have

$$ \begin{aligned} \vec p'(u) &= n .\sum_{i=0}^{n} \Big( B^{n-1}_{i-1} - B^{n-1}_i \Big ) \vec p_i \end{aligned} $$

When developing the above sum we notice that the first and a last terms $\purple{ B^{n-1}_{-1} }$ and $ \green{ B^{n-1}_{n}}$ are null:

$$ \begin{array}{lllll} \vec p'(u) =& \Big( \purple{ B^{n-1}_{0-1} } & - & B^{n-1}_0 \Big ) & n.\vec p_0 & + & \annotation{ i = 0 }\\ & & \cdots & & & & \\ & \Big( B^{n-1}_{i-1} & - & B^{n-1}_{i} \Big ) & n. \vec p_{i} & + & \annotation{ i }\\ & \Big( B^{n-1}_{(i+1)-1} & - & B^{n-1}_{i+1} \Big ) & n. \vec p_{i+1} & + & \annotation{ i + 1 } \\ & & \cdots & & & \\ & \Big( B^{n-1}_{n-1} & - & \green{ B^{n-1}_{n} } \Big ) & n. \vec p_{n} & & \annotation{ i = n } \end{array} $$
$$ \purple{B^{n-1}_{0-1}} = B^{n-1}_{-1} = {n \choose i} \ \underbrace{i}_{0}.u^{i-1} \ (1-u)^{n-i} = 0 $$ $$ \green{ B^{n-1}_{n} } = {n \choose i} \ u^i \ \underbrace{(n-i)}_{0}.(1-u)^{n-i-1} = 0 $$

Factoring back with the remaining factors, the sum will give us:

$$ \vec v(u) = \vec p'(u) = n \ \sum_{i=0}^{n-1} B_i^{n-1}(u) \ \Delta \vec p_i \quad \text{ with } \Delta \vec p_i = \vec p_{i+1} - \vec p_i $$

2nd derivative

Because the first derivative gives us back the formula of a curve we can re-apply the formula recursively to it itself:

$$ \vec a(u) = \vec p''(u) = n \ (n-1) \ \sum_{i=0}^{n-2} B_i^{n-2}(u) \ \Delta \vec v_i \quad \text{ with } \Delta \vec v_i = \Delta \vec p_{i+1} - \Delta \vec p_i $$

r-th derivative

The r-th derivative of a Bézier curve is:

$$ \vec p^{(r)}(u) = \frac{n!}{(n − r)!} \sum_{i=0}^{n−r} \Delta^r \vec p_i B_i^{n−r}(u) $$

With:

$$ \Delta^1 \vec p_i = \vec p_{i+1} - \vec p_i \\ \Delta^r \vec p_i = \Delta^{(r-1)} \vec p_{i+1} - \Delta^{(r-1)} \vec p_i $$

Additional resources

- https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/Bezier/bezier-der.html
- A primer on Bézier

No comments

Rodolphe Vaillant's homepage

Research, teaching and more...