Vmtk Vesselness Filtering without seed point?

Tommaso_Di_Noto · September 17, 2019, 12:45pm

Hi everyone!

I was studying tutorials about vmtk and I saw the Vesselness Filtering module. I find it really useful for a single patient, and I manged to make it work by manually specifying the initial seed (fiducial) point.

But what if I have 300 cases? Is there a way to perform a batch filtering (ideally from command line) WITHOUT the manual seed (even less accurate, but automatic)?

Thanks a lot in advance!

lassoan · September 17, 2019, 1:36pm

Seed point is used only to automatically compute min/max diameter and contrast setting. If you work with similar images then you can determine on a few images what values tend to work well and use those in all images.

Tommaso_Di_Noto · September 17, 2019, 1:41pm

Thanks a lot @lassoan!
So to my next question…how can I run the tool from command line so that I can then encapsulate everything as subprocess in a python script looping over all my cases?

In other words, how can I access this slicer module (and in general any module) from command line, knowing the args, etc?

Thanks again!

lassoan · September 17, 2019, 1:48pm

It is a Python scripted module, so instead of old-school and platform-dependent bash/batch scripting, you can write a short Python script or Jupyter notebook to iterate through your data sets.

You can run vessel filtering by calling computeVesselnessVolume method of the VesselnessFiltering module logic:

import VesselnessFiltering
vfl=VesselnessFiltering.VesselnessFilteringLogic()
vfl.computeVesselnessVolume(inputVolumeNode, outputVolumeNode, minimumDiameterMm=0.2, maximumDiameterMm=25, alpha=0.3, beta=0.3, contrastMeasure=150)

Tommaso_Di_Noto · September 17, 2019, 1:52pm

Ok great, that sounds even easier. I’ll give it a try and let you know how it goes.

Thank you very much again!

Tommaso_Di_Noto · September 17, 2019, 2:43pm

So, I’ve opened a New->Slicer4.9 notebook using Binder but when I import VesselnessFiltering, I get:

Traceback (most recent call last):
  File "<string>", line 1, in <module>
ImportError: No module named VesselnessFiltering

I tried running the example 01_Data_Loading_and_Visualization and that works.

lassoan · September 17, 2019, 4:06pm

VMTK extension is not installed in Slicer version that runs in Binder. You could install it using a Python script, but for batch processing, it would not be convenient to upload/download data anyway, so I would recommend to run the processing script on your local computer.

Tommaso_Di_Noto · September 18, 2019, 3:10pm

Dear @lassoan, going through a local script on jupyter worked.

So, I’ve been trying the two following snippets and none of them work.

Option 1:

import numpy as np
import VesselnessFiltering
import SimpleITK as sitk
import sitkUtils as su

one_patient_path = "/home/newuser/.../sub-204_ses-20110825_angio_brain.nii.gz"
InputVolumeNode = slicer.util.loadVolume(one_patient_path, returnNode=True)

# Create output volume
OutputVolume = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode")  # add new node
OutputVolume.SetName("output_volume_trial")  # Set name for output volume
OutputVolume.CreateDefaultDisplayNodes()

vfl = VesselnessFiltering.VesselnessFilteringLogic()
vfl.computeVesselnessVolume(InputVolumeNode, OutputVolume, minimumDiameterMm=0, maximumDiameterMm=25, alpha=0.3, beta=0.3, contrastMeasure=150)

and this gives me the following error:

File "/home/newuser/.config/NA-MIC/Extensions-27931/SlicerVMTK/lib/Slicer-4.10/qt-scripted- 
modules/VesselnessFiltering.py", line 492, in computeVesselnessVolume
inImage.DeepCopy(currentVolumeNode.GetImageData())
AttributeError: 'tuple' object has no attribute 'GetImageData'

Option 2:

# same imports...

# go through sitk
input_volume = sitk.ReadImage(one_patient_path)  # load nifti volume as sitk image
InputVolumeNode = su.PushVolumeToSlicer(input_volume, None)

# Create output volume
OutputVolume = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode")  # add new node
OutputVolume.SetName("output_volume_trial")  # Set name for output volume
OutputVolume.CreateDefaultDisplayNodes()

and in this case the last cell below never ends and the output volume is not created.

vfl = VesselnessFiltering.VesselnessFilteringLogic()
vfl.computeVesselnessVolume(InputVolumeNode, OutputVolume, minimumDiameterMm=0, maximumDiameterMm=25, alpha=0.3, beta=0.3, contrastMeasure=150)

What am I doing wrong?

Thank you very much again!

lassoan · September 18, 2019, 4:27pm

The problem is what the error message tells: currentVolumeNode is a Python tuple. It is not a MRML node. Maybe you used slicer.util.loadVolume in Slicer-4.10 and earlier, which returns a tuple of (success, node) instead of just the node. Check the API documentation.

Tommaso_Di_Noto:

and in this case the last cell below never ends and the output volume is not created.

vfl = VesselnessFiltering.VesselnessFilteringLogic()
vfl.computeVesselnessVolume(InputVolumeNode, OutputVolume, minimumDiameterMm=0, maximumDiameterMm=25, alpha=0.3, beta=0.3, contrastMeasure=150)

Does this work with the same input and output volumes if you use the vesselness filtering module GUI? Is there any error displayed in the console? You can check the application log in Slicer or type the commands in the Python console in Slicer. If Slicer hangs then you can restart it and find logs of the previous sessions in menu: Help / Report a bug.

Tommaso_Di_Noto · September 19, 2019, 8:23am

@lassoan, apologies for not reading the API doc carefully. I now modified as:

_, InputVolumeNode = slicer.util.loadVolume(one_patient_path, returnNode=True)
dim = InputVolumeNode.GetImageData().GetDimensions()
spacing = InputVolumeNode.GetSpacing()

and following this thread, I created the Output Node with:

def createEmptyVolume(imageSize, imageSpacing, nodeName):
    voxelType=vtk.VTK_FLOAT
    # Create an empty image volume
    imageData=vtk.vtkImageData()
    imageData.SetDimensions(imageSize)
    imageData.AllocateScalars(voxelType, 1)
    thresholder=vtk.vtkImageThreshold()
    thresholder.SetInputData(imageData)
    thresholder.SetInValue(0)
    thresholder.SetOutValue(0)
    thresholder.Update()
    # Create volume node
    volumeNode=slicer.mrmlScene.AddNewNodeByClass("vtkMRMLScalarVolumeNode", nodeName)
    volumeNode.SetSpacing(imageSpacing)
    volumeNode.SetAndObserveImageData(thresholder.GetOutput())
    volumeNode.CreateDefaultDisplayNodes()
    volumeNode.CreateDefaultStorageNode()
    return volumeNode  

# creat output node with same dim and spacing of InputNode
OutputNode = createEmptyVolume(dim, spacing, 'VesselnessFiltered')  # invoke function

Now, if I try to apply the VesselNessFiltering from the GUI with these two volumes (and with the manual seed point), the “VesselnessFiltered” volume is created correctly.

However, when I run in Jupyter:

vfl = VesselnessFiltering.VesselnessFilteringLogic()
vfl.computeVesselnessVolume(InputVolumeNode, OutputVolume, minimumDiameterMm=0, maximumDiameterMm=25, alpha=0.3, beta=0.3, contrastMeasure=150)

I get the following log in Slicer (from Help/Report a bug):

[DEBUG][Python] 19.09.2019 10:04:22 [Python] (/home/newuser/.config/NA-MIC/Extensions-27931
SlicerVMTK/lib/Slicer-4.10/qt-scripted-modules/VesselnessFiltering.py:465) - Vesselness filtering 
started: diameter min=0, max=25, alpha=0.3, beta=0.3, contrastMeasure=150
[CRITICAL][Stream] 19.09.2019 10:04:23 [] (unknown:0) - ERROR: during execute_request
[CRITICAL][Stream] 19.09.2019 10:04:23 [] (unknown:0) - /work/Stable/Slicer-4101- 
build/ITK/Modules/Filtering/Smoothing/include/itkRecursiveGaussianImageFilter.hxx:344:
[CRITICAL][Stream] 19.09.2019 10:04:23 [] (unknown:0) - itk::ERROR: 
RecursiveGaussianImageFilter(0xb7ded90): Sigma must be greater than zero.

The [DEBUG] tells us that the VesselnessFiltering started, but then the itkRecursiveGaussianImageFilter raises the CRITICAL error about the Sigma. What I don’t understand is where, inside computeVesselnessVolume, this itk function is invoked.

Pardon me if the question is trivial, but how could I find this out? Or, in other words, how can I debug in cases like this to understand what calls itkRecursiveGaussianImageFilter?

Thanks again for your time!

lassoan · September 19, 2019, 12:12pm

This error message suggests that the minimum vessel diameter must be greater than zero. Default 0.0 value in computeVesselnessVolume method definition is misleading, it should be probably set to 0.2 (this value is used in the GUI).

Tommaso_Di_Noto · September 23, 2019, 12:11pm

Dear @lassoan, thanks again for the tip!! It worked.

Finally, following several threads and examples, I think I managed to obtain what I wanted. Below, a code snippet for whomever might encounter a similar task (i.e. batch brain vessel extraction from MRA using VMTK and Jupyter Notebook). As the code is now, everything works, but please let me know in case you notice any imprecisions or avoidable passages (e.g. I’m not quite sure that the passage LabelNode–>SegNode–>LabelNodeAgain is very efficient):

# data_path is the folder where I have all my subfolders containing the MRA volumes

for subdir, dirs, files in os.walk(data_path):  # loop over folders and sub-folders
    for file in files:

        # load input volume
        input_volume_path = os.path.join(subdir,file)
        _, InputVolumeNode = slicer.util.loadVolume(input_volume_path, returnNode=True)  # load input volume
        dim_in = InputVolumeNode.GetImageData().GetDimensions()  # save input volume shape
        spacing_in = InputVolumeNode.GetSpacing()  # save input volume spacing
        
        # Create output volume with same dim and spacing of InputVolume
        OutputNode = createEmptyVolume(dim_in, spacing_in, 'VesselnessFiltered')  # invoke function
        
        # create vesselness filtering logic
        vfl = VesselnessFiltering.VesselnessFilteringLogic()  
        # compute VesselnessFiltering
        vfl.computeVesselnessVolume(InputVolumeNode, OutputNode, minimumDiameterMm=0.5, maximumDiameterMm=30, alpha=0.3, beta=0.3, contrastMeasure=130)
        
        # Create temporary labelmap from output volume
        labelVolumeNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLLabelMapVolumeNode")
        slicer.vtkSlicerVolumesLogic().CreateLabelVolumeFromVolume(slicer.mrmlScene, labelVolumeNode, OutputNode)
        
        # Fill temporary labelmap by thresholding
        thresholdValue = 0.08  # set a threshold value
        voxelArray = slicer.util.arrayFromVolume(OutputNode)  # extract numpy array from output volume
        labelVoxelArray = slicer.util.arrayFromVolume(labelVolumeNode)  # extract numpy array from labelmap volume
        labelVoxelArray[voxelArray < thresholdValue] = 0 # set all voxels below thr to 0
        labelVoxelArray[voxelArray >= thresholdValue] = 1 # set all voxels above thr to 1
        
        # Import labelmap to segmentation
        segmentationNode = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLSegmentationNode')  # create segmentation node
        slicer.modules.segmentations.logic().ImportLabelmapToSegmentationNode(labelVolumeNode, segmentationNode)
        segmentationNode.CreateClosedSurfaceRepresentation()  # make segmentation visible in 3D
        
        # Create segment editor to get access to effects
        segmentEditorWidget = slicer.qMRMLSegmentEditorWidget()
        segmentEditorWidget.setMRMLScene(slicer.mrmlScene)
        segmentEditorNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLSegmentEditorNode")
        segmentEditorWidget.setMRMLSegmentEditorNode(segmentEditorNode)
        segmentEditorWidget.setSegmentationNode(segmentationNode)
        segmentEditorWidget.setMasterVolumeNode(OutputNode)
        
        # Remove small islands from segmentation
        segmentEditorWidget.setActiveEffectByName("Islands")
        effect = segmentEditorWidget.activeEffect()
        effect.setParameter('Operation','REMOVE_SMALL_ISLANDS')  # remove islands smaller than a certain dimension
        effect.setParameter('MinimumSize', 170)
        effect.self().onApply()
        
        # Clean up
        segmentEditorWidget = None
        slicer.mrmlScene.RemoveNode(segmentEditorNode)  # remove segmentation editor node
        slicer.mrmlScene.RemoveNode(labelVolumeNode)  # remove temporary LabelMapVolume but keep segmentation node
        
        # Convert again from segmentation node to labelnode matching input volume dimensions
        reference = getNode("VesselnessFiltered") # this will be the volume the segmentation was drawn on
        labelmapVolumeNodeClosed = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLLabelMapVolumeNode')
        seg = getNode("Segmentation")  # set segmentation node            
        ids = vtk.vtkStringArray()
        seg.GetDisplayNode().GetVisibleSegmentIDs(ids)
        smsl = slicer.modules.segmentations.logic()
        smsl.ExportSegmentsToLabelmapNode(segmentationNode,ids,labelmapVolumeNodeClosed,reference)
        
        dim_out_lab = labelmapVolumeNodeClosed.GetImageData().GetDimensions()  # save output shape
        spacing_out_lab = labelmapVolumeNodeClosed.GetSpacing()  # save output spacing
        
        # make sure that shape and spacing match from input to final output nodes
        assert(dim_in == dim_out_lab)
        assert(spacing_in == spacing_out_lab)
        
        # save final node
        file_name = os.path.basename(file)
        index_of_dot = file_name.index('.')
        file_name_without_extensions = file_name[:index_of_dot]  # remove extensions
        output_file_path = os.path.join(subdir, file_name_without_extensions) + "_vessels.nii.gz"
        
        slicer.util.saveNode(labelmapVolumeNodeClosed, output_file_path)  # save output node
        print("Vessel mask saved for {0}".format(file))
        
        # clear Slicer scene (it removes all volumes)
        slicer.mrmlScene.Clear(0)

P.S. if you modify the minimum diameter to 0.2 in your first answer, I think we can (I can) mark it as “Solution” answer

Tommaso_Di_Noto · February 20, 2020, 4:18pm

Dear @lassoan,

Is there, to your knowledge, a way to obtain a smoothed (i.e. probabilistic) vesselness filtering rather than a binary one as in this post, either in VMTK or with another software?

Thanks a lot in advance!

lassoan · February 20, 2020, 6:58pm

VMTK’s vesselness filter computes a probablilistic image with values is between 0.0 and 1.0. By default we display the computed vesselness image by applying a threshold value in display options, so it may appear binary, but it is actually not.

Tommaso_Di_Noto · February 21, 2020, 10:53am

Oh, you are right!

The OutputNode in the example is indeed probabilistic after applying vfl.computeVesselnessVolume, so I can just save that.

Thanks a lot, as usual, for the super quick answer

Topic		Replies	Views
Vmtk Vesselness Filtering: select seedpoint in Python Support segmentation , python , vmtk , vessels	1	582	July 3, 2022
Questions about the vesselness filtering function Development vmtk	0	313	August 9, 2022
SlicerExtension-VMTK - VesselnessFiltering - Code Development segmentation , extensions-manager , python	2	437	November 26, 2019
VMTK toolkit is too different to the available tutorial Support	8	720	March 28, 2019
Vessel Segmentation using VMTK in Slicer Support	1	387	September 17, 2019

Vmtk Vesselness Filtering without seed point?

Related topics