Magnetic Resonance Imaging (MRI) is one of the important medical imaging modalities in modern medicine for both clinical practice and medical research. An MRI system measures samples of the Fourier transform of a cross-sectional plane of interest. The measured samples in the Fourier plane typically lie on a trajectory that is determined by the physical constraints of the imaging system. In most systems, the measured samples are transformed into a set of samples on a cartesian grid via interpolation, such that the required real space image is obtained by simple inverse discrete Fourier transform. However, this type of image reconstruction method leads to serious artifacts, when the measured Fourier samples are not sufficiently dense. A dense sampling is impractical in cases where it is required to reduce the scan time. In such cases, the method of compressive sensing, pioneered by Donoho and others, becomes a better alternative, as it can reconstruct good quality images from sparse samples. Unlike direct inversion, the method of compressive sensing formulates the reconstruction problem as an inverse problem and solves it by employing numerical optimization. However, it was observed that the quality of reconstruction yielded by compressive sensing methods falls abruptly when the under-sampling and/or measurement noise goes beyond a certain threshold. Recently developed learning-based methods are believed to outperform the compressive sensing methods without a steep drop in the reconstruction quality under such imaging conditions. However, the need for training data limits their applicability. We develop a novel spatially adaptive regularization method that outperforms compressive sensing methods as well as selected learning-based methods without any need for training data in the context of reconstructing images from sparse Fourier samples. As the proposed method recovers good quality reconstruction from highly under-sampled Fourier data, it will have significant impact on speeding-up magnetic resonance imaging (MRI).

The proposed regularization method adaptively combines Eigenvalues of the Hessian and first-order TV (TV-1) functionals with spatially varying weights. This allows the functional to restore fine image structures through directional weighting, in terms of the local Eigenvalues. The proposed regularization functional R_hsa is given by

where r=[x,y] represents the 2D image grid, s(r) is the candidate image, are the spatially varying weights, represents the TV-1 term, and denote the Eigenvalues of the image Hessian at pixel location r. We propose to obtain the restored image s^* and optimal spatial weights by joint minimization of the combined cost function given by

where h is the imaging system point spread function (psf), m is the given measurement, F(s,h,m) is the data fitting term, B(s) is bounding constraint on the image pixel values, and is a weight penalty term to avoid rapid spurious changes in spatial weights for insignificant differences in the values of three terms: TV-1 and two Hessian Eigenvalues.

The cost function J_hsa is non-convex jointly with respect to image s and spatial weights, but it is convex with respect to either image or spatial weights individually. In this regard, we propose a Block Coordinate Descent (BCD) scheme that alternates between the spatially varying weight estimation and image computation. We solve the first problem of weight estimation by constructing a novel iterative method, and the second problem of image computation by deriving a novel proximal operator for regularization involving unequally weighted Hessian Eigenvalues. Since joint minimization is non-convex, we adopt a multi-resolution approach to efficiently initialize the BCD method. We call our method the Hessian Combined Order Regularization with Optimal Spatial Adaptation (H-COROSA). We experimentally compare H-COROSA with well-known regularization methods and selected learning-based methods for MRI reconstruction from under-sampled Fourier data. The comparison of reconstruction results with a recent PlugnPlay ADMM based deep learning method (referred to as PNP) [1] is shown in Figure 1.

[1] Ryu E, Liu J, Wang S, Chen X, Wang Z, Yin W. Plug-and-play methods provably converge with properly trained denoisers. In International Conference on Machine Learning 2019 May 24 (pp. 5546-5557).