This test runs a two-stages registration process including following stages:
1) Rigid
2) Affine
The main difference between this test and the test 1 is the way that registration process is initialized and the different stages are connected to eachother. Here we use SetMovingInitialTransform() method to initialize the second stage using the output of the firts stage. Also, we use Inplace=False, so at each stage the output transform object is different than the input transform object. We use a composite transform as container to concatenate the initial transform and the outputs of different stages together.
Also, here we use fixed and moving images masks to do the registration only on regions of interest.
We evaluate the results of the current registration process based on the results of BRIANSFit program that is run with a consistent command set. BRIANSFit is a 3D registration tool included in BRAINSTools. Get BRAINSTools from: https://github.com/BRAINSia/BRAINSTools
In [1]:
## Boiler plate code common to many notebooks. See the TestFilesCommonCode.ipynb for details
from __future__ import print_function
%run TestFilesCommonCode.ipynb
In [2]:
######## WARNING:
# You will need to replace the following path with your local BRAINSFit path.
BRIANSToolsPath="/Users/aghayoor/WorkSpace/BS/release/bin"
BRAINSFitPath=os.path.join(BRIANSToolsPath,"BRAINSFit")
if not os.path.isfile(BRAINSFitPath):
print("ERROR: This script requires that BRAINSFit can be found")
print("{0} does not exists".format(BRAINSFitPath))
In [3]:
fixedImageFilename = fdata('BRAINSTools/test.nii.gz')
movingImageFilename = fdata('BRAINSTools/rotation.test.nii.gz')
fixed = sitk.ReadImage(fixedImageFilename,sitk.sitkFloat32)
moving = sitk.ReadImage(movingImageFilename, sitk.sitkFloat32)
#referenceBaselineImageFilename = fdata('BRAINSTools/BRAINSFitTest_AffineRotationMasks.result.nii.gz')
#baseline = sitk.ReadImage(referenceBaselineImageFilename, sitk.sitkFloat32)
fixedMaskFilename = fdata("BRAINSTools/test_mask.nii.gz")
movingMaskFilename = fdata("BRAINSTools/rotation.test_mask.nii.gz")
fixedMask = sitk.ReadImage(fixedMaskFilename, sitk.sitkFloat32)
movingMask = sitk.ReadImage(movingMaskFilename, sitk.sitkFloat32)
myshow3d(fixed, xslices=[], yslices=[], zslices=[], title="fixed")
myshow3d(moving, xslices=[], yslices=[], zslices=[], title="moving")
#myshow3d(baseline, xslices=[], yslices=[], zslices=[], title="baseline")
myshow3d(fixedMask, xslices=[], yslices=[], zslices=[], title="fixed mask")
myshow3d(movingMask, xslices=[], yslices=[], zslices=[], title="moving mask")
In [4]:
cshow3d(fixed,moving)
In [5]:
# Read the initial transform
initialMovingTxFilename = fdata('BRAINSTools/BRAINSFitTest_Initializer_RigidRotationNoMasks.h5')
initialMovingTx = sitk.ReadTransform(initialMovingTxFilename)
print("Initial moving transform: ")
print(initialMovingTx)
# Define a composite transform that includes the initial transform.
# This composite tranform has the container role.
ctx = sitk.Transform(initialMovingTx)
In [6]:
BFitOut=os.path.realpath('Data/BFitOut.nii.gz')
# RUN THE BRAINSFIT FIRST ######################
!{BRIANSToolsPath}/BRAINSFit \
--costMetric MMI \
--numberOfIterations 250,250 \
--numberOfHistogramBins 200 \
--samplingPercentage 1.0 \
--translationScale 250 \
--maximumStepLength 0.5 \
--minimumStepLength 0.001,0.001 \
--outputVolumePixelType float \
--initialTransform {initialMovingTxFilename} \
--transformType Rigid,Affine \
--fixedVolume {fixedImageFilename} \
--movingVolume {movingImageFilename} \
--outputVolume {BFitOut}
In [7]:
outputBFit = sitk.ReadImage(BFitOut,sitk.sitkFloat32)
myshow3d(outputBFit, xslices=[], yslices=[], zslices=[], title="OutputVolume of BRANSFit")
In [8]:
# Stage 1
# Set the registration filter
R = sitk.ImageRegistrationMethod()
# Transform of the first stage
tx1 = sitk.VersorRigid3DTransform()
tx1.SetFixedParameters(initialMovingTx.GetFixedParameters())
R.SetInitialTransform(tx1,inPlace=False)
# Set the initial transform of this stage
R.SetMovingInitialTransform(ctx)
# Set the other components of the registration filter
R.SetMetricFixedMask(fixedMask)
R.SetMetricMovingMask(movingMask)
R.SetMetricAsMattesMutualInformation( 200 )
R.SetOptimizerAsRegularStepGradientDescent(learningRate =1.0,
minStep=1e-3,
numberOfIterations=250,
gradientMagnitudeTolerance=1e-4,
estimateLearningRate=R.Never)
R.SetOptimizerScales([1, 1, 1, 1.0/250, 1.0/250, 1.0/250])
R.SetInterpolator(sitk.sitkLinear)
R.SetMetricSamplingPercentage(0.5)
R.SetMetricSamplingStrategy(R.RANDOM)
R.SetShrinkFactorsPerLevel([1])
R.SetSmoothingSigmasPerLevel([0])
Out[8]:
In [9]:
R.RemoveAllCommands()
R.AddCommand( sitk.sitkIterationEvent, lambda: command_iteration(R) )
outTx = R.Execute(sitk.Cast(fixed,sitk.sitkFloat32), sitk.Cast(moving,sitk.sitkFloat32))
print("-------")
print(outTx)
print("Optimizer stop condition: {0}".format(R.GetOptimizerStopConditionDescription()))
print(" Iteration: {0}".format(R.GetOptimizerIteration()))
print(" Metric value: {0}".format(R.GetMetricValue()))
In [10]:
# Update the composite transform with the output of the first stage
ctx.AddTransform(outTx)
cshow3d(fixed,moving,ctx)
In [11]:
# Stage 2
tx2 = sitk.AffineTransform(fixed.GetDimension())
tx2.SetFixedParameters(tx1.GetFixedParameters())
R.SetInitialTransform(tx2,inPlace=False)
R.SetMovingInitialTransform(ctx)
# Set a new optimizer
R.SetOptimizerAsConjugateGradientLineSearch(learningRate=1.0,
numberOfIterations=250,
convergenceMinimumValue=1e-6,
convergenceWindowSize=10,
lineSearchLowerLimit=0,
lineSearchUpperLimit=2,
lineSearchEpsilon=0.2,
lineSearchMaximumIterations=40,
estimateLearningRate=R.EachIteration,
maximumStepSizeInPhysicalUnits=1.0)
R.SetOptimizerScalesFromPhysicalShift()
affineOut = R.Execute(fixed,moving)
print("-------")
print(affineOut)
print("Optimizer stop condition: {0}".format(R.GetOptimizerStopConditionDescription()))
print(" Iteration: {0}".format(R.GetOptimizerIteration()))
print(" Metric value: {0}".format(R.GetMetricValue()))
In [12]:
# Update the composite transform with the output of the second stage
ctx.AddTransform(affineOut)
cshow3d(fixed,moving,ctx)
In [13]:
compare_with_baseline(fixed,moving,outputBFit,transform=ctx,numberOfPixelsTolerance=300,radiusTolerance=1,intensityTolerance=10)