# Curvature of a 2D image

49 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?
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.

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
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

Image Analyst on 11 Dec 2022
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.

### Categories

Find more on Time Series Events 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!