Extracting Vessel Dimensions from 2D Fluorescence Microscopy Images in MATLAB

I am working with fluorescence microscopy images of blood vessels in the brain, stored as HDF5 files (.h5). My goal is to extract vessel dimensions (diameters) from each frame in a time series and analyze how they change over time.

Currently, I am using MATLAB to:

1. Load the .h5 file and extract individual frames.
2. Apply Gaussian filtering and contrast enhancement to improve vessel visibility.
3. Perform edge detection (Canny) and skeletonization to extract vessel centerlines.
4. Compute vessel diameters using the Euclidean distance transform from the skeleton.

However, I am running into several issues:

• Edge detection is inconsistent and picks up artifacts inside the vessel.
• The skeleton sometimes extends beyond the vessels or is not continuous.
• Extracted diameters appear constant, even when the vessel visibly changes size.

Current Approach in MATLAB

Here’s the code I am using to manually select vessels and compute diameters:

vars = dloc.saved_vars;

img = dloc.load_var("reference_img");
img = squeeze(img(:,:,1));

frame_idx = 1000;
frame = double(h5read(dloc.path, '/recording', [1,1,1,frame_idx], [512, 512,1,1]));

frame_smooth = imgaussfilt(frame, 2);
% Create a normalized & contrast-enhanced frame for Canny
frame_enhanced = mat2gray(frame_smooth);
frame_enhanced = imadjust(frame_enhanced, stretchlim(frame_enhanced, [0.01, 0.99]));

%% --- Frangi Vesselness Filter ---
[I_vessel, ~] = FrangiFilter2D(frame_enhanced, struct( ...
    'FrangiScaleRange', [1 5], ...
    'FrangiScaleRatio', 2, ...
    'FrangiBetaOne', 0.5, ...
    'FrangiBetaTwo', 15, ...
    'BlackWhite', true));

bw_vessel = I_vessel > 0.1;
bw_vessel = bwareaopen(bw_vessel, 100);
bw_vessel = imclose(bw_vessel, strel('disk', 3));

skeleton = bwskel(bw_vessel);
skeleton = bwmorph(skeleton, 'spur', 30);

distance_map = bwdist(~bw_vessel);
diameter_map = distance_map * 2;

% Extract vessel diameters along skeleton
[y, x] = find(skeleton);
diameters = arrayfun(@(i) diameter_map(y(i), x(i)), 1:length(x));

%% --- Edge Detection for Overlay ---

edge_thresholds = [0.1, 0.4]; 
edges = edge(frame_enhanced, 'Canny', edge_thresholds);

%% --- Visualization: Frangi + Edges + Skeleton + Diameters ---
figure;
imshow(frame_enhanced, []); 
hold on;

contour(edges, [1,1], 'g', 'LineWidth', 1);
plot(x, y, 'r.', 'MarkerSize', 5);

scatter(x, y, 15, diameters, 'filled');
colormap(parula); 
c = colorbar;
c.Label.String = 'Diameter (pixels)';
title(['Edges + Skeleton + Diameters (Frame ', num2str(frame_idx), ')']);

%% --- Plot Diameter Profile ---
figure;
plot(1:length(diameters), diameters, 'b.-');
xlabel('Skeleton Point Index');
ylabel('Vessel Diameter (pixels)');
title('Vessel Diameter Profile Along Skeleton');

Resulting edges and vessel diameter detection

As you can see, the edges are picked up but not always, and the diameters are calculated on no real basis. They don’t seem to be particularly telling.

My question is:

How can I ensure accurate edge detection that captures the actual vessel boundaries and calculate their diameters?

Answer 1

[William Rose]

@Alessio,

Consider Better Skeletonization in the File Exchange.

Please upload a frame of the raw .h5 video image.

You write "As you can see, the edges are picked up but not always, and the diameters are calculated on no real basis. They don’t seem to be particularly telling." What you say I can see is not obvious to me. The grayscale and color images in your post do not match either of the figures made by the code you posted. Please post the images that were generated by the code that you posted.

The vessel diameter is a function of both position along the vessel and time: d(x,t), where x is a measure of distance along the axis of the vessel, and t is time (frame number). When you say "Extracted diameters appear constant, even when the vessel visibly changes size", do you mean that d(x) appears constant or d(t) appears constant, or both?

The code you posted does not include a measure of position along the axis of the vessel ("x", in my comment above). This seems like something you may want to do. If you assign an "x" value to each pixel along the centerline, then you will be able to plot diameter versus x, which may be useful.

A number of years ago I wrote code in Labview to extract time-varying vessel diameter from ultrasound videos of human arteries. It worked pretty well. The output, d(t), reqired filtering to remove glitches and gaps, due to occasional changes in insonation (illumination) that occur with ultrasound.

[Alessio]

Hello @William Rose and thank you for the answer, the frame of the raw .h5 video image is the one I posted above, it's an output of the stabilized image that is obtained by the h5 file and then retrieved with

img = dloc.load_var("reference_img");

I don't understand what you mean by the grayscale and color images do not match the figures, that is the code I used and the results it gave me. I have tweaked it a bit in between yesterday and now, and I will update the question to match the changes.

For hte extracted diameters appear constant, what I meant was that it didn't seem to matter that the vessels branches out and/or becomes bigger. Please let me know if this is clearer, I just started working with this data a few weeks ago and I admit I get lost still.

[William Rose]

@Alessio,

If you make a posting that contains mistakes or inconsistencies, please don't edit the original post, if that renders the responses of others non-sensical. We all make mistakes. If you go back and change your post, then people trying to help you cannot figure out what really happened and the discussion makes no sense. The better approach is to post a comment with your revisions. In this case, you went back and changed your code and changed the figures - at least the second figure. In your original posting, the titles and other features of the images did not match the code you posted.

In your original posting, the second image had a skeleton "line" (not a straight line) running down the approximate middle of the central vessel. The skeleton line split at the Y in the vessel, and a branch of the skeleton continued along each of the two vessel branches. Now the skeleton line in the middle is not visible. I think it is because nyour code includes:

plot(x, y, 'r.', 'MarkerSize', 5);          % plot skeleton as red dots
scatter(x, y, 15, diameters, 'filled');     % plot circles with diameter 15

The circles plotted in the second line of code above cover up the skeleton. The second line plots circles, all with size 15. The circle colors are given by "diamter" in your code. Is that your intention?

Please post the files and data necessary to allow others to run your code. That includes the file reference_img, and the contents of variable dloc, and the file which is referenced by dloc.path.

[William Rose]

I see that you extract the red component of the image, with

img = squeeze(img(:,:,1));

That makes sense, if the vessels are labelled with DiI and you use a rhodamine filter, or a similar setup that produces a good signal in the red channel. Otherwise, you could experiment with other approaches such as summing all three colors to make a grayscale image:

img=sum(img,3);

Maybe that will have improved contrast.