diff --git a/Modules/Registration/example/register.par b/Examples/register.cfg similarity index 100% rename from Modules/Registration/example/register.par rename to Examples/register.cfg diff --git a/Examples/register_brain_images_using_symmetric_svffd.py b/Examples/register_brain_images_using_symmetric_svffd.py new file mode 100755 index 000000000..302016818 --- /dev/null +++ b/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})