# Numerical way to solve for the curvature of a curve

I am aware of the equation to solve for the curvature of a curve. It is $\kappa = \displaystyle\frac{|r' \times r''|}{|r'|^3}$, where $r$ is the parametrization of the curve.

In the absence of the parametric equation, where only data points are given, how can you solve/approximate the curvature with error about $O(h^2)$?

• curvature is also the inverse of the radius of the osculating circle; the curvature at $(x_k, y_k)$ can be approximated by inverting the radius of the unique circle through $(x_{k-1},y_{k-1})$, $(x_k, y_k)$, and $(x_{k+1},y_{k+1})$.
– dbx
Nov 6 '17 at 14:35
• Note that the convergence results about any notion of discrete curvature can be pretty subtle. For example, if $\gamma$ is a smooth plane curve that traces out the unit circle, one can easily construct a sequence of increasingly oscillatory discrete curves that converge pointwise to $\gamma$. Any notion of discrete curvature that I've seen does not converge to the underlying smooth distribution of curvature under just pointwise convergence. Nov 6 '17 at 17:31
• On the other hand, if you consider convergence under the Sobolev norm $d(\Gamma^i,\gamma)=\int\|\Gamma^i(s)-\gamma(s)\|+\|T_j(s)-\gamma'(s)\|~ds$, then you'll get convergence (I'm not sure of the rates exactly rn) of every discretization of curvature. The key is that you need both positions and tangent vectors to converge. In particular, if all of your discrete curves are sampled exactly from $\gamma$ then you should be okay :) I'll write an answer later tonight describing some of these different notions of discrete curvature, and I'll talk a bit more about convergence. Nov 6 '17 at 17:32
• What data points do you have? Do you have $r(t_j)$ for $t_j = a + jh$ for a fixed time step $h$ and an initial value $a$? Alternatively, are the points $r_j$ equidistant when measured along the curve? Moreover, is the curve smooth? If not, then what is the highest derivative known to exists? Nov 13 '17 at 8:54

As dbx mentioned in the comments, you can use the radius of a circle passing through three points. The curvature is the inverse of the radius. The radius can be found through $$r = \frac{abc}{4k},$$ where $$a$$, $$b$$ and $$c$$ denote the distances between the three points and $$k$$ denotes the area of the triangle formed by the three points. Obviously, the curvature is the reciprocal of this, thus $$\kappa = \frac{4k}{abc}$$ I happened to code this in the past in python. Below is the code. I used Heron's formula for computing the area of the triangle $$k$$.

# points_utm is a 3-by-2 array, containing the easting and northing coordinates of 3 points
# Compute distance to each point
a = np.hypot(points_utm[0, 0] - points_utm[1, 0], points_utm[0, 1] - points_utm[1, 1])
b = np.hypot(points_utm[1, 0] - points_utm[2, 0], points_utm[1, 1] - points_utm[2, 1])
c = np.hypot(points_utm[2, 0] - points_utm[0, 0], points_utm[2, 1] - points_utm[0, 1])

# Compute inverse radius of circle using surface of triangle (for which Heron's formula is used)
k = np.sqrt((a+(b+c))*(c-(a-b))*(c+(a-b))*(a+(b-c))) / 4    # Heron's formula for triangle's surface
den = a*b*c  # Denumerator; make sure there is no division by zero.
if den == 0.0:  # Very unlikely, but just to be sure
curvature = 0.0
else:
curvature = 4*k / den


Additional note: this might not be optimal as only three points are used to compute the curvature. If more points are available, it would be beneficial to use them. Still, one can use the same approach: fit a circle through the datapoints and compute the reciprocal of the radius. However, it is not trivial to find the best circle fit for your data. For more information, see http://scipy-cookbook.readthedocs.io/items/Least_Squares_Circle.html.

Based on above's link, I wrote a small python script that shows how you can compute the curvature using a linear least squares fit. Here's the code:

from scipy import optimize
import numpy as np
import matplotlib.pyplot as plt

class ComputeCurvature:
def __init__(self):
""" Initialize some variables """
self.xc = 0  # X-coordinate of circle center
self.yc = 0  # Y-coordinate of circle center
self.r = 0   # Radius of the circle
self.xx = np.array([])  # Data points
self.yy = np.array([])  # Data points

def calc_r(self, xc, yc):
""" calculate the distance of each 2D points from the center (xc, yc) """
return np.sqrt((self.xx-xc)**2 + (self.yy-yc)**2)

def f(self, c):
""" calculate the algebraic distance between the data points and the mean circle centered at c=(xc, yc) """
ri = self.calc_r(*c)
return ri - ri.mean()

def df(self, c):
""" Jacobian of f_2b
The axis corresponding to derivatives must be coherent with the col_deriv option of leastsq"""
xc, yc = c
df_dc = np.empty((len(c), self.xx.size))

ri = self.calc_r(xc, yc)
df_dc = (xc - self.xx)/ri                   # dR/dxc
df_dc = (yc - self.yy)/ri                   # dR/dyc
df_dc = df_dc - df_dc.mean(axis=1)[:, np.newaxis]
return df_dc

def fit(self, xx, yy):
self.xx = xx
self.yy = yy
center_estimate = np.r_[np.mean(xx), np.mean(yy)]
center = optimize.leastsq(self.f, center_estimate, Dfun=self.df, col_deriv=True)

self.xc, self.yc = center
ri = self.calc_r(*center)
self.r = ri.mean()

return 1 / self.r  # Return the curvature

# Apply code for an example
x = np.r_[36, 36, 19, 18, 33, 26]
y = np.r_[14, 10, 28, 31, 18, 26]
comp_curv = ComputeCurvature()
curvature = comp_curv.fit(x, y)

# Plot the result
theta_fit = np.linspace(-np.pi, np.pi, 180)
x_fit = comp_curv.xc + comp_curv.r*np.cos(theta_fit)
y_fit = comp_curv.yc + comp_curv.r*np.sin(theta_fit)
plt.plot(x_fit, y_fit, 'k--', label='fit', lw=2)
plt.plot(x, y, 'ro', label='data', ms=8, mec='b', mew=1)
plt.xlabel('x')
plt.ylabel('y')
plt.title('curvature = {:.3e}'.format(curvature))
plt.show()

• I actually have more than just 3 points. Please elaborate on how this can be used to make a more 'precise' circle.
– cgo
Nov 9 '17 at 7:19
• I added an example how you can do this. Let me know if it is unclear.
– EdG
Nov 9 '17 at 8:08
• Where is the example? Nov 10 '17 at 18:21
• It is the python code. If you scroll down a bit, you'll see "# Apply code for an example". If you run the python code, the resulting curvature is shown in a figure.
– EdG
Nov 11 '17 at 2:04
• Beware of blindly applying Heron's formula. Some care in writing the formula is necessary (see this as well). Nov 13 '17 at 7:19