- Email: [email protected]
3D image registration using a fast noniterative algorithm
Magnetic Resonance Imaging 18 (2000) 1143–1150
3D image registration using a fast noniterative algorithm P. Zhilkin*, M. E. Alexander National Research Council Canada, Institute for Biodiagnostics, 435 Ellice Avenue, Winnipeg, Manitoba R3B 1Y6, Canada Received 20 July 1999; accepted 3 September 2000
Abstract This note describes the implementation of a three-dimensional (3D) registration algorithm, generalizing a previous 2D version [Alexander, Int J Imaging Systems and Technology 1999;10:242–57]. The algorithm solves an integrated form of linearized image matching equation over a set of 3D rectangular sub-volumes (‘patches’) in the image domain. This integrated form avoids numerical instabilities due to differentiation of a noisy image over a lattice, and in addition renders the algorithm robustness to noise. Registration is implemented by first convolving the unregistered images with a set of computationally fast [O(N)] filters, providing four bandpass images for each input image, and integrating the image matching equation over the given patch. Each filter and each patch together provide an independent set of constraints on the displacement field derived by solving a set of linear regression equations. Furthermore, the filters are implemented at a variety of spatial scales, enabling registration parameters at one scale to be used as an input approximation for deriving refined values of those parameters at a finer scale of resolution. This hierarchical procedure is necessary to avoid false matches occurring. Both downsampled and oversampled (undecimating) filtering is implemented. Although the former is computationally fast, it lacks the translation invariance of the latter. Oversampling is required for accurate interpolation that is used in intermediate stages of the algorithm to reconstruct the partially registered from the unregistered image. However, downsampling is useful, and computationally efficient, for preliminary stages of registration when large mismatches are present. The 3D registration algorithm was implemented using a 12-parameter affine model for the displacement: u(x) ⴝ Ax ⴙ b. Linear interpolation was used throughout. Accuracy and timing results for registering various multislice images, obtained by scanning a melon and human volunteers in various stationary positions, is described. The algorithm may be generalized to more general models of the displacement field, and is also well suited to parallel processing. © 2000 Elsevier Science Inc. All rights reserved. Keywords: 3D image processing; Multiscale filtering; Parallel image processing
1. Introduction Motion artifacts are a major problem in interpretation of MR imaging data. Even for short TR, such as in EPI (TR ⬍ 100 ms), inter-image motion can plague the analysis of a time series of such images, that are acquired during functional MR Imaging (fMRI) experiments. For 3D imaging, the time between successive acquisitions is substantially longer than for single-slice (2D) imaging, so motion artifacts are potentially more serious. Several algorithms for realignment of fMRI images are available. Friston [1] considered the problem of realigning and adjusting fMRI time-series. The method of removing the confounding effects of subject movement was proposed. Woods [2] described an automated image registration * Corresponding author. Tel.: ⫹1-204-984-4540; fax: ⫹1-204-9845472. E-mail address: [email protected] (P. Zhilkin).
method based on matching of voxel intensities. In this note, we describe an efficient method of registration that is noniterative, robust to noise, and has a high degree of computational parallelism. It can, in principle, handle arbitrary distortions– global or local– depending on the generality (essentially, the number of parameters) of the displacement model. The method employs a ‘patch’ algorithm, which employs an integrated form of a linearized image matching equation. By using the integrated form, numerical instabilities– caused by differentiation of a noisy function (image) on a discrete lattice– can be avoided, and robustness to noise is greatly improved. The image matching equation has its origins in the “brightness constancy equation” of optical flow [3,4], except that (i) the optical flow field is replaced by a displacement field u that maps features occurring in one of a pair of images to be registered into the corresponding features in the other image; (ii) the spatial mean intensity is allowed to vary between the images. The latter can be important since
0730-725X/00/$ – see front matter © 2000 Elsevier Science Inc. All rights reserved. PII: S 0 7 3 0 - 7 2 5 X ( 0 0 ) 0 0 2 0 9 - 5
1144
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
the acquisition hardware can cause image-to-image fluctuations in the power of the received signal during a long scanning sequence. Prior to registration, the images are convolved with a set of predefined filters. Each filtered image pair is assumed to satisfy, at a given scale of spatial resolution, the same image matching equation (i.e., the equation ‘sees’ the same unknown displacement field u, regardless of the filter used). Each filter will highlight a different set of features in the images, and provides an independent constraint on the displacement field u. In addition, each such filter can be generalized to its multiscale version, allowing features at different scales to be compared. The displacement u will, in general, ‘look’ different at different scales of resolution. However, the field computed at one scale can be used to provide an approximation to the field at another (finer) scale, giving rise to a hierarchical method of registration. In order to simplify the computations, a linearized form of the image matching equation is considered. The hierarchical approach then becomes necessary to avoid the problem of the image matching procedure being trapped in false matches, as explained in the Method section below. The present method is a generalization to 3D of a previous, 2D-patch registration algorithm [4], and is described in the Method section. Accuracy and timing results for testing the algorithm on multislice images of a human brain and of a melon, are given in the Results section. Here, the object was scanned in two or more different stationary positions, and the images were registered against a chosen reference image. Some outstanding problems with the present version of the algorithm, as well as its possible generalizations, are mentioned in the Discussion section.
2. Method Let I(x) and I⬘(x) be two 3D images acquired at different time instances, where x denotes a point in the image. I⬘(x) is the result of changes in I(x) caused by displacement field u(x) and a mean intensity change by factor (1 ⫹ e), i.e., I⬘(x) ⫽ (1 ⫹ e)I(x ⫹ u). The first order approximation to this image matching equation is ⌬I ⬅ I⬘ ⫺ I ⬃ eI ⫹ u 䡠 …I The approximation is valid if u is smaller than the scale of the smallest features appearing in the images. For larger displacements the approximation, generally speaking, is not valid. If so, there is the danger that the equations will converge to a false (local) match instead of the true (global) one [4]. However, by filtering at a sufficiently coarse scale, all fine-scale features smaller than 兩u兩 in size are removed from both images, which in turn reduces the number of features in the image, thus alleviating the problem. If the images I and I⬘ are each convolved with some filter h, a
post-filtered approximation of the above equation can be written ⌬I h ⫽ eI h ⫹ u 䡠 ⵜI h ⫽ eI h ⫹ ⵜ 䡠 共I hu兲 ⫺ I hⵜ 䡠 u,
(1)
where Ih ⫽ I*h, I⬘h ⫽ I⬘*h, and (*) denotes convolution. We are making the tacit assumption that the same displacement field u is present in the original and the filtered imagematching equation. In general, u will appear different at different spatial scales, since different features will be present at each scale. We will assume [4] in Eq. (1) u to be the same for all image pairs Ih, Ih⬘ filtered by a pre-chosen set of filters h at a given scale. Solving Eq. (1) gives an initial estimate of the displacement field u. Using this estimate an approximate distortion correction is applied to the image I⬘, reducing disparity with corresponding features in image I. The residual distortion after such correction is smaller than the original distortion. This allows using a finer-scale filter retaining more details in the images, and guarantees the validity of Eq. (1) for images with finer-scale features present. The corrected image I⬘ is then used as a starting point for a further correction. In this way, image I⬘ is incrementally registered against I. Eq. (1) contains derivatives of the noise-corrupted image Ih whose values are known only at the image lattice points. Numerical differentiation of such a function is ill-posed and leads to unstable solutions [5]. We use a ‘patch’ algorithm [4] to overcome this problem. One can integrate Eq. (1) over a pre-chosen 3D patch V (e.g., a rectangular prism) in the image domain. Using Gauss’s theorem, the divergence term 兰Vⵜ 䡠 (Ihu)dV may be expressed as a flow of the vector Ihu through the surface Sv of the patch V, resulting in the following equation
冕
V
⌬I hdV ⫽ e
冕
I hdV ⫹
V
⫺
冕
冕
I huds
SV
I hⵜ 䡠 udV
(2)
V
In this form, all derivatives have been moved from the image onto displacement field u, which we model as a continuously differentiable function. Solving Eq. (2) with respect to displacement u and mean intensity change e gives the required correction. The optimal placement of the patches is not known, but in general smooth regions provide less reliable estimates than regions with contrast [4]. The choice of filters h is an important part of the method. In particular, in the previous paper [4] it was found that Simoncelli’s matched pair of filters [6] has some advantages over the Daubechies wavelet filters [7]. Therefore, we use a five-tap 1D Simoncelli low-pass (L), high-pass (H) filter pair. The low-pass filter acts also as an interpolating function, and the high-pass filter as a first-derivative of this interpolating function [6]. 2D filters were constructed by
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
1145
Fig. 1. An example of application of Simoncelli downsampling filters to a brain image. A single slice is shown. Image (a) is original; image (b) is the result of applying filters Eq. (3) to image (a); image (c) is obtained by applying the same set of filters to LL image in (b). The images (b) and (c) were modified to improve contrast for the purpose of presentation.
forming the tensor product of these basic 1D filters. If C(x) denotes the 1D low-pass filter, then four 2D filters can be constructed as follows. LL共 x, y兲 ⫽ C共 x兲C共 y兲 LH共 x, y兲 ⫽ C共 x兲
⭸C共 y兲 ⭸y
HL共 x, y兲 ⫽
⭸C共 x兲 C共 y兲 ⭸x
HH共 x, y兲 ⫽
⭸C共 x兲 ⭸C共 y兲 ⭸x ⭸y (3)
These four filters are applied to image slices (i.e., in the x,y-plane), with no filtering in the z-direction, since usually the slice thickness is several times the x,y-pixel size, resulting in the image having a coarser resolution in that direction than in the x,y-plane. We also assume that the displacement in the z-direction is smaller than the thickness of a slice. For larger displacements solving Eq. (2) may result (that doesn’t mean always results) in an incorrect displacement field. It
will be shown in the result section that such simplification of the algorithm is justified for such images since common displacements in fMRI are usually in the range of several mm [1] which is typically of the order of slice thickness. Fig. 1 shows the application of the set of filters (Eq. (3)) to an image slice (a), resulting in four images (b) denoted LL, LH, HL and HH according to the applied filters. Each image is a quarter of the size of the original image (a). Such filtering is called filtering with downsampling, obtained by applying the filter set (Eq. (3)), then decimating the output by a factor 2 in each direction x and y. We refer to this procedure as filtering at level j ⫽ 1. Image details smaller than the width of the filter are washed out in the resulting images. If coarser images are required one can apply the same set of filters to the LL image in (b); the resulting level j ⫽ 2 filtering gives four images each a quarter the size of the LL image (c). This procedure may continue up to some maximum level jmax.
Fig. 2. An example of application of Simoncelli oversampling filters at levels j ⫽ 1,2. Only the LL image is shown. Filtered images become coarser with increasing j. The profile of image intensity along the horizontal white line is shown below the images.
1146
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
and each patch V generates a single equation. The number of patches should be such that the number of equations, nfiltnpatch, exceeds the number of unknowns in u (and e). Using these parameters and a suitable interpolation procedure, the image of interest I⬘ is transformed into image I⬘jmax. This image in turn becomes an input image at the next finer registration level, and so on until the image is registered at the finest level j ⫽ jmin. The value of jmin depends on the signal-to-noise ratio (SNR) and should be chosen smaller (allowing finer details to be registered) for larger SNR. In principle, such a scheme allows handling arbitrary distortions. A linear set of regression equations results if we choose u to have the parameterized form
冘 c ⌽ 共x兲 M
u⫽
k
k
k⫽1
Fig. 3. Generalized registration scheme. Image I⬘ is being registered at levels j ⫽ jmax,. . . jmin. Output at the present registration level is used as input for the subsequent level. I is the reference image, Ih and I⬘h are prefiltered images, u is the displacement, and e the mean intensity change.
In Fig. 2 an example is shown of so-called oversampling filtering, when the resulting images have the same size as the original image. Only the LL images are shown. The scale of the finest detail after j levels of filtering is ⬃2j pixels. The image intensity profile along the horizontal white line is shown below the images in Fig. 2. The higher the level of filtering (j), the smoother is the profile. The noise is also filtered out, and, together with the integration over patches (Eq. (2)), this makes the algorithm robust to noise [4]. The solution of Eq. (2) requires less computation for downsampled images than for the oversampled (fullsize) images. Downsampling is useful in the initial stages if there are large displacements, but does not preserve the translation invariance of the images so that interpolation procedures produce large errors. The redundancy inherent in the oversampled images is necessary to achieve accurate interpolation of the images required between different levels of the registration hierarchy [4] (see below). The generalized registration scheme is shown in Fig. 3. The reference image I and target image I⬘ are first filtered at the coarsest level j ⫽ jmax. The value of jmax depends on the maximum displacement expected (jmax ⬇ log2(兩u兩max), if 兩u兩max is expressed in pixels). A number of patches npatch is selected and placed across the image (overlapping is possible). Eq. (2) is solved for all patches and each of the nfilt filtered image pairs (where nfilt ⫽ number of filters; here, nfilt ⫽ 4), with respect to displacement u and mean intensity change e. Each filter h (of the same spatial scale)
where {⌽1,. . . ,⌽M} is a set of prescribed basis functions over the image domain, and the ck are regression variables (3-component vectors). In general, passing from one level to the next requires image interpolation, since a feature at image grid point x is mapped to the interstitial point x ⫹ u. Every image registered at a given level of resolution is obtained by interpolation of the image using registration parameters computed at a preceding (coarser-scale) level. In this way interpolation errors accumulate, and the quality of the resulting image I⬘reg may deteriorate significantly, depending on the number of registration levels and interpolation method used. In general, a high quality interpolation procedure is required. As a result the registration procedure becomes slower. However, if we limit ourselves to a linear model for displacement–i.e., the affine transformation u ⫽ Ax ⫹ b–we can use the group property of this transformation [4] to incrementally update the parameters A, b and e instead of the images themselves (Fig. 4). In other words, at every registration level j ⫽ jmax,. . . ,jmin, the registered image is obtained from the original image of interest I⬘ using an update of the affine parameters A, b provided by the group transformation rule, and not from the image registered at the previous level. This allows the use of a simple and fast (linear) interpolation procedure to achieve reasonable accuracy and speed up the process of registration. Unfortunately, this method cannot be generalized to handle non-linear distortions. Existing methods [1,2,8] use 6-parameter rigid body model for displacement. The affine model can handle more general motion (shearing) than rigid body model and still is linear. It may be useful in a situation when slices of a 3D image are acquired sequentially (as pointed out by a reviewer), since translation motion between acquisition of the first and last slices causes a shearing effect through the 3D image. In this paper only the affine model for displacement is considered. In this case Eq. (2) becomes linear with respect to the unknown displacement parameters A, b, e. Image
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
Fig. 4. Registration scheme using affine model for the displacement. Image I⬘ is being registered at levels j ⫽ jmax,. . . jmin. I is the reference image, Ih and I⬘h are prefiltered images, A, b are the parameters of the affine transformation, and e is the mean intensity change. Parameters of the affine transformation are updated after registration on each level. The registered image I⬘ at every level is obtained by interpolation of original image I⬘.
registration requires computing integrals in Eq. (2) over a sufficient number of patches and solving the resulting linear regression equations with respect to the 13 parameters A, b, e. Image I⬘ is corrected for distortion using the computed estimates of the parameters and a trilinear interpolation procedure. In the present method, the entire image is simply ‘tiled’ with rectangular patches of uniform size; and images are filtered using 2D oversampling filters, with no filtering in z-direction. At any given registration level, the regression equations are formed from all nfilt ⫽ 4 filtered images, giving rise to 4npatch regression equations. The solution of linear regression equations is computed using Singular Value Decomposition [9].
3. Results T1-weighted spin echo MR brain images of a volunteer were taken using a 1.5T system, with FOV ⫽ 22 cm. After the reference image was taken, the volunteer rotated inplane and tilted his head several degrees (introducing maximal in-plane and out-of-plane motion of approximately 1–2 mm), and another image was taken. Each image consists of 35 slices (slice thickness approximately 5 mm, no spacing between slices), and each slice has 256 ⫻ 256 pixels. The proposed method was applied to these images. In order to avoid boundary-related problems as explained below, only slices 3–32 (count starts from 0) were used for registration. When a transformation is applied to an image, regions containing “no image” usually appear near the boundary of
1147
the image. Contribution of such regions to the MSE between the original and transformed images is small if these regions in both images don’t contain features (brain). In the considered 3D images it may happen (and often does) that such regions contain features in one image and don’t in another. Contribution of such regions to MSE may be large, depending on feature location, displacements, etc. To avoid such artifact the border slices were excluded from consideration in estimating the parameters of the affine transformation, although all slices were used in image transformation. During registration the images were filtered at descending levels of smoothing (increasing levels of resolution): j ⫽ 1,0,0. Linear interpolation was used for all image transformations. A single slice of the images was chosen for presentation, and the result is shown in Fig. 5. Image (a) is the reference image, image (b) is the image of interest, and image (c) is their difference. Image (b) was then registered against (a) as reference. The registered image and the difference between the reference and registered images are shown in (d) and (e), respectively. It is clearly visible that head motion was corrected. Some distortion is still visible in difference image (e). Such residual distortion is to be expected: first, the affine model for displacement cannot compensate for possible nonlinear distortions in this real-life example; and second, linear interpolation itself introduces errors. The mean squared error (MSE), defined as MSE ⫽
储I ⫺ I⬘储
冑储I储储I⬘储
where 储䡠储 denotes the L2 norm of the image function over its domain, was computed (over all slices used in estimating displacement) between images (a) and (b), and (a) and (d). The error decreased almost 50% for the registered image. Table 1 shows how the MSE decreases after each iteration. It is apparent that the image was practically registered after the first step. In another example the T2*-weighted EPI MR brain images (256 ⫻ 256 pixels, 67 slices, slice thickness ⬇3 mm, FOV ⫽ 25 cm) of a volunteer were taken using a 1.5T system. The volunteer was asked to tilt her head several degrees introducing out-of-plane motion approximately 5–10 mm that is 2–3 times greater than the slice thickness. The described method was applied to the pair of images (slices 4 –59, level of smoothness j ⫽ 4,1,1,0 were used for registration). A single slice of the images was chosen for presentation, and the result is shown in Fig. 6. The error (between the presented slices) decreased more than 50% for the registered image. Table 1 shows how the MSE decreases after each iteration. Again, most of the correction is achieved after the first step. A slight increase of the MSE occurred when nonfiltered images were registered (j ⫽ 0). This increase may be explained by the influence of noise, which corrupts the MSE error measure. Although the described method (without filtering in z-direction) in general is not applicable to the images with such large (i.e., several
1148
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
Fig. 5. An example of registering an image using the affine model with linear interpolation (out-of-plane motion smaller than slice thickness). A single slice of the image is shown. (a)–reference image, (b)–image of interest before registration, (c)–their difference, (d)–image obtained after registration of (b) is completed, (d)– difference image between registered and reference images. The MSE (over the slices 3–32 used for registration) before and after registration is also presented.
slice thickness) out-of-plane motion, it works well in this particular example. The displacements obtained for the same pair of images depends on what levels of filtering are used. Sometimes using one set of levels one obtains a more accurate solution (gets certain accuracy faster) than using another set. We tried several sets of levels in the described examples and chose the best one. In another example, EPI MR images of a melon and a human brain were each registered. Images of the melon consist of 64 slices, of 256 ⫻ 256 pixels each (the slice thickness is approximately 2.8 mm); and the brain images contain 20 slices of 128 ⫻ 128 pixels each. The melon was moved by hand (in such way that out-of-plane motion of several mm was introduced) after the reference image was taken. For brain images motion occurred naturally. The melon images were registered at levels j ⫽ 3,2,1,0; the brain at levels j ⫽ 1,0,0; linear interpolation was used in both cases. Similar to previous examples, slices 1–18 for brain, and 4 –59 for melon images were Table 1 MSE between reference and original/registered images, shown in Figures 5 (line 2) and 6 (line 4), as obtained at different levels of registration 1 2 3 4
Registration level MSE Registration level MSE
Original image 0.266 Original image 0.472
j⫽1 0.150 j⫽4 0.314
j⫽0 0.149 j⫽1 .270
j⫽0 0.149 j⫽ 1 0.245
j⫽0 0.247
MSE is computed over all slices used in estimating displacement. Linear interpolation was used for image transformations.
used for displacement estimation stage, and a single slice of the images was selected for presentation. At the top part of Fig. 7 is shown how the difference between the reference and the original/registered melon image evolves after each iteration. At the bottom of Fig. 7 the same is shown for the brain image. Significant residual differences in the ventricles after registration of the brain images may be explained by circulating cerebrospinal fluid, which is more subject to pulsatile circulating changes from image to image than tissue. The MSE was computed for the slices participating in displacement estimation. From the difference images and MSE, one can see a successive improvement after each step. Again, most of the correction is achieved after the first iteration. It was found that starting registration of the melon image at level j ⫽ 3 does not improve the accuracy compared to starting at j ⫽ 2. The results justify the application of the proposed algorithm (without filtering in z-direction) to the images with thick slices and consequently small out-of-plane motion. For the images with thinner slices (or larger out-of-plane motion) the proposed algorithm may arrive at a correct solution (but doesn’t guarantee it, because the solution found may correspond to a local, not global minimum). In general, filtering in z-direction may be required if motion exceeds the voxel dimension in z-direction. The brain images in Fig. 7 were originally obtained for a functional MRI experiment. The purpose of the experiment was to locate activation regions in the brain that participate in performing certain tasks. It is difficult to locate these regions by just looking at the difference image before registration because image misalignment masks the activation
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
1149
Fig. 6. An example of registering a brain image using the affine model with linear interpolation (out-of-plane motion greater than slice thickness). A single slice of the image is shown. (a)–reference image, (b)–image of interest before registration, (c)–their difference, (d)–image obtained after registration of (b) is completed, (d)– difference image between registered and reference images. The MSE (for this slice only) before and after registration is also presented.
zones. However, after the images were registered two activation zones became clearly visible. The CPU time required for registration of images in Fig. 6 was measured on an SGI娀 R10000 computer with 180 MHz clock. Cubic convolution [10] and sinc [8] interpolation procedures were also employed for comparison with linear interpolated results. The MSE after each iteration was also computed, and the results are shown in Table 2. It is evident that most of the correction was achieved after the first step. Also, the cubic convolution and sinc methods produced comparable errors, even though they are computationally more expensive than the linear method.
4. Discussion An algorithm for efficient registration of 3D images has been described. It employs a multiscale filtering approach to registering pairs of images at several levels of spatial resolution, starting from a coarse level of resolution and proceeding to a fine level. This hierarchical approach is necessary to avoid being trapped in a false match caused by displacements being larger than the scale of the features in the image, as discussed in the method section. The proposed method also circumvents the problem of instability caused by numerical differentiation of a noisy image over a lattice
Fig. 7. Registering brain and melon images using a linear displacement model with linear interpolation. A single slice of difference image between reference and original/registered images at different levels is shown. The MSE (for all slices used for registration) is also shown. Images are presented in the same intensity scale.
1150
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
Table 2 Efficiency of the registration method applied to melon and brain images shown in Figure 6 Interpolation method Melon linear Cubic conv. sinc Brain linear Cubic conv. sinc
CPU (sec.)
MSE
300 1400 12000
0.309 3 0.202 3 0.08 3 0.073 3 0.073 0.309 3 0.204 3 0.081 3 0.076 3 0.077 0.309 3 0.208 3 0.083 3 0.077 3 0.077
10 70 600
0.321 3 0.098 3 0.085 3 0.083 0.321 3 0.088 3 0.083 3 0.081 0.321 3 0.099 3 0.086 3 0.083
The melon image of 256 ⫻ 256 ⫻ 64 was registered at levels j ⫽ 3, 2, 1, 0; the brain image of 128 ⫻ 128 ⫻ 20 was registered at levels j ⫽ 1, 0, 0. Slices 4-59 and 1-18 were used in estimating displacement. The MSE, defined over the above slices, was computed after each iteration. Initial MSE between the reference and target images is printed in bold. Computations were done on an SGI Origin 2000 server with R10000 CPU at 180 MHz.
of points, by considering an integrated (‘patch’) form of the image matching equation. This integrated form in turn provides robustness of the registration parameters against noise. In this paper, the method was implemented in its simplest form, namely, by allowing affine displacements of features between images. However, the method can be extended to more general types of distortion, such as global quadratic [4] or higher-order nonlinear global distortions, or local distortions. These more general models involve more parameters than the affine case (e.g., the quadratic model is defined by 31 parameters, including the mean intensity correction e), and thus require more computations. The number of levels of registration is determined by the size of the displacements (which determines an appropriate jmax) and the SNR (which determines an appropriate jmin). If smaller displacements are expected, then the parameter jmax may be chosen smaller, thereby reducing computational time. Most of the computation time is spent in filtering the images (though the computations are still only O(N) in time) and computing integrals in Eq. (2). An outstanding problem is the choice of optimal size and placement of the patches within the image domain. Ideally, these should be determined by the image data, with the size determined by the scale of structural detail of the image features, and the patches located in regions having the greatest contrast. In this way, it may be possible to achieve the same level of accuracy of registration using fewer
patches than the method used in this paper–namely, the uniform ‘titling’ of the image domain, without regard to image structure–and result in saving on computational cost. An important characteristic of the patch algorithm is that it is well suited to a parallel processing implementation. The computation of the integrals in Eq. (2) can be carried out independently for each patch, and these computations involve only the image data local to that patch. Moreover, if the patches have the same volumes, regardless of shape, the number of computations for each patch is the same. Applying oversampling filters Eq. (3) to the images can be also done independently. This property could become important for fast registration when sophisticated nonlinear local distortions have to be modeled—for example, in 3D breast imaging. Acknowledgments We thank Drs W. Richter and L. Ryner, and Ms. C. Sweetland, of the Institute for Biodiagnostics for providing the 3D images. We also thank anonymous reviewers for constructive remarks and suggestions. References [1] Friston KJ, Williams S, Howard R, Frackowiak RSJ, Turner R. Movement-related effects in fMRI time-series. MRM 1996;35:346 – 55. [2] Woods RP, Grafton ST, Holmes CJ, Cherry SR, Mazziotta JC. Automated image registration. I. General methods and intrasubject, intramodality validation. J Comput Assist Tomo 1998;22(1):139 –52. [3] Horn BKP, Schunck BG. Determining optical flow. Artificial Intelligence 1981;17:185–203. [4] Alexander ME. Fast hierarchical noniterative registration algorithm. Int J Imaging Syst Technol 1999;10:242–57. [5] Gupta N, Kanal L. Gradient based image motion estimation without computing gradients. Int J Comput Vision 1997;22:81–101. [6] Simoncelli EP. Design of multidimensional derivative filters. Proc ICIP-94, Austin, TX, 1994. p. 790 – 4. [7] Daubechies I. Ten lectures on wavelets. Philadelphia, PA: Society for Industrial and Applied Mathematics, 1992. [8] 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 Comp Assist Tomogr 1995;19(2): 289 –96. [9] Press WH, Flanney BP, Teukolsky SA, Vetterling WT. Numerical recipes in C. 2nd ed. Cambridge, MA: Cambridge University Press, 1992. [10] Park SK, Schowengerdt RA. Image reconstruction by parametric cubic convolution. Computer Vision, Graphics, and Image Processing 1983;23:258 –72.
3D image registration using a fast noniterative algorithm P. Zhilkin*, M. E. Alexander National Research Council Canada, Institute for Biodiagnostics, 435 Ellice Avenue, Winnipeg, Manitoba R3B 1Y6, Canada Received 20 July 1999; accepted 3 September 2000
Abstract This note describes the implementation of a three-dimensional (3D) registration algorithm, generalizing a previous 2D version [Alexander, Int J Imaging Systems and Technology 1999;10:242–57]. The algorithm solves an integrated form of linearized image matching equation over a set of 3D rectangular sub-volumes (‘patches’) in the image domain. This integrated form avoids numerical instabilities due to differentiation of a noisy image over a lattice, and in addition renders the algorithm robustness to noise. Registration is implemented by first convolving the unregistered images with a set of computationally fast [O(N)] filters, providing four bandpass images for each input image, and integrating the image matching equation over the given patch. Each filter and each patch together provide an independent set of constraints on the displacement field derived by solving a set of linear regression equations. Furthermore, the filters are implemented at a variety of spatial scales, enabling registration parameters at one scale to be used as an input approximation for deriving refined values of those parameters at a finer scale of resolution. This hierarchical procedure is necessary to avoid false matches occurring. Both downsampled and oversampled (undecimating) filtering is implemented. Although the former is computationally fast, it lacks the translation invariance of the latter. Oversampling is required for accurate interpolation that is used in intermediate stages of the algorithm to reconstruct the partially registered from the unregistered image. However, downsampling is useful, and computationally efficient, for preliminary stages of registration when large mismatches are present. The 3D registration algorithm was implemented using a 12-parameter affine model for the displacement: u(x) ⴝ Ax ⴙ b. Linear interpolation was used throughout. Accuracy and timing results for registering various multislice images, obtained by scanning a melon and human volunteers in various stationary positions, is described. The algorithm may be generalized to more general models of the displacement field, and is also well suited to parallel processing. © 2000 Elsevier Science Inc. All rights reserved. Keywords: 3D image processing; Multiscale filtering; Parallel image processing
1. Introduction Motion artifacts are a major problem in interpretation of MR imaging data. Even for short TR, such as in EPI (TR ⬍ 100 ms), inter-image motion can plague the analysis of a time series of such images, that are acquired during functional MR Imaging (fMRI) experiments. For 3D imaging, the time between successive acquisitions is substantially longer than for single-slice (2D) imaging, so motion artifacts are potentially more serious. Several algorithms for realignment of fMRI images are available. Friston [1] considered the problem of realigning and adjusting fMRI time-series. The method of removing the confounding effects of subject movement was proposed. Woods [2] described an automated image registration * Corresponding author. Tel.: ⫹1-204-984-4540; fax: ⫹1-204-9845472. E-mail address: [email protected] (P. Zhilkin).
method based on matching of voxel intensities. In this note, we describe an efficient method of registration that is noniterative, robust to noise, and has a high degree of computational parallelism. It can, in principle, handle arbitrary distortions– global or local– depending on the generality (essentially, the number of parameters) of the displacement model. The method employs a ‘patch’ algorithm, which employs an integrated form of a linearized image matching equation. By using the integrated form, numerical instabilities– caused by differentiation of a noisy function (image) on a discrete lattice– can be avoided, and robustness to noise is greatly improved. The image matching equation has its origins in the “brightness constancy equation” of optical flow [3,4], except that (i) the optical flow field is replaced by a displacement field u that maps features occurring in one of a pair of images to be registered into the corresponding features in the other image; (ii) the spatial mean intensity is allowed to vary between the images. The latter can be important since
0730-725X/00/$ – see front matter © 2000 Elsevier Science Inc. All rights reserved. PII: S 0 7 3 0 - 7 2 5 X ( 0 0 ) 0 0 2 0 9 - 5
1144
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
the acquisition hardware can cause image-to-image fluctuations in the power of the received signal during a long scanning sequence. Prior to registration, the images are convolved with a set of predefined filters. Each filtered image pair is assumed to satisfy, at a given scale of spatial resolution, the same image matching equation (i.e., the equation ‘sees’ the same unknown displacement field u, regardless of the filter used). Each filter will highlight a different set of features in the images, and provides an independent constraint on the displacement field u. In addition, each such filter can be generalized to its multiscale version, allowing features at different scales to be compared. The displacement u will, in general, ‘look’ different at different scales of resolution. However, the field computed at one scale can be used to provide an approximation to the field at another (finer) scale, giving rise to a hierarchical method of registration. In order to simplify the computations, a linearized form of the image matching equation is considered. The hierarchical approach then becomes necessary to avoid the problem of the image matching procedure being trapped in false matches, as explained in the Method section below. The present method is a generalization to 3D of a previous, 2D-patch registration algorithm [4], and is described in the Method section. Accuracy and timing results for testing the algorithm on multislice images of a human brain and of a melon, are given in the Results section. Here, the object was scanned in two or more different stationary positions, and the images were registered against a chosen reference image. Some outstanding problems with the present version of the algorithm, as well as its possible generalizations, are mentioned in the Discussion section.
2. Method Let I(x) and I⬘(x) be two 3D images acquired at different time instances, where x denotes a point in the image. I⬘(x) is the result of changes in I(x) caused by displacement field u(x) and a mean intensity change by factor (1 ⫹ e), i.e., I⬘(x) ⫽ (1 ⫹ e)I(x ⫹ u). The first order approximation to this image matching equation is ⌬I ⬅ I⬘ ⫺ I ⬃ eI ⫹ u 䡠 …I The approximation is valid if u is smaller than the scale of the smallest features appearing in the images. For larger displacements the approximation, generally speaking, is not valid. If so, there is the danger that the equations will converge to a false (local) match instead of the true (global) one [4]. However, by filtering at a sufficiently coarse scale, all fine-scale features smaller than 兩u兩 in size are removed from both images, which in turn reduces the number of features in the image, thus alleviating the problem. If the images I and I⬘ are each convolved with some filter h, a
post-filtered approximation of the above equation can be written ⌬I h ⫽ eI h ⫹ u 䡠 ⵜI h ⫽ eI h ⫹ ⵜ 䡠 共I hu兲 ⫺ I hⵜ 䡠 u,
(1)
where Ih ⫽ I*h, I⬘h ⫽ I⬘*h, and (*) denotes convolution. We are making the tacit assumption that the same displacement field u is present in the original and the filtered imagematching equation. In general, u will appear different at different spatial scales, since different features will be present at each scale. We will assume [4] in Eq. (1) u to be the same for all image pairs Ih, Ih⬘ filtered by a pre-chosen set of filters h at a given scale. Solving Eq. (1) gives an initial estimate of the displacement field u. Using this estimate an approximate distortion correction is applied to the image I⬘, reducing disparity with corresponding features in image I. The residual distortion after such correction is smaller than the original distortion. This allows using a finer-scale filter retaining more details in the images, and guarantees the validity of Eq. (1) for images with finer-scale features present. The corrected image I⬘ is then used as a starting point for a further correction. In this way, image I⬘ is incrementally registered against I. Eq. (1) contains derivatives of the noise-corrupted image Ih whose values are known only at the image lattice points. Numerical differentiation of such a function is ill-posed and leads to unstable solutions [5]. We use a ‘patch’ algorithm [4] to overcome this problem. One can integrate Eq. (1) over a pre-chosen 3D patch V (e.g., a rectangular prism) in the image domain. Using Gauss’s theorem, the divergence term 兰Vⵜ 䡠 (Ihu)dV may be expressed as a flow of the vector Ihu through the surface Sv of the patch V, resulting in the following equation
冕
V
⌬I hdV ⫽ e
冕
I hdV ⫹
V
⫺
冕
冕
I huds
SV
I hⵜ 䡠 udV
(2)
V
In this form, all derivatives have been moved from the image onto displacement field u, which we model as a continuously differentiable function. Solving Eq. (2) with respect to displacement u and mean intensity change e gives the required correction. The optimal placement of the patches is not known, but in general smooth regions provide less reliable estimates than regions with contrast [4]. The choice of filters h is an important part of the method. In particular, in the previous paper [4] it was found that Simoncelli’s matched pair of filters [6] has some advantages over the Daubechies wavelet filters [7]. Therefore, we use a five-tap 1D Simoncelli low-pass (L), high-pass (H) filter pair. The low-pass filter acts also as an interpolating function, and the high-pass filter as a first-derivative of this interpolating function [6]. 2D filters were constructed by
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
1145
Fig. 1. An example of application of Simoncelli downsampling filters to a brain image. A single slice is shown. Image (a) is original; image (b) is the result of applying filters Eq. (3) to image (a); image (c) is obtained by applying the same set of filters to LL image in (b). The images (b) and (c) were modified to improve contrast for the purpose of presentation.
forming the tensor product of these basic 1D filters. If C(x) denotes the 1D low-pass filter, then four 2D filters can be constructed as follows. LL共 x, y兲 ⫽ C共 x兲C共 y兲 LH共 x, y兲 ⫽ C共 x兲
⭸C共 y兲 ⭸y
HL共 x, y兲 ⫽
⭸C共 x兲 C共 y兲 ⭸x
HH共 x, y兲 ⫽
⭸C共 x兲 ⭸C共 y兲 ⭸x ⭸y (3)
These four filters are applied to image slices (i.e., in the x,y-plane), with no filtering in the z-direction, since usually the slice thickness is several times the x,y-pixel size, resulting in the image having a coarser resolution in that direction than in the x,y-plane. We also assume that the displacement in the z-direction is smaller than the thickness of a slice. For larger displacements solving Eq. (2) may result (that doesn’t mean always results) in an incorrect displacement field. It
will be shown in the result section that such simplification of the algorithm is justified for such images since common displacements in fMRI are usually in the range of several mm [1] which is typically of the order of slice thickness. Fig. 1 shows the application of the set of filters (Eq. (3)) to an image slice (a), resulting in four images (b) denoted LL, LH, HL and HH according to the applied filters. Each image is a quarter of the size of the original image (a). Such filtering is called filtering with downsampling, obtained by applying the filter set (Eq. (3)), then decimating the output by a factor 2 in each direction x and y. We refer to this procedure as filtering at level j ⫽ 1. Image details smaller than the width of the filter are washed out in the resulting images. If coarser images are required one can apply the same set of filters to the LL image in (b); the resulting level j ⫽ 2 filtering gives four images each a quarter the size of the LL image (c). This procedure may continue up to some maximum level jmax.
Fig. 2. An example of application of Simoncelli oversampling filters at levels j ⫽ 1,2. Only the LL image is shown. Filtered images become coarser with increasing j. The profile of image intensity along the horizontal white line is shown below the images.
1146
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
and each patch V generates a single equation. The number of patches should be such that the number of equations, nfiltnpatch, exceeds the number of unknowns in u (and e). Using these parameters and a suitable interpolation procedure, the image of interest I⬘ is transformed into image I⬘jmax. This image in turn becomes an input image at the next finer registration level, and so on until the image is registered at the finest level j ⫽ jmin. The value of jmin depends on the signal-to-noise ratio (SNR) and should be chosen smaller (allowing finer details to be registered) for larger SNR. In principle, such a scheme allows handling arbitrary distortions. A linear set of regression equations results if we choose u to have the parameterized form
冘 c ⌽ 共x兲 M
u⫽
k
k
k⫽1
Fig. 3. Generalized registration scheme. Image I⬘ is being registered at levels j ⫽ jmax,. . . jmin. Output at the present registration level is used as input for the subsequent level. I is the reference image, Ih and I⬘h are prefiltered images, u is the displacement, and e the mean intensity change.
In Fig. 2 an example is shown of so-called oversampling filtering, when the resulting images have the same size as the original image. Only the LL images are shown. The scale of the finest detail after j levels of filtering is ⬃2j pixels. The image intensity profile along the horizontal white line is shown below the images in Fig. 2. The higher the level of filtering (j), the smoother is the profile. The noise is also filtered out, and, together with the integration over patches (Eq. (2)), this makes the algorithm robust to noise [4]. The solution of Eq. (2) requires less computation for downsampled images than for the oversampled (fullsize) images. Downsampling is useful in the initial stages if there are large displacements, but does not preserve the translation invariance of the images so that interpolation procedures produce large errors. The redundancy inherent in the oversampled images is necessary to achieve accurate interpolation of the images required between different levels of the registration hierarchy [4] (see below). The generalized registration scheme is shown in Fig. 3. The reference image I and target image I⬘ are first filtered at the coarsest level j ⫽ jmax. The value of jmax depends on the maximum displacement expected (jmax ⬇ log2(兩u兩max), if 兩u兩max is expressed in pixels). A number of patches npatch is selected and placed across the image (overlapping is possible). Eq. (2) is solved for all patches and each of the nfilt filtered image pairs (where nfilt ⫽ number of filters; here, nfilt ⫽ 4), with respect to displacement u and mean intensity change e. Each filter h (of the same spatial scale)
where {⌽1,. . . ,⌽M} is a set of prescribed basis functions over the image domain, and the ck are regression variables (3-component vectors). In general, passing from one level to the next requires image interpolation, since a feature at image grid point x is mapped to the interstitial point x ⫹ u. Every image registered at a given level of resolution is obtained by interpolation of the image using registration parameters computed at a preceding (coarser-scale) level. In this way interpolation errors accumulate, and the quality of the resulting image I⬘reg may deteriorate significantly, depending on the number of registration levels and interpolation method used. In general, a high quality interpolation procedure is required. As a result the registration procedure becomes slower. However, if we limit ourselves to a linear model for displacement–i.e., the affine transformation u ⫽ Ax ⫹ b–we can use the group property of this transformation [4] to incrementally update the parameters A, b and e instead of the images themselves (Fig. 4). In other words, at every registration level j ⫽ jmax,. . . ,jmin, the registered image is obtained from the original image of interest I⬘ using an update of the affine parameters A, b provided by the group transformation rule, and not from the image registered at the previous level. This allows the use of a simple and fast (linear) interpolation procedure to achieve reasonable accuracy and speed up the process of registration. Unfortunately, this method cannot be generalized to handle non-linear distortions. Existing methods [1,2,8] use 6-parameter rigid body model for displacement. The affine model can handle more general motion (shearing) than rigid body model and still is linear. It may be useful in a situation when slices of a 3D image are acquired sequentially (as pointed out by a reviewer), since translation motion between acquisition of the first and last slices causes a shearing effect through the 3D image. In this paper only the affine model for displacement is considered. In this case Eq. (2) becomes linear with respect to the unknown displacement parameters A, b, e. Image
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
Fig. 4. Registration scheme using affine model for the displacement. Image I⬘ is being registered at levels j ⫽ jmax,. . . jmin. I is the reference image, Ih and I⬘h are prefiltered images, A, b are the parameters of the affine transformation, and e is the mean intensity change. Parameters of the affine transformation are updated after registration on each level. The registered image I⬘ at every level is obtained by interpolation of original image I⬘.
registration requires computing integrals in Eq. (2) over a sufficient number of patches and solving the resulting linear regression equations with respect to the 13 parameters A, b, e. Image I⬘ is corrected for distortion using the computed estimates of the parameters and a trilinear interpolation procedure. In the present method, the entire image is simply ‘tiled’ with rectangular patches of uniform size; and images are filtered using 2D oversampling filters, with no filtering in z-direction. At any given registration level, the regression equations are formed from all nfilt ⫽ 4 filtered images, giving rise to 4npatch regression equations. The solution of linear regression equations is computed using Singular Value Decomposition [9].
3. Results T1-weighted spin echo MR brain images of a volunteer were taken using a 1.5T system, with FOV ⫽ 22 cm. After the reference image was taken, the volunteer rotated inplane and tilted his head several degrees (introducing maximal in-plane and out-of-plane motion of approximately 1–2 mm), and another image was taken. Each image consists of 35 slices (slice thickness approximately 5 mm, no spacing between slices), and each slice has 256 ⫻ 256 pixels. The proposed method was applied to these images. In order to avoid boundary-related problems as explained below, only slices 3–32 (count starts from 0) were used for registration. When a transformation is applied to an image, regions containing “no image” usually appear near the boundary of
1147
the image. Contribution of such regions to the MSE between the original and transformed images is small if these regions in both images don’t contain features (brain). In the considered 3D images it may happen (and often does) that such regions contain features in one image and don’t in another. Contribution of such regions to MSE may be large, depending on feature location, displacements, etc. To avoid such artifact the border slices were excluded from consideration in estimating the parameters of the affine transformation, although all slices were used in image transformation. During registration the images were filtered at descending levels of smoothing (increasing levels of resolution): j ⫽ 1,0,0. Linear interpolation was used for all image transformations. A single slice of the images was chosen for presentation, and the result is shown in Fig. 5. Image (a) is the reference image, image (b) is the image of interest, and image (c) is their difference. Image (b) was then registered against (a) as reference. The registered image and the difference between the reference and registered images are shown in (d) and (e), respectively. It is clearly visible that head motion was corrected. Some distortion is still visible in difference image (e). Such residual distortion is to be expected: first, the affine model for displacement cannot compensate for possible nonlinear distortions in this real-life example; and second, linear interpolation itself introduces errors. The mean squared error (MSE), defined as MSE ⫽
储I ⫺ I⬘储
冑储I储储I⬘储
where 储䡠储 denotes the L2 norm of the image function over its domain, was computed (over all slices used in estimating displacement) between images (a) and (b), and (a) and (d). The error decreased almost 50% for the registered image. Table 1 shows how the MSE decreases after each iteration. It is apparent that the image was practically registered after the first step. In another example the T2*-weighted EPI MR brain images (256 ⫻ 256 pixels, 67 slices, slice thickness ⬇3 mm, FOV ⫽ 25 cm) of a volunteer were taken using a 1.5T system. The volunteer was asked to tilt her head several degrees introducing out-of-plane motion approximately 5–10 mm that is 2–3 times greater than the slice thickness. The described method was applied to the pair of images (slices 4 –59, level of smoothness j ⫽ 4,1,1,0 were used for registration). A single slice of the images was chosen for presentation, and the result is shown in Fig. 6. The error (between the presented slices) decreased more than 50% for the registered image. Table 1 shows how the MSE decreases after each iteration. Again, most of the correction is achieved after the first step. A slight increase of the MSE occurred when nonfiltered images were registered (j ⫽ 0). This increase may be explained by the influence of noise, which corrupts the MSE error measure. Although the described method (without filtering in z-direction) in general is not applicable to the images with such large (i.e., several
1148
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
Fig. 5. An example of registering an image using the affine model with linear interpolation (out-of-plane motion smaller than slice thickness). A single slice of the image is shown. (a)–reference image, (b)–image of interest before registration, (c)–their difference, (d)–image obtained after registration of (b) is completed, (d)– difference image between registered and reference images. The MSE (over the slices 3–32 used for registration) before and after registration is also presented.
slice thickness) out-of-plane motion, it works well in this particular example. The displacements obtained for the same pair of images depends on what levels of filtering are used. Sometimes using one set of levels one obtains a more accurate solution (gets certain accuracy faster) than using another set. We tried several sets of levels in the described examples and chose the best one. In another example, EPI MR images of a melon and a human brain were each registered. Images of the melon consist of 64 slices, of 256 ⫻ 256 pixels each (the slice thickness is approximately 2.8 mm); and the brain images contain 20 slices of 128 ⫻ 128 pixels each. The melon was moved by hand (in such way that out-of-plane motion of several mm was introduced) after the reference image was taken. For brain images motion occurred naturally. The melon images were registered at levels j ⫽ 3,2,1,0; the brain at levels j ⫽ 1,0,0; linear interpolation was used in both cases. Similar to previous examples, slices 1–18 for brain, and 4 –59 for melon images were Table 1 MSE between reference and original/registered images, shown in Figures 5 (line 2) and 6 (line 4), as obtained at different levels of registration 1 2 3 4
Registration level MSE Registration level MSE
Original image 0.266 Original image 0.472
j⫽1 0.150 j⫽4 0.314
j⫽0 0.149 j⫽1 .270
j⫽0 0.149 j⫽ 1 0.245
j⫽0 0.247
MSE is computed over all slices used in estimating displacement. Linear interpolation was used for image transformations.
used for displacement estimation stage, and a single slice of the images was selected for presentation. At the top part of Fig. 7 is shown how the difference between the reference and the original/registered melon image evolves after each iteration. At the bottom of Fig. 7 the same is shown for the brain image. Significant residual differences in the ventricles after registration of the brain images may be explained by circulating cerebrospinal fluid, which is more subject to pulsatile circulating changes from image to image than tissue. The MSE was computed for the slices participating in displacement estimation. From the difference images and MSE, one can see a successive improvement after each step. Again, most of the correction is achieved after the first iteration. It was found that starting registration of the melon image at level j ⫽ 3 does not improve the accuracy compared to starting at j ⫽ 2. The results justify the application of the proposed algorithm (without filtering in z-direction) to the images with thick slices and consequently small out-of-plane motion. For the images with thinner slices (or larger out-of-plane motion) the proposed algorithm may arrive at a correct solution (but doesn’t guarantee it, because the solution found may correspond to a local, not global minimum). In general, filtering in z-direction may be required if motion exceeds the voxel dimension in z-direction. The brain images in Fig. 7 were originally obtained for a functional MRI experiment. The purpose of the experiment was to locate activation regions in the brain that participate in performing certain tasks. It is difficult to locate these regions by just looking at the difference image before registration because image misalignment masks the activation
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
1149
Fig. 6. An example of registering a brain image using the affine model with linear interpolation (out-of-plane motion greater than slice thickness). A single slice of the image is shown. (a)–reference image, (b)–image of interest before registration, (c)–their difference, (d)–image obtained after registration of (b) is completed, (d)– difference image between registered and reference images. The MSE (for this slice only) before and after registration is also presented.
zones. However, after the images were registered two activation zones became clearly visible. The CPU time required for registration of images in Fig. 6 was measured on an SGI娀 R10000 computer with 180 MHz clock. Cubic convolution [10] and sinc [8] interpolation procedures were also employed for comparison with linear interpolated results. The MSE after each iteration was also computed, and the results are shown in Table 2. It is evident that most of the correction was achieved after the first step. Also, the cubic convolution and sinc methods produced comparable errors, even though they are computationally more expensive than the linear method.
4. Discussion An algorithm for efficient registration of 3D images has been described. It employs a multiscale filtering approach to registering pairs of images at several levels of spatial resolution, starting from a coarse level of resolution and proceeding to a fine level. This hierarchical approach is necessary to avoid being trapped in a false match caused by displacements being larger than the scale of the features in the image, as discussed in the method section. The proposed method also circumvents the problem of instability caused by numerical differentiation of a noisy image over a lattice
Fig. 7. Registering brain and melon images using a linear displacement model with linear interpolation. A single slice of difference image between reference and original/registered images at different levels is shown. The MSE (for all slices used for registration) is also shown. Images are presented in the same intensity scale.
1150
P. Zhilkin, M.E. Alexander / Magnetic Resonance Imaging 18 (2000) 1143–1150
Table 2 Efficiency of the registration method applied to melon and brain images shown in Figure 6 Interpolation method Melon linear Cubic conv. sinc Brain linear Cubic conv. sinc
CPU (sec.)
MSE
300 1400 12000
0.309 3 0.202 3 0.08 3 0.073 3 0.073 0.309 3 0.204 3 0.081 3 0.076 3 0.077 0.309 3 0.208 3 0.083 3 0.077 3 0.077
10 70 600
0.321 3 0.098 3 0.085 3 0.083 0.321 3 0.088 3 0.083 3 0.081 0.321 3 0.099 3 0.086 3 0.083
The melon image of 256 ⫻ 256 ⫻ 64 was registered at levels j ⫽ 3, 2, 1, 0; the brain image of 128 ⫻ 128 ⫻ 20 was registered at levels j ⫽ 1, 0, 0. Slices 4-59 and 1-18 were used in estimating displacement. The MSE, defined over the above slices, was computed after each iteration. Initial MSE between the reference and target images is printed in bold. Computations were done on an SGI Origin 2000 server with R10000 CPU at 180 MHz.
of points, by considering an integrated (‘patch’) form of the image matching equation. This integrated form in turn provides robustness of the registration parameters against noise. In this paper, the method was implemented in its simplest form, namely, by allowing affine displacements of features between images. However, the method can be extended to more general types of distortion, such as global quadratic [4] or higher-order nonlinear global distortions, or local distortions. These more general models involve more parameters than the affine case (e.g., the quadratic model is defined by 31 parameters, including the mean intensity correction e), and thus require more computations. The number of levels of registration is determined by the size of the displacements (which determines an appropriate jmax) and the SNR (which determines an appropriate jmin). If smaller displacements are expected, then the parameter jmax may be chosen smaller, thereby reducing computational time. Most of the computation time is spent in filtering the images (though the computations are still only O(N) in time) and computing integrals in Eq. (2). An outstanding problem is the choice of optimal size and placement of the patches within the image domain. Ideally, these should be determined by the image data, with the size determined by the scale of structural detail of the image features, and the patches located in regions having the greatest contrast. In this way, it may be possible to achieve the same level of accuracy of registration using fewer
patches than the method used in this paper–namely, the uniform ‘titling’ of the image domain, without regard to image structure–and result in saving on computational cost. An important characteristic of the patch algorithm is that it is well suited to a parallel processing implementation. The computation of the integrals in Eq. (2) can be carried out independently for each patch, and these computations involve only the image data local to that patch. Moreover, if the patches have the same volumes, regardless of shape, the number of computations for each patch is the same. Applying oversampling filters Eq. (3) to the images can be also done independently. This property could become important for fast registration when sophisticated nonlinear local distortions have to be modeled—for example, in 3D breast imaging. Acknowledgments We thank Drs W. Richter and L. Ryner, and Ms. C. Sweetland, of the Institute for Biodiagnostics for providing the 3D images. We also thank anonymous reviewers for constructive remarks and suggestions. References [1] Friston KJ, Williams S, Howard R, Frackowiak RSJ, Turner R. Movement-related effects in fMRI time-series. MRM 1996;35:346 – 55. [2] Woods RP, Grafton ST, Holmes CJ, Cherry SR, Mazziotta JC. Automated image registration. I. General methods and intrasubject, intramodality validation. J Comput Assist Tomo 1998;22(1):139 –52. [3] Horn BKP, Schunck BG. Determining optical flow. Artificial Intelligence 1981;17:185–203. [4] Alexander ME. Fast hierarchical noniterative registration algorithm. Int J Imaging Syst Technol 1999;10:242–57. [5] Gupta N, Kanal L. Gradient based image motion estimation without computing gradients. Int J Comput Vision 1997;22:81–101. [6] Simoncelli EP. Design of multidimensional derivative filters. Proc ICIP-94, Austin, TX, 1994. p. 790 – 4. [7] Daubechies I. Ten lectures on wavelets. Philadelphia, PA: Society for Industrial and Applied Mathematics, 1992. [8] 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 Comp Assist Tomogr 1995;19(2): 289 –96. [9] Press WH, Flanney BP, Teukolsky SA, Vetterling WT. Numerical recipes in C. 2nd ed. Cambridge, MA: Cambridge University Press, 1992. [10] Park SK, Schowengerdt RA. Image reconstruction by parametric cubic convolution. Computer Vision, Graphics, and Image Processing 1983;23:258 –72.