Elementary ITK usage learned in Python
Notification
I wrote this article as a memo of my study. The reference source is this notebook. https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization
ITK A library that is indispensable for image processing in the biomedical field, which specializes in image visualization technology. It is used in Slicer and ImageJ.
environment
Google Colab (Launch any notebook.)
Sample data
First download the sample data from GitHub using subversion.
!sudo apt install subversion
#Get the data
#Included in Github URL/tree/master/Replace with trunk
!svn checkout https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization/trunk/data
ITK installation
#Install ITK
!pip install itk
Image Filtering
What to learn here
--Get to know the image processing pipeline model used by ITK. --Images and Meshes are used to represent data objects. --A new data object is created within the process of manipulating the data object. --Processes are classified as process objects and have source, filter and mapper objects. --source (reader, etc.) gets the data and filters the objects to generate new data. --mapper accepts data for output to a file or other system.
The term filter is widely used in three types:
--Data pipeline Data and process objects are typically converted to each other using the SetInput () and GetOutput () methods. No new output or new pixel data for any processing will occur until the Update () method is called at the end of the pipeline (process or data object). --Pipeline updates Multidimensional image data is large and can be made even larger. ITK addresses this issue through its data streaming capabilities. --Streaming Streaming is done by splitting the image into non-overlapping areas at the end of the pipeline. Then propagate the pipeline as a Requested Region.
There are three types of Image Regions in ITK.
--BufferedRegion: Pixel area stored in memory --LargestPossibleRegion: Maximum possible area of the image --RequestedRegion: The area requested by a single processing path during streaming (BufferedRegion and LargestPossibleRegion are the same as or larger than RequestedRegion)
Try.
import numpy as np
import itk
import matplotlib.pyplot as plt
#Practice spatial filtering with Pac-Man
fileName = 'data/PacMan.png'
reader = itk.ImageFileReader.New(FileName=fileName)
#Try using a Gaussian filter
smoother = itk.RecursiveGaussianImageFilter.New(Input=reader.GetOutput())
print("reader's Output: %s" % reader.GetOutput())
print("smoother's Output: %s" % smoother.GetOutput())
#Update to generate processing results()Don't forget.
smoother.Update()
print("reader's Output: %s" % reader.GetOutput())
print("smoother's Output: %s" % smoother.GetOutput())
image = reader.GetOutput()#Confirmation of the original
plt.subplot(1,2,1),plt.title("original"),plt.imshow(itk.GetArrayFromImage(image))
smoothed = smoother.GetOutput()#After smoothing
plt.subplot(1,2,2),plt.title("smoothing"),plt.imshow(itk.GetArrayFromImage(smoothed))
If you update with a pipeline object (although you have already done so), the processing result will be the latest
smoother.Update()
#However, if you change the process, you need to update again.
smoother.SetSigma(10.0)
smoother.Update()
plt.imshow(itk.GetArrayFromImage(smoothed))
Also, I don't think there are many, but if you make any changes to the reader in the middle, don't forget to update the reader.
# reader.Modified()
#You can try different filters.
# https://itk.org/Doxygen/html/group__IntensityImageFilters.html
fileName = 'data/PacMan.png'
reader = itk.ImageFileReader.New(FileName=fileName)
smoother = itk.MedianImageFilter.New(Input=reader.GetOutput())#For example, median
smoother.SetRadius(5)
smoother.Update()
plt.imshow(itk.GetArrayFromImage(smoother.GetOutput()))
Supplement
#Place StreamingImageFilter at the end of the pipeline to stream the pipeline.
#smoother generates output multiple times for each streaming division of the image area.
#The reader can't stream, so the output is only generated once.
streamer = itk.StreamingImageFilter.New(Input=smoother.GetOutput())
streamer.SetNumberOfStreamDivisions(4)
reader.Modified()
streamer.Update()
#You only have to specify the writing including the extension
writer = itk.ImageFileWriter.New(Input=smoother.GetOutput())
writer.SetFileName('my_output.mha')
writer.SetNumberOfStreamDivisions(4)
writer.Update()
Image Segmentation
What to learn here
--Know the segmentation methods available in ITK (First, ConfidenceConnectedImageFilter) --Understand how region glowing and its parameters change behavior --Understand how level sets and their parameters change their behavior
Types of segmentation available in ITK
- Level Set --Watershed (gradient descent of image features and analysis along domain boundaries) --Connected component (an abstract technique for working with connected objects used in binarization and labeling) --Label voting (for counting pixel values as labels and using them in the process)
- Region Growing --Classifier (machine learning algorithm)
Here, we will introduce a processing example for Region Growing and Level Set.
Region growing
The basic approach of the region expansion algorithm is to start expansion from the seed area (usually one or more pixels) that is considered to be inside the segmented object. Pixels adjacent to this segmented object area are evaluated to determine if they should be considered part of the object. If segmented, they are added to the area and the process continues as long as new pixels are added to the area. The region expansion algorithm is the criteria used to determine whether a region contains pixels, the type of connection used to determine neighbors, and the process (or process) used to search for neighboring pixels. It depends on the strategy).
Here, the following method is used as an example.
- ConfidenceConnected
- ConnectedThreshold
- IsolatedConnected
First, check how to use "Confidence Connected". ConfidenceConnected is a method of determining the upper and lower limits of the domain boundary, finding the standard deviation from the values in between, and setting the allowable range of this standard deviation (like leverage. SD * Multiplier).
# BrainProtonDensitySlice.Use png
# 'unsigned char'The reason for loading with is to show the data conversion method later.
input_image = itk.imread('data/BrainProtonDensitySlice.png', itk.ctype('unsigned char'))
# View the input image
plt.imshow(itk.GetArrayFromImage(input_image))
Check how to use the "ConfidenceConnectedImageFilter" class. If you don't know how to use it, look it up.
confidence_connected = itk.ConfidenceConnectedImageFilter.New(input_image)
#Check how to use using the help function of python
help(confidence_connected)
#Set confidence connected filter parameter arbitrarily
confidence_connected.SetMultiplier(2.3)#Tolerance, ± SD in this example*2.3
confidence_connected.SetNumberOfIterations(5)
confidence_connected.SetInitialNeighborhoodRadius(3)
confidence_connected.SetReplaceValue(255)
If you don't know how to use the function, look it up.
# *Multiplier*?
confidence_connected.SetMultiplier?
Set the starting point "seed point" of area expansion. Note that the origin of the ITK Image object is at the bottom left of the image coordinates. (NumPy is in the upper left)
# Set the seed points
confidence_connected.SetSeed([100, 110])
#Run and view
confidence_connected.Update()
plt.imshow(itk.GetArrayFromImage(confidence_connected))
#Try to improve the processing result. Apply a median filter to the original image.
median_filtered_image = itk.median_image_filter(input_image, radius=1)
#Update the parameters again. Input the image after median filter.
confidence_connected_image = itk.confidence_connected_image_filter(median_filtered_image,
seed=[100,110],
replace_value=255,
multiplier=3.0,
number_of_iterations=5,
initial_neighborhood_radius=3)
plt.imshow(itk.GetArrayFromImage(confidence_connected_image))
Next, let's compare the three region glowing processes. ConnectedThresholdImageFilter vs ConfidenceConnectedmageFilter vs IsolatedConnectedImageFilter
# ConnectedThresholdImageFilter vs ConfidenceConnectedmageFilter vs IsolatedConnectedImageFilter
input_image = itk.imread('data/BrainProtonDensitySlice.png', itk.ctype('unsigned char'))
connected_threshold = itk.ConnectedThresholdImageFilter.New(input_image)
connected_threshold.SetLower(190)
connected_threshold.SetUpper(255)
connected_threshold.SetReplaceValue( 255 )
connected_threshold.SetSeed( [100, 110] );
connected_threshold.Update()
confidence_connected = itk.ConfidenceConnectedImageFilter.New(input_image)
confidence_connected.SetMultiplier(2.3)
confidence_connected.SetSeed( [100, 110] );
confidence_connected.AddSeed( [80, 126] );
confidence_connected.SetNumberOfIterations(5)
confidence_connected.SetInitialNeighborhoodRadius(3)
confidence_connected.SetReplaceValue(255)
confidence_connected.Update()
isolated_connected = itk.IsolatedConnectedImageFilter.New(input_image)
isolated_connected.SetSeed1([98, 112])
isolated_connected.SetSeed2([98, 136])
isolated_connected.SetUpperValueLimit( 245 )
isolated_connected.FindUpperThresholdOff();
isolated_connected.Update()
plt.subplot(131)
plt.title("ConnectedThreshold")
plt.imshow(itk.GetArrayFromImage(connected_threshold))
plt.subplot(132)
plt.title("ConfidenceConnected")
plt.imshow(itk.GetArrayFromImage(confidence_connected))
plt.subplot(133)
plt.title("IsolatedConnected")
plt.imshow(itk.GetArrayFromImage(isolated_connected))
plt.tight_layout()
plt.show()
Fast Marching Segmentation(LevelSet) Next, try the Fast Marching Segmentation of the LevelSet method. If the target area of segmentation is a simple structure, a fast expansion algorithm called fast marching can be used. Simplify the image with a spatial filter before trying.
#Find out which image data types the smoothing filter can handle
itk.CurvatureAnisotropicDiffusionImageFilter.GetTypes()
#How to cast image data to supported data types
Dimension = input_image.GetImageDimension()
FloatImageType = itk.Image[itk.ctype('float'), Dimension]
caster = itk.CastImageFilter[input_image, FloatImageType].New()
caster.SetInput(input_image)
#Apply a smoothing filter.
smoothed_image = itk.curvature_anisotropic_diffusion_image_filter(caster,
time_step=0.125,
number_of_iterations=5,
conductance_parameter=3.0)
plt.imshow(itk.GetArrayFromImage(smoothed_image))
#Execute edge detection filter as it is
gradient = itk.gradient_magnitude_recursive_gaussian_image_filter(smoothed_image, sigma=1.0)
plt.imshow(itk.GetArrayFromImage(gradient))
#Continue to run the sigmoid filter
basin_minimum = 2.25
border_minimum = 3.75
alpha = - (border_minimum - basin_minimum) / 3.0
beta = (border_minimum + basin_minimum) / 2.0
sigmoid = itk.sigmoid_image_filter(gradient, output_minimum=0.0, output_maximum=1.0, alpha=alpha, beta=beta)
plt.imshow(itk.GetArrayFromImage(sigmoid))
#High speed marching
fast_marching = itk.FastMarchingImageFilter.New(sigmoid)
seed_value = 0.0
NodeType = itk.LevelSetNode[itk.ctype('float'), Dimension]
node = NodeType()
node.SetIndex([100,110])
node.SetValue(seed_value)
NodeContainerType = itk.VectorContainer[itk.ctype('unsigned int'), NodeType]
nodes = NodeContainerType.New()
nodes.Initialize()
nodes.InsertElement(0, node)
fast_marching_image = itk.fast_marching_image_filter(sigmoid, trial_points=nodes, stopping_value=80)
plt.imshow(itk.GetArrayFromImage(fast_marching_image))
Binarize the image.
time_threshold = 100 #Set threshold
thresholded_image = itk.binary_threshold_image_filter(fast_marching_image,
lower_threshold = 0,
upper_threshold=time_threshold,
outside_value=0,
inside_value=255)
plt.imshow(itk.GetArrayFromImage(thresholded_image))
From here, try the Shape Detection Level Set. The processing speed is improved, but the curvature is not good.
The pipeline looks like this:
First, create a Balanced image in the lower part of the figure.
#Binary image is White Matter.Use png
binary_mask = itk.imread('data/WhiteMatter.png', itk.ctype('float'))
smoothed_image = itk.smoothing_recursive_gaussian_image_filter(binary_mask)
#Gradation processing
rescaler = itk.IntensityWindowingImageFilter[FloatImageType,FloatImageType].New(smoothed_image)
rescaler.SetWindowMinimum(0)
rescaler.SetWindowMaximum(255)
#Invert the input level to have a negative value internally
rescaler.SetOutputMinimum(5.0)
rescaler.SetOutputMaximum(-5.0)
rescaler.Update()
plt.imshow(itk.GetArrayFromImage(rescaler))
Next, create the feature image in the figure.
#Pre-process using the loaded image in the above procedure
gradient = itk.gradient_magnitude_recursive_gaussian_image_filter(caster, sigma=1.0)
basin_minimum = 1
border_minimum = 5.0
alpha = - (border_minimum - basin_minimum) / 3.0
beta = (border_minimum + basin_minimum) / 2.0
#Feature image in the figure
sigmoid = itk.sigmoid_image_filter(gradient, output_minimum=0.0, output_maximum=1.0, alpha=alpha, beta=beta)
plt.imshow(itk.GetArrayFromImage(sigmoid))
#shape_detection_level_Execute set processing
shape_detection_image = itk.shape_detection_level_set_image_filter(rescaler,
feature_image=sigmoid,
maximum_r_m_s_error=0.001,
number_of_iterations=100,
propagation_scaling=1.0,
curvature_scaling=2.0)
#Binarization
thresholded_image = itk.binary_threshold_image_filter(shape_detection_image,
lower_threshold=-1e7,
# Cut at the zero set
upper_threshold=0.0,
outside_value=0,
inside_value=255)
plt.imshow(itk.GetArrayFromImage(thresholded_image))
Let's introduce some parameters related to the conditions of curvature (Curvature) and propagation (Propagation) of the level set into the shape detection level set. Min Basin Min Boundary Propagation Scaling Curvature Scaling Number of Iterations
def explore_shape_detection_image_parameters(basin_minimum, boundary_minimum, propagation_scaling,
curvature_scaling, number_of_iterations):
input_image = itk.imread('data/BrainProtonDensitySlice.png', itk.ctype('unsigned char'))
# Prepare the initial level set image
binary_mask = itk.imread('data/WhiteMatter.png', itk.ctype('float'))
smoothed_image = itk.smoothing_recursive_gaussian_image_filter(binary_mask)
Dimension = input_image.GetImageDimension()
FloatImageType = itk.Image[itk.ctype('float'), Dimension]
caster = itk.CastImageFilter[input_image, FloatImageType].New()
caster.SetInput(input_image)
# Prepare the speed image
gradient = itk.gradient_magnitude_recursive_gaussian_image_filter(caster, sigma=1.0)
alpha = - (boundary_minimum - basin_minimum) / 3.0
beta = (boundary_minimum + basin_minimum) / 2.0
sigmoid = itk.sigmoid_image_filter(gradient, output_minimum=0.0, output_maximum=1.0, alpha=alpha, beta=beta)
rescaler = itk.IntensityWindowingImageFilter[FloatImageType,FloatImageType].New(smoothed_image)
rescaler.SetWindowMinimum(0)
rescaler.SetWindowMaximum(255)
# Invert the input level set to have negative values inside
rescaler.SetOutputMinimum(5.0)
rescaler.SetOutputMaximum(-5.0)
rescaler.Update()
shape_detection_image = itk.shape_detection_level_set_image_filter(rescaler,
feature_image=sigmoid,
maximum_r_m_s_error=0.001,
number_of_iterations=number_of_iterations,
propagation_scaling=propagation_scaling,
curvature_scaling=curvature_scaling)
thresholded_image = itk.binary_threshold_image_filter(shape_detection_image,
lower_threshold=-1e7,
# Cut at the zero set
upper_threshold=0.0,
outside_value=0,
inside_value=255)
plt.imshow(itk.GetArrayFromImage(thresholded_image))
plt.title(str(basin_minimum)+"-"+str(boundary_minimum)+"-"+str(propagation_scaling)+"-"+str(curvature_scaling)+"-"+str(number_of_iterations))
plt.show()
return
explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=5.,
propagation_scaling=1., curvature_scaling=2., number_of_iterations=100)
explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=5.,
propagation_scaling=1., curvature_scaling=2., number_of_iterations=500)
explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=5.,
propagation_scaling=1., curvature_scaling=4., number_of_iterations=100)
explore_shape_detection_image_parameters(basin_minimum=1, boundary_minimum=20.,
propagation_scaling=1., curvature_scaling=1., number_of_iterations=100)
Image Registration
What to learn here
--Know how ITK performs calculations in physical space --Understanding why ITK performs registration in physical space instead of pixel space --Touch the components of the ITK registration framework and know the available values
#Are the target data calibrated in the same geometric space before alignment? If you need resampling, do it.
#I wonder if I will do it with an IJ sample. future work
# https://imagej.nih.gov/ij/images/
# https://imagej.nih.gov/ij/images/pet-series/
# https://imagej.nih.gov/ij/images/ct-scan.zip
#Load the data to be used (reference image)
PixelType = itk.ctype('float')
fixedImage = itk.imread('data/BrainProtonDensitySliceBorder20.png', PixelType)
plt.imshow(itk.array_view_from_image(fixedImage))
#Image that is out of position
movingImage = itk.imread('data/BrainProtonDensitySliceShifted13x17y.png', PixelType)
plt.imshow(itk.array_view_from_image(movingImage))
#Gather the information needed for the conversion
#Number of dimensions
Dimension = fixedImage.GetImageDimension()
#Image type
FixedImageType = type(fixedImage)
MovingImageType = type(movingImage)
#Type of data to convert
TransformType = itk.TranslationTransform[itk.D, Dimension]
#Initialization of saucer
initialTransform = TransformType.New()
#Define a function that optimizes to minimize alignment error
optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
LearningRate=4,
MinimumStepLength=0.001,
RelaxationFactor=0.5,
NumberOfIterations=200)
#Set alignment algorithm (mean squared error)
metric = itk.MeanSquaresImageToImageMetricv4[
FixedImageType, MovingImageType].New()
#Initialize the alignment object
registration = itk.ImageRegistrationMethodv4.New(FixedImage=fixedImage,
MovingImage=movingImage,
Metric=metric,
Optimizer=optimizer,
InitialTransform=initialTransform)
#Alignment settings
#Conversion type of moving image
movingInitialTransform = TransformType.New()
initialParameters = movingInitialTransform.GetParameters()
initialParameters[0] = 0 #x
initialParameters[1] = 0 #y
movingInitialTransform.SetParameters(initialParameters)
registration.SetMovingInitialTransform(movingInitialTransform)
#Criteria conversion type
identityTransform = TransformType.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)
registration.SetNumberOfLevels(1) # number of multi-resolution levels
registration.SetSmoothingSigmasPerLevel([0]) #Smoothing level
registration.SetShrinkFactorsPerLevel([1]) #Shrink level
#Update the alignment settings and confirm
registration.Update()
#Get the result of conversion process
transform = registration.GetTransform()
finalParameters = transform.GetParameters()
#Moving distance
translationAlongX = finalParameters.GetElement(0)
translationAlongY = finalParameters.GetElement(1)
#Number of repetitions
numberOfIterations = optimizer.GetCurrentIteration()
#Method accuracy (smaller is better)
bestValue = optimizer.GetValue()
print("Result = ")
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Iterations = " + str(numberOfIterations))
print(" Metric value = " + str(bestValue))
#Image composition
CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(movingInitialTransform)
outputCompositeTransform.AddTransform(registration.GetModifiableTransform())
#Resampling
resampler = itk.ResampleImageFilter.New(Input=movingImage,
Transform=outputCompositeTransform,
UseReferenceImage=True,
ReferenceImage=fixedImage)
#Result display
resampler.SetDefaultPixelValue(100)
OutputPixelType = itk.ctype('unsigned char')
OutputImageType = itk.Image[OutputPixelType, Dimension]
resampler.Update()
result = resampler.GetOutput()
plt.subplot(1,2,1)
plt.title("registered")
plt.imshow(itk.array_view_from_image(result))
plt.subplot(1,2,2)
plt.title("moving_org")
plt.imshow(itk.array_view_from_image(movingImage))
Let's see the misalignment by the difference
difference = itk.SubtractImageFilter.New(Input1=fixedImage,
Input2=resampler)
resampler.SetDefaultPixelValue(1)
difference.Update()
plt.title("difference")
plt.imshow(itk.array_view_from_image(difference.GetOutput()))
#Misalignment before and after registration
resampler.SetTransform(identityTransform)
difference.Update()
plt.imshow(itk.array_view_from_image(difference.GetOutput()))
Next, different types of images (multimodal) are aligned. Use Mattes Mutual Information, which is often used for alignment in multimodal cases. Use rigid transformation.
#Reference image
t1_image = itk.imread('data/DzZ_T1.nrrd', itk.F)
print("Pixel spacing:", t1_image.GetSpacing())
plt.imshow(itk.array_view_from_image(t1_image))
#Rotating T2-weighted image
t2_image = itk.imread('data/DzZ_T2.nrrd', itk.F)
print("Pixel spacing:", t2_image.GetSpacing())
plt.imshow(itk.array_view_from_image(t2_image))
#Supplement: Get direction
pyDirM=itk.GetArrayFromVnlMatrix(t2_image.GetDirection().GetVnlMatrix().as_matrix())
print("Direction matrix:\n", pyDirM)
Alignment is performed using the correlation by Mattes Mutual Information, which is often used for alignment in multimodal cases.
T1ImageType = type(t1_image)
T2ImageType = type(t2_image)
assert T1ImageType==T2ImageType, "Images must have the same pixel type!"
#Initialize transformation type as rigid transformation
TransformTypeR = itk.Rigid2DTransform[itk.D]
initialTransform = TransformTypeR.New()
MetricType = itk.MattesMutualInformationImageToImageMetricv4[
T1ImageType, T2ImageType]
metric = MetricType.New()
scales = itk.OptimizerParameters[itk.D](initialTransform.GetNumberOfParameters())
translation_scale = 1.0 / 1000
scales.SetElement(0, 1.0)
scales.SetElement(1, translation_scale)
scales.SetElement(2, translation_scale)
#It can automatically estimate the scale, but it doesn't seem to work with Mattes Mutual Information.
#ScalesEstimatorType = itk.RegistrationParameterScalesFromPhysicalShift[MetricType]
#scalesEstimator = ScalesEstimatorType.New(Metric=metric)
optimizer = itk.RegularStepGradientDescentOptimizerv4.New(
Scales=scales,
#ScalesEstimator=scalesEstimator,
#LearningRate=1.0,
MinimumStepLength=0.001,
# RelaxationFactor=0.8,
NumberOfIterations=300)
registration = itk.ImageRegistrationMethodv4.New(FixedImage=t1_image,
MovingImage=t2_image,
Metric=metric,
MetricSamplingPercentage=0.1, # 10%
Optimizer=optimizer,
InitialTransform=initialTransform)
initialParameters = initialTransform.GetParameters()
#Offset setting
initialParameters[0] = 0.0 # rotation in radians
initialParameters[1] = -40.0 # x translation in millimeters
initialParameters[2] = +30.0 # y translation in millimeters
initialTransform.SetParameters(initialParameters)
resampler = itk.ResampleImageFilter.New(Input=t2_image,
Transform=initialTransform,
UseReferenceImage=True,
ReferenceImage=t1_image)
resampler.SetDefaultPixelValue(0)
resampler.Update()
plt.subplot(131)
plt.imshow(itk.GetArrayFromImage(resampler.GetOutput()))
plt.subplot(132)
plt.imshow(itk.GetArrayFromImage(t1_image))
It's quite off. Maybe I should adjust the parameters. Don't be afraid to repeat it once more using the aligned image.
registration.SetMovingInitialTransform(initialTransform)
Dimension = t1_image.GetImageDimension()
identityTransform = TransformTypeR.New()
identityTransform.SetIdentity()
registration.SetFixedInitialTransform(identityTransform)
registration.SetNumberOfLevels(1)
registration.SetSmoothingSigmasPerLevel([0])
registration.SetShrinkFactorsPerLevel([1])
try:
registration.Update()
except RuntimeError as exc:
print("Exception ocurred:\n", exc, "\n\n")
finally:
print("Registration finished")
CompositeTransformType = itk.CompositeTransform[itk.D, Dimension]
outputCompositeTransform = CompositeTransformType.New()
outputCompositeTransform.AddTransform(registration.GetTransform())
outputCompositeTransform.AddTransform(initialTransform)
resamplerResult = itk.ResampleImageFilter.New(Input=t2_image,
Transform=outputCompositeTransform,
UseReferenceImage=True,
ReferenceImage=t1_image)
resamplerResult.SetDefaultPixelValue(0)
resamplerResult.Update()
# final position of T2 image when resampled into T1 image's pixel grid
plt.imshow(itk.GetArrayFromImage(resamplerResult.GetOutput()))
Outputs alignment information.
transform = registration.GetTransform()
finalParameters = transform.GetParameters()
rotationInRadians = finalParameters.GetElement(0)
translationAlongX = finalParameters.GetElement(1)
translationAlongY = finalParameters.GetElement(2)
rotationInDegrees = rotationInRadians * 360 / 3.141592 # radian to degree
numberOfIterations = optimizer.GetCurrentIteration()
bestValue = optimizer.GetValue()
print("Result = ")
print(" Rotation degr = " + str(rotationInDegrees))
print(" Translation X = " + str(translationAlongX))
print(" Translation Y = " + str(translationAlongY))
print(" Rotation rad. = " + str(rotationInRadians))
print(" Iterations = " + str(numberOfIterations))
print(" Metric value = " + str(bestValue))
Let's look at the difference.
difference = itk.AddImageFilter.New(Input1=t1_image,
Input2=resamplerResult)
difference.Update()
plt.imshow(itk.GetArrayFromImage(difference.GetOutput()))
that's all.
Reference
- https://github.com/KitwareMedical/2019-03-13-KRSCourseInBiomedicalImageAnalysisAndVisualization
- https://itk.org/ITKSoftwareGuide/html/Book2/ITKSoftwareGuide-Book2ch4.html#x37-1770004
Recommended Posts