VMTK center lines extraction: how to accelerate it and how to extract center lines from small vessel?

VMTK is really a good tool. I am using VMTK to extract center lines from vessel. And I have two questions:

  1. Currently, it would cost 8s for one calculation. Is there any method to accelerate it? Can we use gpu in VMTK?

  2. I find it failed for VMTK to extract center line from small vessel. Just like the following picture:

VMTK can not extract center line for the two points (source and target point). Is there any method to extract center line for small vessel?

Thanks for any suggestion!

For reference, the question was posted and discussed on VMTK project page - https://github.com/vmtk/vmtk/issues/354

Unrelated to VMTK, but you can try to use SGEXT to get a graph representation of the airway, nodes correspond to end points and branch points, and edges correspond to center lines of the vessels. SGEXT has been used for to to quantitative analysis in airways, cardiovascular and polymer networks by myself with good results.
Python wrappings for the whole library are provided, and you can install it with pip install sgext. If you want to give it a try, and have questions, please open issues there.
Hope it helps.

Thank you for providing SGEXT. It seems to have similar features as VMTK, but I guess there are internal differences that may be important for some applications. Do you have any comparison to VMTK? What are the main differences, advantages/disadvantages? What features are available without restrictive (commmercial/GPL) license?

We could probably make SGEXT available with a nice GUI in 3D Slicer but we would only invest time into this if there are clear advantages.

I haven’t made the comparison. VMTK is built around meshes, the goal of SGEXT is to take a spatial graph representation of the object (from an image) and work with it as a quantitative analytical tool.
It uses a thinning algorithm, and this one is LGPL from DGtal ( I contributed it from recent research of Digital Topology). Thinning is at the core of SGEXT, but it’s the only part that uses DGtal.

SGEXT includes modules doing reconstructions of graphs, graph comparisons, analysis of branching points (generations), even simulations with pair forces between nodes.

Up to you to decide if it is well suited for Slicer, but its integration won’t take long at all, the only extra dependency would be DGtal.
There are functions to convert the spatial graph to a vtu representation already, and to voxelize the graph into a label image. The library is 95% wrapped in python.

EDIT: The other dependency, which I guess is not that light is Boost, with components: graph, serialization and filesystem.

Thank you for the information, it is useful. Working directly from labelmap input can be useful.

What interfaces does it have with VTK (or ITK)? How the graphs, meshes, etc. are available in Python?

The graphs are from the boost graph library (BGL), and they are internally wrapped in python.

The interaction with ITK is provided in the module sgext.itk to convert back and forth from ITK images.
The images used for SGEXT are internally wrapped, and there are only two types, Binary (unsigned char), and Float (float), both 3D.
If you want to mix pipelines with ITK and sgext, you need to use the member function in the sgext images to_pyarray and from_pyarray and use the analogous functions in ITK. Basically using numpy as a bridge.
This interoperability doesn’t exist with VTK, VTK is used internally for some modules (graph comparisons and location of points), to visualize, and IO operations with the filesystem.

Oh, I will edit the previous post, the other dependency is boost: graph, serialization, and filesystem.

It is not easy to use sgext.

I think the readme of sgext should be update. The example code is:

import sgext

input_filename="/path/to/binary_image.nrrd" # or any format that ITK can read
out_folder="/path/to/output_folder" # or any format that ITK can read

sgext.scripts.thin(input=input_filename,
        out_folder=out_folder,
        foreground="black",
        skel_type="end",
        select_type="first",
        persistence=2,
        visualize=False,
        verbose=True)

However, the method thin accept:

    input: BinaryImageType
        input binary image.

    table_folder: str
        Location of the DGtal look-up-tables for simplicity and isthmusicity,
        for example simplicity_table26_6.zlib.
        These tables are distributed with the sgext package. Use the variable
        'sgext.tables_folder'.

Then, my code is:

itkImg = itk.GetImageFromArray(mask)
sgext.scripts.thin(input=itkImg, 
                   table_folder= './',
                   foreground="black",
                   skel_type="end",
                   select_type="first",
                   persistence=2,
                   visualize=False,
                   verbose=True
                   )

The mask is a numpy matrix (500500400, uint8). And a bug is reported:

TypeError: thin(): incompatible function arguments. The following argument types are supported:
    1. (input: itk::SmartPointer<itk::Image<unsigned char,3> >, skel_type: str, select_type: str, tables_folder: str, persistence: int = 0, input_distance_map_image: itk::SmartPointer<itk::Image<float,3> > = None, profile: bool = False, verbose: bool = False, visualize: bool = False) -> itk::SmartPointer<itk::Image<unsigned char,3> >

Invoked with: kwargs: input=<itkImagePython.itkImageUC3; proxy of <Swig Object of type 'itkImageUC3 *' at 0x00000241E1D999F0> >, table_folder='./', foreground='black', skel_type='end', select_type='first', persistence=2, visualize=False, verbose=True

Two points:

  • If you read the table_folder help, it states that you need to use table_folder=sgext.tables_folder
  • You need to convert the numpy array into a sgext image. Use:
sgextImg = sgext.itk.IUC3P()
sgextImg.from_pyarray(mask)

So:

sgextImg = sgext.itk.IUC3P()
sgextImg.from_pyarray(mask)
thin_image = sgext.scripts.thin(input=sgextImg, 
                   table_folder= sgext.tables_folder,
                   skel_type="end",
                   select_type="first",
                   persistence=2,
                   visualize=False,
                   verbose=True
                   )

If your mask has black foreground, first you need to convert foreground to white (255) and blackground to (0) for this case.

Readme has been updated.

Sorry, but there is still a bug:

my code is:

sgextImg = sgext.itk.IUC3P()
sgextImg.from_pyarray(mask*255)
thin_image = sgext.scripts.thin(input=sgextImg,
                   table_folder= sgext.tables_folder,
                   skel_type="end",
                   select_type="first",
                   persistence=2,
                   visualize=False,
                   verbose=True
                   )

And the bug is:

TypeError: thin(): incompatible function arguments. The following argument types are supported:
    1. (input: itk::SmartPointer<itk::Image<unsigned char,3> >, skel_type: str, select_type: str, tables_folder: str, persistence: int = 0, input_distance_map_image: itk::SmartPointer<itk::Image<float,3> > = None, profile: bool = False, verbose: bool = False, visualize: bool = False) -> itk::SmartPointer<itk::Image<unsigned char,3> >

Invoked with: kwargs: input=Dimension: 3
LargestPossibleRegion: 
  Index: [0, 0, 0]
  Size (i,j,k) (c_array): [440, 501, 501]
Origin: [0, 0, 0]
Spacing: [1, 1, 1]
Direction: 
1 0 0
0 1 0
0 0 1
Buffer: 
ImportImageContainer (000001B85A5370F0)
  RTTI typeinfo:   class itk::ImportImageContainer<unsigned __int64,unsigned char>
  Reference Count: 1
  Modified Time: 56
  Debug: Off
  Object Name: 
  Observers: 
    none
  Pointer: 000001B864783040
  Container manages memory: false
  Size: 110440440
  Capacity: 110440440
, table_folder='D:\\ProgramData\\Anaconda3\\envs\\general\\lib\\site-packages\\sgext\\tables', skel_type='end', select_type='first', persistence=2, visualize=False, verbose=True

There is a bug indeed, I will update sgext meanwhile, and it will solved for the next version, try with:
Change select_type="dmax" and create the distance map as follow:

sgextImg = sgext.itk.IUC3P()
sgextImg.from_pyarray(mask*255)
dmap_image = sgext.scripts.create_distance_map(input=sgextImg, verbose=True)
thin_image = sgext.scripts.thin(input=sgextImg,
                   table_folder= sgext.tables_folder,
                   skel_type="end",
                   select_type="dmax",
                   input_distance_map_image=dmap_image,
                   persistence=2,
                   visualize=False,
                   verbose=True
                   )
dmap_image = sgext.scripts.create_distance_map(input=sgextImg, verbose=True)

This code causes another bug:

unhandled win32 exception occurred in python.exe

My environment is:

win10
Python 3.7.5 (default, Oct 31 2019, 15:18:51) [MSC v.1916 64 bit (AMD64)] :: Anaconda, Inc. on win32
sgext: I don't know the version. I think it is the latest version.

That’s unfortunate, I am sorry. I will have a look. Is there a way you could get a more verbose error about that exception?
Issue created: https://github.com/phcerdan/SGEXT/issues/51

Meanwhile, you can use the following, use select_type="first" and provide an empty distance map image. This will be solved in the next release. Thanks for your time, and please report anything in the issues tab of SGEXT.

sgextImg = sgext.itk.IUC3P()
sgextImg.from_pyarray(mask*255)
dmap_image = sgext.itk.IF3P()
thin_image = sgext.scripts.thin(input=sgextImg,
                   table_folder= sgext.tables_folder,
                   skel_type="end",
                   select_type="first",
                   input_distance_map_image=dmap_image,
                   persistence=2,
                   visualize=False,
                   verbose=True
                   )