SimpleITK 图像配准
在网上找的资源，效果不佳，等清楚了函数和原理再细改，调试效果。
1 # -*- coding : UTF-8 -*-
2 # @file : regist.py
3 # @Time : 2021-11-12 17:00
4 # @Author : wmz
5
6 import SimpleITK as sitk
7
8 # Utility method that either downloads data from the MIDAS repository or
9 # if already downloaded returns the file name for reading from disk (cached data).
10 # %run update_path_to_download_script
11 # from downloaddata import fetch_data as fdata
12
13 # Always write output to a separate directory, we don't want to pollute the source directory.
14 import os
15 OUTPUT_DIR = 'Output'
16
17 import matplotlib.pyplot as plt
18 # % matplotlib
19 # inline
20
21 from ipywidgets import interact, fixed
22 from IPython.display import clear_output
23
24
25 # Callback invoked by the interact IPython method for scrolling through the image stacks of
26 # the two images (moving and fixed).
27 def display_images(fixed_image_z, moving_image_z, fixed_npa, moving_npa):
28 # Create a figure with two subplots and the specified size.
29 plt.subplots(1, 2, figsize=(10, 8))
30
31 # Draw the fixed image in the first subplot.
32 plt.subplot(1, 2, 1)
33 plt.imshow(fixed_npa[fixed_image_z, :, :], cmap=plt.cm.Greys_r);
34 plt.title('fixed image')
35 plt.axis('off')
36
37 # Draw the moving image in the second subplot.
38 plt.subplot(1, 2, 2)
39 plt.imshow(moving_npa[moving_image_z, :, :], cmap=plt.cm.Greys_r);
40 plt.title('moving image')
41 plt.axis('off')
42
43 plt.show()
44
45
46 # Callback invoked by the IPython interact method for scrolling and modifying the alpha blending
47 # of an image stack of two images that occupy the same physical space.
48 def display_images_with_alpha(image_z, alpha, fixed, moving):
49 img = (1.0 - alpha) * fixed[:, :, image_z] + alpha * moving[:, :, image_z]
50 plt.imshow(sitk.GetArrayViewFromImage(img), cmap=plt.cm.Greys_r);
51 plt.axis('off')
52 plt.show()
53
54
55 # Callback invoked when the StartEvent happens, sets up our new data.
56 def start_plot():
57 global metric_values, multires_iterations
58
59 metric_values = []
60 multires_iterations = []
61
62
63 # Callback invoked when the EndEvent happens, do cleanup of data and figure.
64 def end_plot():
65 global metric_values, multires_iterations
66
67 del metric_values
68 del multires_iterations
69 # Close figure, we don't want to get a duplicate of the plot latter on.
70 plt.close()
71
72
73 # Callback invoked when the IterationEvent happens, update our data and display new figure.
74 def plot_values(registration_method):
75 global metric_values, multires_iterations
76
77 metric_values.append(registration_method.GetMetricValue())
78 # Clear the output area (wait=True, to reduce flickering), and plot current data
79 clear_output(wait=True)
80 # Plot the similarity metric values
81 plt.plot(metric_values, 'r')
82 plt.plot(multires_iterations, [metric_values[index] for index in multires_iterations], 'b*')
83 plt.xlabel('Iteration Number', fontsize=12)
84 plt.ylabel('Metric Value', fontsize=12)
85 plt.show()
86
87
88 # Callback invoked when the sitkMultiResolutionIterationEvent happens, update the index into the
89 # metric_values list.
90 def update_multires_iterations():
91 global metric_values, multires_iterations
92 multires_iterations.append(len(metric_values))
93
94
95 fixed_image = sitk.ReadImage(r"E:\Data\right_knee\01-YCQ-cl-L.nrrd", sitk.sitkFloat32)
96 moving_image = sitk.ReadImage(r"E:\Data\right_knee\02-LAXI-r-cl lv.nrrd", sitk.sitkFloat32)
97
98 interact(display_images, fixed_image_z=(0,fixed_image.GetSize()[2]-1), moving_image_z=(0,moving_image.GetSize()[2]-1), fixed_npa = fixed(sitk.GetArrayViewFromImage(fixed_image)), moving_npa=fixed(sitk.GetArrayViewFromImage(moving_image)));
99
100
101 initial_transform = sitk.CenteredTransformInitializer(fixed_image,
102 moving_image,
103 sitk.Euler3DTransform(),
104 sitk.CenteredTransformInitializerFilter.GEOMETRY)
105
106 moving_resampled = sitk.Resample(moving_image, fixed_image, initial_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID())
107
108 interact(display_images_with_alpha, image_z=(0,fixed_image.GetSize()[2]), alpha=(0.0,1.0,0.05), fixed = fixed(fixed_image), moving=fixed(moving_resampled));
109
110 registration_method = sitk.ImageRegistrationMethod()
111
112 # Similarity metric settings.
113 registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
114 registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
115 registration_method.SetMetricSamplingPercentage(0.01)
116
117 registration_method.SetInterpolator(sitk.sitkLinear)
118
119 # Optimizer settings.
120 registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, convergenceMinimumValue=1e-6, convergenceWindowSize=10)
121 registration_method.SetOptimizerScalesFromPhysicalShift()
122
123 # Setup for the multi-resolution framework.
124 registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])
125 registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,0])
126 registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
127
128 # Don't optimize in-place, we would possibly like to run this cell multiple times.
129 registration_method.SetInitialTransform(initial_transform, inPlace=False)
130
131 # Connect all of the observers so that we can perform plotting during registration.
132 registration_method.AddCommand(sitk.sitkStartEvent, start_plot)
133 registration_method.AddCommand(sitk.sitkEndEvent, end_plot)
134 registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, update_multires_iterations)
135 registration_method.AddCommand(sitk.sitkIterationEvent, lambda: plot_values(registration_method))
136
137 final_transform = registration_method.Execute(sitk.Cast(fixed_image, sitk.sitkFloat32),
138 sitk.Cast(moving_image, sitk.sitkFloat32))
139
140 print('Final metric value: {0}'.format(registration_method.GetMetricValue()))
141 print('Optimizer\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))
142
143 moving_resampled = sitk.Resample(moving_image, fixed_image, final_transform, sitk.sitkLinear, 0.0, moving_image.GetPixelID())
144
145 interact(display_images_with_alpha, image_z=(0,fixed_image.GetSize()[2]), alpha=(0.0,1.0,0.05), fixed = fixed(fixed_image), moving=fixed(moving_resampled));
146
147 sitk.WriteImage(moving_resampled, os.path.join(OUTPUT_DIR, 'RIRE_training_001_mr_T1_resampled.mha'))
148 sitk.WriteTransform(final_transform, os.path.join(OUTPUT_DIR, 'RIRE_training_001_CT_2_mr_T1.tfm'))
