Multi-contrast brain magnetic resonance image super-resolution using the local weight similarity

Background Low-resolution images may be acquired in magnetic resonance imaging (MRI) due to limited data acquisition time or other physical constraints, and their resolutions can be improved with super-resolution methods. Since MRI can offer images of an object with different contrasts, e.g., T1-weighted or T2-weighted, the shared information between inter-contrast images can be used to benefit super-resolution. Methods In this study, an MRI image super-resolution approach to enhance in-plane resolution is proposed by exploring the statistical information estimated from another contrast MRI image that shares similar anatomical structures. We assume some edge structures are shown both in T1-weighted and T2-weighted MRI brain images acquired of the same subject, and the proposed approach aims to recover such kind of structures to generate a high-resolution image from its low-resolution counterpart. Results The statistical information produces a local weight of image that are found to be nearly invariant to the image contrast and thus this weight can be used to transfer the shared information from one contrast to another. We analyze this property with comprehensive mathematics as well as numerical experiments. Conclusion Experimental results demonstrate that the image quality of low-resolution images can be remarkably improved with the proposed method if this weight is borrowed from a high resolution image with another contrast. Graphical Abstract Multi-contrast MRI Image Super-resolution with Contrast-invariant Regression Weights Electronic supplementary material The online version of this article (doi:10.1186/s12880-016-0176-2) contains supplementary material, which is available to authorized users.


Background
In MRI, low-resolution (LR) images may be acquired in applications, e.g., functional MRI [1,2] and diffusion tensor imaging [3,4], due to limited data acquisition time or other physical constraints. High-resolution (HR) images appear favorable to perform subsequent posterior image processing and visualization [5]. Super-resolution methods are widely utilized to improve image resolution [6][7][8][9][10]. Typical methods include sparse representations [6][7][8], projection onto convex sets (POCS) [9], tensor frames [10], etc. However, these methods need numerous iterations to accomplish super-resolution, thus they inevitably lead to high computational costs. For MRI, since a great number of images have to be processed, fast and stable methods are desired. Recently, the prior information of MRI has been explored in superresolution. For example, (a) redundant information produced by sub-pixel spatial shifts between multiple images [3], (b) space homogeneity constraint from orthogonal anisotropic acquisitions [2], and (c) the learned dictionary with a nature of the orthogonality [11] have been employed to refine structural details and edges. Besides, image contrast can also be utilized to produce sharper images [12]. However, these methods may not lead to faithful super-resolution results when multiple-shifted images are inapplicable or the information is very limited within a single image. Thus, one may expect other prior information beyond a single image.
Multi-contrast images are frequently acquired in MRI experiments [13]. For example, plentiful edge structures are visible both in T1-weighted and T2-weighted brain images of the same subject. According to the principles of MRI [14], we pick up T1 or T2 weighted signal denoted by SI and take the form where ρ(H) refers to the proton density, TR is the repetition time and TE is the echo time. There are different TR value and TE value within a section of medical tissue that would result in multiple contrast images. Yet, these images share the proton density of the subject so that they largely share similar anatomical structures but with different contrasts in regions. The shared information between inter-contrast images can be considered to benefit super-resolution. Therefore, it is possible to improve the LR image resolution by incorporating prior information from the different contrast image in HR. Rousseau proposed a patch-based iterative framework combining with non-local similarity to share information among multiple contrast images in [15], and later many more detailed analysis was studied in [16]. A constraint that the downsampled version of the reconstructed LR data must be equal to the original LR data is imposed in the iterative framework [5]. The non-local similarity is also measured with both voxel intensity and gradient intensity in super-resolution [17]. However, these methods require training sets or time-consuming iteration processing.
New edge-directed interpolation (NEDI) [18] is a fast and statistical super-resolution method for a single image. It estimates local covariance coefficients from a LR image and assumes that this statistical information is also valid for the corresponding HR image. A pixel of the HR image is interpolated by performing the linear regression of neighboring pixels, which originate from the LR image. This regression process is based on noniterative operations, thus the super-resolution can be performed fast. The NEDI provides a nice way of analyzing statistical information in the image super-resolution. Some recent methods [19][20][21] also use regressions to improve the image resolution and achieve remarkable performances. However, these methods train hundreds of external images prior to recovering structural details, and require plenty of computations. Due to the nice statistical property and low computation time of NEDI, in this work, we extend it into the multi-contrast image super-resolution and demonstrate its superior performance on MRI images.
We will explore how to incorporate the statistics from one image into another contrast image. Regression weights, estimated from a HR image in one contrast, and neighboring pixels around the interpolated location in the LR image of another contrast work together to generate a new pixel value. The fact that neighbors are provided by the LR image itself can offer a guarantee and support for the consistent contrast between the LR one and the interpolated result. Mathematical . a-f share the same structure but have different intensities analysis and experimental evidence will be presented to address a fundamental question of why these weights between two contrast images constitute faithful criteria. Then, the proposed approach probes the information both from a LR image and its corresponding HR image in another contrast. Our method will be compared with the classic bicubic method, NEDI method [18], and the state-of-the-art contrast-guided interpolation (CGI) method [12] in terms of objective-evaluation criteria and visual perceptions.
The remainder of this article is organized as follows: In section II, we briefly review basic concepts of NEDI. In section III, we derive conditions that must be satisfied in our method. Experimental results and discussions will be presented in sections IV. Finally, concluding remarks are made in section V.

Brief review of NEDI
In NEDI, regression weights are estimated in a local region then target pixels are calculated as a linear regression of neighbors [18]. Thus, it is crucial to determine the regression weights in the interpolation. Within a neighborhood, four neighbors are commonly used in NEDI, and consequently there are four regression weights for one pixel interpolation.
The interpolation process is shown in Fig. 1. The NEDI uses patches in the local region to estimate regression weights b j (j = 1, 2, 3, 4) (Fig. 1a). The variable n (i = 1, ⋅ ⋅ ⋅, n) denotes the number of patches and each patch is composed of one pixel y i and its four neighbors x i,j along diagonal directions. Then, the target pixel γ is obtained by multiplying neighbors and their weights (Fig. 1b). The basic regression model (Fig. 1a) applied in our work is where ε i is the residual error. By continually sampling in a 9 × 9 region, a vector y = [y 1 , ⋯, y 49 ] T ∈ ℝ 49 is formed to represent pixels in this region and meanwhile a matrix X = [x 1 , ⋯, x 49 ] ∈ ℝ 49 × 4 , whose column x i contains four neighbors of y i , is formed to represent all neighboring pixels around those pixels of y.
Assuming the image pixel values in a local region satisfy a locally stationary Gaussian process [18], the regression and its solution is The above analysis can be also interpreted from the classical Wiener filtering theory. Let R = (X T X) − 1 ∈ ℝ 4 × 4 represents a covariance between two arbitrary members of the four nearest neighbors, r = X T y ∈ ℝ 4 represents a covariance between the center-pixel and the one of the four nearest neighbors around it, the optimal coefficients can be found by Multi-contrast image super-resolution In the proposed method, a HR image of one contrast is assumed to be available for interpolating a LR image of another contrast. This assumption is reasonable since multi-contrast images are always available in MRI experiments [5,7,13]. The regression weights b i for the i th pixel, borrowed from one contrast HR image according to Eq. (4), is incorporated into the interpolation of the LR image in another contrast. Interpolated pixels ỹ i of an expected HR image are given bỹ where the vector s i includes four pixels of the LR image that are the nearest neighbors along diagonal directions of the i th pixel in the center. This means we assume that the HR image in Fig. 1a is in one contrast and the LR image in Fig. 1b is in another contrast. Then b i is estimated from Fig. 1a and s i comes from    To facilitate following discussion, intensities of images are all normalized between 0 and 1. Furthermore, we assume that multi-contrast images are well registered before super-resolution.

Weights in multi-contrast images
For example, multi-contrast images (Fig. 2) share similar anatomical structures but are with different intensities in sub-regions.
An interesting phenomenon is that, regression weights for different contrast images in Fig. 2 are nearly the same ( Table 1). The same observation is also found ( Table 2) for MRI images generated from the BrainWeb [22] that embody more complex structures (Figs. 3a-d). However, regression weights (Table 3) will be totally different if images do not share the similar anatomical structures (Figs. 3a, e-g). These observations convey important information: The regression weights obtained using the least square estimation is nearly invariant to image contrasts. If this is possible, one may easily employ the information from another contrast image by making use of these weights.
Besides, one may find that the sum of weights in each vector is approximately 1 (Tables 1, 2 and 3). We will analyze this property with comprehensive mathematics and empirical tests on MRI images. This property will be an important foundation to derive similar regression weights for multi-contrast images.

Sum of weights is approximately equal to 1
Suppose there are n central pixels, by adding n operations in a local region, Eq. (2) is written as X n Here, ε i is assumed to satisfy the normal distribution, i.e., ε i~N (μ i , σ 2 ). The variable μ i is the mean and σ 2 is the variance associated with ε i . Then we can easily have where there exists ε i ′ ∼ N(0, σ 2 ). Next, according to the principle of the law of large number, meaning that sufficient central pixels are sampled, one has X n where ∑ i = 1 n μ i is a fairly small constant. Then, given that 4 and ∑ i = 1 n y i are equal to one another, one can obtain that the sum of weights follows We verify this property that sum of weights is approximately equal to 1 on MRI images. Statistical analysis in Fig. 4 show that most of ∑ j = 1 4 b j are very close to 1 for  tested images. In each image, the range of ∑ j = 1 4 b j lies between 0.95 and 1.05 can cover above 95% pixels of local regions.
An explanation on why sum of weights is nearly 1 is given. As shown in Fig. 5a, the red solid wireframe indicates the local region of size 9 × 9. Inside this region, all upper left pixels x j (j = 1) come from the pixels in the marked region X 1 in Fig. 5c. In the same way, the upper right x j (j = 2), bottom left x j (j = 3) and bottom right x j (j = 4) will be from  Fig. 3a and b). c-d list the frequency of ∑ j = 1 4 b j for real images (i.e., Fig. 9a and b) Fig. 5 An illustration of ∑ j = 1 4 b j ≈ 1. a A synthetic image of size 256 × 256 in which the red solid wireframe draws out a local region of size 9 × 9; b Repeated pixels of each x j and y are indicated by an arrow; Collections of all of pixels from x j and y were displayed in (c-g) respectively X 2 , X 3 and X 4 , respectively. Meanwhile, the central pixels y are extracted from the marked region in Fig. 5g. Thus, we can see abundantly repeated pixels (suggestion of an arrow in Fig. 5b) are in these vectors. When the repeated pixels account for a big proportion in the region with a sufficiently large size, the sum of pixel value in each vector comes near to one another, implying that 4 . Then, one can infer that sum of weights can be nearly 1 in Eq. (9).

Shared weights in multi-contrast images
In this section, the case where the weights in one image are close to those of another contrast image will be analyzed.
Regression weights within a small region are determined mostly by the main edge direction in it. These weights are mainly estimated from similar image patches located on edges. In the sense of least square, the influence of contrast on weights regression is very limited since multiplication of a linear system of equations by a constant factor does not change its solution. For example, in Fig. 6a and b, one can see that corresponding regions in the T1 image (Fig. 6a) and T2 images (Fig. 6b) generate similar weights ( Table 4).
The mathematical analysis on weights is simplified as listed below: Weights error meets the following equation (see the derivation of Eq. (A.6) in the Additional file 1: Appendix Regression weights are estimated by continually sampling 3 × 3 patches in a 9 × 9 region, and each patch is composed of one pixel y i and its 4 neighbors x i,j (j = 1, 2, 3, 4) along diagonal directions. Consequently, the vector y = [y 1 , ⋯, y i , ⋯, y 49 ] T ∈ ℝ 49 denotes pixels in this region and the matrix X = [x 1 , ⋯, x i , ⋯, x 49 ] T ∈ ℝ 49 × 4 stands for all neighboring pixels around those pixels of y. Here, X (orX ) is the column-full-rank matrix and their generalized inversions are represented by X + andX þ , respectively. In addition, there are the vector d = ỹ − y ∈ ℝ 49 and the matrix C ¼X−X∈ℝ 49Â4 . We measure the right hand of Eq. (10) on real MRI images at different regions and observations are summarized in Fig. 7. First, most ofX þ y−Xb ð Þ 2 are very close to 0 (Fig. 7a). Besides, most of is close to 0 (Fig. 7b). Therefore, the left hand of Eq. (10) approaches to 0 in most regions, implying that b≈b . This conclusion is confirmed in Fig. 7c, showing that almost 84% ofb−b 2 lies in small values (in the range [0, 0.25]) for the tested multi-contrast MRI images.

Results and discussions
In experiments, we verify our approach on realistic T1-weighted and T2-weighted brain MRI images. 256 × 256 T1 and T2 HR images in Fig. 9 are from Philips Company. The T1 (TR = 170 ms, TE = 3.9 ms) and T2 (TR = 3000 ms, TE = 80 ms) datasets are acquired with Fast Field Echo (FFE) sequence (FOV = Fig. 6 Regression weights within local regions of T1-weighted and T2-weighted MRI images. a is the T1-weighted image; b is the T2-weighted image. Two pairs of image region of size 9 × 9 (enclosed in wireframes, marked as S1 and S2) are extracted from (a) and (b). Note: The data are acquired on a 3 T SIMENS scanner  Figure 10 and Fig. 11 are acquired at a 3 T Siemens Trio Tim MRI scanner using a turbo spin echo sequence (FOV = 230 × 187 mm 2 , slice thickness = 5.0 mm) and the matrix size of T1 (TR = 2000 ms, TE = 9.7 ms) and T2 (TR = 5000 ms, TE = 97 ms) HR images is 384 × 324.

Super-resolution experiments
Before conducting the interpolation simulation, HR images are first blurred by 3 × 3 Gaussian smooth filter with standard deviation 0.5 and then down-sampled by a factor of 2 to obtain their LR versions as listed in Fig. 8.
The LR image will be expanded as large as the HR reference by using the basic nearest neighbor interpolation.
Then these interpolated pixels will be updated using the proposed approach. The proposed method aims to recover edge details of LR brain image. We only borrow the weight from another HR contrast image if a pixel in the expanded LR image is located on an edge. In our work, a pixel is declared to be an edge pixel if the local variance within the nearest neighbors is above a given threshold (=0.0001, under the condition of intensities of images are all normalized between 0 and 1). We set the same value of the threshold in all experiments. Although, in some locations, it is not enough to satisfy the property of weights similarity, they only take a very small proportion of the total and are not processed specially in the proposed method.
The proposed approach is compared with the bicubic method, NEDI [18], and CGI [12]. The CGI method is used to guide the interpolation process by conducting directional filtering and achieves superior Signal-to-Noise Ratio (PSNR), the Structural Similarity (SSIM) [23] and the relative l 2 norm error (RLNE), are used to quantitatively measure the supper-resolution performance. The higher PSNR indicates that the reconstructed pixel value is more consistent to the original HR image and the higher SSIM implies better image structures are preserved. Also, the lower RLNE implies better consistency to the original HR image. For the proposed method, we set the region size as 9 × 9. Within each region, 3 × 3 size patches with 1pixel-width overlap between adjacent patches is set to maximally explore the statics in the local region. These are typical settings in the original NEDI method and works well for tested images. For CGI, default parameters are used in the shared source code.
First pair of images in Fig. 9 clearly show the advantage of employing the statistical information from a HR image in another contrast. Blocky artifacts in Fig. 9c are obviously generated using the classic bicubic method. The NEDI method outperforms the bicubic method since sharper edges are observed in Fig. 9d. The CGI method recovers brain boundaries in Fig. 9e much better than NEDI. Most promising edges (Fig. 9f ) are produced by the proposed approach.
For another two pairs of images acquired on a 3 T MRI scanner in Figs. 10 and 11, it can also be observed that there are many artifacts around some edges (seeing arrows) by the bicubic method. Such artifacts can be reduced by interpolation of using NEDI and CGI, and the proposed method still produces most faithful edges.
The CGI obtains higher PSNR and SSIM and lower RLNE than both NEDI and the classic bicubic. The best objective criteria are achieved by the proposed approach as listed in Table 5. These criteria are consistent to the image quality analyzed above.

Sensitivity to the Misregistration
To evaluate how the misalignment affects the accuracy of the reconstruction result, we shift reference images   along different directions (e.g., slant, anti-slant, vertical and horizontal) by a certain amount of pixels [5]. First, we compute the evaluation criteria of CGI and the proposed method using the ground truth HR image and the interpolated HR images; Second, each number in Table 6 is obtained by subtracting the evaluation criteria of the CGI from of the proposed method and is referred as "the improvement of the PSNR or SSIM or RLNE". The positive number means that the proposed method outperforms CGI method, implying better tolerance of image misregistration. From Table 6, one can see that, under 1 to 2-pixelshift, the proposed method holds advantage over CGI.
One slice of the brain image in Fig. 11 is used in simulation.

Structural distinction in T1 and T2
In MRI, T1 and T2 images have some distinct signal intensity that may cause structural distinctions appeared. For example, a structure can be visible clearly in the T2 image and is embodied too little in the T1 image ( Fig. 10a and b, arrow B), or, in turn, a structure can be visible in T1 image and is embodied too little in the T2 image ( Fig. 10a and b, arrow A). These distinct structures may be lesions or normal organisms but are not ghosts. This is normal phenomenon in MRI.
As discussed in Super-resolution experiments, we know estimated weights are nearly invariant to image contrasts. Therefore, the super-resolution still can work decently. Fig. 12 demonstrates the proposed method produces structures consistent with the ground-truth. For example, if a structure is observed on the reference but not on the ground-truth HR image, the proposed approach will not introduce the structure into the reconstruction (Fig. 12, arrow A). Other structures, which are found on the groundtruth image but not on the reference, can be recovered faithfully (Fig. 12, arrow B). These recovered structures are not reproduced correctly as well as in the ground-truth image, and appear blurrier than its vision in the ground-truth image.

Image denoising
We agree that the noise is not obviously presented in the tested brain imaging datasets. But the proposed method has the ability to suppress noise since regression weights are estimated according to the least square rule, which intrinsically has the ability to suppress noise.
To further elaborate the noise removal, the noise at common levels (1, 3 and 5% of the maximum intensity) [24,25] is added into the ground-truth image. Results of 3% noise in Fig. 13 imply that reducing the region size to 5 × 5 or increase to be larger than 9 × 9 will reduce the PSNR, SSIM and increase the reconstruction error, RLNE. Therefore, a region size of 7 × 7 or 9 × 9 is suggested to optimally suppress the noise. For other noise levels, trend curves of objective criteria are similar with Fig. 13 and come to the same conclusion.
We also comment that if the serious noise that may injury the interpolation result, noise removal before the  interpolation should be accomplished. This is beyond the scope of this work and we leave this as the future work.

Computation time
Our method is implemented with MATLAB on a personal computer with Dual-Core CPU 3.00GHz and 2GB memory. The computation time of the proposed method is very close to NEDI, and costs around 10 s.

Conclusions
An MRI image super-resolution approach is proposed to employ the statistical information retrieved from another contrast MRI image that shares similar anatomical structures. It is found that local regression weights are very similar among multi-contrast MRI images. This property is analyzed with comprehensive mathematics and experimental evidence. Experiment results demonstrate that the image quality of the low-resolution image can be truly improved if the contrast-invariant weight is borrowed from the high resolution image of another contrast. In the future, we plan to further improve the sharpness of edges and textures by utilizing sparse representation [26][27][28][29] and local geometric directions [30][31][32]. The code of this work is available at http:// www.quxiaobo.org/project/MultiContrastMRI/Toolbox _MultiContrastMRI_Superresolution.zip.

Highlights
Multi-contrast MR images share similar anatomical structures, e.g., the T1-weighted and the T2-weighted images. Regression weights are found to be similar among multi-contrast images. Comprehensive mathematics and numerical experiments are presented trying to analyze the weights-similarity property. Regression weights are learnt from another contrast high-resolution MRI image. An MRI image super-resolution approach using local regression weights is proposed. Compared with classic state-of-the-art interpolation techniques, the performance of the proposed method is remarkably improved.