Train Vascular Model Notebook

This notebook contains an example of how to train a vascular model for automatic vessel segmentation from a CTA. For this example we are using CTA of the head. This example requires an expert mask as an input. Expert mask is a binary image volume, where vessels are marked 1 and rest is 0. It also requires one more mask, which serves as mask of the brain region with in the head CTA.

If ITK is not installed in your python environment, you need to define the environment variable ITK_BUILD_DIR that contains the path to where ITK was built.

We need to find the directory in which TubeTK was build. This is required to find the path to the testing data, and may be also required to find the TubeTK library paths if your python environment does not include it. The environment variable TubeTK_BUILD_DIR needs to be defined.



In [6]:

    
import os
import sys
import numpy



In [7]:

    
# Path for TubeTK libs and bin
#Values takend from TubeTK launcher

#sys.path.append("C:/src/TubeTK_Python_ITK/SlicerExecutionModel-build/GenerateCLP/")
#sys.path.append("C:/src/TubeTK_Python_ITK/SlicerExecutionModel-build/GenerateCLP/Release")

#sys.path.append("C:/src/TubeTK_Python_ITK/ITK-build/bin/")
#sys.path.append("C:/src/TubeTK_Python_ITK/ITK-build/bin/Release")

#sys.path.append("C:/src/TubeTK_Python_ITK/TubeTK-build/bin/")
#sys.path.append("C:/src/TubeTK_Python_ITK/TubeTK-build/bin/Release")
sys.path.append("C:/src/TubeTK_Python_ITK/TubeTK-build/lib/")
sys.path.append("C:/src/TubeTK_Python_ITK/TubeTK-build/lib/Release")

#sys.path.append("C:/src/TubeTK_Python_ITK/VTK-build/bin/")
#sys.path.append("C:/src/TubeTK_Python_ITK/VTK-build/bin/Release")



In [8]:

    
# Setting TubeTK Build Directory

TubeTK_BUILD_DIR=None
if 'TubeTK_BUILD_DIR' in os.environ:
    TubeTK_BUILD_DIR = os.environ['TubeTK_BUILD_DIR']
else:
    print('TubeTK_BUILD_DIR not found!')
    print('  Set environment variable')
    os.environ["TubeTK_BUILD_DIR"] = "C:/src/TubeTK_Python_ITK/TubeTK-build"
    TubeTK_BUILD_DIR = os.environ["TubeTK_BUILD_DIR"]
    #sys.exit( 1 )
    
if not os.path.exists(TubeTK_BUILD_DIR):
    print('TubeTK_BUILD_DIR set by directory not found!')
    print('  TubeTK_BUILD_DIR = ' + TubeTK_BUILD_DIR )
    sys.exit(1)



In [9]:

    
# Setting ITK Build Directory and importing ITK

try:
    import itk
except:
    ITK_BUILD_DIR = None
    if 'ITK_BUILD_DIR' in os.environ:
        ITK_BUILD_DIR = os.environ['ITK_BUILD_DIR']
    else:
        print('ITK_BUILD_DIR not found!')
        print('  Set environment variable')
        os.environ["ITK_BUILD_DIR"] = "C:/src/TubeTK_Python_R/ITK-build"
        ITK_BUILD_DIR = os.environ["ITK_BUILD_DIR"]
        #sys.exit( 1 )

    if not os.path.exists(ITK_BUILD_DIR):
        print('ITK_BUILD_DIR set by directory not found!')
        print('  ITK_BUIDL_DIR = ' + ITK_BUILD_DIR )
        sys.exit(1)
    # Append ITK libs
    sys.path.append("C:/src/TubeTK_Python_ITK/ITK-build/Wrapping/Generators/Python/Release")
    sys.path.append("C:/src/TubeTK_Python_ITK/ITK-build/lib/Release")
    sys.path.append("C:/src/TubeTK_Python_ITK/ITK-build/lib")
    
    # Append TubeTK libs
    sys.path.append("C:/src/TubeTK_Python_ITK/TubeTK-build/ITKModules/TubeTKITK-build/Wrapping/Generators/Python/Release")
    import itk



In [10]:

    
from itk import TubeTKITK as itktube

Initialization



In [13]:

    
Dimension = 3
PixelType = itk.F

CTImageFileName = os.path.join(TubeTK_BUILD_DIR, 'MIDAS_Data\inputCTA.mha')
ExpertMaskImageFileName = os.path.join(TubeTK_BUILD_DIR, 'MIDAS_Data\inputExpertMask.mha')
MaskImageFileName = os.path.join(TubeTK_BUILD_DIR, 'MIDAS_Data\inputMask.mha')


SpatialObjectType = itk.SpatialObject[Dimension]

Read the input images



In [14]:

    
ImageType = itk.Image[PixelType, Dimension]
ImageReaderType = itk.ImageFileReader[ImageType]

imageReader1 = ImageReaderType.New()
imageReader1.SetFileName(CTImageFileName)
imageReader1.Update()

CTImage = imageReader1.GetOutput()

imageReader2 = ImageReaderType.New()
imageReader2.SetFileName(ExpertMaskImageFileName)
imageReader2.Update()

ExpertMaskImage = imageReader2.GetOutput()

imageReader3 = ImageReaderType.New()
imageReader3.SetFileName(MaskImageFileName)
imageReader3.Update()

MaskImage = imageReader3.GetOutput()









    



itkImageFileReaderIF3: 0.000000itkImageFileReaderIF3: 1.000000itkImageFileReaderIF3: 0.000000itkImageFileReaderIF3: 1.000000itkImageFileReaderIF3: 0.000000itkImageFileReaderIF3: 1.000000

STEP 1: Crop the input volumes to make them of same size as MaskImage



In [27]:

    
boundary = itk.Index[3]()
boundary.Fill(10)

#Create the crop image filter
CropImageFilterType = itktube.CropImage[ImageType, ImageType]
cropImageFilter = CropImageFilterType.New()
cropImageFilter.SetBoundary(boundary)
#cropImageFilter.SetMatchVolume(MaskImage)  #Giving error

#Crop Input CTA
cropImageFilter.SetInput(CTImage)
cropImageFilter.Update()

croppedCTImage = cropImageFilter.GetOutput()

#Crop Expert Mask
cropImageFilter.SetInput(ExpertMaskImage)
cropImageFilter.Update()

croppedExpertMaskImage = cropImageFilter.GetOutput()

#Crop Mask
cropImageFilter.SetInput(MaskImage)
cropImageFilter.Update()

croppedMaskImage = cropImageFilter.GetOutput()









    



---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-27-6049da1a2d19> in <module>()
     31 imageWriter.SetInput(croppedMaskImage)
     32 imageWriter.SetFileName(MaskImageFileName)
---> 33 imageWriter.Update()

RuntimeError: c:\src\tubetk_python_itk\itk\modules\io\imagebase\include\itkImageFileWriter.hxx:290:
itk::ERROR: ImageFileWriter(000000001700C090): Largest possible region does not fully contain requested paste IO regionPaste IO region: ImageIORegion (000000000021CA00)
  Dimension: 3
  Index: 0 0 0 
  Size: 0 0 0 
Largest possible region: ImageRegion (000000000021CAC8)
  Dimension: 3
  Index: [0, 0, 0]
  Size: [0, 0, 0]

STEP 2: Resample the cropped images.



In [ ]:

    
boundary = itk.Index[3]()
boundary.Fill(10)

#Create the resample image filter
ResampleImageFilterType = itktube.ResampleImage[ImageType, ImageType]

#Resample Input CTA
resampleImageFilter1 = ResampleImageFilterType.New()
resampleImageFilter1.SetInput(croppedCTImage)
resampleImageFilter1.SetMakeIsotropic(True)
resampleImageFilter1.SetInterpolator("Sinc")

resampleCTImage = resampleImageFilter1.GetOutput()

#Resample Expert Mask
resampleImageFilter2 = ResampleImageFilterType.New()
resampleImageFilter2.SetInput(croppedExpertMaskImage)
resampleImageFilter2.SetMakeIsotropic(True)
resampleImageFilter2.SetInterpolator("NearestNeighbor")

resampleExpertMaskImage = resampleImageFilter2.GetOutput()

#Resample Mask
resampleImageFilter3 = ResampleImageFilterType.New()
resampleImageFilter3.SetInput(croppedMaskImage)
resampleImageFilter3.SetMakeIsotropic(True)
resampleImageFilter3.SetInterpolator("NearestNeighbor")

resampleMaskImage = resampleImageFilter3.GetOutput()

STEP 3: Create Mask-only images. this step required Image Math



In [ ]:

    
# create resampleMaskImage -> erodedResampleMaskImage
# resampleExpertMaskImage -> erodedResampleExpertMaskImage
# resampleCTImage -> maskedResampleCTImage

STEP 4: Compute Training Mask



In [ ]:

    
# Create image to save not-vessel mask.
ShortImageType = itk.Image[itk.S, Dimension]
notVesselMaskImage = ShortImageType.New()

#Create Compute Training Mask Filter
ComputeTrainingMaskFilterType = itktube.ComputeTrainingMask[ImageType]

computeTrainingMaskFilter = ComputeTrainingMaskFilterType.New()
computeTrainingMaskFilter.SetInput(erodedResampleExpertMaskImage)
computeTrainingMaskFilter.SetNotVesselMask(notVesselMaskImage)
computeTrainingMaskFilter.SetGap(0.5)
computeTrainingMaskFilter.SetNotVesselWidth(2)
computeTrainingMaskFilter.Update()

expertTrainMaskImage = computeTrainingMaskFilter.GetOutput()

STEP 5: Enhance Vessels in maskedResampleCTImage



In [ ]:

    
DiscriminantInfoFileName = os.path.join(TubeTK_BUILD_DIR, 'Temporary\\vascularModel.mrs')
enhancedScalesExpertMaskImage = ImageType.New()

# Create EnhanceTubesUsingDiscriminantAnalysis Filter
EnhanceTubesUsingDiscriminantAnalysisFilterType = itktube.EnhanceTubesUsingDiscriminantAnalysis[ImageType, ImageType]

ETUDAFilter = EnhanceTubesUsingDiscriminantAnalysisFilterType.New()
ETUDAFilter.SetInput(maskedResampleCTImage)
ETUDAFilter.SetLabelMap(expertTrainMaskImage)
ETUDAFilter.SetTubeId(255)
ETUDAFilter.SetBackgroundId(127)
ETUDAFilter.SetSaveDiscriminantInfo(DiscriminantInfoFileName)
ETUDAFilter.SetOutputSeedScaleImage(enhancedScalesExpertMaskImage)
ETUDAFilter.SetTubeScales(0.4,0.8,1.2,1.6)

enhancedExpertMaskImage = ETUDAFilter.GetOutput()

STEP 6: Compute Segment tubes Parameters



In [ ]:

    
vasculaModelParameterFileName = os.path.join(TubeTK_BUILD_DIR, 'Temporary\\vascularModel.mtp')

# Create SegmentTubesParameters Filter
ComputeSegmentTubesParametersFilterType = itktube.ComputeSegmentTubesParameters[ImageType]

CSTPFilter = ComputeSegmentTubesParametersFilterType.New()
CSTPFilter.SetInput(maskedResampleCTImage)
CSTPFilter.SetMaskImage(expertTrainMaskImage)
CSTPFilter.SetScaleImage(enhancedScalesExpertMaskImage)
CSTPFilter.SetParametersFileName(vasculaModelParameterFileName)
CSTPFilter.Update()

STEP 7: This step requires Image Math



In [ ]:

    
# enhancedExpertMaskImage -> vesselEnhancedExpertMaskImage

STEP 8: Compute Seeds



In [ ]:

    
# Create SegmentBinaryImageSkeleton Filter
SegmentBinaryImageSkeletonFilterType = itktube.SegmentBinaryImageSkeleton[Imagetype]

SBISFilter = SegmentBinaryImageSkeletonFilterType.New()
SBISFilter.SetInput(vesselEnhancedExpertMaskImage)
SBISFilter.Update()

seedsVesselEnhancedExpertMaskImage = SBISFilter.GetOutput()

STEP 9: Segment Tubes



In [ ]:

    
outputVesselsFileName = os.path.join(TubeTK_BUILD_DIR, 'Temporary\\outputVessels.tre')



# Create SegmentTubes Filter
SegmentTubesFilterType = itktube.SegmentTubes[ImageType]

SegmenttubesFilter = SegmentTubesFilterType.New()
SegmenttubesFilter.SetInput(maskedResampleCTImage)