Jing Cheng^{1}, Haifeng Wang^{1}, Yuchou Chang^{2}, and Dong Liang^{1,3}

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.

Introduction

Recently, some schemes of applying SEM fieldsMethods

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 algorithm^{13}, 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 operators^{14}, 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.

Conclusion

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.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.