# Analytical expression

Let $$s: \mathbb R \rightarrow \mathbb R^3$$ be a vector valued function representing a parametric curve $$s(t)$$ with $$t \in \mathbb R$$ the curve's parameter; the curvature of $$s$$ is given by $$\kappa: \mathbb R \rightarrow \mathbb R$$ as follows:

$\kappa(t) = \frac{ \| s'(t) \times s''(t) \| }{ \| s'(t) \|^3 }$

Where '×' is the cross product, $$s'(t)$$ speed (velocity) at $$t$$, $$s''(t)$$ the acceleration, and $$\| \ \|$$ the Euclidean norm of a vector.

# 2D version

The 2D cross product is not a vector like the 3D cross product but a real value $$\mathbb R$$. In addition, the norm of the 2D cross product is the value of the cross product itself: $$\| s' \times s'' \| = s_x' s''_y - s_y's_x''$$. Finally the 2D cross product can be expressed as the determinant $$\left | \phantom{x} \right |$$ of the 2D matrix below:

$\kappa(t) = \frac{ \left | \begin{matrix} s_x'(t) & s_x''(t) \\ s_y'(t) & s_y''(t) \\ \end{matrix} \right | }{ \| s'(t) \|^3 }$

# Numerical computation

If you don't have the formula (i.e. analytical expression) of $$s(t)$$, you can numerically compute it using finite differences:

$$s'(t) = \frac{ s(t+h)-s(t-h) } {2h}$$ $$s''(t) = \frac{ s'(t+h)-s'(t-h) } {2h}$$

with $$h$$ "small" (ex $$h < 0.0001$$ )

Glsl code to visualize the curvature:

// Under MIT License

// Computes the curvature of a parametric curve f(x) as
// c(f) = | f' x f''| / |f'|^3

vec3 a, b, c, m, n;

// parametric curve value s(t):
vec3 mapD0(float t){
return 0.25 + a*cos(t+m)*(b+c*cos(t*7.0+n));
}

// curve derivative (velocity) s'(t)
vec3 mapD1(float t){
return -7.0*a*c*cos(t+m)*sin(7.0*t+n) - a*sin(t+m)*(b+c*cos(7.0*t+n));
}

// curve second derivative (acceleration) s''(t)
vec3 mapD2(float t){
return 14.0*a*c*sin(t+m)*sin(7.0*t+n) - a*cos(t+m)*(b+c*cos(7.0*t+n)) - 49.0*a*c*cos(t+m)*cos(7.0*t+n);
}

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

float curvature( float t ){
vec3 r1 = mapD1(t); // first derivative
vec3 r2 = mapD2(t); // second derivative
return length(cross(r1,r2)) / pow(length(r1),3.0);
}

float curvature_reciprocal( float t ){
vec3 r1 = mapD1(t); // first derivative
vec3 r2 = mapD2(t); // second derivative
return pow(length(r1),3.0) / length(cross(r1,r2));
}