|Home | About | Journals | Submit | Contact Us | Français|
We show that image registration using conventional interpolation and summation approximations of continuous integrals can generally fail because of resampling artifacts. These artifacts negatively affect the accuracy of registration by producing local optima, altering the gradient, shifting the global optimum, and making rigid registration asymmetric. In this paper, after an extensive literature review, we demonstrate the causes of the artifacts by comparing inclusion and avoidance of resampling analytically. We show the sum-of-squared-differences cost function formulated as an integral to be more accurate compared with its traditional sum form in a simple case of image registration. We then discuss aliasing that occurs in rotation, which is due to the fact that an image represented in the Cartesian grid is sampled with different rates in different directions, and propose the use of oscillatory isotropic interpolation kernels, which allow better recovery of true global optima by overcoming this type of aliasing. Through our experiments on brain, fingerprint, and white noise images, we illustrate the superior performance of the integral registration cost function in both the Cartesian and spherical coordinates, and also validate the introduced radial interpolation kernel by demonstrating the improvement in registration.
Image resampling is indispensable in many image processing applications. Given that digital images are represented as discrete values on a regular grid, virtually any image transformation, such as translation, rotation, zooming, or nonlinear warping, requires mapping of the grid points (pixels) to a new set of coordinate points that do not generally correspond to the original grid. To present the transformed image on the same regular grid, it is thus necessary to estimate, i.e. interpolate, the image values between the grid points where no sampled data are available.
Accurately speaking, interpolation is fitting a continuous function passing through the image pixels or voxels (see – for reviews). Interpolation is required in advance of resampling this continuous function on the desired points . In practice, these two processes are commonly combined, sometimes simply referred to as resampling. It has been known that both interpolation and resampling can produce artifacts in the resulting image (see Sec. II for a literature review). Interpolation with non-ideal kernels not only blurs the image, but also creates artificial high-frequency components. In addition, aliasing distortions can appear during resampling of such a non-ideally interpolated image (and even an ideally interpolated one, if rotation is involved).
In addition to helping to represent transformed images on a regular grid, interpolation and resampling are used in other procedures, such as image alignment or registration, as discussed in this paper. Registration provides a transformation maximizing a similarity measure between the transformed version of an image and a second reference image (see – for reviews). Interpolation and resampling are commonly used in the transformation step of the sub-pixel registration, and in turn produce the aforementioned artifacts–notably patterns of local extrema in the similarity measure–making the optimization challenging. The artifacts may in some cases cause the registration to fail completely, for example by lowering the signal-to-noise ratio of the image gradient leading the algorithm to wrong optima, by shifting the global optimum significantly, or by making a false stronger global optimum appear in another region of the image. Additionally, resampling generally breaks the inherent symmetry (inverse-consistency ) of rigid image registration , since the artifacts alter only one of the images. The fact that the transformed image is merely used to compute the similarity measure raises the following question: Are interpolation and resampling both necessary in registration, or would interpolation alone be sufficient? As we will see in Sec. III, the latter is true, in which case aliasing artifacts due to resampling can be avoided by using an integral in the registration cost function instead of a sum, without resorting to costlier closer-to-ideal interpolation kernels. Another question addressed in this paper is whether using ideal interpolants can eliminate resampling artifacts and make the two cases of including and avoiding the resampling stage equivalent.
In this work, we first review the literature extensively and provide a comprehensive list of references where resampling artifacts have been observed in image registration, and also point out several studies in which an integral cost function has been used to avoid resampling. Although these artifacts have been explored in a number of papers and there are instances in the literature where they are mostly prevented by avoiding resampling, to the best of our knowledge, the two cases of keeping and avoiding the resampling step have not been analytically compared and contrasted. We (a) juxtapose inclusion and avoidance of resampling by showing the resulting cost functions side by side in a simple case of registration, and compare their corresponding errors analytically. In particular, we show that the resampling errors are asymmetric, resulting in the inverse-inconsistency of registration, and that by avoiding resampling, the rigid image registration is kept symmetric. Next, we recall that contrary to the popular belief that employing the ideal (sinc) interpolator would eliminate the artifacts in rigid registration, aliasing still occurs when rotation is involved, and show that using multiple ideal low-pass filters to avert it biases the registration. We then (b) propose to use alternative interpolation kernels corresponding to disk/ball-shaped anti-aliasing filters that preserve as much signal energy as possible without biasing the registration, while preventing aliasing. With this approach such filters are implicitly incorporated in the already existing step of interpolation in image registration. We show experimental results on several images, including examples of triangular-mesh image representation on the sphere, while taking two different analytical and stochastic implementation approaches to compute the integral similarity measure. We demonstrate the improvement in registration accuracy achieved by both eschewing resampling and employing the new interpolation kernel, and in registration symmetry in the former case.
We continue by reviewing the relevant prior work in Sec. II. Section III compares the two cases of including and avoiding the resampling step, and Sec. IV demonstrates rotation-induced aliasing while introducing an unbiased kernel to prevent it. Practical considerations are mentioned in Sec. V, and experimental results are provided in Sec. VI. Section VII concludes the paper with a few final remarks.
We review four categories of the related prior work here, presenting them mostly in a chronological order in each subsection.
Artifacts arising from signal interpolation and resampling have been initially studied outside the context of registration , , , , and further exploited in exposing digital forgeries . Ref.  clearly distinguishes resampling from interpolation, and is the first paper that we are aware of to mention that resampling aliases the higher frequencies of the non-ideal interpolant into the lower frequencies (see Sec. III), a concept which was subsequently restated by the authors of , . We will see in Sec. III why this is the major cause of the artifacts in image registration.
Sub-pixel image registration artifacts appear mostly as scalloping patterns with the period of one pixel in the rigid registration objective function. One of the first attempts to alleviate this issue was to remove the periodic elements with such a frequency, and was applied in realignment of functional magnetic resonance imaging (fMRI) time series . These artifacts were primarily noticed in the entropy-based mutual information (MI, –) methods. In , the authors studied the MI metric and noticed sudden changes in the metric for grid-aligning transformations. It was shown that the two commonly used linear and partial-volume  interpolation methods produced dissimilar patterns of artifact, and that a slight resampling of one of the images might make these patterns smoother. Furthermore, they combined gradient information in the match metric to compensate for the artifacts . Random resampling and inclusion of the prior joint probability were shown in  to reduce such artifacts. The authors of  presented a generalized partial volume estimation that by properly choosing the kernel function, reduces the artifacts in MI-based registration. Various interpolators for MI were compared in  and several strategies to reduce the artifacts were proposed. Ref.  presented an estimation error cancel method which reduces the peak estimation error when the sum of squared differences (SSD) or the cross-correlation (CC) is computed discretely and then interpolated. The extent and the nature of artifactual displacements produced by non-rigid SSD-based registration techniques for different interpolators were compared in . In , it was shown that quasi-random sampling based on Halton sequences alleviates the grid effect in MI registration. The authors of  observed reduced interpolation artifacts when using a full-image mutual histogram as compared to only a sub-volume one. In , the noise reduction filtering that occurs when image samples are interpolated was hypothesized to be the reason behind the artifacts, and a constant-variance filter for linear interpolation was presented to remove them in CC- and MI-based registration. The authors of  proposed the use of Euclidean distance in partial volume interpolation to improve the artifacts in MI-based registration. Refs , ,  further report observing these artifacts in registration.
Virtually all of the above references propose alterations in the registration algorithm that lessen the artifacts, or at least make them less noticeable by reducing their regularity. Next, we mention a rather different category in the literature, analyzing the artifacts from the resampling view point.
Formulating the registration cost function using an integral over the continuous interpolated images–instead of a sum over discrete resampled values–eliminates the resampling step (see Sec. III). The integral SSD has been employed outside the registration framework to assess interpolation, for instance when there is a jump in the original signal at a non-grid location , or in parameter optimization for cubic convolution interpolation . Conversely, the resampling error when using the discrete SSD has been shown to depend on the lattice location .
As for sub-pixel registration, omitting the resampling step was first suggested yet again in MI maximization , where computing the joint histogram from continuous interpolated images was argued to be, in principal, better than doing so from resampled images, because it is free of sampling effects. The authors also showed that the partial volume interpolator  is a special case of the continuous evaluation of the MI from interpolated images, when a boxcar kernel is used. The continuous approximation of MI was performed by computing the MI discretely from over-sampled images. Alternative implementations of this method were later provided via nonparametric windows , autocorrelated kernels , and numerical approximation of the entropy , .
Similarly to , it was demonstrated in  that resampling the image using finite-length interpolants, such as linear and cubic convolution, acts as a low-pass filter which affects the registration differently depending on the grid location. In , the authors defined the CC similarity measure as the integral of the product of the interpolated images and computed it analytically in translation-only registration, thereby circumventing resampling. Later, they pointed out in  that the energy of the resampled image appearing in registration objective functions (SSD, CC, or MI) can have local optima, and quantified the oscillation artifacts similarly to . They showed that these artifacts disappear if the resampling step is skipped, which they proposed to do in rigid registration by computing the objective function using stochastic integration (instead of summation) .
Thorough mathematical analysis of the resampling artifacts in MI-based registration would be a complex task, especially because the computation of the image histogram is involved. Concerning the SSD- and CC-based registration algorithms, noteworthy investigation has been conducted in the above references. Nonetheless, we have not encountered a clear analytical comparison between the errors introduced in the two cases of performing resampling (using the sum cost function) and leaving it out (using the integral cost function), which we aim to provide in Sec. III. We also discuss the implementation of the more accurate resampling-free registration in Sections V and VI-A.
As we will explain in Sec. IV, aliasing occurs in resampling of a rotated image, even when the ideal interpolant is employed. To our knowledge, this was originally demonstrated in , and afterwards mentioned again in –. Low-pass filtering the image has been suggested to eliminate the resampling artifacts ,  (outside the registration context), and similar ideas have been employed in the steerable pyramids framework – as well. Multi-pass rotation algorithms, such as those decomposing the rotation into sequences of 1D translations through shearing , –, have been shown to introduce even further aliasing in the image – (although negligible ).
Regarding the image registration application, aliasing in the image (as opposed to the common frequency) domain due to rotation was discussed in  for phase-correlation methods. In a similar context, the authors of  addressed the issue of ringing artifacts in rigid registration of fMRI datasets by post-registration filtering. In Sec. IV, we will propose the use of specific interpolation kernels in image registration to avoid aliasing in the frequency domain.
In this section, we consider the sum of squared differences (SSD) as the registration metric to analyze the effect of resampling in the simple case of translation-only registration. For the mutual information (MI) metric, we refer the interested reader to –.
Let f, g: d → be two continuous d-dimensional images, the axes of which can be aligned via translating one of them, f, to minimize the following SSD cost function:
The subscript id in the optimal shift id stands for ideal, meaning that we consider the registration results of the original unsampled images as the ground truth. This is justified by the fact that acquiring higher-resolution images is expected to increase the precision of the registration, and at the limit, the original continuous version of the image (with infinite resolution) produces the ideal solution that can be considered as the ground truth while comparing the two approximate (integral and sum) cost functions.
By expanding the integral in Eq. (1) and taking into account that the integral of the shifted f2(·) over the entire space is independent of the shift, the equation can be simplified to:
which is equivalent to maximizing the cross correlation (CC) of the two images. Using Parseval’s theorem and the shift property of the Fourier transform, f(x − Δ) ↔ (v)e−i2πv · Δ, this can also be written as:
with * being the complex conjugate, and and ĝ the Fourier transforms of f and g, which for the reasons that we will see below are assumed to be nonzero only in [−1/2, 1/2]d.
Now, suppose that the two images have been sampled on discrete points with the rate of one sample per unit length in each spatial direction. We try to approximate the same integral cost function as (1) by interpolating the two sampled images, i.e. convolving them with the kernel h. Note that the polynomial (e.g. cubic) spline interpolation is also applicable here, since it is equivalent to a convolution with the corresponding cardinal spline function .
Sampling a continuous function s(x) (1D for simplicity), in mathematical terms, is multiplying s by the impulse train (Dirac comb), Σn δ(x − n). This multiplication can be viewed in the frequency domain as the convolution of ŝ Fig. 1a) with the Fourier transform of the impulse train, also an impulse train, resulting in replicas of ŝ placed at the frequencies n , as in Σnŝ(ν − n) (Fig. 1b). According to the sampling theorem , if s is bandlimited, i.e., ŝ(ν) = 0 for |ν| > 1/2 (as in Fig. 1a), it can be reconstructed exactly from its sampled points using the ideal low-pass filter (Fig. 1b, red dashed line), i.e. by interpolation with the sinc function. Otherwise, the aliasing effects make it impossible to recover the original signal , negatively impacting the accuracy of the alignment , . For simplicity, here we assume f and g to be bandlimited, meaning that they only have frequency components in [−1/2, 1/2]d.
The infinite impulse response (IIR) of the sinc function makes the ideal low-pass filter difficult to implement (although some IIR filters can be efficiently evaluated using, e.g., recursive filtering), and the use of finite impulse response (FIR) interpolants (e.g. linear interpolation; the green dotted curve in Fig. 1b) more attractive. Nevertheless, in contrast to their low computational cost, such non-ideal kernels have two undesirable properties: attenuation of low-frequency components which are supposed to remain intact (ν1 in Fig. 1c), and leaving behind some high-frequency components which are supposed to be removed (ν2 in Fig. 1c). The latter is the major cause of the resampling artifacts in general, as explained here and in , , .
Let h be the kernel that we use to interpolate the two sampled images. Convolving the sampled versions of images f and g with h results in continuous images that we use in place of f and g in the cost function of Eq. (1). As in the previous case, the integrals of such continuous functions in the entire space are independent of the shift , and therefore SSD minimization reduces again to CC maximization (similarly to Eqs. (2,3)). Convolution of the sampled version of f with h can be seen in the frequency domain as multiplication of the replicas of by the frequency response of the kernel, ĥ, resulting in ĥ(v) Σnd (v − n). The CC in this case becomes the following integral, which we simplify by dividing it into finite-length segments (see also ):
where we used the bandlimitedness of f and g to reach line 4 from line 3 of this equation. Thus, the optimum shift using the integral cost function, l, becomes:
By comparing Eqs. (3) and (5), it becomes apparent that the error factor v (Δ): = Σkd |ĥ(v + k)|2 ei2πk·Δ for v [−1/2, 1/2]d, which is periodic with the period of 1 pixel in every dimension, is introduced as a result of non-ideal interpolation. v(Δ) can be seen as the discrete-time Fourier transform (DTFT) of hv[k]: = |ĥ(v − k)|2. The error-free case corresponds to , meaning that , or:
which is the ideal low-pass filter. Note that there is no restriction on the phase of ĥef (v), since the filter is applied to both images and the phase is cancelled in the CC. The intrinsic symmetry of rigid registration also becomes apparent when images are treated as continuous, i.e. in Eqs. (3) and (5), meaning that swapping f and g and negating Δ does not change the cost function.
It is common in image registration to compute the SSD merely on the grid points, while interpolating and resampling only one of the images, f, as needed:
where fΔ[·] is the image f after being sampled, interpolated, shifted with Δ, and resampled on the original grid, and g[·] is the image g sampled on the reference grid. The optimum shift of the above sum cost function, s, can be expanded as:
Contrary to the previous case, we cannot trivially remove the energy of the resampled image, since as also demonstrated in , , depends on the position of the transformed image relative to the grid, and therefore on Δ. The discrete case of Parseval’s theorem states that the inner product of two discrete signals can be computed in the frequency domain as the inner product of their DTFTs in [−1/2, 1/2]d. Image g is bandlimited and not interpolated, therefore according to the sampling theorem, DTFT (g[p]) = ĝ(v) for ν [−1/2, 1/2]d. Regarding fΔ [·], we saw in Sec. III-A that the Fourier transform of f after being sampled, interpolated, and shifted, is ĥ(v) Σnd(v − n)e−i2π · Δ (Fig. 1c). Resampling, results in replicas of this function at k d (Fig. 1d), i.e.,
for ν [−1/2, 1/2]d, where the bandlimitedness of f was once again exploited. Finally, using the DTFTs of f[·] and g[·] and applying Parseval’s theorem to Eq. (8) leads to:
The CC (first) integral of the above equation is similar to that of the previous case, Eq. (5), with the exception that ĥ*, as opposed to |ĥ|2, appears in the periodic error term, meaning that for a given interpolation kernel, the CC after resampling has more aliasing artifacts, but less smoothing effects. This is because the larger magnitude of ĥ* (compared to |ĥ|2) results in larger unwanted aliased frequency components (ν2 in Fig. 1c), but at the same time prevents too much attenuation of the desired frequencies (ν1 in Fig. 1c). However, as shown in , , the second integral of Eq. (10) is an additional periodic error term that makes the sum cost function (Eq. (10)) deviate from the ground truth (Eq. (3)) even further.1 This error in the energy of the resampled image, which for instance breaks the symmetry inherent to rigid image registration, does not appear when resampling is avoided (Eq. (5)). Note that the aliased frequencies (ν2 in Fig. 1c) are results of using a non-ideal interpolant; aliasing happens even though the original images are assumed bandlimited.
With similar analysis as in Sec. III-A, both sources of error in Eq. (10) can be seen to disappear (become independent of Δ) if the ideal interpolant (sinc) is used. This raises the question of whether employing ideal (or more feasibly, close to ideal) kernels generally eliminates the resampling artifacts in image registration. Although we just proved this to be true in translation-only registration, we will see in the next section that this is in general not the case in image registration, when rotation is involved. We will adopt a rather intuitive analysis, as rigorous comparison of the integral and the sum cost functions with non-ideal interpolation in the presence of rotation would be tedious and outside the scope of this paper.
Images are often represented on radially non-symmetric Cartesian grids. An instance of such asymmetry, as illustrated in Fig. 2 (see also ), would be the fact that a d-dimensional Cartesian grid provides a resolution times higher in the diagonal directions than in the axis directions, which can also be perceived as higher frequency components in the corners of the DTFT domain. Rotation of an interpolated image may align its high-resolution diagonal features with the axes of the grid, which do not have the capacity of representing such high-frequency components. Thus, if the rotated image is resampled on the original grid, aliasing will inherently occur (Fig. 3, left), even if the ideal interpolation has been used –.
This problem can be particularly serious in image registration, where the dependence of the aliasing error in the (sum) cost function on the rotation angle may bias the registration. The area of the aliased regions in the frequency domain (dashed triangles on Fig. 3, left) for a rotation of angle θ can be calculated geometrically and shown to be A(θ) = (secθ − 1) (cscθ − 1), plotted in Fig. 4 for θ [0°, 90°]. Although A(θ) in this interval is maximal at θ = 45°, we cannot simply deduce that the bias is worst at this angle (and at 135°, 225°, etc.), since the minimum, and not the zeros, of the cost function is searched for. This minimum happens at a zero of the derivative of the cost function, in which the error would depend on both and A(θ) and A′(θ), with the magnitude of the latter being largest when the grids are aligned, i.e. at θ = 0°, ± 90°, 180°.
Note again that computing the integral cost function instead of the sum eliminates the resampling step and consequently the resulting aliasing artifacts. Nonetheless, with rotation involved, such an analytical computation becomes much more complicated than the translation-only case of Sec. III-A, as various overlap configurations for the two meshes will have to be considered. Monte-Carlo methods may be employed to approximate the integral, yet a good such approximation requires oversampling the image, thereby increasing the computational cost. (An experimental comparison can be seen in [40, Fig. 12] for small rotation angles, although the authors do not mention the type of aliasing discussed here). Two alternatives to using the integral cost function to avoid aliasing would be: 1) applying the ideal low-pass filter a second time on the rotated interpolated image to clip off the corners , or 2) upsampling the d–dimensional images times in each direction, thus increasing the total number of sample points by a factor of dd/2. Besides the additional computational cost imposed by these approaches to avoid aliasing, one can see that the cost function would still be biased, since the intersection area of the two registered images (in the frequency domain), which is their non-aliased overlapping regions, 1 − A(θ), still depends on the rotation angle. Normalizing the cost function by this area would not in general be helpful, either, since the frequency components are not necessarily distributed uniformly.
It has been suggested in  (outside the registration context) to filter an image with a circular disk (Fig. 3, right) before rotation is applied, to avoid resampling-induced aliasing for any arbitrary angle. Additionally, the corners outside of such a disk have been ignored in phase-correlation registration in , which would be theoretically equivalent to applying such a filter. The authors of  have also used this filter, however as a post-processing step to remove “ringing artifacts” (as opposed to aliasing).
In this work, we propose to make use of the disk-shaped filter in rigid image registration by incorporating it in the interpolation step. We start by identifying the interpolation kernel h with the Fourier transform . Such functions have been computed in ,  for different dimensions; in 2D and 3D spaces they are:
with r: = ||x||, and J1 the first order Bessel function. Such oscillatory interpolation kernels have been considered in  as radial basis functions for interpolation, however, neither for convolution-based interpolation, nor in the registration framework.
Interpolating the image with the IIR kernels in (11)–instead of with the ideal interpolant (sinc)–filters the image with the disk/ball-shaped kernel of a diameter equal to the sampling rate. Such a filter keeps the maximal frequency components in the image, while preventing the occurrence of aliasing at every arbitrary rotation in an unbiased fashion, meaning that the same amount of information is preserved for all the possible rotations. Therefore, with no aliasing coming from resampling, the two sum and integral cost functions theoretically produce the same results.
Nevertheless, for the same practical considerations as mentioned in Sec. III-A, the use of radial windows such as the radial Hanning window becomes indispensable to make FIR filters out of kernels (11). The radial Hanning window, for instance, is defined as , where N is the diameter of the window, 1 the characteristic function, and the distance from the center of the window. The use of the window function, however, makes the kernels non-ideal which may consequently cause aliasing artifacts again as described in Sec. III, making the sum and integral cost functions in practice dissimilar for smaller windows. This is a general disadvantage of the FIR filters; the fact that they cannot exactly produce an arbitrary desired (especially bandlimited) frequency response.
Analytical computation of the integral cost function requires solving definite integrals in irregular d-dimensional regions defined by the intersections of the pixels (or voxels) of the two meshes. This can be complicated particularly when rotation is involved, however, a simple Monte-Carlo approximation of the integral cost function is achieved by oversampling the image on random  or quasi-random  sample points.
Nevertheless, the particular case of translation-only registration with the integral cost function discussed in Sec. III-A (Eq. (5)) can be shown to be equivalent to maximizing the CC after interpolating its discrete version with the self-convolution of the kernel, h*h (see also ). This can be done quickly by computing the CC using the fast Fourier transform approaches (e.g. as in , ), interpolating and upsampling it to the desired resolution, and finding the maximum of the upsampled CC. Otherwise, computing the integral cost function would require dividing each voxel into 2d subvoxels depending on the subvoxel shift between the two images, and subsequently calculating the integral analytically in each subregion. Once the analytical formula for a subregion is derived, the formula for the rest of the subregions can be computed by swapping the two images and shifting them. For instance, if the integral cost function is desired to be computed in the 1D case with 0 ≤ Δ ≤ 1, the pixel [0, 1] needs to be divided into the two subpixel intervals of [0, Δ] and [Δ, 1]. Once the analytical formula I(f[n], g[n], Δ) is derived for the integral of the first interval, the integral of the second interval can be seen to be I(g[n + 1], f[n], 1 − Δ), and therefore does not need to be computed from scratch.
For our experiments on the sphere, we estimated the integral cost function via the Monte-Carlo integral approximation method. We generated spatially-uniform quasi-random Halton sample points on the 2D Cartesian interval of (z, ) [−1, 1]×[0, 2π], using the haltonset command of Matlab. We then transferred the points onto the sphere by computing θ = acos z. One can see that the surface element |sinθ dθd| = |d (cosθ) d| = |dzd| remains unchanged, preserving the uniformity of the sample points.
Regarding the disk-shaped kernels (11) in Sec. IV, a drawback of their employment is the fact that they cannot be implemented in a separable fashion, i.e. as a product of 1D functions, thereby making them computationally more expensive. In fact, besides the Gaussian kernel, no other radially symmetric kernel is separable . To make up for the computational complexity, one might want to pre-filter the image with the disk-shaped kernel offline to remove the corner frequency components and subsequently interpolate it at each iteration using a separable kernel such as bicubic or windowed sinc, the downside of which, yet, would be the introduction of undesirable artifacts (aliasing, blurring, ringing, etc.) twice over, instead of only once.
In Sec. IV, we made the general assumption that images have nonzero components in the corners of the frequency domain, since they are represented on the Cartesian lattice, and that they have not been upsampled to avoid diagonal aliasing. In case the interpolated image is known not to contain any components in the corners of the frequency domain, such as in MRI with spiral trajectories in the k-space, aliasing due to rotation will naturally not occur and the use of disk-shaped kernels is not expected to bring about any improvement.
In this section, we first compare the integral and sum cost functions by visualizing them in 2D self-registration of a synthesized brain magnetic resonance image (MRI) slice taken from the BrainWeb simulated brain database , .2 A T1 image of a normal brain with isotropic 1-mm3 voxels was generated, and a sagittal slice of it with the dimensions 129×129 was taken for the task of 2D translation-only registration. The ideal low-pass filter was then applied to the image which was subsequently downsampled 10 times (to 13×13). Both high-res and low-res versions of the images are shown on Fig. 5(a). Next, the high-res image was registered to itself with the resolution of 1 pixel, the cost function of which is illustrated on Fig. 5(b) with respect to the 2D shift. This cost function is considered as the ground truth for the low-res downsampled data, which were next registered with the resolution of 0.1 pixels using linear interpolation. The sub-pixel registration of the low-res data was performed once using the sum cost function (Fig. 5(c)), and again using the integral cost function by analytically computing the integral of the square of the difference of the continuous interpolated images on the sub-pixel regions created by the overlapping pixels of the fixed and the shifted meshes (Fig. 5(d)). Resampling artifacts can be clearly seen in the sum cost function as periodic piecewise-convex regions. Their convex nature is due to the quadratic convexity of the square of the linearly interpolated image, appearing in the SSD. Such artifacts can, for instance, alter the direction of the gradient in gradient-based optimization approaches, potentially leading to local optima and slower convergence. The artifacts, however, are much smaller and almost invisible in the smooth integral cost function. Furthermore, when we computed the correlation coefficient between each low-res cost function and the ground truth, it was slightly higher (0.3%) for the integral case (r = 0.996) than for the sum case (r = 0.993).
We also looked into the patterns of artifacts in spherical image registration, which has applications in, e.g., registration of the cerebral cortex (see  and the references therein). We generated two separate images of independent and identically distributed white noise on the sphere with a triangular mesh based on the icosahedron subdivision , with 2562 vertices and 5120 triangles. We then computed the sum SSD registration cost function of the two images with barycentric (linear) interpolation while rotating one of them around a fixed axis. In the absence of artifacts, the cost function is expected to be roughly flat superimposed with an irregular noise pattern. However, structured artifacts with a period of 4.6°–about the angular size of the edges–were observed in the cost function (Fig. 6, blue solid curve), most strongly affecting the origin where the two triangular meshes are entirely aligned. Rotation of a sphere acts as rotation close to the poles and as translation close to the equator, which explains why the cost function appears to be a combination of patterns resulted by both rotation with a bias at the grid aligning origin (see for instance [40, Fig. 12]), and translation with periodic artifacts (similar to Fig. 5(c)).
To compute the integral cost function, we resorted to the Monte Carlo approach of stochastic estimation of an integral , and approximated the integral as a sum over randomly distributed points on the sphere. Moreover, since the uniform distribution produces sample points that do not fill the surface as uniformly as expected (due to the clustering effects), we followed the suggestion in  and used Halton sampling  to produce a set of quasi-random but more uniform set of points. Figure 6 shows the results with as many Halton sample points as the number of vertices (1×, green dotted curve) and with five times more points (5×, red dashed curve). As can be seen, there are significantly less artifacts compared to the case with the sum cost function. Additionally, a better approximation of the integral (more sample points) results in a smoother cost function. The standard deviations of the integral–1× and 5×–and the sum cost functions over the rotation angle were 3.3, 2.5, and 16.7, respectively.
We then used the cortical sulcal maps of the five subjects included in the Spherical Demons package , originally from the OASIS public database , to test the spherical registration on real brain MRI data. The maps were projected on icosahedron meshes of orders four (low-res Ic4, 2562 vertices) and six (high-res Ic6, 40962 vertices), for both of which we computed the sum and integral cost functions similarly to the previous case, except that this time images were rotated about two axes, resulting in 2D cost functions in [−10°, 10°]2. Every pair in the five subjects was forward and backward registered, totaling 20 experiments. We found the location of the minimum of each cost function, and used the high-res Ic6 minima as the ground truths for comparing the low-res Ic4 results of the different techniques. The error for each technique was computed as the L2 distance between the locations of cost function minima in Ic4 and Ic6, the average of which across experiments is listed in Table I. As can be seen, low-res registration with Halton instead of regular sampling produces minima that are more compatible with those of the high-res image. In addition, a fivefold increase in the number of Halton samples further improves the accuracy of registration via better approximation of the integral. One can also see that when Halton sampling is used in the ground truth Ic6 registration, the low- and high-res experiments result in closer minima. This is expected, since the artifacts exist even in high-res registration, albeit with lower magnitude, and reducing them via Halton sampling eliminates one source of perturbation in the distance between the minima.
We also measured the inverse-consistency of registration  by computing the L1 distance between the forward and backward registration cost functions corresponding to each pair of subjects, and observed that on average, compared to regular sampling, the asymmetry error is reduced by 36% and 56% using Halton sampling with the same number of and five times more sample points, respectively.
Next, we performed rotation-only registration experiments to validate the interpolation kernel introduced in Sec. IV. We used a high-res 760×760 planar fingerprint image  (Fig. 7, top left) as the ground truth, and each time rotated it with a specific angle (Fig. 7, top right). Subsequently, we applied the ideal low-pass filter to both the original and rotated images and downsampled them 20 times (Fig. 7, bottom). Then we registered the low-res images by an exhaustive search in the vicinity of the known rotation angle with a resolution of 0.1°, using the sum cost function. We performed the experiment three times, making use of the Bessel (h2D, Eq. (11)) and the sinc kernels (both with a 15-tap Hanning window) and also the bicubic interpolation implemented in Matlab (imrotate). We repeated the experiment for 181 rotation angles from 0° to 180°, and noted that the Bessel kernel with a mean error of 0.18° outperformed the bicubic and the sinc interpolations with mean errors of 0.20° and 0.29°, respectively. The corresponding p-values derived using Wilcoxon’s paired signed rank test for the equality of the median of errors were 3 × 10−19, 3×10−2, and 3×10−7, for Bessel vs. sinc, Bessel vs. bicubic, and bicubic vs. sinc, respectively. The histogram of the signed errors is plotted in Fig. 8, where the Bessel kernel can clearly be seen to have a higher concentration of error around zero than the other two methods do.
Bicubic interpolation seems to be more robust than the windowed sinc, which may be explained by the fact that it is equivalent to convolution with the cardinal cubic kernel which approximates the IIR sinc , as opposed to the FIR 15-tap windowed sinc kernel used here. Also note that no matter what the rotation angle is, the low-res images (Fig. 7, bottom) resolve the high frequencies better in the diagonal than in the horizontal and vertical directions, in accordance with our discussion in Sec. IV and Fig. 2.
In this paper, we discussed the artifacts in subpixel registration from a sampling point of view, and showed that resampling the interpolated image is the major cause of these artifacts. We compared the two cases of including and excluding the resampling step in image registration analytically, and also through experiments by calculating the integral cost function via both analytical and stochastic approaches. Computing the objective function in a continuous manner dramatically reduced the amplitude of the induced artifacts, thereby increasing the accuracy of the gradient and the global optimum of the cost function, and also the inverse-consistency of registration. We also demonstrated aliasing errors due to rotation, and proposed to use radially symmetric interpolation kernels to avoid them. We observed significant improvement in registration accuracy by choosing the proposed interpolation kernels over the traditional ones.
Applying the derived methods to functional connectivity MRI registration is part of an ongoing project, where we have noted that not addressing the resampling artifacts guaranteed that almost every registration would converge to the wrong optimum. Future work consists of devising practical implementations with lower computational complexity for the integral cost function.
This work was supported in part by the National Center for Research Resources (P41-RR14075 and the NCRR BIRN Morphometric Project BIRN002, U24RR021382), the National Institute for Biomedical Imaging and Bioengineering (R01EB006758), the National Institute on Aging (AG022381), the National Center for Alternative Medicine (RC1 AT005728-01), the National Institute for Neurological Disorders and Stroke (R01NS052585-01, 1R21NS072652-01, 1R01NS070963), the Shared Instrumentation Grant 1S10RR023401, Grant 1S10RR019307, and Grant 1S10RR023043, the Autism and Dyslexia Project funded by the Ellison Medical Foundation, and the NIH Blueprint for Neuroscience Research (5U01-MH093765), part of the multi-institutional Human Connectome Project. The work of M. R. Sabuncu was supported in part by Harvard Catalyst, under the KL2 Medical Research Investigator Training (MeRIT) Grant, the Harvard Clinical and Translational Science Center under NIH Grant 1KL2RR025757-01 and financial contributions from Harvard University and its affiliated academic health care centers, and the NIH K25 Grant NIBIB 1K25EB013649-01. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Marios S. Pattichis.
Iman Aganj received the B.S. degree in computer science from the École Polytechnique, Paris, France, in 2005, and the Ph.D. in electrical engineering under the guidance of Prof. G. Sapiro from the University of Minnesota, Minneapolis, in 2010. His Ph.D. research was on various problems in medical image processing and analysis, focusing on diffusion-weighted MRI, anatomical MRI, and electron microscopy.
He has been a Research Fellow with the Martinos Center for Biomedical Imaging, Radiology Department, Massachusetts General Hospital, Harvard Medical School, Charlestown, MA, and a Research Affiliate with the Laboratory for Information and Decision Systems, Electrical Engineering and Computer Science Department, Massachusetts Institute of Technology, Cambridge, since 2011, where he is currently involved in research with Prof. B. Fischl’s Laboratory for Computational Neuroimaging and focuses on anatomical and functional MRI image registration.
Boon Thye Thomas Yeo received the B.S. and M.S. degrees from Stanford University, Stanford, CA, in 2002, and the Ph.D. degree from the Massachusetts Institute of Technology, Cambridge, in 2010.
He is currently a Research Fellow with the Duke-NUS Graduate Medical School and is engaged in research on various problems in medical image analysis, including image registration and segmentation. His current research interests include emerging approaches to measure the structural and functional architectures of human brains. By characterizing the brain systems of healthy human subjects, he seeks to understand how these systems support cognition and how they are disrupted as a result of psychiatric disorders.
Dr. Yeo was a recipient of several awards for his research, including the MICCAI Young Scientist Award in 2007, the MICCAI Young Scientist Runner-Up Award in 2008, and the MICCAI Young Investigator Publication Impact Award in 2011.
Mert R. Sabuncu received the Ph.D. degree in electrical engineering, with a focus on image processing and its applications to various medical imaging problems, from Princeton University, Princeton, NJ, in 2006.
He was a Post-Doctoral Researcher with the Massachusetts Institute of Technology (MIT), Cam-bridge, where he was with Prof. P. Golland’s Laboratory for three years and was involved in research on medical image analysis problems in collaboration with researchers from the Harvard Medical School (HMS), Charlestown, MA. He joined HMS as a Faculty Member in 2009, where he is currently an Assistant Professor with the Martinos Center for Biomedical Imaging and leads a research group, which focuses on developing computational algorithms for examining the relationships between medical image data, genotype data, and clinical outcome.
Bruce Fischl is currently involved in research on the development of novel techniques for the generation of models of neuroanatomical structures using MRI, including the automatic construction of geometrically accurate and topologically correct models of the human cerebral cortex, as well as segmentation of over 35 subcortical and ventricular structures. The cortical models have been useful in a number of domains, including the development of a technique that exploits the correlation between cortical folding patterns and function to generate a more accurate mapping across different brains. This high-dimensional nonlinear registration procedure results in a substantial increase in statistical power over more standard methods of intersubject averaging, and allows the automatic labeling of many anatomical features of the cortex. The surface models also yield measures of the thickness of the cortical ribbon, which is of great clinical and research significance as many neurodegenerative diseases result in progressive regionally specific atrophy of the cortical gray matter. The confluence of these two techniques—the construction of highly accurate and topologically correct models of the cortex, and the capability to generate high-resolution mappings across subjects—allows the comparison of cortical atrophy in patient and control population. His research group freely distributes the tools via the FreeSurfer analysis suite to a community of over 12 000 researchers worldwide, which has resulted in hundreds of publications of clinically and neuroscientifically significant findings. His current research interests include a number of studies utilizing these tools to understand the pattern of cortical thinning in healthy aging, development, and disorders such as Huntington’s disease, Alzheimer’s disease, multiple sclerosis, dyslexia, autism, and schizophrenia.
1This error term is however constant for shiftable transforms .
2To avoid algorithmic biases in the results, we only change the interpolation strategy in each experiment while using the same registration algorithm.
Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.
Iman Aganj, Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Harvard Medical School, Charlestown, MA 02129 USA, and also with the Laboratory of Information and Decision Systems, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139 USA.
Boon Thye Thomas Yeo, Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Harvard Medical School, Charlestown, MA 02129 USA, and also with the Department of Neuroscience and Behavioral Disorders, Duke-NUS Graduate Medical School, Singapore.
Mert R. Sabuncu, Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Harvard Medical School, Charlestown, MA 02129 USA, and also with the Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139 USA.
Bruce Fischl, Athinoula A. Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Harvard Medical School, Charlestown, MA 02129 USA, also with the Computer Science and Artificial Intelligence Laboratory, Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA 02139 USA, and also with the Division of Health Sciences and Technology, Massachusetts Institute of Technology, Cambridge, MA 02139 USA.