ITK  5.4.0 Insight Toolkit
SphinxExamples/src/Registration/Metricsv4/RegisterTwoPointSets/Code.py
1 #!/usr/bin/env python3
2
3 # =========================================================================
4 #
6 #
8 # you may not use this file except in compliance with the License.
9 # You may obtain a copy of the License at
10 #
12 #
13 # Unless required by applicable law or agreed to in writing, software
15 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
16 # See the License for the specific language governing permissions and
17 # limitations under the License.
18 #
19 # =========================================================================*/
20
21 # Adapted from ITK itkJensenHavrdaCharvatTsallisPointSetMetricRegistrationTest.cxx
22
23 import sys
24 from math import pi, sin, cos
25
26 import itk
27
28
29 # Generate two circles with a small offset
30 def make_circles(l_dimension: int = 2):
31  PointSetType = itk.PointSet[itk.F, l_dimension]
32
34  offset = [2.0] * l_dimension
35
36  fixed_points = PointSetType.New()
37  moving_points = PointSetType.New()
38  fixed_points.Initialize()
39  moving_points.Initialize()
40
41  step = 0.1
42  for count in range(0, int(2 * pi / step) + 1):
43
44  theta = count * step
45
46  fixed_point = list()
48  for dim in range(1, l_dimension):
50  fixed_points.SetPoint(count, fixed_point)
51
52  moving_point = [fixed_point[dim] + offset[dim] for dim in range(0, l_dimension)]
53  moving_points.SetPoint(count, moving_point)
54
55  return fixed_points, moving_points
56
57
58 def test_registration(l_dimension: int = 2):
59  # Define test parameters
60
61  num_iterations = 10
62
63  passed = True
64  tolerance = 0.05
65
66  # Define types
67  PointSetType = itk.PointSet[itk.F, l_dimension]
68  AffineTransformType = itk.AffineTransform[itk.D, l_dimension]
70  PointSetType
71  ]
73  PointSetMetricType
74  ]
76
77  # Make point sets
78  fixed_set, moving_set = make_circles(l_dimension)
79
80  transform = AffineTransformType.New()
81  transform.SetIdentity()
82
83  metric = PointSetMetricType.New(
84  FixedPointSet=fixed_set,
85  MovingPointSet=moving_set,
86  PointSetSigma=1.0,
87  KernelSigma=10.0,
88  UseAnisotropicCovariances=False,
89  CovarianceKNeighborhood=5,
90  EvaluationKNeighborhood=10,
91  MovingTransform=transform,
92  Alpha=1.1,
93  )
94  metric.Initialize()
95
96  shift_scale_estimator = ShiftScalesType.New(
97  Metric=metric, VirtualDomainPointSet=metric.GetVirtualTransformedPointSet()
98  )
99
100  optimizer = OptimizerType.New(
101  Metric=metric,
102  NumberOfIterations=num_iterations,
103  ScalesEstimator=shift_scale_estimator,
104  MaximumStepSizeInPhysicalUnits=3.0,
105  MinimumConvergenceValue=0.0,
106  ConvergenceWindowSize=10,
107  )
108
109  def print_iteration():
110  print(
111  f"It: {optimizer.GetCurrentIteration()}"
112  f" metric value: {optimizer.GetCurrentMetricValue():.6f} "
113  )
114
116
117  # Run optimization to align the point sets
118  optimizer.StartOptimization()
119
120  print(f"Number of iterations: {num_iterations}")
121  print(f"Moving-source final value: {optimizer.GetCurrentMetricValue()}")
122  print(f"Moving-source final position: {list(optimizer.GetCurrentPosition())}")
123  print(f"Optimizer scales: {list(optimizer.GetScales())}")
124  print(f"Optimizer learning rate: {optimizer.GetLearningRate()}")
125
126  # applying the resultant transform to moving points and verify result
127  print("Fixed\tMoving\tMovingTransformed\tFixedTransformed\tDiff")
128
129  moving_inverse = metric.GetMovingTransform().GetInverseTransform()
130  fixed_inverse = metric.GetFixedTransform().GetInverseTransform()
131
132  def print_point(vals: list) -> str:
133  return f'[{",".join(f"{x:.4f}" for x in vals)}]'
134
135  for n in range(0, metric.GetNumberOfComponents()):
136  transformed_moving_point = moving_inverse.TransformPoint(moving_set.GetPoint(n))
137  transformed_fixed_point = fixed_inverse.TransformPoint(fixed_set.GetPoint(n))
138
139  difference = [
140  transformed_moving_point[dim] - transformed_fixed_point[dim]
141  for dim in range(0, l_dimension)
142  ]
143
144  print(
145  f"{print_point(fixed_set.GetPoint(n))}"
146  f"\t{print_point(moving_set.GetPoint(n))}"
147  f"\t{print_point(transformed_moving_point)}"
148  f"\t{print_point(transformed_fixed_point)}"
149  f"\t{print_point(difference)}"
150  )
151
152  if any(abs(difference[dim]) > tolerance for dim in range(0, l_dimension)):
153  passed = False
154
155  if not passed:
156  raise Exception("Transform outside of allowable tolerance")
157  else:
158  print("Transform is within allowable tolerance.")
159
160
161 if __name__ == "__main__":
162  if len(sys.argv) == 2:
163  dimension = int(sys.argv[1])
164  test_registration(dimension)
165  else:
166  test_registration()
itk::PointSet
A superclass of the N-dimensional mesh structure; supports point (geometric coordinate and attribute)...
Definition: itkPointSet.h:82
itk::AffineTransform
Definition: itkAffineTransform.h:101
itk::Math::abs
bool abs(bool x)
Definition: itkMath.h:856