Skeletonize 3d cylinders

I'm having trouble figuring out a reliable way to draw lines marking the midpoints of mostly hollow cylinders in a 3d volume. Demo volume attached. I've been trying mostly morpholigical operations (dilate/erode/close/open) to try to finagle out a 'rod' from the cylinder's extent. It's not terrible, but a few get lost and a few get smushed together into a bit of nonsense every time. It's also inflexible and needs handholding with different sizes. Is there a clever way of doing this, hopefully that also avoids spending 95% of the time on super slow dilate operations?

The cylinders are open-ended tubes, so I can't just fill them then erode. I would love to be able to fill them, though. The cylinders also sometimes have junk in them. This makes it nonviable to try to subtract an expanded cylinder from a filled one to generate the rod. I've also tried 3d convolution/template matching and various clustering/watershed methods of breaking smushed regions apart, but it didn't get me any further.

Is there some neat trick that can fill the cylinders, or otherwise identify them with very high accuracy and specificity?

[William Rose]

@Carson Purnell, tubeplot in the file exchange was useful to me in making 3D images of blood vessel networks.

[William Rose]

@Carson Purnell, what am I supposed to do with democylinders.mat? I tried

load('democylinders.mat');
imshow(raw)

but imshow() produced a bunch of error messages. I tried opening democylinders.mat in the Matlab volume viewer, but that produced an "invalid file format" error.

[DGM]

raw appears to be a volumetric image, but it's 'double', scaled [0 1725], so I don't know how volume viewer would treat it. I'm not used to doing anything technical with volumetric images, but it looks like this:

... which would have the following mean projections:

[Matt J]

With volshow(raw>0),

Answer 1

[Matt J]

The key problem is that some of the cylinders are touching other cylinders. For those that don't, we can easily separate them with regionprops3 and fit cylinders to them using this FEX download:

https://www.mathworks.com/matlabcentral/fileexchange/87584-object-oriented-tools-for-fitting-conics-and-quadrics

Once you have a cylinder fit, it is easy to derive the medial line. Maybe you can refine the code below with your watershed code to get a better separation of the cylinders

load('democylinders.mat');
BW=imopen(raw>0,strel('sphere',3));


cv=regionprops3(BW,'VoxelList');
cv=cellfun(@transpose, cv{:,1}','un',0);


fobj=cellfun(@cylindricalFit,cv);
fobj=fobj(abs([fobj.radius]-9)<2); %discard cylinder fits giving outlier radii

for i=1:N
  surf(fobj(i)); hold on
  medialLine(fobj(i));
end; hold off; grid on
axis equal

function medialLine(cyl)
%Plot medial line of a cylinder fit

   d=cyl.R(:,1)'*cyl.height/2;
   p1=cyl.center+d;
   p2=cyl.center-d;
   XYZ=num2cell([p1;p2],1);
   plot3(XYZ{:},':b','LineWidth',3);

end