Evaluation of Interpolation Effects on Upsampling and Accuracy of Cost Functions-Based Optimized Automatic Image Registration

Interpolation has become a default operation in image processing and medical imaging and is one of the important factors in the success of an intensity-based registration method. Interpolation is needed if the fractional unit of motion is not matched and located on the high resolution (HR) grid. The purpose of this work is to present a systematic evaluation of eight standard interpolation techniques (trilinear, nearest neighbor, cubic Lagrangian, quintic Lagrangian, hepatic Lagrangian, windowed Sinc, B-spline 3rd order, and B-spline 4th order) and to compare the effect of cost functions (least squares (LS), normalized mutual information (NMI), normalized cross correlation (NCC), and correlation ratio (CR)) for optimized automatic image registration (OAIR) on 3D spoiled gradient recalled (SPGR) magnetic resonance images (MRI) of the brain acquired using a 3T GE MR scanner. Subsampling was performed in the axial, sagittal, and coronal directions to emulate three low resolution datasets. Afterwards, the low resolution datasets were upsampled using different interpolation methods, and they were then compared to the high resolution data. The mean squared error, peak signal to noise, joint entropy, and cost functions were computed for quantitative assessment of the method. Magnetic resonance image scans and joint histogram were used for qualitative assessment of the method.


Introduction
1.1. Interpolation. One of the most important parts of designing a registration algorithm is choosing a good interpolation function in order to increase the accuracy of registration. Also, Interpolation is required if the fractional unit of the motion is not matched and located on high resolution (HR) grid. One of the ways by which we can help physicians in coming up with a better diagnosis and treatment is improving the resolution of images. One scheme for the interpolation step is shown in Figure 1. Here, a circle shows the reference HR image, and a diamond and a triangle represent a shifted HR pixel. For instance, if the image is downsampled by a factor of 4, a diamond has (0.25, 0.25) subpixel shift for the vertical and horizontal directions and a triangle has a shift that is less than (0.25, 0.25). In Figure 1, a triangle is not placed on the HR grid and it needs interpolation, but a diamond does not need interpolation. Therefore, some interpolation approaches are proposed to overcome the problem of low resolution in medical imaging. Magnetic resonance imaging (MRI) is an invaluable modality in the medical field. Particularly, neuroimaging with MRI helps physicians to study the internal structure and functionality of the human brain. In these cases, high resolution and isotropic images are important because higher isotropic resolution could theoretically reduce partial volume artifacts, leading to better accuracy/precision in deriving volumetric measurement and decreasing considerable errors in the registration [1]. Clinically, acquiring a fully isotropic 3D image set is not feasible because of time, motion artifacts, and PSNR factors. Thus, typically, in 3D MR data, the inplane direction has higher resolution than the slice direction ( -axis). In this case, invaluable information will be lost in the latter direction. Our objective is to recover and fill in this missing information in order to enable the physicians to have a more accurate perspective of the underlying structure 2 International Journal of Biomedical Imaging The study of interpolation approaches dates back to the 1980s [2]. In which a great diversity of techniques can be found in the literature . For example, B-splines were sometimes referred to as cubic splines [3], whereas cubic interpolation was also known as cubic convolution [4][5][6] and as high resolution spline interpolation [2]. Eight interpolation algorithms are reviewed in the following sections. We first present cubic Lagrangian, quintic Lagrangian, and heptic Lagrangian. Then, we explain a nearest neighbor interpolation approach which is associated with strong aliasing and blurring effect. Next, discussions of the trilinear interpolation approach as well as B-spline 3rd order, B-spline 4th order, and windowed Sinc are explained. Finally, we discuss and evaluate the performance of these interpolation algorithms in order to find the best interpolation method for upsampling of 3D MR images. Different 2D interpolation approaches exist in medical imaging [4]. However, in this paper, we compare the performance (quality and quantity) of eight common interpolation approaches on 3D data.

Nearest Neighbor Interpolation.
Nearest neighbor interpolation (also known as zero-order interpolation) is the simplest method, and strong aliasing and blurring effects are associated with this interpolation [14]. The local 1-point Lagrange interpolation is equivalent to the nearest neighbor interpolation, defined by The images when scaled up in size may look very blocky. Likewise, the local 2-point Lagrange interpolation is equivalent to the linear interpolation, defined by ( , ) = { 1 − | − | , for − 1 ≤ < + 1, 0, otherwise.

Trilinear Interpolation.
Trilinear interpolation calculates values placed between existing voxel values by linearly weighting the eight closest neighboring values. In other words, trilinear is the name given to the process of linearly interpolating points within a 3D box, given the values at the vertices of the box (see Figure 2) [15].
The known values at each vertex are indicated as 000, 100, 010, . . . , 111, and the unknown value is calculated by merging the known corner values weighted by their distance from the point ( , , ) within the cube.

B-Spline Interpolation. B-spline interpolation uses
weighted voxel values in a wider neighborhood compared to trilinear interpolation, but both the B-spline and trilinear kernels are symmetrical and separable. The place of the neighboring points as control points relates to B-spline interpolation and combines the intensity values at these places using a set of polynomial basis according to (5) [16]. Equation (5) shows -order B-spline with + 1 control points ( 1, 2, . . . , ), In (5), , are the polynomial functions of order (degree − 1), and is the number of control points; must be at least 2 (linear) and less than + 1.
( ) is validly defined for min ≤ < max , where min = and max = +2 . A knot vector ( 1 , 2 , . . . , +( +1) ) must be determined. This specifies the values of at which the pieces of curve join, like knots joining bits of string. It is important to note that the degree of the weighting polynomial (the order of the curve) is not dependent on the number of control points, [17]. The weighting polynomial can be recursively defined by the following equation [18]: In (6), ( ) represents an index that refers to the control points and ( ) are generally referred to as knot points (see Figure 3). The series of control point is defined as a control surface. This indexing scheme allows one to weight different control points more than other control points by using it once during the computation. Typically, the first and last control points are weighted more heavily than the internal points to give a smooth interpolating curve. Generally, the shape of the curve ( , ) is specified by the relative spacing between the knots ( 0 , 1 , . . . ,  [3,5,9,[16][17][18][19] = 0 , < , for example, [0, 0, 0, 1, 2, 3, 4, 4, 4] (here = 3, = 5).

Windowed Sinc
Interpolation. This interpolation function has minimum aliasing artifacts in contrast to linear interpolation. Sinc function can be windowed more generally to yield [5,10] the following 4 International Journal of Biomedical Imaging Think of an image data set comprising a 3D matrix voxel with intensities ( , , ), specified by integer position coordinates ( , , ). If one wants to calculate the intensity value at an interior point defined by noninteger coordinates ( , , ), this can be obtained by the following equation [5]: For satisfying (9), two limiting conditions are required: (i) ( , , ) must be band limited. Put differently, the function must have Fourier transform { ( , , )} = ( ) = 0 for | | > for some maximum frequency > 0; (ii) the sampling rate-must be greater than twice the bandwidth, for example, > 2 .
The following section discusses the registration process and explains how interpolation is involved with registration.

Image Registration Algorithms.
Image registration methods in medical imaging seek to align two or more images and can be applied in the same modality on the same patient for the purpose of monitoring and quantifying disease progression over time. Registration can also be applied across different modalities, which is useful for correction of different patient positions across scans, for instance, aligning positron emission tomography (PET) data to an MRI image. Also, image registration can be used on the different patients, which is useful for studies of variability between subjects. Image registration is classified into the following categories and depends on several factors: image modalities (MRI, PET, CT, etc.), the subject of registration (a single person or different persons), the object of registration (head or heart), the image dimensionality (e.g., 2D, 3D, and 4D), and geometrical transformation (affine, rigid, projective, etc.). This study examines 3D affine registration of brain images using voxel intensities similarity measures such as normalized mutual information (NMI), normalized cross correlation (NCC), least squares (LS), and correlation ratio (CR). More explicitly, if a target image is resampled to match a reference image, the image intensities at each voxel should be similar in the two images. In fact, when utilizing an intensity-based cost function, it is essential to repeatedly resample one of the images to match the others at several various resolutions, while searching for the min cost function. This resampling process requires interpolation during the registration process [19]. In OAIR method, interpolation involves resampling of anisotropic voxels in the -direction into isotropic cubic voxels. Also, it is important to note that in the the optimized automatic image registration (OAIR) method, the interpolation technique utilized for registration does not necessarily need to be the same interpolation technique used during registration to compute a final image using the optimal parameters.
In this paper, we are focusing on the effect of interpolation technique and cost function used for intensity-based registration. The following sections are organized as follows, we first give some background and provide a means of defining the critical components involved in image registration and establish a theoretical framework.

Geometric Transformation.
When registering images, one should specify a geometric transformation that specially aligns one image to another. The common transformations can be classified as rigid, affine, and projective. Rigid transformation can be defined as a simple transformation that includes only translation and rotation. The projective transformation is the most general transformation and maps lines to lines (but does not necessarily preserve parallelism). An affine transformation includes scaling, rotation, translation, shearing, and reflection. There are several scanner-produced errors that can result in skewing or scaling terms, and affine transformations are applied to overcome these problems.
An affine transformation maps straight lines to straight lines and keeps the parallelism of lines, but not their lengths or their angles. Changing scaling and shearing factors for each image dimension will extend the degree of freedom (DOF, the number of independent pieces of information that go into the estimate of a parameter) of the rigid transformation [20][21][22][23][24][25][26][27]. Figure 4 shows the five basic components of affine transformations. The following matrices constitute the basic affine transforms in 3D, addressed in homogeneous form.
Translation. Translate a point in the -plane to a new place by adding a vector ( , , ). x = + , y = + , and z = + . P represents scaled matrices: P = TP, where Rotation. If a point ( , , ) is rotated an angle about the coordinate origin to become a new point (x ,y , z ), the three basic rotations in 3D can be defined as follows: International Journal of Biomedical Imaging Rotation about the -axis: ] .
Rotation about the -axis: Rotation about the -axis: ] .
There are several reasons for using homogeneous coordinates, including the ability to apply all four transformations multiplicatively. In view of the fact that transformation combinations (shearing, scaling, and rotation) are all multiplicative transforms, only translation is an additive transform. Next, the following cost functions are defined and described: LS, NCC, CR, and NMI. Furthermore, an overview of the literature on their use in registration for medical applications is included. Among these cost functions, the NMI-based registration has become commonplace in many medical applications [28].

Cost Functions.
The cost function or similarity measure evaluates the similarity between two images. In this section, the behaviors of four commonly used cost functions will be examined.
Least Squares (LS). The least squares method measures the average of the squared difference in image intensities [29]: where is the reference image, is the input image, is the number of values over which the sum is performed, and is the least square. When two images differ only by Gaussian noise, the least squares will be the optimum cost function. Images of two different modalities such as MRI and PET will never differ by only Gaussian noise. Due to patient motion, even two images of the same modality, such as two MRI images, will rarely only differ by Gaussian noise. The effectiveness of LS will be extremely decreased by a small number of voxels having considerable intensity differences.
Correlation Ratio (CR). The main principle of the correlation ratio method is to calculate a "similarity measure" between a reference image and an input image and search for a spatial transformation and an intensity mapping such that by dis-replacing and remapping its intensities, the resulting image ( × ) can be seen as equivalent as possible to . This can be obtained by minimizing the following CR function [30]: (17) which integrates over the voxel positions in the image . The minimum and maximum values for the CR are 0 and 1, respectively. The CR can be applied in multimodal image registration involving positron emission tomography (PET), MRI, and computed tomography (CT) images, providing a good tradeoff between accuracy and robustness [31].
Normalized Cross-Correlation (NCC). The cross-correlation function works very well for aligning images of the same modality. Cross-correlation function is defined by the following equation: where is the reference image intensity, is the input image intensity, and and represent the partials of images and in and -directions, respectively. The summation is taken over the region ( , V), where and overlap. When ( , ) best matches ( , ), CrossCorr( , V) shows the maximum value.
Normalized Mutual Information (NMI). The algorithms of mutual information (MI) have been the most investigated measure for registration of medical image to date. The mutual information of images and is defined by the following [32,33]: where , is the joint probability distribution image of and and and are the marginal probability distribution function of and , respectively. The minimum and maximum values for normalized mutual information are 0 and 1, respectively.
When images are correctly registered and aligned, there is maximal dependence between the gray values of the images, meaning that the amount of mutual information would be high. Misregistration will cause a decrease in the MI measure [34]. NMI has been used with success for a wide variety of combinations, including MR, CT, SPET, PET, and also time series images [35]. NMI can be found in a large number of studies [34,36].

Optimized Automatic Image Registration 3D.
OAIR is a robust image registration algorithm based on FLIRT (FLIRT stands for FMRIB's Linear Registration tool 1.3) [26,[37][38][39]. The OAIR technique specifies a transformation that minimizes a cost function, which represents the quality of alignment between two images. The method assesses the cost function at the number of different image resolutions, starting with the lowest resolution. Each step of increasing resolution uses the previously specified optimal transformation as the starting point and further refines its values. OAIR method usually works very well with the image of the same modality (e.g., MRI-MRI, CT-CT, and PET-PET). During the OAIR registration, the resampling process will influence the computed value of the cost function; therefore, choosing the best interpolation is important.
Outline of the OAIR Method. (1) The registration algorithm specifies the minimum resolution for each dimension of the target and reference images (they are subsampled by factors two, four, and eight).
(4) For each resampled image (which is 8, 4, 2, and 1 times) specifies the transform that minimizes the cost function.
Optimization Steps. The theoretical registration problem is completely determined by an interpolation method, a cost function, and a transformation space. However, in practice, an optimization method is needed to find the transformation that minimizes the cost function [26]. In general, all cost functions require global optimization. As a part of the transformation optimization process, the images are subsampled by several factors (e.g., eight, four, and two times) [38].
Levels Eight, Four, Two, and One Optimization. Reference and target images are interpolated and subsampled by eight, so each image is eight times smaller. The parameters corresponding to the minimum cost function are specified and used as the initial transformation. For the next level (level four) in the optimization, the reference and target images are interpolated and subsampled by four and, like in level eight, the transformation parameters corresponding to the minimum cost function are specified and used as the initial transformation for the next level (level two) in the optimization. For level two optimization, the process repeats, except that the reference and target images are first interpolated and subsampled by factor two. As mentioned above, the parameters of the transformation are systematically varied, and the cost function is assessed for each setting. For level one optimization, 1 mm interpolated images are used and the transformation is generalized to contain 12-DOF. The merit of this multiresolution technique is that the initial optimization, at large , has a noticeably reduced computational load, since the number of sample points is considerably less. Additionally, a large subsampling ( = 8) uses the lowest resolution image and coarse rotation angle, in which the large features of the image are dominate, and so the overall alignment is easier to find. The following sections explain in detail how each resampled image (which is 8, 4, 2, and 1 times) specifies the transform that minimizes the cost function.
(1) Level Eight Optimization. One of the most difficult tasks in image registration is finding the right orientation or rotation, International Journal of Biomedical Imaging where the values represent the amount of rotation in degrees about , , and axes, respectively. In this case, there would be 125 possible angle configurations, since there are three rotation angles and each angle can contain five different values. A 4-DOF local optimization is also applied to find optimal translation and global scale for each angle configuration. The best 20% of the cost values and corresponding angle configurations (candidate local minima) are stored in vector of minima that is used as starting point for a further optimization, which uses a smaller step size over a narrow range of angles. For each parameter setting related to the top 20% of the cost function minima, the algorithm performs the minima over rotation as well as global scale and translation (previously, the algorithm had not optimized over rotation). For each of these sets of parameters, a 7-DOF optimization is then performed, storing the results of the transformation and cost before and after optimization in a vector of minima. A vector of parameters and top 20% of the min cost function values are considered for the next higher resolution (level four optimization) stage, because the relative costs of each candidate solution may change at higher resolutions. The algorithm also uses interpolation to transform images to this new orientation [26,[37][38][39].
(2) Level Four Optimization. The algorithm now calls level four with the interpolated images subsampled by 4 to specify the transformation that minimizes the cost function starting with the transformations determined in level eight. The optimization specifies a 7-DOF transformation that corresponds to minimum value of the cost function. This transformation is then perturbed and the cost function is calculated for these new settings. The perturbations correspond to six degrees for each rotation parameter and a global scaling factor of 0.8, 0.9, 1.0, 1.1, and 1.2. A vector of parameters and the top 20% of the cost function minima values are considered for the next step, which involves images subsampled by 2.
(3) Level Two Optimization. The algorithm uses the images interpolated and subsampled by 2 and computes the value of the cost function for each parameter setting obtained from the level four optimization. It finds the best minimum, and then optimizes it with the 7-DOF, then 9-DOF, and 12-DOF. The algorithm then returns the best minimum after optimization.
(4) Level One Optimization. The algorithm now level one to use the unsubsampled interpolated images and computes the value of the cost function for each parameter setting obtained from the level two optimization. In this step, one optimization run is performed, with the maximum allowable DOF, as determined by the user (max 12-DOF). The best answer is returned from level one and gives us the minimum cost of differences between the images.

Materials and Methods
A 3D spoiled gradient recalled (SPGR). MRI of the brain was acquired at Nationwide Children's Hospital of Columbus, Ohio, USA, using a 3T GE MR scanner from a 34-yearold participant. Interpolation techniques were performed on brain scans. Relevant imaging parameters are listed in Table 1. The first initial reference, also the HR, image dimensions were 512 × 512 × 120 (this is native scanner output) with voxel size of 0.5 × 0.5 × 1.3 mm 3 and with slice thickness and spacing between slices of 1.3 mm (acquiring a fully isotropic 3D scan 8 International Journal of Biomedical Imaging was not feasible because of time, motion artifact, and SNR factors). Because of interpolation and registration time, we simulated the new 3D HR images (simulated reference 2) with a resolution of 256 × 256 × 120 and with a voxel size of 1×1× 1.3 mm 3 and with slice thickness and spacing between slices of 1.3 mm. In the absence of gold standards, simulations are sometimes utilized to assess registration accuracy. A common tactic is to take real data and deform it using appropriate spatial transformation model (affine, rigid, and projective) and other factors that are thought to be relevant in limiting registration accuracy such as simulating the addition of noise and blurring. The first low resolution (LR) images were generated from the reference one, and the resolution was decreased (512 × 512 × 60 and with a voxel size of 0.5 × 0.5 × 2.6 mm 3 ) along the slice direction by subsampling by a factor of 2. The second, third, and fourth LR images were generated from simulated reference 2, and they were subsampled by a factor of 2 in the -, -, and directions. The second LR images were generated with a resolution of 256 × 256 × 60 (axial plane) and with a voxel size of 1 × 1 × 2.6 mm 3 . The third LR images were generated with resolution 256 × 128 × 120 (sagittal plane) and with voxel size of 1 × 2 × 1.3 mm 3 , and finally, the fourth LR images were generated with a resolution of 128 × 256 × 120 (coronal plane) and with a voxel size of 2 × 1 × 1.3 mm 3 . We rotated the LR images in -direction by 5 degrees. Then, we translated the rotated image above in by 2 mm and in by 3 mm. The LR images are corrupted by Gaussian noise (10 standard deviation) and Gaussian blurring (horizontal 5 radius).
Afterward, we used these LR images as input to our interpolation algorithms (trilinear, cubic Lagrangian, quintic Lagrangian, heptic Lagrangian, windowed Sinc, B-spline 3rd order, and B-spline 4th order) to remap to a common size. They were upsampled and back to their original dimension (256 × 256 × 120), and then we compared them to the reference images in order to find the minimum interpolation error during upsampling. Image restoration (adaptive noise reduction and blind deconvolution techniques) was implemented upon the upsampled images to reduce blurring and noise. Adaptive noise reduction algorithm reduces noise without blurring the edges by replacing a pixel value with a weighted sum of all local pixels reached by following a path with small pixel intensity values between neighboring pixels, and blind deconvolution is a method, which allows recovering of the target object from a set of blurred images in the presence or a poorly specified or unknown point spread function (PSF) [40]. Restoration can be implemented by applying any deconvolution method that considers the presence of noise and blurring.
OAIR was applied on high resolution data set (simulated reference 2) with a resolution of 256 × 256 × 120 and with a voxel size of 1 × 1 × 1.3 mm 3 , and the transformed image with a resolution of 256 × 256 × 120 and with a voxel size of 1 × 1 × 1.3 mm 3 . Throughout the OAIR, when an optimal fit was achieved, the target image was reformatted using the transformation function and interpolations described above to match the reference image. For achieving a good registration (intensity-based cost function) between the fixed image (reference image) and the moving image (target image), the resampling was essential because the moving image did not necessarily have the same origin, spacing, and number of pixels as the fixed image. Therefore, the resampling process helped us to have the moving image in the grid of the fixed image.
The intensity-based registration method looked for the transformation that would give the smallest value of the cost function, which we assumed was the transformation that also gave the best alignment. During this registration for analyzing the effect of interpolation and cost function, we applied and tested various interpolations and cost functions. The cost functions which were performed in this method included (1) normalized mutual information (NMI); (2) normalized cross correlation (NCC); (3) least squares (LS); (4) correlation ratio (CR).

Image Assessment.
There are various ways to evaluate the accuracy of registration technique. They can be divided into qualitative and quantitative methods. For qualitative and quantitative assessment of registered images, we proposed five ways to compare and evaluate the new transformation with the old; we needed to quantify the difference between the geometrically transformed source images with the target image.

Quantitative Assessment.
For the quantitative assessment, we considered an a mean squared error (MSE), peak signal to noise ratio (PSNR), and entropy. The MSE and PSNR measures are estimates of the quality of registration images, and entropy is also a suitable choice for quantitative assessment of the accuracy of registration method.
(1) Mean Square Error. MSE was computed between the original image (reference) and reconstructed image in order to measure the average of the squared difference in image intensities: where , , and represent the direct comparison of each coordinate location, is the reference image, and is the International Journal of Biomedical Imaging reconstructed image. The MSE was computed for 3D brain image in order to assign a value and compare the results.
where , , and are the number of points in the -, -,directions, respectively, for the reconstructed volume.
(2) Peak Signal to Noise. PSNR in decibels (dB) between the original image and the registered image is defined by [41] PSNR = 20 × log 10 ( where MAX is the maximum pixel value of the image and RMSE is the square root of the MSE. (3) Entropy. The desire for a measure of information (commonly termed entropy) of a message stems from communication theory [42]. Shannon introduced an adapted measure in 1984 [43], which weights the information per outcome by the probability of that outcome occurring. Given the events occurring with the probabilities, the Shannon entropy is defined as where = (histogram count in bin)/total count. The Shannon entropy can be applied and computed for an image and be used on the distribution of the gray values of the image. An image with a low entropy value has almost a single intensity; it contains very little information. An image with a high entropy value has more or less equal quantities of numerous different intensities; it contains a lot of information [42]. For instance, blurring an image reduces noise and high frequency and thus sharpens the images histogram, resulting in reduced entropy.

Qualitative Assessment.
One way for qualitative assessment is the subtraction of the reference and registered images. Image subtraction techniques can be used to detect subtle changes that may reflect clinically important disease progression [21]. Another way to conduct a qualitative assessment is to create a joint histogram. The joint histogram is a functional tool for visualizing the relationship between the intensities of corresponding voxels in two or more images. Visual assessment is also considered for qualitative assessment.
(1) Joint Histogram. The joint histogram is two-dimensional for two grayscale images A and B and is created by plotting the intensity of each voxel in image A against the intensity of the corresponding voxel in image B. When two images of different modalities are produced, the spatial resolution is likely to be different (see Figure 6). Therefore, before calculating a joint histogram, it is essential to rescale the range of data of the first image to the range of data of the second image. When two images are perfectly aligned, the corresponding anatomical areas overlap, and their joint histogram is highly focused. In misaligned images, anatomical areas are not matched, and they are mixed up and their joint histogram is scattered, for example, the cerebrum region of one image overlaid onto the skull region of another image causes a more dispersed joint histogram. Note, joint histograms are commonly used for different modalities like MR-CT and PET-MR [7,44]. We implemented the joint histogram technique on the registered images of the same modality for checking the effect of interpolations and cost functions on the accuracy of OAIR.

Results and Discussion
In Figure 7, the interpolated images in axial, coronal, and sagittal views are shown for the first LR images. Part (a) illustrates the 3D image interpolated by the trilinear method; in part (b), the image generated by heptic Lagrangian and quintic Lagrangian interpolation images appears in part (c). Part (d), part (e), and part (f) show that 3D images are generated by windowed Sinc, cubic Lagrangian, and nearest neighbor interpolations, respectively. The B-spline 3rd and B-spline 4th interpolations are shown in parts (g) and (h), respectively. The 3D MSEs for 512-sized MRI images were computed in all three planes. To further test the interpolations in 3D, three typical matrix sizes were simulated, namely, 64, 128, and 256. The 3D MSEs of these matrix sizes were tabulated in Table 2. The MSE is inversely proportional to the 3D MRI images size. As a result, the trilinear method yielded more accurate (lower MSE) values than the discussed interpolations. Also, the interpolated images were quantitatively evaluated by computing the PSNR, which is widely used in the evaluation of reconstructed images [11]. The 3D PSNRs for 64-512 sized MRI images were computed. The 3D PSNRs of these matrix sizes were tabulated in Table 3. The PSNR results for trilinear interpolation for matrix size of 64-512 were approximately 91 (dB), 100 (dB), 111 (dB), and 121 (dB), respectively. Based on Table 3, trilinear interpolation shows PSNR superiority against the other interpolation. In addition, the PSNR was found to slowly increase as the matrix size increased. The second LR images, the matrix size of 256 × 256 × 60 (axial view), the third LR images, the matrix size of 256 × 128 × 120 (sagittal view), and the fourth LR images, the matrix size of 128 × 256 × 120 (coronal view) were simply interpolated separately, and MSE was computed. The results are tabulated in Table 4.
As a result, in Table 4, the interpolated matrix size of 256× 256×60 (axial view) yields more accurate (lower MSE) results than both matrix sizes of 256× 128 × 120 (sagittal view), 128 × 256 × 120 (coronal view). In other words, the interpolated images in -direction have more quality (lower MSE) than the interpolated images in -and -directions. However, the perceived quality of 128×256×120 (coronal view) was nearly as good as that of 256 × 128 × 120 (sagittal view). Figure 8, we applied the trilinear algorithm on the second LR images, and the resolution in the axial view did not improve because the resolution in -and -directions was already high (256 × 256 × 60) with voxel size of 1 × 1 × 2.6 mm 3 , and the resolution in axial direction was constant. In contrast, the resolutions in sagittal and coronal views showed different results, and we saw enough improvement in both planes (interpolated images). As was mentioned before, the resolution in the slice-select direction is lower than plane direction, and  we would like to improve the resolution in the slice-selection direction. In Figure 9, we applied a trilinear algorithm on the third LR images with a resolution of 256 × 128 × 120 and a voxel size of 1 × 2 × 1.3 mm 3 ; the resolution in coronal view did not improve because the resolution in -direction was high. However, the resolutions in sagittal and axial views were changed, and we saw enough improvement in their resolutions. The fourth LR images had a resolution of 128×256×120 and a voxel size of 2 × 1 × 1.3 mm 3 ; those results are shown in Figure 10. The resolution in the sagittal view was constant because the resolution in was high and just the sagittal resolution in axial and coronal views was changed. The downsampled results in Figure 8 (middle row) in the sagittal and the coronal views show significant jagged-edge distortion. The trilinear interpolation results in Figure 8 (bottom row) and the coronal views have smoother edges but somewhat blurred appearance overall. Also, the downsampled results in Figure 9 (middle row) in the axial and sagittal views were almost equivalent to the downsampled results in the sagittal and coronal views in Figure 9, and they showed noticeable jagged-edge distortion. The trilinear interpolation results in Figure 9 (bottom row) in the axial and sagittal views had smoother edges but blurred appearance slightly. One can see enough improvement in resolution of both planes (axial and sagittal), but the resolution of the coronal view was constant. The downsampled results in Figure 10 (middle row) were similar to the downsampled results in Figures 8 and 9, but with different perspective (axial and coronal). The trilinear interpolation results in Figure 10 (bottom row) in axial and coronal views showed enough improvement in their planes, but the sagittal view was constant.

Runtime Measurement.
The runtimes of the various interpolation schemes were computed on MR images with a resolution of 256 × 256 × 60 (axial view). In the axial view, the nearest neighbor was the fastest interpolation, with a run time of 16.8 s, and the trilinear was a bit slower than the nearest neighbor with a runtime of 20.4 s. Cubic Lagrangian was fairly fast (35.7 s) and required less time than the quintic Lagrangian, heptic Lagrangian, and windowed Sinc, with 75.2 s, 159.6 s and 182.3 s, respectively. Interpolation with the B-spline 3rd order and B-spline 4th order took about 48 and 185 times as long as nearest neighbor interpolation. This weak performance was caused by evaluation of the exponential function necessary to specify the weights and increasing the order of B-spline showed the interpolation drastically. The results of the run times are presented in Figure 11.
Among the interpolations techniques discussed, the trilinear method was one of the fastest techniques and had the smallest interpolation error. The nearest neighbor had a strong point in which the original voxel intensities were preserved, but the resulting image was degraded significantly and had a blocky appearance. Our experiments showed that the heptic Lagrangian technique had smaller error than the quintic Lagrangian and the cubic Lagragian. The windowed Sinc had a smaller error than the nearest neighbor, B-spline 3rd order, and B-spline 4th order. The main drawback of windowed Sinc interpolation was that, it generated significant ripple artifacts in the surrounding of the images edges. The B-spline 3rd order and B-spline 4th order were the slowest techniques in this study, and B-spline 3rd order produced one of the worst results in terms of similarity to the original image. These results demonstrated that the increment of the order in B-spline will not significantly improve the interpolation quality, and this will just magnify the edge effects and the degree of blurriness, which already noticeable when compared to trilinear and Lagrangian methods. The theory and application of B-spline were analyzed by researchers [45,46], and they found the third-order B-spline interpolator to be sufficient for some specific practical applications [47]. Currently, we believe that the trilinear can offer the best compromise between speed and accuracy in upsampling.

Analyzing the Effect of Interpolation Techniques on Accuracy of Cost Functions-Based OAIR Algorithm.
We implemented OAIR 3D described in Section 1 to perform registration between images, and seven interpolation techniques using similarity measures NCC, LS, CR, and NMI. We computed MSE and PSNR of our results, and the experimental results are listed in Table 5. It is important to note that the interpolation error during upsampling (before registration) is different than the interpolation error of geometric transformation (during registration). For instance, the interpolation algorithm, which has remarkable performance in upsampling process, may have insufficient performance in geometric transformation [48]. Statistical analysis of Table 5 showed that there was insignificant difference between the sets of image registered using CR, LS, NCC, and NMI ( value > 0.9994 for all cost functions). However, the effect of interpolation was considerable, and we observed significant difference between the sets of image registered using different interpolations ( -value < 0.0001 for all interpolations). For instance, sets of images registered using windowed Sinc interpolation were significantly better than sets of images registered using B-spline 3rd-order interpolation with similar cost functions (lower MSE and higher PSNR). For qualitative assessment, we investigated the accuracy of registered results using intensity-based cost functions (CR, LS, NCC, NMI). Windowed Sinc and B-spline 3rd-order interpolations were used during registration (other interpolations schemes can also be used if more investigation is desired). Figure 12 shows axial slices from two registered 3D MRI volumes with their subtractions. The panels show axial slices from two data sets (3D simulated images with a resolution of 256 × 256 × 120 and with a voxel size of 1 × 1 × 1.3 mm 3 ) after registration of the three-dimensional volumes using an intensity-based CR (first column), LS (second column), NCC (third column),    and NMI (fourth column). Windowed Sinc (panel one) and B-spline 3rd-order (panel two) interpolations were used as resampling. In the second row of the panels, after registration, the pixel intensities of the reference and target images were roughly identical and different images were considerably smooth. In both panels, although differences of registered 3D MRI volumes using an intensity-based CR (first column), LS (second column), NCC (third column), and NMI (fourth column) were difficult to see by visual inspection, changing the interpolation showed that significant differences between two panels exist. In panel two, where B-spline 3rd order was used during registration, the boundary of the skull, which was masked out of the images during registration, could still be observed easily in the difference images, whereas in the panel one, where windowed Sinc was used during registration, the skull in different images was not easily detected by eye. A possible reason for checking the effect of interpolations and cost functions for registration can be seen by visual inspection of the joint histogram in Figure 13, which contains several histograms of 3D MRI using an intensity-based CR (first column), LS (second column), NCC (third column), and NMI (fourth column). In the top and bottom rows, windowed Sinc and B-spline 3rd-order interpolations were used during registration, respectively. The top row in Figure 13 showed the joint histogram for 3D images at registration using windowed interpolation and with small amount of misregistration, and there was a diagonal in the distribution with the small dispersion. However, in the bottom row, the B-spline 3rdorder interpolation led to large mis-registration of the image, and increased off-diagonal entries started to appear, and the distribution became more dispersed. The distribution of the bottom row is nonsymmetric, and so the linear relationship is not preserved. In general, the intense inhomogeneity will noticeably change for different interpolations; this is one of the reasons that will induce nonsymmetric dispersion and a nonlinear relationship between intensities. However, the change in the appearance of the histograms for these 3D MRI volumes using CR, LS, NCC, and NMI is insignificant. We used these joint histograms to better understand the effect of different interpolations and cost functions during  Figure 13: Joint histogram for registered 3D MRI volumes using an intensity-based cost function (CR (first column), LS (second column), NCC (third column), and NMI (fourth column)) are shown. The top row is generated from the images when registered using windowed Sinc interpolation, and the bottom row is generated from the images when registered using B-spline 3rd-order interpolation.
registration. We also measured joint entropy of the registered images using various interpolations and cost functions. Because joint entropy is directly related, the joint probability distribution described the statistical relationship of corresponding voxel intensities. Entropy increased with increasing mis-registration as can be seen in visual appearance of the joint histogram (see Figure 13). High dispersion of the joint histogram is equivalent to high joint entropy [49]. The joint entropy results are shown in Table 6, and there are no significant differences between the entropies of registered images using CR, LS, NCC, and NMI ( -value ≥ 0.9999 for all cost functions), which means that they have the similar dispersion in the joint histogram, but there were significant differences between the entropies of registered images using different interpolations ( -value < 0.0001 for all interpolations). Also, we computed costs for different voxel similarity cost functions that were used in registration with different interpolations and cost functions. We used inverse cost functions; thus, the minimum cost function corresponded to better registration. The results are shown in Table 7, and the registered images using CR, LS, NCCI, and NMI yielded very close results (value ≥ 0.9999 for all cost functions), whereas the registered images using different interpolations yielded different results ( -value ≤ 0.0041). Statistical analysis of Tables 5, 6, and 7 showed that interpolations had a significant effect on the registration accuracy, whereas cost functions had no effect on the registration accuracy.

Conclusion
Interpolation techniques play a critical role in the improvement and deterioration of the quality of the image as the resolution changes. Thus, the interpolation error is crucial in assessing the interpolation techniques. The interpolation error depends on features such as geometric deformation and the content of the image; therefore, only one evaluation method would not be adequate to evaluate all properties of an algorithm, and a variety of methods should be applied. The comparison is performed by visual quality assessment, quantitative interpolation error determination, and run time measurement. In this study, the results of the algorithms showed that the trilinear method had the smallest interpolation error and the highest PSNR and was one of the fastest techniques, making it appropriate for upsampling in 3D MR images and super resolution. Although super computers are able to compute a huge amount of data in real time, fast methods might be required for online resampling of image sequences or films [50]. The resulting images for trilinear interpolation were less smooth and blocky than other interpolated images. Nevertheless, trilinear interpolation has the effect of losing some high frequency information from the image [51].
Also, the effect of cost functions (LS, NMI, NCC, and CR), and interpolations (trilinear, cubic Lagrange, quintic Lagrange, heptic Lagrange, windowed Sinc, B-spline 3rd order, and B-spline 4th order) for OAIR of 3D brain images was examined, and our experimental results showed that interpolations can effectively decrease or increase the failure possibility of the registration algorithm, and the robustness of method was not due to the choice of cost function, but the choice of interpolation was critical on the robustness of registration. In addition, each component of the optimization method was also necessary to achieve the accurate registrations. Studying the precise effect of transformation on registrations is an active research area [52] but is beyond the scope of this paper.
We noted that the study presented in the tables and figures relied on the same modality (MRI); in addition, these results are not representative of different modality combinations. Perhaps other conclusions would be obtained by the use of different modality combinations or transformations, for instance, MRI-PET or CT-MRI. Also, OAIR was explained in detail, and it was a powerful image registration algorithm. This algorithm was a fully automated algorithm and proposed various resampling interpolation methods combined with CR, LS, NCC, and NMI as a cost function. One of the advantages of this method was using a feature detector (corners are used as the features) to automatically choose a large number of potentially matchable feature points in both images. The algorithm is able to detect identical features in all projections of the scene regardless of the particular image deformation. This method International Journal of Biomedical Imaging has direct potential for registering clinical MRI images. We have validated this method quantitatively and qualitatively, on the simulated and real data, respectively. Also, there are many excellent sources for more in-depth discussions of image registration that the reader may wish to read and learn [53][54][55][56].

Future Work
We would like to study the precise effect of transformation on registration, and we are also interested in combining the information of three MRI plane orientations using brain images in order to increase the resolution of 3D brain image based on super resolution reconstruction (SRR) technique.