Tortuosity of blood vessels

Hello,

I am trying to calculate the Tortuosity of blood vessels. (Vessel tortuosity is calculated as the sum of branch lengths divided by the sum of the euclidean distance between their end points).

I have the following questions:

1. How to calculate the lengths of both actual branches and the imaginary straight lines between nodes
2. How do I mark the vessel branches (orange) and the branch nodes (yellow) as shown in the picture.
3. I am calculating "spinelength" as the sum of all pixels. How do I calculate individual branch lengths?
4. Any suggestions in preprocessing would be appreciated.

I followed https://www.mathworks.com/videos/medical-image-processing-with-matlab-81890.html video and read the https://www.mathworks.com/matlabcentral/answers/522378-how-to-measure-the-diameter-along-the-length-of-an-object#answer_429784 post. Any help would be highly appreciated. Thank you.

clc;
clear;
close all

% Read the image
I=imread('VAD.png');
figure,imshow(I)

%convert it to gray scale
I_gray=rgb2gray(I);
%Sharpen the image
b = imsharpen(I_gray,'Amount',8);

h = fspecial('average', [3 3]);
b = imfilter(b, h);

%choose brighter objects
Bina=b>150
figure,imshow(Bina);

se = strel('cube',3)
erodedBW = imerode(Bina,se);

%Remove small objects from binary image
BW2 = bwareaopen(Bina,100)
figure,imshow(BW2);

skelImage = bwskel(BW2, 'MinBranchLength', 10);

MinBranchLength = round(sum(skelImage(:))/2)
skelImage = bwskel(BW2,'MinBranchLength',MinBranchLength);
figure,imshow(skelImage)

endpointImage = bwmorph(skelImage, 'endpoints');
[rows, columns] = find(endpointImage)

spineLength = sum(skelImage(:))
straightLineDistance = sqrt((columns(2) - columns(1))^2 + (rows(2) - rows(1))^2)
tortuosity = spineLength / straightLineDistance

Answer 1

[Image Analyst]

Call bwmorph() to get the branchpoints then use that to erase the branch points. Then label each curve and call bwmorph() to get the endpoints to get the straight line distance. Pretty easy, in fact you've probably already done it by now. You might have gotten something like (untested):

bp = bwmorph(mask, 'branchpoints');
mask(bp) = false; % Erase branchpoints.
[labeledImage, numRegions] = bwlabel(mask);
for k = 1 : numRegions
    thisRegion = ismember(labeledImage, k);
    endpoints = bwmorph(thisRegion, 'endpoints');
    % Get coordinate
    [r, c] = find(endpoints);
    euclideanDistance(k) = sqrt((r(end)-r(1)).^2 + (c(end)-c(1)).^2);
    area(k) = sum(thisRegion);
    tortuosity(k) = area(k) / euclideanDistance(k);
end

How does that compare to your code?

Answer 2

[darova]

Since you already have a skeletonized image

• dilate your image at nodes to separate each lines

BW1 = imdilate(BW,zeros(3));        % dilate image with matrix 3x3

• use bwlabel to separate each line