Minimal registration example ¶

This tutorial will walk you through a very minimalistic example for image registration via mermaid.

Contents

Minimal registration example

Introduction ¶

There are essentially three different ways of how to do registration with pytorch. Below are these options in increasing level of difficulty, but in order of increasing flexibility.

Use the simple registration interface (see example_simple_interface.py in the demos directory.
Manually instantiating the mermaid optimizer and setting up the registration parameters (which we are going to do here).
Instantiating the mermaid registration models and writing your own optimizer for them.

Importing modules ¶

We start by importing some popular modules. The used mermaid modules are:

mermaid.multiscale_optimizer provides single scale and multi-scale optimizing functionality for all the registration models.
mermaid.example_generation allows to generate simply synthetic data and some real image pairs to test the registration algorithms
mermaid.module_parameters allows generating mermaid parameter structures which are used to keep track of all the parameters

import torch
from mermaid.data_wrapper import AdaptVal
import numpy as np

import mermaid.example_generation as eg
import mermaid.module_parameters as pars
import mermaid.multiscale_optimizer as MO

Specifying registration model properties ¶

Next we specify some settings which will need to be specified when setting up the registration model. These are the meanings of the following variables:

use_map: Most image registration models in mermaid (and all the ones that are based on integration through velocity fields) support solutions by directly solving the associated equations on the images or alternatively via evolution of the transformation map (if use_map=True). In general, a map-based solution is preferred as it avoids numerical discretization issues when, for example, directly advecting an image and because it directly results in the desired transformation map. Note that within mermaid a map phi denotes a map in target space (i.e., to map from source to target space) and its inverse map is a map defined in the coordinates of the source image (to map from the target space to the source space).
model_name: Simply the name of the desired registration model. For example, lddmm_shooting is an LDDMM implementation which optimizes over the initial vector momentum and operates either on maps or directly on the image (depening on use_map).
map_low_res_factor: Especially when computing in 3D memory sometimes becomes an issue. Hence, for map-based solutions mermaid supports computing solutions of the evolution equations at a lower resolution. map_low_res_factor specifies how low this resoultion should be. None or 1.0 is the original resolution. 0.5, for example, uses half the resolution in each direction. In any case, the similarity measure is always evaluated at the original resolution via upsampling of the map.

use_map = True
model_name = 'lddmm_shooting'
map_low_res_factor = 0.5

if use_map:
    model_name = model_name + '_map'
else:
    model_name = model_name + '_image'

mermaid support mostly sgd or lbfgs_ls (an LBFGS method with linesearch). Other optimizers can be used, but this is not well supported at the moment. As each of these optimizers iterates the number of desired iterations needs to be specified. The optimizer also suppport visualizing intermediate output at given iteration intervals.

optimizer_name = 'sgd'
nr_of_iterations = 50
visualize = True
visualize_step = 10

Creating the parameter structure ¶

AS typical in mermaid, we create a parameters structure (params) that will keep track of all the parameters used during a run.

# keep track of general parameters
params = pars.ParameterDict()

Creating an example ¶

Now we are ready to create some example data. We create squares for source and target images in two dimensions.

dim = 2
example_img_len = 64

szEx = np.tile( example_img_len, dim )         # size of the desired images: (sz)^dim
I0,I1,spacing= eg.CreateSquares(dim).create_image_pair(szEx,params) # create a default image size with two sample squares
sz = np.array(I0.shape)

Out:

Creating new category: root.square_example_images
Using default value = 10 for key = len_s of category = root.square_example_images
Using default value = 16 for key = len_l of category = root.square_example_images

Moving from numpy to pytorch ¶

The example generation produces numpy arrays. As mermaid uses pytorch these need to be converted to pytorch arrays. Also, we support running mermaid on the CPU and the GPU. The convenience function AdaptVal takes care of this by moving an array either to the GPU or leaving it on the CPU.

# create the source and target image as pyTorch variables
ISource = AdaptVal(torch.from_numpy(I0.copy()))
ITarget = AdaptVal(torch.from_numpy(I1))

Instantiating the optimizer ¶

Now we are ready to instantiate the optimizer with the parameters defined above. We instantiate a single-scale optimizer here, but instantiating a multi-scale optimizer proceeds similarly. We then start the optimization (via so.optimizer()).

so = MO.SingleScaleRegistrationOptimizer(sz,spacing,use_map,map_low_res_factor,params)
so.set_model(model_name)
so.set_optimizer_by_name( optimizer_name )
so.set_visualization( visualize )
so.set_visualize_step( visualize_step )

so.set_number_of_iterations(nr_of_iterations)

so.set_source_image(ISource)
so.set_target_image(ITarget)

# and now do the optimization
so.optimize()

../_images/sphx_glr_example_minimal_registration_without_simple_interface_001.png

Out:

Creating new category: root.optimizer
Creating new category: root.model
Creating new category: root.model.deformation
Creating new category: root.model.registration_model
Creating key = use_map; category = root.model.deformation; value = True
Creating key = map_low_res_factor; category = root.model.deformation; value = 0.5
Using default value = False for key = compute_similarity_measure_at_low_res of category = root.model.deformation
Creating new category: root.optimizer.single_scale
Using default value = 0.0001 for key = rel_ftol of category = root.optimizer.single_scale
Using default value = 1 for key = spline_order of category = root.model.registration_model
Using default value = none for key = weight_clipping_type of category = root.optimizer
Using default value = 1.0 for key = weight_clipping_value of category = root.optimizer
Creating new category: root.optimizer.gradient_clipping
Using default value = True for key = clip_display of category = root.optimizer.gradient_clipping
Using default value = False for key = clip_individual_gradient of category = root.optimizer.gradient_clipping
Using default value = 1.0158730158730158 for key = clip_individual_gradient_value of category = root.optimizer.gradient_clipping
Using default value = True for key = clip_shared_gradient of category = root.optimizer.gradient_clipping
Using default value = 1.0 for key = clip_shared_gradient_value of category = root.optimizer.gradient_clipping
Creating key = type; category = root.model.registration_model; value = lddmm_shooting_map
Overwriting key = type; category = root.model.registration_model; value =  lddmm_shooting_map -> lddmm_shooting_map
Creating new category: root.model.registration_model.env
Using default value = False for key = get_momentum_from_external_network of category = root.model.registration_model.env
Using map-based lddmm_shooting_map model
Using default value = True for key = use_CFL_clamping of category = root.model.registration_model
Using default value = True for key = use_odeint of category = root.model.registration_model.env
Using default value = False for key = use_ode_tuple of category = root.model.registration_model.env
Creating new category: root.model.registration_model.forward_model
Creating new category: root.model.registration_model.forward_model.smoother
Using default value = multiGaussian for key = type of category = root.model.registration_model.forward_model.smoother
Using default value = [0.05, 0.1, 0.15, 0.2, 0.25] for key = multi_gaussian_stds of category = root.model.registration_model.forward_model.smoother
Using default value = [0.06666666666666667, 0.13333333333333333, 0.19999999999999998, 0.26666666666666666, 0.3333333333333333] for key = multi_gaussian_weights of category = root.model.registration_model.forward_model.smoother
the param of smoother is <mermaid.smoother_factory.MultiGaussianFourierSmoother object at 0x7f30e2d376d8>
Creating new category: root.model.registration_model.shooting_vector_momentum
Using default value = False for key = use_velocity_mask_on_boundary of category = root.model.registration_model.shooting_vector_momentum
works in mermaid iter mode
Using default value = 1.0 for key = reg_factor of category = root.model.registration_model.env
Creating new category: root.model.registration_model.loss
Using default value = False for key = display_max_displacement of category = root.model.registration_model.loss
Using default value = False for key = limit_displacement of category = root.model.registration_model.loss
Using default value = 0.05 for key = max_displacement of category = root.model.registration_model.loss
LDDMMShootingVectorMomentumMapNet(
  (integrator): ODEWrapBlock()
)
Creating key = name; category = root.optimizer; value = sgd
Creating key = nr_of_iterations; category = root.optimizer.single_scale; value = 50
Creating new category: root.optimizer.sgd
Creating new category: root.optimizer.sgd.individual
Using default value = 0.01 for key = lr of category = root.optimizer.sgd.individual
Using default value = 0.9 for key = momentum of category = root.optimizer.sgd.individual
Using default value = 0.0 for key = dampening of category = root.optimizer.sgd.individual
Using default value = 0.0 for key = weight_decay of category = root.optimizer.sgd.individual
Using default value = True for key = nesterov of category = root.optimizer.sgd.individual
Creating new category: root.optimizer.sgd.shared
Using default value = 0.01 for key = lr of category = root.optimizer.sgd.shared
Using default value = 0.9 for key = momentum of category = root.optimizer.sgd.shared
Using default value = 0.0 for key = dampening of category = root.optimizer.sgd.shared
Using default value = 0.0 for key = weight_decay of category = root.optimizer.sgd.shared
Using default value = True for key = nesterov of category = root.optimizer.sgd.shared
Using default value = True for key = use_step_size_scheduler of category = root.optimizer
Creating new category: root.optimizer.scheduler
Using default value = True for key = verbose of category = root.optimizer.scheduler
Using default value = 0.5 for key = factor of category = root.optimizer.scheduler
Using default value = 10 for key = patience of category = root.optimizer.scheduler
Using default value = rk4 for key = solver of category = root.model.registration_model.forward_model
Using default value = True for key = adjoin_on of category = root.model.registration_model.forward_model
Using default value = 1e-05 for key = rtol of category = root.model.registration_model.forward_model
Using default value = 1e-05 for key = atol of category = root.model.registration_model.forward_model
Using default value = 20 for key = number_of_time_steps of category = root.model.registration_model.forward_model
Creating new category: root.model.registration_model.similarity_measure
Using default value = ssd for key = type of category = root.model.registration_model.similarity_measure
Using ssd similarity measure
Using default value = 0.1 for key = sigma of category = root.model.registration_model.similarity_measure
[31m    0-Tot: E=015.7218 | simE=015.7218 | regE=000.0000 | optParE=000.0000 | relF=  n/a    | [0m
[34m    0-Img: E=015.7218 | simE=015.7218 | regE=000.0000 |[0m
[31m    1-Tot: E=015.4231 | simE=015.4231 | regE=000.0000 | optParE=000.0000 | relF=000.0182 | [0m
[34m    1-Img: E=015.4231 | simE=015.4231 | regE=000.0000 |[0m
[31m    2-Tot: E=015.0197 | simE=015.0197 | regE=000.0000 | optParE=000.0000 | relF=000.0252 | [0m
[34m    2-Img: E=015.0197 | simE=015.0197 | regE=000.0000 |[0m
[31m    3-Tot: E=014.6717 | simE=014.6716 | regE=000.0000 | optParE=000.0000 | relF=000.0222 | [0m
[34m    3-Img: E=014.6717 | simE=014.6716 | regE=000.0000 |[0m
[31m    4-Tot: E=014.4240 | simE=014.4239 | regE=000.0001 | optParE=000.0000 | relF=000.0161 | [0m
[34m    4-Img: E=014.4240 | simE=014.4239 | regE=000.0001 |[0m
[31m    5-Tot: E=014.2155 | simE=014.2153 | regE=000.0002 | optParE=000.0000 | relF=000.0137 | [0m
[34m    5-Img: E=014.2155 | simE=014.2153 | regE=000.0002 |[0m
[31m    6-Tot: E=013.9136 | simE=013.9133 | regE=000.0002 | optParE=000.0000 | relF=000.0202 | [0m
[34m    6-Img: E=013.9136 | simE=013.9133 | regE=000.0002 |[0m
[31m    7-Tot: E=013.5529 | simE=013.5526 | regE=000.0004 | optParE=000.0000 | relF=000.0248 | [0m
[34m    7-Img: E=013.5529 | simE=013.5526 | regE=000.0004 |[0m
[31m    8-Tot: E=013.1096 | simE=013.1091 | regE=000.0005 | optParE=000.0000 | relF=000.0314 | [0m
[34m    8-Img: E=013.1096 | simE=013.1091 | regE=000.0005 |[0m
[31m    9-Tot: E=012.6062 | simE=012.6055 | regE=000.0007 | optParE=000.0000 | relF=000.0370 | [0m
[34m    9-Img: E=012.6062 | simE=012.6055 | regE=000.0007 |[0m
[31m   10-Tot: E=012.1933 | simE=012.1923 | regE=000.0009 | optParE=000.0000 | relF=000.0313 | [0m
[34m   10-Img: E=012.1933 | simE=012.1923 | regE=000.0009 |[0m
[31m   11-Tot: E=011.8421 | simE=011.8409 | regE=000.0012 | optParE=000.0000 | relF=000.0273 | [0m
[34m   11-Img: E=011.8421 | simE=011.8409 | regE=000.0012 |[0m
[31m   12-Tot: E=011.2298 | simE=011.2283 | regE=000.0015 | optParE=000.0000 | relF=000.0501 | [0m
[34m   12-Img: E=011.2298 | simE=011.2283 | regE=000.0015 |[0m
[31m   13-Tot: E=010.6379 | simE=010.6360 | regE=000.0019 | optParE=000.0000 | relF=000.0509 | [0m
[34m   13-Img: E=010.6379 | simE=010.6360 | regE=000.0019 |[0m
[31m   14-Tot: E=010.1932 | simE=010.1909 | regE=000.0023 | optParE=000.0000 | relF=000.0397 | [0m
[34m   14-Img: E=010.1932 | simE=010.1909 | regE=000.0023 |[0m
[31m   15-Tot: E=009.5512 | simE=009.5484 | regE=000.0028 | optParE=000.0000 | relF=000.0608 | [0m
[34m   15-Img: E=009.5512 | simE=009.5484 | regE=000.0028 |[0m
[31m   16-Tot: E=008.9230 | simE=008.9197 | regE=000.0033 | optParE=000.0000 | relF=000.0633 | [0m
[34m   16-Img: E=008.9230 | simE=008.9197 | regE=000.0033 |[0m
[31m   17-Tot: E=008.2878 | simE=008.2839 | regE=000.0039 | optParE=000.0000 | relF=000.0684 | [0m
[34m   17-Img: E=008.2878 | simE=008.2839 | regE=000.0039 |[0m
[31m   18-Tot: E=007.6788 | simE=007.6742 | regE=000.0046 | optParE=000.0000 | relF=000.0702 | [0m
[34m   18-Img: E=007.6788 | simE=007.6742 | regE=000.0046 |[0m
[31m   19-Tot: E=007.1398 | simE=007.1344 | regE=000.0053 | optParE=000.0000 | relF=000.0662 | [0m
[34m   19-Img: E=007.1398 | simE=007.1344 | regE=000.0053 |[0m
[31m   20-Tot: E=006.7816 | simE=006.7756 | regE=000.0061 | optParE=000.0000 | relF=000.0460 | [0m
[34m   20-Img: E=006.7816 | simE=006.7756 | regE=000.0061 |[0m
[31m   21-Tot: E=006.2190 | simE=006.2123 | regE=000.0067 | optParE=000.0000 | relF=000.0779 | [0m
[34m   21-Img: E=006.2190 | simE=006.2123 | regE=000.0067 |[0m
[31m   22-Tot: E=005.5830 | simE=005.5756 | regE=000.0073 | optParE=000.0000 | relF=000.0966 | [0m
[34m   22-Img: E=005.5830 | simE=005.5756 | regE=000.0073 |[0m
[31m   23-Tot: E=005.0480 | simE=005.0401 | regE=000.0079 | optParE=000.0000 | relF=000.0885 | [0m
[34m   23-Img: E=005.0480 | simE=005.0401 | regE=000.0079 |[0m
[31m   24-Tot: E=004.5942 | simE=004.5857 | regE=000.0085 | optParE=000.0000 | relF=000.0811 | [0m
[34m   24-Img: E=004.5942 | simE=004.5857 | regE=000.0085 |[0m
[31m   25-Tot: E=004.0458 | simE=004.0367 | regE=000.0091 | optParE=000.0000 | relF=000.1087 | [0m
[34m   25-Img: E=004.0458 | simE=004.0367 | regE=000.0091 |[0m
[31m   26-Tot: E=003.7798 | simE=003.7703 | regE=000.0095 | optParE=000.0000 | relF=000.0556 | [0m
[34m   26-Img: E=003.7798 | simE=003.7703 | regE=000.0095 |[0m
[31m   27-Tot: E=003.5297 | simE=003.5197 | regE=000.0100 | optParE=000.0000 | relF=000.0552 | [0m
[34m   27-Img: E=003.5297 | simE=003.5197 | regE=000.0100 |[0m
[31m   28-Tot: E=003.2569 | simE=003.2465 | regE=000.0104 | optParE=000.0000 | relF=000.0641 | [0m
[34m   28-Img: E=003.2569 | simE=003.2465 | regE=000.0104 |[0m
[31m   29-Tot: E=003.0882 | simE=003.0775 | regE=000.0107 | optParE=000.0000 | relF=000.0413 | [0m
[34m   29-Img: E=003.0882 | simE=003.0775 | regE=000.0107 |[0m
[31m   30-Tot: E=002.7677 | simE=002.7567 | regE=000.0110 | optParE=000.0000 | relF=000.0851 | [0m
[34m   30-Img: E=002.7677 | simE=002.7567 | regE=000.0110 |[0m
[31m   31-Tot: E=002.5098 | simE=002.4985 | regE=000.0113 | optParE=000.0000 | relF=000.0735 | [0m
[34m   31-Img: E=002.5098 | simE=002.4985 | regE=000.0113 |[0m
[31m   32-Tot: E=002.0954 | simE=002.0838 | regE=000.0116 | optParE=000.0000 | relF=000.1339 | [0m
[34m   32-Img: E=002.0954 | simE=002.0838 | regE=000.0116 |[0m
[31m   33-Tot: E=001.7128 | simE=001.7009 | regE=000.0119 | optParE=000.0000 | relF=000.1410 | [0m
[34m   33-Img: E=001.7128 | simE=001.7009 | regE=000.0119 |[0m
[31m   34-Tot: E=001.2203 | simE=001.2081 | regE=000.0122 | optParE=000.0000 | relF=000.2218 | [0m
[34m   34-Img: E=001.2203 | simE=001.2081 | regE=000.0122 |[0m
[31m   35-Tot: E=000.8146 | simE=000.8021 | regE=000.0125 | optParE=000.0000 | relF=000.2236 | [0m
[34m   35-Img: E=000.8146 | simE=000.8021 | regE=000.0125 |[0m
[31m   36-Tot: E=000.5843 | simE=000.5714 | regE=000.0129 | optParE=000.0000 | relF=000.1454 | [0m
[34m   36-Img: E=000.5843 | simE=000.5714 | regE=000.0129 |[0m
[31m   37-Tot: E=000.5364 | simE=000.5232 | regE=000.0132 | optParE=000.0000 | relF=000.0312 | [0m
[34m   37-Img: E=000.5364 | simE=000.5232 | regE=000.0132 |[0m
[31m   38-Tot: E=000.5772 | simE=000.5638 | regE=000.0134 | optParE=000.0000 | relF=000.0259 | [0m
[34m   38-Img: E=000.5772 | simE=000.5638 | regE=000.0134 |[0m
[31m   39-Tot: E=000.6588 | simE=000.6452 | regE=000.0136 | optParE=000.0000 | relF=000.0492 | [0m
[34m   39-Img: E=000.6588 | simE=000.6452 | regE=000.0136 |[0m
[31m   40-Tot: E=000.7195 | simE=000.7057 | regE=000.0137 | optParE=000.0000 | relF=000.0353 | [0m
[34m   40-Img: E=000.7195 | simE=000.7057 | regE=000.0137 |[0m
[31m   41-Tot: E=000.7643 | simE=000.7504 | regE=000.0138 | optParE=000.0000 | relF=000.0254 | [0m
[34m   41-Img: E=000.7643 | simE=000.7504 | regE=000.0138 |[0m
[31m   42-Tot: E=000.7926 | simE=000.7787 | regE=000.0139 | optParE=000.0000 | relF=000.0158 | [0m
[34m   42-Img: E=000.7926 | simE=000.7787 | regE=000.0139 |[0m
[31m   43-Tot: E=000.7945 | simE=000.7805 | regE=000.0139 | optParE=000.0000 | relF=000.0010 | [0m
[34m   43-Img: E=000.7945 | simE=000.7805 | regE=000.0139 |[0m
[31m   44-Tot: E=000.7576 | simE=000.7437 | regE=000.0139 | optParE=000.0000 | relF=000.0210 | [0m
[34m   44-Img: E=000.7576 | simE=000.7437 | regE=000.0139 |[0m
[31m   45-Tot: E=000.6660 | simE=000.6521 | regE=000.0139 | optParE=000.0000 | relF=000.0550 | [0m
[34m   45-Img: E=000.6660 | simE=000.6521 | regE=000.0139 |[0m
[31m   46-Tot: E=000.5532 | simE=000.5393 | regE=000.0138 | optParE=000.0000 | relF=000.0727 | [0m
[34m   46-Img: E=000.5532 | simE=000.5393 | regE=000.0138 |[0m
[31m   47-Tot: E=000.4480 | simE=000.4342 | regE=000.0138 | optParE=000.0000 | relF=000.0727 | [0m
[34m   47-Img: E=000.4480 | simE=000.4342 | regE=000.0138 |[0m
[31m   48-Tot: E=000.3575 | simE=000.3439 | regE=000.0137 | optParE=000.0000 | relF=000.0666 | [0m
[34m   48-Img: E=000.3575 | simE=000.3439 | regE=000.0137 |[0m
[31m   49-Tot: E=000.2884 | simE=000.2748 | regE=000.0136 | optParE=000.0000 | relF=000.0537 | [0m
[34m   49-Img: E=000.2884 | simE=000.2748 | regE=000.0136 |[0m
[32m-->Elapsed time 49.81070[s][0m

Writing out the parameter configuration file ¶

To inspect what settings were in fact used during the run we can write them out using json. We can write out both the actual settings as well as a json file describing what the individual settings mean. If desired one can then modify these settings and load them back in as the new desired setting for a new optimization run.

params.write_JSON( 'test_minimal_registration_' + model_name + '_settings_clean.json')
params.write_JSON_comments( 'test_minimal_registration_' + model_name + '_settings_comments.json')

Out:

Writing parameter file = test_minimal_registration_lddmm_shooting_map_settings_clean.json
Writing parameter file = test_minimal_registration_lddmm_shooting_map_settings_comments.json

Total running time of the script: ( 0 minutes 49.975 seconds)

Gallery generated by Sphinx-Gallery