# Curvature of a 2D image

50 views (last 30 days)

Show older comments

Francesco Pignatelli
on 9 Dec 2022

Edited: J. Alex Lee
on 28 Dec 2022

Hello all,

I would like to plot the Probability Density Function of the curvature values of a list of 2D image. Basically I would like to apply the following formula for the curvature:

k = (x'(s)y''(s) - x''(s)y'(s)) / (x'(s)^2 + y'(s)^2)^2/3

where x and y are the transversal and longitudinal coordinates, s is the arc length of my edge and ' and '' denote the first and second derivative.

What I have done so far, is to binarize and extracting the edge, as you can see here:

however, I am quite puzzled about how to get the curvaure based on the previous formula.

Does anyone has any hints about how to get the curvaure of such edge?

##### 4 Comments

J. Alex Lee
on 23 Dec 2022

Edited: J. Alex Lee
on 23 Dec 2022

@Francesco Pignatelli, can you also post the code you used to identify the edges you are interested in analyzing?

To your questions:

- I'm not sure if by "3rd order" you mean specifically that versus maybe mis-reading that i was considering spline differentiation as a 3rd method in a list of methods I could think of...anyway, there's a matlab command called "spline", or you could also call it through "interp1". I believe they all default to cubic (3rd order) pieces. You'd want at least that to ensure non-zero second derivatives. Aside: a quick exploration taught me that out of the 3 available piecewise interpolation options "pchip","makima","spline", only "spline" ensures continuous second derivative across the pieces. At this stage in solving your problem, I would not get too hung up on this (see second bullet below) - no matter what you do here, you will have a problem because of the fact that your starting set of coordinates come from a binarized image.
- To be clear, the problem isn't about the image, it's about once you convert an image (however ideal it is) into boundary coordinates, which can only have single-pixel resolution, assuming you trace boundaries from a binarized image. So, consider, in a perfect circle of constant curvature, you will inevitably have edges at N,S,E,W with long stretches of completely straight lines. It's hard to imagine succeeding to recover the correct radius from such set of coordinates. You would have to sample pixels far beyond the range of the straight regions to get any curvature at all (straight edges have zero curvature). This problem persists no matter what method you are using to compute curvature, circle fitting or derivatives.

### Accepted Answer

J. Alex Lee
on 22 Dec 2022

Edited: J. Alex Lee
on 22 Dec 2022

The calculation of arc length and curvature using arc length parameterization and its derivatives, given "perfect coordinates", is straightforward application of pythagoreas and basic trig. Given a curve defined by (x,y), define arc length parameterization

So

where s is arc length ϕ is inclination angle.

Implementation wise, if you choose any black box differentiator to evaluate the derivatives, then you can say

The formulat given by @Francesco Pignatelli in the original question (if it is correct - i have not seen that form before) - can be gotten from the above equations. But, easier to just compute the inclination angle ϕ since we need to be taking derivatives anyway.

There are 2 steps:

- define arc length coordinate for each xy coordinate using pythgoreas. This is implemented in my attached function "arclen"
- determine inclination and curvature using derivatives of parameters w.r.t arc length, implemented in my attached function "plancurv_arclen"

You can choose to take derivatives any way you want (finite difference on various stensiles, diff(p)./diff(s), etc. I chose convolutions for its convenience in dealing with closed loop (such as boundaries of images). But the wrinkle is to be "legit", you will need the xy coordinates sampled uniformly in arc length, which means the coordinates need to be resampled, and in a way that you don't know a priori will result in uniform sampling once you re-compute the arc length based on the newly sampled points. So you can just resample a few times and those iterations should converge - practically this seems to work and i dont have reason to believe it would ever fail, but i dont know how to prove that or anything.

To demonstrate on a "perfect" set of edge coordinates, run the below, which compares a straightforward implementation by me (plancurv_arclen) compared to modified version of @Image Analyst's circle fitting and the FEX contribution by Manjunatha (removing the preprocessing for dealing with images and/or freehand ROIs). I did not survey other FEX contributions or ANSWERS on this topic, but the Manjunatha method is also in the other thread between you two, so I figured to throw it in. I don't think it is fundamentally different from the circle fitting, though I did not dive into the details. I did play around with parameters for other methods to make them closer.

Some take-aways

- IA's circle fit method is blind to concavity
- Both IA and Driscoll-Manjunatha can only accomodate a single fitting window, but this seems inadequate because you would want narrower windows to better resolve regions of higher curvature and vice versa - this manifests in IA's method, but not as much in Driscoll-Manjunatha.
- The parallelization overhead for Driscoll-Manjunatha outweights benefits in this case - without parfor, my laptop ran it in no more than 100 times the runtime of the arc length method, which came out to 0.26s vs .003s.

Now, the application of interest does not involve perfect coordinates, but rather severly aliased coordinates coming from an image. I will make a comment to this answer to address.

Edit: it occurs to me that since all of us are using splines anyways, we should be able to take derivatives directly of the spline and make things even cleaner - i just don't have a lot of experience working with splines to implement.

demo_perfect_coordinates

##### 7 Comments

Image Analyst
on 28 Dec 2022

J. Alex Lee
on 28 Dec 2022

### More Answers (1)

Image Analyst
on 11 Dec 2022

See my answer here, where I compute curvature by fitting a circle to a sliding window:

##### 15 Comments

J. Alex Lee
on 23 Dec 2022

@Francesco Pignatelli, thanks for the link, Bruno has the missing link to provide that 3rd method based on direct spline differentiation I mention in below answer.

In any of the cases for computing curvature (either based on circle fits or arc length parameter differentiation, which in turn can be accomplished in many ways), the main problem for you and for others on this forum, deal with boundaries defined by pixel coordinates of images. So the zero-th order problem is actually not the curvature-seeking problem, but the image/coordinate pre-processing problem to avoid pixel aliasing to corrupt whatever measure you are seeking.

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!