File Exchange

image thumbnail

Skel2Graph 3D

version 1.22.0.1 (10.8 KB) by Philip Kollmannsberger
Calculate the network graph of a 3D voxel skeleton.

24 Downloads

Updated 05 Dec 2018

GitHub view license on GitHub

MATLAB code to derive the network graph of a 3D voxel skeleton

This function converts a 3D binary voxel skeleton into a network graph described by nodes and edges. The input is a 3D binary image containing a one-dimensional voxel skeleton, generated e.g. using the "Skeleton3D" thinning function available on MFEX. The output is the adjacency matrix of the graph, and the nodes and links of the network as MATLAB structure.

Usage:

[A,node,link] = Skel2Graph(skel,THR)

where "skel" is the input 3D binary image, and "THR" is a threshold for the minimum length of branches, to filter out skeletonization artifacts.

A is the adjacency matrix with the length of the links as matrix entries, and node/link are the structures describing node and link properties.

A node has the following properties:

- idx List of voxel indices of this node
- links List of links connected to this node
- conn List of destinations of links of this node
- comX,comY,comZ Center of mass of all voxels of this node
- ep 1 if node is endpoint (degree 1), 0 otherwise

A link has the following properties:

- n1 Node where link starts
- n2 Node where link ends
- point List of voxel indices of this link

A second function, "Graph2Skel3D.m", converts the network graph back into a cleaned-up voxel skeleton image.

An example of how to use these functions is given in the script "Test_Skel2Graph3D.m", including a test image. In this example, it is also demonstrated how to iteratively combine both conversion functions in order to obtain a completely cleaned skeleton graph.

Any comments, corrections or suggestions are highly welcome. If you include this in your own work, please cite our publicaton [1].

Philip Kollmannsberger 09/2013, 01/2016
philipk@gmx.net

[1] Kollmannsberger, Kerschnitzki et al.,
"The small world of osteocytes: connectomics of the lacuno-canalicular network in bone."
New Journal of Physics 19:073019, 2017.

Cite As

Kollmannsberger, Kerschnitzki et al., "The small world of osteocytes: connectomics of the lacuno-canalicular network in bone." New Journal of Physics 19:073019, 2017.

Comments and Ratings (30)

hao

Hi Julien,
I use it to generate many points.How to connect them to lines.
Thank you!

Can it exactly compute branchpoints and endpoints? The present bwmorph3 gives them, but it is unclear how they are connected? Thanks

alphin j

Can I use it for a 3D point cloud?

Cody Crosby

Hi Julien,
Disclaimer: I did not write the code. Generally, voxels are considered to be the 3D representations of pixels (i.e. volume as opposed to area). Hope that helps!

Me back again, I have a question about your algorithm:
Is there a difference between a "voxel" and "pixels"? If yes, which is it exactly?

Thank you!

Hello, I would like to implement a colormap on my 3D graph where the color depend on the size of the vessel. Would that be possible you think?

I have found this program very helpful in our study of trabecular bone, but we have been sidetracked by other projects and haven't been able to explore our research using this program...hopefully we can get back to it soon.

But, I also wanted to comment on Jianhao question on converting stl files/data to 3D binary images. You can do a search in MatLab's file exchange for 'stl to 3D volume' which will give you a bunch of possible tools to use. I found the 'Mesh voxelization' tool very useful

Hello.I am trying to use your function to do analysis for 3d printing model.But I meet trouble on converting the stl document into the required 3d binary image.Any tools that can be used?Thanks.

I made some bug fixes and improvements to the code on github. This should address some of the issues mentioned below:

- empty link and conn arrays resulting from short branches should now be fixed
- nodes and links now have a "label" variable to identify multiple skeletons in a volume

Remaining issues or "features":

- node degree 2 occurs in case of multi-voxel nodes with two links emanating (usually skeletonization artifacts)
- XY coordinates are swapped, this is by design but can easily be converted.

Nodes can have multiple "indexes", if multiple voxels belong to that node. Regarding STL files, these are polygon mesh files and need to be converted to binary images first.

Hi all,
I would like to use this code or the other skeletonize code. We are using vessels segmented in an STL file type, but we keep getting the error array exceeds max array size. I am wondering if anyone else has had this problem and if there is an easy fix?

Thanks for your help~!

Hello. This function is extremely useful. Thank you very much for putting it together.

I have a question regarding my output from the "node" struct.
It seems that some "idx" entries contain different numbers of indexes (anywhere from 1 to 9 different indexes, it seems), and some of them appear to be right-justified in the variables workspace, instead of left as per the majority, for some reason. Some entries in "links" are simply empty arrays, as are the "conn" entries of that particular node.

Is there more information you could provide as to why these issues are occurring?

I've copied the full "node" struct of my skeleton below:

ep: comz: comy: comx: conn: links: idx:
----------------------------------------------------------------------------------------------------------------------
0 9.50000000000000 32.5000000000000 23.5000000000000 [31,3] [1,2] [62591;69918]
0 16 30 37 [32,3] [3,4] 114283
0 17.2500000000000 31.2500000000000 43 [1,2,4,35] [2,4,5,6] [121813;121912;121913;129138]
0 21 32 57 [3,5,8] [5,7,8] 151626
0 22.5000000000000 27.5000000000000 49.5000000000000 [4,37,40,39] [7,9,10,11] [158648;165973]
0 24 27 56 [41,7,43] [12,13,14] 173405
0 26 27 59 [6,8,45] [13,15,16] 188258
0 26.5000000000000 30.5000000000000 67.5000000000000 [4,7,46,11] [8,15,17,18] [188663;195988]
0 28 25.5000000000000 64.5000000000000 [34,11,44,47] [19,20,21,22] [202916;203014]
0 30 22 68 [36,11] [23,24] 217472
0 30.8571428571429 25.8571428571429 68 [8,9,10,13,51,42] [18,20,24,25,26,27] [217867;218067;225094;225194;225293;225392;232620]
0 32 6 69 [49,13] [28,29] 230739
0 34 28.4444444444444 73.6666666666667 [11,12,52,54,53,14] [25,29,30,31,32,33] [232919;240246;240345;240445;247772;255197;255297;255397;262723]
0 38 29 78 [13,15] [33,34] 277575
0 42.6666666666667 31.3333333333333 80 [14,58,17] [34,35,36] [307475;314900;314999]
0 51.5000000000000 59 84.7500000000000 [56,17,18] [37,38,39] [377077;377176;384403;384501]
0 53 39 82 [15,16,19] [36,38,40] 389944
0 53 64.3333333333333 79.3333333333333 [16,59,57] [39,41,42] [392416;392417;392515]
0 61 42 84 [17,20,23] [40,43,44] 449643
0 67 55 81 [19,21,22] [43,45,46] 495477
0 67 61 78 [20,60,61] [45,47,48] 496068
0 73 62 84 [20,64] [46,49] 540723
0 76 38.2500000000000 82.5000000000000 [19,27,62,25] [44,50,51,52] [553295;560620;560720;567946]
0 83 17.5000000000000 93.2500000000000 [65,26,67,70] [53,54,55,56] [603102;610626;610627;617952]
0 85 30 80.4000000000000 [23,66,28,68,75] [52,57,58,59,60] [619226;626651;626652;626752;633976]
0 86 24 88.5000000000000 [24,69,71,28,29] [54,61,62,63,64] [626065;633392;633491;641014]
0 86 53 81 [23,74,73] [50,65,66] 636354
0 88 29 85 [25,26,77] [58,63,67] 648832
0 91 16 92 [26,72,30] [64,68,69] 669827
0 95.5000000000000 13.5000000000000 91 [29,76] [69,70] [699328;706654]
1 7 33 18 1 1 47736
1 14 29 33 2 3 99330
1 16 29 37 [] [] 114184
1 17 13 76 9 19 120064
1 19 27 39 3 6 136263
1 21 7 74 10 23 149168
1 21 25 43 5 9 150919
1 21 27 49 [] [] 151123
1 23 19 45 5 11 165177
1 23 25 47 5 10 165773
1 23 27 53 6 12 165977
1 25 9 57 11 27 179049
1 25 19 49 6 14 180031
1 27 17 53 9 21 194687
1 27 23 55 7 16 195283
1 28 37 74 8 17 204113
1 29 17 57 9 22 209541
1 29 21 67 [] [] 209947
1 31 1 69 12 28 222819
1 31 5 69 [] [] 223215
1 31 17 61 11 26 224395
1 31 34 75 13 30 226092
1 33 36 77 13 32 241142
1 38 25 75 13 31 277176
1 38 28 77 [] [] 277475
1 45 62 86 16 37 332825
1 53 68 78 18 42 392811
1 54 15 94 15 35 395005
1 55 59 74 18 41 406766
1 67 58 75 21 47 495768
1 69 67 77 21 48 511511
1 71 26 75 23 51 522300
1 73 63 85 [] [] 540823
1 74 64 82 22 49 548344
1 76 14 93 24 53 558255
1 77 24 86 25 57 566663
1 79 13 96 24 55 580434
1 81 17 73 25 59 595657
1 81 19 83 26 61 595865
1 82 13 93 24 56 602706
1 87 19 91 26 62 640423
1 91 11 95 29 68 669335
1 91 59 80 27 66 674072
1 93 47 81 27 65 687735
1 94 31 68 25 60 693563
1 99 15 89 30 70 729125
1 117 46 77 28 67 865832

Avik Mondal

Hi there,
I've been using skel2graph with trabecular networks. I remember reading in one of your papers (The small world of osteocytes) that there are nodes of degree 1 which are endpoints, and nodes of degree 3 and up by definition. However, I have been occasionally finding nodes of degree 2. I'm analyzing a large dataset and so to keep things computationally efficient I usually cut the image in to much smaller VOI's so I assumed that nodes of degree 2 were a side effect of that. I was recently wondering however if it was possible that nodes of degree 2 occur when a node consists of a cluster of voxels (rather than the 1 voxel as is expected from skeletonization) and that cluster has two links emerging from it?

Avik Mondal

Nevermind my previous comment, I simply didn't use a skeletonized binary matrix

Avik Mondal

I was wondering about Test_Skel2Graph3D. I've tried to convert a few skeleton's into graphs and I occasionally run into the error "Index exceeds matrix dimensions." The big for loop in the code is bounded by the length of node2. However, for my skeletons, my link2 has less elements than node2. Do you know how I can fix the code to handle this?

Thanks for reporting - this should be fixed now (pull request by @kevinbbb87 merged today), at least at the two obvious locations.

J C R

Very nice tools!
There seems to be a bug at "Skel2Graph"
at
"s2n=zeros(w*l*h,1);
s2n(nhi(:,14))=1:length(nhi);
"
length(nhi) should give the number of rows, but this fails when the number of rows is less than the number of columns. Simply changing it to:
"
s2n(nhi(:,14))=1:length(nhi(:, 1));
"
does the trick.

Kind regards

Yoann Atlas

After some tests, I have finally just changed the loop condition "for i=1:num_realnodes" by "for i=1:length(node)" in order to consider all the nodes and to keep the single links. It works for the moment.

Yoann

Yoann Atlas

Hello

Thank you for your code, it's a really good job!

It seems there is still an issue. When the skeleton is only composed of one link and two end-points (what you called "ep"), you skip visiting these nodes (the loop "for i=1:numrealnodes"). So the structure is losing all single links.
Maybe it's a little bit dirty but I copy/paste the loop and change the condition "for i=1:cc3.NumObjects" to get link beetween two end-points.

Hope it could help!
Yoann

Hi Philip!
Thanks for this really cool function. I am using this function to find the end points of the branches in a 3D binary image of dendrite. I have two queries: (1) Is there a way to decide a threshold automatically for removing the artifacts since skeleton3D gives too many branches than required. Or perhaps a threshold inside skeleton3D to not get too many branches along small rough surfaces. (2) After getting nodes from skel, is there a way to decide the nodes that lie on end points of the branches. Nodes with ep=1 are too many, I just need the ones at ends of big branches and not along the small rough surfaces.
Thanks
Sahil

The issues with the "while" loop mentioned below should now be fixed, and some documentation of the resulting structures has been added.

To adhere to common practice, the endpoints of branches are now also included as nodes with degree = 1. They can easily be filtered out to reproduce the behavior of the previous version.

Philip

Helton Leal

Hi Phillip, I'm working with neuron segmentation in 3D images, your code is helping alot, but I'm stuck understanding the Node and Link structure in order to advance.

If I have the index of a pixel, how do I get which Node is the one who represents that pixel? I tried iterating along all nodes and test membership in idx, it didnt work for all pixels.

Hi Philip,

I'm currently working with your submission and I'm stuck with something I can't understand. Would you help me? It is related to the while loop in Test_Skel2Graph3D that filters spurious nodes. My question is why generating the graph from the skel and the skel from the graph alternatively is supposed to remove those nodes? I also found out that sometimes repeating this filtering operation do not remove any nodes after a few iterations, until it suddenly does.

I would really appreciate your help with this.
Thanks a lot!

José Ignacio

Carlos

Hi Philip, thanks a lot for the script, it works perfectly!
I just wonder why in some cases (with own skel volumes) the Test_Skel2Graph3D gets stuck in the while loop ((min(cellfun('length',{node2.conn}))<3)). When this happens, I skip the loop and it shows the figure properly (or at least it seems so)

Carlos

I'm referring to Test_Skel2Graph3D file

Sébastien

Thanks, now it is working fine! Great contribution! One question: In case multiple unconnected skeletons are detected could each of their constituent nodes/branches be assigned an unique skeleton label?

Nicolas

Hi Sébastien,

thanks for reporting this. The error appears whenever the skeleton volume has non-zero values on its boundary. I updated the function accordingly. Now, if there are positive voxels on the boundary, they are deleted (and lost).

This should not be a big problem however, since the skeleton is topologically not defined on its borders anyway (the boundary has no neighbourhood).

Please let us know if it works now.

Philip

Sébastien

Hi,

This function does not work properly. I keep on getting the same error:

"Error using sub2ind (line 56)
Out of range subscript."

You can for instance trigger this error with this example input skeleton:

S=logical(zeros(3,3,3));
S(2,2,1)=1;
S(2,2,2)=1;
S(2,2,3)=1;

but it essentially failed with same error message for most of the larger skeletons I generated by calling Skeleton3D.

Updates

1.22.0.1

updated description

1.22.0.0

Linked to Skeleton3D

1.21.0.0

Linked to github

1.2.0.0

fixed typo

1.2.0.0

Fixed typo

1.2.0.0

- zero padding added, outermost layer is now included
- branch ends are included as 1-nodes
- added documentation of the node/link structure
- more intuitive test script, removed bugs (infinite loops)

1.1.0.0

Function returned an error if the skeleton volume was not padded with zeros. This is now corrected.

MATLAB Release Compatibility
Created with R2010b
Compatible with any release
Platform Compatibility
Windows macOS Linux
Acknowledgements

Inspired by: Skeleton3D

Discover Live Editor

Create scripts with code, output, and formatted text in a single executable document.


Learn About Live Editor