add: Example Python scripts for image registration

Examples/register_brain_images_using_symmetric_svffd.py

@@ -0,0 +1,157 @@
+#!/usr/bin/env python
+
+"""
+Python script illustrating the use of a symmetric pairwise registration of two
+brain scans using the SVFFD algorithm based on the exponential map of a stationary
+velocity field (SVF) represented as multi-variate cubic B-spline function.
+
+For illustration, we use N4 bias corrected T1-weighted adult brain MR images of
+the OASIS dataset for which manual annotations are provided by Neuromorphometrics Inc.
+(see MICCAI 2012 Workshop on Multi-Atlas Labeling). These images are not included
+and need to be obtained elsewhere. Alternatively, use other images available to you.
+
+After the registration is done, we transform each image by the half transformation
+to illustrate how to obtain the flow of spatial mappings corresponding to the
+one-parameter subgroup of diffeomorphic mappings generated by a SVF, i.e., exp(tv)
+where exponent t corresponds to the length of the integration interval.
+By default, the length of the integration interval for the full transformation from
+target to source image is given by t=1.
+
+When 'mirtk' Python module is not found, ensure to add the absolute path of the
+lib/python directory of your MIRTK installation to the PYTHONPATH environment variable.
+"""
+
+import os
+import mirtk
+
+
+# ------------------------------------------------------------------------------
+# global settings of mirtk.subprocess module
+mirtk.subprocess.showcmd = True   # whether to print executed commands with arguments
+mirtk.subprocess.threads = 0      # 0: use all available cores
+mirtk.subprocess.onerror = 'exit'
+
+
+# ------------------------------------------------------------------------------
+# auxiliary functions
+def check(path, msg=None):
+    """Check if file does not exist."""
+    if os.path.exists(path):
+        if msg: print("Skip: " + msg)
+        return False
+    else:
+        if msg: print(msg)
+        return True
+
+
+def remove(path):
+    """Delete file if it exists."""
+    if os.path.exists(path):
+        os.remove(path)
+
+
+# ------------------------------------------------------------------------------
+# settings for this example
+tgtuid = '1000'
+srcuid = '1001'
+
+imgdir = os.path.expanduser(os.path.join('~', 'Datasets', 'MAL35', 'N4'))
+imgfmt = os.path.join(imgdir, '{id}.nii.gz')
+tgtimg = imgfmt.format(id=tgtuid)
+srcimg = imgfmt.format(id=srcuid)
+
+affdof = 'aff.dof'
+svffd  = 'svffd.dof.gz'
+disp   = 'disp.nii.gz'
+velo   = 'velo.nii.gz'
+outfmt = '{id}-warped.nii.gz'
+tgtout = outfmt.format(id=tgtuid)
+srcout = outfmt.format(id=srcuid)
+
+params = dict(sim='NMI', bins=64, bg=0, ds=5, be=.001)
+energy = 'SIM[Image dissimilarity](I(1) o T^-0.5, I(2) o T^0.5) + BE[Bending energy](T)'
+mffd   = 'None'  # required for inverse consistent/symmetric energy using SVFFD model
+
+view_sequence = True # open image viewer with source deformation sequence when done
+
+
+# ------------------------------------------------------------------------------
+# affine pre-alignment
+if check(affdof, "Perform rigid and affine registration with 9 DOFs"):
+    remove(svffd)
+    # No shearing allowed to avoid warning with SVFFD registration using -dof_i option
+    # instead of explicitly resampling the input source image in target image space.
+    # The "register" command below may still display a warning because the inverse
+    # homogeneous transformation matrix may be decomposed with small shearing values
+    # due to numerical inaccuracies.
+    mirtk.register(tgtimg, srcimg, '-par', 'Allow shearing', False,
+                   model='Rigid+Affine', dofout=affdof, **params)
+
+
+# ------------------------------------------------------------------------------
+# compute SVFFD using symmetric energy function
+if check(svffd, "Compute SVFFD using symmetric energy function"):
+    remove(tgtout)
+    remove(srcout)
+    remove(disp)
+    remove(velo)
+    # '-image' argument needed before srcimg because we use '-dof_i'
+    mirtk.register(tgtimg, '-image', srcimg, '-dof_i', affdof,
+                   '-par', 'Energy function', energy,
+                   '-par', 'Multi-level transformation', mffd,
+                   model='SVFFD', dofin=affdof, dofout=svffd, **params)
+
+
+# ------------------------------------------------------------------------------
+# illustrate conversion of register command output to other file formats
+
+# convert SVFFD to dense diplacement field saved as NIfTI image
+if check(disp, "Convert SVFFD to dense displacement field"):
+    mirtk.convert_dof(svffd, disp, output_format='disp')
+
+# convert SVFFD to dense velocity field saved as NIfTI image
+if check(velo, "Convert SVFFD to dense velocity field"):
+    mirtk.convert_dof(svffd, velo, output_format='svf')
+
+# apply SVFFD to transform each image half-way
+do_evaluate_similarity = False
+if check(tgtout, "Transforming target image using exp(-0.5 v)"):
+    mirtk.transform_image(tgtimg, tgtout, dofin=svffd, Tt=0, St=-0.5)
+    do_evaluate_similarity = True
+if check(srcout, "Transforming source image using exp(+0.5 v)"):
+    mirtk.transform_image(srcimg, srcout, source_invdof=affdof, dofin=svffd, Tt=0, St=+0.5, target=tgtimg)
+    do_evaluate_similarity = True
+
+
+# ------------------------------------------------------------------------------
+# evaluate similarity of transformed images
+if do_evaluate_similarity:
+    print("Evaluate similarity of transformed images:")
+    output = mirtk.evaluate_similarity(tgtout, srcout, '-noid', onexit='output')
+    print(output.replace(", ", "\n"))
+
+
+# ------------------------------------------------------------------------------
+# gradually deform source image towards target anatomy
+# may be used to create a movie of the deformation
+print("\nGradually deform source image using exp(tv) for t in [0, 1]:")
+steps = 10
+dt = 1. / steps
+frames = []
+for i in range(steps + 1):
+    t = i * dt
+    frame = '{id}-t{t}.nii.gz'.format(id=srcuid, t=t)
+    if check(frame, "Deform source image with t={}".format(t)):
+        dofin = 'Id' if t == 0. else svffd
+        mirtk.transform_image(srcimg, frame, source_invdof=affdof, dofin=dofin, Tt=0, St=t, target=tgtimg)
+    frames.append(frame)
+
+# illustrate creation of image sequence (3D+t) NIfTI from 3D volumes
+srcseq = '{id}-sequence.nii.gz'.format(id=srcuid)
+if check(srcseq, "Make image sequence from individual time steps"):
+    mirtk.combine_images(*frames, output=srcseq)
+
+# open image viewer
+if view_sequence and mirtk.path("view"):
+    print("Opening image viewer with target image and deformed source image sequence...")
+    mirtk.view(target=tgtimg, source=frames, opts={'diff': None, 'xy': None, 'res': 2})