Researchers interested in fast and accurate reconstruction of the susceptibility distributionsAccelerating the selection of optimal regularization weights, for L1-regularized, magnitude-weighted QSM

Quantitative Susceptibility Mapping reconstructions may benefit from L1-regularization and magnitude weighing [1], however these iterative reconstruction methods are time-consuming. Recently Split Bregman iterations have been proposed[2] for L1-regularized, magnitude-weighted QSM, minimizing:

$$ \frac{1}{2}||F^{-1}DF\chi-\phi||_2^2 + \lambda_1||WG\chi||_1 $$

with F the Fourier operators, $$$ D = \frac{1}{3}-\left(\frac{k_z}{k}\right)^2 $$$ the dipole function φ is the normalized field shift, W is a mask on edges in the magnitude image G the gradient operator in three dimensions.

Prior to the final reconstruction, the optimal regularization weight λ1 balancing data-consistency and smoothness should be selected (fig1). To this end, a series of reconstructions with different regularization weights are performed, of which only the norms of the residual and the gradient of the reconstruction are used to create an L-curve [3], in which the points with the highest curvature corresponds to the optimal regularization parameter. Only the final reconstruction is carried out with magnitude weighting. For fast convergence of the L1 regularized reconstructions, also the optimal value λ2 for the L2 regularized problem should be determined with an additional L-curve: The general framework was followed as described in [1] and illustrated in fig 1. It was investigated 1) if an anti-aliasing buffer is needed in the construction of the L-curves, as it was hypothesized that both coordinates of the L-curve, the norm of the data-consistency and the regularization norm are almost completely determined by the voxels in the Region of interest. 2) If the complete L2-curve could be constructed in the frequency domain (Parseval’s theorem) based on only one FFT, and using only half the frequency components of the phase data (conjugate symmetry). 3) if it was possible for L1-curve to reconstruct two χ distributions with different regularization weights simultaneously for the cost of one, in the imaginary component of the phase distribution. This is justified by the purely real character of the field, the dipole and the susceptibility distribution, and the linearity of the Fourier and the gradient operation. Only for the single non-linear operation of the each iteration (soft thresholding) the datasets were temporarily separated 4) if for the Split Bregman iterations, gradient operations in the spatial domain were faster than those in the frequency domain.1) The point of maximal curvature in the L-curve was found to be invariant to the buffer size, making a buffer of four voxels sufficient, reducing calculation times almost twofold. Without any accuracy loss, L-curves constructions were further accelerated 2) twenty fold for the L2-curve using the single FFT approach and 3) almost twofold for the L1-curve by processing the additional reconstruction in the imaginary channel. 4) Switching the order of the Fourier and gradient operations gave 30% reduction in the final magnitude weighted reconstruction without accuracy loss. The total processing time was reduced more than threefold, without accuracy loss.Efficient reconstruction algorithms make QSM more suitable for (real-time) clinical usage, and may allow the incorporation of more sophisticated/time-consuming regularization tools. The here proposed accelerations serve these goals, yielding numerically equivalent results while reducing the total computational load more than threefold. As mentioned before [1] the efficiency of the L1-sweep can be further optimized by parallel computing.The complete pipeline of optimally L1-regularized QSM with magnitude-weighting can be accelerated more than threefold without accuracy loss.