%matplotlib inline
import matplotlib.pyplot as plt

import SimpleITK as sitk

print(sitk.Version())
from myshow import myshow

# Download data to work on
%run update_path_to_download_script
from downloaddata import fetch_data as fdata

OUTPUT_DIR = "Output"

SimpleITK Version: 2.4.0 (ITK 5.4)
Compiled: Aug 15 2024 01:21:37

fixed_rgb = sitk.ReadImage(fdata("vm_head_rgb.mha"))
fixed_rgb = fixed_rgb[735:1330, 204:975, :]
fixed_rgb = sitk.BinShrink(fixed_rgb, [3, 3, 1])

Fetching vm_head_rgb.mha

moving = sitk.ReadImage(fdata("vm_head_mri.mha"))

Fetching vm_head_mri.mha

myshow(moving)

# Segment blue ice
seeds = [[10, 10, 10]]
fixed_mask = sitk.VectorConfidenceConnected(
    fixed_rgb,
    seedList=seeds,
    initialNeighborhoodRadius=5,
    numberOfIterations=4,
    multiplier=8,
)

# Invert the segment and choose largest component
fixed_mask = sitk.RelabelComponent(sitk.ConnectedComponent(fixed_mask == 0)) == 1

myshow(sitk.Mask(fixed_rgb, fixed_mask));

# pick red channel
fixed = sitk.VectorIndexSelectionCast(fixed_rgb, 0)

fixed = sitk.Cast(fixed, sitk.sitkFloat32)
moving = sitk.Cast(moving, sitk.sitkFloat32)

initialTransform = sitk.Euler3DTransform()
initialTransform = sitk.CenteredTransformInitializer(
    sitk.Cast(fixed_mask, moving.GetPixelID()),
    moving,
    initialTransform,
    sitk.CenteredTransformInitializerFilter.MOMENTS,
)
print(initialTransform)

itk::simple::Euler3DTransform
 Euler3DTransform (0x7f8df87dcfc0)
   RTTI typeinfo:   itk::Euler3DTransform<double>
   Reference Count: 1
   Modified Time: 3728
   Debug: Off
   Object Name: 
   Observers: 
     none
   Matrix: 
     1 0 0 
     0 1 0 
     0 0 1 
   Offset: [-0.305777, 2.55523, 135.031]
   Center: [0.853078, -32.3853, -144.027]
   Translation: [-0.305777, 2.55523, 135.031]
   Inverse: 
     1 0 0 
     0 1 0 
     0 0 1 
   Singular: 0
   AngleX: 0
   AngleY: 0
   AngleZ: 0
   ComputeZYX: Off

def command_iteration(method):
    print(
        f"{method.GetOptimizerIteration()} = {method.GetMetricValue()} : {method.GetOptimizerPosition()}",
        end="\n",
    )
    sys.stdout.flush();

tx = initialTransform
R = sitk.ImageRegistrationMethod()
R.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
R.SetOptimizerAsGradientDescentLineSearch(learningRate=1, numberOfIterations=100)
R.SetOptimizerScalesFromIndexShift()
R.SetShrinkFactorsPerLevel([4, 2, 1])
R.SetSmoothingSigmasPerLevel([8, 4, 2])
R.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
R.SetMetricSamplingStrategy(R.RANDOM)
# specifying a seed eliminates registration variability associated with the
# random sampling
R.SetMetricSamplingPercentage(percentage=0.1, seed=42)
R.SetInitialTransform(tx)
R.SetInterpolator(sitk.sitkLinear)

import sys

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(tx)
print(f"Optimizer stop condition: {R.GetOptimizerStopConditionDescription()}")
print(f" Iteration: {R.GetOptimizerIteration()}")
print(f" Metric value: {R.GetMetricValue()}")

0 = -0.7633063775373199 : (-0.0007251069363556011, -0.0004616650065359306, -0.00029872387211812005, -0.06285427361028123, 2.7316499067940554, 109.7703328160343)

1 = -1.1826751740226433 : (-0.0011684068679123868, -0.0008655481622207755, -0.00046279630880065466, -0.30863411589667533, 2.940196913878562, 107.77275330345998)

2 = -1.1860220026357704 : (-0.002246494426982868, -0.0017899175822320085, -0.000786754457434835, -0.7467617667792382, 3.3673346459401876, 109.49425731262846)

3 = -1.2009357450679656 : (-0.002341740740537487, -0.0018407596235358388, -0.0008253641732948884, -0.753240492974779, 3.3904028961788812, 109.08275356856258)

4 = -1.2015891475478269 : (-0.0024923686885144573, -0.0019244612147116405, -0.0008860834831289548, -0.766393621470172, 3.4269550263805386, 108.59996629816861)

5 = -1.2015458441069766 : (-0.002720248924479156, -0.002055953250176367, -0.000977055382048611, -0.7888867189246096, 3.4828483959313936, 108.19536711229534)

6 = -1.2011280524310775 : (-0.003731092760915537, -0.002649496200497969, -0.0013736111964535093, -0.8887844499139288, 3.7179029406428747, 107.8555022013922)

7 = -1.2025209799957448 : (-0.005267615447186065, -0.003448821792621109, -0.002000225362384655, -0.9480401025275874, 3.963201740079521, 109.11659258148521)

8 = -1.2078375789568716 : (-0.005395885489474075, -0.00350538757747309, -0.002058325699255972, -0.939383650637257, 3.9726453861356177, 108.58393446133363)

9 = -1.20680147640917 : (-0.005679610734367443, -0.0036401125632585336, -0.0021854480940744683, -0.928426945514153, 3.9961813547998197, 107.96676271497321)

tx = sitk.CompositeTransform([initialTransform, sitk.AffineTransform(3)])

R.SetOptimizerAsGradientDescentLineSearch(learningRate=1, numberOfIterations=100)
R.SetOptimizerScalesFromIndexShift()
R.SetShrinkFactorsPerLevel([2, 1])
R.SetSmoothingSigmasPerLevel([4, 1])
R.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
R.SetInitialTransform(tx)

outTx = R.Execute(
    sitk.Cast(fixed, sitk.sitkFloat32), sitk.Cast(moving, sitk.sitkFloat32)
)
R.GetOptimizerStopConditionDescription()

0 = -0.9246489027264059 : (0.9974089908558552, 0.00032056228165761844, 5.297830597443564e-05, 0.00024604716550388116, 0.9982222039714761, -2.614438023337558e-05, -0.0022798356115667894, 0.014886279113618249, 0.9963704151924245, -0.03139888493098891, 0.11756924709526634, -0.41062839425054326)

1 = -0.9329251612157431 : (0.9961749786784537, 0.0004609968345551488, 5.613374670103203e-05, 0.00035762348818118856, 0.9974467421133562, 6.204854435219327e-05, -0.0022037096800858755, 0.01574145770293372, 0.9923283190620598, -0.033951733435116344, 0.11078474400363453, 0.8216753604241602)

2 = -0.9425675484883225 : (0.9958706287141689, 0.0004935470848067009, 4.203973409404129e-05, 0.0003834253051336116, 0.9972271541239995, 0.00011006747372080706, -0.0024773068663882427, 0.01737741151722013, 0.9921159311240098, -0.028980145238754352, 0.10034515001213463, 0.48552444026831476)

3 = -0.9439829114314803 : (0.9950607387816928, 0.0005841741387997406, 1.181229515171041e-05, 0.000457327414101372, 0.9966498410979091, 0.00022906224832368763, -0.0029882093296480614, 0.02026541582293507, 0.9908417279521354, -0.018718005621909874, 0.07282199758671677, 0.11589828606011865)

4 = -0.9447162295820674 : (0.9913125678083762, 0.0009898387952096911, -0.00010690574384879314, 0.0007809755631358855, 0.9940465899678642, 0.0007672792593900695, -0.0038067523461829647, 0.025127700414040267, 0.9816201419309787, 0.023578336142760052, -0.07851669428640713, 0.9097848754926282)

5 = -0.9534880135550752 : (0.9910101204131582, 0.0010244093217783494, -0.00012817374271016267, 0.0008014833536019125, 0.993801612503512, 0.0008433963450818441, -0.004151410396904882, 0.02687404814784744, 0.9818370626352215, 0.031570976340459814, -0.10428417193052229, 0.20175019122840498)

6 = -0.9564247354763359 : (0.9903912900478714, 0.0010982845231765487, -0.0001608324278799042, 0.0008494578842833573, 0.9933143472450819, 0.0009698417844007319, -0.004432321128259562, 0.027818154356902134, 0.9808336543534639, 0.04406450391953347, -0.15022204068584988, -0.1645013256975954)

7 = -0.9571498786128856 : (0.9872597398582902, 0.0014688806542792496, -0.00029016742322651777, 0.0010814149860426767, 0.9908796471023696, 0.0015060193827100632, -0.00493872000692705, 0.028336570981587123, 0.9732761836959808, 0.09412695543876856, -0.3468273803055635, 0.04978144712258992)

8 = -0.9638760405748776 : (0.9870235117160736, 0.0015063989538739542, -0.00030108477528093897, 0.0010986947924285283, 0.9906506307045011, 0.00154840762910178, -0.005134636906454853, 0.02968317778666504, 0.9733496544056116, 0.09790117281994681, -0.3590765246833331, -0.39943152225587164)

9 = -0.9651801825477209 : (0.9864179893153423, 0.00160021903801876, -0.00031890129865016337, 0.0011413713939034524, 0.990095075602002, 0.0016431511294565369, -0.005399728010554285, 0.031096572871528474, 0.9724605918455429, 0.10405636863818245, -0.3891055602793314, -0.7306461277580022)

resample = sitk.ResampleImageFilter()
resample.SetReferenceImage(fixed_rgb)
resample.SetInterpolator(sitk.sitkBSpline3)
resample.SetTransform(outTx)
resample.AddCommand(
    sitk.sitkProgressEvent,
    lambda: print(f"\rProgress: {100*resample.GetProgress():03.1f}%...", end=""),
)
resample.AddCommand(sitk.sitkProgressEvent, lambda: sys.stdout.flush())
resample.AddCommand(sitk.sitkEndEvent, lambda: print("Done"))
out = resample.Execute(moving)

Progress: 0.0%...

Progress: 1.0%...

Progress: 17.8%...

Progress: 1.8%...

Progress: 4.3%...

Progress: 3.6%...

Progress: 5.5%...

Progress: 7.3%...

Progress: 9.1%...

Progress: 10.9%...

out_rgb = sitk.Cast(
    sitk.Compose([sitk.RescaleIntensity(out)] * 3), sitk.sitkVectorUInt8
)
vis_xy = sitk.CheckerBoard(fixed_rgb, out_rgb, checkerPattern=[8, 8, 1])
vis_xz = sitk.CheckerBoard(fixed_rgb, out_rgb, checkerPattern=[8, 1, 8])
vis_xz = sitk.PermuteAxes(vis_xz, [0, 2, 1])

myshow(vis_xz, dpi=30)

import os

sitk.WriteImage(out, os.path.join(OUTPUT_DIR, "example_registration.mha"))
sitk.WriteImage(vis_xy, os.path.join(OUTPUT_DIR, "example_registration_xy.mha"))
sitk.WriteImage(vis_xz, os.path.join(OUTPUT_DIR, "example_registration_xz.mha"))

Registration Example¶