# Curvature of a 2D image

54 views (last 30 days)
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 CommentsShow 2 older commentsHide 2 older 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.
Francesco Pignatelli on 24 Dec 2022
Hi @J. Alex Lee, yes of course! My code is:
clusterdata = extractSingleShotOH(namel,i); %%Here I import a specific frame from the data acquired
image=clusterdata.imageData;
image=flip(image,2);
image=flip(image,1);
%% Filtering and Edge Extraction
BW = medfilt2(image,[2,2]);
BW = mat2gray(BW);
level = graythresh(BW);
binaryImage = BW > level;
binaryImage = imfill(binaryImage, 'holes');
binaryImage = bwareafilt(binaryImage, 10);
BWmean = rescale(binaryImage);
boundaries = bwboundaries(BWmean);
Regarding the first bullet point, yes, sorry for the misunderstaning, with "3rd order" I mean cubic interpolation and using the interpolation options you mentioned, a continuous second derivative is ensured. This is the reason why I am currently using spline function.

Sign in to comment.

### 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:
1. define arc length coordinate for each xy coordinate using pythgoreas. This is implemented in my attached function "arclen"
2. 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
Arc length method took 0.0922560 s IA circle fit method 0.0111560 s Starting parallel pool (parpool) using the 'Processes' profile ... Connected to the parallel pool (number of workers: 2). pared Driscoll-Manjunatha took 27.5067510 s
##### 7 CommentsShow 5 older commentsHide 5 older comments
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

Sign in to comment.

### More Answers (1)

Image Analyst on 11 Dec 2022
##### 15 CommentsShow 13 older commentsHide 13 older 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.
Francesco Pignatelli on 23 Dec 2022
@J. Alex Lee thanks for your reply. I am not 100% sure if I have understood correct but yes I am quite puzzled about the 3rd order method based on direct spline interpolation and I am trying to ffix it somehow. Anyway I am having a look inyto your functions and I am trying to take all the different ingredients trying to implement the formula I am trying to compute. Thus, thanks for your scripts :) I will attach mine soon.

Sign in to comment.

### Categories

Find more on Geometric Transformation and Image Registration in Help Center and File Exchange

R2018b

### Community Treasure Hunt

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

Start Hunting!