Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC5741191

Formats

Article sections

Authors

Related links

IEEE Trans Med Imaging. Author manuscript; available in PMC 2017 December 22.

Published in final edited form as:

Published online 2014 June 12. doi: 10.1109/TMI.2014.2329271

PMCID: PMC5741191

NIHMSID: NIHMS925584

Kourosh Jafari-Khouzani, Athinoula A. Martinos Center for Biomedical Imaging, Department of Radiology, Massachusetts General Hospital and Harvard Medical School, Boston, MA, USA;

The publisher's final edited version of this article is available at IEEE Trans Med Imaging

In magnetic resonance imaging (MRI), spatial resolution is limited by several factors such as acquisition time, short physiological phenomena, and organ motion. The acquired image usually has higher resolution in two dimensions (the acquisition plane) in comparison with the third dimension, resulting in highly anisotropic voxel size. Interpolation of these low resolution (LR) images using standard techniques, such as linear or spline interpolation, results in distorted edges in the planes perpendicular to the acquisition plane. This poses limitation on conducting quantitative analyses of LR images, particularly on their voxel-wise analysis and registration. We have proposed a new non-local means feature-based technique that uses structural information of a high resolution (HR) image with a different contrast and interpolates the LR image. In this approach, the similarity between voxels is estimated using a feature vector that characterizes the laminar pattern of the brain structures, resulting in a more accurate similarity measure in comparison with conventional patch-based approach. This technique can be applied to LR images with both anisotropic and isotropic voxel sizes. Experimental results conducted on brain MRI scans of patients with brain tumors, multiple sclerosis, epilepsy, as well as schizophrenic patients and normal controls show that the proposed method is more accurate, requires fewer computations, and thus is significantly faster than a previous state-of-the-art patch-based technique. We also show how the proposed method may be used to upsample regions of interest drawn on LR images.

Image spatial resolution is a limiting factor in quantitative analyses of medical images. In magnetic resonance imaging (MRI), the acquisition time, short physiological phenomena, and organ motion may limit the image resolution. For example, in dynamic contrast enhanced (DCE) and dynamic susceptibility contrast (DSC) MRI of brain, a tracer is injected into the vein during continuous rapid brain imaging. As the arrival of the bolus, its passage through the brain, and its wash-out occur in a short period of time, rapid imaging is required to achieve sufficient temporal resolution to capture the dynamics of the contrast. This poses limitations on the spatial resolution. In other cases, spatial resolution is limited by image acquisition time. For example, fluid attenuated inversion recovery (FLAIR) MRI requires longer acquisition time compared to T1-weighted (T1W) MRI. Thus, spatial resolution may be compromised to achieve reasonable scan time and minimize the likelihood of motion artifact. This is usually done by the acquisition of relatively small number of 2D slices resulting in large slice thickness and spacing between slices. As a result, the images have lower resolution in the slice dimension compared to the other two dimensions (acquisition plane), resulting in highly anisotropic voxels (e.g., 1mm×1mm×6mm).

Quantitative analyses of these low resolution (LR) images are further distorted when these images are transformed into a different space. For instance, voxel-wise analysis usually requires all image volumes of one subject to be aligned and analyzed in one space. Thus, if slices of different image volumes are not already acquired exactly from the same location, they need to be coregistered, which involves interpolation. However, standard interpolation techniques are not accurate and result in distorted edges in the planes perpendicular to the acquisition plane. Fig. 1 shows the result of nearest neighbor (NN), linear, and spline interpolation techniques applied to an image with voxel size 1mm×1mm×6mm. The reason for an inaccurate interpolation is that the neighboring voxels, from which the intensity of an unknown voxel is estimated, may not have the same tissue type as the unknown voxel. This results in an incorrect estimation of the voxel intensity using standard interpolation techniques.

A HR T2W image with 1mm×1mm×1mm resolution was downsampled to 1mm×1mm×6 mm by averaging 6 consecutive slices to simulate partial volume effect and then interpolated using (a) NN, (b) linear, and (c) cubic techniques and **...**

Several methods have been proposed to improve the interpolation of LR images. In one category of these methods, a single LR image from the subject is available for upsampling. Adjacent slices of the LR image may be registered using a non-rigid method to reconstruct the slices between them [1–6]. These registration-based methods best work when each slice is acquired from a thin slab of the organ as blurred edges between thick slices lead to poor coregistration. In [7], a family of adaptive 3D interpolation filters are designed and conditioned on different spatial contexts in order to reconstruct a slice based on four neighboring slices. This method, however, requires a training set. A diffusion-based technique has been used in [8] to reconstruct a slice using neighboring slices. In [9], a sparse super-resolution approach was proposed using overcomplete dictionaries. This method also requires a training set. Finally, in [10], a structure-preserving constraint is enforced by a patch-based reconstruction approach using a variant of non-local means (NLM) filter [11]. This approach is based on the self-similarity concept and this patch-based approach has been also applied successfully in [12] to enhance the resolution of 4D-CT images.

Our proposed method is related to another category of methods, in which a high resolution (HR) image with a different contrast is used to upsample the LR image [13–15]. Conventionally, a HR image, usually a T1W image with isotropic voxel size, is acquired to reveal the brain structure. A LR image may be interpolated with the help of this HR image. This is referred to as brain hallucination by Rousseau in [13] and later the method was elaborated in [14]. Rousseau uses patch-based image interpolation to upsample the LR image, which is based on NLM denoising [11, 16–18]. In NLM denoising, the intensity of a voxel is approximated by the linear combination of all intensities within the image [11]. Rousseau performs this approximation within a neighborhood of each voxel. In other words:

$$\mathbf{x}\left(\mathbf{v}\right)\approx \sum _{\mathbf{k}\in \mathrm{\Omega}\left(\mathbf{v}\right)}w\left(\mathbf{v},\mathbf{k}\right)\mathbf{x}\left(\mathbf{k}\right)$$

(1)

where **x** is the image intensity function and Ω(**v**) is a neighborhood of voxel **v**. The weights *w*(**v**, **k**) are calculated based on the similarity of 3D patches centered in voxels **v** and **k** in the HR image registered to the LR image. The same weights are used for interpolation of the LR image based on the assumption that the local patterns in HR and LR images are similar. This assumption may not be valid – for example when a lesion is clearly visible in the LR image but not in the HR image due to the difference between their contrasts. To resolve this issue, a combination of weights calculated from both HR and reconstructed images is used in an iterative algorithm in [14]. In each iteration, a correlation map showing the inconsistencies between HR and reconstructed images is generated to adjust the enforcement of the weights of the two images. The enforcement of the HR image pattern is lowered when the correlation is poor. In this paper, we refer to this method as the SRIP (super-resolution using intermodality priors) method.

In [15], a similar strategy is used with the difference that the similarity between two voxels of the HR image is based on the similarity of their intensities, as opposed to the similarity of patches centered in those voxels. However, to resolve HR and LR pattern inconsistencies mentioned above as well as any slight misalignment between HR and LR images, a self-similarity strategy is applied by taking into account the similarity between patches centered in two voxels in the reconstructed LR image. The weights of these two similarities (calculated from HR and LR images, respectively) are controlled by adjusting some parameters. In this paper, we refer to this method as the SSIP (self-similarity and image priors) method. SSIP has been shown to outperform SRIP in upsampling synthetic MRI from the BrainWeb database [15, 19].

To perform an accurate interpolation, it is critical to estimate the correct tissue type similarity for the LR image. A patch-based approach to measure the similarity between voxels may be suitable for images with directional textural patterns, which is not an accurate assumption for MR images. The brain gyri and sulci form laminar structures where voxels with similar tissue types form layers (as shown in Fig. 2). A patch-based approach is not rotation invariant and thus two close voxels of one tissue type placed at curved structures may be interpreted dissimilar. This is better shown in Figs. 3 (a)–(b) where the performance of a patch-based approach is illustrated by presenting the similarity map between one pixel (shown by a cross sign) and all other pixels in the image. The pixel is located at the boundary of white matter (WM) and gray matter (GM). As shown, the laminar shape of the brain structure is ignored by a patch-based approach, whereas our approach specifically incorporates this anatomy as a prior in the upsampling procedure. Not only is the assumption of laminar structures valid for normal brain structures, as will be shown by experiments, it also results in more accurate interpolation within lesions compared to a patch-based approach.

Laminar structure of cortex, courtesy of David H. Hubel (http://hubel.med.harvard.edu/). As shown, the cortex is made of layers of cells that form a laminar structure.

Comparison of the performances of similarity measurements by similarity maps. Higher values show more similarity. (a) original image with a point of interest shown by a cross sign, (b) similarity map showing similarity of image pixels with the point of **...**

In this paper, we use the structural information of a HR image to guide the interpolation of a LR image. We propose a novel similarity measure that takes the brain laminar structure into account by assuming that voxels that have equal distances to a strong edge are more likely to belong to the same tissue type. This approach results in higher interpolation accuracy compared to SSIP, a state-of-the-art approach [15], while also requiring fewer computations and being significantly faster. Experiments are conducted on MRI scans of patients with brain tumors, multiple sclerosis, and epilepsy, MRI dataset of schizophrenic patients and normal controls, as well as BrainWeb MRI with multiple sclerosis model [19]. We also test the performance of the algorithm in the presence of noise, intensity non-uniformity, and motion artifact in the HR image, as well as noise and susceptibility artifacts in the LR image. We will also show how the proposed algorithm may be used to upsample regions of interest (ROIs) manually drawn on LR images.

In a standard interpolation technique, the continuous function (image intensity in our problem)–from which a limited number of samples is available–is locally modeled, and the value of the function in one point (voxel in our problem) is estimated by fitting the model to its neighboring points and calculating the value of the model at that point. For example, in linear interpolation, the intensities are assumed to change as a linear combination of the values at some neighboring points, similar to the approximation in equation (1). The weights *w*(**v**, **k**) should be proportional to the similarity between **x**(**v**) and **x**(**k**). Conventionally, *w*(**v**, **k**) correlates with the distance between points **v** and **k**. This is based on the assumption that points closer to each other are more likely to have similar values. This is not always true, particularly in images with highly anisotropic voxels. In the NLM approach, this distance is not the spatial distance between the two points. The similarity between two points may be calculated by a patch-based approach as:

$${w}_{\mathbf{x}}\left(\mathbf{v},\mathbf{k}\right)=\frac{1}{{Z}_{\mathbf{v}}}exp\left(-\frac{{\Vert {P}_{3\mathrm{D}}\left(\mathbf{x}\left(\mathbf{v}\right)\right)-{P}_{3\mathrm{D}}\left(\mathbf{x}\left(\mathbf{k}\right)\right)\Vert}_{2}^{2}}{2{N}_{i}\beta {\widehat{\sigma}}^{2}}\right)$$

(2)

where *P*_{3D}(**x**(**v**)) and *P*_{3D}(**x**(**k**)) are 3D patches of **x** centered in voxels **v** and **k**, respectively [16]. *Z*** _{v}** is a constant of normalization,

$${w}_{\mathbf{x}}\left(\mathbf{v},\mathbf{k}\right)=\frac{1}{N}P\left(\mathbf{v},\mathbf{k}\right)$$

(3)

where *P*(**v**, **k**) is the probability that voxels **v** and **k** have similar tissue types in terms of tissue properties specified by intensity function **x**, and *N* is a normalization factor such that Σ_{k}*w*** _{x}**(

Furthermore, a LR image may not be simply considered as a down-sampled version of a higher resolution image. Partial volume effect and geometric transformation are usually present in the LR images. We consider the following degradation model to take these into account:

$$\mathbf{y}=H\mathbf{x}+\mathbf{n}$$

(4)

where **y** is the LR image, **x** is the original HR image of the same contrast, *H* is a matrix incorporating geometric transformation, partial volume effect, and sub-sampling, and **n** is the observation noise. An approach to estimate the LR image **x** is to minimize a least-square cost function such as:

$$\widehat{\mathbf{x}}=\underset{\mathbf{x}}{argmin}{\Vert \mathbf{y}-H\mathbf{x}\Vert}^{2}$$

(5)

As this problem does not have a unique solution, some prior information is needed to constrain the solutions. This prior information may be applied in the form of equation (1), for example by introducing a regularization term *R*(**x**) and solving:

$$\widehat{\mathbf{x}}=\underset{\mathbf{x}}{argmin}{\Vert \mathbf{y}-H\mathbf{x}\Vert}^{2}+\lambda R\left(\mathbf{x}\right)$$

(6)

where *R*(**x**) incorporates the constraint in (1). Since the assumption is that *H* is known, we would need to select an appropriate *R*(**x**) and a minimization algorithm. We consider the following regularization term:

$$R\left(\mathbf{x}\right)=\sum _{\mathbf{v}}{\left|\mathbf{x}\left(\mathbf{v}\right)-\sum _{\mathbf{k}\in \mathrm{\Omega}\left(\mathbf{v}\right)}w\left(\mathbf{v},\mathbf{k}\right)\mathbf{x}\left(\mathbf{k}\right)\right|}^{2}$$

(7)

This ensures the condition in equation (1) is met [20]. An overview of the proposed algorithm is presented in Fig. 4. An initial estimation of the upsampled image is generated by NN interpolation of the LR image. The HR image is then coregistered to our initial estimate, and the weights *w*** _{x}**(

In this paper, to obtain the values of *w*** _{x}**(

$$P\left(\mathbf{v},\mathbf{k}\right)\propto exp\left(-{\Vert F\left(\mathbf{v}\right)-F\left(\mathbf{k}\right)\Vert}^{2}\right)$$

(8)

where *F* is a feature vector representing tissue property of a voxel. Assuming that **z** represents the intensity of the HR image with a different contrast coregistered to the LR image, *F* could be set to **z** or a vector of **z** values in a neighborhood of the input voxel, e.g., a patch [11, 13–15]. However, such feature vector does not take the brain laminar structure into account and is not rotation-invariant. To capture the laminar structure, voxels having equal distances to the main edges should have similar *F* values. Equivalently, *F* contour maps should have isolines parallel to the main edges. The key point is to propagate the edge information perpendicular to the edge direction. The gradient operator can be used to extract edges. Fig. 5 (b) shows the contour map of a gradient operator applied to the sample image in Fig. 5 (a) (before applying the gradient operator, the image is slightly smoothed). The edge information, however, stays local in a gradient operator. This causes the information of sharp curves to be retained but the edge information is lost as one deviates from the edges.

Performance of the proposed features on a sample image with laminar structure (a), showing contour lines for gradient operator (b), and the same image after Gaussian filtering (c–d) with two different standard deviations with (d) having larger **...**

To further propagate the edge information, we propose using a multi-resolution Gaussian filtering approach. In this approach, the image is filtered by a set of Gaussian kernels with different standard deviations. Figs. 5 (c)–(d) show the contour maps of the filtered image in Fig. 5 (a) using Gaussian kernels with two different standard deviations. As shown, by increasing the standard deviation, the edge information is further propagated but detailed information of smaller structures is lost. A combination of the above operators in the similarity vector *F* may be used in order to retain both high and low resolution information. We propose to use the following feature vector:

$${F}_{\mathbf{z}}\left(\mathbf{v}\right)={\left(\mathbf{z}\left(\mathbf{v}\right),|\nabla \mathbf{z}\left(\mathbf{v}\right){|,\mathbf{z}\ast {\mathbf{g}}_{1}{|}_{\mathbf{v}},\mathbf{z}\ast {\mathbf{g}}_{\mathbf{2}}{|}_{\mathbf{v}},\dots ,\mathbf{z}\ast {\mathbf{g}}_{k}|}_{\mathbf{v}}\right)}^{T}$$

(9)

where is gradient operator, * is convolution operator, **g_{i}**,

An iterative method is used for the minimization in (6). In this method, in iteration *t*+1, the image
${\widehat{\mathbf{x}}}_{\text{rec}}^{t+1}$ is reconstructed by (1) from **$\widehat{x}$*** ^{t}* and the condition in (4) is enforced by:

$${\widehat{\mathbf{x}}}^{t+1}={\widehat{\mathbf{x}}}_{\text{rec}}^{t+1}-\text{NN}\left(H{\widehat{\mathbf{x}}}_{\text{rec}}^{t+1}-\mathbf{y}\right)$$

(10)

where NN is the nearest neighbor interpolation operator. This constraint has been referred to as subsampling consistency [21]. To take into account the inconsistencies between the LR and HR weights, we apply a self-similarity approach based on the concepts from image denoising [16–18]. The values of *w*** _{x}**(

$${w}_{\mathbf{x}}\left(\mathbf{v},\mathbf{k}\right)=\frac{1}{N}exp\left(-{\alpha}_{\mathbf{z}}{\Vert {F}_{\mathbf{z}}\left(\mathbf{v}\right)-{F}_{\mathbf{z}}\left(\mathbf{k}\right)\Vert}^{2}-{\alpha}_{\mathbf{x}}{\Vert {F}_{\mathbf{x}}\left(\mathbf{v}\right)-{F}_{\mathbf{x}}\left(\mathbf{k}\right)\Vert}^{2}\right)$$

(11)

where *α*** _{x}** and

The most time-consuming part of the algorithm is the generation of *w*** _{x}**(

An overview of the algorithm is presented below. First a NN operator interpolates the image to generate an initial estimation of the upsampled image. This is equivalent to repeating each slice *n-*1 times when we would like to divide the slice thickness by *n*. We may want to insert enough number of slices to make the voxels of the resulting image near-isotropic. The HR image is then coregistered to the resulting image and the values of *w*** _{x}**(

- Use NN operator to generate initial estimation
**$\widehat{x}$**^{0} - Coregister the HR image to
**$\widehat{x}$**^{0}and call it**z** - Set
*α*= 0_{x} - Estimate the weighs
*w*(_{x}**v**,**k**) using**z**and**$\widehat{x}$**^{t} - Repeat until ||
**$\widehat{x}$**^{t}^{+1}−**$\widehat{x}$**||^{t}^{2}<*ε*- Correct
**$\widehat{x}$**^{t}^{+1}for downsampling and partial volume effects by applying constraint (10)

- Repeat steps 4–5 with a non-zero
*α*_{x}

Experiments were conducted on artificially downsampled images and clinically acquired LR datasets. Assuming that partial volume effect can be modeled by averaging, downsampling by a factor of *n* in one direction was performed by averaging every *n* slices of the HR image **x**. This is equivalent to the application of operator *H* in (4). The downsampled image was then reconstructed by the HR image **z** and evaluated by two quality measures, namely Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index (SSIM). PSNR is defined by:

$$\text{PSNR}\left(\widehat{\mathbf{x}}\right)=10{\text{log}}_{10}\left(\frac{{d}^{2}}{\text{MSE}\left(\widehat{\mathbf{x}}\right)}\right)$$

(12)

where **$\widehat{x}$** is the reconstructed image,
$\text{MSE}\left(\widehat{\mathbf{x}}\right)=1/\mid \mathrm{\Omega}\mid \xb7{\displaystyle \sum _{\mathbf{v}\in \mathrm{\Omega}}}{\left(\mathbf{x}\left(\mathbf{v}\right)-\widehat{\mathbf{x}}\left(\mathbf{v}\right)\right)}^{2}$, Ω is the region covering the organ, i.e., the whole image excluding background, and *d* is the actual dynamic range of the image, i.e., *d* = max(**x**) − min(**x**). SSIM [22] is a measure more consistent with the human visual system. For images **x** and **$\widehat{x}$**, SSIM is defined as:

$$\text{SSIM}\left(\mathbf{x},\widehat{\mathbf{x}}\right)=\frac{\left(2{\mu}_{\mathbf{x}}{\mu}_{\widehat{\mathbf{x}}}\right)\left(2{\sigma}_{\mathbf{x}\widehat{\mathbf{x}}}+{c}_{2}\right)}{\left({\mu}_{\mathbf{x}}^{2}+{\mu}_{\widehat{\mathbf{x}}}^{2}+{c}_{1}\right)\left({\sigma}_{\mathbf{x}}^{2}+{\sigma}_{\widehat{\mathbf{x}}}^{2}+{c}_{2}\right)}$$

(13)

where *μ*** _{x}** and

Experimentally, *N _{v}* was set to 10 to maximize PSNR and minimize computation time and memory usage for one sample MR image (shown later in the experiments). The effect of parameter

In the following experiments, we evaluated performance of the proposed method in comparison with SSIP (also with multi-threading implementation), a state-of-the-art method. In order to compare the methods in different imaging conditions and in the presence of lesions with different sizes, we selected images with different contrasts and resolutions from three datasets of patients with brain tumors, multiple sclerosis, and epilepsy. The Institutional Review Board approved this study. We also tested the algorithm on publicly available MRI dataset of schizophrenic patients and normal controls provided by National Alliance for Medical Image Computing (NAMIC) as well as BrainWeb MRI. A fixed set of upsampling parameters was used for the proposed method in all experiments. These parameters have been optimized for a single image in the brain tumor dataset. We also evaluated the effect of upsampling parameters, slice thickness, noise, intensity non-uniformity, as well as presence of motion and susceptibility artifacts, and showed the application of the proposed method in upsampling of ROIs. We limited the number of iterations in SSIP method to 30 in all experiments to avoid convergence to a suboptimal solution. This number was selected empirically by a compromise between accuracy and computation time.

We selected 5 patients with brain tumors from our dataset of glioblastoma patients [23, 24]. Each patient had sagittal T2-weighted (T2W) and Gadolinium-enhanced MPRAGE images, both with 1mm isotropic voxel size (matrix size 256×256×176). For each subject, the T2W image was downsampled axially to generate a LR image with the slice thickness of 6 mm, which is a common slice thickness for clinically acquired LR images. The downsampled image was then upsampled using the Gadolinium-enhanced MPRAGE image. We also upsampled the T2W image using NN, linear, cubic, and spline interpolation techniques, as well as SSIP.

Fig. 6 shows the results of upsampling using four different methods, namely NN interpolation, spline interpolation, SSIP, and our proposed method. For a better comparison, a smaller field of view of each image is displayed. As shown, the standard interpolation techniques result in distorted edges. However, SSIP and the proposed method better reconstruct the edges. Fig. 7 shows a magnified section of two images presented in Fig. 6 to better illustrate the difference between the performances of SSIP and the proposed method. As shown, the proposed method has visually superior performance, particularly near the edges. This is mainly because of the better performance of the proposed feature-based similarity measure compared to the patch-based method, which was illustrated in Fig. 3.

Upsampling of T2W images of five subjects from the brain tumor dataset using NN interpolation (row 1), spline interpolation (row 2), SSIP (row 3), and the proposed method (row 4), along with the original T2W image (ground truth, row 5) and HR MPRAGE image **...**

Table I compares the PSNR and SSIM values of different upsampling methods applied to the above images. The table also presents the computation times. The proposed method was significantly faster and more accurate than SSIP. Average computation time for the proposed method was 6 minutes whereas the average computation time for SSIP was 84 minutes. Nevertheless, both SSIP and the proposed methods outperformed standard interpolation techniques.

Comparison of PSNR and SSIM values of different upsampling methods applied to T2W images of brain tumor patients. The proposed method is compared with four standard interpolation techniques as well as SSIP. The computation times have been reported for **...**

Fig. 8 (a) presents the PSNR curve for the image presented in the first column of Fig. 6 (subject 1 in Table I). To better show the trend of changes in PSNR values of the SSIP method, we have presented the PSNR values for 60 iterations. The proposed method calculates the weights only twice, which correspond to the two exponential-like parts of the curve of the proposed method in Fig. 8. The proposed method outperformed SSIP even if the weights were extracted only from the HR image. As shown, SSIP has decreasing performance after seven iterations while the termination condition is not met. This may be due to the problem with iteratively taking the LR image into account explained in the Methods section. SSIP had similar decreasing performance in the rest of images in Fig. 6. Note that each iteration of the SSIP method takes significantly longer than each iteration of the proposed method as SSIP requires the weights to be calculated iteratively. In the case presented in Fig. 8 (a), maximum PSNR for SSIP method was achieved in seven iterations and took about 19 minutes, which is still longer than the 6 minutes using the proposed method.

Comparison of the PSNR values of SSIP and the proposed method. PSNR curves of upsampling the (a) T2W image of subject 1 of brain tumor patients, (b) FLAIR image of subject 1 of multiple sclerosis patients, (c) FLAIR image of subject 1 of epilepsy patients, **...**

To test whether the assumption of laminar structures of the MR images is also more accurate than the patch-based similarity assumption within a lesion, we calculated the values of PSNR and SSIM within the tumors. Table II reports the values. As shown, the proposed method has higher values of PSNR and SSIM in all subjects.

We selected 5 patients with multiple sclerosis. Each patient had a sagittal FLAIR image with 0.5mm×0.5mm×1mm voxel size (matrix size 384×384×160) and a sagittal MPRAGE image with 0.9 mm isotropic voxel size (matrix size 256×256×192). The FLAIR image was downsampled sagittally to generate a LR image with the slice thickness of 6 mm (i.e., 0.5mm×0.5mm×6mm voxel size). Similar to the previous experiment, the downsampled FLAIR image was upsampled by MPRAGE image using the proposed method as well as SSIP and four standard interpolation techniques. Table III compares the PSNR and SSIM values and computation times of different upsampling methods. Both SSIP and the proposed methods outperformed standard interpolation techniques. As shown, the proposed method is significantly faster and is more accurate than SSIP. Fig. 9 compares a sample slice of the FLAIR image of subject 1 in Table 3 upsampled using SSIP and the proposed method. As shown, small lesions that are not visible in the images upsampled by NN and spline techniques are visible in the upsampled image using the proposed method as well as SSIP method, even though the lesions in the MPRAGE and FLAIR images do not always have similar patterns of contrasts. Fig. 8 (b) presents the PSNR curve for subject 1 in Table 3 for SSIP and the proposed method.

Upsampling of FLAIR image of subject 1 from the multiple sclerosis dataset using HR MPRAGE image. First row (left to right): upsampled FLAIR image using NN, spline, and SSIP techniques. Second row (left to right): upsampled image using the proposed method, **...**

We selected 5 patients with epilepsy. Each patient had a sagittal FLAIR image (matrix size 192×192×176), a sagittal T2W image (matrix size 192×192×160), and a sagittal MPRAGE image (matrix size 256×256×176) – all with 1 mm isotropic voxel size. In one experiment, the FLAIR image was downsampled coronally to generate a LR image with slice thickness of 6 mm and then was upsampled using the T2W image. In another experiment, the downsampled FLAIR image was upsampled using the MPRAGE image. Table IV compares the PSNR and SSIM values and computation times of different upsampling methods applied to the above images. Similar to the previous experiments, both SSIP and the proposed method outperformed standard interpolation techniques. However, the proposed method was significantly faster and more accurate than SSIP. In this dataset, we did not find a significant difference in using MPRAGE and T2W images for upsampling the FLAIR image based on PSNR. However, the images upsampled using T2W images had generally higher values of SSIM. Fig. 10 presents the upsampling results using NN interpolation, spline interpolation, SSIP, and the proposed method for subject 1 in Table IV. Although the MPRAGE image had better overall performance compared to T2W image in this subject, T2W image had better reconstruction in the substantia nigra and red nucleus region since it had better contrast in this region. Fig. 8 (c) presents the PSNR curves for this subject using SSIP and the proposed method upsampled by T2W and MPRAGE images.

Upsampling of FLAIR image of subject 1 from the epilepsy dataset using different methods. Row 1: NN interpolation, SSIP using T2W image and SSIP using MPRAGE image. Row 2: spline interpolation, the proposed method using T2W image, and the proposed method **...**

NAMIC Brain Multimodality (http://hdl.handle.net/1926/1687) is a publicly available MRI dataset of 10 schizophrenic patients and 10 normal controls. Each subject has T1W and T2W images with 1 mm isotropic resolution (matrix size 256×256×176). In this experiment, we were interested in evaluating the effects of noise and intensity non-uniformity on upsampling. The T2W image was first downsampled axially to generate a LR image with the slice thickness of 9 mm. We added Rician noise with levels of 0% and 5% (relative to the brain mean intensity) to both LR and HR images. We also applied intensity non-uniformity to the HR image, modeled by a 3D Gaussian function. The downsampled T2W image was then upsampled by the T1W image using the proposed method as well as SSIP and four standard interpolation techniques. We tested the upsampling algorithm with and without noise and intensity non-uniformity, comprising 80 upsampling operations. In each experiment, similar levels of noise were applied to both LR and HR images. For a fair comparison, we denoised the LR images in all methods in the same way the LR images were denoised in SSIP method. Table V presents the mean and standard deviation values of PSNR and SSIM over all 20 subjects. The upsampled images were compared with noise-free images. As shown, the proposed method outperformed all other methods. It also outperformed SSIP in computation time. Furthermore, as expected, noise and intensity non-uniformity generally decreased both PSNR and SSIM values. However, the proposed method was less affected by noise and intensity non-uniformity compared to SSIP. Fig. 11 presents a sample upsampling result with %5 noise level in the LR and HR images and intensity non-uniformity applied to the HR image. Fig. 8 (d) presents the PSNR curves for this image.

Upsampling of a sample T2W image from NAMIC Brain Multimodality dataset. The LR T2W image had 5% noise. The HR T1W image had intensity non-uniformity and 5% noise. First row (from left to right): upsampling using NN, spline, and SSIP methods. Second row **...**

BrainWeb [19] (http://brainweb.bic.mni.mcgill.ca/brainweb/) is a simulated brain MRI database providing two anatomical models: normal and multiple sclerosis (MS). MRI scans have been simulated using T1W, T2W, and proton-density-weighted (PDW) sequences. We were interested in evaluating the performance of the algorithm in upsampling along all three dimensions. We downsampled the T2W and PDW images of multiple sclerosis model to 3 mm, 6 mm, and 9 mm isotropic resolutions and then upsampled by T1W image to 1 mm isotropic resolution. Table VI presents the values of PSNR and SSIM for standard interpolation methods as well as the proposed method. The proposed method outperformed all standard methods. Fig. 12 displays the results of upsampling the T2W and PD images with 9 mm isotropic resolution.

Upsampling of BrainWeb T2W and PD images with 9 mm isotropic resolution to 1 mm isotropic resolution using HR T1W image. First row (from left to right): upsampled T2W images using NN, spline, and the proposed method, the original HR T2W image (ground **...**

Comparison of PSNR and SSIM values of different standard upsampling methods and the proposed method applied to T2W and PDW images of BrainWeb database with epilepsy model. Upsampling was performed along all three directions. The computation times have **...**

In all the above experiments, the LR and HR images had different contrasts (e.g., LR T2W and HR T1W images). However, the method works on images with similar contrasts as well, and indeed due to their contrast similarity, it should result in higher PSNR. In practice, it is uncommon to have LR and HR images from the same contrast. Nevertheless, this may find application in upsampling dynamic images given a static HR image with similar contrast. In one experiment, we upsampled the BrainWeb T2W and PD images with 6 mm isotropic resolution using the T2W and PD images with 1 mm isotropic resolution and obtained PSNR (SSIM) values of 21.78 (0.893) and 21.99 (0.927), respectively, compared to the values of 19.59 (0.807) and 20.01 (0.853) in Table VI.

For the image presented in the first column of Fig. 6 (subject 1 in Table I), we evaluated the effects of the number of voxels *N _{v}*, neighborhood size, number of Gaussian kernels and their standard deviations, parameter

Effect of (a) the number of selected voxels *N*_{v}, (b) neighborhood size, (c) standard deviation of Gaussian kernels, (d) number of Gaussian kernels, (e) threshold *ε*, and (f) slice thickness on upsampling the LR T2W image of subject 1 from the brain **...**

The effect of neighborhood size on PSNR is shown in Fig. 13 (b) for *N _{v}*=10. As shown, when increasing the neighborhood size, the PSNR ratio saturates, although the computation time increases as well. We chose to use a neighborhood size of 7 mm (i.e., a cube of size 7×7×7 around the voxel of interest in the image with 1mm×1mm×1mm resolution). Similar to

To evaluate the effect of the number of Gaussian kernels and their standard deviations, we evaluated six filters with standard deviations of *σ _{n}* = (2

In another experiment on subject 1 from the brain tumor dataset, we changed the threshold *ε* to alter the start point of the first weight update and also repeated step 6 of the algorithm twice, resulting in three weight updates. The PSNR curves with two weight updates are presented in Fig. 13 (e) (PSNR after the third weight update was not shown for better visualization). The final PSNR slightly decreased with increasing *ε*, and the computation time increased due to an additional weight update. However, the third weight update did not notably change the PSNR (less than 0.08 change) and actually slightly decreased PSNR for *ε*= 0.05 and *ε*= 0.03. Thus two weight updates were sufficient to achieve reasonable PSNR.

Fig. 13 (f) shows how performance of the proposed method changes with different slice thicknesses in comparison with SSIP and spline interpolation techniques. As expected, PSNR decreases with increased slice thickness. The difference between the proposed method and spline technique increases at larger slice thicknesses. The same is true for the SSIP method, though it has lower PSNR compared to the proposed method.

We were interested in the performance of our method in the presence of two common artifacts, motion artifact in the HR image and susceptibility artifact causing signal drop and geometric distortion in the LR image. In one experiment (Fig. 14), we selected a HR MPRAGE image with 1mm×1mm×1mm resolution with motion artifact and used it to upsample a LR (downsampled) T2W image with a 1mm×1mm×6mm voxel size. As shown, the result is still visually acceptable, though the structural information in the HR image was limited. It also provided higher values of PSNR and SSIM (PSNR=27.93, SSIM=0.773) compared to NN interpolation (PSNR=26.88, SSIM=0.735), linear (PSNR=27.41, SSIM=0.738), cubic (PSNR=27.72, SSIM=0.755), spline (PSNR=27.77, SSIM=0.757), and SSIP (PSNR=26.25, SSIM=0.710) techniques.

Effect of motion artifact of the HR image on upsampling. First row (left to right): T2W image upsampled using NN (PSNR=26.88, SSIM=0.735), spline (PSNR=27.77, SSIM=0.757), and SSIP (PSNR=26.25, SSIM=0.710) techniques. Second row (left to right): T2W image **...**

In another experiment, we upsampled spin echo (SE) and gradient echo (GE) baseline DSC images (median of 4D images over time) using a HR T2W image (1mm isotropic). The DSC images, and in particular GE DSC, had susceptibility artifact appearing as signal drop and geometric distortion. The resolution of baseline DSC images was 1.2mm×1.2mm×6.5mm, which was upsampled to 1.2mm×1.2mm×1.08mm. The results obtained by NN interpolation, SSIP method, and the proposed method are presented in Fig. 15. The upsampled SE DSC images using SSIP and the proposed method have visually more details than the ones upsampled using NN method. While the upsampled GE DSC images have more details compared to the upsampled ones using NN method, they have some artifacts near the cortex edges due to higher geometric distortion in this image (resulting in inaccurate coregistration), particularly near the temporal lobe. On the other hand, there is inter-slice intensity variation appearing as dark and bright slices in the LR image, which is also present in the upsampled images using both SSIP and the proposed method.

We modeled the partial volume effect with simple averaging. In practice, this assumption may not be valid. Fig. 15 presented the results for upsampling of real clinical images (i.e., not artificially downsampled). In another experiment, we used the proposed algorithm to upsample a FLAIR image with 0.43mm×0.43mm×6mm resolution using a MPRAGE image as well as a T2W image, both with 1mm×1mm×1mm resolution. The upsampled images have 0.43mm×0.43mm×1mm resolution and are displayed in Fig. 16. As shown, the proposed method has a visually superior performance compared to the standard interpolation methods. In particular, some details in the GM and WM boundaries within the cortex are now visible in the upsampled image. Note that the FLAIR image had higher in-plane resolution (0.43mm×0.43mm) compared with the MPRAGE and T2W images (1mm×1mm). Nevertheless, the upsampling algorithm had an acceptable performance. As shown, the choice of HR images had little impact on the final outcome. However, a comparison of the two images’ details reveals that each upsampled image has its own advantages and disadvantages.

The proposed method can be used to upsample ROIs manually drawn on LR images. Note that due to the large number of slices, manual drawing of ROIs in a HR image is a tedious task. On the other hand, a reliable automated tool to outline the regions of interest may not exist. Thus, in practice, the outlines may have to be manually drawn on the LR images. The proposed method could be used to upsample ROIs using a HR image with a similar contrast. The algorithm to upsample ROIs is similar to the one introduced above with two differences: (i) partial volume condition (i.e., (10)) is not enforced as the ROI values are all constant, (ii) the weights are calculated only once using the HR image (i.e., step 6 of the algorithm is eliminated). The resulting image has a range of values which needs to be thresholded at 0.5 (if ROI values are 1) in order to get the final upsampled ROI.

Fig. 17 shows the result of ROI upsampling on two subjects with brain tumors. In one subject, the tumor ROI was drawn on post-contrast T1W image with 0.43mm×0.43mm×6mm resolution and then upsampled using a HR post-contrast MPRAGE image (1mm×1mm×1mm resolution). In another subject, the edema ROI was drawn on FLAIR image with 0.43mm×0.43mm×6mm resolution and then upsampled using a HR T2W image (1mm×1mm×1mm resolution). The results are overlaid on the HR images.

Upsampling ROIs using the proposed method. (a) Sagittal view of a HR MPRAGE postcontrast image showing enhancing tumor coregistered to a LR T1W postcotrast image. The ROIs drawn on the LR image are overlaid on the MPRAGE image. (b) The result of 3D upsampling **...**

Note that segmentation algorithms that work based on the intensity of region may not provide good results here as the region may have varying intensity. However, the proposed upsampling works even if the outlined regions present varying intensity.

A new method was proposed to upsample LR images. The proposed method outperformed SSIP, a state-of-the-art upsampling method, in PSNR, SSIM, and computation time. The entire upsampling takes about 6 minutes for a 256×256×176 image on a quad-core Xeon 5472 3.0GHz CPU. The upsampling parameters were optimized based on the observations in a sample image and were fixed in all other experiments. This showed the method had a robust performance when image resolution or contrast changed. The method successfully worked on the LR images, both with anisotropic and isotropic voxel sizes. In addition, we showed how the method could be applied to upsample ROIs.

We used features based on image intensity, gradient, and multi-resolution Gaussian filtering. Other feature vectors may be proposed that can better estimate the similarity of the voxels on a LR image with laminar structures. However, the proposed features are simple and produce reasonably accurate results in a reasonable time. Although the assumption of having laminar structures is valid in healthy brain tissues, the results showed the method had superior performance in the abnormal tissues as well. The proposed method had decreased but still acceptable performance in the presence of noise, intensity non-uniformity and motion artifact in the HR image, and susceptibility artifact in the LR image. It was also observed that intensity artifacts in the LR image will remain in the upsampled image.

For the majority of subjects, the SSIP method had decreased PSNR with increasing iterations and thus we limited the number of iterations to 30, resulting in an average of 84 minutes of computation time for patients with brain tumors. Our algorithm is faster than the SSIP algorithm because it only requires two updates to the weights instead of each iteration, and this calculation is the most time consuming part. One may want to keep SSIP weights constant to speed up the algorithm, but since the weight data, which is a large amount of data, would need to be saved in the memory, its access may be very slow, resulting in an undesirable decrease in speed.

It is required for the HR image to have some structural similarity to the LR image. In the presence of multiple HR images, we may choose the one that is closest in contrast to the LR image. However, the proposed method may also be extended in future work to combine several HR images to generate more accurate weights. The weights may be generated by multiplying the *P*(**v**, **k**) values of the two HR images.

In the presence of large geometric distortion in the LR image due to susceptibility artifact, the performance of the proposed method decreases. This is because the calculated weights using the HR image do not correspond to the correct voxels in the LR image. To reduce this effect, we may use a non-rigid registration technique to register the HR image to the LR image [25]. The registration needs to be smooth enough for the HR image to avoid following the blocky edges of the LR image interpolated by the NN algorithm (step 2 of the proposed algorithm). After upsampling, the inverse registration may be applied to correct for geometric distortion of the upsampled image. This, however, does not correct the signal drop due to susceptibility artifact. We also showed the performance of the proposed method on the baseline image of LR dynamic MRI scans. The method may also be applied to each time point of the dynamic image but must be fast to reduce the total reconstruction time.

Another application of the proposed method is in the coregistration of LR images. The LR images may be upsampled by a HR image separately and then coregistered to each other. This may have application for the accurate correction of motion artifacts in dynamic MRI scans.

This work was supported in part by the National Institutes of Health under Advanced Multimodal Neuroimaging Training Program Grant R90-DA023427, Grant 5U01CA154601-03, Grant 2P41EB015896-16, Grant R01CA129371, and Grant K24CA125440A.

We would like to thank Tracy Batchelor for providing the MRI data of glioblastoma patients, Caterina Mainero and Sindhuja Tirumalai Govindarajan for providing the MRI data of multiple sclerosis patients, Steven M. Stufflebeam and Naoaki Tanaka for providing the MRI data of epilepsy patients, and Elizabeth Gerstner for providing ROIs in the ROI upsampling experiments. We would like to thank José V. Manjón and colleagues for providing their upsampling codes and helping with using the codes (SSIP method) and Alireza Akhondi-Asl, Ali Gholipour, and Jayashree Kalpathy-Cramer for their useful comments.

1. Frakes DH, Dasi LP, Pekkan K, Kitajima HD, Sundareswaran K, Yoganathan AP, Smith MJ. A new method for registration-based medical image interpolation. IEEE Trans Med Imag. 2008 Mar;27(3):370–7. [PubMed]

2. Penney GP, Schnabel JA, Rueckert D, Viergever MA, Niessen WJ. Registration-based interpolation. IEEE Trans Med Imag. 2004 Jul;23(7):922–926. [PubMed]

3. Ólafsdóttir H, Pedersen H, Hansen MS, Lyksborg M, Hansen MF, Darkner S, Larsen R. Registration-based interpolation applied to cardiac MRI. SPIE Medical Imaging. 2010:762336-762336–8.

4. Laursen LF, Ólafsdóttir H, Bærentzen JA, Hansen MS, Ersbøll BK. Registration-based interpolation real-time volume visualization. Proceedings of the 28th Spring Conference on Computer Graphics; 2012. pp. 15–21.

5. Leng J, Xu G, Zhang Y. Medical image interpolation based on multi-resolution registration. Computers & Mathematics with Applications. 2013 Aug;66(1):1–18.

6. Ólafsdóttir H, Pedersen H, Hansen MS, Larsson H, Larsen R. Improving image registration by correspondence interpolation. Biomedical Imaging: From Nano to Macro, 2011 IEEE International Symposium on; 2011. pp. 1524–1527.

7. Alipour S, Wu X, Shirani S. Context-based interpolation of 3-D medical images. Conference proceedings: Annual International Conference of the IEEE Engineering in Medicine and Biology Society. IEEE Engineering in Medicine and Biology Society. Conference; 2010. pp. 4112–5. [PubMed]

8. Neubert A, Salvado O, Acosta O, Bourgeat P, Fripp J. Constrained reverse diffusion for thick slice interpolation of 3D volumetric MRI images. Computerized medical imaging and graphics : the official journal of the Computerized Medical Imaging Society. 2012 Mar;36:130–8. [PubMed]

9. Rueda A, Malpica N, Romero E. Single-image super-resolution of brain MR images using overcomplete dictionaries. Medical image analysis. 2013 Jan;17:113–32. [PubMed]

10. Manjon JV, Coupe P, Buades A, Fonov V, Louis Collins D, Robles M. Non-local MRI upsampling. Medical image analysis. 2010 Dec;14:784–92. [PubMed]

11. Buades A, Coll B, Morel J-M. A review of image denoising algorithms, with a new one. Multiscale Modeling & Simulation. 2005;4:490–530.

12. Zhang Y, Wu G, Yap P, Feng Q, Lian J, Chen W, Shen D. Hierarchical Patch-Based Sparse Representation-A New Approach for Resolution Enhancement of 4D-CT Lung Data. IEEE Trans Med Imag. 2012 Nov;31(11) [PubMed]

13. Rousseau F. Springer, editor. Computer Vision–ECCV 2008. 2008. Brain hallucination; pp. 497–508.

14. Rousseau F. A non-local approach for image super-resolution using intermodality priors. Medical image analysis. 2010 Aug;14:594–605. [PMC free article] [PubMed]

15. Manjon JV, Coupe P, Buades A, Collins DL, Robles M. MRI superresolution using self-similarity and image priors. International journal of biomedical imaging. 2010;2010:425891. [PMC free article] [PubMed]

16. Coupe P, Yger P, Prima S, Hellier P, Kervrann C, Barillot C. An optimized blockwise nonlocal means denoising filter for 3-D magnetic resonance images. IEEE Trans Med Imag. 2008 Apr;27(4):425–41. [PMC free article] [PubMed]

17. Manjon JV, Carbonell-Caballero J, Lull JJ, Garcia-Marti G, Marti-Bonmati L, Robles M. MRI denoising using non-local means. Medical image analysis. 2008 Aug;12:514–23. [PubMed]

18. Protter M, Elad M, Takeda H, Milanfar P. Generalizing the nonlocal-means to super-resolution reconstruction. IEEE Trans Image Process. 2009 Jan;18(1):36–51. [PubMed]

19. Cocosco CA, Kollokian V, Kwan RK-S, Evans AC. BrainWeb: Online interface to a 3D MRI simulated brain database. NeuroImage, Proceedings of 3-rd International Conference on Functional Mapping of the Human Brain; Copenhagen. 1997.

20. Mignotte M. A non-local regularization strategy for image deconvolution. Pattern Recognition Letters. 2008;29:2206–2212.

21. Banerjee J, Jawahar C. Super-resolution of text images using edge-directed tangent field. Document Analysis Systems, 2008. DAS’08. The Eighth IAPR International Workshop on; 2008. pp. 76–83.

22. Wang Z, Bovik AC, Sheikh HR, Simoncelli EP. Image quality assessment: from error visibility to structural similarity. IEEE Trans Image Process. 2004 Apr;13(4):600–612. [PubMed]

23. Batchelor TT, Sorensen AG, di Tomaso E, Zhang WT, Duda DG, Cohen KS, Kozak KR, Cahill DP, Chen PJ, Zhu M, Ancukiewicz M, Mrugala MM, Plotkin S, Drappatz J, Louis DN, Ivy P, Scadden DT, Benner T, Loeffler JS, Wen PY, Jain RK. AZD2171, a pan-VEGF receptor tyrosine kinase inhibitor, normalizes tumor vasculature and alleviates edema in glioblastoma patients. Cancer cell. 2007 Jan;11:83–95. [PMC free article] [PubMed]

24. Batchelor TT, Duda DG, di Tomaso E, Ancukiewicz M, Plotkin SR, Gerstner E, Eichler AF, Drappatz J, Hochberg FH, Benner T, Louis DN, Cohen KS, Chea H, Exarhopoulos A, Loeffler JS, Moses MA, Ivy P, Sorensen AG, Wen PY, Jain RK. Phase II study of cediranib, an oral pan-vascular endothelial growth factor receptor tyrosine kinase inhibitor, in patients with recurrent glioblastoma. Journal of clinical oncology : official journal of the American Society of Clinical Oncology. 2010 Jun 10;28:2817–23. [PMC free article] [PubMed]

25. Heinrich MP, Jenkinson M, Brady M, Schnabel JA. Robust Super-Resolution Reconstruction with Multi-Modal Registration and Guidance. Proceedings of the Medical Image Understanding and Analysis; Swansea, UK. 2012. pp. 81–86.