itk中的配准整理

Perform 2D Translation Registration With Mean Squares

原文地址: https://examples.itk.org/src/registration/common/perform2dtranslationregistrationwithmeansquares/documentation

效果:

源码:

import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

import itk


# Perform 2D Translation Registration With Mean Squares
baseDir = r'D:\learn\itk\itk5.3.0\ITK-5.3.0\Examples\Data'
fixed_input_image = baseDir + r"\BrainProtonDensitySlice.png"
moving_input_image = baseDir + r"\BrainProtonDensitySliceShifted13x17y.png"

PixelType = itk.ctype("float")
fixedImage = itk.imread(fixed_input_image, PixelType)
movingImage = itk.imread(moving_input_image, PixelType)

Dimension = fixedImage.GetImageDimension()
FixedImageType = itk.Image[PixelType, Dimension]
MovingImageType = itk.Image[PixelType, Dimension]

TransformType = itk.TranslationTransform[itk.D, Dimension]
initialTransform = TransformType.New()

optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
    LearningRate=4,
    MinimumStepLength=0.001,
    RelaxationFactor=0.5,
    NumberOfIterations=200,
)

metric = itk.MeanSquaresImageToImageMetricv4[FixedImageType, MovingImageType].New()

registration = itk.ImageRegistrationMethodv4.New(
    FixedImage=fixedImage,
    MovingImage=movingImage,
    Metric=metric,
    Optimizer=optimizer,
    InitialTransform=initialTransform,
)

movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()
initialParameters[0] = 0
initialParameters[1] = 0
movingInitialTransform.SetParameters(initialParameters)
registration.SetMovingInitialTransform(movingInitialTransform)

identityTransform = TransformType.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)

registration.SetNumberOfLevels(1)
registration.SetSmoothingSigmasPerLevel([0])
registration.SetShrinkFactorsPerLevel([1])
registration.Update()

transform = registration.GetTransform()
finalParameters = transform.GetParameters()
translationAlongX = finalParameters.GetElement(0)
translationAlongY = finalParameters.GetElement(1)
numberOfIterations = optimizer.GetCurrentIteration()
bestValue = optimizer.GetValue()

print("Result = ")
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Iterations    = " + str(numberOfIterations))
print(" Metric value  = " + str(bestValue))

CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(movingInitialTransform)
outputCompositeTransform.AddTransform(registration.GetModifiableTransform())

resampler = itk.ResampleImageFilter.New(
    Input=movingImage,
    Transform=outputCompositeTransform,
    UseReferenceImage=True,
    ReferenceImage=fixedImage,
)
resampler.SetDefaultPixelValue(100)

OutputPixelType = itk.ctype("unsigned char")
OutputImageType = itk.Image[OutputPixelType, Dimension]
caster = itk.CastImageFilter[FixedImageType, OutputImageType].New(Input=resampler)
output_image_arr = itk.GetArrayFromImage(caster)  # 配准后的结果

difference = itk.SubtractImageFilter.New(Input1=fixedImage, Input2=resampler)
intensityRescaler = itk.RescaleIntensityImageFilter[
    FixedImageType, OutputImageType
].New(
    Input=difference,
    OutputMinimum=itk.NumericTraits[OutputPixelType].min(),
    OutputMaximum=itk.NumericTraits[OutputPixelType].max(),
)
resampler.SetDefaultPixelValue(1)
intensityRescaler.Update()
difference_image_after_arr = itk.GetArrayFromImage(intensityRescaler.GetOutput())  # 配准后与原图像的差别

resampler.SetTransform(identityTransform)
intensityRescaler.Update()
difference_image_before_arr = itk.GetArrayFromImage(intensityRescaler.GetOutput())  # 配准前与原图像的差别

# output_image  difference_image_after  difference_image_before
plt.subplot(231), plt.imshow(itk.GetArrayFromImage(fixedImage), 'gray'), plt.title('fixedImage')
plt.subplot(232), plt.imshow(itk.GetArrayFromImage(movingImage), 'gray'), plt.title('movingImage')
plt.subplot(233), plt.imshow(output_image_arr, 'gray'), plt.title('output_image')
plt.subplot(234), plt.imshow(difference_image_after_arr, 'gray'), plt.title('difference_image_after')
plt.subplot(235), plt.imshow(difference_image_before_arr, 'gray'), plt.title('difference_image_before')
# plt.subplot(236), plt.imshow(itk.GetArrayFromImage(moving_smoothed_image), 'gray'), plt.title('moving_smoothed_image')
plt.show()

多模态互信息配准 Perform Multi Modality Registration With Viola Wells Mutual Information

原文地址: https://examples.itk.org/src/registration/common/performmultimodalityregistrationwithmutualinformation/documentation

效果图

源码:

import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

import itk


dim = 2
ImageType = itk.Image[itk.F, dim]
FixedImageType = ImageType
MovingImageType = ImageType

baseDir = r'D:\learn\itk\itk5.3.0\ITK-5.3.0\Examples\Data'

fixed_img_path = baseDir + r"\BrainT1SliceBorder20.png"
moving_img_path = baseDir + r"\BrainProtonDensitySliceShifted13x17y.png"

fixed_img = itk.imread(fixed_img_path, itk.F)
moving_img = itk.imread(moving_img_path, itk.F)

fixed_normalized_image = itk.normalize_image_filter(fixed_img)
fixed_smoothed_image = itk.discrete_gaussian_image_filter(fixed_normalized_image, variance=2.0)

moving_normalized_image = itk.normalize_image_filter(moving_img)
moving_smoothed_image = itk.discrete_gaussian_image_filter(moving_normalized_image, variance=2.0)

TransformType = itk.TranslationTransform[itk.D, dim]
OptimizerType = itk.GradientDescentOptimizer
MetricType = itk.MutualInformationImageToImageMetric[ImageType, ImageType]
RegistrationType = itk.ImageRegistrationMethod[ImageType, ImageType]
InterpolatorType = itk.LinearInterpolateImageFunction[ImageType, itk.D]
n_iterations = 200
transform = TransformType.New()
metric = MetricType.New()
optimizer = OptimizerType.New()
registrar = RegistrationType.New()
interpolator = InterpolatorType.New()

registrar.SetFixedImage(fixed_smoothed_image)
registrar.SetMovingImage(moving_smoothed_image)
registrar.SetOptimizer(optimizer)
registrar.SetTransform(transform)
registrar.SetInterpolator(interpolator)
registrar.SetMetric(metric)
registrar.SetFixedImageRegion(fixed_img.GetBufferedRegion())
registrar.SetInitialTransformParameters(transform.GetParameters())
numberOfSamples = int(fixed_img.GetBufferedRegion().GetNumberOfPixels() * 0.01)
print('numberOfSamples:', numberOfSamples)
metric.SetNumberOfSpatialSamples(numberOfSamples)
metric.SetFixedImageStandardDeviation(0.4)
metric.SetMovingImageStandardDeviation(0.4)
metric.ReinitializeSeed(121212)

optimizer.SetLearningRate(15.0)
optimizer.SetNumberOfIterations(n_iterations)
optimizer.MaximizeOn()


def log_iteration():
    print('Iteration:', optimizer.GetCurrentIteration())
    print('Value:', optimizer.GetValue())
    print('CurrentPosition:', list(optimizer.GetCurrentPosition()))


optimizer.AddObserver(itk.IterationEvent(), log_iteration)
registrar.Update()
print('-'*50)
print(f"Its: {optimizer.GetCurrentIteration()}")
print(f"Final Value: {optimizer.GetValue()}")
print(f"Final Position: {list(registrar.GetLastTransformParameters())}")

ResampleFilterType = itk.ResampleImageFilter[MovingImageType, FixedImageType]
resample = ResampleFilterType.New(
    Transform=transform,
    Input=moving_img,
    Size=fixed_img.GetLargestPossibleRegion().GetSize(),
    OutputOrigin=fixed_img.GetOrigin(),
    OutputSpacing=fixed_img.GetSpacing(),
    OutputDirection=fixed_img.GetDirection(),
    DefaultPixelValue=100,
)
resample.Update()
output_image_arr = itk.GetArrayFromImage(resample.GetOutput())

CheckerBoardFilterType = itk.CheckerBoardImageFilter[FixedImageType]
checker = CheckerBoardFilterType.New()
checker.SetInput1(fixed_img)
checker.SetInput2(resample.GetOutput())
checker.Update()
difference_image_after_arr = itk.GetArrayFromImage(checker.GetOutput())

identityTransform = TransformType.New()
identityTransform.SetIdentity()
resample.SetTransform(identityTransform)
checker.Update()
difference_image_before_arr = itk.GetArrayFromImage(checker.GetOutput())

# output_image  difference_image_after  difference_image_before
plt.subplot(231), plt.imshow(itk.GetArrayFromImage(fixed_img), 'gray'), plt.title('fixed_img')
plt.subplot(232), plt.imshow(itk.GetArrayFromImage(moving_img), 'gray'), plt.title('moving_img')
plt.subplot(233), plt.imshow(output_image_arr, 'gray'), plt.title('output_image')
plt.subplot(234), plt.imshow(difference_image_after_arr, 'gray'), plt.title('difference_image_after')
plt.subplot(235), plt.imshow(difference_image_before_arr, 'gray'), plt.title('difference_image_before')
# plt.subplot(236), plt.imshow(itk.GetArrayFromImage(moving_smoothed_image), 'gray'), plt.title('moving_smoothed_image')
plt.show()

Register Image to Another Using Landmarks 通过标记点配准图像

原文地址:https://examples.itk.org/src/registration/common/registerimagetoanotherusinglandmarks/documentation

效果图

源码

import os
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

import itk


# Perform 2D Translation Registration With Mean Squares
baseDir = r'D:\learn\itk\itk5.3.0\ITK-5.3.0\Examples\Data'
fixed_input_image = baseDir + r"\BrainProtonDensitySlice.png"
moving_input_image = baseDir + r"\BrainProtonDensitySliceShifted13x17y.png"

PixelType = itk.ctype("float")
fixedImage = itk.imread(fixed_input_image, PixelType)
movingImage = itk.imread(moving_input_image, PixelType)

Dimension = fixedImage.GetImageDimension()
ImageType = itk.Image[PixelType, Dimension]

LandmarkPointType = itk.Point[itk.D, Dimension]
LandmarkContainerType = itk.vector[LandmarkPointType]

fixed_landmarks = LandmarkContainerType()
moving_landmarks = LandmarkContainerType()
fixed_point = LandmarkPointType()
moving_point = LandmarkPointType()
fixed_point[0] = 66
fixed_point[1] = 79
moving_point[0] = 99
moving_point[1] = 116
fixed_landmarks.push_back(fixed_point)
moving_landmarks.push_back(moving_point)
fixed_point[0] = 115
fixed_point[1] = 81
moving_point[0] = 148
moving_point[1] = 117
fixed_landmarks.push_back(fixed_point)
moving_landmarks.push_back(moving_point)
fixed_point[0] = 75
fixed_point[1] = 139
moving_point[0] = 108
moving_point[1] = 176
fixed_landmarks.push_back(fixed_point)
moving_landmarks.push_back(moving_point)

TransformInitializerType = itk.LandmarkBasedTransformInitializer[
    itk.Transform[itk.D, Dimension, Dimension]
]
transform_initializer = TransformInitializerType.New()
transform_initializer.SetFixedLandmarks(fixed_landmarks)
transform_initializer.SetMovingLandmarks(moving_landmarks)

transform = itk.Rigid2DTransform[itk.D].New()
transform_initializer.SetTransform(transform)
transform_initializer.InitializeTransform()

output = itk.resample_image_filter(
    movingImage,
    transform=transform,
    use_reference_image=True,
    reference_image=fixedImage,
    default_pixel_value=50,
)
print('transform:', transform)

output_image_arr = itk.GetArrayFromImage(output)  # 配准后的结果

OutputPixelType = itk.ctype("unsigned char")
OutputImageType = itk.Image[OutputPixelType, Dimension]

difference = itk.SubtractImageFilter.New(Input1=fixedImage, Input2=output)
intensityRescaler = itk.RescaleIntensityImageFilter[
    ImageType, OutputImageType
].New(
    Input=difference,
    OutputMinimum=itk.NumericTraits[OutputPixelType].min(),
    OutputMaximum=itk.NumericTraits[OutputPixelType].max(),
)
intensityRescaler.Update()
difference_image_after_arr = itk.GetArrayFromImage(intensityRescaler.GetOutput())  # 配准后与原图像的差别

# output_image  difference_image_after  difference_image_before
plt.subplot(221), plt.imshow(itk.GetArrayFromImage(fixedImage), 'gray'), plt.title('fixedImage')
plt.subplot(222), plt.imshow(itk.GetArrayFromImage(movingImage), 'gray'), plt.title('movingImage')
plt.subplot(223), plt.imshow(output_image_arr, 'gray'), plt.title('output_image')
plt.subplot(224), plt.imshow(difference_image_after_arr, 'gray'), plt.title('difference_image_after')
plt.show()