Curvature of a 2D image
Show older comments
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 10 Dec 2022
Edited: J. Alex Lee
on 23 Dec 2022
First you need to order your edge location (xy) coordinates along the arc of the curve; this can be achieved with the image processing toolbox in-built "bwboundaries", so you can do that before you extract the edge (because the edge image is not going to be all that helpful to get your ordered coordinates along perimeter).
In principle, you can then calculate the arc coordinate for each xy using pythagoras on each distance between each consecutive pair: 
% generate a curve
x = transpose(linspace(0,2,7));
y = x.^2;
% pack into a coordinates array
xy = [x,y];
% compute [pseudo] arc-length
d = diff(xy); % nearest neighbor component diff
r = sqrt(sum(d.^2,2)); % nearest neighbor distance
s = cumsum([0;r]); % arc length
You can then think about several ways to achieve the derivatives: finite differences, splines, etc., to implement your formula above.
Here's a simple way with 1st order finite differences using diff(), but you'd have to modify to deal with edges (you keep shortening the coordinate lists by 1 each application of diff)
dxds = diff(xy(:,1))./diff(s); % this also equals cos(phi) where phi is inclination angle
dyds = diff(xy(:,2))./diff(s);
phi = unwrap(atan2(dyds,dxds);
curv = diff(phi)./diff(s(2:end)); % curvature is just the derivative of inclination angle
ddxdss = diff(dxds)./diff(s(2:end));
% etc
You've already linked to another answer involving spline interpolation and differentiating on the spline.
You could imagine other ways to take derivatives, like higher order finite differences (can do on unevenly sampled arc length) or convolutions (need evenly sampeld arc length).
In practice, the problem is that coordinates you get from images will be aliased so derivatives of your curve, parameterized or not, are going to be really bad - useless if you are only using a few neighboring pixels.
So your real problem is actually about how to pre-process your pixel coordinates and/or how to sample your coordinates to avoid aliasing problems.
Francesco Pignatelli
on 23 Dec 2022
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.
Francesco Pignatelli
on 24 Dec 2022
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
Francesco Pignatelli
on 23 Dec 2022
J. Alex Lee
on 23 Dec 2022
@Francesco Pignatelli, glad that you have already tried - what I wanted to say above was that my implementation if applied directly to bwboundaries coordinates of an image will almost certainly fail and give horrible results. There is a need to properly pre-process your image coordinates before applying any method to calculate curvature based on derivatives.
@Francesco Pignatelli, @Image Analyst, Here's a complete package with different methods to compute curvature from xy coordinates. There are several downsampling and upsampling options to aid with smoothing of pixel aliasing when coordinates are extracted from images.
The basic idea behind downsampling is to calculate a pseudo-curvature (arbitrary units) and use as weights to determine sampling concentrations: sample more where absolute curvature is high. Upsampling from this downsampled and smoothed description can help avoid some of the oscillations you get when upsampling from aliased coordinates.
I also modified Image Analyst's circle fit method to return a signed curvature. The sign is determined from the cross product of the tangent along the curve, and the normal directed toward the center of the radius of curvature at that point.
Demonstrated on the same lobed shape, which was converted to an image. All methods do a decent job if you pick the right set of parameters and don't scrutinize too much.
demo_birdie_b
edit: added histograms
And here is an application to the original image posted by @Francesco Pignatelli
demo_pignatelli
Francesco Pignatelli
on 28 Dec 2022
Image Analyst
on 28 Dec 2022
Just to finish up and answer the original question can you "plot the Probability Density Function of the curvature values"?
J. Alex Lee
on 28 Dec 2022
The histogram of curvature values will change depending on how values are sampled along the curve, but sure, I'll edit the previous comment with the histograms
More Answers (1)
Image Analyst
on 11 Dec 2022
0 votes
See my answer here, where I compute curvature by fitting a circle to a sliding window:
15 Comments
Francesco Pignatelli
on 12 Dec 2022
Image Analyst
on 12 Dec 2022
I did not use that formula. I simply fit the data in a sliding window to a circle and gave you the radius of the circle.
Francesco Pignatelli
on 12 Dec 2022
Image Analyst
on 12 Dec 2022
If the formula is y=ax^2 + bx + c, the derivative is 2ax+b.
But here the equation has two independent variables: x and y and I wasn't sure how to fit something like a polynomial when both x and y vary. I guess you could do partial derivatives with respect to x and y, and maybe sum them in quadrature or something, but it was far easier to just take a simple circle fitting function and use that.
Francesco Pignatelli
on 13 Dec 2022
Image Analyst
on 13 Dec 2022
Sorry, no. It's been 45 years since I took differential calculus. Why won't the fitting it to a circle method work for you?
Francesco Pignatelli
on 13 Dec 2022
J. Alex Lee
on 14 Dec 2022
you can define s using pythagoreas's theorem on pairs of points: dx^2+dy^2 = ds^2 and cumsum on ds.
Image Analyst
on 14 Dec 2022
@J. Alex Lee can you take my demo and build that in so that we have a complete demo?
Francesco Pignatelli
on 19 Dec 2022
J. Alex Lee
on 19 Dec 2022
@Image Analyst, Thanks for the invitation, I did start taking a deeper look at your demo and think there's a good way to put everythign together. @Francesco Pignatelli, I will create another answer when I've worked through an example and include Image Analyst's demo.
I will just briefly say that in studying the examples and application I realized that while the arc length idea as I have laid out is still correct, it is not sufficient to apply to coordinates extracted from an image, because of aliasing. Majority of the work is in smoothing out and resampling such pixel coordinates, not computing the curvature (which is nearly trivial if you have a good description of the xy coordinates).
Image Analyst
on 20 Dec 2022
@Francesco Pignatelli, did "someone else's methodology from a scientific article" give the formula or process they went through to determine the curvature at each (x,y) location?
Francesco Pignatelli
on 23 Dec 2022
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.
Francesco Pignatelli
on 23 Dec 2022
Categories
Find more on Gears in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
