Affine registration: a comparison of several programs

Magnetic Resonance Imaging 22 (2004) 55– 66

Affine registration: a comparison of several programs Peter Zhilkin*, Murray E. Alexander Institute for Biodiagnostics, National Research Council Canada, Manitoba, Canada Received 30 November 2002; accepted 21 May 2003

Abstract Several registration programs with an affine model for the displacement field were tested on various 2D and 3D MRI of the same modality. The following programs were considered: AIR 3.0 (Woods, J. Comp. Assist. Tomogr, 22(1): 139 –152, 1998), COCGV (Ostuni, JMRI, 7(2): 410 – 415, 1997), FLIRT (Jenkinson, Med. Image Analysis, 5(2): 143–156, 2001), Intramodal Registration (Thevenaz, IEEE Trans. Image Proc., 7(1): 27– 41, 1998), SPM (Friston, Human Brain Mapping, 2: 165–189, 1995), and Patch Algorithm (Zhilkin, MRI 18(9): 1143–1150, 2000). Although some of these programs can perform multimodal registration, none was used in such a mode. This paper attempts a fair comparison of the performance of the Patch Algorithm with other programs. However, different settings of the programs’ parameters may further improve the quality of the registration and/or change execution speed. The registered images, the CPU time required to perform the registration, and the error between the registered and reference images, are presented. Most of the programs give comparable accuracies of registration, but their execution times vary considerably. In general the AIR and Patch Algorithm require the least time. The Patch Algorithm can be easily parallelizable on a multi-processor computer. © 2004 Elsevier Inc. All rights reserved. Keywords: Image processing; Affine registration; Unimodal registration; Comparison

1. Introduction A variety of programs exists for correcting motion artifacts between images acquired over time. In this communication we consider registration of images of the same modality and compare the accuracy and computational performance of the Patch Algorithm (PA) [1] with several published registration programs that use an affine model for the motion displacement field. They include Automated Image Registration (AIR), version 3.08, [2]; Correspondence of Closest Gradient Voxels (COCGV), version 5.10, [3]; FMRIB’s Linear Image Registration Tool (FLIRT), version 3.1, [4]; Intramodal Registration (IR), 1997, [5]; Statistical Parametric Mapping (SPM), version SPM99b [6]. Although some of these programs can perform multimodal registration, such features are not considered in this study, and all programs were run in a unimodal setting. The programs are written in C/C⫹⫹ (except SPM, which is written in Matlab and C), and are freely available from the Internet. This significantly simplifies comparison of computational

* Corresponding author. Tel.: ⫹1-204-984-4540; fax: ⫹1(204)-9845472. E-mail address: [email protected] (P. Zhilkin). 0730-725X/04/$ – see front matter © 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.mri.2003.05.004

performance because the programs may be compiled and run on the same hardware. In AIR, a cost function (for which several options are available) is iteratively minimized. The latest version of the program is characterized by “increased speed and accuracy and ability to address a broader range of registration problems” [2]. COCGV registers images by aligning areas of intensity transition. This algorithm can “successfully perform registrations under conditions of unrelated intervolume voxel intensities, significant object displacements, and/or significant amounts of missing data” [3]. FLIRT uses a global optimization method (combining local optimization with several search strategies) to minimize a cost function. It is “more reliable at finding the global minimum than several of the currently available registration packages in common use” [4]. IR uses a robust multiresolution pyramid strategy and a nonlinear optimization algorithm. Most of the iterations are performed at coarse resolution, where the number of pixels is small, so that computational cost is significantly reduced. The estimates are then propagated to finer resolution where fewer iterations are required because of near-optimal initial conditions. This algorithm “has been used successfully for registration of a whole variety of biomedical images” [5]. SPM uses a least squares fitting procedure, based on spatial and intensity transformations

56

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

applied to one image to match the other [6]. PA reduces the registration problem to a linear regression problem based on the linearization of an “image matching” equation [1,7]. This equation is integrated over a set of patches (subvolumes in 3D images) in the image domain, to avoid the instabilities caused by differentiation over a noisy image lattice. The integrations are carried out over a hierarchy of coarse-to-fine resolutions, to reduce the risk of false matches. All programs were compiled and run on an Origin 2000 SGI© computer. The same set of 2D and 3D Magnetic Resonance images (MRI) was presented to each of the tested programs. Most of the programs have various parameters that influence the registration process, and therefore, also the program’s computational performance and the quality of the registered images obtained. For each program we performed a limited number of runs on the same data, using different values for the most important parameters, and present the results that seem best to us. When comparing programs, we tried to use them under approximately the same conditions (if possible), in order that none gets an advantage. For example, we used the same interpolation method (linear) for the final image reconstruction for each program. We also tried to demonstrate the best abilities of the programs by using their ‘key’ features. For example, IR strongly suggests using a cubic B-spline interpolation for image reconstruction. The effect of noise on each of the program’s robustness and quality of the registration is also briefly addressed.

2. Methods The software packages AIR, COCGV, FLIRT, IR, and SPM, their underlying algorithms, advantages and limitations, and recommendations on how to use them, are described by their respective authors in the corresponding references. In this section, the PA, its features, limitations and recommendations on application are briefly described. The PA uses a linearized image matching equation integrated over a given patch P in the image domain [1,7]

冕

P

⌬I hds ⫽ e

冕

I hds ⫹

P

冕

I hu · de ⫺

⭸P

冕

I hⵜ · uds

P

(1) where Ih denotes an image I convolved with a filter h, ⌬Ih is the difference between unregistered target and reference images, u the displacement field, e the mean intensity change between the two images, and ⭸P is the boundary of the patch P. In this note, an affine model for the displacement field is considered that has 6 and 12 unknown param1 The maximal expected displacement umax determines the maximal level of filtering jmax ⬇ log2(兩umax兩), where 兩umax兩 is expressed in pixels.

Table 1 Registration of an affinely distorted 2D MR image of a human brain against the original image Program, cost function, interpolation method AIR, standard deviation of ratio image, linear Patch Algorithm, linear Intramodal Registration, linear

WI (pixels)

CPU (sec)

0.0098

0.7

0.0183 0.0104

0.7 1.9

The warping index (WI), and CPU time required for registration.

eters (plus one extra parameter, an intensity change factor e) for 2D and 3D registration, respectively. The above integral equation is linear with respect to these parameters. Selecting a sufficient number of patches and filters, one can cast registration as a linear regression problem. After the image has been registered at a coarse resolution, the disparity between corresponding image features is reduced, and the image may then be registered at successively finer resolutions. Such a hierarchical approach reduces the risk of false matches. Furthermore, by performing integration over a patch, the algorithm becomes robust to noise. The PA updates registration parameters after each step, and the final registered image is obtained from the image that is being registered (also referred as the target image) by using a single reconstruction (image interpolation). In turn, this allows one to achieve good results using only (fast, low accuracy) linear interpolation. The method has a number of parameters that drive the registration process. The most important parameter that should be known prior registration 1is the maximal expected displacement between two images. Other parameters, such as levels of filtering, number of patches, their size and location, may be heuristically estimated. In the current version of the program, the patches (all of the same size) simply cover the entire image in a uniform way. To avoid “edge” effects, the patches adjacent to image boundaries are not included in the linear system. The reduction of a registration problem to a linear regression problem used in PA obliges users to carefully set certain parameters. Improper choice of maximal filtering level may invalidate the underlying linearized image matching equation, even though linear regression will always yield a solution. Choosing the proper patch size is also important. The rule of thumb is to select the patch size to be large enough that matching features in the images are located in the same patch, but small enough to have the number of patches (which determines the number of equations) sufficiently greater than the number of unknowns. For a dozen unknowns there are many ways to accommodate such requirements. The PA doesn’t minimize the difference between two images explicitly. Instead, it fits the values of certain integrals computed over patches (left-hand side of Eq. 1) to a linear model (right-hand side of Eq. 1). Choosing different patch sizes or different location of the patches will

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

57

Fig. 1. Registration of a 2D MR brain image. Row 1: target image and images registered by AIR and PA; row 2: images registered by IR, FLIRT, and SPM. Rows 3 and 4: the absolute value of the difference between the reference and target/registered images (in the same order).

result in slightly different solutions that are expected to be in the vicinity of a ‘true’ solution. To compare the quality of registration, a single measure

should be chosen. The Mean Square Error (MSE) between an image pair, defined as err ⫽ 兩I'⫺I兩/公(兩I'兩 兩I兩), was selected for this purpose. To avoid ‘edge’ effects, the MSE is

58

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

Table 2 The MSE between the reference and registered images, and CPU time required for registration of a slice of a 3D MR image of a human brain (Fig. 1) Program, cost function, interpolation method

MSE

1. Unregistered image pair 2. AIR, standard deviation of ratio image, linear 3. Patch Algorithm, linear 4. Intramodal Registration, linear 5. FLIRT (‘rigid body’ model), correlation ratio, linear 6. SPM (‘rigid body’ model), linear 7. AIR, least squares, 6-point semi-support sinc 8. Patch Algorithm, spline 9. Intramodal Registration, spline 10. FLIRT (‘rigid body’ model), least squares, linear 11. SPM (‘rigid body’ model), 9-point support sinc

CPU (sec)

0.2614 0.0836

0.7

0.0845 0.0835

0.3 4.3

0.0900

12.6

0.0950

5.9

0.0821

7.9

0.0828 0.0808

0.45 7.4

0.0899

12.4

0.0953

6.2

3.1. Registration of a 2D image with a known affine transformation A 256 ⫻ 256-pixel artificially reconstructed MR image of a human brain was registered against the original one. COCGV was not used because it can register only 3D images. Neither were FLIRT and SPM, which support only the ‘rigid body’ model for displacements in 2D images. The expected (correct) affine matrix (the inverse of the applied affine matrix) and translation2 are 0.97081 冋0.00693

A⫽

usually computed over a central area, avoiding boundaries of the images. In the case of registering image pairs with known ‘true’ transformation, the following geometric warping index (WI) was also computed [5]: WI ⫽

1 N

冘储共A

⫺1

duce the final registered image (unless stated otherwise). All difference images within each figure are shown with the same intensity scale. Same intensity scale implies that the darker the difference image is, the closer is the match between two images.

⫺ A r兲 · x k ⫹ 共 ⫺ A ⫺1 · b ⫺ b r兲储

k

(2) where {A,b} denote the matrix and vector representing the known ‘true’ transformation parameters, and {Ar,br} those restored by registration. Index k ranges over all N pixels/ voxels of the images. The WI is computed in pixel units for isotropic pixels (in 2D example below) or millimeters for anisotropic voxels (in 3D examples below). All programs were compiled (with optimization enabled) and run on an Origin 2000 SGI© computer with 8 CPUs running at 250 MHz clock rate. None of the tested programs took advantage of several CPU availability, though some of the algorithms, e.g., PA, are suitable for parallel processing. The MIPS C/C⫹⫹ compiler was used to compile the programs, Matlab 6.1 was used for SPM.

3. Results The programs were tested on several artificially generated (cases 1,4,9) and real (cases 2,3,5,6,7) MRI. Artificial images were obtained by reconstruction of real images according to a known affine transformation using sinc interpolation [9]. All programs used linear interpolation to pro-

册 冋

册

⫺0.00971 ⫺1.13002 ; b⫽ 1.02034 ⫺3.68154

The affine matrices and translations obtained after registration are very close to the expected ones, e.g., the maximal difference between the expected and obtained elements of the affine matrices are 0.00011, 0.00013, and 0.00011 for AIR, PA, and IR, respectively. The WI (expressed in pixel units), and CPU time required to perform the registration are presented in Table 1. The CPU time includes the processor time required to perform registration and image reconstruction, and does not include input/output. Small values of the WI indicate high quality of registration. 3.2. Registration of a 2D MR brain image A T1-weighted spin echo 3D MRI (256 ⫻ 256 pixels, 35 slices) of a human brain was acquired. The subject then slightly turned his head, introducing an in-plane motion, and another image was acquired. Knowing the purpose of the study the subject attempted to prevent any out-of-plane motion, which is considered small and negligible. A pair of corresponding slices of the 3D images was selected for 2D registration. The target image and images registered by AIR, PA, IR, FLIRT, and SPM (the last two using a ‘rigid body’ model), in that order, are shown in rows 1 and 2 of Fig. 1. The absolute value of the difference between the reference and the target/registered images are shown in rows 3 and 4. The MSE between the reference and registered images was computed over a central part of the images (a 10-pixel wide border was excluded to avoid image border effects). The initial (before registration) MSE is shown in Table 2, line 1. The MSE and CPU time required for the registration are shown in Table 2, lines 2– 6. The same image pair was also registered using different cost functions and/or image interpolation methods. The registered images 2 Translation depends on origin location and is recomputed for the origin located at the image corner.

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

59

Fig. 2. Registration of a pair of 2D contrast-enhanced MRI of a breast. Row 1: target image and images registered by AIR and PA (from left to right); row 2: images registered by IR, FLIRT, and SPM. Rows 3 and 4: the absolute value of the difference between the reference and target/registered images (in the same order).

60

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

Table 3 The MSE between the reference and registered images, and CPU time required for registration of a 2D contrast-enhanced breast image (Fig. 2)

Table 4 Registration of an affinely distorted 3D MR image of a human brain against the original image

Program, cost function, interpolation method

Program, cost function, interpolation method

Unregistered image pair AIR, standard deviation of ratio image, linear Patch Algorithm, linear Intramodal Registration, linear FLIRT (‘rigid body’ model), correlation ratio, linear SPM (‘rigid body’ model), linear

MSE

CPU (sec)

0.3134 0.1239

1.0

0.1210 0.1222 0.1887

0.3 1.4 12.7

0.1924

5.3

were similar to those shown in Fig. 1, and the MSE and CPU time required for registration are presented in Table 2, lines 7–11. All programs registered the image with small, comparable MSEs (0.0808 ⫺ 0.0953 vs. initial MSE ⫽ 0.2614). Visual inspection of the obtained images indicates that registration was successful. The PA has the smallest computation time (0.3 s), while other programs took 0.7 – 12.6 s to perform the registration. The actual computational time of SPM is smaller than stated because the program also prepares and displays various plots during registration. 3.3. Registration of a 2D MR breast image A pair of contrast-enhanced images of a breast (256 ⫻ 256 pixels) was registered. The target image and the images registered by AIR and PA are shown in row 1 of Fig 2; the images registered by IR, FLIRT and SPM are shown in row 2. The absolute value of the difference between the reference and target images before registration and after registration by AIR, PA, IR, FLIRT, and SPM are shown in rows 3 and 4. All programs, except FLIRT and SPM, registered the images properly; the ‘rigid body’ model for displacement field used by those two does not adequately describe the existing displacement field between the images. The MSE between the reference and registered images computed over a central area of the images, and the CPU times required for the registration are presented in Table 3. While the MSEs obtained by AIR, PA, and IR are comparable (0.1239, 0.1210, 0.1222, respectively, vs. 0.3134 for the unregistered image pair), it is evident that in this example the PA is significantly faster (0.3 s.) than AIR and IR (1.0 and 1.4 s.). 3.4. Registration of a 3D image with a known affine transformation Similarly to the 2D case, a sinc-interpolated reconstructed image, transformed according to a prescribed affine transformation, was registered against the original one. The expected (correct) affine matrix (the inverse of the applied affine matrix) and translation are

1. AIR, least squares, linear 2. Patch Algorithm, linear, filtering levels 1(0), 0(0) along X&Y (Z) 3. Patch Algorithm, linear, filtering levels 2(1), 1(0), 0(0) along X&Y (Z) 4. Intramodal Registration, linear 5. SPM, linear 6. COCGV, linear 7. FLIRT, correlation ratio, linear

MSE

WI (mm)

CPU (sec)

0.0534 0.0536

0.0396 0.0570

90 31

0.0534

0.0186

180

0.0533 0.2028 0.1929 0.0552

0.0137 n/a n/a n/a

190 110 1020 550

The MSE, WI (in mm), and CPU time required for the registration.

A⫽

冋 冋 册

0.9774 0.0091 ⫺ 0.0104

⫺ 0.0100 1.0204 ⫺ 0.0071

册

0.0051 ⫺ 0.0154 ; b 1.0102

⫺ 1.0 ⫽ ⫺ 2.5 ⫺ 0.1 The MSE, WI (expressed in mm) and CPU time required to perform the registration, are presented in Table 4. The WI was computed only for the programs with origin location known to the authors. Most of the tested programs performed a proper registration and reconstructed the affine transformation to a high degree of accuracy. For example, the maximal difference between the expected and obtained elements of the affine matrices was 0.0005, 0.0007, 0.0003, 0.0002, 0.0541, 0.0224, and 0.0128 for AIR, PA, PA (second run), IR, SPM, COCGV, and FLIRT, respectively. However, we had difficulties to properly register the image using SPM and COCGV, and the MSEs between the reference and registered images are large (⬃0.2). Comparable low MSEs obtained by the other programs (⬃0.05) and low values of the WI (in the range 0.01 – 0.06 mm) indicate the high quality of the registration (images are not provided). Computational time varies, with PA and AIR being the fastest (0.5 and 1.5 min), and FLIRT and COCGV the slowest (9 and 17 min). In this example, it is clearly demonstrated that the computational performance of the PA depends heavily on the filtering levels used (compare lines 2 and 3). A consequence of this is that it may require more CPU time to register an image pair with a larger displacement between images, than with a smaller one. 3.5. Registration of a 3D MR image of a brain with a small displacement A T1-weighted spin echo image of human brain (256 ⫻ 256 pixels, FOV ⫽ 22 cm, 35 slices, slice thickness approximately 5 mm, no spacing between slices) was acquired in

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

61

Fig. 3. Registration of a 3D spin echo MRI of a human brain with a small displacement. A single slice is shown. Row 1: target image and images registered by AIR, PA, IR (from left to right); row 2: images registered by SPM, COCGV, and FLIRT. Rows 3 and 4: the absolute value of the difference between the target/registered images and the reference image (in the same order).

a 1.5T system (reference image). After the subject slightly moved his head, introducing a small in-plane and out-ofplane motion, another (target) image was acquired. The target image was registered against the reference image. A slice of the 3D image was selected for presentation, and the results obtained (all using linear interpolation) are shown in Fig. 3. The reference image and the images registered by AIR, PA, and IR are shown in row 1 (from left to right). The images registered by SPM, COCGV, and FLIRT are shown in row 2. The corresponding absolute value of the differences between the reference image and the target/registered images are shown in rows 3 and 4. The MSE between the target/registered images and the reference image, computed

over a central part of the images, and CPU time required for registration, are presented in Table 5, lines 1–7. All programs registered the target image properly, with MSE ⬃0.14 (vs. 0.24 for the unregistered pair). It took approximately 4 and 7 s for PA and AIR to do the registration, 40 s – 10 min for other programs. The same images were also registered using different cost functions and/or interpolation methods. The MSE and CPU time are presented in Table 5, lines 8 –12. Most of the programs managed to reduce MSE at the expense of (significantly) increased CPU time, although COCGV and FLIRT registration took approximately the same and moderately increased CPU time, respectively.

62

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

Table 5 The MSE between the reference and registered images, and CPU time required for registration of a 3D brain image with a small displacement (Fig. 3) Program, cost function, interpolation method 1. Unregistered image pair 2. AIR, standard deviation of ratio image, linear 3. Patch Algorithm, linear 4. Intramodal Registration, linear 5. SPM, linear 6. COCGV, linear 7. FLIRT, correlation ratio, linear 8. AIR, least squares, linear 9. Patch Algorithm, spline 10. Intramodal Registration, spline 11. COCGV, spline 12. FLIRT, correlation ratio, sinc, 7-point support

MSE

CPU (sec)

0.2356 0.1414

7

0.1381 0.1401

4 150

0.1408 0.1416 0.1401

39 600 450

0.1401 0.1304 0.1307

23 27 520

0.1374 0.1354

610 560

3.7. Registration of a 3D MR image of a melon with a large displacement An EPI image (256 ⫻ 256, 64 slices, with a slice thickness of approximately 3 mm) of a melon was acquired. Then the melon was moved by hand (thus, introducing both in-plane and out-of-plane motion), and another image acquired. The second image was registered against the first one, and the results are shown in Fig. 5. The reference, target images and the images registered by AIR and PA (from left to right) are shown in row 1. The images registered by IR, SPM, COCGV, and FLIRT are shown in row 2. The absolute values of the difference between the reference and target/registered images are shown in rows 3 and 4. The MSE between the target/registered and reference images and CPU time required for registration are shown in Table 7. All programs successfully registered the target image with a comparable MSEs in the range 0.05 – 0.06. AIR and PA were faster (38 and 35 sec) than other programs, which required from 87 s for SPM up to 66 min for COCGV. 3.8. Effect of noise on quality of registration

3.6. Registration of a 3D MR brain image with a large displacement In another example, a pair of T2*-weighted EPI images of a human brain (256 ⫻ 256, 67 slices, slice thickness approximately 2 mm), this time with large in-plane and out-of-plane displacements, was registered. The results are shown in Fig. 4. The reference and target images, and the images registered by AIR and PA, are presented in row 1. The images registered by IR, SPM, COCGV, and FLIRT are shown in row 2. The absolute value of the difference between the target/registered images and the reference image are shown in the rows 3 and 4. The MSE between the registered and reference images, and CPU time required for registration are shown in Table 6. AIR, IR, and FLIRT did the best, achieving MSE ⬃0.09 (vs. ⬃0.41 for unregistered pair). The PA had difficulty with this image pair. As mentioned above, PA reduces the registration problem to a linear regression problem. For some reason the distribution of the points (the values of the integrals computed over the patches) in the parameter space is such that the solution for linear regression doesn’t correspond to a good match between images. Starting with only relatively coarse filters, with levels 4(3), 1(0), 0(0), 0(0) along X&Y (Z) directions, PA managed to register the image fairly well (with MSE ⬃0.1) at the expense of increased CPU time. SPM and COCGV also had difficulty registering this image pair, and achieved MSE in the range 0.15 – 0.17, which is still significantly smaller than initial MSE ⬃0.41. Overall, AIR had the best computational performance (less than 1 min) in this example, while PA, IR, and FLIRT required 4.3, 6.3, and 10.1 min, respectively. SPM and COCGV registration took slightly less than 2 and approximately 20 min, respectively.

A slice of the 3D MR images of a human brain was used to examine the effect of noise on registration accuracy of the programs (except COCGV that can’t register 2D images). For this purpose, Rician noise with amplitudes of 5, 10, 15, and 20% of the maximal intensity of the images was independently added to the original reference and target images, and the images thus obtained were registered. All programs participating in this test registered the images properly, regardless of noise, as was confirmed by visual inspection of the registered images and the difference between the reference and registered images. The absolute values of the difference between the reference and target/registered images for the 10% noise case are presented in Fig. 6 as an example. 3.9. Registration of a series of 3D images with known affine transformations To get statistically more significant results, a series of 3D images was generated. First, 50 affine matrices A, translation vectors b and residual intensity changes e were generated by small perturbation of the unit matrix, zero translation vector and zero residual intensity change. The perturbations were uniformly distributed over the ranges [⫺0.01, 0.01], [⫺1, 1] pixels, [⫺0.1, 0.1] pixels, and [⫺0.025, 0.025] for each component of the affine matrix, xand y-components of translation, z-component of translation, and intensity change, respectively. Such perturbations yield up to several mm displacements. A T1-weighted spin echo image of human brain (256 ⫻ 256 ⫻ 35, with voxel size 0.86 ⫻ 0.86 ⫻ 5.0 mm, intensity range ⬃0 –25000) was reconstructed according to the generated transformations using sinc interpolation. Then the generated images were

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

63

Fig. 4. Registration of a 3D EPI MRI of a human brain with a large displacement. A single slice is shown. Row 1: reference, target and the images registered by AIR and PA (from left to right); row 2: images registered by IR, SPM, COCGV, and FLIRT. Row 3 and 4: the absolute value of the difference between the target/registered images and the reference image (in the same order).

registered against the original image by all tested programs. AIR was used with SD of ratio image as an objective function (other objective functions were also tried) and default threshold. The image pairs that had significant MSE after registration were registered again with a lower threshold (which yields more voxels for registration), and the best results were kept. Approximately one third of all image pairs required registration with a smaller threshold. IR used the recommended choices for most of the parameters, and all voxels were used in registration. Approximately onequarter of the image pairs were registered using only brighter voxels (selected by the program). PA used filtering levels {1,0,0} with no filtering in the z-direction, patch size

Table 6 MSE and CPU time required for registering 3D image of a human brain with a large displacement (Fig. 4) Program, cost function, interpolation method Unregistered image pair AIR, standard deviation of ratio image, linear Patch Algorithm, linear Intramodal Registration, linear SPM, linear COCGV, linear FLIRT, correlation ratio, linear

MSE

CPU (sec)

0.4120 0.0873

52

0.1044 0.0872 0.1679 0.1523 0.0886

260 380 115 1200 610

64

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

Fig. 5. Registration of a 3D EPI MRI of a melon. A single slice is shown. Row 1: reference, target and images registered by AIR and PA (from left to right); row 2: images registered by IR, SPM, COCGV, and FLIRT. Row 3 and 4: the absolute value of the difference between the reference and target/registered images (in the same order).

was 16 ⫻ 16 ⫻ 5 voxels, and image border patches were not used in registration. Another PA run used filtering levels {2,1,0,0}, keeping other parameters as in the first run. COCGV, SPM, and FLIRT were used with default parameters. The results are provided in Table 8. They include median MSE and WI (expressed in mm), CPU time required to register the entire series, and median CPU time required to register a single image pair. The WI was computed only for the program with known origin location. The results indicate that for the considered image series AIR, PA, IR, and FLIRT on average demonstrated higher registration accuracy (in terms of MSE) than COCGV and SPM. The WI for AIR, PA, and IR all had comparable values. AIR and PA were the most computationally efficient (26 and 15 min

Table 7 The MSE between the reference and registered images, and CPU time required for registering 3D image of a melon with a large displacement (Fig. 5) Program, cost function, interpolation method Unregistered image pair AIR, standard deviation of ratio image, linear Patch Algorithm, linear SPM, linear Intramodal Registration, linear COCGV, linear FLIRT, correlation ratio, linear

MSE

CPU (sec)

0.2861 0.0546

38

0.0530 0.0604 0.0541 0.0580 0.0540

35 87 210 4000 730

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

65

Fig. 6. Registration of a single slice of an MRI corrupted by Rician noise with intensity of 10% of the maximal intensity of the image. Shown are the images of the absolute value of the difference between the reference and target images, and the differences for the registered images by AIR, PA (row 1), IR FLIRT, and SPM (row 2).

for the whole series) followed by SPM and IR (42 and 80 min). FLIRT and COCGV were found to be significantly slower (6.5 and more than 11.5 h). The median CPU time required to register a single image pair was the smallest for AIR (14 s), followed by PA (18 and 37 s), SPM (50 s), IR (100 s), FLIRT (slightly less than 8 min), and COCGV (more than 12 min). Similar results were obtained for the same series of images, but with Rician noise (with intensity Table 8 The median MSE, WI and CPU time required to register a series of 50 3D images of human brain Program

Median Median WI CPU (sec) MSE (mm) total/median

AIR Patch Algorithm, levels {1,0,0} Patch Algorithm, levels {2,1,0,0} Intramodal Registration SPM COCGV FLIRT

0.1032 0.1033 0.1007 0.1021 0.2408 0.2382 0.1336

0.0958 0.0924 0.0537 0.0700 n/a n/a n/a

1570/14.2 910/18.3* 1855/37.1* 4900/102 2520/50 41600**/740** 23400/470

* Mean time, median time is approximately same because CPU time doesn’t depend on image content, given that the same number and size of patches and filtering levels are used. **For 48 registrations, CPU time was not measured for two registrations because of time counter overflow.

up to 5% of the maximal image intensity) added independently to each image, Table 9.

4. Discussion A detailed comparison of the performance of several available registration programs, for both 2D and 3D image registration, using an affine model for the displacement, was Table 9 The median MSE, WI, and CPU time required to register a series of 50 3D images of human brain with Rician noise added independently to each image Program

Median Median WI CPU (sec) MSE (mm) total/median

AIR Patch Algorithm, levels {1,0,0} Patch Algorithm, levels {2,1,0,0} Intramodal Registration SPM COCGV FLIRT

0.1751 0.1738 0.1731 0.1761 0.2426 0.2361 0.1800

0.0970 0.1439 0.0961 0.0858 n/a n/a n/a

981/17.6 949/19* 1853/37* 4100/82 2370/40 43000**/813** 31000/620

* See note in Table 8. **For 49 registrations, CPU time was not measured for one registration because of time counter overflow.

66

Zhilkin, Alexander / Magnetic Resonance Imaging 22 (2004) 55– 66

carried out using a set of MRI of the same modality. Errors and CPU execution times were determined in each case. In three cases, images were artificially transformed and the programs were compared for their accuracy in recovering the (known) affine transformation. Although in many cases the registration accuracy of the programs was comparable, there was considerable variation in their execution times. In general, the execution times for AIR and PA were the shortest among all the programs tested. It is important to note that for registration of image series using AIR and IR, the program parameters were chosen separately for each pair of images, whereas for PA and the other programs the same parameters were used for the entire series. Therefore, the results obtained may be more favorable for AIR and IR in these cases. Some programs had problems registering certain images. In addition, Rician noise of different amplitudes was added independently to the images before registration; all programs demonstrated robustness on these noisy images. The MSE between image pairs was used as a measure to assess the quality of registration. Another measure (WI) was also used for images with known correct transformation. Although the second measure is more appropriate, we would like to note that often the images with unusually large MSE values are registered more poorly and have larger WI values. Therefore, the MSE may be used as a rough measure of the quality of registration, especially since many programs use MSE as their objective function for minimization. It may seem that MSE as a performance measure favors the programs that explicitly minimize intensity differences. Although this may be so, from the results above one may see that programs minimizing objective functions different from MSE have MSEs comparable to (sometimes even smaller than) those obtained by the programs that explicitly minimize MSE. Computational performance of the PA (for a given image size) depends on filtering levels and patch size, and on their number, and does not depend on the image content. The maximal level of filtering is assigned according to the maximal expected displacement. From our practical experience, this maximal level of filtering often is not required, resulting in less computation. However, as shown in one example (case 6), under-filtering may result in a poor match; more filtering (and CPU time) is usually required in such cases. In

PA the affine transformation that registers the target image is a solution of a linear regression system. The accuracy of PA may be increased if robust regression (that eliminates outliers) is used instead. Using multi-processor computers may further increase the computational performance of the PA, because the evaluation of the integrals in Eq. 1 for different patches can be done independently. Acknowledgments We thank Drs. J. Ashburner of Functional Imaging Lab, M. Jenkinson of University of Oxford, J. Ostuni of National Institutes of Health, P. Thevenaz of Swiss Federal Institute of Technology Lausanne, and R. Woods of UCLA School of Medicine for their help with software packages. We also thank Drs. L. Ryner and W. Richter of National Research Council Canada for providing the MRI used in this study. References [1] Zhilkin P, Alexander ME. 3D image registration using a fast noniterative algorithm. Magn Reson Imaging 2000;18(9):1143–50. [2] Woods RP, Grafton ST, Holmes CJ, Cherry SR, Mazziotta JC. Automated image registration: I. General methods and intrasubject, intramodality validation. J Comput Asst Tomogr 1998;22(1):139 –52, http://bishopw.loni.ucla.edu/AIR3/. [3] Ostuni JL, Levin RL, Frank JA, DeCarli C. Correspondence of closest gradient voxels - a robust registration algorithm, J of Magn Reson Imaging, 1997, 7(2):410 – 415, http://www.cc.nih.gov/ldrr/staff/ostuni/ cocgv/cocgv.html. [4] Jenkinson M, Smith S. A global optimisation method for robust affine registration of brain images. Med Image Anal 2001;5(2):143–56, http://www.fmrib.ox.ac.uk/fsl/flirt. [5] Thevenaz P, Ruttimann UE, Unser M. A pyramid approach to subpixel registration based on intensity. IEEE Transac on Image Proc 1998; 7(1):27– 41, http://bigwww.epfl.ch/thevenaz/registration/. [6] Friston KJ, Ashburner J, Frith CD, Poline J-B, Heather JD, Frackowiak RSJ. Spatial registration and normalization of images. Human Brain Mapping 1995;2:165– 89, http://www.fil.ion.ucl.ac.uk/spm/. [7] Alexander ME. Fast hierarchical noniterative registration algorithm. Int J Imaging Sys Technolo 1999;10:242–57. [8] Gupta N, Kanal L. Gradient based image motion estimation without computing gradients. Int J Comput Vision 1997;22(1):81–101. [9] Hajnal JV, Saeed N, Soar EJ, Oatridge A, Young IR, Bydder GM. A registration and interpolation procedure for subvoxel matching of serially acquired MR images. J Comput Asst Tomogr 1995;19(2):289 – 296.