Introduction

The choice of the optimal parameters is a major problem of optimization based approaches to solve inverse problems such as the local phase to susceptibility deconvolution for QSM. However, in real clinical scenarios these parameters have always to be selected automatically to avoid time consuming and unfended trials. In this work noise is assumed to be a random vector and the method of unbiased predictive risk estimator (UPRE) is used to select suitable parameters. The proposed algorithm is tested on the simulated multi-echo data provided at the 2019 QSM Reconstruction Challenge (RC2)¹. This work is a further development of the algorithm presented at the ISMRM Meeting 2021² and its purpose is to show that the method of UPRE can be applied advantageously to a shearlet /TGV based susceptibility reconstruction.

Theory

The inverse problem, typically solved to obtain the susceptibility map $$$\chi$$$ from the local phase $$$\psi$$$, can be formulated in Fourier space as : $$ \mathcal{F}(\psi) (\textbf{k}) = D(\textbf{k}) \mathcal{F}(\chi) (\textbf{k}) , \quad \textbf{k} = (k_x,k_y,k_z) \in \mathbb{R}^3. $$ Here $$$\mathcal{F}$$$ denotes the Fourier transform, and $$$D$$$ the dipole kernel: $$ D(\textbf{k}) := \frac{1}{3}- \frac{k_z^2}{|\textbf{k}|^2}. $$ To achieve a unique solution without artifacts, the shearlet transform³ of the to be determined $$$\chi$$$ and a TGV-term are used as regularizers². The method of alternating direction of multipliers (ADMM) is applied to solve the resulting optimization problem containing two data-fit terms and the two regularization terms^2,4. Therefore, auxillary variables $$$w,d,t,b^w,b^d,b^t$$$ and three new parameters $$$\mu =(\mu_1,\mu_2,\mu_3) > 0$$$ are introduced and the problem is split into multiple subproblems which are easier to solve^2,5,6,7 :$$ \begin{split} \chi_{k+1} &= \text{arg}\min_ {\chi}\:\frac{\beta_1}{2} ||E\chi-\psi||_2^2 + \frac{\beta_2}{2} || \chi_{well}- Mx||_2^2 + \frac{\mu_1}{2} ||w_k - \Psi \chi – b^w_k||_2^2 + \frac{\mu_2}{2} ||d_k – (\nabla^f \chi – v) – b^d_k||_2^ +\frac{\mu_3}{2} ||t_k -\zeta^bv – b^t_k||_2^2 2, \quad (1) \\ w_{k+1} + &= \text{arg}\min_w \: ||w||_1 + \frac{\mu_1}{2} ||w - \Psi \chi_{k+1} – b^w_k||_2^2,\\ d_{k+1} + &= \text{arg}\min_d \: ||d||_1 + \frac{\mu_2}{2} ||d – (\nabla^f \chi_{k+1} – v_{k+1}) – b^d_k) – b^d_k||_2^2,\\ t_{k+1} + &=\text{arg}\min_t: ||t||_1 + \frac{\mu_3}{2} ||t – \zeta^bv_{k+1} – b^t_k ||_2^2 ,\\ b^w_{k+1} &= b^w_k + \Psi \chi_{k+1} - w_{k+1}, \\ b^d_{k+1} &= b^d_k + (\nabla^f \chi_{k+1} - v_{k+1} )- d_{k+1}, \\ b^t_{k+1} &= b^t_k + \zeta v_{k+1} - t_{k+1}, \end{split}. $$ To determine suitable parameters $$$\alpha = \beta_{1,2},\mu_{2,3}$$$ one by one (whereby the other unknown parameters are set to zero and $$$\mu_1$$$ is set as an upper limit of the shearlet coefficients of a LSQR solution $$$\chi_{LSQR}$$$), the predictive risk $$$\frac{1}{n}||p_\alpha||_2^2 = \frac{1}{n}||K\chi_{\alpha}-K \chi_{true})||_2^2$$$ with the true susceptibility map $$$\chi_{true}$$$ and the regularized solution $$$\chi_{\alpha}$$$ is estimated⁸ and paramters were chosen to minimize the unbiased predictive risk estimator $$$U$$$.
Here $$$K = \mathcal{F}^*\tilde{D}(\textbf{k}) \mathcal{F}$$$ and $$$\tilde{D}$$$ is the pointspread function of the respective data term, e.g. $$$\tilde{D} = D$$$ for the determination of $$$\beta_1$$$.

Methods

The RC2 provided simulated local phase data related to two different susceptibility ground truths (Sim1/2) and different error metrics¹. Phase data from different echoes were combined after weighting with the inverse standard deviation of the noise of the phase data and a laplacian based suppression of phase errors for each echo. Noise was assumed to be a random vector such that related quantities like $$$\psi$$$ and $$$\chi_{\alpha}$$$ and the predictive risk were also random vectors or variables. As equation (1) can be solved by diagonalization in Fourier space the regularized solution $$$\chi_{\alpha} = R_{\alpha} d$$$, where $$$R_{\alpha}$$$ is a linear regularization operator and $$$d = \psi$$$ for all but the determination of $$$\beta_2$$$ - in that case $$$d = \chi_{well}$$$^2,4. With $$$A_{\alpha}=KR_{\alpha}$$$ the predictive error can be split into a deterministic and a stochastic term: $$p_{\alpha}=(A_\alpha-I)K\chi_{true} + A_{\alpha}\eta.$$ Here $$$\eta$$$ is the noise vector. Since $$$A_{\alpha}$$$ is selfadjoint and noise is approximately Gaussian distributed (for high SNR in the magnitude image), neglecting covariance terms in the case of determining $$$\beta_2$$$ and by approximating the voxelwise variance $$$\sigma^2$$$ by the mean of the VOI, the UPRE $$$U$$$ can be expressed via the regularized residuum $$$r_{\alpha}$$$⁸: $$U(\alpha) = \frac{1}{n}||r_{\alpha}||_2^2 + \frac{2 \text{mean}(\sigma^2)}{n}\text{trace}(A_{\alpha})-\text{mean}(\sigma^2).$$ The suitable regularization parameter is then obtained by finding the roots of the derivative of $$$U$$$ with the Matlab function fzero. To determine $$$\mu_{2,3}$$$ only the data term with the local phase $$$\psi$$$ was considered. Finally, by solving $$ \tilde{s} = \text{arg}\min_s U_1(s \beta_1 )+U_2((1-s)\beta_2)$$ the scale factor $$$\tilde{s}$$$ was determined and parameters were updated to:
$$\begin{align} \tilde{\beta}_1 &= \tilde{s} \beta_1,\\ \tilde{\beta}_2 &= (1-\tilde{s}) \beta_2\\ \tilde{\mu}_{2,3} &= \tilde{s} \mu_{2,3}. \end{align}$$.

Results and Discussion

Table 1 shows the achieved image quality measures (IQM) for the susceptibility maps of the proposed algorithm (PA) and the winner algorithms of RC2 in the Robustness/RMSE category (mcTFI⁹ /wL1L2TV¹⁰).The PA achieves a top ranking except for the streaking error. The just mentioned susceptibility maps (related to the susceptibility model Sim2) are shown in Figure 1. The good IQM achieved by the PA and the visual appearance of the obtained QSM verify that the UPRE method is suitable for automatic parameter selection. The purpose of this work is to show that the UPRE method can be applied within a susceptibility map reconstruction algorithm, using the benefits of shearlets and TGV, to find suitable regularization parameters.

Acknowledgements

This work was supported by a grant from Deutsche Forschungsgemeinschaft (Grant number: DFG STR 1480/2-1). The use of the data from the 2019 Reconstruction Challenge is kindly acknowledged. The use of code for the shearlet reweighing provided by Dr. Jackie Ma is kindly acknowledged as well as the use of the code for the shearlet transform from ShearLab3D.