Computational and Mathematical Methods in Medicine

Volume 2016, Article ID 3647202, 9 pages

http://dx.doi.org/10.1155/2016/3647202

## Diffusion-Weighted Images Superresolution Using High-Order SVD

^{1}Department of Computer Science, Chengdu University of Information Technology, Chengdu 610225, China^{2}College of Electronic Engineering, Sichuan University, Chengdu 610065, China^{3}Department of Computer Science, Xihua University, Chengdu 610039, China

Received 27 January 2016; Revised 9 July 2016; Accepted 28 July 2016

Academic Editor: Po-Hsiang Tsui

Copyright © 2016 Xi Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

The spatial resolution of diffusion-weighted imaging (DWI) is limited by several physical and clinical considerations, such as practical scanning times. Interpolation methods, which are widely used to enhance resolution, often result in blurred edges. Advanced superresolution scanning acquires images with specific protocols and long acquisition times. In this paper, we propose a novel single image superresolution (SR) method which introduces high-order SVD (HOSVD) to regularize the patch-based SR framework on DWI datasets. The proposed method was implemented on an adaptive basis which ensured a more accurate reconstruction of high-resolution DWI datasets. Meanwhile, the intrinsic dimensional decreasing property of HOSVD is also beneficial for reducing the computational burden. Experimental results from both synthetic and real DWI datasets demonstrate that the proposed method enhances the details in reconstructed high-resolution DWI datasets and outperforms conventional techniques such as interpolation methods and nonlocal upsampling.

#### 1. Introduction

Diffusion-weighted imaging (DWI) is a noninvasive magnetic resonance modality which can be used to infer features of local tissue anatomy, composition, and microstructure from water displacement measurements [1]. Water does not diffuse equally throughout the brain and this property has been applied widely for* in vivo* analysis of white matter architecture and neuronal diseases [2, 3]. Despite the rapid development of this technology and the broad application of anisotropic diffusion properties, DWI exhibits an inherently low signal-to-noise ratio (SNR) compared with other imaging modalities. In addition, because DWI implements EPI scanning in multiple directions, the spatial resolution is relatively poor under limited physical and clinical considerations such as durable scanning time of patients.

It has been shown that the limited resolution of DWI introduces a partial volume effect (PVE) which results in bias during DWI imaging analysis [4]. The improvement of DWI spatial resolution with high SNR provides a better sensitivity for the analysis of brain structure and clinical disease [5, 6]. Moreover, high-resolution DWI could improve the estimation accuracy of diffusion tensor imaging, thus proving to be beneficial to fiber tractography and finer bundle analysis [7].

Several methods have been proposed in the literature to enhance the spatial resolution of DWI. During the acquisition stage, long acquisition times remain a primary obstacle preventing the method from being of real interest clinically. For example, Miller et al. [8] implemented five days of scanning to obtain postmortem high-resolution DWI with high SNR. In order to avoid the long scanning times, superresolution (SR) acquisition emerged as an effective technology, initially proposed for MRI, but soon adapted to DWI. Subpixel shifting in the in-plane dimension was proposed to obtain multiple low-resolution images for reconstruction into higher resolution images [9]. Anisotropic scanning was another strategy used to obtain low-resolution images for reconstruction. Scherrer et al. [10] employed a maximum of a posteriori estimation from anisotropic orthogonal acquisition to reconstruct isotropic high-resolution DWI images.

Compared with the SR acquisition implemented in specific scanning protocols, the SR algorithm in the postprocessing stage has been involved with the scene image SR reconstruction. This method category is independent of acquisition protocol and was previously implemented on MRI. The most intuitive methods for increasing resolution are interpolation methods, such as bicubic and B-spline interpolation [11]. These methods estimate the new voxels according to some smoothness assumptions which are not valid in inhomogeneous areas. Hence, interpolation methods usually result in blurred edges and artifacts in lines. Nowadays, advanced superresolution algorithms for scene images have been proposed for MRI and the main idea is reconstructing high-resolution information from the image content. Manjón et al. [12] used a nonlocal estimator to reconstruct high-resolution MRI using a single low-resolution dataset. Rousseau [13] involved multimodality MRI to improve the SR quality. Coupé et al. [7] implemented a nonlocal estimator in DWI and also incorporated information to enhance the reconstruction results. Sparse representation, a recent trend in signal and image processing, has also been effectively implemented in MRI. Rueda et al. [14] reconstructed single MRI datasets based on prior knowledge with a pretrained overcomplete dictionary. Trinh et al. [15] extended the sparse representation with nonnegative expressions in order to remove noise and superresolve the two sets together.

Most of the currently used SR methods are based on the structure MRI modality and, when applied to DWI, joint information should be considered for better reconstruction [16]. Joint information comes from the redundancy acquired from adjacent scanning directions. This extra information is beneficial for enhancing the spatial resolution of DWI. Single value decomposition (SVD) plays a central role in reducing high-dimensional data to lower-dimensional data and is a classical method involved in inverse problems such as denoising [17] and restoration [18]. Recently, high-order single value decomposition (HOSVD) was used to generalize the SVD of a matrix into a high-order matrix and offered a simple yet elegant method for handling similar patches [19]. Additionally, the HOSVD basis was adapted from image content and may achieve a more sparse representation than the fixed basis. In this paper, we propose a novel SR method for DWI datasets using HOSVD. Similar to the nonlocal patch-based SR approaches successfully implemented on both MRI and DWI [7, 12, 13], HOSVD was used to construct the regularization framework in the proposed method. The merit of the HOSVD SR method stems from the adaptive HOSVD basis which results in a more accurate reconstruction. Also, the HOSVD is only implemented over similar patches in a stack, which effectively decreases the computational complexity simultaneously. This is especially useful for DWI, since the involvement of joint information from adjacent directions in DWI datasets dramatically increases the computation burden [16].

The remainder of this paper is organized as follows: we first describe the proposed method in detail and then apply it to both synthetic and* in vivo* DWI datasets for experimental evaluation. Experimental results and computational efficiency are demonstrated in Section 4 and concluding remarks are given in Section 5.

#### 2. Methods

Image SR leads to an ill-posed inverse problem which is related to the LR image and HR image ; the general model can be expressed as follows:where represents acquisition noise, represents the decimator operator, and represents the degradation function [7, 12, 13].

Based on this model, the SR image can be estimated by minimizing a least-square cost function as follows:

For such inverse problems, a regularization term should be added to stabilize the convergence; thus, the HR image can be estimated from the LR observation using the following equation:where is the regularization term, is a fidelity term, and is a balancing parameter. As demonstrated by Coupé et al. [7], nonlocal patches methods can be an efficient way to define the regularization term. Instead of using a nonlocal mean estimator, we propose the implementation of a high-order SVD to be used as the estimator in this study, owing to its simple application and promising performance [19].

The HOSVD estimator clusters similar patches into a stack, in a manner similar to other patch-based methods [7, 12], and then performs an HOSVD transformation to obtain the HOSVD basis and coefficients. After the truncation of the coefficients, the patches are then reconstructed by an inverse HOSVD transform.

With this in mind, the regularization term for the superresolution process in (3) can be defined as follows:where is the HOSVD base estimator.

Given an patch centered in , we define such similar patches (including ) as , where , and the similar patches are obtained as follows.

Let us denote as the stack ; the HOSVD of the stack can then be defined as [20]where is the set of coefficient matrices for a three-order tensor with , stands for the th mode tensor product defined in [20], and , , and are orthonormal unitary matrices.

After applying the HOSVD transform, the patches can be estimated by nullifying the coefficients under the assumption that the coefficients of the clean image have a sparse distribution. As indicated in [19], the coefficients can be truncated using hard thresholding as follows:where denotes the hard threshold defined by for the stack with patches of size . As noted in [20], the coefficients in tensor** S** are not necessarily positive and the hard thresholding is defined as the absolute value of the coefficient array:where denotes the th element of tensor .

After truncation, the stack is reconstructed using an inverting transform with truncated coefficients to obtain the final HOSVD estimator :

Since the DWI datasets are three-dimensional, the above method should be extended to include fourth-order HOSVD transforms for stacks with 3D similarity patches. Besides this, the threshold should be modified as for a stack of size .

In [12], a mean consistency correction followed the estimator to ensure coherence with the physical acquisition model. This was implemented in the fidelity term:

This was done for the entirety of location in the LR dataset while subsampling consistency was imposed on the reconstructed patches. Finally, the iteration process is summarized by (8) and (9) and is applied until convergence:where NN is the nearest neighbor interpolation and is the iteration number.

In order to further improve the SR performance, the proposed HOSVD SR method can be augmented using joint information from adjacent directions in the DWI dataset [16]. For each patch , the corresponding stack was constructed with similar patches which were determined as follows: the distance threshold selected all patches for which was chosen as , where is the variance of the noise. This threshold is balanced between the estimation accuracy and the computational speed as indicated in [19]. The joint information was introduced by enlarging the search window into the adjacent DWI datasets, where is defined as and denotes the directions before and after it. In this paper, the HOSVD superresolution method which uses joint information in multiple directions is referred to as HOSVD-M.

#### 3. Experiments

In order to quantitatively evaluate the quality of the reconstruction, B-spline interpolation, which has been introduced for DWI resolution enhancement in the literature [21, 22], is used for comparison. In addition to this, a nonlocal approach for image SR [12] is also involved as an effective nonlocal patch-based SR method for comparison purposes. In this section, both synthetic and* in vivo* datasets were implemented for evaluation. The patches size was empirically set to 5 as suggested in [23], which was for denoising purposes primarily and also demonstrated robust results in this work. The balance parameter was set to 0.01 in all experiments. Since a sensitivity analysis for this parameter showed that the values between 0.001 and 0.2 only generated less than 0.1 dB variations of the PSNR, this means that the reconstruction has little dependence of this parameter which was also observed in the literature [14].

The simulated dataset without noise was chosen as ground truth, which consists of the 3D structure field presented at the 2012 HARDI Reconstruction Challenge [24] and occupies a 16 × 16 × 5 volume, mimicking a realistic 3D tract configuration. As shown in Figure 1(a), this dataset is comprised of five different fiber bundles which give rise to the nonplanar configurations of bending, crossing, and kissing tracts. All fiber tracts were characterized with a fractional anisotropy between 0.75 and 0.90. To better explore the proposed method, this synthetic dataset was also corrupted by Rician noise (SNR = 30) as demonstrated in Figure 1(e). Both the original dataset and the noisy set were downsampled by factor 2 using nearest neighbor interpolation along each axis. Afterwards, the LR datasets were superresolved using the B-spline method, the nonlocal method, and the proposed method, respectively. In addition to the visual comparison demonstrated in Figure 1, the angular accuracy was also measured for quantitative evaluation purposes [24]. The angular accuracy in the orientation of the estimated fiber compartments was assessed by medians of the average error (in degrees) between the estimated fiber direction and the true direction present in each voxel:where the unitary vectors and are a true fiber population in the voxel and the closest of the estimated directions. For further analysis, the average error in all voxels was calculated and demonstrated in box-and-whisker diagrams. The upper and lower edges of the boxes are the 75 and 25 percentile, respectively; the smallest and biggest observations are the two ends of the whisker. The mean and median are demonstrated using red dot and line, and for each reconstruction dataset, 2% of the worst results were selected as outliers to eliminate the anomaly results.