Skeleton3D: Parallel medial axis thinning of a 3D binary volume
This code calculates the 3D medial axis skeleton of an arbitrary 3d binary volume. It is an optimized MATLAB implementation of the homotopic thinning algorithm described in . We developed it to quantify the network of cell processes in bone , but it should work on images of any tubular or filamentous structures. An example volume (testvol.mat) is included, along with an example script (Test_Skeleton3D.m). Any comments, corrections or suggestions are highly welcome.
skel = Skeleton3D(bin)
where "bin" is a 3D binary image, and "skel" the resulting image containing only the skeleton voxels, or
skel = Skeleton3D(bin,mask)
to mask all foreground voxels in 'mask' from skeletonization, e.g. to preserve certain structures in a image volume.
For additional cleanup, e.g. pruning of short branches, please use my Skel2Graph3D package on MATLAB File Exchange.
This code is inspired by the ITK implementation by Hanno Homann  and the Fiji/ImageJ plugin by Ignacio Arganda-Carreras . If you include this in your own work, please cite our original publicaton .
Philip Kollmannsberger 09/2013
 Ta-Chih Lee, Rangasami L. Kashyap and Chong-Nam Chu
"Building skeleton models via 3-D medial surface/axis thinning algorithms."
Computer Vision, Graphics, and Image Processing, 56(6):462–478, 1994.
 Kollmannsberger, Kerschnitzki et al.,
"The small world of osteocytes: connectomics of the lacuno-canalicular network in bone."
New Journal of Physics 19:073019, 2017.
Kollmannsberger, Kerschnitzki et al., "The small world of osteocytes: connectomics of the lacuno-canalicular network in bone." New Journal of Physics 19:073019, 2017.
How can I use it for getting curves from a point cloud(set of dense curves)?
in the file p_oct_label.m
you checked the adjacency of octant 8 according to 17 twice.
idx = cube(:,17) == 1;
cube(idx,17) = label(idx);
cube(idx,:) = p_oct_label(4,label(idx),cube(idx,:));
I think you can delete one of them.
How much effort will it cost to generate surface skeletons like in ?
Thanks for answering me Philip, I got it! Now, I wanna do something more. I have another file with the diameter size value of the vessel in it (in bmp format) and which is of the same matrix size as the 3D-Skeleton that I obtain. My idea is to color the points of the skeleton with a sort of colormap where hot values means big vessel size and blue small vessel size. Any idea to help me in this challenge would be greatfull for me. Thanks!
Thank you, Phillip!
@Alex: the easiest way is to compute the distance transform on the original volume using bwdist() and to multiply the result with the final skeleton. This should give you the skeleton with a local thickness estimate stored as pixel intensity.
@Julien: this sounds like "skel" is not of type logical - try "skel = logical(skel)", or replace line 42 by "cands=zeros(width,height,depth);".
I have the following error appearing using a 3D matrix as input and full of 0 and 1... Do you have any idea on how I can fix that? I really not see why this message appear since my matrix is already supposed to be logical and comes from slices of a CT-scan...
"Error using false
Input following 'like' is not a logical array.
Error in Skeleton3D (line 42)
Is there a way to have the code output an estimate of local thickness at each point along the skeleton?
great function, thanks!
I updated the code by merging some bug fixes and performance improvements made by David Dreher some time ago.
@Eric Chadwick: can you send me a small test dataset or code how you generate the data?
Works well, except that I have some strange artifacts when I run this on a volume with intersecting vessels. I am working on processing a much larger vasculature scan, but have created 3 intersecting cylinders to practice on. At the intersection, the function preserves the outline of the cylinders as well as the centreline. Do you know why this might be?
This seems like a very useful function! Can it deal with shapes like a torus (doughnut), that have holes?
To remove short branches, please try my Skel2Graph3D script also here on FEX.
A slight directional bias can result from the fact that the medial axis is not uniquely defined in a voxel structure, but the error should never be more than +/- 1 voxel.
In the last update, a H1 line has been added and the code was cleaned up a bit.
Be useful if minimum branch length (noise removal) option was included
Seems somehow biased to the left direction. That is, Skeleton3D(BW) and fliplr(Skeleton3D(fliplr(BW))) do not produce the same result
I have an MRA vessel that I'm trying to find the centre-line for, running it through the Skeleton3D code doesn't seem to isolate anything from the matrix. Any thoughts on whats going wrong?
Almost 5 stars. Wonderfully useful tool, but the help is a little sparse (it would be useful to have a H1 line so we know what the parameters and outputs are), and it looks like the function once had global variables which are now (thankfully) replaced but their global designation remains.
Since this is a port from other code I bet it can be vectorised in places for speed but I've found it useful so far - thanks.
Linked to github
Removed designation of global variables, and added H1 line for help. For cleaning noise such as short artifacts, please use my Skel2Graph3D package.