This notebook is intended to demonstrate how vessel segmentation methods of ITKTubeTK can be applied to multi-channel MRI (MRA + T1, T2, etc).
In [1]:
import itk
from itk import TubeTK as ttk
from itkwidgets import view
import numpy as np
In [2]:
ImageType = itk.Image[itk.F, 3]
ReaderType = itk.ImageFileReader[ImageType]
ResampleType = ttk.ResampleImage[ImageType]
reader1 = ReaderType.New(FileName="data/mra-Brain.mha")
reader1.Update()
im1 = reader1.GetOutput()
res = ResampleType.New(Input = im1) # Vessel extraction requires Isotropic voxel spacing
res.SetMakeHighResIso(True)
res.Update()
im1iso = res.GetOutput()
reader2 = ReaderType.New(FileName="data/mri-t1-Brain.mha")
reader2.Update()
im2 = reader2.GetOutput()
res = ResampleType.New(Input = im2)
res.SetMakeHighResIso(True)
res.Update()
im2iso = res.GetOutput()
reader3 = ReaderType.New(FileName="data/mri-t2-Brain.mha")
reader3.Update()
im3 = reader3.GetOutput()
res = ResampleType.New(Input = im3)
res.SetMakeHighResIso(True)
res.Update()
im3iso = res.GetOutput()
In [ ]:
# Manually extract a few vessels to form an image-specific training set
xp=[[-3.9,-54.8,12.3],
[25.1,-20.8,-22.7],
[-27.0,-25.8,-23.9]]
xs=[0.333,2,2]
vSeg = ttk.SegmentTubes[ImageType].New()
vSeg.SetInput(im1iso)
vSeg.SetVerbose(True)
vSeg.SetMinRoundness(0.1)
vSeg.SetMinCurvature(0.05)
vSeg.SetRadiusInObjectSpace( xs[0] )
vSeg.ExtractTubeInObjectSpace( xp[0], 1 )
vSeg.SetRadiusInObjectSpace( xs[1] )
vSeg.ExtractTubeInObjectSpace( xp[1], 2 )
vSeg.SetRadiusInObjectSpace( xs[2] )
vSeg.ExtractTubeInObjectSpace( xp[2], 3 )
tubeMaskImage = vSeg.GetTubeMaskImage()
In [ ]:
view(tubeMaskImage)
In [5]:
LabelMapType = itk.Image[itk.UC,3]
trMask = ttk.ComputeTrainingMask[ImageType,LabelMapType].New()
trMask.SetInput( tubeMaskImage )
trMask.SetGap( 3 )
trMask.SetNotObjectWidth( 1 )
trMask.Update()
fgMask = trMask.GetOutput()
In [6]:
im1Math = ttk.ImageMath[ImageType,LabelMapType].New(Input=im1iso)
im1Math.Threshold(0,0,63,0)
bkgMask = im1Math.GetOutput()
imCombMath = ttk.ImageMath[LabelMapType,LabelMapType].New(Input=bkgMask)
imCombMath.AddImages(fgMask, 1, 1)
mask = imCombMath.GetOutput()
In [7]:
imWriter = itk.ImageFileWriter[LabelMapType].New(Input=fgMask)
imWriter.SetFileName("data/im1iso-trainMask.mha")
imWriter.Update()
In [13]:
enhancer = ttk.EnhanceTubesUsingDiscriminantAnalysis[ImageType,LabelMapType].New()
enhancer.SetInput( im1iso )
enhancer.AddInput( im2iso )
enhancer.AddInput( im3iso )
enhancer.SetLabelMap( fgMask )
enhancer.SetRidgeId( 255 )
enhancer.SetBackgroundId( 127 )
enhancer.SetUnknownId( 0 )
enhancer.SetTrainClassifier(True)
enhancer.SetUseIntensityOnly(True)
enhancer.SetScales([0.3333,1,2.5])
enhancer.Update()
enhancer.ClassifyImages()
In [16]:
#view(enhancer.GetOutput())
view(enhancer.GetClassProbabilityImage(0))
#imMath = ttk.ImageMath[ImageType,ImageType].New(Input = segmenter.GetClassProbabilityImage(0))
#imMath.AddImages( segmenter.GetClassProbabilityImage(1), 1, -1 )
#view(imMath.GetOutput())
In [18]:
writer = itk.ImageFileWriter[ImageType].New(Input=enhancer.GetClassProbabilityImage(0))
writer.SetFileName("Class0Prob.mha")
writer.Update()
In [ ]: