Numerical Approximation for 2D Curvature

8.1k Views Asked by At

I have a list of points (x, y) that are taken from an unknown 2D parametric curve $\vec{f}(t)$. These points are monotonically increasing in t (ie: they're a "connect the dots" version of $\vec{f}(t)$). I would like to numerically approximate the curvature at all of these points.

I know that

$\kappa = |\frac{dT}{ds}| = \frac{|\vec{T}'(t)|}{|\vec{r}'(t)|} = \frac{|\vec{r}'(t) \times \vec{r}'(t)|}{|\vec{r}'(t)|^3}$

and I was planning on using the last formula, approximating the derivative as:

$\frac{dx}{dt} \approx \frac{x(t + \Delta t) - x(t - \Delta t)}{2\Delta t}$

Is there a better way of doing this? (Maybe by fitting a curve to the points?)

4

There are 4 best solutions below

5
On

The problem is that although $t$ is monotone increasing, you don't know that the sequential $t$ values are in equal steps.

So the best thing to do is to try to fit a quadric to your data, and then compute the curvature of that quadric.

Details: Suppose the points are $x_i, y_i$. You want to take a few points to either side of $i$...say $x_{i-2}, y_{i-2}$ up through $x_{i+2}, y_{i+2}$. I'm going to shift indices and call these $(x_1, y_1)$ up to $(x_5, y_5)$.

The hope is that there's a curve, implicitly defined by something like $$ A x^2 + 2B xy + Cy^2 + 2Dx + 2E y + F = 0 $$ that all these points satisfy; if there is, then you can compute the curvature of the quadric. But how do you find $A... F$? You write a matrix equation: for each $x_i, y_i$ pair, say $x_1, y_1$, you have $$ x_1^2 A + 2x_1 y_1 B + y_1^2 C + 2x_1 D + 2y_1 E + 1 F = 0 $$ which you can write as $$ \begin{bmatrix} x_1^2 & 2x_1 y_1 & y_1^2 & 2x_1& 2y_1 &1 \end{bmatrix} \begin{bmatrix} A\\B\\C\\D\\E\\F\end{bmatrix} = 0. $$ If you do this for each $xy$ pair, perhaps for seven different points, you get seven row vectors on the left and seven zeroes on the right, and a matrix equation that looks like $$ M v = 0 $$ where $M$ is $7 \times 6$ and $v$ is the vector of unknown coefficients. You find a least-squares solution to this, and you've got the coefficients you need. (One problem: all zeros is a solution...so you have to restrict to the "nonzero" case, perhaps by insisting that the $F$ coefficient be 1 or something).

When you do that, and look for the curvature of the quadric, I believe that something like $B^2 - AC$ ends up being the curvature, but it's been a while since I worked it out.  

2
On

A very general method is to reconstruct (for example using splines) the parametric curve from the discrete data, and then compute the curvature of the reconstructed curve. The end answer can be understood as some finite difference operator on the discrete points. Different choice of interpolation will give different finite difference operators.

For example, suppose you are using a second order approximation, with the three points $ (x_i, y_i) $ for $ i = -1, 0, 1 $. A parabolic curve passing through the three points is: $$ x(t) = \frac{1}{2} (x_1 - 2 x_{0} + x_{-1}) t^2 + \frac{1}{2} (x_1 - x_{-1}) t + x_0 \\ y(t) = \frac{1}{2} (y_1 - 2 y_{0} + y_{-1}) t^2 + \frac{1}{2} (y_1 - y_{-1}) t + y_0 $$ You can already recognize the coefficients of $ t $ as finite second and first derivative operators: $$ D_t^2 x_i = \frac{1}{2} (x_{i-1} - 2 x_{i} + x_{i+1}) \\ D_t x_i = \frac{1}{2} (x_{i+1} - x_{i-1}) $$ Now you can find the curvature of the curve, for example using the formula: $$ \kappa = \frac{| x' y'' - y' x''|}{(x' x' + y' y')^{3/2}} $$ In fact, you should find: $$ \kappa(0) = \frac{| D_t x_0 \; D_t^2 y_0 - D_t y_0 \; D_t^2 x_0 |}{((D_t x_0)^2 + (D_t y_0)^2)^{3/2}} $$

2
On

A method for computing the curvature is described for example in:

Driscoll MK, McCann C, Kopace R, Homan T, Fourkas JT, et al. (2012) "Cell Shape Dynamics: From Waves to Migration." PLoS Comput Biol 8(3): e1002392. doi:10.1371/journal.pcbi.1002392

They have 400 points (they call them boundary points) and:

At each boundary point, we calculate the boundary curvature by fitting a circle to that boundary point and the two points that are 10 boundary points away from it. The magnitude of the boundary curvature is then defined as the reciprocal of the radius of that circle.

See also this video from the above paper.

0
On

Have a look at chapter "10.4 The Curvature of a Planar Digital Curve" in KLETTE, Reinhard; ROSENFELD, Azriel. Digital geometry: Geometric methods for digital picture analysis. Morgan Kaufmann, 2004.

See also Adina CIOMAGA, Pascal MONASSE, and Jean-Michel MOREL, The Image Curvature Microscope: Accurate Curvature Computation at Subpixel Resolution, Image Processing On Line, 7 (2017), pp. 197–217. https://doi.org/10.5201/ipol.2017.212