Skeletonize a 3D geometry with N-furcation, obtain branch ids and vertices of terminal and N-furcation points, compute MISR

[Natasha]

Hi there,

I'd like to skeletonize a 3D geometry i.e convert to a graph with nodes and edges. Nodes will be the terminal points or junctions where N-furcation occurs. Edges will be the segments that lie between two nodes.

I started by doing the following,

Using the Vessel-Tree.stl file available here, MISR could be obtained via the following code

# vmtkCommand = '''vmtkcenterlines -ifile vesseltree.stl -ofile foo_centerlines.vtp'''

    # p = pypes.PypeRun(vmtkCommand)

    centerlineReader = vmtkscripts.vmtkSurfaceReader()

    centerlineReader.InputFileName = 'foo_centerlines.vtp'

    centerlineReader.Execute()

    clNumpyAdaptor = vmtkscripts.vmtkCenterlinesToNumpy()

    clNumpyAdaptor.Centerlines = centerlineReader.Surface

    clNumpyAdaptor.Execute()

    numpyCenterlines = clNumpyAdaptor.ArrayDict

    #print('numpyCenterlines["Points"] shape = ', numpyCenterlines['Points'].shape)

    print('Point Data Keys: ', numpyCenterlines['PointData'].keys())

    print('Cell Data Keys: ', numpyCenterlines['CellData'].keys())

    print('Point Data Shape: ', numpyCenterlines['PointData']['MaximumInscribedSphereRadius'].shape,' = Number of Points: ', numpyCenterlines['Points'].shape[0])

    print('Point Data Shape: ', numpyCenterlines['PointData']['MaximumInscribedSphereRadius']) 

    print('Point Data EdgeArray: ', numpyCenterlines['PointData']['EdgeArray'].shape)

    print('Point Data EdgePCoordArray: ', numpyCenterlines['PointData']['EdgePCoordArray'].shape)

    print('Point Data CellPointIds: ', numpyCenterlines['CellData']['CellPointIds'])

I could obtain MISR from 

Point Data MISR:  [ 1.63092614  1.6309224   1.64535931 ...,  3.11330085  3.11331409

  3.11300885]

But, I am not sure how to obtain the corresponding vertices/coordinates at which MISR has been computed. but I am not sure how to retrieve the brand id / edge ids and terminal/n-furcation points from CellData.

I would like to generate an output file that would contain branch id, MISR of the point data/ vertices in the respective branch and a connectivity/topology nd-array

that will define the adjacency/connectivity of terminal nodes/vertex in each branch/edge

Any suggestions on how to proceed will be really helpful

Many thanks,

Deepa

[Natasha]

I'd a chance to look at this class, vmtknetworkextraction . Can this feature be used?

[Natasha]

Hi All,

I tried the following and the network skeleton has been generated

 

vmtkCommand = '''vmtknetworkextraction -ifile vesseltree.stl --pipe vmtksurfaceviewer -i @.o'''
p = pypes.PypeRun(vmtkCommand)  

 

I would like to know how to obtain the branch ids and MISR from `p`.Any suggestions?





 
I

[Natasha]

I tried generating a vtkPolyData object as an output and retrieve the branch id and topology information from Point and cell data numpy arrays.

 

vmtkCommand = '''vmtknetworkextraction -ifile vesseltree.stl --pipe vmtksurfaceviewer -i @.o -ofile network.vtp'''

Somehow the output vtkPolyData object network.vtp isn't generated.

Any suggestions will be really helpful.

Deepa 

[Evan Kao]

Hello Deepa, 

Keep in mind that vmtkNetworkExtraction uses a different algorithm for calculating centerlines and will not give you the MISR, but an estimated radius that is different (the details are in this paper).

But, I am not sure how to obtain the corresponding vertices/coordinates at which MISR has been computed.  

The indexes of the points will correspond to the indexes of the point data.  For example, assuming you are using vmtkNetworkExtraction so that MaximumInscribedSphereRadius is replaced by Radius, then:

point_0 = numpyCenterlines['Points'][0,]

radius_at_point_0 = numpyCenterlines['PointData']['Radius'][0]

I tried generating a vtkPolyData object as an output and retrieve the branch id and topology information from Point and cell data numpy arrays.  

I would either skip vmtkSurfaceViewer so that

vmtkCommand = "vmtknetworkextraction -ifile vesseltree.stl -ofile network.vtp"

Or I would just use VMTK's Python API directly.

- Evan

On Wed, Apr 22, 2020 at 8:29 AM Natasha <deepama...@gmail.com> wrote:

I'd a chance to look at this class, vmtknetworkextraction . Can this feature be used?

--
You received this message because you are subscribed to the Google Groups "vmtk-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to vmtk-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/vmtk-users/6895a7e6-4129-492f-8497-d9bb304c8bd0%40googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "vmtk-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to vmtk-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/vmtk-users/5f7d893c-7685-4b66-8614-81ec1574c7ab%40googlegroups.com.

[Evan Kao]

I forgot to mention that another advantage with using VMTK's Python API is that you'll have access to multiple outputs from vmtkNetworkExtraction, the skeletonization and the graph:

They can be accessed like this (assuming surface is the output vmtkSurfaceReader):

networkExtraction = vmtkscripts.vmtkNetworkExtraction()

networkExtraction.Surface = surface

networkExtraction.Execute()

network = networkExtraction.Network

graph = networkExtraction.GraphLayout

- Evan

[Natasha]

Hello Evan,

I tried the following in python,

vmtkCommand = '''vmtksurfacereader -ifile vesseltree.stl
                 --pipe vmtksurfacewriter -ofile foo_net.vtp'''

p = pypes.PypeRun(vmtkCommand)

networkExtraction = vmtkscripts.vmtkNetworkExtraction()
networkExtraction.Surface = 'foo_net.vtp'

networkExtraction.Execute()
network = networkExtraction.Network
graph = networkExtraction.GraphLayout

 I get the following error,

line 23, in <module>

    networkExtraction.Execute()

 networkExtraction.SetInputData(self.Surface)

TypeError: SetInputData argument 1: method requires a VTK object

I'm not sure what's wrong. Could you please have a look?

And is `graph` a vmtkPolyData object that contains cell data and point data?

Also, could you please explain how the skeleton and network image attached in the previous email were created ?

[Deepa]

I tried to use vmtkmarchingcubes before vmtksurfacereader, I still get the same error unfortunately

--
You received this message because you are subscribed to the Google Groups "vmtk-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to vmtk-users+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/vmtk-users/c7f1d567-0baf-4631-8b70-df6ab3c8b53c%40googlegroups.com.

[Evan Kao]

Hello Deepa,

I may not have been clear enough earlier.  The surface needs to be an actual vtkPolyData object that comes after using vmtkSurfaceReader, not the filename.  So something like this:

reader = vmtkscripts.vmtkSurfaceReader()

reader.InputFileName = 'vesseltree.stl'

reader.Execute()

networkExtraction = vmtkscripts.vmtkNetworkExtraction()

networkExtraction.Surface = reader.Surface

networkExtraction.Execute()

network = networkExtraction.Network

graph = networkExtraction.GraphLayout

And is `graph` a vmtkPolyData object that contains cell data and point data?  

Yes, the graph is a vtkPolyData made of vtkLines, which have only two points, while the network is made of vtkPolyLines that are made of multiple points.

Also, could you please explain how the skeleton and network image attached in the previous email were created ?  

I used the latest version of ParaView (installed from binaries). 

- Evan

To view this discussion on the web visit https://groups.google.com/d/msgid/vmtk-users/CAGchuN7hG9Zp%3DG8QLWq6YeDv9QhBBXjYdu7BJpCrOXrkUcM57A%40mail.gmail.com.

[Deepa]

Hi Evan,

Thanks a lot. I could obtain the ''Radius' and 'Topology'.

pprint(network_dsa.CellData['Topology'])

gives

VTKArray([[0, 1],
          [0, 2],
          [0, 2],
          [0, 1],
          [0, 3],
          [0, 3],
          [0, 4],
          [0, 5],
          [0, 5],
          [1, 6],
          [4, 6],
          [1, 5],
          [2, 6],
          [3, 4]], dtype=int64)

pprint(graph_dsa.CellData['Radius'])

VTKArray([ 1.45944192,  1.30663604,  1.30884438,  2.47456279,  1.10002427,
           1.06081075,  1.81926081,  1.3795525 ,  1.42647552,  2.88554481,
           2.29184943,  1.79042397,  1.49456846,  1.21113038])

However, I couldn't understand the information stored in the topology VTKarray. I actually want to use this data to obtain incidence matrix of the `graph`  .

Both, the above-mentioned  VTKArray's have 14 values that correspond to 14 edges in the input file, vesseltree.stl.

t tried to process the CellData like the following. But I am not sure if I understand the output returned and how to derive the incidence matrix from the following. Any suggestions will be helpful.

polyline1 = network_dsa.GetCell(2)  # 1 is the cell id
pprint(polyline1)
for i in range(polyline1.GetNumberOfPoints()):
    pointId = polyline1.GetPointId(i)  # pointId is the point id of the i-th point along the polyline
    pprint(pointId)

To view this discussion on the web visit https://groups.google.com/d/msgid/vmtk-users/CAEhu9hCq6bXHhW5jHGUuyxgi5ozKU8m1OrENAFCspztwYcrLiw%40mail.gmail.com.

[Evan Kao]

Yes, unfortunately the data from the Topology array is not that helpful.  Each cell of the network polydata is a branch.  Each branch is represented as a line/edge on the graph polydata.  The Topology array contains the node/point Ids of the graph.  So for example:

# look at branch 7

topology = network_dsa.CellData['Topology']

print(topology[7]) # gives [1, 5]

node1 = graph_dsa.Points[1, ...]

node5 = graph_dsa.Points[5, ...]

If you are building an incidence matrix B with rows as node ids and columns as edge ids then from the above you can fill two entries as:

B[1, 7] = 1

B[5, 7] = 1

However, the output of the Topology array treats all the end points as 0, which is problematic.  In the end, I suggest ignoring the given topology array and building your own topology array from the graph polydata, i.e., go through every cell/edge/branch and get the point ids, something like this:

for edgeId in range(graph.GetNumberOfCells()):

    edge = graph.GetCell(edgeId)  # this will be a vtkLine, which only has two points

    edgeNode0 = edge.GetPointId(0) 

    edgeNode1 = edge.GetPointId(1) 

Also, you should use vtkCleanPolyData first (or the Clean filter in ParaView) to merge duplicate points before you build the topology array, because vmtkNetworkExtraction does not combine end points for each branch.

- Evan

[Deepa]

Thanks a lot for the answer.

I'm sorry again for the noobish questions

I tried

 print(topology[7]) # gives [1, 5]

This gives [0,5] and not [1,5]. But think this is what you meant by "the output of the Topology array treats all the end points as 0, which is problematic." So I can't really build the icidence matrix like 

B[1, 7] = 1

B[5, 7] = 1.

Next, I tried

for edgeId in range(graph.GetNumberOfCells()):

    edge = graph.GetCell(edgeId)  # this will be a vtkLine, which only has two points

    edgeNode0 = edge.GetPointId(0) 

    edgeNode1 = edge.GetPointId(1) 

returns

'edgeId : 0'
'edgeNode0: 0'
'edgeNode1: 1'

'edgeId : 1'
'edgeNode0: 2'
'edgeNode1: 3'

'edgeId : 2'
'edgeNode0: 4'
'edgeNode1: 5'

'edgeId : 3'
'edgeNode0: 6'
'edgeNode1: 7'

'edgeId : 4'
'edgeNode0: 8'
'edgeNode1: 9'

'edgeId : 5'
'edgeNode0: 10'
'edgeNode1: 11'

and 

so on

:

:

'edgeId : 13'
'edgeNode0: 26'
'edgeNode1: 27'

From what I see here there are 27 nodes in the network. But the actual network doesn't have 27 nodes. Also to create an incidence matrix it's important to have the connection between the edges.

For example,

I will be nice if I can get something like

'edgeId : 0'
'edgeNode0: 0'
'edgeNode1: 1'

'edgeId : 1'
'edgeNode0: 1'
'edgeNode1: 2'

This actually shows the connection between edgeId0 and  edgeId1 (i.e  'edgeNode1: 1' in edgeId 0 is  'edgeNode0: 1' in edgeid1)

And I think this is what you meant by saying , "Also, you should use vtkCleanPolyData first (or the Clean filter in ParaView) to merge duplicate points before you build the topology array, 

because vmtkNetworkExtraction does not combine endpoints for each branch."

So I went ahead and tried

graph = networkExtraction.GraphLayout

clean_graph = vtk.vtkCleanPolyData(graph)

Error:

 clean = vtk.vtkCleanPolyData(graph)
TypeError: method requires a string argumen

I don't really get why a string argument has to be passed. I am trying to pass the ployData in the graph as an input argument.  

Please help me with this

Deepa

[Evan Kao]

Hello Deepa,

From what I see here there are 27 nodes in the network. But the actual network doesn't have 27 nodes. 

 

And I think this is what you meant by saying , "Also, you should use vtkCleanPolyData first (or the Clean filter in ParaView) to merge duplicate points before you build the topology array,   

Yes, exactly.  It actually shows 28 nodes (since indexing starts from 0), 2 for each line. 

Error:
 clean = vtk.vtkCleanPolyData(graph)
TypeError: method requires a string argumen
I don't really get why a string argument has to be passed. I am trying to pass the ployData in the graph as an input argument. 

vtkCleanPolyData is a class, not a function.  I'm not familiar with the technical details over why it admits a string, but in general initializing the class requires no input arguments.  vtkCleanPolyData (and other VTK algorithms) works like this:

cleaner = vtk.vtkCleanPolyData()

cleaner.SetInputData(graph) # set the data to processed

cleaner.Update() # execute the algorithm

cleanedGraph = cleaner.GetOutput() # get the processed data

vtkCleanPolyData has other optional parameters (like setting the tolerance), that you can check in the documentation.

- Evan 

[Deepa]

perfect, thanks!

I did try cleaner.SetInput(graph) after looking at the mention in this post. cleaner.SetInputData(graph) works

Deepa

[Deepa]

Hi Evan,

I want to retrieve the coordinates of the nodes

So, I tried

edgeNode0 = edge.GetPointId(0)
edgeNode1 = edge.GetPointId(1)

edgeNode0_pos = cleanedGraph_dsa.Points[edgeNode0]
edgeNode1_pos = cleanedGraph_dsa.Points[edgeNode1]

edgeNode0_pos is a VTKArray .To convert this to a numpy array, I tried  vtk_to_numpy. But this returns an error:  AttributeError: 'VTKArray' object has no attribute 'GetDataType'. Could you please suggest what's wrong here?

from vtk.util.numpy_support  import vtk_to_numpy

vtk_to_numpy(edgeNode0_pos)

[Deepa]

Hi Evan,

I just noticed 

edgeNode0_pos = cleanedGraph_dsa.Points[edgeNode0]
edgeNode1_pos = cleanedGraph_dsa.Points[edgeNode1]

might not be giving the right coordinates of the nodes.(edgeNode0 and edgeNode1).

pprint(cleanedGraph_dsa.Points.shape) 

still gives (28,3) i.e positions corresponding to 28 nodes. 

Could you please clarify why this happens even after we use vtkCleanPolyData. In the documentation, it is mentioned PointMerging is On by default. It's not clear why cleanedGraph

still has 28 nodes.

Am I missing something here?

Thanks a lot for your support

Deepa

[Deepa]

Hi Evan,

I also have another question. I had a chance to read through one of the paper by Luca, "Automatic Generation of Glomerular Capillary

Topological Organization" . Figure 3 of the paper shows, how the distance between two junctions ( segment length) L is measured.

I'd like to know how these lengths are computed. Is there an in-build class in VMTK? Or is it okay to compute the euclidean distance between edgenode0 and edgenode1?

Thanks again,

Deepa

[Evan Kao]

Hi Deepa,

 

edgeNode0_pos is a VTKArray .To convert this to a numpy array, I tried  vtk_to_numpy. But this returns an error:  AttributeError: 'VTKArray' object has no attribute 'GetDataType'. Could you please suggest what's wrong here?  

VTKArray is derived from a numpy array so there should be no need to convert it.

Could you please clarify why this happens even after we use vtkCleanPolyData  

 

Maybe change the tolerance?  I can't really tell without looking at the code.

I'd like to know how these lengths are computed. Is there an in-build class in VMTK? Or is it okay to compute the euclidean distance between edgenode0 and edgenode1?  

Euclidean distance is probably not what you want.  Try vmtkCenterlineAttributes (see the VMTK tutorial "Geometric Analysis").   

- Evan

[Deepa]

Hi Evan,

Sorry for the delay in response.

Please find the code below

    reader = vmtkscripts.vmtkSurfaceReader()
    reader.InputFileName = 'vesseltree.stl'
    reader.Execute()

    networkExtraction = vmtkscripts.vmtkNetworkExtraction()
    networkExtraction.Surface = reader.Surface
    networkExtraction.Execute()
    network = networkExtraction.Network
    graph = networkExtraction.GraphLayout

    network_dsa = dsa.WrapDataObject(network)

    cleaner = vtk.vtkCleanPolyData()
    cleaner.SetInputData(graph)  # set the data to processed
    cleaner.Update()  # execute the algorithm
    cleanedGraph = cleaner.GetOutput()

    cleanedGraph_dsa = dsa.WrapDataObject(graph)

    tail = []
    head = []
    radius = []

    for edgeId in range(cleanedGraph.GetNumberOfCells()):
        edge = cleanedGraph.GetCell(edgeId)  # this will be a vtkLine, which only has two points

        edgeNode0 = edge.GetPointId(0)
        edgeNode1 = edge.GetPointId(1)

        tail.append(edgeNode0)
        head.append(edgeNode1)

        radius.append((cleanedGraph_dsa.CellData['Radius'][edgeId]))

The input file, ` 'vesseltree.stl ` is attached in the first post of this thread. I could obtain the radius of each branch. What I couldn't obtain is the centerlne distance using vmtkCenterlineAttributes .

I could approximate the lengths using euclidean distance and not using vmtkCenterlineAttributes .

Could you please offer suggestions on how to integrate  vmtkCenterlineAttributes in the above code?

Thank you very much

Deepa

[Deepa]

Hi Evan, 

When I tried to extract the centerline distance, the output of  numpyCenterlines['PointData']. is empty

 reader = vmtkscripts.vmtkSurfaceReader()
    reader.InputFileName = 'vesseltree.stl'
    reader.Execute()

    # centerlines
    # ref: https://mail.google.com/mail/u/0/?tab=rm&ogbl#search/evan/FMfcgxwHMsSvqTVwpztSNlcSFTvfBnvb
    clNumpyAdaptor = vmtkscripts.vmtkCenterlinesToNumpy()
    clNumpyAdaptor.Centerlines = reader.Surface
    clNumpyAdaptor.Execute()
    numpyCenterlines = clNumpyAdaptor.ArrayDict
    pprint(numpyCenterlines['PointData'].keys())

I am not sure how to compute the length of centerlines of each edge of the network from here.

Should one use  numpyCenterlines['Points']?

Thanks,

Deepa

[Evan Kao]

Hi Deepa,

You can use vmtkCenterlineAttributes on any centerline object (e,g, the output of vmtkNetworkExtractor).  It will create an array called "Abscissas".  The Abscissas is the length along each cell/branch.  So the length of the branch will be the abscissas value at the last endpoint of the branch.  See the VMTK tutorial "Geometric Analysis" for more information.

- Evan

[Deepa]

Hi Evan,

Thanks a lot for the response. I could also get to understand `Abscissas` from the tutorial.

I could find only the following fields associated with  networkExtraction = vmtkscripts.vmtkNetworkExtraction().

I also check the fields associated with `network` obtained using,   network = networkExtraction.Network. I couldn't find vmtkCenterline or vmtkCenterlineAttributes .

If there is an example that shows how to retrieve vmtkCenterlineAttributes  from the output of vmtkNetworkExtractor that will be of great help.

Many thanks for your time and attention,

Deepa

[Deepa]

Hi Evan,

This is a polite reminder.

Deepa

[Evan Kao]

Hello Deepa,

Sorry I forgot to answer this.  vmtkCenterlineAttributes is another class, so you shouldn't be looking for it as an attribute of networkExtraction but instead apply it.  Actually, I think I misspoke previously. You should use vmtkCenterlineGeometry instead, which will give you a "Length" for each branch:

network = networkExtraction.Network

# I would also do some resampling and smoothing of the network here using vmtkCenterlineResampling and vmtkCenterlineSmoothing

centerlinesGeometry = vmtkscripts.vmtkCenterlineGeometry()

centerlinesGeometry.Centerlines = network

centerlinesGeometry.Execute()

Then you can use vmtkCenterlinesToNumpy to extract the data.  The length of each branch is stored in the cell data attributes as "Length".

Here the branches are colored by length.  The crude red arrow points to cell 0, which is highlighted in the spreadsheet view and has a length of 26.7205.

- Evan

[Deepa]

Hi Evan,

Thanks a lot. 

May I know how you saved vesseltree_network_abscissas.vtp as output? I'd also like to know if Paraview is used for visualizing the contents of the vtp file.

Deepa

[Evan Kao]

Yes, I used ParaView to visualize the centerlines. Centerlines are technically vtkPolyData, the same as surfaces, so you can use vmtkSurfaceWriter to save them to file.

[Deepa]

Thanks you very much, Evan

I tried this

reader = vmtkscripts.vmtkSurfaceReader()
reader.InputFileName = 'vesseltree.stl'
reader.Execute()

# network extraction

networkExtraction = vmtkscripts.vmtkNetworkExtraction()
networkExtraction.Surface = reader.Surface
networkExtraction.Execute()

network = networkExtraction.Network

centerlinesGeometry = vmtkscripts.vmtkCenterlineGeometry()
centerlinesGeometry.Centerlines = network
centerlinesGeometry.Execute()

surfacewriter = vmtkscripts.vmtkSurfaceWriter()
surfacewriter.Surface = centerlinesGeometry.Centerlines
surfacewriter.OutputFileName = "vesseltree_network.vtp"
surfacewriter.Execute()

After loading  vesseltree_network.vtp in Paraview

May I know how to view the table that you shared in the second last email and color branches? 

[Evan Kao]

In the top right corner of the render window are some buttons for manipulating the window (red box).  You can split the view horizontally or vertically:

You can then choose what type of viewer you want.  "Render View" is the default, used for looking at data in 3D.  To view the numerical data, click on the last option, Spreadsheet View (red box):

The top of the Spreadsheet View will have drop down menus to choose which dataset you want to view (red box) and which data field (point, cell, or field data) you want to examine (green box):

Highlighting the cell/branch you want in the Spreadsheet View (green arrow) will also highlight the corresponding branch in the Render View in pink (red arrow):

This can be tough to see though because the centerlines are pretty thin, which is why I used the Tube filter to make the centerlines easier to see.  But the Tube filter creates multiple new cells for each polyline of the network, so highlighting each branch requires selecting multiple cells:

- Evan

[Deepa]

Thank a ton for your tremendous support, Evan

I really like the spreadsheet view.

However, I find the values of `Length` are slightly off from the values displayed in the snapshot shared by you. I've generated the output file using the code

shared in the second last email. 

Also, could you please let me know how the legend has to be enabled? I used the search tab in the color map editor to look for Length. But for some reason the 

the length field is not accessible.

Also, while using the tube filter may I know how the branch ids are displayed in the spreadsheet?

 

Deepa

[Evan Kao]

However, I find the values of `Length` are slightly off from the values displayed in the snapshot shared by you. I've generated the output file using the code
shared in the second last email.

That's probably because I resampled and smoothed the centerlines.  Whether it's better to use the raw or smoothed network is up to you.

Also, could you please let me know how the legend has to be enabled? I used the search tab in the color map editor to look for Length. But for some reason the 
the length field is not accessible.

For Paraview questions, it's probably best to read the Paraview guide or watch a tutorial.   As for coloring, you can find the Coloring section in the properties tab.  There is a dropdown menu to choose which array you want to color with (red box).  When you pick it, the Legend should automatically pop up.

Also, while using the tube filter may I know how the branch ids are displayed in the spreadsheet?  

I used the Generate Ids filter and under Cell Ids Array Name I entered "BranchId".

- Evan