{"id":3506,"date":"2012-04-27T08:11:08","date_gmt":"2012-04-27T13:11:08","guid":{"rendered":"https:\/\/blogs.mathworks.com\/pick\/?p=3506"},"modified":"2022-05-29T20:18:00","modified_gmt":"2022-05-30T00:18:00","slug":"calculating-arclengths-made-easy","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/pick\/2012\/04\/27\/calculating-arclengths-made-easy\/","title":{"rendered":"Calculating arclengths…made easy!"},"content":{"rendered":"
\r\n <\/introduction>\r\n

Brett<\/a>'s Pick this week is Arclength,<\/a> by John D'Errico<\/a>.\r\n <\/p>\r\n

First, a nod (and some MATLAB swag!) to Frank Engel<\/a>, who steered us to John's awesome code! We recently asked users to nominate their favorite File Exchange contributions<\/a> and Frank jumped in quickly to steer us to Arclength.\r\n <\/p>\r\n

How timely! For a medical image processing seminar I recently put together, I wanted to measure the tortuosity<\/a> of blood vessels. Defining tortuosity as the ratio of arclength to endpoint distance, I segmented an image of the retinal\r\n vasculature, skeletonized<\/a> the image, and broke the vessels into sub-units after detecting branch points<\/a>, and then sought to measure the length of the sub-units.\r\n <\/p>\r\n

After working on the problem for a while, I came up with two reliable methods--after a few misfires. For my first attempt,\r\n I sought to reorient vessel segments so that there were no repeated \"x-values,\" and to fit and calculate the length of splines.\r\n That approach was unwieldy, and yielded poor results. After playing around some more, I found a couple of reliable, robust\r\n approaches. First, I used regionprops<\/tt><\/a> to measure the perimeters of the segments, and divided that value by two--works like a charm, since the vessels were skeletonized.\r\n Next, I isolated vessel segments using bwlabel<\/tt><\/a>, and then calculated the maximum of the bwdistgeodesic<\/tt><\/a> (geodesic distance transform). Another success!\r\n <\/p>\r\n

<\/p>\r\n

John's approach helped me to see where I went astray on my first misguided attempt, and made the solution easy:<\/p>

img = imread('seg10.png'<\/span>);\r\n%imshow(img);<\/span><\/pre>
 <\/p>\r\n   
Here are three approaches:<\/p>
% Perimeter:<\/span>\r\nperim = regionprops(img,'Perimeter'<\/span>);\r\nperim = perim.Perimeter\/2\r\n\r\n% Geodesic Distance Transform:<\/span>\r\n[r,c] = find(bwmorph(img,'endpoints'<\/span>));\r\ngdt = max(max(bwdistgeodesic(img,c(1),r(1),'quasi-euclidean'<\/span>)))\r\n\r\n% John's arclength:<\/span>\r\npts = regionprops(img,'pixellist'<\/span>);\r\npts = [pts.PixelList];\r\nvessellength = arclength(pts(:,1),pts(:,2))<\/pre>
perim =\r\n       74.527\r\ngdt =\r\n       74.527\r\nvessellength =\r\n       74.527\r\n<\/pre>
So it appears that all of these approaches work, and isn't it nice to have options? John's code is particularly useful even\r\n      if your x-y- coordinates aren't necessarily extracted from an image. And because he breaks his path \"into a series of integrals\r\n      between each pair of breaks on the curve,\" he avoids the problems I had trying to fit a continuous spline to a rotated set\r\n      of coordinates.\r\n   <\/p>\r\n   
Very nice, John. And thanks for the note, Frank. Swag on the way to both of you!<\/p>\r\n   
Keep those suggestions coming. This makes my job a lot easier! :)<\/p>\r\n   
As always, comments to this blog post<\/a> are welcome. Or leave a comment for John here<\/a>.\r\n   <\/p>