Curve curvature

The curvature of a curve corresponds to the variation of its tangent : if the tangent stays the same, the curve is flat. When the curve bends, its tangent vector varies. Since the tangent is related to the derivative of the curve, the curvature is therefore related to its second derivative.

%matplotlib inline

#These are the python libraries used to illustrate this course

import matplotlib
import matplotlib.pyplot as plt
import matplotlib.animation
import numpy as np
import io
import math
from IPython.display import HTML
import ipywidgets as widgets
from mpl_toolkits.mplot3d import Axes3D

Let's first plot some Bézier curves and their tangent vector.

def bezier(p0, p1, p2, p3):

    B0 = lambda t : (1-t)**3
    B1 = lambda t : 3*t*((1-t)**2)
    B2 = lambda t : 3*(t**2)*(1-t)
    B3 = lambda t : t**3

    return lambda t : p0 * B0(t) + p1 * B1(t) + p2 * B2(t) + p3 * B3(t)

def bezier_tangent(p0, p1, p2, p3):
    B0t = lambda t : -3*((1-t)**2)
    B1t = lambda t : 3*(1-t)*(1-3*t)
    B2t = lambda t : 3*t*(2-3*t)
    B3t = lambda t : 3*(t**2)

    return lambda t : p0 * B0t(t) + p1 * B1t(t) + p2 * B2t(t) + p3 * B3t(t)

def bezier_tangent_update(cx = 0.2, t = 0.5):
    p = np.array([
        [0,0],
        [cx, 1],
        [0.8, -1],
        [1,0]
    ])

    bx = bezier(*p[:,0])
    by = bezier(*p[:,1])
    ts = np.linspace(0,1,100)

    fig, ax = plt.subplots(figsize = (8,8))
    
    ax.plot(bx(ts), by(ts))
    ax.plot(p[:2,0], p[:2,1], '-o')
    ax.plot(p[2:,0], p[2:,1], '-o')
    
    pt = np.array([bezier(*p[:,0])(t), bezier(*p[:,1])(t)])
    ax.scatter([pt[0]], [pt[1]], color = 'tomato')
    
    dpt = np.array([bezier_tangent(*p[:,0])(t), bezier_tangent(*p[:,1])(t)])
    ax.quiver([pt[0]], [pt[1]], [dpt[0]], [dpt[1]], scale = 20, width = 0.008, angles = 'xy', color = 'black', zorder = 2)
    
    plt.show()
    
_ = widgets.interact(bezier_tangent_update, cx = (0,1,0.1), t = (0,1,0.02))

The tangent vector varies in two ways: its length and its direction. The variation in the length of the tangent vector is due to the parameterization. In its current form, when the parameter $t$ varies linearly, the corresponding point on the curve does not move at constant speed along the curve. Using a different parameterization for the same curve will generate different variations in length. When the length of the tangent vector is constant, the curve is parameterized by arclength. The second derivative of a parameterisation by arclength yields the curvature. With other parameterizations, we have to take that into account, since a change in tangent length does not mean the curve bends. Formally, for a curve $\mathbf{c}(t)$, its tangent is given by

$$\mathbf{\dot{c}} = \frac{\mathrm{d}\mathbf{c}}{\mathrm{d}t}.$$

If the curve were parameterized by arclength $s$ , its tangent would be

$$\frac{\mathrm{d}\mathbf{c}}{\mathrm{d}s} = \frac{\mathbf{\dot{c}}}{\lVert\mathbf{\dot{c}}\rVert}.$$

Due to the chain tule, we have

$$\frac{\mathrm{d}\mathbf{c}}{\mathrm{d}s} = \frac{\mathrm{d}\mathbf{c}}{\mathrm{d}t}\frac{\mathrm{d}t}{\mathrm{d}s} = \mathbf{\dot{c}}\frac{\mathrm{d}t}{\mathrm{d}s}.$$

Identifying with the previous equation this yields

$$\frac{\mathrm{d}t}{\mathrm{d}s} = \frac{1}{{\lVert\mathbf{\dot{c}}\rVert}}.$$

The second derivative of the curve parameterized by artlength is therefore

$$ \begin{aligned} \frac{\mathrm{d}^2\mathbf{c}}{\mathrm{d}^2s} &= \frac{\mathrm{d}}{\mathrm{d}t} \left(\frac{\mathrm{d}\mathbf{c}}{\mathrm{d}s}\right) \frac{1}{\lVert\mathbf{\dot{c}}\rVert}\\ &= \frac{\mathrm{d}}{\mathrm{d}t} \left(\frac{\mathbf{\dot{c}}}{\lVert\mathbf{\dot{c}}\rVert}\right) \frac{1}{\lVert\mathbf{\dot{c}}\rVert}\\ &= \frac{1}{\lVert\mathbf{\dot{c}}\rVert} \frac{ \mathbf{\ddot{c}}\lVert\mathbf{\dot{c}} \rVert - \mathbf{\dot{c}}\frac{\mathrm{d}}{\mathrm{d}t}\lVert\mathbf{\dot{c}} \rVert }{\lVert\mathbf{\dot{c}}\rVert^2}\\ \end{aligned} $$

Using the fact that $\lVert \mathbf{u} \rVert = (\mathbf{u}.\mathbf{u})^{\frac{1}{2}}$ we have that

$$ \frac{\mathrm{d}}{\mathrm{d}t}\lVert \mathbf{u} \rVert = \frac{\mathbf{\dot{u}}.\mathbf{u}}{\lVert \mathbf{u} \rVert} $$

and therefore

$$ \begin{aligned} \frac{\mathrm{d}^2\mathbf{c}}{\mathrm{d}^2s} &= \frac{1}{\lVert\mathbf{\dot{c}}\rVert^2}\left( \mathbf{\ddot{c}} - \mathbf{\ddot{c}}.\mathbf{\dot{c}}\frac{\mathbf{\dot{c}}}{\lVert\mathbf{\dot{c}} \rVert^2} \right)\\ \end{aligned} $$

Intuitively, the $\frac{1}{\lVert\mathbf{\dot{c}}\rVert^2}$ factor corresponds to two compensations of the fact that we derive wrt. $t$ rather than $s$. The $\mathbf{\ddot{c}} - \mathbf{\ddot{c}}.\mathbf{\dot{c}}\frac{\mathbf{\dot{c}}}{\lVert\mathbf{\dot{c}} \rVert^2}$ corresponds to the orthogonal projection of $\mathbf{\ddot{c}}$ on the plane orthogonal to the tangent : we substract the component of $\mathbf{\ddot{c}}$ that is colinear with $\mathbf{\dot{c}}$.

The curvature $\kappa$ corresponds to the length of that projection. When it is zero, the tangent does not change direction. The radius of curvature is the radius of the osculating circle at $t$ and has value $\frac{1}{\kappa}$.

def bezier_second(p0, p1, p2, p3):
    B0tt = lambda t : 6*(1-t)
    B1tt = lambda t : 18*t - 12
    B2tt = lambda t : 6 - 18*t
    B3tt = lambda t : 6*t

    return lambda t : p0 * B0tt(t) + p1 * B1tt(t) + p2 * B2tt(t) + p3 * B3tt(t)

def bezier_curv_update(cx = 0.2, t = 0.3):
    p = np.array([
        [0,0],
        [cx, 1],
        [0.8, -1],
        [1,0]
    ])

    bx = bezier(*p[:,0])
    by = bezier(*p[:,1])
    ts = np.linspace(0,1,100)

    fig, ax = plt.subplots(figsize = (8,8))
    ax.set_xlim((-0.2, 1.2))
    ax.set_ylim((-0.7, 0.7))

    
    ax.plot(bx(ts), by(ts))
    ax.plot(p[:2,0], p[:2,1], '-o')
    ax.plot(p[2:,0], p[2:,1], '-o')
    
    pt = np.array([bezier(*p[:,0])(t), bezier(*p[:,1])(t)])
    ax.scatter([pt[0]], [pt[1]], color = 'tomato')
    
    dpt = np.array([bezier_tangent(*p[:,0])(t), bezier_tangent(*p[:,1])(t)])
    
    ax.quiver([pt[0]], [pt[1]], [dpt[0]], [dpt[1]], scale = 20, width = 0.008, angles = 'xy', color = 'black', zorder = 2)
    
    dptt = np.array([bezier_second(*p[:,0])(t), bezier_second(*p[:,1])(t)])

    n = dptt - dptt.dot(dpt) * dpt / (dpt.dot(dpt))
    ax.quiver([pt[0]], [pt[1]], [n[0]], [n[1]], width = 0.008, angles = 'xy', color = 'black', zorder = 2)

    curv = np.linalg.norm(n) / np.linalg.norm(dpt)**2
    if abs(curv) > 1e-5:
        radius = abs(1 /curv)
        center = pt + n * radius / np.linalg.norm(n)
        ax.scatter([center[0]], [center[1]], color = 'tomato')
        ax.add_patch(matplotlib.patches.Circle(center, radius, fill = None))
    
    plt.show()

def bezier_curv_update(cx = 0.2, t = 0.3):
    p = np.array([
        [0,0],
        [cx, 1],
        [0.8, -1],
        [1,0]
    ])

    bx = bezier(*p[:,0])
    by = bezier(*p[:,1])
    ts = np.linspace(0,1,100)

    fig, ax = plt.subplots(figsize = (8,8))
    ax.set_xlim((-0.2, 1.2))
    ax.set_ylim((-0.7, 0.7))

    
    ax.plot(bx(ts), by(ts))
    ax.plot(p[:2,0], p[:2,1], '-o')
    ax.plot(p[2:,0], p[2:,1], '-o')
    
    pt = np.array([bezier(*p[:,0])(t), bezier(*p[:,1])(t)])
    ax.scatter([pt[0]], [pt[1]], color = 'tomato')
    
    dpt = np.array([bezier_tangent(*p[:,0])(t), bezier_tangent(*p[:,1])(t)])
    
    ax.quiver([pt[0]], [pt[1]], [dpt[0]], [dpt[1]], scale = 20, width = 0.008, angles = 'xy', color = 'black', zorder = 2)
    
    dptt = np.array([bezier_second(*p[:,0])(t), bezier_second(*p[:,1])(t)])

    n = dptt - dptt.dot(dpt) * dpt / (dpt.dot(dpt))
    ax.quiver([pt[0]], [pt[1]], [n[0]], [n[1]], width = 0.008, angles = 'xy', color = 'black', zorder = 2)

    curv = np.linalg.norm(n) / np.linalg.norm(dpt)**2
    if abs(curv) > 1e-5:
        radius = abs(1 /curv)
        center = pt + n * radius / np.linalg.norm(n)
        ax.scatter([center[0]], [center[1]], color = 'tomato')
        ax.add_patch(matplotlib.patches.Circle(center, radius, fill = None))
    
    plt.show()
    
_ = widgets.interact(bezier_curv_update, cx = (0,1,0.1), t = (0,1,0.02))