In [ ]:

    

## Boiler plate code common to many notebooks.  See the TestFilesCommonCode.ipynb for details
from __future__ import print_function
%run TestFilesCommonCode.ipynb


    

In [ ]:

    

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

    

fixedImageFilename = fdata('BRAINSTools/test.nii.gz')
movingImageFilename = fdata('BRAINSTools/rotation.rescale.rigid.nii.gz')
fixed = sitk.ReadImage(fixedImageFilename,sitk.sitkFloat32)
moving = sitk.ReadImage(movingImageFilename, sitk.sitkFloat32)

myshow3d(fixed, xslices=[], yslices=[], zslices=[], title="fixed")
myshow3d(moving, xslices=[], yslices=[], zslices=[], title="moving")


    

In [ ]:

    

cshow3d(fixed,moving)


    

In [ ]:

    

BFitOut=os.path.realpath('Data/BFitOut.nii.gz')
# RUN THE BRAINSFIT FIRST ######################
!{BRAINSFitPath} \
--costMetric MMI \
--numberOfIterations 250,250,250 \
--numberOfHistogramBins 200 \
--samplingPercentage 0.5 \
--translationScale 250 \
--maximumStepLength 0.5 \
--minimumStepLength 0.01,0.003,0.001 \
--outputVolumePixelType float \
--initializeTransformMode useMomentsAlign \
--transformType Rigid,ScaleVersor3D,Affine \
--fixedVolume {fixedImageFilename} \
--movingVolume {movingImageFilename} \
--outputVolume {BFitOut}


    

In [ ]:

    

outputBFit = sitk.ReadImage(BFitOut,sitk.sitkFloat32)
myshow3d(outputBFit, xslices=[], yslices=[], zslices=[], title="OutputVolume of BRANSFit")


    

In [ ]:

    

# First transform is passed through an initializer
tx1 = sitk.CenteredVersorTransformInitializer(fixed, moving, sitk.VersorRigid3DTransform())


    

In [ ]:

    

ctx = sitk.Transform(tx1) # Set composite transform
                          # This composite transform is update at each stage
                          # and finally will be considered as the output of 
                          # the registration process since the InPlace is TRUE.
ctx.SetFixedParameters(ctx.GetFixedParameters()) # hack to force deep copy, as registion is done in place..
# Set the registration filter
R = sitk.ImageRegistrationMethod()
R.SetMetricAsMattesMutualInformation( 200 )
R.SetOptimizerAsRegularStepGradientDescent(learningRate =0.5,
                                           minStep=1e-2,
                                           numberOfIterations=250,
                                           gradientMagnitudeTolerance=1e-4,
                                           estimateLearningRate=R.Never)
R.SetOptimizerScales([1, 1, 1, 1.0/250, 1.0/250, 1.0/250])
R.SetShrinkFactorsPerLevel([1])
R.SetSmoothingSigmasPerLevel([0])
R.SetInitialTransform(ctx)
R.SetInterpolator(sitk.sitkLinear)
R.SetMetricSamplingPercentage(0.5)
R.SetMetricSamplingStrategy(R.RANDOM)


    

In [ ]:

    

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

    

cshow3d(fixed,moving,ctx)


    

In [ ]:

    

# Set the second stage transform
tx2 = sitk.ScaleVersor3DTransform()
# Fixed parameters of the next transform are taken from the last transform.
tx2.SetFixedParameters(tx1.GetFixedParameters())


    

In [ ]:

    

ctx.AddTransform(tx2) # Update the composite transform
ctx.SetFixedParameters(ctx.GetFixedParameters())
R.SetInitialTransform(ctx)
R.SetOptimizerAsRegularStepGradientDescent(learningRate =0.5,
                                           minStep=3e-3,
                                           numberOfIterations=250,
                                           gradientMagnitudeTolerance=1e-4,
                                           estimateLearningRate=R.Never)
R.SetOptimizerScales([1, 1, 1, 1.0/250, 1.0/250, 1.0/250, 1.0/25, 1.0/25, 1.0/25])

scaleTxOut = R.Execute(fixed,moving)
print("-------")
print(scaleTxOut)
print("Optimizer stop condition: {0}".format(R.GetOptimizerStopConditionDescription()))
print(" Iteration: {0}".format(R.GetOptimizerIteration()))
print(" Metric value: {0}".format(R.GetMetricValue()))


    

In [ ]:

    

cshow3d(fixed,moving,ctx)


    

In [ ]:

    

# Set the third stage
tx3 = sitk.AffineTransform(fixed.GetDimension())
tx3.SetFixedParameters(tx2.GetFixedParameters())


    

In [ ]:

    

ctx.AddTransform(tx3) # Update the composite transform
ctx.SetFixedParameters(ctx.GetFixedParameters())
R.SetInitialTransform(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=20,
                                            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 [ ]:

    

cshow3d(fixed,moving,ctx)


    

In [ ]:

    

compare_with_baseline(fixed,moving,outputBFit,transform=ctx,numberOfPixelsTolerance=100,radiusTolerance=0,intensityTolerance=9)