Computational and Mathematical Methods in Medicine

Volume 2017, Article ID 6462832, 12 pages

https://doi.org/10.1155/2017/6462832

## Second-Order Regression-Based MR Image Upsampling

Department of Computer Science, Chengdu University of Information Technology, Chengdu 610225, China

Correspondence should be addressed to Jiliu Zhou; nc.ude.tiuc@uilijuohz

Received 4 January 2017; Revised 9 March 2017; Accepted 15 March 2017; Published 30 March 2017

Academic Editor: Giancarlo Ferrigno

Copyright © 2017 Jing Hu 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 magnetic resonance imaging (MRI) is often limited due to several reasons, including a short data acquisition time. Several advanced interpolation-based image upsampling algorithms have been developed to increase the resolution of MR images. These methods estimate the voxel intensity in a high-resolution (HR) image by a weighted combination of voxels in the original low-resolution (LR) MR image. As these methods fall into the zero-order point estimation framework, they only include a local constant approximation of the image voxel and hence cannot fully represent the underlying image structure(s). To this end, we extend the existing zero-order point estimation to higher orders of regression, allowing us to approximate a mapping function between local LR-HR image patches by a polynomial function. Extensive experiments on open-access MR image datasets and actual clinical MR images demonstrate that our algorithm can maintain sharp edges and preserve fine details, while the current state-of-the-art algorithms remain prone to some visual artifacts such as blurring and staircasing artifacts.

#### 1. Background

Magnetic resonance (MR) imaging is widely used to assess brain diseases, spinal disorders, cardiac function, and musculoskeletal injuries. Compared with computed tomography, MR imaging requires a longer acquisition time [1]. Hence, in order to minimize involuntary patient motion in MR imaging, the scan time is often shortened, thereby obtaining MR images with fewer slices and larger spacing [2]. In other words, MR images are usually highly anisotropic (e.g., 1 × 1 × 6 mm^{3}), with a lower resolution in the slice-selection direction than in the in-plane direction [3]. However, in many medical applications, an isotropic MR image is required [4]. In addition, a higher resolution of MR image is also essential for more detailed understanding of human anatomy, facilitating early detection of abnormalities, and improving clinical assessment accuracy [4].

A possible approach to increase MR image resolution is an interpolation-based image upsampling [2]. Traditional interpolation methods adopted in natural image processing field such as spline interpolation can be directly employed. But these methods use fixed interpolation coefficients and only select spatially nearby sampling voxels, thereby producing images with blurred edges and staircasing artifacts. To reduce these unwanted artifacts, some sophisticated interpolation methods [2, 5–11] have been recently proposed in biomedical image processing. Largely, these advanced interpolation methods are derived from a nonlocal redundancy concept [12]. That is, they specifically select sampled voxels or adapt interpolation coefficients. For instance, Manjón et al. [5] proposed determining the interpolation coefficients through the similarity of intensities between two 3D image patches around the unknown voxel and sampled voxels. Observing that image structures in a slice-selection direction of a 3D MR image also exist in an in-plane direction, while the resolution of the latter is usually higher than the former’s, Plenge et al. [6] proposed reconstructing a high-resolution (HR) version in a slice-selection direction by leveraging cross-scale self-similarity and cross-scale resolution discrepancy from patches in in-plane direction. Furthermore, Qu et al. [7] extended Plenge’s work by using the sparsity of similar image patches. Besides the methods proposed in [5–7] which use only the input low-resolution (LR) image, several other algorithms have also been advocated to leverage HR images of the same subjects in other imaging modalities, with the aim of reconstructing more high-frequency image details. For example, a combined interpolation weight proposed in [8] was calculated from both an example HR image and an LR input image. Further, Manjón et al. [9, 10] developed a correlation map which shows the inconsistencies between HR image and LR input so as to adjust the relative importance of these two images in the combinational weight computation. On noticing that the voxel similarity comparison via Euclidean distance of intensity was not rotationally invariant, Jafari-Khouzani et al. [2, 11] proposed to concatenate several image features computed from each voxel, such as gradients and Gaussian-blurred intensity values, into a vector and measure the voxel similarity from this feature vector therein.

Intuitively, whether using a different imaging modality or not, the above-mentioned advanced interpolation algorithms all focus on the refinement of interpolation weights. Nevertheless, using a weighted averaging scheme, these algorithms essentially can be described as a zero-order regression problem [13]. Unfortunately, the zero-order regression prototype risks blur the laminar shape of brain structure, and therefore high-order image details cannot be well retrieved [2, 13]. To the best of our knowledge, few efforts have been made in MR image processing to promote an interpolation algorithm that preserves higher-order details.

To address the fine-detail limitation in current interpolation algorithms, we propose a regression-inspired interpolation method using second-order polynomials. Under a kernel regression framework, we attempted to transform an LR image space into the expected HR image space by establishing a set of high-order prototypes so as to locally approximate the ground truth image space. Moreover, to further strengthen local consistency in the proposed method, a patch-based image reconstruction scheme is utilized rather than voxel-by-voxel scheme as other interpolation methods do. With these two salient features, the proposed method successfully enhances high-frequency details in final results.

In the rest of this paper, the proposed method will be specified in Section 2. Extensive experimental results and discussions are demonstrated in Section 3. Finally, a conclusion is described in Section 4.

#### 2. Methods

In this section, we firstly describe a MR image acquisition model and mathematical notations therein and then present the specific implementation of our regression-based image upsampling algorithm.

##### 2.1. The Acquisition Model and Notations

In the context of MR imaging, the most accepted image acquisition model assumes that an acquired LR image is generated from an HR counterpart through a sequence of degradation operations, such as motion blur, field inhomogeneity, and noise, and can be represented as [1]where denotes an observed LR image, is a downsampling operator, is a blurring operator, is an HR image, and is acquisition noise, which is often regarded as Rician-distributed in MR imaging. To infer from , image upsampling is basically an inverse process of MR imaging. Moreover, owing to the fact that the specific form of the blurring operation is unknown, this blurring operation could be regarded as an unspecified mapping function. In light of this, the HR image could then be formulated as where is a simple interpolation operator (i.e., bicubic interpolation) which generates the initial estimation of the upsampled image and is the denoised version of the LR image ; the mapping function denotes the inversed blurring operation, representing the relations from an LR space to the target HR space.

Equation (2) translates the image upsampling problem into a regression problem, with the particular form of regression function remaining unspecified. To estimate , we propose the development of a generic local expansion of about one image patch in image. This particular method will be discussed in Section 2.2.

As shown in Figure 1, the goal of the proposed method is to construct an HR version from a given LR image (which is in fact image in (2)). is a smoothed version of (where represents smoothing) and is an upscaled version of . The bolded lowercase and denote the column vectors of two image patches which are extracted from and , respectively; and are the column vectors of two image patches taken from and , respectively. Among all patches in image , is the most similar one to ; patches and have the same coordinates for the center pixel. In the proposed method, constitutes LR-HR training patch pairs and constitutes LR-HR testing patch pairs.