Find the endpoints of 3D skeleton

I have images of vessels which are of size 316x316x856. I have used the Skeleton3D function which is available from FileExchange to extract the skeleton, and so the result is a binary volume of the same size. Now however I don't know how to find the endpoints, seeing as the bwmorph command only works for 2D. I tried calculating the Euclidean distance of every point in the skeleton to the other, in a double loop, and keep the maximum distance, but the results I get I think are wrong. Are there any other ways to correctly compute the endpoints?

Answer 1

[Sean de Wolski]

Convolve the 3d skeleton using convn() with a rubik's cube of ones. Anywhere in this new image that is equal to one or two should be an endpoint.

endpts = ismember(convn(double(skel3d),ones(3,3,3),'same'),[1 2]);

[Iason Solomos]

Thank you for the prompt answer. I do get quite a large number of endpoints though, ~8million out of the 85 million total in the image. The number of nonzero elements in the binary vessel volume is around 17m. And the nonzero elements in the skeleton is 7 million. Does it make sense for the nonzero elements to be 1 million more in the endpoint image than in the skeleton image?

[Sean de Wolski]

This is also going to get spur points so depending on the skeleton algorithm it might be picking up lots of spurs. You could repeat the above 3-4 times which would shorten the skeleton a little but remove the spurs.

Alternatively: There is another option that is more involved that uses the geodesic distance to give you two end points for each vessel that maximize the length of it. If you're interested in this(?) I can try and put some code together. It would involve grabbing all voxels for a specific object. Selecting all of the end/spur points for it and then calculating the geodesic distance between each of them and selecting the two with the largest distance.

[Iason Solomos]

Yes you are right, there are probably a lot of spurs in my image. The geodesic distance idea is interesting. Here is my code:

    try
    load red
catch
    warning('File, red.mat is not present. Currently generating file and will save it for future uses...')
    second = 'D:\Data\(20141210_11_09_27)__Perivascular_drainage_XYTZ_ch_2.tif';
    channel = resizeChannel(second,'red',0.5);
end

red = channel;
red = padarray(red, [30 30 30], 0);
redOne = double(red)/double(max(red(:)));
vessel_volume = (redOne>0.01);
L = bwlabeln(vessel_volume);
props = regionprops(L, 'area');
areas=cat(1,props.Area);
[~, index]=max(areas);
vessel_volume =(L==index);
skel = Skeleton3D(vessel_volume);
endpts = ismember(convn(double(skel),ones(3,3,3),'same'),[1 2]);
ind = find(endpts);
[i1,i2,i3]=ind2sub(size(endpts),ind);
maxdist = 0;
for i = 1:i1
    for j = 1:i1
        if sqrt(((i1(i)-i1(j))^2+(i2(i)-i2(j))^2+(i3(i)-i3(j))^2)) > maxdist
        maxdist = sqrt(((i1(i)-i1(j))^2+(i2(i)-i2(j))^2+(i3(i)-i3(j))^2));
        Point1 = i;
        Point2 = j;
        end
    end
end

startp = [i1(Point1),i2(Point1),i3(Point1)];
endp = [i1(Point2),i2(Point2),i3(Point2)];
  end

[Sean de Wolski]

Do you have an example zip file you can attach with the MAT file?

[Sean de Wolski]

Also, right idea and good first pass!

You're calculating regular distance, not geodesic. Look at the 'SubImage' and 'PixelList' option in regionprops for pulling out each object, and then bwdistgeodesic for a selected end point.

[Sean de Wolski]

Actually, can you attach a MAT file in a zip file with just skel stored.

[Iason Solomos]

Sorry for coming back to you late, yes I can do that. Will send my short code as well.