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()


---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
<ipython-input-2-5e9f3d8180df> in <module>()
      1 ImageType = itk.Image[itk.F, 3]
      2 ReaderType = itk.ImageFileReader[ImageType]
----> 3 ResampleType = ttk.ResampleImage[ImageType]
      4 
      5 reader1 = ReaderType.New(FileName="data/mra-Brain.mha")

/home/stephen.aylward/src/ITK-Release/Wrapping/Generators/Python/itkLazy.py in __getattribute__(self, attr)
     50             module = self.__belong_lazy_attributes[attr]
     51             namespace = {}
---> 52             itkBase.LoadModule(module, namespace)
     53             self.loaded_lazy_modules.add(module)
     54             for k, v in namespace.items():

/home/stephen.aylward/src/ITK-Release/Wrapping/Generators/Python/itkBase.py in LoadModule(name, namespace)
     91     loader = LibraryLoader()
     92     if not swigModuleName in sys.modules:
---> 93         module = loader.load(swigModuleName)
     94 
     95     # OK, now the modules on which this one depends are loaded and

/home/stephen.aylward/src/ITK-Release/Wrapping/Generators/Python/itkBase.py in load(self, name)
    209         self.setup()
    210         try:
--> 211             return importlib.import_module(name)
    212         finally:
    213             self.cleanup()

/usr/lib/python3.6/importlib/__init__.py in import_module(name, package)
    124                 break
    125             level += 1
--> 126     return _bootstrap._gcd_import(name[level:], package, level)
    127 
    128 

/usr/lib/python3.6/importlib/_bootstrap.py in _gcd_import(name, package, level)

/usr/lib/python3.6/importlib/_bootstrap.py in _find_and_load(name, import_)

/usr/lib/python3.6/importlib/_bootstrap.py in _find_and_load_unlocked(name, import_)

/usr/lib/python3.6/importlib/_bootstrap.py in _load_unlocked(spec)

/usr/lib/python3.6/importlib/_bootstrap_external.py in exec_module(self, module)

/usr/lib/python3.6/importlib/_bootstrap.py in _call_with_frames_removed(f, *args, **kwds)

/home/stephen.aylward/src/ITK-Release/lib/TubeTKPython.py in <module>()
     30                 fp.close()
     31             return _mod
---> 32     _TubeTKPython = swig_import_helper()
     33     del swig_import_helper
     34 else:

/home/stephen.aylward/src/ITK-Release/lib/TubeTKPython.py in swig_import_helper()
     26         if fp is not None:
     27             try:
---> 28                 _mod = imp.load_module('_TubeTKPython', fp, pathname, description)
     29             finally:
     30                 fp.close()

/usr/lib/python3.6/imp.py in load_module(name, file, filename, details)
    241                 return load_dynamic(name, filename, opened_file)
    242         else:
--> 243             return load_dynamic(name, filename, file)
    244     elif type_ == PKG_DIRECTORY:
    245         return load_package(name, filename)

/usr/lib/python3.6/imp.py in load_dynamic(name, path, file)
    341         spec = importlib.machinery.ModuleSpec(
    342             name=name, loader=loader, origin=path)
--> 343         return _load(spec)
    344 
    345 else:

/usr/lib/python3.6/importlib/_bootstrap.py in _load(spec)

/usr/lib/python3.6/importlib/_bootstrap.py in _load_unlocked(spec)

/usr/lib/python3.6/importlib/_bootstrap.py in module_from_spec(spec)

/usr/lib/python3.6/importlib/_bootstrap_external.py in create_module(self, spec)

/usr/lib/python3.6/importlib/_bootstrap.py in _call_with_frames_removed(f, *args, **kwds)

ImportError: /home/stephen.aylward/src/ITK-Release/lib/_TubeTKPython.so: undefined symbol: _ZN2af5arrayD1Ev

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 [ ]: