How to calculate the distance between blood vessel edges?

The project is to calculate the diametre of a blood vessel in multiple points and to plot the vessel lenght (x-axis) by diameter (y-axis), basically to calculate the diameter of the vessel in multiple points based on the length between the two edges. In the code it only calculates it for two points, and also the 2 lines for the vessel edges aren't correctly defined in the "Skeleton Image". The image I worked with is also attached.

How can I finish my project?

% Read in mat file.
storedStructure = load('imagine_binară.mat');
binaryImage = storedStructure.BWfinal;

% Fill the image.
binaryImage = imclose(binaryImage, strel('disk', 30));
binaryImage(1,:) = true;
binaryImage(end,:) = true;
binaryImage = imfill(binaryImage, 'holes');

% Get the Euclidean Distance Transform.
binaryImage(1,:) = false;
binaryImage(end,:) = false;
edtImage = bwdist(~binaryImage);

% Skeletonize the image.
skelImage = bwmorph(binaryImage, 'skel', inf);
skelImage = bwmorph(skelImage, 'spur', 40);

% Utilizează regionprops pentru a obține proprietățile regiunilor conectate.
stats = regionprops(skelImage, 'Centroid', 'MajorAxisLength');

% Extrage coordonatele centrelor și lungimile majore ale regiunilor.
centroids = cat(1, stats.Centroid);
majorAxisLengths = [stats.MajorAxisLength];

% Calculează diametrul pentru fiecare regiune.
diameters = majorAxisLengths;
% Utilizarea regionprops pentru a calcula proprietățile regiunilor conectate
stats = regionprops(skelImage, 'MajorAxisLength');
% Extrage lungimea maximă a segmentului
maxLength = max([stats.MajorAxisLength]);
% Afișează skeletonul imaginii.
subplot(1,2,1)
imshow(skelImage, []);
impixelinfo;
title('Skeleton Image');

% Plotează distribuția diametrului.
figure;
plot(maxLength, diameters, 'o');
xlabel('Length');
ylabel('Diameter');
title('Diameter Distribution along Skeleton');

Answer 1

[Image Analyst]

You need to invert your binary image. The vessel should be white:

binaryImage = ~storedStructure.BWfinal;
binaryImage = bwareafilt(binaryImage, 1); % Take largest blob only.

An alternate method would be to scan across the image getting the top row and bottom row of the image then subtracting then and averaging.

[Ana Brostic]

Thank you for you answer @Image Analyst. My region of interest is the vessel, but I only need its walls in order to measure the diameter along the vessel from the top wall to the bottom wall in multiple places, that is why I worked with them inverted. Can you help me with this please?

[Image Analyst]

Yes. But you only want the interior opening ("lumen") of the vessel, right? You don't want to find the walls and then get the center line of the walls, right? I'll be out a lot today so I might not get to it right away. I think the original binary image is not good so can you attach the original gray scale image and your segmentation algorithm?

[Ana Brostic]

Yes, I only want the interior. I will attach a picture to better describe my task. So I want to define the walls, to only have the contour that are marked with blue, which is what I tried to do with the skeletization. And then I need to find the diameter, so basically the distance between the walls in multiple points which is marked with the red lines. At the end, I need to do a graph with the distances depending on the wall length.I hope now it is more clear. Thank you so much for your help!

[Image Analyst]

Take the two largest blobs with bwareafilt. Then label them with bwlabel. Then use ismember to extract the upper one and lower one individually. Then scan column by column finding the lower edge of the upper blob and thetop edge of the lower blob. Then subtract them. Something like

mask = bwareafilt(mask, 2); % Take 2 largest blobs.
labeledImage = bwlabel(mask);
% Find upper blob
props = regionprops(mask, 'Centroid');
xy = vertcat(props.Centroid)
[~, minBlobNumber] = min(xy(:, 1)) % minBlobNumber will be either 1 or 2 depending on which blob is left-most.
upperBlob = ismember(labeledImage, minBlobNumber)
lowerBlob = mask & ~upperBlob; % Lower blob is the other one.
[rows, columns] = size(mask);
% Preallocate arrays for top and bottom edges.
topEdge = nan(1, columns);
bottomEdge = nan(1, columns);
for col = 1 : columns
    thisCol = upperBlob(:, col);
    t1 = find(thisCol, 1, 'last');
    if ~isempty(t1)
        topEdge(col) = t1;
    end

    thisCol = lowerBlob(:, col);
    t2 = find(thisCol, 1, 'last');
    if ~isempty(t2)
        bottomEdge(col) = t2;
    end
end
widths = bottomEdge - topEdge;
x = 1 : columns;
% Get rid of nans
badColumns = isnan(widths);
x(badColumns) = [];
widths(badColumns) = [];
% Plot
plot(x, widths, 'b-', 'LineWidth', 2);
grid on;
xlabel('x')
ylabel('Widths')
title('Width vs. x')

[Ana Brostic]

@Image Analyst Thank you for your quick response. I will try this method and get back to you with my result. Happy holidays!