In compressed sensing the optimization problem typically has a cost function of the form

where Ψ is an appropriate sparsifying (often wavelet) transform, first described for MRI by Lustig et al. (2007).

As is well known, any solution that minimizes the least squares function

satisfies the normal Equation 11 as well. We can therefore use this form to compare the approaches. Here, we include the full expression for the space-limited and band-limited discrete Fourier transform matrix A(= S_kFS_χ) to illustrate its similarity with the compressed sensing framework. It should be noted that, though masking is often applied in iterative QSM reconstructions, this is the first time that the mask-based background field removal unique to QSM is leveraged to generate conditions through which the masking turns the problem into a well-posed integral equation.

To compare QSM reconstruction performance between the compressed sensing and incomplete spectrum approaches, we added the additional sparsity-promoting regularization term λ||Ψχ||₁ to our incomplete spectrum optimization procedure (see Equation 12) in a “regularized incomplete spectrum” reconstruction. This was implemented in the Julia programming language using the “RegularizedLeastSquares” package, which is closely related to the MRI-specific work by Knopp and Grosser (2021). The sparsifying transform used in the “regularized incomplete spectrum approach” was a Daubechies wavelet with 2 vanishing moments (db2) (Vonesch et al., 2007), which is commonly used for this purpose (Majumdar and Ward, 2012).

To validate the incomplete spectrum QSM reconstruction method, we used a numerical phantom. This means that there was a known ground-truth susceptibility distribution so the QSM reconstruction error could be computed. Rather than using a reconstructed susceptibility map as a ground truth (which could lead to an “inverse crime” (Marques et al., 2021), we used the QSM reconstruction challenge 2.0 dataset from the QSM Challenge 2.0 Organization Committee et al. (2021) (see Figure 2), simulated using a comprehensive model.

The input to this method was the Sim2 dataset's local field map with signal to noise ratio SNR1. As the unwrapped, local field map was available, no additional pre-processing (i.e., background field removal) was necessary for susceptibility calculation. The brain mask used was the ground-truth mask provided with the dataset.

Since the simulated phantom is a high-resolution, high SNR data set, we also tested the performance of the incomplete spectrum method using brain images acquired in vivo. The in vivo dataset used is from Karsa et al. (2019), a gradient recalled echo acquisition at 3 Tesla with 1 mm isotropic resolution in five healthy volunteers. This was a 5 echo acquisition with echo times TE₁ = 3ms, 5.4 ms echo spacing, 20° flip angle, TR = 29 ms, and pixel bandwidth = 270 Hz. Before susceptibility calculation, these data were processed according to the pipeline described in Karsa et al. (2019), described briefly here: the total field map was calculated from the multi-echo data using non-linear complex fitting (Liu et al., 2013). The total field map was then unwrapped using Laplacian unwrapping (Schweser et al., 2013), followed by background field removal with projection onto dipole fields (PDF) (Liu et al., 2011). The brain mask for background field removal was generated by combining a mask from the FMRIB Software Library's brain extraction tool (FSL BET) (Smith, 2002) (applied to the last-echo magnitude image) with a mask obtained by thresholding the inverse noise map derived from the non-linear fitting (Liu et al., 2013) as proposed by Karsa et al. (2020).

The simulated numerical phantom dataset was used to investigate the effect of the choice of mask S_χ and band-limit S_k on the incomplete spectrum QSM reconstruction.

To illustrate the effect of the band limit on incomplete spectrum reconstruction, in Figure 1 we present a sagittal “slice” in k-space of: the k-space data (ν) input into conventional QSM algorithms, the band limit (S_k), and the k-space of the incomplete spectrum QSM reconstruction (Fχ_recon) corresponding to this band-limit. The reconstructed spectrum shows that the single streak in the input spectrum, which typically results in similar streaks in image space, is essentially split into two streaks at the boundaries of the discarded k-space, after being “filled-in” by the incomplete spectrum approach.

The incomplete spectrum QSM reconstruction was applied in vivo with the threshold value optimized on the QSM Challenge dataset. In both datasets the novel incomplete spectrum method (as well as a compressed sensing regularized version of it) were compared to four different QSM reconstruction methods chosen from different categories of susceptibility calculation algorithms:

In the numerical phantom the parameters of these reconstruction methods were tuned as follows, and the results can be found in Table 1.

The results of the band limit analysis can be found in Figures 3A, C with the sagittal slice of the PSNR-optimal and XSIM-optimal incomplete spectrum QSM reconstructions shown in Figures 3B, D for comparison. These show the effect on the reconstruction of changing the threshold t_well, which changes the size of the well-posed k-space region. The PSNR-optimal reconstruction shows slightly less pronounced streaking artifacts than the XSIM-optimal reconstruction: See for example, the orange arrow in Figures 3B, D. However, the PSNR-optimal reconstruction is smoother and has much lower contrast than the XSIM-optimal reconstruction (Figures 3B, C). This loss in contrast is highlighted when comparing the two reconstructions in the six brain regions of interest, as shown in Figure 4. Based on these results, the XSIM-optimal threshold t_well = 0.25 was used for the reconstructions of the in-vivo datasets. The results of the space-limit S_χ investigation, where the effect of erosion and dilation of the mask are analyzed, are shown in Figure 5.

The PSNR-optimal regularization weight for FANSI was 1·10⁻⁴ but this gave over-regularized reconstructions (i.e. smoothed and with loss of contrast) for the in-vivo dataset, and was therefore reduced to the XSIM-optimal weight of 1·10⁻⁵ which gave acceptable reconstructions with minimal streaking artifacts. Both regularization weights are included in Table 1 for reference. The regularization parameters giving optimal PSNR for each method in the numerical phantom are all shown in Table 1. These regularization parameters, tuned in the numerical phantom, were used to reconstruct all the in-vivo volunteer data.

The XSIM and PSNR metrics for all reconstruction methods on the challenge phantom data as well as their parameters can be found in Table 1. These metrics show that the incomplete spectrum approach is more accurate, with fewer streaking artifacts than the direct TKD approach. However, the incomplete spectrum approach performs slightly worse than the regularized iterative FANSI and CS methods. The performance of the algorithms are slightly different according to the XSIM and PSNR metrics. Adding regularization to the IS method offers minor improvement to XSIM at the cost of PSNR.

Figure 6 shows a comparison of ROI mean susceptibility values in the simulated phantom for the different QSM reconstruction methods. For reference, ground truth ROI susceptibility values are given by the horizontal red lines.

The XSIM-optimal threshold of t_well = 0.25 was used in the in vivo reconstructions, as it provided higher contrast and less smooth susceptibility maps than the PSNR optimal threshold. All six QSM reconstruction methods are visually compared in Figure 7, where we have displayed sagittal and coronal slices to emphasize QSM contrast between brain structures typically of interest as well as the level of streaking artifacts. Figure 8 shows difference maps between our proposed incomplete spectrum QSM reconstruction method and the conventional QSM reconstruction methods investigated. The difference maps show that most of the differences between the conventional QSM reconstructions and the incomplete spectrum reconstruction seem to be residual streaking susceptibility differences and not anatomical artifacts although some deep-brain gray-matter regions appear brighter in the IS QSM, particularly compared with NDI, highlighted by orange arrows in Figures 8F, L.

In Figure 9 the ROI segmentation of a representative healthy volunteer is shown. And Figure 10 shows a comparison of ROI mean susceptibility values, averaged over all five volunteers, for the different QSM reconstruction methods. For reference, literature ROI mean susceptibility values are given by the horizontal red lines.

Investigations in the numerical phantom from the QSM challenge showed that there is a frequency-space band limit S_k that is optimal for this dataset Figure 3, which is relatively robust (a 50% reduction of the threshold resulted in a 7.1 % decrease in the PSNR metric, and decreased the XSIM metric by 2.8%). This XSIM optimal threshold value of t_well = 0.25 seems high, as it fills in almost two thirds of k-space. The results, however, do not indicate over-regularization (through the loss of anatomical contrast or smoothing), and the incomplete spectrum method provided high-quality reconstructions in vivo without requiring additional tuning, further supporting its robustness.

Masking, although non-trivial (Smith, 2002), is already an integral part of most QSM pipelines i.e., for background field removal (Schweser et al., 2017). This means that finding an image support S_χ suitable for this incomplete spectrum approach is straightforward. Investigating the effect of eroding and dilating the mask in the numerical phantom showed that using a mask that was slightly larger than the brain region of interest does not negatively affect the QSM reconstruction (actually increasing PSNR in Figure 5). It should be noted, however, that a larger mask slows down the convergence of the method. The reconstruction converged more than twice as fast for 8 voxels of erosion with respect to the original mask: 7.3 seconds compared to 16.2 seconds with no erosion. The decrease in the XSIM metric compared to the PSNR metric as the mask is dilated highlights XSIM's increased sensitivity to local variance error (in this case outside of the original brain mask). Visually, neither the eroded nor the dilated reconstructions feature artifacts related to the choice of mask (see Figures 5C, D), which demonstrates the IS method to be relatively robust to masking (provided the background field removal was successful for the given mask).

Comparing the reconstructions in the numerical phantom (see Table 1 and Figure 6) we find that the incomplete spectrum performed slightly better than the direct TKD method but worse than the conventional state-of-the-art methods. Since the compressed sensing approach used such a small ℓ₁ penalty term, adding the same regularization term to the incomplete spectrum reconstruction did not improve the metrics, although there is a visual difference between the regularized and unregularized IS reconstructions in vivo (see Figures 7A, B, 8). The ROI based analysis (Figure 6) shows a slight underestimation by the incomplete spectrum approach of the mean susceptibility in the globus pallidus (GP), caudate, red nucleus, thalamus and substantia nigra, compared to conventional methods except for the compressed sensing reconstruction which consistently underestimated the mean susceptibility in these ROIs. Overall, NDI performed the best on the challenge phantom as it provides ROI mean susceptibility values closest to the ground truth values except in the substantia nigra (Figure 4). This is not reflected in the XSIM and PSNR metrics (Table 1) nor in the in-vivo dataset ROIs (Figure 6), where NDI generally gave the lowest and least accurate reconstructed mean susceptibilities of the algorithms compared.

The “doubling” of the streaking artifact as observed, for example, around the calcification in the top of the brain in the simulated dataset, seems to be a side-effect of the way the incomplete spectrum method reconstructs the missing frequency domain information as illustrated in Figure 1. It is not as apparent in vivo, most likely because this dataset has a lower overall SNR. Since this IS method fills in k-space information between cones at two angles centered on the origin it replaces more high frequency than low frequency information, leading to a minor denoizing effect (observable as smoothing, specifically when comparing the in vivo TKD and incomplete spectrum reconstructions).

Since the incomplete spectrum method presented here does not rescale regions of the frequency domain (like TKD), it does not directly change the point spread function of the operator, and no rescaling of the output susceptibility map should be necessary [as was proposed for TKD by Schweser et al. (2013)]. This is supported by the incomplete spectrum results on the simulated dataset where no scaling in contrast was observed for higher levels of (intrinsic) regularization (although a slight loss in contrast due to smoothing was observed, as in the comparison between XSIM-optimal and PSNR-optimal reconstructions shown in Figures 3, 4).

Figure 10 shows that, in vivo, the incomplete spectrum approach gave ROI mean susceptibility values similar to those from commonly used, state-of-the-art algorithms i.e., NDI and nlTV. This suggests that the incomplete spectrum method did not systematically overestimate or underestimate the susceptibilities in these ROIs. In most of the ROIs, the incomplete spectrum approach gave susceptibility values between those from nlTV and TKD. The compressed sensing reconstruction gave the lowest susceptibility values in all ROIs. The literature values were larger than the values reconstructed by all algorithms in the caudate, putamen, substantia nigra (and globus pallidus) perhaps because the healthy volunteers were relatively young compared to the subjects included in the studies cited here and may, therefore, have had a lower iron content and a lower susceptibility in these deep-brain gray matter ROIs than the subjects in the cited studies (Li et al., 2014, 2023; Zhang et al., 2018).

With additional ℓ₁ wavelet regularization, the incomplete spectrum QSM reconstruction more closely resembles the compressed sensing QSM (Wu et al., 2012) (Figures 7A–C, G–I, 8B, C, H, I). At high levels of regularization the regularized incomplete spectrum and compressed sensing reconstructions become identical, as expected from their similar cost functions (see Equations 12, 13), and illustrated in Supplementary Figure S1. Adding the regularization to the incomplete spectrum approach leads to less pronounced streaking (as illustrated by the orange arrow, Figures 7A, B). However, regularization adds another parameter to the method which then requires tuning, as opposed to the original implementation which is less sensitive to t_well initially tuned on the numerical phantom. Further, it has been shown that band-limiting based on thresholding the dipole kernel as performed in this incomplete spectrum approach leads to correlated artifacts (Wu et al., 2012) thereby violating the assumptions of compressed sensing (i.e. uncorrelated artifacts). Therefore, additional regularization using this compressed sensing term is disadvantageous compared to the original unregularized incomplete spectrum approach.

The biggest limitation of the incomplete spectrum QSM reconstruction approach lies in the input data ν = D⁻¹Fb, which is the directly deconvolved phase information. Streaking artifacts are introduced by multiplication with the inverse dipole kernel D⁻¹ and the incomplete spectrum method then “corrects” for these. Future work will involve applying the incomplete spectrum method to reconstructed susceptibility maps, to provide additional “correction” of the spectrum where applicable. This could lead to frequency domain correction schemes that would work in tandem with conventional iterative methods to improve the robustness to noise, or artifacts that can be identified in specific regions of frequency domain (such as streaking). This technique provides a new tool for filling in missing regions of frequency space, or correcting “ill-posed” or noisy frequency-space regions with a suitable band-limit in QSM maps by using an image space mask.