<jupyter_start><jupyter_text>Registration Example<jupyter_code>%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"<jupyter_output><empty_output><jupyter_text>This section of the Visible Human Male is about 1.5GB. To expedite processing and registration we crop the region of interest, and reduce the resolution. Take note that the physical space is maintained through these operations.<jupyter_code>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])
moving = sitk.ReadImage(fdata("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)
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()}")
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()
resample = sitk.ResampleImageFilter()
resample.SetReferenceImage(fixed_rgb)
resample.SetInterpolator(sitk.sitkBSpline)
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)
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"))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11998154130133826, 'fraction_numerical': 0.02007383479464698, 'mean_word_length': 7.384912959381045, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '179007', 'license': 'Apache-2.0', 'n_tokens_mistral': 1616, 'n_tokens_neox': 1525, 'n_words': 264, 'path': 'Python/56_VH_Registration1.ipynb', 'repo_name': 'chirayukong/SimpleITK-Notebooks', 'size': '4334'}