This notebook is intended to demonstrate how select registration, segmentation, and image mathematical methods of ITKTubeTK can be combined to perform multi-channel brain extraction (aka. skull stripping for patient data containing multiple MRI sequences).

There are many other (probably more effective) brain extraction methods available as open-source software such as BET and BET2 in the FSL package (albeit such methods are only for single channel data). If you need to perform brain extraction for a large collection of scans that do not contain major pathologies, please use one of those packages. This notebook is meant to show off the capabilities of specific ITKTubeTK methods, not to demonstration how to "solve" brain extraction.



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]

reader1 = ReaderType.New(FileName="data/mra.mha")
reader1.Update()
im1 = reader1.GetOutput()

reader2 = ReaderType.New(FileName="data/mri_t1_sag.mha")
reader2.Update()
im2 = reader2.GetOutput()

reader3 = ReaderType.New(FileName="data/mri_t2.mha")
reader3.Update()
im3 = reader3.GetOutput()



In [3]:

    
RegisterImagesType = ttk.RegisterImages[ImageType]

reg2To1 = RegisterImagesType.New(FixedImage=im1, MovingImage=im2)
reg2To1.SetReportProgress(True)
reg2To1.SetExpectedOffsetMagnitude(5)
reg2To1.SetExpectedRotationMagnitude(0.005)
reg2To1.SetRigidMaxIterations(100)
reg2To1.SetRegistration("PIPELINE_AFFINE")
reg2To1.SetMetric("MATTES_MI_METRIC")
reg2To1.Update()
im2Reg = reg2To1.ResampleImage()



In [4]:

    
reg3To1 = RegisterImagesType.New(FixedImage=im1, MovingImage=im3)
reg3To1.SetReportProgress(True)
reg3To1.SetExpectedOffsetMagnitude(5)
reg3To1.SetExpectedRotationMagnitude(0.005)
reg3To1.SetRigidMaxIterations(100)
reg3To1.SetRegistration("PIPELINE_AFFINE")
reg3To1.SetMetric("MATTES_MI_METRIC")
reg3To1.Update()
im3Reg = reg3To1.ResampleImage()



In [5]:

    
CheckerboardImagesType = itk.CheckerBoardImageFilter[ImageType]
cb23 = CheckerboardImagesType.New(Input1=im2Reg, Input2=im3Reg)
cb23.Update()
im23 = ImageType.New()
im23 = cb23.GetOutput()
view(im23)



In [6]:

    
N = 8
readerList = ["003", "010", "026", "034", "045", "056", "063", "071"]

imBase = []
imBaseB = []
for i in range(0,N):
    name = "data/Normal"+readerList[i]+"-FLASH.mha"
    nameB = "data/Normal"+readerList[i]+"-FLASH-Brain.mha"
    reader = ReaderType.New(FileName=name)
    reader.Update()
    imBaseTmp = reader.GetOutput()
    reader = ReaderType.New(FileName=nameB)
    reader.Update()
    imBaseBTmp = reader.GetOutput()
    imBase.append(imBaseTmp)
    imBaseB.append(imBaseBTmp)



In [7]:

    
regB = []
regBB = []
for i in range(0,N):
    regBTo1 = RegisterImagesType.New(FixedImage=im2Reg, MovingImage=imBase[i])
    regBTo1.SetReportProgress(True)
    regBTo1.SetExpectedOffsetMagnitude(5)
    regBTo1.SetExpectedRotationMagnitude(0.005)
    regBTo1.SetRigidMaxIterations(100)
    regBTo1.SetRegistration("PIPELINE_AFFINE")
    regBTo1.SetMetric("MEAN_SQUARED_ERROR_METRIC")
    regBTo1.Update()
    img = regBTo1.ResampleImage()
    regB.append( img )
    img = regBTo1.ResampleImage("LINEAR", imBaseB[i])
    regBB.append( img )



In [8]:

    
regBBT = []
for i in range(0,N):
    imMath = ttk.ImageMath[ImageType,ImageType].New( Input=regBB[i] )
    imMath.Threshold(0,1,0,1)
    img = imMath.GetOutput()
    if i==0:
        imMathSum = ttk.ImageMath[ImageType,ImageType].New( img )
        imMathSum.AddImages( img, 1.0/N, 0 )
        sumBBT = imMathSum.GetOutput()
    else:
        imMathSum = ttk.ImageMath[ImageType,ImageType].New( sumBBT )
        imMathSum.AddImages( img, 1, 1.0/N )
        sumBBT = imMathSum.GetOutput()
        
view(sumBBT)



In [9]:

    
insideMath = ttk.ImageMath[ImageType,ImageType].New( Input = sumBBT )
insideMath.Threshold(1,1,1,0)
insideMath.Erode(5,1,0)
brainInside = insideMath.GetOutput()

outsideMath = ttk.ImageMath[ImageType,ImageType].New( Input = sumBBT )
outsideMath.Threshold(0,0,1,0)
outsideMath.Erode(10,1,0)
brainOutsideAll = outsideMath.GetOutput()
outsideMath.Erode(20,1,0)
outsideMath.AddImages(brainOutsideAll, -1, 1)
brainOutside = outsideMath.GetOutput()

outsideMath.AddImages(brainInside,1,2)
brainCombinedMask = outsideMath.GetOutputUChar()

outsideMath.AddImages(im2Reg, 512, 1)
brainCombinedMaskView = outsideMath.GetOutput()
view(brainCombinedMaskView)



In [10]:

    
LabelMapType = itk.Image[itk.UC,3]

segmenter = ttk.SegmentConnectedComponentsUsingParzenPDFs[ImageType,LabelMapType].New()
segmenter.SetFeatureImage( im1 )
segmenter.AddFeatureImage( im2Reg )
segmenter.AddFeatureImage( im3Reg )
segmenter.SetInputLabelMap( brainCombinedMask )
segmenter.SetObjectId( 2 )
segmenter.AddObjectId( 1 )
segmenter.SetVoidId( 0 )
segmenter.SetErodeDilateRadius( 5 )
segmenter.Update()
segmenter.ClassifyImages()
brainCombinedMaskClassified = segmenter.GetOutputLabelMap()



In [11]:

    
cast = itk.CastImageFilter[LabelMapType, ImageType].New()
cast.SetInput(brainCombinedMaskClassified)
cast.Update()
brainMaskF = cast.GetOutput()

brainMath = ttk.ImageMath[ImageType,ImageType].New(Input = brainMaskF)
brainMath.Threshold(2,2,1,0)
brainMath.Dilate(1,1,0)
brainMaskD = brainMath.GetOutput()
brainMath.SetInput( im1 )
brainMath.ReplaceValuesOutsideMaskRange( brainMaskD, 1, 1, 0)
brain = brainMath.GetOutput()

view(brain)



In [ ]: