From 4bb15b14453a4ea3c53959522e320e526ce3b185 Mon Sep 17 00:00:00 2001 From: Andreas Schuh Date: Wed, 19 Jul 2017 17:34:24 +0100 Subject: [PATCH 1/3] add: Example symmetric SVFFD registration script, conversion to SVF [Registration] --- ...ster_brain_images_using_symmetric_svffd.py | 94 +++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100755 Modules/Registration/example/register_brain_images_using_symmetric_svffd.py diff --git a/Modules/Registration/example/register_brain_images_using_symmetric_svffd.py b/Modules/Registration/example/register_brain_images_using_symmetric_svffd.py new file mode 100755 index 000000000..fc15ce707 --- /dev/null +++ b/Modules/Registration/example/register_brain_images_using_symmetric_svffd.py @@ -0,0 +1,94 @@ +#!/usr/bin/env python + +""" +Python script illustrating the use of a symmetric pairwise registration of two +brain scans using the SVFFD model 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 SVFFD algorithm 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 + +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 + + +# affine pre-alignment +if not os.path.exists(affdof): + print("\nPerform rigid and affine registration with 9 DOFs") + if os.path.exists(svffd): os.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 not os.path.exists(svffd): + print("\nCompute SVFFD using symmetric energy function") + if os.path.exists(tgtout): os.remove(tgtout) + if os.path.exists(srcout): os.remove(srcout) + if os.path.exists(disp): os.remove(disp) + if os.path.exists(velo): os.remove(velo) + # '-image' argument needed before srcimg because of use of '-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) + +# convert SVFFD to dense diplacement field saved as NIfTI image +if not os.path.exists(disp): + print("\nConvert SVFFD to dense displacement field") + mirtk.convert_dof(svffd, disp, output_format='disp') + +# convert SVFFD to dense velocity field saved as NIfTI image +if not os.path.exists(velo): + print("\nConvert SVFFD to dense velocity field") + mirtk.convert_dof(svffd, velo, output_format='svf') + +# apply SVFFD to transform each image half-way +if not os.path.exists(tgtout): + print("\nTransforming target image using exp(-0.5 v)") + mirtk.transform_image(tgtimg, tgtout, dofin=svffd, Tt=0, St=-0.5) +if not os.path.exists(srcout): + print("\nTransforming source image using exp(+0.5 v)") + mirtk.transform_image(srcimg, srcout, source_invdof=affdof, dofin=svffd, Tt=0, St=+0.5, target=tgtimg) + +# evaluate similarity of transformed images +print("\nEvaluate similarity of transformed images:") +output = mirtk.evaluate_similarity(tgtout, srcout, '-noid', onexit='output') +print(output.replace(", ", "\n")) \ No newline at end of file From d4ec7fcf903d8331dd92bb2beb2e89dcaa6ee73e Mon Sep 17 00:00:00 2001 From: Andreas Schuh Date: Fri, 21 Jul 2017 12:08:49 +0100 Subject: [PATCH 2/3] enh: Move registration examples to top-level Examples folder --- .../Registration/example/register.par => Examples/register.cfg | 0 .../register_brain_images_using_symmetric_svffd.py | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename Modules/Registration/example/register.par => Examples/register.cfg (100%) rename {Modules/Registration/example => Examples}/register_brain_images_using_symmetric_svffd.py (100%) 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/Modules/Registration/example/register_brain_images_using_symmetric_svffd.py b/Examples/register_brain_images_using_symmetric_svffd.py similarity index 100% rename from Modules/Registration/example/register_brain_images_using_symmetric_svffd.py rename to Examples/register_brain_images_using_symmetric_svffd.py From f78df6600b3c3d391bda05be5ab83dc7f147c945 Mon Sep 17 00:00:00 2001 From: Andreas Schuh Date: Fri, 21 Jul 2017 12:54:33 +0100 Subject: [PATCH 3/3] enh: Demonstrate creation of source image deformation sequence --- ...ster_brain_images_using_symmetric_svffd.py | 109 ++++++++++++++---- 1 file changed, 86 insertions(+), 23 deletions(-) diff --git a/Examples/register_brain_images_using_symmetric_svffd.py b/Examples/register_brain_images_using_symmetric_svffd.py index fc15ce707..302016818 100755 --- a/Examples/register_brain_images_using_symmetric_svffd.py +++ b/Examples/register_brain_images_using_symmetric_svffd.py @@ -2,7 +2,7 @@ """ Python script illustrating the use of a symmetric pairwise registration of two -brain scans using the SVFFD model based on the exponential map of a stationary +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 @@ -10,7 +10,7 @@ (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 SVFFD algorithm is done, we transform each image by the half transformation +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. @@ -24,6 +24,34 @@ 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' @@ -44,11 +72,13 @@ 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 not os.path.exists(affdof): - print("\nPerform rigid and affine registration with 9 DOFs") - if os.path.exists(svffd): os.remove(svffd) +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 @@ -57,38 +87,71 @@ mirtk.register(tgtimg, srcimg, '-par', 'Allow shearing', False, model='Rigid+Affine', dofout=affdof, **params) + +# ------------------------------------------------------------------------------ # compute SVFFD using symmetric energy function -if not os.path.exists(svffd): - print("\nCompute SVFFD using symmetric energy function") - if os.path.exists(tgtout): os.remove(tgtout) - if os.path.exists(srcout): os.remove(srcout) - if os.path.exists(disp): os.remove(disp) - if os.path.exists(velo): os.remove(velo) - # '-image' argument needed before srcimg because of use of '-dof_i' +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 not os.path.exists(disp): - print("\nConvert SVFFD to dense displacement field") +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 not os.path.exists(velo): - print("\nConvert SVFFD to dense velocity field") +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 -if not os.path.exists(tgtout): - print("\nTransforming target image using exp(-0.5 v)") +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) -if not os.path.exists(srcout): - print("\nTransforming source image using exp(+0.5 v)") + 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 -print("\nEvaluate similarity of transformed images:") -output = mirtk.evaluate_similarity(tgtout, srcout, '-noid', onexit='output') -print(output.replace(", ", "\n")) \ No newline at end of file +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})