Numerical way to solve for the curvature of a curve

Solution 1:

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[0] = (xc - self.xx)/ri                   # dR/dxc
        df_dc[1] = (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)[0]

        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()