Non-smooth Convex Optimization for O-Space Reconstruction
Jing Cheng1, Haifeng Wang1, Yuchou Chang2, and Dong Liang1,3

1Lauterbur Research Center for Biomedical Imaging, Shenzhen Institude of Advanced Technoleoy,Chinese Academy of Sciences, Shenzhen, China, 2Department of Computer Science and Engineering Technology, University of Houston - Downtown, Houston, TX, United States, 3Research Center for Medical AI, Shenzhen Institude of Advanced Technoleoy, Chinese Academy of Sciences, Shenzhen, China


Non-linear spatial encoding magnetic (SEM) fields can accelerate data acquisitions and improve the image quality. O-Space imaging generates a radially varying SEM field for spatial encoding in order to achieve more efficient encoding. In this work, we introduce and evaluate a novel primal dual algorithm which can handle the inverse problems of non-smooth convex optimization with non-linear forward operators to reconstruct O-Space images from undersampled data. The experimental results on simulated data show that the proposed method can achieve better image quality compared with the existing methods.


Recently, some schemes of applying SEM fields1-7 have been studied to accelerate data acquisition. In terms of speed, O-Space outperforms Cartesian SENSE when the effective acceleration factor approaches, equals, or exceeds the number of radiofrequency (RF) coils used. Due to the nonlinear gradient, O-Space data are usually reconstructed using Kaczmarz algorithm8, a L2- norm minimization method also called the algebraic reconstruction technique (ART). Compressed sensing (CS) approaches have been successfully combined with parallel imaging using conventional linear gradient encoding9-11. CS has also been applied to O-Space to reduce the mean squared error of images at high acceleration12. In this work, we applied a novel primal dual algorithm which was suitable for solving non-smooth convex optimization problems with non-linear operators to reconstruct images from highly undersampled O-Space data. Simulations demonstrated that this method could reduce undersampling artifacts and improve image quality for O-Space imaging.


In nonlinear encoding, the magnetic resonance signal from the q-th RF channel should satisfy the following equation: $$s_q(t)= \int_{\Omega} M(X)C_q(X)e^{j\Phi(X,t)}dX (1)$$ where $$$M(X)$$$ is the magnetization at location $$$X$$$;$$$C_q(X)$$$is the sensitivity of coil q; $$$\Omega$$$ is the region of interest; $$$\Phi(X,t)$$$ is the spatially dependent encoding phase, which can be decomposed as, $$\Phi(X,t)=\kappa^T(t)\Psi(X) (2)$$ where $$$\kappa(t)$$$ is the vector of gradient moments, and $$$\Psi(X)$$$ represents the gradient encoding fields at position $$$X=(x,y,z)$$$. More explicitly, for the O-Space echo corresponding to the l-th center placement (CP) of $$$(x_l,y_l)$$$, the Eq. (2) becomes $$\Phi(X,t)=\kappa_x(t)\cdot x+\kappa_y(t)\cdot y-\frac{\kappa_{z2}(t)}{2}\cdot((x-x_l)^2+(y-y_l)^2) (3)$$ where $$$\kappa_x(t)=\gamma\int_{0}^{t} G_x(\tau)d\tau$$$, $$$\kappa_y(t)=\gamma\int_{0}^{t} G_y(\tau)d\tau$$$ and $$$\kappa_{z2}(t)=\gamma\int_{0}^{t} G_{z2}(\tau)d\tau$$$; $$$G_x(t)$$$, $$$G_y(t)$$$ and $$$G_{z2}(t)$$$ are gradients waveforms on X, Y and Z2 directions; $$$\gamma$$$ is the gyromagnetic ratio. For O-Space imaging, $$$G_x(t)$$$ and $$$G_y(t)$$$ are actually the same as those of radial imaging, while $$$G_{z2}(t)$$$ describes a constant amplitude of the circumferentially invariant, radially varying nonlinear SEM. For the case of the l-th CP located at the location of $$$(x_l,y_l)$$$, the O-Space signal Eq. (1) can be written explicitly as$$s_{l,q}(t)= \int_{\Omega} M(X)C_q(X)e^{j\gamma\int_{0}^{t} \left(G_x(\tau)x+G_x(\tau)y+G_{z2}(\tau)\left[\frac{1}{2}((x-x_l)^2+(y-y_l)^2) \right]\right)d\tau}dX (4)$$ In the discrete case, assuming $$$E$$$ is the O-Space encoding matrix and with desired images $$$M\in C^{N\times N}$$$, acquired data $$$s\in C^N$$$, the Eq. (4) simply becomes a matrix equation: $$$s=EM$$$ which includes encodings from all coils and CPs. The desired image recovery problem can be written as the following convex problem,$$\widetilde{M}=argmin\left\{\lambda R(M)+||s-EM||^2_2 \right\}(5)$$ where $$$R(M)$$$ is a penalty function that incorporates the prior knowledge.

The minimization in (5) has been traditional addressed using gradient based methods such as gradient descent or its extensions to higher order derivatives. However, many regularizers of interest result in a non-differentiable objective functional, thus gradient based methods are not applicable. A common solution is to consider a smooth approximation, which however introduces additional parameters and gives non-exact results. A first-order primal dual algorithm13, introduced by Chambolle and Pock, directly works with non-smooth convex optimization problems. This primal dual algorithm (PD) has been extended to problems with non-linear forward operators14, which is suitable for O-Space reconstruction.With PD algorithm, the O-Space images can be reconstructed by repeating the following steps until a convergence criterion is satisfied: $$\begin{cases}h_{n+1}=prox_{\sigma F^*}(h_n+\sigma E(\overline{M}_n)) \\M_{n+1}=prox_{\tau G}(M_n-\tau \left[\triangledown E(M_n)\right]^*(h_{n+1})) \\\overline{M}_{n+1}=M_{n+1}+\theta (M_{n+1}-M_n)\end{cases} (6)$$ where $$$\left[\triangledown E(M_n)\right]^*$$$ is the adjoint of the derivative of $$$E$$$ in point $$$M_n$$$, $$$h$$$ is the dual variable. The proximal mapping of the convex function $$$H(x)$$$ is obtained by the following minimization:$$prox_{\sigma H}(x)=agrmin_{z}\left\{H(z)+\frac{||z-x||^2_2}{2\sigma} \right\} (7)$$

Experiments and Results

We compared the proposed PD method with ART and conjugate gradient (CG) on simulated data. Simulations were performed in MATLAB (The MathWorks, Natick, MA, USA) using one geometric numerical phantom (Resolution: 128×128; FOV: 10cm; Z2: 1.5 mT/m2) and another two brain phantoms (Resolution1: 64×64; Resolution2: 128×128; FOV1: 10cm; FOV2: 20cm; Z2: 0.15 mT/m2). Simulated receiver coil sensitivity profiles of eight element arrays were used to generate O-Space data acquisitions.

Numerical simulations are shown in Figure 1 to compare the reconstructions of undersampled O-Space data with different methods at acceleration factor of 32. As seen in Figure 1, the proposed PD algorithm reconstructs the image more accurately than other methods. And the improvement of the proposed PD algorithm is also indicated by the quantification metrics in Table 1. Figure 2 illustrates the reconstructions on the simulated brain data. The acceleration factor here is 16. The results show that the proposed PD algorithm can remove more artifacts than the previous ART or CG methods.


In this work, we reconstructed O-Space images from highly undersampled data using a proposed primal dual algorithm which is suitable for non-linear operators. Compared with other commonly used methods, the proposed algorithm shows superiority in accuracy according to the simulation experiments. In the future,more in vivo experiments using the proposed method will be studied.


This study is supported by the grant from the National Science Foundation of China (No. 61471350, No. 81729003, No. 61871373), National Science Foundation of Guangdong Province (No. 2018A0303130132), and Guangdong Provincial Key Laboratory of Medical Image Processing (No. 2017A050501026).


1. Hennig J, Welz AM, Schultz G, et al. Parallel imaging in non-bijective, curvilinear magnetic field gradients: a concept study. MAGMA. 2008; 21(1-2): 5-14.

2. Stockmann JP, Ciris PA, Galiana G, et al. O-Space images: highly efficient parallel imaging using second-order nonlinear fields as encoding gradients with no phase encoding. Magn Reson Med. 2010; 64(2): 447-56.

3. Tam LK, Stockmann JP, Galiana G, Constable RT. Null Space Imaging: nonlinear magnetic encoding fields designed complementary to receiver coil sensitivities for improved acceleration in parallel imaging. Magn Reson Med. 2012; 68(4): 1166-75.

4. Gallichan D, Cocosco C, Dewdney A, et al. Simultaneously driven linear and nonlinear spatial encoding fields in MRI. Magn Reson Med. 2011; 65(3): 702-714.

5. Galiana G, Constable RT. Single Echo MRI. PLoS One. 2014; 9(1): e86008.

6. Wang H, Tam LK, Constable RT, Galiana G. Fast rotary nonlinear spatial acquisition (FRONSAC) imaging. Magn Reson Med. 2016; 75(3):1154-65.

7. Wang H, Tam LK, Kopanoglu E, et al. Experimental O-Space Turbo Spin Echo Imaging. Magn Reson Med 2016; 75(4):1654-61.

8. Herman GT, Lent A, Rowland SW. ART: mathematics and applications, A report on the mathematical foundations and on the applicability to real data of the algebraic reconstruction techniques. J Theor Biol 1973; 42(1):1-32.

9. Lustig M, Donoho D, Pauly JM. Sparse MRI: The application of compressed sensing for rapid MR imaging. Magn Reson Med 2007; 58(6):1182–1195.

10. Liang D, Liu B, Wang J, Ying L. Accelerating SENSE using compressed sensing. Magn Reson Med 2009; 62(6):1574–1584.

11. Lustig M, Pauly JP. SPIRiT: Iterative self-consistent parallel imaging reconstruction from arbitary k-space, Magn Reson Med 2010; 64(2):457–471.

12. Tam LK, Galiana G, Stockmann JP, et al. Pseudo-random center placement O-space imaging for improved incoherence compressed sensing parallel MRI. Magn Reson Med 2014; 73(6):2212–2224.

13. Chambolle A, Pock T. A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging. J. Math. Imag. Vis. 2011; 40:120-145.

14. Valkonen T. A primal–dual hybrid gradient method for nonlinear operators with applications to MRI. Inverse Problems2014; 30(5):05501.


Figure 1. Comparison of different reconstruction results with R=32 on numerical phantom. The enclosed part is enlarged for a close-up comparison. The PD algorithm can achieve better image quality.

Table 1. The detailed comparison on performance metrics of different reconstruction methods with R=32 on the simulated numerical phantom.

Figure 2 Comparison of different reconstruction results with R=16 on the simulated brain data. The PD algorithm shows better performance in artifacts removal.

Figure 3. Reconstructions of different methods with R=16 on the brain data. The corresponding error maps and quantification metrics indicate the better performance of the proposed PD method.

Proc. Intl. Soc. Mag. Reson. Med. 27 (2019)