8
$\begingroup$

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)$?

$\endgroup$
4
  • $\begingroup$ 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})$. $\endgroup$
    – dbx
    Nov 6 '17 at 14:35
  • $\begingroup$ 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. $\endgroup$ Nov 6 '17 at 17:31
  • $\begingroup$ 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. $\endgroup$ Nov 6 '17 at 17:32
  • $\begingroup$ 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? $\endgroup$ Nov 13 '17 at 8:54
8
+150
$\begingroup$

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()
$\endgroup$
8
  • $\begingroup$ I actually have more than just 3 points. Please elaborate on how this can be used to make a more 'precise' circle. $\endgroup$
    – cgo
    Nov 9 '17 at 7:19
  • $\begingroup$ I added an example how you can do this. Let me know if it is unclear. $\endgroup$
    – EdG
    Nov 9 '17 at 8:08
  • $\begingroup$ Where is the example? $\endgroup$
    – Narasimham
    Nov 10 '17 at 18:21
  • $\begingroup$ 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. $\endgroup$
    – EdG
    Nov 11 '17 at 2:04
  • 3
    $\begingroup$ Beware of blindly applying Heron's formula. Some care in writing the formula is necessary (see this as well). $\endgroup$ Nov 13 '17 at 7:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.