Interpolation-based super-resolution reconstruction: effects of slice thickness

Abstract. Standard clinical magnetic resonance imaging (MRI) is acquired in two-dimensions where the in-plane resolution is higher than the slice select direction. These acquisitions include axial, coronal, and sagittal planes. To date, there have been few attempts to combine the information of these three orthogonal orientations. This paper aims to take advantage of the different in-plane resolution acquired from each plane orientation and combine them into one volume in order to attain a higher resolution image. This combination of MRI data will allow the detection of smaller areas that would otherwise be missed using only one slice orientation. A comparison of slice thicknesses along with image registration is performed. The mean-squared error and peak signal-to-noise were computed for quantitative assessment. MRI and phantom scans and joint histograms were used for qualitative assessment.


Introduction
Super-resolution reconstruction (SRR) is the process of attaining several low-resolution (LR) images and combining them to achieve a high-resolution (HR) image, 1 as shown in Fig. 1.These LR images are usually shifted in space such that one image contains data that the other does not.In medical imaging, the ability to reconstruct HR images is useful to a physician for making a correct diagnosis by better detecting brain tumors or pathological changes.Acquiring multiple scenes can be through the use of one sensor with several acquisitions or from several sensors placed in different positions.This concept is translated to magnetic resonance imaging (MRI) where we utilize LR images from different scanning planes/perspectives (axial, sagittal, and coronal) with focus on unimodal data combination based on a priori knowledge of voxel resolution. 2,3The major advantage of the combination of MRI data from three planes is that it will allow both visualization and data analysis of smaller areas than by using only one plane of orientation.It will also allow for an increase in signal-to-noise ratio (SNR) as a result of combining thicker voxels.This is achievable because the inplane resolution is commonly higher than the slice select direction.
Commonly, brain views are shown in one of three perspectives.The transverse (axial or x-y) planes slice the patient from top to bottom, the sagittal (y-z) planes slice the patient laterally, and the coronal (x-z) planes slice the patient lengthwise from front to back, as shown in Fig. 2. Different slice selection directions are utilized to obtain an image volume.Axial, sagittal, and coronal planes in the brain MRI volume consist of two-dimensional (2-D) slices.Each 2-D image is considered a slice plane, whereas the slice selection direction is oriented along the z-axis.In these cases, HR and isotropic images are important because higher isotropic resolution could theoretically reduce partial volume artifacts, leading to better accuracy in deriving volumetric measurement, and decreasing considerable errors in registration. 4linically, acquiring a fully isotropic three-dimensional (3-D) image set is not feasible because of time, motion artifacts, and SNR factors.Thus typically, in 3-D MRI data, the in-plane direction has a higher resolution than the slice direction (z-axis).In this case, invaluable information will be lost in the latter direction.The objective is to recover and fill in this missing information in order to enable the physicians to obtain a more accurate perspective of the underlying structure available in the data by optimizing the choice of interpolation techniques.

Observation Model
In the absence of gold standards, simulations are sometimes used to assess SRR accuracy.A common tactic is to take real data and deform it using an appropriate spatial transformation model (affine, rigid, and projective) and other factors thought to be relevant in limiting SRR accuracy, such as simulating the addition of noise and blurring.An observation model describes the process of obtaining an LR image from an HR image.The LR image can be obtained from warping, subsampling, blurring, and noise operators executed on the HR image.The observation model can be defined as 5,6 where X is the ideal undegraded HR image, D k represents a decimation operator for the k'th image (subsampling), G k is the geometric transformation operator for the k'th image, B k represents the blur operator of the k'th image, and E k is the random sensor noise.

Super-Resolution Algorithms
There exists a variety of super-resolution methods; two will be discussed, the maximum a-posteriori (MAP) super-resolution and Irani and Peleg's 7,8 method.The method utilized in this paper will follow.It should be noted that there are many great sources for in-depth discussions of reconstruction of HR image from several LR MRI scans. 9-12

MAP Super-resolution image reconstruction
The MAP super-resolution image reconstruction of two orthogonal planes was investigated by Bai et al. 11 They investigated the combination of two orthogonal planes (axial and coronal) in order to create an HR image based on an MAP super-resolution method with improved resolution of in-plane resolution in all directions as well as improved SNR.] SRR techniques typically start with an observation model [see Eq. (1)].An MAP method to find the best fit of X, denoted by X MAP , which is intended to maximize the a-posteriori probability, P, where K represents of the total number of LR observations and it is assumed that the K observations are independent of each other: 11 X MAP ¼ arg max x ½PfXjY 1 ; : : : ; Y k g: (2) The following equation can be obtained by applying Bayes rule, and through some standard derivations: 11 where log PfY k jXg is the log-likelihood function for the k'th observation.log PfXg is the logarithm of the a-priori probability distribution of X.
If it is assumed that typical noise in MR images is Gaussian, the Markov random field (MRF) model as the prior distribution for PfXg can be modeled by the following equation: 11,12 X MAP ¼ arg min x JðXÞ; (4) where where λ is the temperature parameter for the MRF model and ν X ðsÞ is its local potential function.σ 2 k is the standard deviation of Gaussian noise for the k'th observation.

Irani and Peleg's approach
The approach by Irani and Peleg in reconstructing an HR image treats dynamic scenes and more complex motions than static scenes and pure translational motion in the image plane. 6heir algorithm can create a set of simulated LR images.The image differences between the observed LR actual images and simulated LR images are backprojected.The backprojecting kernel can be used as an initial estimate of the HR image.The observed LR images sequences fg k g are obtained from the HR image.The imaging model can be expressed by the following equation: where g k is the k'th observed image, f is the HR image that the super-resolution algorithm is trying to find and T k is the 2-D geometric transformation that transform f to g k .h is a blurring operator, which is specified by the point spread function (PSF) of the sensor.η k is an additive noise term and α k is a downsampling operator. 6he super-resolution algorithm starts creating a higher resolution image with an initial guess f 0 for the HR image, and then the imaging process is simulated to acquire a set of LR images.This simulated set of LR images fg 0 k g K k¼1 corresponds to the set of observed images fg k g K k¼1 . 6The imaging process of g k at the n'th iteration can be expressed by the following equation: Fig. 1 Combining low-resolution (LR) images to achieve a high-resolution (HR) image in order to achieve super-resolution (SR) image construction.
Fig. 2 An illustration of the standard three perspectives of the brain, where axial is also known as the transversal plane.
where ↓ s is a downsampling by a factor s, * is the convolution operator, and n is the n'th iteration.
If the initial estimate image f 0 is the correct HR image, then the simulated LR images fg 0 k g K k¼1 must be equivalent to the observed LR images fg k g K k¼1 .The difference between the images fg k − g ð0Þ k g K k¼1 is computed and used to improve the initial guess image f 0 by backprojecting in order to acquire an HR image f ð1Þ .Each value in the difference image is backprojected onto its receptive field in the initial guess image.The following equation is repeated iteratively to minimize the error function. 6

Interdependence of interpolation, and registration in our SRR technique
In this section, the foundations for this study are described by discussing the interpolation techniques and image registration processes and their relation to the SRR.The aim is to expand the understanding of interpolation and registration and to discover how they are involved with SRR.
Interpolation.Interpolation has become a default operation in image processing and medical imaging and is of the important factors in the success of SRR.Interpolation is needed if the fractional unit of motion is not matched and located on the HR grid.The study of interpolation approaches date back to the 1980s, 13 for which a great diversity of techniques can be found in the literature.For example, B-splines were sometimes referred to as cubic splines, 14 whereas cubic interpolation was also known as cubic convolution [15][16][17] and as HR spline interpolation. 18In our paper, 19 eight interpolation algorithms have been reviewed in more detail, including trilinear, nearest neighbor, cubic Lagrangian, quintic Lagrangian, heptic Lagrangian, B-spline third order, B-spline fourth order, and windowed Sinc.In Sec. 3, we discuss and evaluate the performance of these interpolation algorithms in order to find the best interpolation method for the high accuracy of super-resolution image reconstruction.
Image registration algorithms.This study focuses on an intensity-based registration method.In this registration, interpolation, geometric transformation, and cost function assessment are essential steps, as they can affect the accuracy of registration.This section examines 3-D affine registration of brain images using voxel intensities similarity measures such as normalized mutual information, normalized cross correlation, least squares, and correlation ratio.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 re-sample one of the images to match the other at several various resolutions, while searching for the min cost function.This re-sampling process requires interpolation during the registration process. 20In the optimized automatic image registration (OAIR) method, interpolation involves the re-sampling of anisotropic voxels in the z-direction into isotropic cubic voxels.Also, it is important to note that in the 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.The OAIR technique has been explained in more details in our paper. 19Materials and Methods This work utilizes the concept of acquiring LR images that are not from the same slice planes.Figure 3 shows three scans with the same in-plane resolution but different orientations.Each LR volume is then mapped onto an HR grid based on a priori knowledge of the in-plane resolution.Afterward, three HR-fitted volumes from different planes are combined.

Main Stages for Our Super-Resolution Method
The SR image reconstruction technique consists of four main stages: up-sampling, restoration, registration, and combination.These stages can be performed separately or simultaneously.
Our scheme for super-resolution is illustrated in Fig. 4. The first LR images were generated from simulated HR SPGR MRI of the brain, and the resolution was decreased (256 × 256 × 60 with a voxel size of 1 × 1 × 2.6 mm 3 ) along the slice direction by subsampling by factor of 2 (axial plane).The second and third LR images were generated from simulated HR SPGR MRI of the brain image, and they were subsampled by a factor of 2 in the x and y directions.The second LR images were generated with a resolution of 256 × 128 × 120 (sagittal plane) and with a voxel size of 1 × 2 × 1.3 mm 3 , and finally, the third LR images were generated with a resolution 128 × 256 × 120 (coronal plane) and with a voxel size of After the LR SPGR MRI of the brain images were created and simulated like LR phantom images, they were rotated in the x direction by 5 deg.Then, we translated the rotated image in x by 2 mm and in y by 3 mm.Each LR image is corrupted by Gaussian noise (10 standard deviations) and Gaussian blurring (5-mm radius).
Afterward, the LR MRI images were used as input to interpolation algorithms (trilinear, cubic Lagrangian, quintic Lagrangian, heptic Lagrangian, windowed Sinc, B-spline third order, and B-spline fourth order) to remap to a common size.They were upsampled and changed back to their original dimensions (256 × 256 × 120), and then we compared them to the reference images (simulated HR images) in order to find the best interpolation (quantitative and qualitative assessments).
Image restorations (adaptive noise reduction and blind deconvolution techniques) were implemented on the upsampled MRI of the brain and phantom images to reduce blurring and noise.The adaptive noise reduction consists of three main phases.The first phase of the scheme efficiently determines whether a pixel is noise or not based on some predefined threshold and calculated values.Once pixels are detected as noise in the previous phase, their new value will be estimated and set in the noise reduction phase.Finally, a conditional image enhancement phase will be conducted for those images which have been corrupted with high density noise to preserve the edges and details of the restored image.This allows the adaptive noise reduction to have a great performance even at a high noise density.The major advantage of the adaptive noise reduction algorithm is reducing noise without blurring the edges by replacing a pixel value with a weighted sum of all the local pixels reached by following a path with small pixel intensity values between neighboring pixels.Blind deconvolution is a method for this process and allows recovering of the target object from a set of blurred images in the presence or a poorly specified or  unknown PSF. 21Restoration can be implemented by applying any deconvolution method that considers the presence of noise and blurring.

Implementation of registration algorithm on 3-D MR images
The OAIR was applied on an HR dataset with a resolution of 256 × 256 × 120 and with a voxel size of 1 × 1 × 1.3 mm 3 , and the transformed image had a resolution of 256 × 256 × 120 and a voxel size of 1 × 1 × 1.3 mm 3 .This was repeated for the ACR phantom data.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.

Implementation of our super-resolution algorithm
In the implementation of our super-resolution algorithm, the best registered images with the minimum error in interpolation were selected.This was performed based on the findings that accurate image registration and interpolation are critical in the super-resolution process.As a result, the super-resolution requires the same size images in three perspectives with the minimum error in interpolation (axial, sagittal, and coronal).Image registration was performed after the three planes were mapped onto a common space, and one volume was chosen as the fixed image (axial) and the other two as the moving images (sagittal and coronal).The best interpolation during registration was then quantitatively determined.Then each registered image was mapped onto an HR grid based on a priori knowledge of the in-plane resolution.Afterward, the three HR fitted volumes from different planes were combined.

Image Assessment
There are various ways to evaluate the accuracy of interpolation techniques, registration, and SRR.They can be divided into qualitative and quantitative methods.

Quantitative assessment
For the quantitative assessment, we considered a mean square error (MSE) and peak-signal-to-noise ratio (PSNR).The MSE and PSNR measures are estimates of the quality of registration images.
Mean square error.The MSE was computed between the original image (reference) and registered image in order to measure the average of the squared difference in image intensities: where I, j, k represent the direct comparison of each coordinate location, R is the reference image, and I is the reconstructed image.The MSE was computed for a 3-D brain image in order to assign a value and compare the results.
where n, m, l are the number of points in the x, y, z directions, respectively, for the reconstructed volume.
Peak signal-to-noise.The PSNR in decibels (dB) between the original image and the registered image is defined by 22 PSNR ¼ 20 × log 10 MAX RMSE where MAX is the maximum pixel value of the image and RMSE is the square root of the MSE.

Qualitative assessment
One way for 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.
Joint histogram.The joint histogram is 2-D 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.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, they are mixed up and their joint histogram is scattered.For example, images of the cerebrum over the skull cause a more dispersed joint histogram.
Example joint histograms for different modalities like magnetic resonance -computed tomography and positron emission tomography -magnetic resonance at different stages have been investigated in some papers. 23,24We implemented the joint histogram technique on the registered images for checking the effect of interpolations on the accuracy of OAIR.

Evaluating Effect of Slice Thickness on Super-Resolution Technique
The 3-D MSE and PSNR of LR datasets were computed.The 3-D MSE and PSNR of the combined (all three planes) SPGR phantom datasets with different slice thicknesses were tabulated in Table 2.The MSE is proportional to the slice thickness.As a result, a thinner slice thickness has a lower MSE value.
The SPGR phantom datasets were also quantitatively evaluated by computing the PSNR, which is widely used in the evaluation of reconstructed images. 20The 3-D PSNRs for different slice thickness phantoms were computed.As can be seen in Table 2, the combined LR SPGR phantom 1 with a smaller slice thickness shows the PSNR superiority against the other combined LR SPGR with a bigger slice thickness.In addition,   the PSNR was found to slowly increase as slice thickness decreased.The smaller slice thickness yielded a lower MSE error and higher PSNR.Statistical analysis of Table 2 showed that there was an insignificant difference between the sets of images registered with the minimum difference slice thickness (P-value >0.78).But it was observed that a significant difference was between the set of images registered with the maximum slice thickness (i.e., combined LR1 and combined LR4) (P-value <0.0001).For instance, sets of images of combined LR SPGR phantom 1 and 2 (smaller slice thickness) were significantly better than the sets of images of the combined LR SPGR phantom 3 and 4 (bigger slice thickness).In addition, the 3-D MSE and PSNR of SPGR MRI of the brain datasets tabulated and used for the validation of the original information preservation method were anatomical volumes in the axial (256 × 256 × 128), sagittal (256 × 120 × 256), and coronal (120 × 256 × 256) planes (see Table 3).Each plane was mapped to a 256 × 256 × 256 grid preserving original data in order to obtain the same size image, then they were combined with each other.Our method yielded more accurate results than separately interpolating each plane.As Table 3 shows

Visual Quality of Reconstructed Images
Combining the datasets would lead to mis-mapping of data.This misalignment was ameliorated by applying image registration.The axial, coronal, and sagittal planes are demonstrated in Figs. 5, 6, and 7 for the LR SPGR phantom dataset and combined LR SPGR phantom dataset and MRI of the brain, respectively.
In Fig. 5, resolution of SPGR phantoms in three perspectives were decreased by increasing the slice thickness (top to bottom).We also computed the joint histograms of axial, sagittal, and coronal planes for the SPGR phantom and MRI of the brain along the 3-D combined datasets in order to evaluate the accuracy of the SRR method (see Figs 5-7).As is clearly seen in Fig. 6, combined LR1 and LR2 datasets (smaller slice thickness) have a lower dispersion in the joint histogram compared to the combined LR3 and LR4 datasets (bigger slice thickness).Also, the combined LR datasets yielded more accurate and more highly resolved results than the LR datasets without the combination (Figs. 5 and 6).
Visually inspecting Fig. 7 for each scan demonstrates higher quality images for each orientation, respectively, and a lower quality for the other two methods.This was the case in all Fig. 6 SPGR phantom images of the axial, coronal, sagittal, and combined datasets (top to bottom).Each column displays the axial, coronal, and sagittal perspectives (left to right).Joint histograms for the combined LR datasets are similar to the HR SPGR phantom dataset and all gray value correspondences lie on the diagonal with minor dispersion; however, in the first and second top rows, all three perspectives have a higher resolution compared to other two datasets (LR3 and LR4).three planes; the quality was better than just in the axial, sagittal, and coronal planes.For instance, the top row (axial plane) demonstrated better quality for the first row (axial perspective) and worse quality for the other two.
We have introduced our method to exploit the original information from MRI data using different planes to increase the resolution and SNR.Thus, by scanning in two or three different planes, we are adding new information, specifically since the in-plane resolution (xy) is usually higher than the slice thickness (z).Note, if slice gaps are not present, the voxels are averaged from either two or three planes, taking advantage of the higher SNR as a result of the slice thickness.This is recommended since small voxels have a lower SNR, so the high spatial resolution image may be too noisy to be diagnostically useful if the voxels are too small for an adequate signal.Also, it is known that small voxels are particularly a problem on low-field imaging systems.This problem is alleviated by our technique since we combined information from the in-plane resolution of another plane and thus do not have to decrease the voxel size at the time of acquisition.

Discussion
Theoretically, acquiring three volumes with thicker slices then generating thin slice images is practical as this will have a gain in SNR due to the thick slices, and thus a gain in the total acquisition time.By scanning in two or three different planes new information is added, specifically since the in-plane resolution (xy) is usually higher than the slice thickness (z).This is recommended since small voxels have a lower SNR so the high spatial resolution image may be too noisy to be diagnostically useful if the voxels are too small for an adequate signal.Also having LR volumes results in the partial volume effect, which arises when the interface between two different tissues occurs within a single voxel. 25It is also known that small voxels are particularly a problem with low field imaging systems.This problem is alleviated by our technique since we combine information from the in-plane resolution of another plane and do not have to decrease the voxel size at the time of acquisition.In SRR, each LR input image focuses on a slightly shifted field of view (FOV) of the HR scene and does increase the resolution. 25However, in our case, the FOV is not shifted, but is completely rotated onto a separate plane allowing for HR voxels from three planes as opposed to one.In turn, this allows for an increase in resolution in all three directions (3-D) as opposed to only the in-plane (2-D) direction.SRR is a process of combining several LR images to create an HR image.Thus, our method is novel in that the images are not purely LR because the in plane resolutions are HR.We are combining HR data from three perspectives or planes.
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.
According to the above discussion, our new method has several advantages when compared to the MAP-Huber method.These advantages are listed below.1.Our method is relatively simple and very intuitive to implement, whereas the MAP-Huber method results in a more complicated optimization problem.
2. The MAP-Huber method creates simple blur and white noise.
3. Our proposed SRR method has a more powerful inclusion of a priori information.
A disadvantage, however, would be the increase in computing time and it is important to develop efficient algorithms to reduce computational costs.

Conclusion
The goal of this paper was to assess the effects of slice thicknesses on the combination of the orthogonal plane MRI SPGR brain and SPGR ACR phantom scans.This combination can be utilized for reducing artifact data loss by utilizing nonisotropic volumes of the same object and combining them to extract the most information from all the volumes.This ability to combine information will enrich the data sets by both enhancing the image for visualization and the data for further computational processing.Future extensions of this study will include the investigation of better registration techniques and the improvement of computational speed.

Fig. 4
Fig. 4 Scheme for super-resolution.Y 1, Y 2, and Y 3 are LR images.X is an HR image.

Fig. 3
Fig. 3 The proposed scheme includes three stages: (a) acquiring LR volumes, (b) mapping LR volumes to a common HR grid, and (c) combining and interpolating the registered LR volumes.

Fig. 5
Fig. 5 Spoiled gradient recalled (SPGR) phantom images of the axial, coronal, and sagittal with different slice thickness (top to bottom).Each column displays the axial, coronal, and sagittal perspectives (left to right).The joint histogram for each LR dataset is dissimilar to the HR SPGR phantom and LR1, LR2, LR3, and LR4 show the histograms with a significant dispersion.
, combining the SPGR MRI of the brain image from the three planes has a lower MSE in comparison with the individual volumes.The runtimes of each plane and the combined SPGR MRI of the brain were computed (run times measured on the Intel Xeon with 2.13 GHz 2 processor).The interpolation and registration of the axial view (256 × 256 × 128) took ∼184 s.The sagittal and coronal views required around 190 and 193 s for 3-D MR images of 256 × 120 × 256 and 120 × 256 × 256, respectively.Compared with the axial view took 4.89% of the coronal view time.The reconstruction of the combined data took about 1.5 times as long as the reconstruction of the axial, sagittal, and coronal views individually.

Fig. 7
Fig.7Brain images of the axial, coronal, sagittal, and combined datasets (top to bottom).Each column displays the axial, sagittal, and coronal perspectives (left to right).Joint histogram for the combined datasets is similar to reference datasets and all gray value correspondences lie on the diagonal with minor dispersion: axial, coronal, and sagittal show the histograms with significant dispersion.

Table 2 3
-D mean square error (MSE) and peak-signal-to-noise (PSNR) for SPGR phantom with different slice thickness.

Table 3 3
-D MSE and PSNR for MR images of 256 for axial, coronal, sagittal, and combined three volumes.