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 cmimport 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 cmimport itkdim = 2
ImageType = itk.Image[itk.F, dim]
FixedImageType = ImageType
MovingImageType = ImageTypebaseDir = 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 cmimport 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()