Calculating arclengths…made easy!
Brett's Pick this week is Arclength, by John D'Errico.
First, a nod (and some MATLAB swag!) to Frank Engel, who steered us to John's awesome code! We recently asked users to nominate their favorite File Exchange contributions and Frank jumped in quickly to steer us to Arclength.
How timely! For a medical image processing seminar I recently put together, I wanted to measure the tortuosity of blood vessels. Defining tortuosity as the ratio of arclength to endpoint distance, I segmented an image of the retinal vasculature, skeletonized the image, and broke the vessels into sub-units after detecting branch points, and then sought to measure the length of the sub-units.
After working on the problem for a while, I came up with two reliable methods--after a few misfires. For my first attempt, I sought to reorient vessel segments so that there were no repeated "x-values," and to fit and calculate the length of splines. That approach was unwieldy, and yielded poor results. After playing around some more, I found a couple of reliable, robust approaches. First, I used regionprops to measure the perimeters of the segments, and divided that value by two--works like a charm, since the vessels were skeletonized. Next, I isolated vessel segments using bwlabel, and then calculated the maximum of the bwdistgeodesic (geodesic distance transform). Another success!
John's approach helped me to see where I went astray on my first misguided attempt, and made the solution easy:
img = imread('seg10.png'); %imshow(img);
Here are three approaches:
% Perimeter: perim = regionprops(img,'Perimeter'); perim = perim.Perimeter/2 % Geodesic Distance Transform: [r,c] = find(bwmorph(img,'endpoints')); gdt = max(max(bwdistgeodesic(img,c(1),r(1),'quasi-euclidean'))) % John's arclength: pts = regionprops(img,'pixellist'); pts = [pts.PixelList]; vessellength = arclength(pts(:,1),pts(:,2))
perim = 74.527 gdt = 74.527 vessellength = 74.527
So it appears that all of these approaches work, and isn't it nice to have options? John's code is particularly useful even if your x-y- coordinates aren't necessarily extracted from an image. And because he breaks his path "into a series of integrals between each pair of breaks on the curve," he avoids the problems I had trying to fit a continuous spline to a rotated set of coordinates.
Very nice, John. And thanks for the note, Frank. Swag on the way to both of you!
Keep those suggestions coming. This makes my job a lot easier! :)
As always, comments to this blog post are welcome. Or leave a comment for John here.
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.