The magnetic field perturbation B→M(r→) measured with MRI, which is caused by the tissue magnetization of the underlying magnetic susceptibility distribution χr→ of the biological tissue, is defined in Eqs. 1a–1c according to Marques and Bowtell (2005)

where M→r→ is the magnetization μ0 the magnetic field constant and r→ the coordinate position in image space. Eq. 1c can be expressed as a convolution between the magnetic field perturbation B→Mr→ and the point dipole response dr→ (Li and Leigh, 2004; Marques and Bowtell, 2005):

where k→ are the coordinates in k-space, kzp the object coordinates projection onto the main magnetic field that is assumed to align along z-direction and θ is the angle between the z-direction and r→ (Liu et al., 2009). Assuming first order approximation for non-ferromagnetic materials χ≪1, separating near- and far-field contributions via the Lorentz sphere (Schweser et al., 2016) and employing the relationship between the susceptibility distribution and the resulting magnetic field variations (Li and Leigh, 2004; Marques and Bowtell, 2005), the relationship between the magnetic susceptibility distribution and the magnetic field perturbation is given by Eq. 3.

Since the dipole kernel dk→ has zeros at the magic angle (θ≈54.7°) and low values in the vicinity of the conical surface, this relationship is ill-posed, posing a challenge in calculating the susceptibility distribution. Therefore, algorithms including proper regularization techniques, such as MEDI (Liu et al., 2012) or HEIDI (Schweser et al., 2012), are required to determine susceptibility maps from a single MRI scan. Alternatively, the approach ‘calculation of susceptibility through multiple orientation sampling’ (COSMOS) (Liu et al., 2009) relies on at least three MRI scans with varying orientations of the object (e.g., the head) to the main magnetic field to shift the undetermined k-space regions and comprehensively fill the entire k-space. Based on Eqs 2, 3 it becomes obvious that the image resolution as well as the orientation of the FoV with respect to the main magnetic field impacts the magnetic field perturbation, where tilting the FoV leads to a change in the alignment of the FoV axes with physical axes describing the main magnetic field. For in-silico experiments and generation of training data, tilting of the FoV was mimicked by tilting the main magnetic field direction.

The concept behind adaptive convolution is the selection of the most appropriate set of convolutional filter weights for specific data sets, by providing additional information (Kang et al., 2020). When the convolutional filter weights are considered as points on a low-dimensional manifold in the high-dimensional filter weight space, the weights move on the manifold as a function of the additional information (Figure 1A). Hence, the convolution filter weights change adaptively as a function of the presented side information. The processing of this additional information as well as the computation of convolutional filter weights is performed by the filter manifold network (FMN), a fully connected feed forward network. When using a tilted FoV and/or anisotropic voxel sizes for MRI data acquisition, the shapes and edges of the anatomical structures in the image change smoothly. Consequently, by correlating these changes in the image with the additional information provided via the FMN, the network model is expected to learn the relationship between the susceptibility map and B-field distribution more easily. While the filter manifold is learned by the FMN during training, in inference the FMN weights remain fixed, however, the corresponding weights of the convolution operation change adaptively. Adaptive convolution is employed to adjust network parameters to data-specific attributes, selecting appropriate parameters for the susceptibility map calculation.

Adaptive convolution layers are built from fully connected feed forward networks as shown in Figure 1B. The FMN receives the side information s→ consisting of six input parameters, i.e., voxel-size v→ and the image orientation o→, at the input layer. The information propagates through the linear layers with ELU activation function f, whereby the dimensions of the extracted feature maps of the hidden layers are gradually increased (12, 48, 196). The number of extracted features at the final layer depends on the placement of the adaptive layer in the network to match the size of the corresponding convolution kernel. Since we placed the adaptive convolution layer in the first encoding block, the final FMN layer encodes 13,824 parameters that are reshaped into the 5D-tensor with shape [16, 32, 3, 3, 3] encoding the filter weights w. Thereafter, the input feature map X from the previous convolution layer is convolved with the adaptive convolutional filter k yielding the output feature map Y. The adaptive layer weights, and consequently the weights of the FMN, are updated and adjusted during network training. This enables the FMN to associate the side information with the filter manifold.

The adaptive U-Net (Figure 1C) computes susceptibility maps from the same-sized local magnetic field perturbation (BM→r→). The network relies on a standard 3D U-Net (Ronneberger et al., 2015) with 16 initial channels, concatenations between encoder and decoder and batch normalization layers (Ioffe and Szegedy, 2015). A skip connection between input and output was included to enable residual learning (Kim et al., 2016), which increases the performance of the network (Jin et al., 2017) and helps to reduce vanishing gradients during training (He et al., 2016). Thus, the network model computes the difference between the magnetic field perturbation and the susceptibility map. The ELU activation function (Clevert et al., 2015) was used to further address vanishing gradients and due to its greater ability to generalize. The adaptive convolution layer was included in the first encoding block. The kernel sizes of the convolutional layers were (3x3x3) and for transposed convolution (2x2x2). As comparison, an adaptive U-Net with the adaptive layer included in the entire encoder (Supplementary Figure S1), in total 10 adaptive convolution layers, and an identical U-Net lacking adaptive convolution were constructed and trained similarly as well as evaluated against the adaptive U-Net. This adaptive U-Net is referred to as adaptive encoder U-Net.

The loss measure (Eq. 4) used for training the models was designed as such that the similarity of the susceptibility maps and magnetic field perturbations is measured, including model-based learning therewith. Lχ compares the ground truth susceptibility χgt with the reconstructed susceptibility map χrec. LB compares the ground truth magnetic field Bgt with the one computed via the forward model (Eq. 3), where the reconstructed susceptibility map χrec was convolved with the unit dipole response dv→,o→ given by the specific voxel-size v→ and image orientation o→. The inclusion of the forward field computation ensures that optimization occurs only towards solutions obeying this convolutional relationship. The weighing factor λ was set to λ=1, to ensure that the physical model is substantially enforced via the loss function.

To circumvent data sparsity of brain MRI data with multiple head orientations, voxel sizes and ground truth susceptibility maps, and consequently allow the network model to learn solving the underlying physical problem of field-to-susceptibility inversion, pre-training on synthetic data sets and transfer learning to in-vivo brain data was performed. Hence, the network models were pre-trained on 1,000 purely synthetic data sets for 500 epochs, detaching the model from developing a possible anatomic bias of brain data therewith. The number of data sets was chosen empirically as trade-off between using a sufficient amount of training data and network training time. The synthetic susceptibility maps (shape: 320 × 320 × 320 voxel) were composed of rectangles, polygons and ellipsoids with susceptibilities drawn from a Gaussian distribution Nμ=0,σ=0.25. The total number of shapes was randomly drawn, ranging from 80 to 150 rectangles and 200 to 300 ellipsoids. The number of polygons was set to 50 empirically. 300 ellipsoids of random dimension were used for the creation of the mask and randomly assembled in image border regions. Each shape was blurred individually with a Gaussian kernel from the normal distribution Nμ=0,σ=U00.8 to smooth the edges. σ was drawn from a uniform distribution U. The susceptibilities of overlaying shapes were averaged to increase the number of edges in the data sets. The magnetic field perturbation was obtained via fast forward convolution in k-space (Eq. 3). During network training, four patches of dimension 160 × 160 × 160 were randomly cut from each data set. Augmentation techniques included random scaling of the susceptibilities, addition of Gaussian noise, random flipping and random 90° rotations for the patches. To obtain magnetic field distributions of random voxel size and image orientations, the dipole kernel for computing the field distributions and the ground truth susceptibility maps were augmented. The voxel-size was randomly drawn from the Gaussian distribution Nμ=1,σ=1.5. Different FoV orientations were simulated by tilting the main field direction of the dipole kernel. The normalized field vector B→0=001T was randomly rotated around the x-, y- and z-axis. The x- and y-values were drawn from the Gaussian distribution Nμ=0,σ=11 and the z-values from Nμ=0,σ=15. Variations in voxel-size and image orientation occurred with probability 0.8, ensuring standard parameters in the training data, and were used to generate the dipole kernel for the forward convolution of susceptibility maps.

In a second step, the network models were tailored specifically to brain data by performing transfer learning to susceptibility maps reconstructed from multi-orientation brain data sets acquired at 3 T (Shi et al., 2022) and 7 T (Deistung et al., 2013). The training data consisted in total of 109 data sets of varying orientations. Since the data sets were originally acquired with isotropic voxel sizes, we applied trilinear interpolation to augment the data sets to anisotropic voxel sizes. The voxel-size was randomly drawn from the distribution Nμ=1,σ=1.5. During transfer learning, it was randomly selected with a probability of 0.5 whether the measured local magnetic field perturbation or whether the magnetic field computed by fast forward convolution of the ground truth COSMOS susceptibility map with the dipole kernel was used as input data. In the latter case, the susceptibilities of the COSMOS maps were scaled by applying a randomly chosen multiplication factors (interval [0, 2]) to further augment the training data. By performing random scaling throughout transfer learning, the network model not only optimizes toward the susceptibility distribution of the provided data, but also maintains the susceptibility scale invariance. Four patches with dimensions of 96 × 96 × 96 were randomly cut from the data sets and padded to 128 × 128 × 128 to ensure border regions in the data. Transfer learning was performed for 30 epochs.

The AdamW optimizer β1=0.9,β2=0.99,ε=1e−08,λ=0.01 (Loshchilov and Hutter, 2017) was used during training and transfer learning. A cosine annealing learning rate schedule with warm restarts (T0=500,Tmult=0.5,ηinit=0.001,ηmin=1e−08) was used for network pre-training (Loshchilov and Hutter, 2016). For the transfer learning, the parameters of the adaptive layer remained fixed, the learning rate for all parameters in the decoding branch of the model was set to η=0.001 and to η=0.00001 for all other parameters. All augmentations, interpolations and dipole convolutions were performed on the GPU. Automated mixed precision together with cuDNN-benchmarking was used during distributed training of the network models on four NVIDIA A100 GPUs (NVIDIA Corporation, Santa Clara, CA, United States). It took approximately 25 h to train the model for 500 epochs and additional 3 h for transfer learning. Training the adaptive U-Net with four patches of voxel-size 160 × 160 × 160 requires around 26 GB of GPU memory per used GPU. The application of the trained adaptive U-Net to a local field map with a dimension of 192 × 224 × 160 voxels requires 18 GB CPU memory and takes approximately 14.3 s on a AMD EPYC 7713 CPU (Advanced Micro Devices Corporation, Santa Clara, CA, United States) or 1.6 s on a NVIDIA A100 GPU.

The adaptive U-Net and the conventional U-Net without and with transfer learning were evaluated on synthetic data sets as well as in vivo brain data sets (Table 1). Gradient-echo data were acquired on a Siemens Magnetom Vida (Siemens Healthcare GmbH, Erlangen, Germany) with a 64-channel head coil. One healthy subject was scanned three times. Two measurements were conducted with the isotropic voxel size of 1 mm × 1 mm × 1 mm, but different FoV orientations, one was oriented purely axial (o→ = [0, 0, 1]^T) and one was tilted by 12-degree from the axial orientation (o→=[–0.15, −0.20, 0.97]^T) (Table 1, this study (1)). A third scan with 0.57 mm × 0.57 mm × 2 mm voxels and a FoV tilted by 25° from the axial plane (o→=[–0.03, –0.42, 0.90]^T, Table 1, this study (2)) was acquired to test the limits of the network models. Since the side information array consisted of voxel sizes vi≥1, this data set serves to assess the network model performance on out-of-distribution data.

Post-processing of the data acquired within this study consisted of unwrapping (Abdul-Rahman et al., 2007) the phase images of each echo, dividing them by B0∗TEi∗γ and averaging them across the various echo times TE to achieve the magnetic field perturbation measured at 1 T. Sophisticated harmonic artefact removal for phase data (Schweser et al., 2011) with 10 different spherical kernels (1 to 10 voxels, regularized with truncated singular value decomposition: 0.1) was applied to reveal the local magnetic field perturbation. These SHARP-processed images were used for field-to-susceptibility inversion with the different deep learning models and homogeneity enabled incremental dipole inversion (HEIDI) (Schweser et al., 2012), while the maps computed using HEIDI served as reference. We referenced all susceptibility maps to the average susceptibility of the brain tissue within the field of view and stated susceptibility values in parts-per-million (ppm).

Data from other studies Deistung et al. (2013); Shi et al. (2022) were also considered for transfer learning and performance evaluations. Scan parameters of these studies are summarized in Table 1. For explicit details on data post-processing the reader is referred to the corresponding articles by Deistung et al. (2013) and Shi et al. (2022).

To evaluate the performance of the adaptive U-Net on pathological data, we simulated a data set containing four arbitrary shaped lesions. To this end, a COSMOS map from Shi et al. (2022) was interpolated to voxel-size 1 mm x 1 mm x 1.5 mm based on which a mask with four lesions was manually drawn. Various susceptibility values (χ=−1.0,−0.5,0.5,1.0) were assigned to the mask regions, smoothed using a Gaussian kernel (μ = 0, σ = 0.5), and incorporated into the COSMOS map. The magnetic field perturbation for the COSMOS data set with incorporated lesions was calculated via fast forward convolution in k-space (Eq. 3).

The reconstructed susceptibility maps of the adaptive U-Net were compared with other field-to-susceptibility methods including the adaptive U-Net without transfer learning, the conventional U-Net with transfer learning, the adaptive encoder U-Net, AFTER-QSM (Xiong et al., 2023) and HEIDI. We downloaded the AFTER-QSM model (implementation and trained network parameters) from the authors GitHub repository (https://github.com/sunhongfu/deepMRI/tree/master/AFTER-QSM, accessed and downloaded on 31.01.2024). Since AFTER-QSM is a DL-based field-to-susceptibility model specifically addressing resolution and orientation invariance and the trained parameters are publicly available, we chose this method for comparison.

Quality assessment of the reconstructed susceptibility maps has been performed on visual inspection in combination with quantitative image metrics. For this, the normalized root mean-squared-error (NRMSE, root mean squared error divided by the L2-norm of the ground truth and multiplied by 100), the structural similarity index (SSIM) and the peak signal-to-noise ratio (PSNR) were employed. To assess the amplitude of the computed susceptibilities, 400 k values were randomly sampled from the computed (c) and ground truth (g) susceptibility maps; ordinary linear least-squares fitting was then applied to determine the functional relationship g(c) (assuming g to be precise values). A mask of positions throughout the whole brain was randomly created and applied to the data sets. Hence, susceptibilities from identical positions were drawn and the random seed was identical for all models. Furthermore, the intensity values of the computed and ground truth susceptibility maps within cortical and deep gray (caudate nucleus, globus pallidus, thalamus, putamen, hippocampus, red nuclei, substantia nigra) matter brain regions were assessed via scatter plots and ordinary least-squares fitting. The gray matter regions were identified using DL-based segmentation (Billot et al., 2023) on the COSMOS map and then reduced in size by applying erosion with a 3 × 3 × 3 matrix of ones as structural element.

To verify the method of adaptive convolution itself and check, whether the adaptive layer has learned the filter manifold, different network activation maps and adaptive filter kernels were investigated for synthetic and brain data (Figure 2). Comparing the activation maps reveals changes in intensities and prominent edges in the respective feature maps (Figure 2, orange arrows). This change of structures is dealt in the adaptive model with rotating the filter kernel itself (Figure 2, orange arrows). On synthetic data and brain data, a constant change in edge sensitivity can be noticed by, e.g., changing the orientation of the edges, changing vertical edges to corners, or flipping filter values.

The computed susceptibility maps of the adaptive 3D U-Net and the conventional 3D U-Net generally resemble the ground truth susceptibility (Figure 3), however, differences in reconstructed object details and susceptibility intensities are observed (orange arrows). Indicated by lower NRMSE and higher SSIM the adaptive model computes susceptibility maps with higher similarity to ground truth data than the conventional U-Net. The plot of the susceptibilities of the adaptive model reconstruction (Figure 3C) against the ground truth susceptibilities shows a point distribution centered on a straight line, having highest point density around zero. The linear fitting showed a slope of 1.022, indicating a high degree of concordance between the calculated and ground truth susceptibility. For the conventional model (Figure 3D), the distribution between computed and ground truth susceptibilities is shifted toward a smaller susceptibility range, which is evident by a slope of 1.273 of the fitting curve.

The computed susceptibility maps of the different network models differ regarding metrics, detail resolution and amplitude of susceptibilities in the evaluation on three different data sets with varying FoV orientation and voxel-size (Figure 4) on 1 mm isotropic data from Shi et al. (2022). The susceptibility maps computed with the adaptive model are the ones mostly similar to the ground truth, visually and in metrics. Differences in the reconstructed maps are present in the globus pallidus (Figure 4A, white arrows) and in the optic radiation (Figure 4A, white circled arrows). While both adaptive models underestimate the susceptibility in the globus pallidus, an overestimation of susceptibilities is visible for the conventional U-Net and AFTER-QSM. The susceptibility map computed by AFTER-QSM (Figure 4A) exhibits a generally increased contrast level. As indicated by the difference maps, the adaptive U-Net and AFTER-QSM delineated most accurately the substantia nigra and red nuclei (Figure 4B, white circled arrows) on the tilted data set, whereas the susceptibility map of the conventional U-Net did not show a clear boundary between these structures. Similar findings are observed on the anisotropic data set, where all models have difficulties in delineating the substantia nigra and red nuclei (Figure 4C, white circled arrows), however, AFTER-QSM achieves the clearest delineation. The adaptive model without transfer learning showed visually the largest deviations to the COSMOS ground truth data set. AFTER-QSM and the adaptive model without transfer learning achieved the lowest image metrics. All visual findings are supported by the respective difference maps. The NRMSE deviates between the models in a greater range as the SSIM.

The scatter plots of randomly drawn susceptibilities (Figure 5) from the maps computed from in-vivo brain data (Figure 4) revealed a point cloud centered on a straight line for all models. On the isotropic and non-tilted data set (Figure 5A), the slope of 1.092 from the fitting line of the adaptive U-Net marks the closest agreement between the computed and ground truth COSMOS susceptibilities. The same model revealed its lowest slope (0.835) on the isotropic tilted data set (Figure 5B). The fitted slopes of the data reconstructed with the conventional U-Net varied between 0.635 and 0.801 (Figures 5G–I). The susceptibility maps reconstructed with the adaptive U-Net without transfer learning produced the broadest point distributions and slopes between 0.416 and 0.579 (Figures 5D–F). For AFTER-QSM (Figures 5J–L), the point distributions are more closely centered on the fitting line compared to those of the adaptive model without transfer learning. The slopes of the fitting line ranged from 0.427 to 0.534.

Scatter plots of susceptibilities within various gray matter regions (Figure 6) obtained from in-vivo brain data (Figure 4) revealed similar slopes like the ones determined across the whole brain (Figure 5). Agreement between the computed and ground truth susceptibilities was closest for the adaptive model as indicated by slopes varying between 0.824 and 1.128 (Figures 6A–C). The point distributions of the adaptive model without transfer learning (Figures 6D–F) and AFTER-QSM (Figures 6J–L) exhibit the largest dispersion, with fitting line slopes of about 0.5 being the lowest.

The performance of two different configurations of the adaptive U-Net differing in the number of included adaptive layers based on two different brain data sets is summarized in Figure 7. The adaptive U-Net with a single adaptive layer in the first encoding block (Adaptive U-Net, Figure 1C) and the adaptive encoder U-Net with multiple adaptive layers (in total 10) included from the first convolution block throughout the entire encoder (Adaptive Encoder U-Net, Supplementary Figure S1) compute susceptibility maps of comparable image quality. The NRMSE decreased slightly from 66.15 to 65.35 on the tilted isotropic data set (Figure 7A) and from 66.21 to 62.26 on the tilted anisotropic data set (Figure 7B). SSIM values of the maps computed with the different models were equal. The difference maps of the adaptive encoder U-Net revealed higher susceptibilities in various regions with high iron content than the single layer adaptive U-Net, e.g., in the right putamen, globus pallidus, and substantia nigra (see arrows and rectangles in Figure 7A) and lower susceptibilities in white matter fiber tracts, e.g., pyramidal tracts, (see white rectangles in Figure 7B). The adaptive encoder U-Net achieves similar demarcation of brain structures like the substantia nigra and red nuclei (Figure 7, white rectangles).

The impact of the side information array s on the adaptive U-Net is presented in Figure 8. If the side information is modified to values that do not match the used acquisition parameters, the image quality deteriorates and deviations to the ground truth increase for both, myelin- and iron-containing structures. The reconstructions with false side information yielded higher NRMSE and lower SSIM compared to maps reconstructed with the correct one (Figure 8). As highlighted by the error maps, the susceptibility maps deviate from the ground truth the more the side information differs from the correct one. The arrows and the rectangles indicate structures that substantially deviate from the susceptibility maps computed with the correct side information.

Susceptibility maps reconstructed from a data set with pathological lesions using different field-to-susceptibility models are presented in Figure 9. The overall susceptibility of the adaptive U-Net was closest to the COSMOS ground truth map, particularly in deep gray matter regions. The appearance of lesions were blurred across their original structural boundaries in all susceptibility maps, with the difference maps indicating largest deviations especially for lesion III (orange, top right hemisphere, susceptibilities of −0.68 ± 0.34 ppm). Lesion II (turquoise, susceptibilities of 0.72 ± 0.31 ppm) is shown at an enlarged scale, revealing subtle differences in the reconstructed lesion susceptibilities. The mean and standard deviations of lesion susceptibilities are summarized in Table 2, where the AFTER-QSM approach revealed closest agreement to the ground truth lesion susceptibilities. The adaptive U-Net had larger susceptibility deviations for lesions with negative susceptibility differences (e.g., calcified lesion III, IV) than for the ones with positive susceptibility differences (e.g., hemorrhagic lesion I, II).

The susceptibility maps from HEIDI and the adaptive U-Net with transfer learning visually have closest agreement in computed susceptibilities on the isotropic non-tilted data set (Figure 10A) and in iron-laden structures such as the putamen (Figure 10A orange arrows) and the dentate nuclei (Figures 10A,B turquoise arrows). Furthermore, the adaptive U-Net and AFTER-QSM achieved improved delineation of large (Figure 10A yellow arrows) and small brain vessels (Figure 10B yellow arrows) than the conventional U-Net. Substantial differences in susceptibilities are also present in the genu of the corpus callosum (Figure 10B orange arrows), with the conventional U-Net reconstructing positive rather than negative susceptibilities. Although the high-resolution dataset (Figure 10C) is out-of-distribution data, since the voxel size of the scan was not part of the network training and hence, not on the learned filter manifold, the contrast between the susceptibility map from HEIDI and the adaptive U-Net are comparable, while less contrast is visible in the susceptibility map from the conventional U-Net. A shadow effect is visible in the vicinity of large veins in the cerebellum of the reconstructed susceptibility maps of the adaptive U-Net that is more pronounced in the map of the adaptive U-Net without transfer learning and less visible in the conventional U-Net. The reason for the shadowing may be attributed to the fact that the local field distribution is sparsely sampled across the slice encoding direction because of the high voxel aspect ratio (voxel size in slice encoding direction divided by the voxel size in-plane) of 3.5. The AFTER-QSM map additionally yields large scale heterogeneities in the cerebellum. All susceptibility maps from the adaptive U-Net without transfer learning substantially deviate from the HEIDI susceptibility and the adaptive U-Net and are prone to artefacts. Like in Figure 4, the susceptibility maps of AFTER-QSM revealed a higher contrast, a more specifically-larger co-domain than those of the adaptive U-Net and HEIDI, especially for deep gray matter regions in the high-resolution data (Figure 10C).