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))

filtering1.png

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))

filtering2.png

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()))

filtering3.png

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

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.

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))

せぐ1.png

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))

セグ2.png

#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))

せぐ3.png

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()

segu1.png

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))

せぐ4.png

#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))

せぐ5.png

#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))

せぐ66.png

#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))

segu7.png

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))

segu8.png

From here, try the Shape Detection Level Set. The processing speed is improved, but the curvature is not good. The pipeline looks like this: shape-detection-pipeline.png

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))

segu9.png

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))

segu11.png

#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))

segu12.png

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)

segu14.png segu15.png segu16.png segu17.png

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))

regi1.png

#Image that is out of position
movingImage = itk.imread('data/BrainProtonDensitySliceShifted13x17y.png', PixelType)
plt.imshow(itk.array_view_from_image(movingImage))

regi2.png

#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))

regi3.png

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()))

regi4.png

#Misalignment before and after registration
resampler.SetTransform(identityTransform)
difference.Update()
plt.imshow(itk.array_view_from_image(difference.GetOutput()))

regi5.png

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))

regi6.png

#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))

regi7.png

#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))

regi8.png

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()))

regi9.png

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()))

regi10.png

that's all.

Reference

Recommended Posts

Elementary ITK usage learned in Python
Refactoring Learned in Python (Basic)
Python classes learned in chemoinformatics
What I learned in Python
Character code learned in Python
Python functions learned in chemoinformatics
I learned about processes in Python
Basic Linear Algebra Learned in Python (Part 1)
When looking at memory usage in Python 3
Non-logical operator usage of or in python
Quadtree in Python --2
Python in optimization
CURL in python
Geocoding in python
SendKeys in Python
Meta-analysis in Python
Unittest in python
Epoch in Python
Discord in Python
Sudoku in Python
DCI in Python
quicksort in python
nCr in python
N-Gram in Python
Programming in python
Plink in Python
Constant in python
FizzBuzz in Python
Sqlite in python
StepAIC in Python
N-gram in python
LINE-Bot [0] in Python
Csv in python
Disassemble in Python
Reflection in Python
Constant in python
nCr in Python.
format in python
Scons in Python3
Puyo Puyo in python
python in virtualenv
PPAP in Python
Quad-tree in Python
Reflection in Python
Chemistry in Python
Hashable in python
DirectLiNGAM in Python
LiNGAM in Python
Flatten in python
flatten in python
Python variables and data types learned in chemoinformatics
Statistical test grade 2 probability distribution learned in Python ②
One liner to 100% CPU1 core usage in Python
From file to graph drawing in Python. Elementary elementary
Survival time analysis learned in Python 2 -Kaplan-Meier estimator
Statistical test grade 2 probability distribution learned in Python ①
TensorFlow: Run data learned in Python on Android
Sorted list in Python
Daily AtCoder # 36 in Python
Clustering text in Python
Daily AtCoder # 2 in Python