next up previous contents
Next: 12. Landmarks, Surfaces and Up: 2 B. Anatomical Image Previous: 10. Linear Registration   Contents

Subsections

11. Non Linear Registration

11.1 Introduction

A nonlinear registration ends up being is a displacement field between two images. In particular it takes a voxel in image 1 to a position in image 2. Some non-linear registration methods simply output directly the actual displacement field. Other methods (such as the one used in BioImage Suite) use a parametric representation for the displacement field based on spline models. These later type of methods enable a more compact representation.

BioImage Suite typically uses a tensor cubic-spline representation. These transformations are stored in a custom .grd file format and can be manipulated using both GUI tools (Brain Register) and command line tools - see Section 11.5 for more details on these.

11.2 Visualizing NonLinear Transformations

Figure 11.1: Left: The Visualize Transformations (Jacobian) Tab. Right: Superposition of a grid showing the effect of the transformation on the warped image. This is achieved using controls in the bottom box (C) of the tab shown on the left.
Image jacobian Image nonlinear_grid

Visualizing and quantifying the effect of a nonlinear transformation is a difficult problem. Unlike linear transformations which are global in nature, nonlinear maps can have varying accuracy over the whole image. For example, in the case of a linear registration between two brains, if the map is accurate at a few dispersed points in the brain, we can be confident that it is equally accurate in the rest of the image. Nonlinear registrations on the other hand, can be highly accurate in some parts of the image and less accurate in others.

11.2.0.0.1 Jacobian Maps:

A key property of nonlinear transformations is the determinant of the Jacobian of the transformation - often referred to as ``Jacobian map''.11.1 This is a measure of the local volume expansion required to map the reference image to the target image. If this number is equal to $1$ then there is no volume change, if it is greater than $1$ it implies that the target image is bigger than the reference image, whereas if it is less than $1$ it implies that the reference image had to be shrunk to match the target. If this value goes below zero (or even close to zero), it is a sign of trouble. It implies that the transformation map has gone singular. In such cases, one needs to discard the transformation and repeat the registration with higher smoothness settings.

To quickly compute a Jacobian map for a transformation use the controls in Figure 11.1(left), Panel A. The two entry-boxes allow one to specify the resolution reduction factor (3.0 is the default) from the original image - this is useful for speeding up the process, the Jacobian maps are usually very smooth, so computing them at full resolution is not really needed, and the Intensity threshold factor (0.0 = do this everywhere). Both of these options can be used to speed up the process.

Three different combinations of results can be computed. The Comp Jacobian button simply computes the Jacobian map (as a percentage). The other two buttons Comp Tensor and Comp Strain store the whole tensor and eigenvalues of the Jacobian tensor respectively.

11.2.0.0.2 Showing Grids over the Warped Image:

The Visualize Transform Box, marked as C in Figure 11.1(left), can be used to overlay a grid on the original image showing the effect of the transformation. The spacing of the grid and the intensity of the grid are controlled by the four option menus in the box. They are strictly used for visualization, so trial and error will often reveal the best settings. The new image which consists of a combination of the warped grid and the warped transform image is placed as ``Results'' in the transform viewer. See Figure 11.1(right) for an example.

11.3 Nonrigid Registration (Intensity Based)

Figure 11.2: The Non-Linear Registration Controls. Above: the most common options under the Simple tab. Below: advanced options.
Image nonlinear_simple

Image nonlinear_advanced

Non-Linear Registration allows for additional degrees of freedom and permits the more accurate mapping of images from different subjects into a common space. The methodology in this control derives from the work of Papademetris et al. [49] which in turn is closely related to work by Rueckert et al. [91]. The computed transformation is parameterized in terms of a tensor b-spline grid with uniform control point spacing (e.g. 20.0 mm in the example on the left).

To perform a nonlinear registration, just as in the case of the linear registration above, first load the reference image in the ``Reference'' Viewer and then the target image in the ``Transform Viewer''. Again, similarly to the linear registration case, there are two sets of controls.

In the Simple Controls (Figure 11.2 top), the user specifies the resolution (as a multiple of the native resolution of the reference image), the control point spacing and whether to automatically save and overwrite the resulting transformation. The registration is invoked by pressing the ``Compute Linear + Non Linear Registration'' which first computes an affine mapping followed by a nonrigid mapping. The resulting transformation when saved in the .grd file format includes both the affine and the non-linear components.

Under the Advanced tab (Figure 11.2 bottom) the user can set additional parameters such as the similarity metric (default=Normalized Mutual Information), the optimization method etc. Perhaps the most important of these parameters is the smoothness (the control for this is highlighted in figure 11.2(bottom)), setting this to a higher value ensures a more rigid mapping at the expense of less accuracy in the registration. Here the registration can be invoked in two ways. The first is as a straight non-linear registration using the ``Compute Non Linear Registration'' button - in this case either the current transformation is used for initialization or the linear component of the registration is assumed to be the identity matrix. The second is as combo Linear+Non Linear Registration in the same way as for the simple controls.

11.4 Distortion Correction (Single Axis Distortion)

Figure 11.3: The Distortion Correction Controls. Above: the most common options under the Simple tab. Below: advanced options.
Image distortion_simple

Image distortion_advanced

The Distortion Correction functionality is essentially a constrained form of the non-linear registration in which the warping is assumed to be only in one direction. This reflects the case of distortion in echoplanar magnetic resonance acquisitions - such as those used for DTI and fMRI. The distortion direction is assumed to be along the ``Y-axis'' of the image (which is typical of the phase encode direction), if this is not the case for the specific images, the axis can be changed using the Phase-encode direction drop menu under the advanced controls (Figure 11.3 bottom). See Figure 11.4 for an example.

Most of the functionality of the distortion correction algorithm is closely related to the nonrigid registration methods and is shown in Figure 11.3. There are two key differences: (a) The resolution can be set to be lower i.e. 2.0 or 4.0 since typically the reference image (anatomical) has voxel size of $\approx$ 1mm and the target (echoplanar) has voxel size $\approx$ 4 mm, in which case running the algorithm at 1 or 1.5 mm resolution is overkill. (b) if the Echo-planar image used is a spin-echo image with no signal loss the adjust drop menu can be used to enable the signal-loss conservation constraint (originally from the work of Studholme et al [107]) - this may result in better registrations.

Figure 11.4: Visualization of a distortion correction transformation. This shows displacements only along the phase-encode (y) axis of the echoplanar image.
Image distortion_grid


11.5 Batch Mode Registration

11.5.1 A quick overview of GNU Make

All Batch Job Generators in BioImage Suite result in a makefile which needs to be processed with the Unix Make Utility (most commonly GNU Make). While ``makefiles'' are primarily used to compile large programs they are a great mechanism for general batch mode processing because they enable:

  1. the definition of job dependencies - i.e. job A must be completed before job B
  2. the use of multiple processors at the same time - i.e. run two or more jobs at once.
  3. the ability to automatically recover from a system crash and to restart at the point that the job failed - i.e. no repeat computations.

The standard name for a makefile is unsurprisingly ``makefile'', although in this context I recommend the use of more descriptive names e.g. ``makefile.controls'' etc. Given a makefile, the batch job is executed using

 
make -f makefile.controls
On Linux systems make=gmake i.e. typing gmake is equivalent to typing make - this may appear in some examples. Additional useful flags include
  1. ``-n'' - do a dry run i.e. simply print a list of commands to be executed without doing anything
    make -n -f makefile.controls
    
  2. ``-j'' - specify how many jobs to run at once - typically equal to the number of processors available e.g. to use 2 processors type
    make -j2 -f makefile.controls
    

In addition makefiles contain a number of ``jobs'' which may be explicitly specified. Typically the first job defined in BioImage Suite batch jobs is the equivalent of ``do everything'' so you need not worry about this - in some cases other specifications will need to be given . Note: If after starting a batch job - it is for some reason terminated (either by a system crash or a reboot or ..) it may be restarted by typing exactly the same command as the one used to originally start the batch job. The job dependency mechanism in make will ensure that no processing that was previously completed is re-run.

Microsoft Windows Info: Whereas make is a standard feature on most unix systems (including Mac OS X), to run batch jobs in MS-Windows you will to download and install GNU Make. A binary (from the UnixUtils distribution ) is included with BioImage Suite.

11.5.2 Single Reference Image Registrations

To compute batch-mode registrations to a single reference (typically inter-subject nonlinear brain mappings for the purpose of computing fMRI composite maps etc.) BioImage Suite provides the pxmultiregister_int.tcl tool. This is a batch-mode generator which will generate a makefile - the batch job itself is executed using the standard make program - see detailed description above. pxmultiregister_int.tcl has a large number of options - with more to be added in the next release. Simply typing pxmultiregister_int.tcl sample lists all options, which might be intimidating for the average user. A basic setup file has the form: (Lines beginning with # are treated as comments and may be added liberally to document the processing.

#Example Setup File
# all lines beginning with # are ignored
#
# Mode of Operation
set intensities_only 1
#
# List all images here
set imglist {
 Az2_2_sm0.8.hdr 
 Az2_3_sm0.8.hdr 
 Az2_4_sm0.8.hdr 
}
# Put reference brain here
#
set refimg Az2_5_sm0.8.hdr
# Linear mode -- choices are "rigid", "similarity", "affine" -- 
# this is used to generate an initial transformation to be refined later.
set linearmode "affine"
# Tweak parameters for intensity based part 
# (the ones below are typical for fMRI composite maps)
# Resolution is a scale factor x native 
#(i.e. if reference image is 1.1x1.1x1.2 setting this to 1.5
#   will resample the images to
# 1.65x1.65x1.7 mm). Integer values are not recommended due to 
#   artifacts in joint histogram computation, 1.5 or 1.2 are good values
set resolution 1.5

# Spacing defines the flexibility of the transformation -- 
#   the gap between control points in the tensor b-spline model 
#   used for the non-linear transformation
# 15.0 mm is a good choice for composite functional maps
# For structural morphometric studies this should be reduced to
#  10 or 8 mm (with a corresponding increase in computational time)
# If the data is a fairly low resolution this
#  can be increased to 20mm or so.
set spacing 15.0

# This is the maximum number of iterations for the
# Conjugate Gradient Minimizer (15 is usually a safe number)
set iterations 15

# Leave this to zero unless other instructed
# Regarding filetype: filetype = 1 includes directory name
#    in output filename 0 does not
set filenametype 0

# Set this higher if data quality is low
set smoothness 0.001

# If linearonly=1 then the resulting transformation will be linear
#         i.e. no warping will be attempted
set linearonly 0

Once the setup file is completed, the next step is to decide where the results (and a host of log files) will be stored. Typically this is in a subdirectory (e.g. results). If the setupfile is called controls_setup and we want the results to go to controls_results type:

pxmultiregister_int.tcl controls_setup control_results
This will check for the existence of all files in the setup file, if a file (image) is missing an error will be given. Once this completes OK, the next step is to generate the makefile using

pxmultiregister_int.tcl controls_setup control_results go > makefile.controls

At this point the batch-job is ready and can be started using the make utility described above.

11.5.3 Pairwise Image Registrations

The case of interest here is where one set of images (e.g. thick-slice ``2D'' conventional images) are to be mapped to another set of images (e.g. high resolution 3D anatomical images). To accomplish this use the pxpairwiseregister.tcl script. This is a batch-mode generator which will generate a makefile - the batch job itself is executed using the standard make program - see above for instructions as to how to use GNU make. pxpairwiseregister.tcl has a large number of options. Simply typing pxpairwiseregister.tcl sample lists all options, which might be intimidating for the average user. A basic setup file has the form: (Lines beginning with # are treated as comments and may be added liberally to document the processing.

#Example Setup File
# all lines beginning with # are ignored
#
#
#
# Put reference list here
#
set reflist {
 /data1/brains/1001/refA.hdr
 /data1/brains/1002/refB.hdr
}
#
# Put the target list here
#
set trglist {
 /data1/brains/1001/targA.hdr
 /data1/brains/1002/targB.hdr
}
#
#
# Type of registration
# mode = rigid,affine,nonlinear
set mode "nonlinear"
#
# Tweak parameters
# filetype = 1 includes directory name in output filename 0 does not
# defaults are for rigid/affine/nonlinear
set resolution 1.5
set spacing 15.0
set iterations 15
set stepsize 1
set smoothness auto

See the previous section for more description of the parameters. Executing pxpairwiseregister.tcl is also identical to pxmultiregister_int.tcl, i.e. first

pxpairwiseregister.tcl setup.txt results
then

pxpairwiseregister.tcl setup.txt results go > makefile.results

and then use the make utility to start the makefile.

11.5.4 Reslicing Images and other command line tool using transformations

The BioImage Suite registration tools (both command line and using the graphical user interface), produce transformation files as outputs, typically with a ``.matr'' or a ``.grd'' extensions. (See the File Formats Chapter for more details.) The command line tools described here can be used to apply these transformations to images (reslicing) or to manipulate the transformations in different ways.

11.5.4.0.1 Image Reslicing

Given a reference image, a target image and an appropriate transformation (or set of transformations) the target image can be resliced in the space of the reference image using the command

pxmat_resliceimage.tcl reference_image target_image output_image
        interpolation_mode xform [xform2] [ xform3]

 interpolation_mode = 0,1,3 none,linear,cubic 
         (Avoid linear if the images are in the range 0.1)
 xform = the .matr or .grd file

Note: The direction of the transformations is a common source of confusion. When computing registrations, the estimated transformations is FROM the reference TO the target. This transformation can be used in pxmat_reslice.tcl to move the target image to the reference image. A good rule of thumb is to remember that images move in the opposite direction as transformations.

If multiple transformations are specified (upto 3) then the concatenation of these transformations will be used. For example, consider the case where we have:

The following command will reslice the functional image to the space of the 3D reference brain:

pxmat_resliceimage.tcl ref_image func_image resliced 3 xf1.grd xf2.matr xf3.matr

11.5.4.0.2 Inverting Transformations

Often the inverse transformation is required. This can be accomplished using the pxmat_inverttransform.tcl script. The syntax is:

pxmat_inverttransform.tcl reference_image output_transformation  xform1 
         [ xform2 ] [ xform3 ]

If multiple transformations are specified (upto 3) then the inverse will be the inverse of the concatenation of these transformations.

11.5.4.0.3 Computing Displacement Fields

Sometimes it is desired to get the displacement field for each voxel. This can be accomplished using the pxmat_displacementfield.tcl command. The syntax is:

pxmat_displacementfield.tcl  reference_image output_field  
         xform1 [ xform2 ] [ xform3 ]

The output_field will be a 4D analyze image file with three frames, storing the x, y and z displacement of each voxel in mm (or nominal units). If more than one transformation is specified the final displacement field will be a result of the concatenation of the transformations specified.

11.6 Example: Co-register reference 3D brain with individual 3D brain

Figure 11.5: Methods used to Co-register an individual 3D wholebrain anatomical with a reference 3D brain image.
Image BrainRegisterNL


Image BrainRegisterNL2

The step-by-step instructions, with reference to Figure 11.5, are:

  1. Choose brainregister from the BioImageSuite main menu. Three windows will appear: a Transform Viewer, a Reference Viewer and a BrainRegister menu bar. In the Reference Viewer choose (File | Standard Images | Colin_1mm_stripped), or the filename of your chosen reference brain.
  2. In the Transform Window choose (File | Load)
  3. Choose the filename that refers to the stripped whole brain acquisition from the brain extraction steps explained above and click Open.
  4. On the BrainRegister menu bar, choose (Registration | Non-Linear Registration). A new Registration/OverlayTool window will appear.
  5. In the new Registration/OverlayTool window click ``Compute Linear + Non Linear Registration''. This will compute a Non-Linear Registration between the two loaded images. If the Auto Save Results box is red, BioImageSuite will automatically save the resulting .grd transformation file in the current working directory.
  6. While the registration is running a PXTkConsole window will appear showing the status of the registration. When the registration is complete a ``For Your Information'' window will appear telling you the name of the saved registration file. If the filename already exists the window will ask if you want to overwrite the existing file.
  7. Check the registration by clicking around in the individual 3D brain image (Transform Viewer) and seeing if the crosshairs match a similar anatomical position in the 3D reference image (Reference Viewer).
    1. Brainstem
    2. Anterior Commissure
    3. Cingulate Gyrus

  8. In a unix window, rename the newly generated .grd file to [studynum]_3DtoRef.grd
    • ex: % mv colin_axi_1mm_stripped_tr3567_3d_stripped.grd tr3567_3DtoRef.grd

Figure 11.6: Methods used to load a previously saved non-linear transformation.
Image BrainTransformation2

11.7 Checking 3D to Reference non-linear registrations

The step-by-step instructions, with reference to Figure 11.6, are:

  1. Choose brainregister from the BioImageSuite main menu. Three windows will appear: a Transform Viewer, a Reference Viewer and a BrainRegister menu bar. In the Reference Viewer choose (File | Standard Images | Colin_1mm_stripped), or the filename of your chosen reference brain.
  2. In the Transform Window choose (File | Load)
  3. Choose the filename that refers to the stripped whole brain acquisition from the brain extraction steps explained above and click Open.
  4. On the BrainRegister menu bar, choose (Registration | Transformation). A new Registration/OverlayTool window will appear.
  5. Under the Transformation block Click Load.
  6. Choose the non-linear transformation for the 3D to reference registration (i.e. [studynumber]_3DtoRef.matr) and Click Open.
  7. Under Reslice Options Click Reslice.
  8. Check the registration by clicking around in the individual 3D brain image (Transform Viewer) and seeing if the crosshairs match a similar anatomical position in the 3D reference image (Reference Viewer).

11.8 Remarks

In the next chapter, we discuss the use of surface based registration methods. These often a useful alternative if one is only interested in making specific structures with high accuracy (higher than what is possible with intensity based methods).

In addition to being useful for computing fMRI composite activations, nonlinear registration methods can be used for morphometric studies by examining the properties of the transformation map itself.


next up previous contents
Next: 12. Landmarks, Surfaces and Up: 2 B. Anatomical Image Previous: 10. Linear Registration   Contents