Introduction
Computing curvature from noisy data points requires smoothing the data first, then computing derivatives of the smooth curve. Direct numerical differentiation amplifies noise, making results unreliable. The standard approach is to fit a smoothing spline to the data points and then compute curvature analytically from the spline's derivatives. In C++, libraries like Eigen, Boost, and ALGLIB provide spline fitting capabilities.
2D Curvature
For a parametric curve (x(t), y(t)):
κ = |x'y'' - y'x''| / (x'² + y'²)^(3/2)
For y = f(x):
κ = |f''| / (1 + f'²)^(3/2)
3D Curvature
For a parametric curve (x(t), y(t), z(t)):
where r' × r'' is the cross product of the first and second derivatives.
Step 1: Smoothing Spline Fitting
Using Eigen (Cubic Spline)
1#include <Eigen/Dense>
2#include <vector>
3#include <cmath>
4
5class CubicSpline {
6 std::vector<double> x_, a_, b_, c_, d_;
7
8public:
9 void fit(const std::vector<double>& x, const std::vector<double>& y) {
10 int n = x.size() - 1;
11 x_ = x;
12 a_ = y;
13 b_.resize(n);
14 c_.resize(n + 1);
15 d_.resize(n);
16
17 std::vector<double> h(n), alpha(n);
18 for (int i = 0; i < n; i++)
19 h[i] = x[i+1] - x[i];
20
21 for (int i = 1; i < n; i++)
22 alpha[i] = 3.0/h[i] * (a_[i+1] - a_[i]) - 3.0/h[i-1] * (a_[i] - a_[i-1]);
23
24 // Solve tridiagonal system for c
25 std::vector<double> l(n+1), mu(n), z(n+1);
26 l[0] = 1; mu[0] = 0; z[0] = 0;
27
28 for (int i = 1; i < n; i++) {
29 l[i] = 2*(x[i+1] - x[i-1]) - h[i-1]*mu[i-1];
30 mu[i] = h[i] / l[i];
31 z[i] = (alpha[i] - h[i-1]*z[i-1]) / l[i];
32 }
33
34 l[n] = 1; z[n] = 0; c_[n] = 0;
35
36 for (int j = n-1; j >= 0; j--) {
37 c_[j] = z[j] - mu[j]*c_[j+1];
38 b_[j] = (a_[j+1] - a_[j])/h[j] - h[j]*(c_[j+1] + 2*c_[j])/3.0;
39 d_[j] = (c_[j+1] - c_[j]) / (3.0*h[j]);
40 }
41 }
42
43 int findSegment(double x) const {
44 int lo = 0, hi = x_.size() - 2;
45 while (lo < hi) {
46 int mid = (lo + hi) / 2;
47 if (x_[mid+1] < x) lo = mid + 1;
48 else hi = mid;
49 }
50 return lo;
51 }
52
53 double eval(double x) const {
54 int i = findSegment(x);
55 double dx = x - x_[i];
56 return a_[i] + b_[i]*dx + c_[i]*dx*dx + d_[i]*dx*dx*dx;
57 }
58
59 double derivative1(double x) const {
60 int i = findSegment(x);
61 double dx = x - x_[i];
62 return b_[i] + 2*c_[i]*dx + 3*d_[i]*dx*dx;
63 }
64
65 double derivative2(double x) const {
66 int i = findSegment(x);
67 double dx = x - x_[i];
68 return 2*c_[i] + 6*d_[i]*dx;
69 }
70};
Step 2: Computing Curvature
2D Curvature from y = f(x)
1double curvature2D(const CubicSpline& spline, double x) {
2 double fp = spline.derivative1(x); // f'
3 double fpp = spline.derivative2(x); // f''
4 return std::abs(fpp) / std::pow(1.0 + fp*fp, 1.5);
5}
6
7// Usage
8CubicSpline spline;
9spline.fit(x_data, y_data);
10
11for (double x = x_data.front(); x <= x_data.back(); x += 0.01) {
12 double k = curvature2D(spline, x);
13 std::cout << x << " " << k << "\n";
14}
2D Parametric Curvature
For curves that cannot be expressed as y=f(x) (e.g., loops):
1struct Point2D { double x, y; };
2
3double parametricCurvature2D(
4 const CubicSpline& sx, const CubicSpline& sy, double t)
5{
6 double xp = sx.derivative1(t);
7 double yp = sy.derivative1(t);
8 double xpp = sx.derivative2(t);
9 double ypp = sy.derivative2(t);
10
11 double num = std::abs(xp * ypp - yp * xpp);
12 double den = std::pow(xp*xp + yp*yp, 1.5);
13
14 return (den > 1e-12) ? num / den : 0.0;
15}
16
17// Fit separate splines for x(t) and y(t)
18std::vector<double> t_param(points.size());
19std::iota(t_param.begin(), t_param.end(), 0); // t = 0, 1, 2, ...
20
21CubicSpline sx, sy;
22sx.fit(t_param, x_coords);
23sy.fit(t_param, y_coords);
3D Curvature
1struct Vec3 { double x, y, z; };
2
3Vec3 cross(const Vec3& a, const Vec3& b) {
4 return {a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x};
5}
6
7double magnitude(const Vec3& v) {
8 return std::sqrt(v.x*v.x + v.y*v.y + v.z*v.z);
9}
10
11double curvature3D(
12 const CubicSpline& sx, const CubicSpline& sy, const CubicSpline& sz,
13 double t)
14{
15 Vec3 r1 = {sx.derivative1(t), sy.derivative1(t), sz.derivative1(t)};
16 Vec3 r2 = {sx.derivative2(t), sy.derivative2(t), sz.derivative2(t)};
17
18 Vec3 cp = cross(r1, r2);
19 double num = magnitude(cp);
20 double den = std::pow(magnitude(r1), 3);
21
22 return (den > 1e-12) ? num / den : 0.0;
23}
Handling Noise: Smoothing Splines
A smoothing spline balances fitting the data closely and maintaining smoothness. The smoothing parameter controls this tradeoff:
1// Using ALGLIB for smoothing splines
2#include "interpolation.h"
3
4void computeCurvatureSmoothed(
5 const std::vector<double>& x,
6 const std::vector<double>& y,
7 double smoothing_param)
8{
9 alglib::real_1d_array xa, ya;
10 xa.setcontent(x.size(), x.data());
11 ya.setcontent(y.size(), y.data());
12
13 alglib::spline1dinterpolant s;
14 alglib::spline1dfitpenalized(xa, ya, x.size(), smoothing_param, s);
15
16 // Evaluate curvature at each point
17 for (size_t i = 0; i < x.size(); i++) {
18 double val, d1, d2;
19 alglib::spline1ddiff(s, x[i], val, d1, d2);
20
21 double curvature = std::abs(d2) / std::pow(1.0 + d1*d1, 1.5);
22 std::cout << x[i] << " " << curvature << "\n";
23 }
24}
Smoothing Parameter Selection
| Parameter | Effect |
| Too low smoothing | Spline overfits noise, curvature is erratic |
| Too high smoothing | Spline over-smooths, real curvature features are lost |
| Optimal | Captures true curvature while suppressing noise |
Methods to choose the smoothing parameter:
Cross-validation: Leave-one-out or k-fold
GCV (Generalized Cross-Validation): Automated selection
Visual inspection: Plot the spline against data points
Common Pitfalls
Direct differentiation of noisy data: Numerical differentiation amplifies noise. First derivatives double the noise level; second derivatives (needed for curvature) quadruple it. Always smooth first.
Parameterization choice: For parametric curves, using arc-length parameterization gives more uniform curvature estimates than using point index as the parameter. Compute cumulative chord length for better results.
Boundary effects: Splines have reduced accuracy at endpoints due to fewer neighboring points. Consider using natural boundary conditions or discarding curvature estimates at the first and last few points.
Non-uniform spacing: If data points are unevenly spaced, some spline implementations produce poor fits. Use libraries that explicitly handle non-uniform knot spacing.
Smoothing parameter sensitivity: Curvature is very sensitive to the smoothing parameter. Too little smoothing gives noisy curvature; too much removes real curvature features. Always validate with cross-validation or domain knowledge.
Summary
Fit a smoothing spline to noisy data before computing curvature — never differentiate noisy data directly
Use the curvature formula κ = |f''| / (1 + f'²)^(3/2) for 2D curves and the cross-product formula for 3D
For parametric curves, fit separate splines for each coordinate (x(t), y(t), z(t))
The smoothing parameter critically affects curvature estimates — tune with cross-validation
Use C++ libraries (ALGLIB, Eigen, Boost) for production-quality spline fitting
