COMPUTE_WARP

compute_warp − Compute a 3D non-linear deformation ("warp") of an MRI volume to match a template.

compute_warp [-h] | [ options-file {keyword=value} ]

compute_warp computes a non-linear deformation for inter-subject normalization of brain scans. The deformation field consists of a set of displacement vectors sampled on a coarser grid. optionsfile contains the parameters for the warp. If ’-’ is specified as optionsfile, no optionsfile is used.

The produced deformation grid is fed to the program warp_reslice which performs the actual reslicing of the volumes.

-h

Display help and a list of recognized keywords with a short description. This output may seem a little cryptical for most people, hence this man-file.

keyword=value

Most major options are controlled using keyword=value assignments in the optionsfile. Options may also be specified using keyword=value pairs on the command line after the optionsfile. Options specified on the command line overrule any definitions in the optionsfile. Note that the keyword=value pair on the command line must be enclosed in quotes if there are spaces in the string, e.g.

compute_warp warp.opt level=info "alpha=0.1 0.05"

to read options from the file warp.opt and to set the verbosity to "info" and to set the alpha parameter. Allowed keyword/value pairs are described the the later section FORMAT OF THE OPTIONSFILE.

This program is only one of a series of programs designed for performing non-linear registration (aka: warp) of a set of Magnetic Resonance (MRI) volumes of the human brain. The program is designed for achieving accurate inter-subject alignment of functional brain images by deforming corresponding co-registered purely anatomical (MR) scans.

The program takes two 8 bit MRI volumes and produces a set of displacement vectors (comprising a displacement field). The displacement field provides 3-dimensional vectors at locations in the template volume. The vector at a given location in the template volume is interpreted as the displacement of the anatomically corresponding location in the subject volume (aka: data volume). The displacement vectors are not computed at every voxel location in the template but on a coarser grid centered on the template volume’s field of view.

A volume’s field of view is the entire collection of space occupied by the voxels in the volume.

Typically (this is entirely controlled in the options file!!!) the displacement vectors are computed on a 4 times coarser grid so that there exists one displacement vector for each 4x4x4 cube of template voxels. This means that if the template volume is 256x256x160 voxels, the displacement field is going to consist of 64x64x40 vectors centered over 4x4x4 segment blocks in the template volume.

When the displacement field is applied (e.g. using warp_reslice or apply_warp ) displacement vectors are needed for locations in between the locations of the displacement vectors. Tri-linear using the 8 surrounding vectors in the grid of displacement vectors resolves this.

This program has been designed and tested for t1 weighted MR images of the human brain. It has been tested on 20ms double spin echo images as well, but t1 is preferable since this modality has the best gray/white matter contrast. The resolution should be better than 1.1mm in at least two dimensions.

The entire procedure for performing non-linear registration of a functional data-set using corresponding anatomical (MR) images consists of the following data steps:

a. Raw MRI’s should be corrected for scanner inhomogeneity.

b. Strip skull & scalp

c. Affine alignment of stripped MRI (b) to Talairach template

d. Non-linear alignment of "affined" subject (c) to template

e. Rigid body (6 parameter) alignment to subjects’ stripped MRI (intra-subject alignment) for each functional scan

f. Reslicing of functional data to template space combining rigid body (e), affine (c) and non-linear (d) transformations.

Ad. a.: Some 3D MR acquisition protocols produces severe inhomogeneity effects, which are clearly visible in the scans as a varying "illumination" of different regions of the scan. This should be corrected for prior to any further processing.

Ad. b.: Brain stripping is quite often done manually. Some reports have been made about accurate semi- or full-automatic procedures. If the MR scanner acquires more than one modality a multi-channel threshold combined with appropriate morphological operations may do the job.

Ad. c.: The affine alignment is computed f.x. by the use of Roger Wood’s AIR package (the align_linear program). The warp_reslice program accepts the AIR format registration files (as well as free style text format). An affine alignment is one that preserves straight lines straight, (but parallel lines do not necessarily remain parallel). An affine transformation is conveniently expressed as a 4x4 matrix. If image coordinates are represented in homogeneous coordinates (4 dimensional representation) the affine transformation from one volume’s coordinate system to the other is represented by the matrix multiplication of the 4x4 transformation matrix with the 4x1 coordinate vector. In order for the transformation to be affine some restrictions apply on the transformation matrix so that the matrix is determined from 12 parameters. Therefore the affine transformation is also labelled the "12 parameter transformation". The 12 parameters take care of rotation, translation and up/down scaling in size in any directions.

Ad. d.: This is the step handled by compute_warp ! The non-linear transformation takes care of what remains to be fixed after the 12-parameter affine alignment has been performed. This means that the subject’s MR volumes need to be resampled according to the initial 12 parameter alignment prior to computing the non-linear warp.

Ad. e.: The rigid-body intra-subject alignment is handled by e.g. Roger Wood’s AIR package. The AIR package may be obtained from

http://bishopw.loni.ucla.edu/AIR3/index.html.

The intra-subject may, similar to the affine transformation, be expressed as a matrix multiplication in homogeneous coordinates only now the matrix is determined from only 6 parameters, 3 for specifying a translation, and 3 for arbitrary rotation.

Ad. f.: The functional scans should always be resampled as few times as necessary since image resolution otherwise may be degraded. The program warp_reslice is designed to perform the entire chain of image transformations in one step combining the 6 parameter functional-to-structural, the 12 parameter and the non-linear transformations in one reslicing operation.

The program uses a hierarchical algorithm for computing the displacement field. The template and data volumes are subsampled to lower resolution and the displacement field is computed on the coarse resolution first. The resolution is increased when the displacement field has converged. The displacement field is made up of fewer vectors in the lower resolution - typically one maintains 1 displacement vector for each 4x4x4 voxels.

For example, for warping a 256x256x160 sized data volume into a template (which has to be same size!) one would typically use three resolution steps: 1/4, 1/2 and full (1/1). In the coarsest resolution step the data and template volumes are subsampled to 64x64x40 voxels and we choose one vector for each 4x4x4 block of voxels (yielding 16x16x10 displacement field vectors).

After 4-8 iterations the displacement field has converged and the program switches to the second resolution step of 128x128x80 voxels. Once again we choose one vector for each 4x4x4 block of voxels now yielding 32x32x20 displacement field vectors. The previous 16x16x10 displacement field is now super-sampled to 32x32x20 using Tri-linear interpolation. And now the program iterates 4-8 times in the second resolution step. Finally in the 3rd resolution step the original 256x256x160 data and template volumes are used and the displacement field is again super-sampled, now to 64x64x40 vectors. The program iterates 4-8 times in the highest resolution and saves the resulting 64x64x40 displacement field.

So what happens in the iterations within a resolution step? Each displacement vector’s block of voxels (4x4x4 in the above example) is shifted locally an integer amount of voxels and the optimal shift is stored as the new displacement vector. The optimal shift is identified by evaluating a cost function at each shift location

C = voxel similarity function + alpha * neighbor-connection

The first term is low when there is good correspondence between voxel values in the template box and the subject box of voxels according to a voxel similarity function. This function may be chosen to be either Least Squares, Modified Least Squares, Correlation, or Mutual Information.

The second term determines how much the displacement vector should be affected by neighboring displacement vectors. This term is low when the displacement vector is similar to the average of the 6 nearest neighboring vectors. This term acts more and more as a smoothness enforcer with increasing values of the parameter alpha.

After each iteration the displacement field is smoothed with a low-pass filter in order to further enforce smoothness. This smoothing is performed independently in the X, Y and Z components of the displacement field vectors using a Gaussian kernel.

As a final regularizer for the warp, the displacement field is reconciled with the inverse mapping. Now, what does this mean ? The program actually computes two warps: One from data to template and one with the role of data and template volumes reversed. Ideally after each iteration, the transformation found from data to template should be the inverse of the transformation from template to data. In practice this is not quite the case. In order to reduce error, these two ideally inverse mappings are corrected using a "reconciling operation" which means that the two displacement fields are adjusted so that they look more like each other’s inverses.

Mutual information is usually a hit. May produce "strange" warps if not regularised enough. It is always a good idea to inspect the warped volumes and compare against the template volumes to ensure that the warped brain still "looks like a normal brain".

Least squares (0) is faster and works well if you are sure to have the same scaling of voxel intensities in all scans. If there are intensity problems present, you may see problems like white matter being mapped to gray (squeezing the gray matter) or similar problems.

Correlation works well if each scan has it’s own global intensity scaling factor.

Local least squares works well if you have problems with local inhomogeneities. It puts more emphasis on matching the edges.

Most of the parameters don’t really need to be set. You may want to chance the resolution steps

There are several ways of evaluating the warp result. One is to compute the warped data volume and compare against the template. The shape of the warped brain should resemble the template and some of the fissures should be matched, some are not. An interactive 3-panel volume browser is recommended for this (e.g. the Matlab program shv.m obtainable from ftp://eivind.imm.dtu.dk/pub/volio/voliomatlab.tar.gz).

Secondly the displacement field may be inspected,- this should look smooth and without fast transitions. The Jacobian volume can be computed from the displacement field. (The Jacobian is the determinant of the derivative matrix of the vector field. If the displacement field is the vector function

F(x,y,z)=(vx(x,y,z), vy(x,y,z), vz(x,y,z))

the Jacobian is

            | dvx/dx, dvy/dx, dvz/dx |
J(x,y,z) =  | dvx/dy, dvy/dy, dvz/dy |
            | dvx/dz, dvy/dz, dvz/dz |

Jacobian values that are either high ( > 4 or so indicating 4 times expansion of the volume) or low ( < 0.25 indicating 4 times compression) are "bad" in the sense that they indicate an unrealistic deformation. The jacobian may be computed using the matlab function below. The program uses the loadvol MEX function obtained from the voliomatlab package (ftp://eivind.imm.dtu.dk/pub/volio/voliomatlab.tar.gz).

function [jac,hjac]=getjac(file)
% [jac,hjac]=getjac(file)
%
% Computes the jacobian from the displacement field in "file"
%
% The cmpix is automatically obtained from the header
% The returned header can be used to save the jacobian volume
htrf=loadhdr(file);
trf=double(loadvol(file));
sz=htrf.size;
x1=1:sz(1)-1; x2=2:sz(1); y1=1:sz(2)-1; y2=2:sz(2);
z1=1:sz(3)-1; z2=2:sz(3);
dtrfx=(trf(x2,:,:,:)-trf(x1,:,:,:))/htrf.cmpix(1);
dtrfy=(trf(:,y2,:,:)-trf(:,y1,:,:))/htrf.cmpix(2);
dtrfz=(trf(:,:,z2,:)-trf(:,:,z1,:))/htrf.cmpix(3);

dd=zeros([htrf.size(1:3)-1,3, 3]);
dd(:,:,:,:,1)=(dtrfx(x1,y1,z1,:)+dtrfx(x1,y2,z1,:)...
   +dtrfx(x1,y1,z2,:)+dtrfx(x1,y2,z2,:))/4;
dd(:,:,:,:,2)=(dtrfy(x1,y1,z1,:)+dtrfy(x2,y1,z1,:)...
   +dtrfy(x1,y1,z2,:)+dtrfy(x2,y1,z2,:))/4;
dd(:,:,:,:,3)=(dtrfz(x1,y1,z1,:)+dtrfz(x2,y1,z1,:)...
   +dtrfz(x1,y2,z1,:)+dtrfz(x2,y2,z1,:))/4;
dd=permute(dd, [4, 5, 1, 2, 3]);
ii = eye(3);
jac=zeros(htrf.size(1:3)-1);
for x=x1,
  for y=y1,
    for z=z1,
      jac(x,y,z)=det(dd(:,:,x,y,z)+ii);
    end,
  end,
  disp(sprintf(’%4.1f’, x/max(x1)*100))
end
hjac = htrf;
hjac.rank = 3;
hjac.size=size(jac);
hjac.cmpix=htrf.cmpix(1:3);
hjac.min=num2str(min(min(min(jac))));
hjac.max=num2str(max(max(max(jac))));
hjac.filename1=’’;

The optionsfile is an ASCII file, each line assigns a value or a list of values to a keyword. Anything after ’#’ on a line is ignored.

Many parameters take a list of values. This is because you may want to use different parameter settings for the different resolution steps. The first resolution step uses the first parameter value specified to all keywords. The program iterates the displacement field in a resolution step until a convergence criteria is met (either maximum number of iterations set via max_repeat or the relative change of the cost function is too small, set via min_rel_change ). Once convergence is reached, the second parameter value of all the keywords is used. You may think of the parameter values in the options file as rows of a matrix, and the program uses a column of the matrix for each resolution step. Note that if a parameter has too few values, the LAST value will be re-used. This means that the keyword with the largest number of values will determine the total number of resolution steps to be applied to the warp.

basedir = string

Base directory to prepend all file names.

templatefile = string list

Name of the template volume. The template volume is the target of the warp. Note that if more than one file is specified the program searches through the list for volumes to use for subsampled versions of the template volume (the original high-res template volume has to be the first in the list).

datafile = string list

Name of the data (aka subject) volume. The data volume is the volume to be warped. Note that if more than one file is specified the program searches through the list for volumes to use for subsampled versions of the data volume (the original high-res template volume has to be the first in the list).

fileformat= string

Determines the default output file format (see volio.h), 0 (default) same as template, 2=VAPET format, 3=ANALYZE format, 4=XPRIME format. The program automatically determines input formats.

outfile_field = string

The name of the generated displacement field data-to-template.

outfile_inv_field = string

The name of the generated reconciled inverse displacement field data-to-template. This is not computed if the reconciled warp feature is disabled (using "no_reconcile=1").

outfile_warped=string

If specified the program will at the end compute the warped data volume and write it to this file. Note that this file is resampled using nearest neighbor sampling which is sub-optimal in most uses. Please use warp_reslice which is capable of Tri-linear and sinc-interpolation.

outfile_inv_warped=string

If specified the program will at the end compute the warped template volume and write it to this file. Note that this file is resampled using nearest neighbor sampling which is sub-optimal in most uses. Please use warp_reslice which is capable of Tri-linear and sinc-interpolation. This is not computed if the reconciled warp feature is disabled.

logfile=string

Technical information about the warp’s inner workings goes into this file.

dumpfile=string

Computational data about the convergence of the warp goes into this text file. The file consists of one row per iteration. The data is

cost_sumsq, cost_sumsqsc, cost_entropy, cost_mutinfo, cost_corr

and 5 more for the inverse mapping if computed.

level=string

Controls amount of output to stdout. Choose from ’silent’, ’std’, ’info’, and ’debug’ for least to most output.

method=integer list

Selects the voxel similarity measure to apply. 0: least square, 1: squared scaled, 2: local mean subtracted least squares, 3: max mutual information, 4: min.entropy (not recommended - mathematically unstable), 5: correlation coefficient.

no_reconcile=integer list

Set this to non-zero to disable the reconciled warp feature. Default is enabled.

autoreg=integer list

Set this to non-zero to enable the automatic on-line setting of the regularizer alpha. Default is disabled.

alpha = float list

Selects the amount of neighborhood connection strength. High values gives more restrained warps. Use values of approximately 0.1 for the Mutual Information, 10-50 for the least squares and correlation coefficient similarity functions.

d_filter_fwhm = float list

Controls amount of brute-force smoothing of vector field in between iterations. The parameter sets the FWHM (measured in vector grid distances) of the Gaussian kernel used for the smoothing.

search = integer list

Set to 1,2,or 3. Controls the size of the local search-region for updating vectors in each iteration. Searches plus minus this many voxels, 1=search 3*3*3=27 places, 2=search 5*5*5=125 places, 3=search 7*7*7=343 locations. 2 is recommended and default.

segx=integer list

segy=integer list

segz=integer list

Sets the number of voxels per displacement vector. Typically 4x4x4.

resx=integer list

resy=integer list

resz=integer list

Sets the resolution to use.

max_repeat = integer list

Sets the maximum number of iterations to run before switching to a higher resolution.

min_rel_change = float

Sets the minimum relative change of the selected voxel similarity measure before switching to a higher resolution.

This is an example options file for warping a 256x256x160 data to a template of identical size producing a displacement field of 64x64x40 vectors using Mutual Information and reconciled warp. Three resolution steps are carried out (64x64x40, 128x128x80 and 256x256x160 for 1/4, 1/2 and full resolution). This example takes approximately 40 minutes on a 266MHz Pentium.

templatefile = b.mri
datafile = a.mri
outfile_field = a2b.trf
outfile_inv_field = b2a.trf
logfile = a2b.log
dumpfile = a2b.dump
level = std
d_filter_fwhm = 0.1
resx = 64 128 256
resy = 64 128 256
resz = 40  80 160
segx = 4
segy = 4
segz = 4
search = 2
method = 3                  # This is Mutual information
alpha =  0.11 0.08 0.06
# above would be something like 40 10 5 if method=0/1/2/5
max_repeat = 8
autoreg = 0
min_rel_change = 0

Failure to make data and template the same size.

Failure to make template of size so that it can be subsampled and divided into an integer number of segments as specified by segx..z.

Failure to select the appropriate voxel similarity measure

Failure to perform the 12 parameter pre-alignment

Failure to correct for scanner inhomogeneities prior to warp

Bad stripping

Scanner/movement artifacts present in scan

Failure to convert the MRI scans to 8 bit prior to warping.

Please tell me if you find any!

Ulrik Kjems <uk@imm.dtu.dk>, see also http://eivind.imm.dtu.dk/staff/kjems

gauss_vol(1), histo_vols(1), subsample(1), file2vapet(1), file2ana(1), file2xprime(1), headerinfo(1), resize_vol(1), convert_vol(1), covar(1), flip(1), region(1), compute_warp(1), warp_reslice(1), warp_reconcile(1),