Weakly-supervised Learning for Single-step Quantitative Susceptibility Mapping

Quantitative susceptibility mapping (QSM) can estimate tissue magnetic susceptibility values from MRI Larmor frequency sensitive phase images to provide novel image contrast [28]. To date, all QSM methods rely on a dipolar convolution that relates susceptibility sources to induced Larmor frequency offsets [21, 18], which is expressed in the k-space as bellow.

where $B(\vec{k})$ is the induced magnetic perturbation along the main magnetic field $B_{0}$ direction, $X(\vec{k})$ is the susceptibility distribution in the k space, $D(\vec{k})$ is the dipole kernel.

The generation of QSM requires solving two challenging ill-posed problems - (1) removal of background field contributions from sources outside the interest, (2) field-to-source inversion by solving the dipole deconvolution. Though existing BFR algorithms demonstrate excellent performance, they have several limitations, including volume erosion, inaccurate BFR close to volume boundary, and residual background leakage [23]. Incorrect BFR often introduces erroneous local field outputs and subsequently affects susceptibility quantification.

The field-to-source inversion faces technical challenges due to the singularities of the dipole kernel. Calculation of susceptibility through multiple orientation sampling (COSMOS) [15] remains the empirical gold-standard of QSM, as the additional field data sufficiently improves the conditioning of this ill-posed inversion. However, multi-orientation data acquisition is time consuming and clinically infeasible. Single-orientation QSM is preferred which is typically computed by either thresholding of the convolution operator [24, 31] or sophisticated regularization methods [6, 20, 14, 1]. In addition, several single-step QSM methods [4, 17, 26] have been proposed to directly estimate QSM from the total field (combined BFR and dipole inversion) to prevent potential error propagation across successive operations. However, existing QSM techniques still struggle to generate reliable QSM estimates in clinical contexts, which greatly hinders QSM clinical translation.

Recently, several deep learning (DL) QSM techniques have been proposed. For dipole inversion, QSMnet [33] and QSMGAN [5] utilized COSMOS estimates as QSM labels, while DeepQSM [2] was trained using purely synthetic data. AutoQSM [29] utilized STAR-QSM [30] estimates after SMV method [32] for BFR as QSM labels for single-step QSM learning. Though these techniques demonstrate promising performance, they have several limitations. Due to the lack of a ground-truth reference, these methods usually used COSMOS maps or synthetic data for training. However, acquiring large amounts of COSMOS data is not only expensive but also time consuming. In addition, COSMOS neglects tissue susceptibility anisotropy [11] and contains errors from BFR and image registration procedures, which compromise its value as a training label. Though synthetic data generated from the physical model provides a reliable and cost-effective way for training, the generalization needs to be addressed. In autoQSM, the robustness and accuracy of STAR-QSM could affect the performance of autoQSM. Moreover, for DL approaches for dipole inversion only, their performance still is affected by the performance of BFR methods.

Here, we propose a weakly-supervised approach for single-step QSM, denoted as wTFI. wTFI utilizes RESHARP [27] local fields as supervision for a multi-task learning of local tissue fields and QSM. For QSM quantitative evaluation, 9 multi-orientation datasets are utilized with comparison to TKD [24], MEDI [12], and STAR-QSM, using COSMOS result as a reference. In addition, the local tissue fields of wTFI are qualitatively compared with BFR methods, SHARP [22], RESHARP [27], PDF [13], and LBV [34]. More qualitative analysis is performed on single-orientation datasets and clinical datasets.

WTFI was trained using a 3D convolutional neural network (CNN) with an encoder-decoder structure, as shown in Fig.1. Since RESHARP produces more accurate local fields when compared with other BFR methods such as SHARP, PDF, and LBV. Local fields of RESHARP were utilized as supervision in the training paradigm. However, RESHARP suffered from brain erosion at the expense of losing information of local field and thus QSM measures in these regions. To address this problem, wTFI utilized a domain adaption technique to recover the lost information at brain edges for whole brain QSM.

During training, wTFI took 4 inputs (2 groups of total fields and brain masks) to get 4 outputs (2 groups of local fields and susceptibility maps). For the total field $f_{T1}$ and brain mask $m_{1}$, the corresponding outputs were $\chi_{1}$ and $f_{L1}$, and $f_{T2}$ and $m_{2}$ with outputs $\chi_{2}$ and $f_{L2}$. Define $m_{2}$ as brain mask obtained using human brain extraction tools from the MR magnitude data; $m_{1}$ as the eroded brain region after RESHARP, which has a smaller region than $m_{2}$ ($m_{2}>m_{1}$); $f_{L_{RESHARP}}$ as RESHARP local field; $f_{T}$ as the total field estimated from phase images; $f_{T1}$ as the total field in brain region $m_{1}$, $f_{T1}=m_{1}f_{T}$; $f_{T2}$ as the total field in brain region $m_{2}$, $f_{T2}=m_{2}f_{T}$.

The loss function consisted of five terms. The first loss term, $L_{\chi_{1}}$, was imposed on $\chi_{1}$. Leveraging the local field results from RESHARP $f_{L_{RESHARP}}$ as a weak supervision, $\chi_{1}$ convoluted with the dipole kernel $d$ should satisfy the well-established QSM inversion physical model.

where $W$ is a data-weighting factor which can be the magnitude image or noise weight matrix, $\ast$ is the convolution operator. Since noise is unknown and spatially variant in the local field measurements, the nonlinear dipole convolution data consistency loss was used to get more robust QSM estimates as conventional QSM methods [16, 19].

Next, we included $L_{f_{L1}}$ on local field output $f_{L1}$,

Next, a data consistency loss was imposed on $\chi_{2}$ and $f_{L2}$.

Furthermore, define the susceptibility consistency loss between $\chi_{1}$ and $\chi_{2}$ inside $m_{1}$.

The total variation (TV) loss $L_{TV}$ on $\chi_{2}$ serves as regularization for edge preserving and denoising on the QSM.

where $G_{x}$, $G_{y}$, and $G_{z}$ are gradient operator in x, y, z directions.

The total loss $L_{total}$ was the weighted sum of 5 loss functions. $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$, and $\lambda_{4}$ were loss weights.

After training, the trained DL model only took the whole brain total field $f_{T2}$ and brain mask $m_{2}$ to get the local field and susceptibility map.

Multi-orientation QSM Data 9 QSM datasets were acquired using 5 head orientations and a 3D GRE scan with voxel size 1x1x1 mm³ from 3T MRI scanners, shared by Dr. Jongho Lee [33]. QSM data processing was implemented as below, offline GRAPPA [7] reconstruction used to reconstruct magnitude and phase images from saved k-space data, BET (FSL, FMRIB, Oxford, UK) [25] for brain mask extraction, the Laplacian method [10] for phase unwrapping, and RESHARP [27] with spherical mean radius 4mm for BFR. COSMOS results were calculated using the 5 orientation data with image registration using FLIRT (FSL, FMRIB, Oxford, UK) [8, 9]. TKD, MEDI, STAR-QSM were performed to get QSM estimates at normal head position.

For wTFI training, leave-one-out cross validation was used. For each dataset, a total of 40 scans from other 8 datasets were used for network training. wTFI was trained using patch-based neural network with patch size 96x96x96. Around 2000 patch pairs of total fields, brain masks, RESHARP local fields with patch size 96x96x96 with an overlapping of 16.6 percent between adjacent patches were cropped. After training, the trained DL model took the whole brain total field and brain mask of the leave-one dataset to get the local field and QSM.

QSM images at the eroded brain region were quantitatively compared with respect to COSMOS maps, using peak signal-to-noise ratio (pSNR), normalized root mean squared error (NRMSE), high-frequency error norm (HFEN), and structure similarity index (SSIM). Due to no ground-truth, the local fields of SHARP, RESHARP, PDF, LBV, and wTFI were qualitatively compared.

Single Orientation QSM Data 200 QSM datasets were collected on a 3T MRI scanner (GE Healthcare MR750) from a susceptibility-weighted software application (SWAN, GE Healthcare). The data acquisition parameters were as follows: in-plane data matrix 320x256, FOV 24 cm, voxel size 0.5x0.5x2.0 mm³, 4 TEs [10.4, 17.4, 24.4, 31.4] ms, TR 58.6 ms, and total acquisition time 4 minutes.

Complex multi-echo images were reconstructed from saved k-space data. The brain masks were obtained using the SPM tool [3]. After BFR using the RESHARP with spherical mean radius 4mm, susceptibility inversion was performed using TKD, STAR-QSM, and MEDI. For wTFI training, 8000 patches of total field map, brain mask and RESHARP results from 100 datasets with patch size 128x128x64 were used for training.

Clinical Data 150 clinical data were acquired using susceptibility-weighted angiography (SWAN, GE Healthcare, Waukesha WI) on a 3T MRI scanner (GE Healthcare MR750) with the following data acquisition parameters: in-plane data matrix 288x224, FOV 22 cm, slice thickness 3 mm, first TE 12.6 ms, echo spacing 4.1 ms, 7 echoes, TR 39.7 ms, pixel bandwidth 244 Hz, and total acquisition time of about 2 minutes.

Complex multi-echo images were reconstructed from raw k-space data. The brain masks were obtained using the SPM tool [3]. RESHARP method with spherical mean radius 6mm were used for BFR. For wTFI training, about 2000 patch pairs of total fields, brain masks and local fields from 100 datasets with patch size 128x128x64 were used.

Multi-orientation Data Table.1 summarized the quantitative metrics from 4 QSM methods on 9 datasets with COSMOS map as a reference. Compared to TKD, STAR-QSM, and MEDI, wTFI achieved the best metric scores in pSNR, NRMSE, and second in HFEN. Fig.2 showed the total fields and local fields from a representative dataset. Residual background fields showed up in SHARP (b), PDF (d), and LBV (e) results. SHARP (b) and RESHARP (c) suffered from brain erosion. PDF and LBV results showed strong shading artifacts and erroneous BFR. wTFI (f) produced RESHARP-like local fields and preserved the whole brain without erosion. Fig.3 displayed the QSM images. TKD (a) and MEDI (b) results showed streaking artifacts (a-b, black arrows). MEDI images lost details due to oversmoothing. STAR-QSM (c) images showed good quality with slight image artifacts. wTFI (d) produced high quality QSM and recovered the susceptibility information of brain edges (d, white arrows).

Single-orientation Data Fig.4 displayed the QSM estimates of a single-orientation dataset. Visually, wTFI (d) outperformed TKD (a), MEDI (b) and STAR-QSM (c) with invisible streaking and shading artifacts. In addition, WTFI was capable of recovering the susceptibility information at brain boundaries.

Clinical Data Fig.5 displayed SWI images and susceptibility maps calculated using wTFI from 6 clinical patients. wTFI results clearly showed the brain hemorrhage, microbleeds, calcification, and vessel malformation.

Ablation Study The ablation study was to investigate the neural network design. 3 neural networks were compared, (1) wTFI, (2) wTFI_m only trained with inputs $f_{T1}$ and $m_{1}$. (3) STAR-Net, which used $f_{T1}$ and $m_{1}$ as inputs and STAR-QSM maps as QSM labels for single-step QSM reconstruction.

Fig.6 showed QSM estimates of a single-orientation dataset. STAR-Net in eroded brain (b) showed STAR-QSM-like quality, but less black shadings (white arrows). Whole brain STAR-Net (c) and wTFI_m (d) showed larger errors close to brain boundary when applying to whole brain $f_{T2}$ and $m_{2}$ for QSM estimates (black arrows). This indicated the domain shift caused susceptibility quantification errors. Visual assessment showed that wTFI (e) produced better whole brain susceptibility maps.

In this work, a weakly-supervised DL method for single-step QSM was proposed. From quantitative evaluation, wTFI achieved high metric scores. However, COSMOS as reference has limitations aforementioned. Based on visual assessment, wTFI outperforms conventional background field removal and dipole inversion methods. Most importantly, wTFI is capable of recovering the QSM information close to brain edges.

The proposed wTFI for QSM reconstruction approach has several advantages. First, wTFI employs RESHARP results as weak supervision and does not require QSM labels. Second, wTFI performs a multi-task learning, which estimates local fields and susceptibility map. Third, wTFI is able to recover magnetic susceptibility of anatomical structures near the edges of the brain. wTFI has a limitation that the wTFI performance is affected by RESHARP.

We thank Professor Jongho Lee for sharing the multi-orientation QSM datasets.