How to measure the diameter along the length of an object?

Hi All,

I have a dataset of many small cells (one cell per image). These cells have different (irrigular) shapes, roundess and orienatioans. I need to measure cells' diamaters at each point along their length. I followed the solutions suggested by "Image Analyst" in these posts:

1- https://www.mathworks.com/matlabcentral/answers/126064-how-to-calculate-the-diameter-of-a-vessel

2- https://www.mathworks.com/matlabcentral/answers/273506-help-in-measuring-diameter-of-vessel

I tested the algorithm on two images as shown below. The top one is a sample of my dataset images and the bottom one is an image that I created to test and understand the algorithm better. However, it seems that I missing something and my reults do not completly make sense to me. In the bottom image, I wanted the algorithm to give me higher values for the dimater of the middle section of the rectangle.However, it seems that bwdist() function gives me lower values for those sections. Could you please help me with that?

Also, how can I get rid of those "big" fluctuations in the results of the top image? I tried to smooth() the original Image and/or apply moving average function to the results. However, non of them worked properly.

Thank you for your help and time!

Answer 1

[Image Analyst]

See my demo to get the average width of a curvy blob. Adapt as needed.

% This demo computes the tortuosity and average width of a snake-like shape by dividing the actual longest spine distance by the straight line distance between the endpoints.
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;

%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
folder = pwd;
baseFileName = 'tortuosity.png';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
	% The file doesn't exist -- didn't find it there in that folder.
	% Check the entire search path (other folders) for the file by stripping off the folder.
	fullFileNameOnSearchPath = baseFileName; % No path this time.
	if ~exist(fullFileNameOnSearchPath, 'file')
		% Still didn't find it.  Alert user.
		errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
		uiwait(warndlg(errorMessage));
		return;
	end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
	% It's not really gray scale like we expected - it's color.
	% Use weighted sum of ALL channels to create a gray scale image.
	grayImage = rgb2gray(grayImage);
end
% Display the image.
hFig = figure;
subplot(1, 4, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;

%--------------------------------------------------------------------------------------------------------
%    SEGMENTATION
% Binarize the image
binaryImage = imbinarize(grayImage);
% Fill holes, if any.
binaryImage = imfill(binaryImage, 'holes');
% Take the largest blob only.
binaryImage = bwareafilt(binaryImage, 1);
% Display the binary image.
subplot(1, 4, 1);
cla;
imshow(binaryImage, []);
title('Binary Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');

%--------------------------------------------------------------------------------------------------------
%    SKELETONIZATION
% Now skeletonize
skelImage = bwskel(binaryImage, 'MinBranchLength', 10); % First guess.  Strangely bwskel() doesn't seem to have any short spurs, no matter what MinBranchLength is.
% skelImage = bwmorph(binaryImage, 'skel', inf); % First guess.  Will have short spurs.  Shouldn't use this, but going to just for the demo.
% Display the original skeleton image, with annoying spurs
subplot(1, 4, 2);
imshow(skelImage, []);
title('Original Skeleton Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');

% First get the area of the skeleton and take half of it to get the MinBranchLength.
% This should give use the longest spine path only, and get rid of all spurs (at least for THIS demo image).
MinBranchLength = round(sum(skelImage(:))/2)
% Now take the skeleton again using that min branch length we just computed.
skelImage = bwskel(binaryImage, 'MinBranchLength', MinBranchLength);
% Display the skeleton image.
subplot(1, 4, 3);
imshow(skelImage, []);
title('Skeleton Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');

% Overlay the skeleton on the original image.
subplot(1, 4, 1);
imshow(imoverlay(binaryImage, skelImage, 'r'));
title('Binary Image with Skeleton Overlaid', 'FontSize', fontSize, 'Interpreter', 'None');

%--------------------------------------------------------------------------------------------------------
%    TORTUOSITY COMPUTATION
% Find the endpoints
endpointImage = bwmorph(skelImage, 'endpoints');
[rows, columns] = find(endpointImage)
% Draw a red line between them.
subplot(1, 4, 3);
hold on;
plot(columns, rows, 'r-', 'LineWidth', 2);

% Find the spine length.
% (Would be more accurate to find the spine length by traversing it.)
spineLength = sum(skelImage(:))
% Compute the straight line distance.
straightLineDistance = sqrt((columns(2) - columns(1))^2 + (rows(2) - rows(1))^2)
% Compute the tortuosity
tortuosity = spineLength / straightLineDistance
caption = sprintf('Skeleton Image.  Tortuosity = %f', tortuosity);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

%--------------------------------------------------------------------------------------------------------
%    AVERAGE WIDTH COMPUTATION
% Compute the Euclidean Distance Transform
edtImage = bwdist(~binaryImage);
subplot(1, 4, 4);
imshow(edtImage, []);
title('Distance Transform', 'FontSize', fontSize, 'Interpreter', 'None');
% Extract the distances only along the skeleton image to extract only radii along spine.
radii = edtImage(skelImage); 
% These are the half widths.  Multiply by 2 to get the full widths.
averageWidth = 2 * mean(radii)
caption = sprintf('Distance Transform.  Mean Width = %.1f', averageWidth);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

[EhsanD]

Hi,

Thank you so much for your time helping me. Actually, I was exactly using the same approach that you mentioned in your answer above:

edtImage = bwdist(~binaryImage);
radii = edtImage(skelImage);

If I understand correctly, bwdist() returns the distance from the nearest non-zero pixel for every pixel in the image. However, it seems that it doesn't give me the correct diameter always. I created this sample image below to explain my question better.

What I want to measure are the "lengths of the blue arrows" along the centerline of this image. These arrows are perpendicular to the centerline and start from one edge and go to the other edge of the shape. However, if you see the point that it is shown with a red asterisk, the code above gives me the length of the red arrow at that point (as it is actually the shortest distance from the asterisk point to the first white pixel at the end of red arrow calculated by bwdist() function). In my application, the correct diameter at that point is the length of the green arrow (which is also perpendicular to the centerline). I would highly appreciate it if you could please help me find the right solution to measure the length of these blue/green arrows. Thank you!

[Image Analyst]

How much of a difference do you think it will make, really? Probably not much. If you want to try though, you can find the skeleton like I did. Then progress along it and fit a stretch to a polynomial, like a quadratic. Then get the slope and take minus one over that slope to get the slope of the perpendicular line. Use improfile() to obtain a profile along that perpendicular cross section. I've done it before when I wanted to get the slope along an edge in the image. It's not easy but if you want to try, go ahead. I can't give you my code because it's to intricately woven into my program and it would take a lot of handholding to talk you through it. But if it's that important to you, you'll figure out how to do it - I already gave you the algorithm so just code it up. Otherwise take the easy route with the code I gave you. It probably won't be that much of a difference anyway.

[EhsanD]

Thanks for your reply and guiding me to find a solution. However, I didn't quite understand why you said the results might not be that different. As part of my research, I am interested in characterizing the morphology of the vessels, especially around the swollen regions. If we have a big swelling, then the relative difference between the length of the green arrow and the red arrow would become bigger. Am I right?

Also, I don't know how to handle situations just like what happened in the bottom image of my original post. Using the solution based on bwdist(~binaryImage) and edtImage(skelImage) functions I get lower diameter values for the swollen regions (as shown in the bottom plot). How can I fix this problem? Thank you!

[Image Analyst]

I told you already. Reread my last comment starting with the third sentence.

[Tri Rowstenkowski]

@Image Analyst Hello! I am trying to calculate the tortuosity of blood vessels- https://www.mathworks.com/matlabcentral/answers/740432-tortuosity-of-blood-vessels. I have a problem with the calculation of individual branch length. Can you please help me?

Answer 2

[Image Analyst]

As you can see, the main axis is not so obvious in general. For some amorphous splat-shaped blob, what is the length or width? See, it's not so clear cut.

1. One way is to use bwskel() to get the main spine and bwdist() to get the distance transform. Then multiply the two and take the mean along the spine and double it to get the width.
2. Another way assumes that the blob is basically a rectangle. So you could use bwferet() to get the maximum caliper diameter, then use regionprops() to get the area. Then get the mean width by dividing the area by the max Feret diameter, since a width of a rectangle is the area divided by its length..

[EhsanD]

Thank you very much for your prompt reply. You are right. In general, my images are more like blood vessels that have a bulging at some point along their lengths. I am interested in measuring the diameter of the vessel at every single point along its centerline (around the swollen region). In order to find the centerline, I used bwmorph(bwin, 'thin', Inf) function. Then, I multiplied the bwdist() output image with the centerline image and got those plots that I showed in my original post above. However, it seems that I am missing something. I think I need something to give me the "maximum" distance (from one edge to the other edge of my vessels) and perpendicular to the centerline. Does that make sense? How can I smooth the output results (reduce those fluctuations)? Thanks!