Help in measuring diameter of vessel

Hello,

I need help in understanding why my code doesn't output quite what I thought it would. I am trying to measure the diameter of a vessel in a 512x512 image by measuring the diameter at several points and averaging the measurements. Here is my code

          segmentLength = round(size(image,1)/segments);
          [rightrow, rightcol] = find(image==3);
          rightwall = [rightcol rightrow];
          for j = 1:segmentLength:(size(image,1)-segmentLength)
              midpt = j+round(segmentLength/2);
              [~, midptcol] = find(image(midpt,:)==2);
              if isempty(midptcol)==0
                  [d, p] = min((rightwall(:,1)-midptcol).^2+(rightwall(:,2)-midpt).^2);
                  %p is the row index of rightwall containing the nearest point
                  d = sqrt(d);
                  diameter = [diameter d];
              end
          end

I changed the value along the right side of the vessel to be 3 in order to identify what wall is the wall on the right and used 20 for my value of segments. When I graph the points to see what this code has done, the points with shortest distance don't seem to actually be those points (as far as I can tell with my eye) and so are not quite the diameter at that point. Can you help me understand what I am missing? Here is the image where the points on the right don't match up with where I would expect them to:

Thanks! Michael

Answer 1

[Image Analyst]

Michael, try this:

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 = 20;

% Read in mat file.
storedStructure = load('binaryimage.mat');
binaryImage = storedStructure.BWfinal;
% Display the image.
subplot(2, 2, 1);
imshow(binaryImage, []);
title('Original Binary Image', 'FontSize', fontSize, 'Interpreter', 'None');

% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') 
drawnow;

% Fill the image.
binaryImage(1,:) = true;
binaryImage(end,:) = true;
binaryImage = imfill(binaryImage, 'holes');
% Display the image.
subplot(2, 2, 2);
imshow(binaryImage, []);
title('Original Binary Image', 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;

% Get the Euclidean Distance Transform.
binaryImage(1,:) = false;
binaryImage(end,:) = false;
edtImage = bwdist(~binaryImage);
% Display the image.
subplot(2, 2, 3);
imshow(edtImage, []);
title('Distance Transform Image', 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;

skelImage = bwmorph(binaryImage, 'skel', inf);
skelImage = bwmorph(skelImage, 'spur', 40);
% There should be just one now.  Let's check
[labeledImage, numLines] = bwlabel(skelImage);
fprintf('Found %d lines\n', numLines);
% Display the image.
subplot(2, 2, 4);
imshow(skelImage, []);
title('Skeleton Image', 'FontSize', fontSize, 'Interpreter', 'None');

% Measure the radius be looking along the skeleton of the distance transform.
meanRadius = mean(edtImage(skelImage))
meanDiameter = 2 * meanRadius
message = sprintf('Mean Radius = %.1f pixels.\nMean Diameter = %.1f pixels',...
  meanRadius, meanDiameter);
uiwait(helpdlg(message));

[Michael J]

This worked fantastic, thank you very much! And so, I will accept it as my answer. However, I am actually writing this code for several images and so the step where you cut down the extra lines using the spur operation may cut down too little or too much. I've tried writing a bit of code to try and get to one longest line without me needing to look at the resulting image and checking to see if it is only one line, but am having an issue with bwlabel. Here is my code:

    numLines = 0;
    while numLines ~=1
        skelImage = bwmorph(skelImage, 'spur', 1);
        [labeledImage, numLines] = bwlabel(skelImage);
    end

The issue here is that the loop ends immediately as the number of lines is already 1 after the creation of the skeleton in skelImage = bwmorph(im, 'skel', inf);

I've also tried subtracting branchpoints to then delete the extra branches, but this leads to false branches containing parts of the main spine. If you have any further insight on how to approach this issue I would greatly appreciate it, thank you again for all of your help!

[Image Analyst]

An alternate way would be to use regionprops to get the orientation. Then use imrotate to rotate the vessel vertical. Then you may have to find a way to chop off the triangles at the top and bottom. Then just go along row by row finding the first and last point

for row = 1 : rows
    thisRow = rotatedImage(row, :);
    leftCol = find(thisRow, 1, 'first');
    rightCol = find(thisRow, 1, 'last');
    width(row) = rightCol - leftCol;
end

[Sandip Sadhukhan]

Can you please send me the "binaryimage.mat" file? I want to test it and then accordingly modify to cater to one of my requirements. Thanks in advance.

[Image Analyst]

@Sandip Sadhukhan, it's attached.

[Sandip Sadhukhan]

Thanks a lot for your quick response.

Answer 2

[Image Analyst]

I think the usual way to do this, and more accurate than your way, is to get a binary image of the vessel, then compute the Euclidean Distance Transform with bwdist(), then to skeletonize the binary image and use it to mask the EDT image to get an image with values only along the centerline of the vessel. This gives you all the radii along the axis of the vessel. Then you simply take the mean or histogram or whatever you want.

If you have an image of just the edges then you need to seal the top and bottom edge, then call imfill, then erase the top and bottom rows.

Something like this:

binaryImage(1,:) = true;
binaryImage(end,:) = true;
binaryImage = imfill(binaryImage, 'holes');
binaryImage(1,:) = false;
binaryImage(end,:) = false;
edtImage = bwdist(binaryImage);
skelImage = bwmorph(binaryImage, 'skel', inf);
meanRadius = mean(edtImage(skelImage));
meanDiameter = 2 * meanRadius;

The skeleton will look like a "Y" near the ends of the vessel, so you can get the main spine of the skeleton by using bwmorph() to find the branchpoints and setting those to false and then use bwareafilt(binaryImage, 1) to get the longest branch.

Try it and let me know how it goes. Attach the initial binary image if you get stuck.

[Michael J]

Thank you for the fast response Image Analyst! I am working with your suggestions but keep getting a meanRadius of 0. I think the issue may be in the skelImage step as the skeleton does indeed look like "Y" at the vessel ends but there are a few "Z" 's in the middle of the skeleton. When trying to get rid of these extra branches, using bwmorph(skelImage, 'branchpoints') and then turning those points false and using them to then get the longest branch with bwareafilt(skelImage,1), I am left with a "Z" as opposed to the main spine. I have attached the binary image in question, it is the result of some imaging segmentation functions to try and isolate the vessel from the background.

[Felipe Suntaxi]

Dear Image analyst I'm doing a similar analysis but in 3D. I would like to calculate the diameter of some vessels in a binary image. Any suggestions?

[Chris Loizou]

@ImageAnalyst, i really need your help with something im working on, is there a way to contact you?

[Image Analyst]

You can post a question on Answers, and if I get time and it doesn't take too much time, I might answer. If it's a major project, I'll probably just refer you to relevant published articles or bibliographies, like this one, or possibly to a university that I know that you could hire to do the work. I usually don't spend more than 5 minutes answering - 30 minutes tops, but that would have to be something that looks really fun, interesting, or challenging, or something I might think would make a nice general purpose demo.

[Jyoti Jethe]

@Image Analyst Hi, I am new to image processing in matlab. Please help to find out the diameter of blood vessel in the given image. There are network of vessels, I tried using skeletonization but still unable to get a proper diameter value and way to plot it (or may be histogram of the value of diameter). Following is the piece of code and the image attached is rerode image:

sktn_img=2*bwdist(~rerode);

%imshow(sktn_img);

skeleton=bwmorph(rerode,'skel',inf);

diameter_img=sktn_img.*double(skeleton);

Please help