Hello,
I want to compare the performance of different optimization algorithms while registering two images. For this purpose, I want to know what are the valid input parameters (step size, function tolerance, number of iteration etc), that each optimizer, available in SimpleITK, can accept. I am using SimpleITK with Python. Where can I get this information?
Hello @debapriya,
This information is found in the ImageRegistrationMethod documentation SetOptimizerAs* methods.
Thank you @zivy . In order to test the Amoeba optimizer, I set the following parameters.
registration_method.SetOptimizerAsAmoeba(simplexDelta=1.00, numberOfIterations=10, withRestarts=False)
However, the optimization continued for 400 iterations, and finally stopped with the following error
Please tell me what is going wrong.
Can this be because of multi-resolution registration? How can I stop this optimization process after a fixed number of iterations?
What is an acceptable value of SimplexDelta?
Hello @debapriya,
The number of iterations is for each run of the optimizer. When using the multi-resolution approach the optimizer is re-run for each level. To see this effect, run the optimizer using a single level and set the number of iterations to a small number (e.g. 10).
The SimplexDelta value depends on the parameter space terrain and how far the initial parameter values are from the final ones (the initial simplex is defined using the initial parameter values and the SimplexDelta).
What is the default value of simplexDelta?
What value of simplexDelta should I set, if I do not know the shape of the parameter space terrain?
Hello @debapriya,
As you can see from the documentation, simplexDelta does not have a default value. This requires domain knowledge so up to the user. The simplex vertices are ([x0[0], x0[1], ..., x0[i]+simplexDelta, ..., x0[d-1]]) where x0 are the initial parameter values. You may need to perform a hyper-parameter search step (try different values for the simplexDelta to identify the one that works best for your case).
A few simplexDelta values, and the corresponding Mean TRE and Max TRE values obtained, are listed below.
However, each time, the code stopped with the following exception
Unexpected exception formatting exception. Falling back to standard exception
Traceback (most recent call last):
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\interactiveshell.py", line 3548, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "C:\Users\debap\AppData\Local\Temp\ipykernel_19716\1123557497.py", line 1, in <module>
get_ipython().run_cell_magic('timeit', '-r1 -n1', '\nglobal tx\n\n# Select the fixed and moving images, valid entries are in [0,9].\n# fixed_image_index = 0\n# moving_image_index = 7\n\n# fixed_image_index = 1\n# moving_image_index = 8\n\n# fixed_image_index = 2\n# moving_image_index = 9\n\n# fixed_image_index = 0\n# moving_image_index = 9\n\n# fixed_image_index = 0\n# moving_image_index = 8\n\n# fixed_image_index = 3\n# moving_image_index = 6\n\n# fixed_image_index = 1\n# moving_image_index = 6\n\n# fixed_image_index = 1\n# moving_image_index = 7\n\n# fixed_image_index = 1\n# moving_image_index = 5\n\nfixed_image_index = 5\nmoving_image_index = 9\n\n\n\ntx = bspline_intra_modal_registration(\n fixed_image=images[fixed_image_index],\n moving_image=images[moving_image_index],\n fixed_image_mask=(masks[fixed_image_index] == lung_label),\n fixed_points=points[fixed_image_index],\n moving_points=points[moving_image_index],\n)\n(\n initial_errors_mean,\n initial_errors_std,\n _,\n initial_errors_max,\n initial_errors,\n) = ru.registration_errors(\n sitk.Euler3DTransform(), points[fixed_image_index], points[moving_image_index]\n)\n(\n final_errors_mean,\n final_errors_std,\n _,\n final_errors_max,\n final_errors,\n) = ru.registration_errors(tx, points[fixed_image_index], points[moving_image_index])\n\nplt.hist(initial_errors, bins=500, alpha=0.5, label="before registration", color="blue")\nplt.hist(final_errors, bins=500, alpha=0.5, label="after registration", color="green")\nplt.legend()\nplt.title("TRE histogram")\nprint(\n f"Initial alignment errors in millimeters, mean(std): {initial_errors_mean:.2f}({initial_errors_std:.2f}), max: {initial_errors_max:.2f}"\n)\nprint(\n f"Final alignment errors in millimeters, mean(std): {final_errors_mean:.2f}({final_errors_std:.2f}), max: {final_errors_max:.2f}"\n)\n\n\n# print(f"{initial_errors_med:.2f}")\n# print(f"{final_errors_med:.2f}")\n')
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\interactiveshell.py", line 2515, in run_cell_magic
result = fn(*args, **kwargs)
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\magics\execution.py", line 1189, in timeit
all_runs = timer.repeat(repeat, number)
File "C:\Users\debap\AppData\Local\Programs\Python\Python310\lib\timeit.py", line 206, in repeat
t = self.timeit(number)
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\magics\execution.py", line 173, in timeit
timing = self.inner(it, self.timer)
File "<magic-timeit>", line 36, in inner
File "C:\Users\debap\AppData\Local\Temp\ipykernel_19716\2841005318.py", line 78, in bspline_intra_modal_registration
return registration_method.Execute(fixed_image, moving_image)
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\SimpleITK\SimpleITK.py", line 10863, in Execute
val = _SimpleITK.ImageRegistrationMethod_Execute(self, fixed, moving)
RuntimeError: Exception thrown in SimpleITK ImageRegistrationMethod_Execute: C:\Users\debap\Desktop\ITK\Modules\Numerics\Optimizersv4\src\itkAmoebaOptimizerv4.cxx:268:
ITK ERROR: AmoebaOptimizerv4(000002A650BF1D80): cost function and simplex delta dimensions mismatch
During handling of the above exception, another exception occurred:
Traceback (most recent call last):
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\interactiveshell.py", line 2142, in showtraceback
stb = self.InteractiveTB.structured_traceback(
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\ultratb.py", line 1435, in structured_traceback
return FormattedTB.structured_traceback(
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\ultratb.py", line 1326, in structured_traceback
return VerboseTB.structured_traceback(
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\ultratb.py", line 1173, in structured_traceback
formatted_exception = self.format_exception_as_a_whole(etype, evalue, etb, number_of_lines_of_context,
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\ultratb.py", line 1063, in format_exception_as_a_whole
self.get_records(etb, number_of_lines_of_context, tb_offset) if etb else []
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\ultratb.py", line 1155, in get_records
FrameInfo(
File "C:\Users\debap\Desktop\Python_Vir_Env\Python\lib\site-packages\IPython\core\ultratb.py", line 780, in __init__
ix = inspect.getsourcelines(frame)
File "C:\Users\debap\AppData\Local\Programs\Python\Python310\lib\inspect.py", line 1129, in getsourcelines
lines, lnum = findsource(object)
File "C:\Users\debap\AppData\Local\Programs\Python\Python310\lib\inspect.py", line 958, in findsource
raise OSError('could not get source code')
OSError: could not get source code
Can you please share an example code which uses the Amoeba optimizer?
My registration function is copied below. Is it possible that Amoeba optimizer is not suitable for non-rigid (BSpline) registration?
def bspline_intra_modal_registration(
fixed_image,
moving_image,
fixed_image_mask=None,
fixed_points=None,
moving_points=None,
):
registration_method = sitk.ImageRegistrationMethod()
# Determine the number of BSpline control points using the physical spacing we want for the control grid.
grid_physical_spacing = [50.0, 50.0, 50.0] # A control point every 50mm
image_physical_size = [
size * spacing
for size, spacing in zip(fixed_image.GetSize(), fixed_image.GetSpacing())
]
mesh_size = [
int(image_size / grid_spacing + 0.5)
for image_size, grid_spacing in zip(image_physical_size, grid_physical_spacing)
]
mesh_size = [int(sz / 4 + 0.5) for sz in mesh_size]
initial_transform = sitk.BSplineTransformInitializer(
image1=fixed_image, transformDomainMeshSize=mesh_size, order=3
)
registration_method.SetInitialTransformAsBSpline(
initial_transform, inPlace=True, scaleFactors=[1, 2, 4]
)
registration_method.SetMetricAsMeanSquares()
# Settings for metric sampling, usage of a mask is optional. When given a mask the sample points will be
# generated inside that region. Also, this implicitly speeds things up as the mask is smaller than the
# whole image.
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(0.01)
if fixed_image_mask:
registration_method.SetMetricFixedMask(fixed_image_mask)
# Multi-resolution framework.
registration_method.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2, 1])
registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1, 0])
registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
registration_method.SetInterpolator(sitk.sitkLinear)
registration_method.SetOptimizerAsAmoeba(numberOfIterations=10, simplexDelta=True)
# If corresponding points in the fixed and moving image are given then we display the similarity metric
# and the TRE during the registration.
if fixed_points and moving_points:
registration_method.AddCommand(
sitk.sitkStartEvent, rc.metric_and_reference_start_plot
)
registration_method.AddCommand(
sitk.sitkEndEvent, rc.metric_and_reference_end_plot
)
registration_method.AddCommand(
sitk.sitkIterationEvent,
lambda: rc.metric_and_reference_plot_values(
registration_method, fixed_points, moving_points
),
)
return registration_method.Execute(fixed_image, moving_image)
Hello @debapriya,
Your observation is generally correct. The BSpline transform usually has a large number of parameters (dim*mesh_size). In this case, the most appropriate optimizer from those supported by ITK/SimpleITK is LBFGSB, Limited memory Broyden Fletcher Goldfarb Shannon minimization with simple Bounds (if no bounds are given it is unbounded). For more about LBFGS and LBFGSB see wikipedia.
The downhill simplex/Nelder Mead/amoeba is usually appropriate for low dimensional parameter spaces (e.g. rigid transformation) and similarity/loss functions that are not differentiable (it doesn’t require the derivative as other algorithms do).