Interpolation artifacts in multimodality image registration based on ...
IEEE Copyright Notice
Copyright © 2003 IEEE. Reprinted from IEEE Transactions on Medical Imaging. This material is posted here with permission of the IEEE. Such permission of the IEEE does not in any way imply IEEE endorsement of any of the authors’ products or services. Internal or personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution must be obtained from the IEEE by writing to [email protected]. By choosing to view this document, you agree to all provisions of the copyright laws protecting it.
854
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
Interpolation Artifacts in Multimodality Image Registration Based on Maximization of Mutual Information Jeffrey Tsao, Member, IEEE
Abstract—Mutual information (MI) is an increasingly popular match metric for multimodality image registration. However, its value is affected by interpolation, which may limit registration accuracy. The purpose of this study was to characterize the artifacts from eight interpolators and to investigate efficient strategies to overcome these artifacts. The interpolators were: 1) nearest neighbor; 2) linear; 3) cubic Catmull–Rom; 4) Hamming-windowed sinc; 5) partial volume; 6) NN with jittered sampling (JIT); 7) NN with histogram blurring (BLUR); and 8) NN with JIT and BLUR. The impact of interpolation on MI was evaluated in two dimensions over different translational and rotational misregistration. Interpolation caused spurious fluctuations in MI whenever the voxel grids had coinciding periodicities and were nearly aligned. The artifacts did not lessen by using intensity interpolators with wider support (e.g., cubic Catmull–Rom, Hamming-windowed sinc). PV could lead to either arch artifacts or inverted-arch artifacts, depending on the relative voxel sizes. Several strategies reduced artifacts and improved registration robustness: JIT, BLUR, avoiding an extreme number of intensity bins, and resampling the images in a rotated orientation with different relative voxel sizes (e.g., 3). These findings also apply to related methods, including normalized MI, joint entropy, and Hill’s third moment. Index Terms—Generalized clustering framework, interpolation artifacts, multimodality image registration, mutual information.
I. INTRODUCTION
M
EDICAL images from different modalities can provide complementary information about anatomy or physiology in a synergistic manner [1]. Since images may be acquired in different poses, considerable effort has been placed on developing methods for multimodality image registration. Over the past decade, several promising methods have emerged [2], including Woods’ [3], moment-based [4], joint entropy [5], and mutual information (MI) [6], [7] methods. These methods are primarily identified by their choice of the match metric Manuscript received June 30, 2002; revised February 5, 2003. This work was supported in part by the National Institutes of Health (NIH) under Grant PHS 5 P41 RR05964 and Grant PHS RO1 CA 51430, in part by the National Science Foundation (NSF) under Science and Technology Center Award NSF DIR 89-20133, in part by the Illinois Department of Commerce and Community Affairs under Grant 91-82128, and in part by the Office of Naval Research (ONR) under Grant N00014-92-J-1160. The Associate Editor responsible for coordinating the review of this paper and recommending its publication was A. A. Amini. The author is currently with the Institute for Biomedical Engineering, Swiss Federal Institute of Technology, Zurich, Building ETF, Room C108, Sternwartstrasse 7, 8092 Zurich, Switzerland (e-mail [email protected]). Digital Object Identifier 10.1109/TMI.2003.815077
(also known as similarity measure), which makes them suitable for certain types of registration tasks. For example, the match metric of Woods’ method [3] assumes that when magnetic resonance (MR) and PET images are registered, the voxel intensities at corresponding positions conform to a many-to-one mapping. In comparison, the match metric of MI does not make any assumptions about the intensity mapping relationship [8] between the images. In general, a match metric with fewer assumptions is more practical, since it can be applied to a wider range of images. However, when the assumptions are valid, they can make the match metric more robust by providing additional constraints. Regardless of the choice of the match metric, these methods register images in a similar fashion [8], [9], which involves four steps: 1) images are superposed in a tentative pose relative to one another; 2) a joint histogram of voxel intensities is constructed from the overlap region of the two images; 3) a match metric is determined from the joint histogram; and 4) the match metric is optimized by reorienting the images iteratively. At the heart of these registration methods lies an interpolation algorithm, which is needed to estimate the voxel intensities at nongrid positions, whenever the voxel grids of the images are not in exact alignment. In general, interpolators with a larger support [4], [10]–[12] (e.g., cubic Catmull–Rom, sinc) lead to smaller intensity errors in the interpolated image. Hence, those interpolators have been argued to be better at achieving higher registration accuracy [4], [10], [13], [14]. However, a better interpolator may not lead to better registration per se. In most registration methods mentioned above [4]–[7], the voxel intensities are typically quantized into a fixed number of intensity bins before calculating the match metric. Therefore, the additional accuracy afforded by interpolators with larger support may not translate directly to improved accuracy of the match metric, due to the inherent rounding of the binning process. Moreover, it is important to note that the requirements for interpolation are not identical to those for registration. In interpolation, experimental nonidealities (e.g., noise and finite sampling) inevitably introduce errors to the interpolated voxel intensities. To minimize errors, it is often desirable to use data-consistent interpolators, so that the voxel intensity at any grid position does not change after interpolation. As a result, errors are only introduced when interpolating at nongrid positions. However, this implies that the amount of interpolation errors will vary depending on the extent of grid alignment. The variation may lead to artifactual fluctuations in the match metric, which may affect registration accuracy [15]–[18]. For some of the registration methods mentioned
0278-0062/03$17.00 © 2003 IEEE
TSAO: INTERPOLATION ARTIFACTS IN MUTUAL INFORMATION IMAGE REGISTRATION
above [4]–[7], their match metrics are particularly susceptible to such fluctuations. Thus, paradoxically, an interpolator that introduces a constant level of errors (i.e., a nondata-consistent interpolator) may in fact be more suitable than one that only introduces errors with misaligned voxel grids [16]. Therefore, for the purpose of registration, the effectiveness of an interpolator cannot be judged solely by the closeness of the interpolated image to the original image, but also by its effect on the match metric. In recent years, MI has become the match metric of choice due to its wide applicability and overall accuracy [19]. Therefore, several publications [15]–[18] have examined the impact of interpolation on MI. However, this impact has not been fully characterized yet. The purpose of this study was to extend the scope of previous work by investigating the patterns of interpolation artifacts for a wider range of interpolators, as well as efficient strategies to overcome these artifacts. Here, the term “interpolation artifact” refers to the effect of interpolation on the match metric. This is different from the conventional sense of the term, which refers to the effect of interpolation on image quality [20]. In this work, the effects of eight interpolation schemes were examined: 1) nearest neighbor (NN); 2) linear; 3) cubic Catmull–Rom; 4) Hamming-windowed sinc; 5) partial volume (PV); 6) NN with jittered sampling (JIT); 7) NN with histogram blurring (BLUR); and 8) NN with JIT and BLUR. This study also investigated the impact of using different numbers of intensity bins on the artifact pattern. Finally, the artifact-reducing effects of image rotation [16] and resampling [17], [18] were analyzed.
855
tional (around the image center) misregistration in the ranges of and , respectively. Such a plot was referred to as a registration curve. The ranges of misregistration were specified relative to the registered pose of the image contents. Thus, for example, a translational misregistration of three voxels for the rotated pair meant that the SPECT image was rotated 30 clockwise (to correct for misalignment introduced in creating the image) and then translated three voxels. Both the rotation and the translation were applied together as a single transformation to avoid accumulating errors in two interpolation steps. The transformations were applied to the SPECT image, which contained less signal energy in the high spatial frequencies, and was therefore less susceptible to degradation in image quality from interpolation. MI was expected to attain a maximum when the images were registered. To evaluate MI, the images were overlaid according to the specified pose relative to one another. A joint histogram was created by plotting an ( ) point for every pair of corresponding voxels in the overlap region, where was the voxel intensity in the first image, and was the interpolated voxel intensity in the second image. The voxel intensities were quantized into a fixed number of discrete bins (referred to as intensity bins hereafter). Thus, if the intensities of each image were quantized into intensity bins, the size of the joint histogram would be . The number of points in each histogram bin was normalized by . the total number of points to give a probabilistic measure MI was determined as
II. METHODS A. Human Head Images Registration experiments were conducted on four pairs of images, denoted as “aligned,” “rotated,” “128/129,” and ” (Fig. 1). The aligned pair consisted of a 128 128 axial “ -weighted MR image of a human head and a corresponding -ECD single photon emission computed 128 128 tomography (SPECT) image. The images were obtained using Fourier interpolation from a central axial slice of three-dimensional (3-D) images, which were preregistered using Woods’ method [3] and verified visually. The other three pairs of images were similar to the aligned pair, except that the SPECT image in the rotated pair was rotated 30 anti-clockwise around the center. The SPECT image in the 128/129 pair was rescaled to a grid size of 129 129 by reducing the voxel size by 128/129. pair was rescaled by increasing The SPECT image in the . These four image pairs were created prior the voxel size by to the registration experiments, using Fourier interpolation to resample the SPECT image in order to minimize degradation of the image contents. Even if some degradation was introduced in this step, it should have little influence on the results of the registration experiments, since the degradation did not vary during each experiment. Therefore, any fluctuation of the match metric resulted solely from subsequent operations. B. Evaluating Mutual Information For each pair of images, MI (i.e., the match metric) was determined as a function of translational (along the axis) and rota-
MI was evaluated for different numbers of intensity bins ( , 16, 32, 64, 128, 192, 256, and 512) and with various interpolation schemes. Interpolation artifacts were quantified by measuring the smoothness of each registration curve as follows (Fig. 2). First, the axis of each registration curve was normalized to a range of one. This allowed the interpolation artifacts to be quantified relative to the variation of the match metric. Then, the registration curve was subtracted from a smoothed version of itself, obtained by median filtering followed by convolution with a Hamming filter. Both filters had the same length, corresponding to a range of 1.5 voxels in translation, or 5 in rotation along axis. To avoid edge effects from filtering, the initial the and final points of the curve, equal to the filter length, were discarded. The difference between the registration curve and the smoothed version contained mostly interpolation artifacts (Fig. 2). Thus, smoothness was quantified as one over the root mean square (rms) of the difference curve. A registration curve with fewer interpolation artifacts was expected to be smoother, and would therefore attain a higher smoothness value. C. Interpolators Eight interpolation schemes were tested. For interpolators 1 to 4, the voxel intensity at a nongrid position was interpolated as a weighted sum of voxel intensities at neighboring grid positions. The total weights summed to unity, and the weighting
856
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
Fig. 1. Four image pairs used in registration experiments, denoted as “aligned,” “rotated,” “128/129,” and “=3.” Each pair consisted of an axial T -weighted MR head image and a corresponding Tc-ECD SPECT image. A different transformation (i.e., rotation or scaling) was applied to the SPECT image of each pair using Fourier interpolation. In this figure, the scaling of the SPECT voxel grid in the “128/129” and “=3” pairs is exaggerated for illustration purposes only.
was determined by the interpolation kernel (Fig. 3 and Table I), which was separable along the spatial dimensions. The width of the interpolation kernel determined the number of neighbors involved in the summation. Interpolators with a smaller support relied on fewer neighbors and were therefore faster. If the interpolation kernel extended beyond the image edge, the voxel intensity to be interpolated was considered to be outside the overlap region between the two images and was discarded
from the match metric. This was to avoid any confounding edge effects that might affect the interpolated intensity and alter the image contents. The number of overlapping voxels changed discretely, as the relative pose of the two images varied [14]. 1) Nearest Neighbor: The voxel at the closest grid position was assigned a weight of one, while all other voxels received zero weights.
TSAO: INTERPOLATION ARTIFACTS IN MUTUAL INFORMATION IMAGE REGISTRATION
857
Fig. 2. Quantification of interpolation artifacts. Left: registration curve showing MI as a function of translational misregistration. Center: smoothed version of registration curve. Right: difference between original registration curve and the smoothed version, revealing mainly interpolation artifacts. Smoothness was defined as one over the rms of the difference curve.
Fig. 3. Interpolation kernels in one dimension. The width of the support is shown below.
TABLE I EXPRESSION FOR INTERPOLATOR IN ONE DIMENSION. THE POINT TO INTERPOLATE IS AT POSITION x, BETWEEN 0 AND 1, WITH NEIGHBORING GRID POINTS LISTED IN ORDER AS . . . ; n ; n ; n ; n ; n ; n ; . . . THE POINT TO INTERPOLATE LIES BETWEEN n AND n , AND IS EQUIVALENT TO n WHEN x = 0, AND TO n WHEN x = 1.
2) Linear: The interpolation kernel was triangular in shape along each spatial dimension, with a support width of two voxels.
3) Cubic Catmull–Rom: The interpolation kernel was determined by the Catmull–Rom [21] cubic spline, with a support width of four voxels.
858
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
6
Fig. 4. Translation results for all interpolators and image pairs, with MI plotted as a function of translational misregistration (within a range of 10 voxels along the y axis). Each column of graphs corresponds to a different interpolator, while each row corresponds to a different image pair. The scales of the graph axes are shown in the lower left panel. In each panel, the curves from bottom to top correspond to 8, 16, 32, 64, 128, 192, 256, and 512 intensity bins, respectively.
4) Hamming-Sinc: The interpolation kernel along each spatial dimension was a sinc function multiplied by a six-voxel-wide Hamming window. The kernel was normalized to ensure that the total weights summed to unity [22]. 5) Partial Volume: This interpolator, as introduced by Collignon et al. [7], [23], is not an interpolator in the traditional sense because it does not estimate the intensities of voxels at nongrid positions. Instead, each voxel at a nongrid position is treated as though it consists of several fractional voxels at neighboring grid positions. Since the fractional voxels are at grid positions, their intensities can be evaluated directly without interpolation. The fractions sum to unity and were determined according to proportions from linear interpolation [7], [23]. In general, the fractions may be determined according to proportions from other interpolators as well [24]. 6) Nearest Neighbor With Jittered Sampling: For every voxel to be interpolated, its coordinates were jittered by , adding a normally distributed random offset ( ). The intensity at the jittered coordinates was interpolated by NN interpolation. 7) Nearest Neighbor With Histogram Blurring: Following NN interpolation, the joint histogram was blurred by convolving it with a filter (coefficients 1/16, 4/16, 6/16, 4/16, and 1/16) along each axis of the histogram. 8) Nearest Neighbor With Jittered Sampling and Histogram Blurring: This interpolator is a combination of interpolators 6 and 7. Therefore, jittering was applied to the coordinates of
every voxel to be interpolated. Then, the joint histogram was blurred by convolution, as described earlier. III. RESULTS Figs. 4 and 5 summarize the results for all interpolators (columns) and image pairs (rows). Each graph panel shows MI as a function of translational (Fig. 4) or rotational (Fig. 5) misregistration for different numbers of intensity bins. In each panel, the curves from bottom to top correspond to increasing number of intensity bins. Fig. 6 summarizes the quantification of interpolation artifacts for the translation (column 1) and rotation (column 2) experiments, and for all image pairs (rows). Each graph panel shows the smoothness of the registration curves as a function of the number of intensity bins in a semi natural-log scale for all interpolators. In general, interpolation artifacts were most evident in the translation results of the aligned pair (Fig. 4, row 1). The artifact pattern was characteristic for each interpolator, and it varied with the number of intensity bins. For most interpolators, the artifacts worsened at the high and/or low ends of the intensity bin range. This trend was reflected as a general decrease in smoothness toward the right and/or left sides of each panel in Fig. 6, respectively. In most cases, increasing the number of intensity bins excessively also had the deleterious effect of changing the match metric around 0 from the global maximum to a local minimum
TSAO: INTERPOLATION ARTIFACTS IN MUTUAL INFORMATION IMAGE REGISTRATION
859
6
Fig. 5. Rotation results for all interpolators and image pairs, with MI plotted as a function of rotational misregistration (within a range of 180 around the image center). Each column of graphs corresponds to a different interpolator, while each row corresponds to a different image pair. The scales of the graph axes are shown in the lower left panel. In each panel, the curves from bottom to top correspond to 8, 16, 32, 64, 128, 192, 256, and 512 intensity bins, respectively.
(e.g., See curves from bottom to top in Fig. 5, row 1, column 2). This would lead to improper registration since the match metric no longer attained a maximum when the images were registered. NN interpolation was found to produce stairs-like artifacts in the translation results of the aligned pair (Fig. 4, row 1, column 1). The stairs-like discontinuities occurred at every half-voxel step, coinciding with the discontinuities of the interpolation kernel (Fig. 3). Additionally, there were subtle discontinuities at every integer-voxel step, which occurred when voxels close to the edges of the overlap region either fell into or outside the overlap region, thus changing the voxels contributing to the match metric. The artifacts did not vary appreciably with different number of intensity bins, as reflected by a smoothness curve that was almost flat (Fig. 6, row 1, column 1). The stairs-like artifacts were reduced substantially for the other image pairs (Fig. 4, rows 2 – 4, column 1), although they remained noticeable for the 128/129 pair. This reduction in artifact was reflected quantitatively by an increase in smoothness (Fig. 6, rows 2–4, column 1). These artifacts were eliminated in the rotation results of all image pairs (Fig. 5, column 1). The artifacts of linear, cubic Catmull–Rom and Hamming-sinc interpolation (Figs. 4 and 5, columns 2 – 4, respectively) were similar in nature. In the translation results of the aligned pair (Fig. 4, row 1, columns 2 – 4), the artifacts approached a stairs-like pattern at the low end of the intensity bin range, with the discontinuities occurring at every integer-voxel step. These stairs-like artifacts were accompanied
by high-frequency fluctuations for cubic Catmull–Rom and Hamming-sinc. By increasing the number of intensity bins, the stairs-like artifacts diminished, but a different type of artifacts with an arch pattern began to emerge. The arch artifacts were accentuated with increasing number of intensity bins, and they bifurcated successively into finer arches. These artifacts were mostly eliminated for the other image pairs, except at the highest numbers of intensity bins (Fig. 4, rows 2 – 4, columns 2 – 4). Artifacts were not seen in the rotation results of all image pairs, except for either small dips or spikes at rotation angles that aligned the voxel grids (Fig. 5, columns 2 – 4). Quantitatively, the smoothness of the registration curves did not differ considerably among the three interpolators (linear, cubic Catmull–Rom and Hamming-sinc), and their smoothness curves overlapped in most cases (Fig. 6). Whenever the smoothness value differed among the three interpolators, it tended to be lower for cubic Catmull–Rom and Hamming-sinc interpolation than for linear interpolation, indicating that the interpolators with wider support did not lead to smoother registration curves. PV interpolation produced inverted arches in the translation results of the aligned pair (Fig. 4, row 1, column 5). The inverted arches were accentuated when there were more intensity bins, as reflected by a drop in smoothness (Fig. 6, row 1, column 1). These inverted arches changed to arches for the 128/129 pair (Fig. 4, row 3, column 5), but were not seen for the rotated and pairs (Fig. 4, rows 2 and 4, column 5), as reflected by the the increased smoothness. In the rotation results (Fig. 5, column
860
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
Fig. 6. Smoothness of the registration curves as a function of the number of intensity bins in a semi natural-log scale. Columns 1 and 2 correspond to the translation and rotation experiments, respectively. Each row corresponds to a different image pair. In each panel, the different curves correspond to the smoothness of the registration curves for different interpolators. Note that the scale on the Y axis is different for each panel.
5), the only artifacts seen for all image pairs consisted of either dips or spikes at rotation angles that aligned the voxel grids. For interpolators that included JIT (i.e., interpolators 6 and 8) (Figs. 4 and 5, columns 6 and 8), the artifacts consisted of small stochastic fluctuations only. Since the jittering degraded the similarity between the images, it decreased the value of the match metric slightly, compared with interpolators without JIT. However, this decrease should have no effect on registration accuracy since it did not depend on the alignment of the image contents. In general, the smoothness values of these interpolators varied over a narrower range than other interpolators (Fig. 6),
due to the absence of any severe interpolation artifacts. At times, these interpolators considerably outperformed the others, such as in the translation results of the aligned pair (Fig. 6, row 1, column 1). However, in situations where the artifacts of other interpolators were suppressed (e.g., in the rotation results of pair, Fig. 6, row 4, column 2), the slight stochastic the fluctuations resulted in a slightly lower overall smoothness. For interpolators that included BLUR (i.e., interpolators 7 and 8) (Figs. 4 and 5, columns 7 and 8), the artifacts were similar to those without BLUR. Nevertheless, BLUR generally improved smoothness (Fig. 6). The blurring also decreased the overall
TSAO: INTERPOLATION ARTIFACTS IN MUTUAL INFORMATION IMAGE REGISTRATION
value of the match metric. The main impact of BLUR was to maintain the maximum of the match metric around 0 when there was a large number of intensity bins (Fig. 5). In other words, the registration curve maintained its shape over widely different number of intensity bins. This was reflected as the increased flatness of the smoothness curves (Fig. 6). IV. DISCUSSION Interpolation artifacts can occur if the voxel grids of the images have coinciding periodicities. Since these artifacts may cause severe fluctuations in the match metric, they can limit the achievable registration accuracy and undermine the robustness of the match metric. The present results showed that the artifacts were reduced substantially when the periodicities of the voxel grids were abolished with strategies such as JIT, or when one image was rotated or resampled relative to the other image. A. Number of Intensity Bins In the present study, the artifact patterns were characterized over widely different number of intensity bins. Overall, the artifact patterns were found to change smoothly with the number of intensity bins (Figs. 4 and 5). Although extreme choices should be avoided, the precise choice for the number of intensity bins may not be critical. In general, it is preferable to err on the low side, since having fewer intensity bins allows for a faster processing time and a smaller size for the joint histogram. Otherwise, if the joint histogram is large, each histogram bin will be populated by a small number of points, and the match metric may suffer from poor statistics. This is similar to the adverse effect seen with image subsampling [17], [18], since both procedures lead to fewer points occupying each histogram bin [25]. If the number of intensity bins is chosen to be even higher, most of the points in the joint histogram may fail to aggregate sufficiently to form clusters. Thus, the only prominent clusters remaining will correspond to the background voxels, due to the sheer number of these voxels [26]. As a result, the match metric will be increasingly influenced by the nonspecific matching of background voxels, and registration accuracy will degrade [26]. This problem was encountered in the present study for most interpolators. When there were too many intensity bins, the match metric at 0 changed from the global maximum to a local minimum, thus leading to improper registration (Fig. 5). This problem may be of particular concern in multiresolution [14], [27] or multiscale [25] matching, since the number of voxels may be small at certain levels of matching, so the number of intensity bins may need to be adjusted accordingly [14], [27]. B. Classification of Interpolators The interpolators examined in this study can be divided into two groups: intensity interpolators (linear, cubic Catmull–Rom, and Hamming-sinc), and nonintensity interpolators (the remaining ones). The main difference between the two groups is that intensity interpolators may introduce new intensity values that are not present in the original image [7], [23]. For example, a binary-valued image may contain fractional intensity values after interpolation. In contrast, nonintensity interpolators only return intensity values from the original image. Therefore, the
861
mechanism of artifact generation is fundamentally different between the two groups. Intensity interpolators generate artifacts by modifying voxel intensities, while nonintensity interpolators generate artifacts by changing the way the match metric is evaluated. As a result, the artifacts from intensity interpolators are inherently more dependent on the image contents, thus making them somewhat difficult to analyze. C. Linear, Cubic Catmull–Rom, & Hamming-Sinc Interpolation Usually, the performance of an intensity interpolator is characterized by its frequency response (or power spectrum) [20]. While this description is suited to determining the amount of interpolation error, it is inadequate in describing the impact of interpolation on the match metric. This is because the match metric is a highly nonlinear function of the voxel intensities. This inadequacy was shown in the present paper by the somewhat unexpected finding that the switch from linear to cubic Catmull–Rom to Hamming-sinc interpolation led to no improvement in the smoothness of the registration curves. Instead, it even worsened the smoothness some of the times. These results suggest that registration may not benefit from using intensity interpolators with wider support, despite their improved frequency response. This finding is important from a practical standpoint, since interpolators with wider support are also computationally more demanding. The present results also suggest that intensity interpolators may not be preferable to nonintensity interpolators, since the former did not appear to be any less susceptible to producing artifacts than the latter. The intensity interpolators being studied (i.e., linear, cubic Catmull–Rom, Hamming-sinc) were found to produce stairs-like artifacts with a small number of intensity bins and arch artifacts with more intensity bins. The stairs-like artifacts were expected since these interpolators would behave in a similar manner to NN interpolation, if the intensities were discretized to very few levels. The stairs-like discontinuities occurred at every integer-voxel step when the interpolation kernel switched to a different set of neighboring voxels for interpolation. For the arch artifacts at higher numbers of intensity bins, their exact shapes vary with image contents. For example, by adding 10% noise pixel-wise (without filtering) to the SPECT image, the artifacts become more pronounced, adopting a rounded appearance (Fig. 7). This agrees with the findings of Pluim et al. [17], [18], who attributed the origin of these artifacts in linear interpolation to the intensity averaging effect of the interpolator. The present study showed that further increases in the number of intensity bins led to bifurcations of these arch artifacts. The onset of arch bifurcation was delayed to a higher number of intensity bins for interpolators with a larger support (i.e., cubic Catmull–Rom, Hamming-sinc). While intensity averaging is a likely cause of these arch artifacts, their complicated shapes (Fig. 4, columns 2–4) suggest that there may be additional contributing mechanisms, which deserve further investigation in the future. D. Partial Volume Interpolation The present study showed that PV interpolation could lead to inverted-arch artifacts (Fig. 4, row 1, column 5) as well as
862
Fig. 7. Artifact patterns of linear, cubic Catmull–Rom, and Hamming-sinc interpolation, after adding 10% noise pixel-wise to the SPECT image of the aligned pair. The noise level indicates the standard deviation of noise relative to the maximum intensity. The noisy SPECT image is shown on the upper left. The registration curves show MI as a function of translational misregistration, for the aligned pair with 64 intensity bins.
arch artifacts (Fig. 4, row 3, column 5). Both types of artifacts originated from the splitting of misaligned voxels into fractional parts, since the fractional parts contributed less in total to the match metric than the whole voxel [8], [28]. Fig. 8 illustrates the origin of these artifacts. Column 1 magnifies one of the registration curves from column 5 of Fig. 4. The grid alignment is illustrated for three particular points of the registration curves ( , , and ) in columns 2–4, respectively. The plots show the distance between each voxel row in the SPECT image and the nearest voxel row in the MR image. This distance is referred to as fractional grid misalignment. For example, a fractional grid misalignment of zero means that the voxel rows in both images are exactly aligned, while a fractional grid mismeans that the SPECT voxel row is right alignment of in-between two nearest voxel rows in the MR image. The inverted-arch artifacts were manifested in translation when the ratio between the voxel sizes of the two images was either one (i.e., both images have identical voxel sizes) or a rational fraction of small numbers (e.g., 3/5) [15]–[18]. The bottom of the inverted arches corresponded to the majority of the voxels being misaligned (Fig. 8, row 1, column 4). In contrast, the arch artifacts were manifested when there were only slight mismatches between the voxel sizes, as found in the 128/129 pair (Fig. 4, row 3, column 5). Due to the small voxel size difference, the amount of grid misalignment differed slightly for each voxel row (Fig. 8, row 2). Thus, some portions of the image had smaller fractional grid misalignment (i.e., were more aligned) than others. The more aligned portions would have higher contribution to the
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
match metric. For example, when the 128/129 pair was registered (Fig. 8, row 2, column 2), voxel rows close to the edges of the image (row 0 or 128; see black arrows) were more aligned than those close to the center (e.g., row 64; see open arrows). Since the voxel rows at the edges corresponded primarily to the image background, their higher contribution to the match metric would lead to an overall decrease in the match metric, resulting in the dips between the arches. The top of the arches occurred at every half-voxel step, when the voxel rows in the image center came into grid alignment, thus contributing more to the match metric (See open arrows in Fig. 8, row 2, column 4). Accordingly, if the image contents were shifted by half a field of view, the top of the arches would occur at every integer-voxel step instead (data not shown). These arch and inverted-arch artifacts can be eliminated by choosing the voxel sizes judiciously. For example, a salient choice was for the voxel sizes of one image to be times the voxel sizes of the other image, since the resultant amount of grid misalignment was considerably different from row to row (Fig. 8, row 3). Thus, all portions of the voxel grids were misaligned by the same amount on average, regardless of the translation between the images (Fig. 8, row 3, columns 2 – 4). With being an irrational number, the amount of misalignment also never repeated itself exactly from row to row, thereby reducing any periodicity. The artifacts of PV interpolation were found to reduce dramatically in the translation results when the relative voxel size was changed from 128/129 (Fig. 4, row 3, (Fig. 4, row 4, column 2). column 2) to E. Jittered Sampling Considerable artifact reduction was achieved in all cases with JIT, which abolished the periodicities of the voxel grids and eliminated the associated artifacts (see Likar and Pernuˇs [25] for a related approach). Even though JIT in the present implementation deliberately introduced error in the interpolation process, it proved to be an effective interpolator for registration, since it introduced approximately the same amount of error regardless of the extent of grid alignment. F. Histogram Blurring Improved registration robustness was found with BLUR, which maintained the global maximum of the match metric at 0 even when there was a large number of intensity bins. By applying BLUR, the points in each histogram bin were spread into neighboring bins, which allowed the amount of clustering in the joint histogram to be estimated accurately, even when there were very few points per histogram bin. In fact, BLUR can be considered an approximation of Parzen windowing [29] and gridding [30], which are two mathematically similar techniques of estimation from discrete samples. In the present context, the blurring kernel served as the window function of Parzen windowing, or the convolution function of gridding. The approximation became more accurate with more intensity bins, due to smaller discretization error, albeit at the expense of increased computation. In practice, for a moderate number
TSAO: INTERPOLATION ARTIFACTS IN MUTUAL INFORMATION IMAGE REGISTRATION
863
Fig. 8. Artifact pattern of PV interpolation depends on the ratio between the voxel sizes of the images. Rows 1 – 3 show the results of the aligned, 128/129, and =3 pairs, with voxel ratios of 1/1, 128/129 and =3, respectively. Column 1 shows a magnified version of the respective registration curves, corresponding to column 5 of Fig. 4. The artifact patterns arose from different grid alignments, as the two images were translated from one another. The alignments at three particular points ( , , and ) are illustrated in columns 2 – 4, respectively. The plots show the distance between each voxel row in the SPECT image and the nearest voxel row in the MR image. This distance is referred to as fractional grid misalignment. Arrows point to relevant features. See text for details.
of intensity bins, BLUR provides an efficient approximation of Parzen windowing to improve registration accuracy. G. Related Match Metrics and Generalization Although this study focused on MI only, the observations can be generalized to other methods such as moment-based [4], and joint entropy [5] methods, since the match metrics of these methods can be expressed in a similar mathematical form as MI [8], [9], [28]. The present observations also apply to normalized .A MI [31], which can be expressed as (MI)/(joint entropy) natural question is whether the present results also apply to other image pairs. Results from previous studies [8], [16], [26], [32] and from other groups [15], [17], [18], [25] suggest that while the exact shapes of the artifacts are dependent on the image contents (especially for intensity interpolators), the general artifact patterns are similar over a wide range of images, including images of different modalities and size scales, and in two and three dimensions (Fig. 9). This is perhaps not too surprising since the artifacts are primarily related to implementation issues of the match metric, which are independent of the image contents. H. Summary In conclusion, registration using MI or similar match metrics can be affected severely by interpolation artifacts. Nevertheless,
Fig. 9. Translation results for an aligned pair of 3-D SPECT and MRI images with 128 intensity bins and identical voxel sizes. MI is plotted as a function of translational misregistration (within a range of 10 voxels along the y axis). Results for the different interpolators are staggered along the vertical axis for illustration purposes.
6
the artifacts can be reduced significantly if the voxel grids are misaligned. In most registration tasks, especially in three dimensions, the voxel grids are unlikely to be in exact alignment, so
864
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
interpolation artifacts are unlikely to be a major source of registration error [15]. In those cases, it may suffice to use NN interpolation, thus allowing a considerable reduction in processing time. Even in those cases, it may be worthwhile to incorporate BLUR, which is computationally inexpensive for a moderate number of histogram bins, and was found to always improve the smoothness of the registration curves. In cases where the voxel grids are almost in alignment (e.g., when a head holder is used to minimize differences in head positioning, such as in fMRI acquisitions [32], [33]), special attention must be paid to reduce interpolation artifacts in order to ensure registration accuracy [15]. This work showed that the following strategies could reduce artifacts considerably and improve registration robustness: 1) avoiding an extreme number of intensity bins; 2) resampling the images in a rotated orientation with unequal voxel sizes (e.g., voxel sizes of one image being times those of the other image); 3) JIT; and 4) BLUR. The first two strategies constitute precomputation steps, and do not increase the amount of processing at each iteration of the registration algorithm. Applying these two strategies together often suffices to eliminate most interpolation artifacts, so NN interpolation combined with BLUR can be used for fast performance, while problems from interpolation artifacts and/or an inappropriately chosen number of intensity bins can be avoided. The rare exception is in pathological situations during the registration process that happen to realign the voxel grids. To ensure small interpolation artifacts at all time, JIT can be incorporated with minimal increase in processing time.
ACKNOWLEDGMENT The author would like to thank Dr. P. C. Lauterbur and Dr. S. R. Teixeira for helpful comments on the manuscript, and Dr. M. Ichise for the SPECT image.
REFERENCES [1] M. W. Vannier and D. E. Gayou, “Automated registration of multimodality images,” Radiology, vol. 169, pp. 860–861, 1988. [2] D. L. Hill, P. G. Batchelor, M. Holden, and D. J. Hawkes, “Medical image registration,” Phys. Med. Biol., vol. 46, pp. R1–45, 2001. [3] R. P. Woods, J. C. Mazziotta, and S. R. Cherry, “MRI-PET registration with automated algorithm,” J. Comput. Assist. Tomogr., vol. 17, pp. 536–546, 1993. [4] D. L. G. Hill, C. Studholme, and D. J. Hawkes, “Voxel similarity measures for automated image registration,” Proc. Vis. Biomed. Comput., pp. 205–216, 1994. [5] A. Collignon, D. Vandermeulen, P. Suetens, and G. Marchal, “3 D multimodality medical image registration using feature space clustering,” in Proc. Int. Conf. Computer Vision Virtual Reality Robot Medicine, 1995, pp. 195–204. [6] W. M. Wells 3rd, P. Viola, H. Atsumi, S. Nakajima, and R. Kikinis, “Multi-modal volume registration by maximization of mutual information,” Med. Image Anal., vol. 1, pp. 35–51, 1996. [7] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, and P. Suetens, “Multimodality image registration by maximization of mutual information,” IEEE Trans Med Imag., vol. 16, pp. 187–198, Apr. 1997. [8] J. Tsao and P. C. Lauterbur, “Generalized clustering-based image registration for multi-modality images,” Proc. IEEE Engineering Medicine Bioloy Soc, pp. 667–670, 1998.
[9] M. Bro-Nielsen, “Rigid registration of CT, MR and cryosection images using a GLCM framework,” in Proc CVRMed/MRCAS, 1997, pp. 171–180. [10] J. V. Hajnal, N. Saeed, E. J. Soar, A. Oatridge, I. R. Young, and G. M. Bydder, “A registration and interpolation procedure for subvoxel matching of serially acquired MR images,” J. Comput. Assist. Tomogr., vol. 19, pp. 289–296, 1995. [11] J. L. Ostuni, A. K. Santha, V. S. Mattay, D. R. Weinberger, R. L. Levin, and J. A. Frank, “Analysis of interpolation effects in the reslicing of functional MR images,” J. Comput. Assist. Tomogr., vol. 21, pp. 803–810, 1997. [12] T. M. Lehmann, C. Gonner, and K. Spitzer, “Survey: Interpolation methods in medical image processing,” IEEE Trans. Med. Imag., vol. 18, pp. 1049–1075, Nov. 1999. [13] W. F. Eddy, M. Fitzgerald, and D. C. Noll, “Improved image registration by using Fourier interpolation,” Magn. Reson. Med., vol. 36, pp. 923–931, 1996. [14] P. Thévenaz and M. Unser, “Optimization of mutual information for multiresolution image registration,” IEEE Trans. Image Processing, vol. 9, pp. 2083–2099, Dec. 2000. [15] F. Maes, “Segmentation and Registration of Multimodal Medical Images: From Theory, Implementation and Validation to a Useful Tool in Clinical Practice,” Catholic University of Leuven, Leuven , Belgium, 1998. [16] J. Tsao, “Efficient interpolation for clustering-based multimodality image registration,” in Proc. Int. Soc. Magnetic Resonence Medicine, 1999, p. 2195. [17] J. P. W. Pluim, J. B. A. Maintz, and M. A. Viergever, “Mutual information matching and interpolation artifacts,” Proc SPIE, pp. 56–65, 1999. , “Interpolation artefacts in mutual information-based image regis[18] tration,” Comput. Vis. Image Understanding, vol. 77, pp. 211–232, 2000. [19] J. West et al., “Comparison and evaluation of retrospective intermodality brain image registration techniques,” J. Comput. Assist. Tomogr., vol. 21, pp. 554–566, 1997. [20] P. Thévenaz, T. Blu, and M. Unser, “Interpolation revisited,” IEEE Trans. Med. Imag., vol. 19, pp. 739–758, 2000. [21] E. Catmull and R. Rom, “A class of local interpolationg splines,” in Computer Aided Geometric Design, R. E. Barnhill and R. F. Riesenfled, Eds. New York: Academic Press, 1974. [22] N. A. Thacker, A. Jackson, D. Moriarty, and E. Vokurka, “Improved quality of re-sliced MR images using re-normalized sinc interpolation,” J. Magn. Reson. Imaging, vol. 10, pp. 582–588, 1999. [23] A. Collignon, F. Maes, D. Delaere, D. Vandermeulen, P. Suetens, and G. Marchal, “Automated multimodality medical image registration using information theory,” in Proc. Int. Conf. Information Processing Medical Imaging: Computational Imaging and Vision, 1995, pp. 263–274. [24] H. M. Chen and P. K. Varshney, “Generalized partial volume interpolation for image registration based on mutual information,” presented at the Proc. IEEE Signal Processing Soc., Western New York, Image Processing Workshop, Univ. Rochester, Ctr. Optoelectron. Imag., 2000. [25] B. Likar and F. Pernuˇs, “A hierarchical approach to elastic registration based on mutual information,” Image Vis. Computing, vol. 19, pp. 33–44, 2001. [26] J. Tsao, “Multimodality image registration by gradient enhanced generalized clustering,” in Proc. Int. Soc. Magn. Reson. Med., 1999, p. 125. [27] M. Jenkinson and S. Smith, “A global optimization method for robust affine registration of brain images,” Med. Image Anal., vol. 5, pp. 143–156, 2001. [28] T. M. Buzug and J. Weese, “Image registration for DSA quality enhancement,” Comput. Med. Imag. Graph., vol. 22, pp. 103–113, 1998. [29] E. Parzen, “On estimation of a probability density function and the mode,” Ann. Math. Statist., vol. 33, pp. 1065–1076, 1962. [30] J. D. O’Sullivan, “A fast sinc function gridding algorithm for Fourier inversion in computer tomography,” IEEE Trans. Med. Imag., vol. MI-4, pp. 200–207, 1985. [31] C. Studholme, D. L. G. Hill, and D. J. Hawkes, “An overlap invariant entropy measure of 3 D medical image alignment,” Pattern Recogn., vol. 32, pp. 71–86, 1999. [32] J. Tsao, J. B. M. Goense, and P. C. Lauterbur, “fMRI image registration using total local density of intensity mapping plot,” in Proc. Int. Soc. Magnetic Resonance Medicine, 1998, p. 2135. [33] B. Kim, J. L. Boes, P. H. Bland, T. L. Chenevert, and C. R. Meyer, “Motion correction in fMRI via registration of individual slices into an anatomical volume,” Magn. Reson. Med., vol. 41, pp. 964–972, 1999.
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
1183
Correspondence________________________________________________________________________ Correction to “Interpolation Artifacts in Multimodality Image Registration Based on Maximization of Mutual Information”
REFERENCES [1] J. Tsao, “Interpolation artifacts in multimodality image registration based on maximization of mutual information,” IEEE Trans. Med. Imag, vol. 22, pp. 854–864, July 2003.
Jeffrey Tsao In the above paper [1], some of the information in Fig 5 was missing. The figure should have appeared as shown here. Manuscript received July 30, 2003. The author is with the Institute for Biomedical Engineering, Swiss Federal Institute of Technology, Zurich, Building ETF, Room C108, Sternwartstrasse 7, 8092 Zurich, Switzerland (e-mail [email protected]). Digital Object Identifier 10.1109/TMI.2003.817660
6
Fig. 5. Rotation results for all interpolators and image pairs, with MI plotted as a function of rotational misregistration (within a range of 180 around the image center). Each column of graphs corresponds to a different interpolator, while each row corresponds to a different image pair. The scales of the graph axes are shown in the lower left panel. In each panel, the curves from bottom to top correspond to 8, 16, 32, 64, 128, 192, 256, and 512 intensity bins, respectively. 0278-0062/03$17.00 © 2003 IEEE
Copyright © 2003 IEEE. Reprinted from IEEE Transactions on Medical Imaging. This material is posted here with permission of the IEEE. Such permission of the IEEE does not in any way imply IEEE endorsement of any of the authors’ products or services. Internal or personal use of this material is permitted. However, permission to reprint/republish this material for advertising or promotional purposes or for creating new collective works for resale or redistribution must be obtained from the IEEE by writing to [email protected]. By choosing to view this document, you agree to all provisions of the copyright laws protecting it.
854
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
Interpolation Artifacts in Multimodality Image Registration Based on Maximization of Mutual Information Jeffrey Tsao, Member, IEEE
Abstract—Mutual information (MI) is an increasingly popular match metric for multimodality image registration. However, its value is affected by interpolation, which may limit registration accuracy. The purpose of this study was to characterize the artifacts from eight interpolators and to investigate efficient strategies to overcome these artifacts. The interpolators were: 1) nearest neighbor; 2) linear; 3) cubic Catmull–Rom; 4) Hamming-windowed sinc; 5) partial volume; 6) NN with jittered sampling (JIT); 7) NN with histogram blurring (BLUR); and 8) NN with JIT and BLUR. The impact of interpolation on MI was evaluated in two dimensions over different translational and rotational misregistration. Interpolation caused spurious fluctuations in MI whenever the voxel grids had coinciding periodicities and were nearly aligned. The artifacts did not lessen by using intensity interpolators with wider support (e.g., cubic Catmull–Rom, Hamming-windowed sinc). PV could lead to either arch artifacts or inverted-arch artifacts, depending on the relative voxel sizes. Several strategies reduced artifacts and improved registration robustness: JIT, BLUR, avoiding an extreme number of intensity bins, and resampling the images in a rotated orientation with different relative voxel sizes (e.g., 3). These findings also apply to related methods, including normalized MI, joint entropy, and Hill’s third moment. Index Terms—Generalized clustering framework, interpolation artifacts, multimodality image registration, mutual information.
I. INTRODUCTION
M
EDICAL images from different modalities can provide complementary information about anatomy or physiology in a synergistic manner [1]. Since images may be acquired in different poses, considerable effort has been placed on developing methods for multimodality image registration. Over the past decade, several promising methods have emerged [2], including Woods’ [3], moment-based [4], joint entropy [5], and mutual information (MI) [6], [7] methods. These methods are primarily identified by their choice of the match metric Manuscript received June 30, 2002; revised February 5, 2003. This work was supported in part by the National Institutes of Health (NIH) under Grant PHS 5 P41 RR05964 and Grant PHS RO1 CA 51430, in part by the National Science Foundation (NSF) under Science and Technology Center Award NSF DIR 89-20133, in part by the Illinois Department of Commerce and Community Affairs under Grant 91-82128, and in part by the Office of Naval Research (ONR) under Grant N00014-92-J-1160. The Associate Editor responsible for coordinating the review of this paper and recommending its publication was A. A. Amini. The author is currently with the Institute for Biomedical Engineering, Swiss Federal Institute of Technology, Zurich, Building ETF, Room C108, Sternwartstrasse 7, 8092 Zurich, Switzerland (e-mail [email protected]). Digital Object Identifier 10.1109/TMI.2003.815077
(also known as similarity measure), which makes them suitable for certain types of registration tasks. For example, the match metric of Woods’ method [3] assumes that when magnetic resonance (MR) and PET images are registered, the voxel intensities at corresponding positions conform to a many-to-one mapping. In comparison, the match metric of MI does not make any assumptions about the intensity mapping relationship [8] between the images. In general, a match metric with fewer assumptions is more practical, since it can be applied to a wider range of images. However, when the assumptions are valid, they can make the match metric more robust by providing additional constraints. Regardless of the choice of the match metric, these methods register images in a similar fashion [8], [9], which involves four steps: 1) images are superposed in a tentative pose relative to one another; 2) a joint histogram of voxel intensities is constructed from the overlap region of the two images; 3) a match metric is determined from the joint histogram; and 4) the match metric is optimized by reorienting the images iteratively. At the heart of these registration methods lies an interpolation algorithm, which is needed to estimate the voxel intensities at nongrid positions, whenever the voxel grids of the images are not in exact alignment. In general, interpolators with a larger support [4], [10]–[12] (e.g., cubic Catmull–Rom, sinc) lead to smaller intensity errors in the interpolated image. Hence, those interpolators have been argued to be better at achieving higher registration accuracy [4], [10], [13], [14]. However, a better interpolator may not lead to better registration per se. In most registration methods mentioned above [4]–[7], the voxel intensities are typically quantized into a fixed number of intensity bins before calculating the match metric. Therefore, the additional accuracy afforded by interpolators with larger support may not translate directly to improved accuracy of the match metric, due to the inherent rounding of the binning process. Moreover, it is important to note that the requirements for interpolation are not identical to those for registration. In interpolation, experimental nonidealities (e.g., noise and finite sampling) inevitably introduce errors to the interpolated voxel intensities. To minimize errors, it is often desirable to use data-consistent interpolators, so that the voxel intensity at any grid position does not change after interpolation. As a result, errors are only introduced when interpolating at nongrid positions. However, this implies that the amount of interpolation errors will vary depending on the extent of grid alignment. The variation may lead to artifactual fluctuations in the match metric, which may affect registration accuracy [15]–[18]. For some of the registration methods mentioned
0278-0062/03$17.00 © 2003 IEEE
TSAO: INTERPOLATION ARTIFACTS IN MUTUAL INFORMATION IMAGE REGISTRATION
above [4]–[7], their match metrics are particularly susceptible to such fluctuations. Thus, paradoxically, an interpolator that introduces a constant level of errors (i.e., a nondata-consistent interpolator) may in fact be more suitable than one that only introduces errors with misaligned voxel grids [16]. Therefore, for the purpose of registration, the effectiveness of an interpolator cannot be judged solely by the closeness of the interpolated image to the original image, but also by its effect on the match metric. In recent years, MI has become the match metric of choice due to its wide applicability and overall accuracy [19]. Therefore, several publications [15]–[18] have examined the impact of interpolation on MI. However, this impact has not been fully characterized yet. The purpose of this study was to extend the scope of previous work by investigating the patterns of interpolation artifacts for a wider range of interpolators, as well as efficient strategies to overcome these artifacts. Here, the term “interpolation artifact” refers to the effect of interpolation on the match metric. This is different from the conventional sense of the term, which refers to the effect of interpolation on image quality [20]. In this work, the effects of eight interpolation schemes were examined: 1) nearest neighbor (NN); 2) linear; 3) cubic Catmull–Rom; 4) Hamming-windowed sinc; 5) partial volume (PV); 6) NN with jittered sampling (JIT); 7) NN with histogram blurring (BLUR); and 8) NN with JIT and BLUR. This study also investigated the impact of using different numbers of intensity bins on the artifact pattern. Finally, the artifact-reducing effects of image rotation [16] and resampling [17], [18] were analyzed.
855
tional (around the image center) misregistration in the ranges of and , respectively. Such a plot was referred to as a registration curve. The ranges of misregistration were specified relative to the registered pose of the image contents. Thus, for example, a translational misregistration of three voxels for the rotated pair meant that the SPECT image was rotated 30 clockwise (to correct for misalignment introduced in creating the image) and then translated three voxels. Both the rotation and the translation were applied together as a single transformation to avoid accumulating errors in two interpolation steps. The transformations were applied to the SPECT image, which contained less signal energy in the high spatial frequencies, and was therefore less susceptible to degradation in image quality from interpolation. MI was expected to attain a maximum when the images were registered. To evaluate MI, the images were overlaid according to the specified pose relative to one another. A joint histogram was created by plotting an ( ) point for every pair of corresponding voxels in the overlap region, where was the voxel intensity in the first image, and was the interpolated voxel intensity in the second image. The voxel intensities were quantized into a fixed number of discrete bins (referred to as intensity bins hereafter). Thus, if the intensities of each image were quantized into intensity bins, the size of the joint histogram would be . The number of points in each histogram bin was normalized by . the total number of points to give a probabilistic measure MI was determined as
II. METHODS A. Human Head Images Registration experiments were conducted on four pairs of images, denoted as “aligned,” “rotated,” “128/129,” and ” (Fig. 1). The aligned pair consisted of a 128 128 axial “ -weighted MR image of a human head and a corresponding -ECD single photon emission computed 128 128 tomography (SPECT) image. The images were obtained using Fourier interpolation from a central axial slice of three-dimensional (3-D) images, which were preregistered using Woods’ method [3] and verified visually. The other three pairs of images were similar to the aligned pair, except that the SPECT image in the rotated pair was rotated 30 anti-clockwise around the center. The SPECT image in the 128/129 pair was rescaled to a grid size of 129 129 by reducing the voxel size by 128/129. pair was rescaled by increasing The SPECT image in the . These four image pairs were created prior the voxel size by to the registration experiments, using Fourier interpolation to resample the SPECT image in order to minimize degradation of the image contents. Even if some degradation was introduced in this step, it should have little influence on the results of the registration experiments, since the degradation did not vary during each experiment. Therefore, any fluctuation of the match metric resulted solely from subsequent operations. B. Evaluating Mutual Information For each pair of images, MI (i.e., the match metric) was determined as a function of translational (along the axis) and rota-
MI was evaluated for different numbers of intensity bins ( , 16, 32, 64, 128, 192, 256, and 512) and with various interpolation schemes. Interpolation artifacts were quantified by measuring the smoothness of each registration curve as follows (Fig. 2). First, the axis of each registration curve was normalized to a range of one. This allowed the interpolation artifacts to be quantified relative to the variation of the match metric. Then, the registration curve was subtracted from a smoothed version of itself, obtained by median filtering followed by convolution with a Hamming filter. Both filters had the same length, corresponding to a range of 1.5 voxels in translation, or 5 in rotation along axis. To avoid edge effects from filtering, the initial the and final points of the curve, equal to the filter length, were discarded. The difference between the registration curve and the smoothed version contained mostly interpolation artifacts (Fig. 2). Thus, smoothness was quantified as one over the root mean square (rms) of the difference curve. A registration curve with fewer interpolation artifacts was expected to be smoother, and would therefore attain a higher smoothness value. C. Interpolators Eight interpolation schemes were tested. For interpolators 1 to 4, the voxel intensity at a nongrid position was interpolated as a weighted sum of voxel intensities at neighboring grid positions. The total weights summed to unity, and the weighting
856
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
Fig. 1. Four image pairs used in registration experiments, denoted as “aligned,” “rotated,” “128/129,” and “=3.” Each pair consisted of an axial T -weighted MR head image and a corresponding Tc-ECD SPECT image. A different transformation (i.e., rotation or scaling) was applied to the SPECT image of each pair using Fourier interpolation. In this figure, the scaling of the SPECT voxel grid in the “128/129” and “=3” pairs is exaggerated for illustration purposes only.
was determined by the interpolation kernel (Fig. 3 and Table I), which was separable along the spatial dimensions. The width of the interpolation kernel determined the number of neighbors involved in the summation. Interpolators with a smaller support relied on fewer neighbors and were therefore faster. If the interpolation kernel extended beyond the image edge, the voxel intensity to be interpolated was considered to be outside the overlap region between the two images and was discarded
from the match metric. This was to avoid any confounding edge effects that might affect the interpolated intensity and alter the image contents. The number of overlapping voxels changed discretely, as the relative pose of the two images varied [14]. 1) Nearest Neighbor: The voxel at the closest grid position was assigned a weight of one, while all other voxels received zero weights.
TSAO: INTERPOLATION ARTIFACTS IN MUTUAL INFORMATION IMAGE REGISTRATION
857
Fig. 2. Quantification of interpolation artifacts. Left: registration curve showing MI as a function of translational misregistration. Center: smoothed version of registration curve. Right: difference between original registration curve and the smoothed version, revealing mainly interpolation artifacts. Smoothness was defined as one over the rms of the difference curve.
Fig. 3. Interpolation kernels in one dimension. The width of the support is shown below.
TABLE I EXPRESSION FOR INTERPOLATOR IN ONE DIMENSION. THE POINT TO INTERPOLATE IS AT POSITION x, BETWEEN 0 AND 1, WITH NEIGHBORING GRID POINTS LISTED IN ORDER AS . . . ; n ; n ; n ; n ; n ; n ; . . . THE POINT TO INTERPOLATE LIES BETWEEN n AND n , AND IS EQUIVALENT TO n WHEN x = 0, AND TO n WHEN x = 1.
2) Linear: The interpolation kernel was triangular in shape along each spatial dimension, with a support width of two voxels.
3) Cubic Catmull–Rom: The interpolation kernel was determined by the Catmull–Rom [21] cubic spline, with a support width of four voxels.
858
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
6
Fig. 4. Translation results for all interpolators and image pairs, with MI plotted as a function of translational misregistration (within a range of 10 voxels along the y axis). Each column of graphs corresponds to a different interpolator, while each row corresponds to a different image pair. The scales of the graph axes are shown in the lower left panel. In each panel, the curves from bottom to top correspond to 8, 16, 32, 64, 128, 192, 256, and 512 intensity bins, respectively.
4) Hamming-Sinc: The interpolation kernel along each spatial dimension was a sinc function multiplied by a six-voxel-wide Hamming window. The kernel was normalized to ensure that the total weights summed to unity [22]. 5) Partial Volume: This interpolator, as introduced by Collignon et al. [7], [23], is not an interpolator in the traditional sense because it does not estimate the intensities of voxels at nongrid positions. Instead, each voxel at a nongrid position is treated as though it consists of several fractional voxels at neighboring grid positions. Since the fractional voxels are at grid positions, their intensities can be evaluated directly without interpolation. The fractions sum to unity and were determined according to proportions from linear interpolation [7], [23]. In general, the fractions may be determined according to proportions from other interpolators as well [24]. 6) Nearest Neighbor With Jittered Sampling: For every voxel to be interpolated, its coordinates were jittered by , adding a normally distributed random offset ( ). The intensity at the jittered coordinates was interpolated by NN interpolation. 7) Nearest Neighbor With Histogram Blurring: Following NN interpolation, the joint histogram was blurred by convolving it with a filter (coefficients 1/16, 4/16, 6/16, 4/16, and 1/16) along each axis of the histogram. 8) Nearest Neighbor With Jittered Sampling and Histogram Blurring: This interpolator is a combination of interpolators 6 and 7. Therefore, jittering was applied to the coordinates of
every voxel to be interpolated. Then, the joint histogram was blurred by convolution, as described earlier. III. RESULTS Figs. 4 and 5 summarize the results for all interpolators (columns) and image pairs (rows). Each graph panel shows MI as a function of translational (Fig. 4) or rotational (Fig. 5) misregistration for different numbers of intensity bins. In each panel, the curves from bottom to top correspond to increasing number of intensity bins. Fig. 6 summarizes the quantification of interpolation artifacts for the translation (column 1) and rotation (column 2) experiments, and for all image pairs (rows). Each graph panel shows the smoothness of the registration curves as a function of the number of intensity bins in a semi natural-log scale for all interpolators. In general, interpolation artifacts were most evident in the translation results of the aligned pair (Fig. 4, row 1). The artifact pattern was characteristic for each interpolator, and it varied with the number of intensity bins. For most interpolators, the artifacts worsened at the high and/or low ends of the intensity bin range. This trend was reflected as a general decrease in smoothness toward the right and/or left sides of each panel in Fig. 6, respectively. In most cases, increasing the number of intensity bins excessively also had the deleterious effect of changing the match metric around 0 from the global maximum to a local minimum
TSAO: INTERPOLATION ARTIFACTS IN MUTUAL INFORMATION IMAGE REGISTRATION
859
6
Fig. 5. Rotation results for all interpolators and image pairs, with MI plotted as a function of rotational misregistration (within a range of 180 around the image center). Each column of graphs corresponds to a different interpolator, while each row corresponds to a different image pair. The scales of the graph axes are shown in the lower left panel. In each panel, the curves from bottom to top correspond to 8, 16, 32, 64, 128, 192, 256, and 512 intensity bins, respectively.
(e.g., See curves from bottom to top in Fig. 5, row 1, column 2). This would lead to improper registration since the match metric no longer attained a maximum when the images were registered. NN interpolation was found to produce stairs-like artifacts in the translation results of the aligned pair (Fig. 4, row 1, column 1). The stairs-like discontinuities occurred at every half-voxel step, coinciding with the discontinuities of the interpolation kernel (Fig. 3). Additionally, there were subtle discontinuities at every integer-voxel step, which occurred when voxels close to the edges of the overlap region either fell into or outside the overlap region, thus changing the voxels contributing to the match metric. The artifacts did not vary appreciably with different number of intensity bins, as reflected by a smoothness curve that was almost flat (Fig. 6, row 1, column 1). The stairs-like artifacts were reduced substantially for the other image pairs (Fig. 4, rows 2 – 4, column 1), although they remained noticeable for the 128/129 pair. This reduction in artifact was reflected quantitatively by an increase in smoothness (Fig. 6, rows 2–4, column 1). These artifacts were eliminated in the rotation results of all image pairs (Fig. 5, column 1). The artifacts of linear, cubic Catmull–Rom and Hamming-sinc interpolation (Figs. 4 and 5, columns 2 – 4, respectively) were similar in nature. In the translation results of the aligned pair (Fig. 4, row 1, columns 2 – 4), the artifacts approached a stairs-like pattern at the low end of the intensity bin range, with the discontinuities occurring at every integer-voxel step. These stairs-like artifacts were accompanied
by high-frequency fluctuations for cubic Catmull–Rom and Hamming-sinc. By increasing the number of intensity bins, the stairs-like artifacts diminished, but a different type of artifacts with an arch pattern began to emerge. The arch artifacts were accentuated with increasing number of intensity bins, and they bifurcated successively into finer arches. These artifacts were mostly eliminated for the other image pairs, except at the highest numbers of intensity bins (Fig. 4, rows 2 – 4, columns 2 – 4). Artifacts were not seen in the rotation results of all image pairs, except for either small dips or spikes at rotation angles that aligned the voxel grids (Fig. 5, columns 2 – 4). Quantitatively, the smoothness of the registration curves did not differ considerably among the three interpolators (linear, cubic Catmull–Rom and Hamming-sinc), and their smoothness curves overlapped in most cases (Fig. 6). Whenever the smoothness value differed among the three interpolators, it tended to be lower for cubic Catmull–Rom and Hamming-sinc interpolation than for linear interpolation, indicating that the interpolators with wider support did not lead to smoother registration curves. PV interpolation produced inverted arches in the translation results of the aligned pair (Fig. 4, row 1, column 5). The inverted arches were accentuated when there were more intensity bins, as reflected by a drop in smoothness (Fig. 6, row 1, column 1). These inverted arches changed to arches for the 128/129 pair (Fig. 4, row 3, column 5), but were not seen for the rotated and pairs (Fig. 4, rows 2 and 4, column 5), as reflected by the the increased smoothness. In the rotation results (Fig. 5, column
860
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
Fig. 6. Smoothness of the registration curves as a function of the number of intensity bins in a semi natural-log scale. Columns 1 and 2 correspond to the translation and rotation experiments, respectively. Each row corresponds to a different image pair. In each panel, the different curves correspond to the smoothness of the registration curves for different interpolators. Note that the scale on the Y axis is different for each panel.
5), the only artifacts seen for all image pairs consisted of either dips or spikes at rotation angles that aligned the voxel grids. For interpolators that included JIT (i.e., interpolators 6 and 8) (Figs. 4 and 5, columns 6 and 8), the artifacts consisted of small stochastic fluctuations only. Since the jittering degraded the similarity between the images, it decreased the value of the match metric slightly, compared with interpolators without JIT. However, this decrease should have no effect on registration accuracy since it did not depend on the alignment of the image contents. In general, the smoothness values of these interpolators varied over a narrower range than other interpolators (Fig. 6),
due to the absence of any severe interpolation artifacts. At times, these interpolators considerably outperformed the others, such as in the translation results of the aligned pair (Fig. 6, row 1, column 1). However, in situations where the artifacts of other interpolators were suppressed (e.g., in the rotation results of pair, Fig. 6, row 4, column 2), the slight stochastic the fluctuations resulted in a slightly lower overall smoothness. For interpolators that included BLUR (i.e., interpolators 7 and 8) (Figs. 4 and 5, columns 7 and 8), the artifacts were similar to those without BLUR. Nevertheless, BLUR generally improved smoothness (Fig. 6). The blurring also decreased the overall
TSAO: INTERPOLATION ARTIFACTS IN MUTUAL INFORMATION IMAGE REGISTRATION
value of the match metric. The main impact of BLUR was to maintain the maximum of the match metric around 0 when there was a large number of intensity bins (Fig. 5). In other words, the registration curve maintained its shape over widely different number of intensity bins. This was reflected as the increased flatness of the smoothness curves (Fig. 6). IV. DISCUSSION Interpolation artifacts can occur if the voxel grids of the images have coinciding periodicities. Since these artifacts may cause severe fluctuations in the match metric, they can limit the achievable registration accuracy and undermine the robustness of the match metric. The present results showed that the artifacts were reduced substantially when the periodicities of the voxel grids were abolished with strategies such as JIT, or when one image was rotated or resampled relative to the other image. A. Number of Intensity Bins In the present study, the artifact patterns were characterized over widely different number of intensity bins. Overall, the artifact patterns were found to change smoothly with the number of intensity bins (Figs. 4 and 5). Although extreme choices should be avoided, the precise choice for the number of intensity bins may not be critical. In general, it is preferable to err on the low side, since having fewer intensity bins allows for a faster processing time and a smaller size for the joint histogram. Otherwise, if the joint histogram is large, each histogram bin will be populated by a small number of points, and the match metric may suffer from poor statistics. This is similar to the adverse effect seen with image subsampling [17], [18], since both procedures lead to fewer points occupying each histogram bin [25]. If the number of intensity bins is chosen to be even higher, most of the points in the joint histogram may fail to aggregate sufficiently to form clusters. Thus, the only prominent clusters remaining will correspond to the background voxels, due to the sheer number of these voxels [26]. As a result, the match metric will be increasingly influenced by the nonspecific matching of background voxels, and registration accuracy will degrade [26]. This problem was encountered in the present study for most interpolators. When there were too many intensity bins, the match metric at 0 changed from the global maximum to a local minimum, thus leading to improper registration (Fig. 5). This problem may be of particular concern in multiresolution [14], [27] or multiscale [25] matching, since the number of voxels may be small at certain levels of matching, so the number of intensity bins may need to be adjusted accordingly [14], [27]. B. Classification of Interpolators The interpolators examined in this study can be divided into two groups: intensity interpolators (linear, cubic Catmull–Rom, and Hamming-sinc), and nonintensity interpolators (the remaining ones). The main difference between the two groups is that intensity interpolators may introduce new intensity values that are not present in the original image [7], [23]. For example, a binary-valued image may contain fractional intensity values after interpolation. In contrast, nonintensity interpolators only return intensity values from the original image. Therefore, the
861
mechanism of artifact generation is fundamentally different between the two groups. Intensity interpolators generate artifacts by modifying voxel intensities, while nonintensity interpolators generate artifacts by changing the way the match metric is evaluated. As a result, the artifacts from intensity interpolators are inherently more dependent on the image contents, thus making them somewhat difficult to analyze. C. Linear, Cubic Catmull–Rom, & Hamming-Sinc Interpolation Usually, the performance of an intensity interpolator is characterized by its frequency response (or power spectrum) [20]. While this description is suited to determining the amount of interpolation error, it is inadequate in describing the impact of interpolation on the match metric. This is because the match metric is a highly nonlinear function of the voxel intensities. This inadequacy was shown in the present paper by the somewhat unexpected finding that the switch from linear to cubic Catmull–Rom to Hamming-sinc interpolation led to no improvement in the smoothness of the registration curves. Instead, it even worsened the smoothness some of the times. These results suggest that registration may not benefit from using intensity interpolators with wider support, despite their improved frequency response. This finding is important from a practical standpoint, since interpolators with wider support are also computationally more demanding. The present results also suggest that intensity interpolators may not be preferable to nonintensity interpolators, since the former did not appear to be any less susceptible to producing artifacts than the latter. The intensity interpolators being studied (i.e., linear, cubic Catmull–Rom, Hamming-sinc) were found to produce stairs-like artifacts with a small number of intensity bins and arch artifacts with more intensity bins. The stairs-like artifacts were expected since these interpolators would behave in a similar manner to NN interpolation, if the intensities were discretized to very few levels. The stairs-like discontinuities occurred at every integer-voxel step when the interpolation kernel switched to a different set of neighboring voxels for interpolation. For the arch artifacts at higher numbers of intensity bins, their exact shapes vary with image contents. For example, by adding 10% noise pixel-wise (without filtering) to the SPECT image, the artifacts become more pronounced, adopting a rounded appearance (Fig. 7). This agrees with the findings of Pluim et al. [17], [18], who attributed the origin of these artifacts in linear interpolation to the intensity averaging effect of the interpolator. The present study showed that further increases in the number of intensity bins led to bifurcations of these arch artifacts. The onset of arch bifurcation was delayed to a higher number of intensity bins for interpolators with a larger support (i.e., cubic Catmull–Rom, Hamming-sinc). While intensity averaging is a likely cause of these arch artifacts, their complicated shapes (Fig. 4, columns 2–4) suggest that there may be additional contributing mechanisms, which deserve further investigation in the future. D. Partial Volume Interpolation The present study showed that PV interpolation could lead to inverted-arch artifacts (Fig. 4, row 1, column 5) as well as
862
Fig. 7. Artifact patterns of linear, cubic Catmull–Rom, and Hamming-sinc interpolation, after adding 10% noise pixel-wise to the SPECT image of the aligned pair. The noise level indicates the standard deviation of noise relative to the maximum intensity. The noisy SPECT image is shown on the upper left. The registration curves show MI as a function of translational misregistration, for the aligned pair with 64 intensity bins.
arch artifacts (Fig. 4, row 3, column 5). Both types of artifacts originated from the splitting of misaligned voxels into fractional parts, since the fractional parts contributed less in total to the match metric than the whole voxel [8], [28]. Fig. 8 illustrates the origin of these artifacts. Column 1 magnifies one of the registration curves from column 5 of Fig. 4. The grid alignment is illustrated for three particular points of the registration curves ( , , and ) in columns 2–4, respectively. The plots show the distance between each voxel row in the SPECT image and the nearest voxel row in the MR image. This distance is referred to as fractional grid misalignment. For example, a fractional grid misalignment of zero means that the voxel rows in both images are exactly aligned, while a fractional grid mismeans that the SPECT voxel row is right alignment of in-between two nearest voxel rows in the MR image. The inverted-arch artifacts were manifested in translation when the ratio between the voxel sizes of the two images was either one (i.e., both images have identical voxel sizes) or a rational fraction of small numbers (e.g., 3/5) [15]–[18]. The bottom of the inverted arches corresponded to the majority of the voxels being misaligned (Fig. 8, row 1, column 4). In contrast, the arch artifacts were manifested when there were only slight mismatches between the voxel sizes, as found in the 128/129 pair (Fig. 4, row 3, column 5). Due to the small voxel size difference, the amount of grid misalignment differed slightly for each voxel row (Fig. 8, row 2). Thus, some portions of the image had smaller fractional grid misalignment (i.e., were more aligned) than others. The more aligned portions would have higher contribution to the
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
match metric. For example, when the 128/129 pair was registered (Fig. 8, row 2, column 2), voxel rows close to the edges of the image (row 0 or 128; see black arrows) were more aligned than those close to the center (e.g., row 64; see open arrows). Since the voxel rows at the edges corresponded primarily to the image background, their higher contribution to the match metric would lead to an overall decrease in the match metric, resulting in the dips between the arches. The top of the arches occurred at every half-voxel step, when the voxel rows in the image center came into grid alignment, thus contributing more to the match metric (See open arrows in Fig. 8, row 2, column 4). Accordingly, if the image contents were shifted by half a field of view, the top of the arches would occur at every integer-voxel step instead (data not shown). These arch and inverted-arch artifacts can be eliminated by choosing the voxel sizes judiciously. For example, a salient choice was for the voxel sizes of one image to be times the voxel sizes of the other image, since the resultant amount of grid misalignment was considerably different from row to row (Fig. 8, row 3). Thus, all portions of the voxel grids were misaligned by the same amount on average, regardless of the translation between the images (Fig. 8, row 3, columns 2 – 4). With being an irrational number, the amount of misalignment also never repeated itself exactly from row to row, thereby reducing any periodicity. The artifacts of PV interpolation were found to reduce dramatically in the translation results when the relative voxel size was changed from 128/129 (Fig. 4, row 3, (Fig. 4, row 4, column 2). column 2) to E. Jittered Sampling Considerable artifact reduction was achieved in all cases with JIT, which abolished the periodicities of the voxel grids and eliminated the associated artifacts (see Likar and Pernuˇs [25] for a related approach). Even though JIT in the present implementation deliberately introduced error in the interpolation process, it proved to be an effective interpolator for registration, since it introduced approximately the same amount of error regardless of the extent of grid alignment. F. Histogram Blurring Improved registration robustness was found with BLUR, which maintained the global maximum of the match metric at 0 even when there was a large number of intensity bins. By applying BLUR, the points in each histogram bin were spread into neighboring bins, which allowed the amount of clustering in the joint histogram to be estimated accurately, even when there were very few points per histogram bin. In fact, BLUR can be considered an approximation of Parzen windowing [29] and gridding [30], which are two mathematically similar techniques of estimation from discrete samples. In the present context, the blurring kernel served as the window function of Parzen windowing, or the convolution function of gridding. The approximation became more accurate with more intensity bins, due to smaller discretization error, albeit at the expense of increased computation. In practice, for a moderate number
TSAO: INTERPOLATION ARTIFACTS IN MUTUAL INFORMATION IMAGE REGISTRATION
863
Fig. 8. Artifact pattern of PV interpolation depends on the ratio between the voxel sizes of the images. Rows 1 – 3 show the results of the aligned, 128/129, and =3 pairs, with voxel ratios of 1/1, 128/129 and =3, respectively. Column 1 shows a magnified version of the respective registration curves, corresponding to column 5 of Fig. 4. The artifact patterns arose from different grid alignments, as the two images were translated from one another. The alignments at three particular points ( , , and ) are illustrated in columns 2 – 4, respectively. The plots show the distance between each voxel row in the SPECT image and the nearest voxel row in the MR image. This distance is referred to as fractional grid misalignment. Arrows point to relevant features. See text for details.
of intensity bins, BLUR provides an efficient approximation of Parzen windowing to improve registration accuracy. G. Related Match Metrics and Generalization Although this study focused on MI only, the observations can be generalized to other methods such as moment-based [4], and joint entropy [5] methods, since the match metrics of these methods can be expressed in a similar mathematical form as MI [8], [9], [28]. The present observations also apply to normalized .A MI [31], which can be expressed as (MI)/(joint entropy) natural question is whether the present results also apply to other image pairs. Results from previous studies [8], [16], [26], [32] and from other groups [15], [17], [18], [25] suggest that while the exact shapes of the artifacts are dependent on the image contents (especially for intensity interpolators), the general artifact patterns are similar over a wide range of images, including images of different modalities and size scales, and in two and three dimensions (Fig. 9). This is perhaps not too surprising since the artifacts are primarily related to implementation issues of the match metric, which are independent of the image contents. H. Summary In conclusion, registration using MI or similar match metrics can be affected severely by interpolation artifacts. Nevertheless,
Fig. 9. Translation results for an aligned pair of 3-D SPECT and MRI images with 128 intensity bins and identical voxel sizes. MI is plotted as a function of translational misregistration (within a range of 10 voxels along the y axis). Results for the different interpolators are staggered along the vertical axis for illustration purposes.
6
the artifacts can be reduced significantly if the voxel grids are misaligned. In most registration tasks, especially in three dimensions, the voxel grids are unlikely to be in exact alignment, so
864
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
interpolation artifacts are unlikely to be a major source of registration error [15]. In those cases, it may suffice to use NN interpolation, thus allowing a considerable reduction in processing time. Even in those cases, it may be worthwhile to incorporate BLUR, which is computationally inexpensive for a moderate number of histogram bins, and was found to always improve the smoothness of the registration curves. In cases where the voxel grids are almost in alignment (e.g., when a head holder is used to minimize differences in head positioning, such as in fMRI acquisitions [32], [33]), special attention must be paid to reduce interpolation artifacts in order to ensure registration accuracy [15]. This work showed that the following strategies could reduce artifacts considerably and improve registration robustness: 1) avoiding an extreme number of intensity bins; 2) resampling the images in a rotated orientation with unequal voxel sizes (e.g., voxel sizes of one image being times those of the other image); 3) JIT; and 4) BLUR. The first two strategies constitute precomputation steps, and do not increase the amount of processing at each iteration of the registration algorithm. Applying these two strategies together often suffices to eliminate most interpolation artifacts, so NN interpolation combined with BLUR can be used for fast performance, while problems from interpolation artifacts and/or an inappropriately chosen number of intensity bins can be avoided. The rare exception is in pathological situations during the registration process that happen to realign the voxel grids. To ensure small interpolation artifacts at all time, JIT can be incorporated with minimal increase in processing time.
ACKNOWLEDGMENT The author would like to thank Dr. P. C. Lauterbur and Dr. S. R. Teixeira for helpful comments on the manuscript, and Dr. M. Ichise for the SPECT image.
REFERENCES [1] M. W. Vannier and D. E. Gayou, “Automated registration of multimodality images,” Radiology, vol. 169, pp. 860–861, 1988. [2] D. L. Hill, P. G. Batchelor, M. Holden, and D. J. Hawkes, “Medical image registration,” Phys. Med. Biol., vol. 46, pp. R1–45, 2001. [3] R. P. Woods, J. C. Mazziotta, and S. R. Cherry, “MRI-PET registration with automated algorithm,” J. Comput. Assist. Tomogr., vol. 17, pp. 536–546, 1993. [4] D. L. G. Hill, C. Studholme, and D. J. Hawkes, “Voxel similarity measures for automated image registration,” Proc. Vis. Biomed. Comput., pp. 205–216, 1994. [5] A. Collignon, D. Vandermeulen, P. Suetens, and G. Marchal, “3 D multimodality medical image registration using feature space clustering,” in Proc. Int. Conf. Computer Vision Virtual Reality Robot Medicine, 1995, pp. 195–204. [6] W. M. Wells 3rd, P. Viola, H. Atsumi, S. Nakajima, and R. Kikinis, “Multi-modal volume registration by maximization of mutual information,” Med. Image Anal., vol. 1, pp. 35–51, 1996. [7] F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, and P. Suetens, “Multimodality image registration by maximization of mutual information,” IEEE Trans Med Imag., vol. 16, pp. 187–198, Apr. 1997. [8] J. Tsao and P. C. Lauterbur, “Generalized clustering-based image registration for multi-modality images,” Proc. IEEE Engineering Medicine Bioloy Soc, pp. 667–670, 1998.
[9] M. Bro-Nielsen, “Rigid registration of CT, MR and cryosection images using a GLCM framework,” in Proc CVRMed/MRCAS, 1997, pp. 171–180. [10] J. V. Hajnal, N. Saeed, E. J. Soar, A. Oatridge, I. R. Young, and G. M. Bydder, “A registration and interpolation procedure for subvoxel matching of serially acquired MR images,” J. Comput. Assist. Tomogr., vol. 19, pp. 289–296, 1995. [11] J. L. Ostuni, A. K. Santha, V. S. Mattay, D. R. Weinberger, R. L. Levin, and J. A. Frank, “Analysis of interpolation effects in the reslicing of functional MR images,” J. Comput. Assist. Tomogr., vol. 21, pp. 803–810, 1997. [12] T. M. Lehmann, C. Gonner, and K. Spitzer, “Survey: Interpolation methods in medical image processing,” IEEE Trans. Med. Imag., vol. 18, pp. 1049–1075, Nov. 1999. [13] W. F. Eddy, M. Fitzgerald, and D. C. Noll, “Improved image registration by using Fourier interpolation,” Magn. Reson. Med., vol. 36, pp. 923–931, 1996. [14] P. Thévenaz and M. Unser, “Optimization of mutual information for multiresolution image registration,” IEEE Trans. Image Processing, vol. 9, pp. 2083–2099, Dec. 2000. [15] F. Maes, “Segmentation and Registration of Multimodal Medical Images: From Theory, Implementation and Validation to a Useful Tool in Clinical Practice,” Catholic University of Leuven, Leuven , Belgium, 1998. [16] J. Tsao, “Efficient interpolation for clustering-based multimodality image registration,” in Proc. Int. Soc. Magnetic Resonence Medicine, 1999, p. 2195. [17] J. P. W. Pluim, J. B. A. Maintz, and M. A. Viergever, “Mutual information matching and interpolation artifacts,” Proc SPIE, pp. 56–65, 1999. , “Interpolation artefacts in mutual information-based image regis[18] tration,” Comput. Vis. Image Understanding, vol. 77, pp. 211–232, 2000. [19] J. West et al., “Comparison and evaluation of retrospective intermodality brain image registration techniques,” J. Comput. Assist. Tomogr., vol. 21, pp. 554–566, 1997. [20] P. Thévenaz, T. Blu, and M. Unser, “Interpolation revisited,” IEEE Trans. Med. Imag., vol. 19, pp. 739–758, 2000. [21] E. Catmull and R. Rom, “A class of local interpolationg splines,” in Computer Aided Geometric Design, R. E. Barnhill and R. F. Riesenfled, Eds. New York: Academic Press, 1974. [22] N. A. Thacker, A. Jackson, D. Moriarty, and E. Vokurka, “Improved quality of re-sliced MR images using re-normalized sinc interpolation,” J. Magn. Reson. Imaging, vol. 10, pp. 582–588, 1999. [23] A. Collignon, F. Maes, D. Delaere, D. Vandermeulen, P. Suetens, and G. Marchal, “Automated multimodality medical image registration using information theory,” in Proc. Int. Conf. Information Processing Medical Imaging: Computational Imaging and Vision, 1995, pp. 263–274. [24] H. M. Chen and P. K. Varshney, “Generalized partial volume interpolation for image registration based on mutual information,” presented at the Proc. IEEE Signal Processing Soc., Western New York, Image Processing Workshop, Univ. Rochester, Ctr. Optoelectron. Imag., 2000. [25] B. Likar and F. Pernuˇs, “A hierarchical approach to elastic registration based on mutual information,” Image Vis. Computing, vol. 19, pp. 33–44, 2001. [26] J. Tsao, “Multimodality image registration by gradient enhanced generalized clustering,” in Proc. Int. Soc. Magn. Reson. Med., 1999, p. 125. [27] M. Jenkinson and S. Smith, “A global optimization method for robust affine registration of brain images,” Med. Image Anal., vol. 5, pp. 143–156, 2001. [28] T. M. Buzug and J. Weese, “Image registration for DSA quality enhancement,” Comput. Med. Imag. Graph., vol. 22, pp. 103–113, 1998. [29] E. Parzen, “On estimation of a probability density function and the mode,” Ann. Math. Statist., vol. 33, pp. 1065–1076, 1962. [30] J. D. O’Sullivan, “A fast sinc function gridding algorithm for Fourier inversion in computer tomography,” IEEE Trans. Med. Imag., vol. MI-4, pp. 200–207, 1985. [31] C. Studholme, D. L. G. Hill, and D. J. Hawkes, “An overlap invariant entropy measure of 3 D medical image alignment,” Pattern Recogn., vol. 32, pp. 71–86, 1999. [32] J. Tsao, J. B. M. Goense, and P. C. Lauterbur, “fMRI image registration using total local density of intensity mapping plot,” in Proc. Int. Soc. Magnetic Resonance Medicine, 1998, p. 2135. [33] B. Kim, J. L. Boes, P. H. Bland, T. L. Chenevert, and C. R. Meyer, “Motion correction in fMRI via registration of individual slices into an anatomical volume,” Magn. Reson. Med., vol. 41, pp. 964–972, 1999.
IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 22, NO. 7, JULY 2003
1183
Correspondence________________________________________________________________________ Correction to “Interpolation Artifacts in Multimodality Image Registration Based on Maximization of Mutual Information”
REFERENCES [1] J. Tsao, “Interpolation artifacts in multimodality image registration based on maximization of mutual information,” IEEE Trans. Med. Imag, vol. 22, pp. 854–864, July 2003.
Jeffrey Tsao In the above paper [1], some of the information in Fig 5 was missing. The figure should have appeared as shown here. Manuscript received July 30, 2003. The author is with the Institute for Biomedical Engineering, Swiss Federal Institute of Technology, Zurich, Building ETF, Room C108, Sternwartstrasse 7, 8092 Zurich, Switzerland (e-mail [email protected]). Digital Object Identifier 10.1109/TMI.2003.817660
6
Fig. 5. Rotation results for all interpolators and image pairs, with MI plotted as a function of rotational misregistration (within a range of 180 around the image center). Each column of graphs corresponds to a different interpolator, while each row corresponds to a different image pair. The scales of the graph axes are shown in the lower left panel. In each panel, the curves from bottom to top correspond to 8, 16, 32, 64, 128, 192, 256, and 512 intensity bins, respectively. 0278-0062/03$17.00 © 2003 IEEE
