Investigating the effect of oblique image acquisition on the accuracy of QSM and a robust tilt correction method

Purpose Quantitative susceptibility mapping (QSM) is used increasingly for clinical research where oblique image acquisition is commonplace, but its effects on QSM accuracy are not well understood. Theory and Methods The QSM processing pipeline involves defining the unit magnetic dipole kernel, which requires knowledge of the direction of the main magnetic field B^0 with respect to the acquired image volume axes. The direction of B^0 is dependent on the axis and angle of rotation in oblique acquisition. Using both a numerical brain phantom and in vivo acquisitions in 5 healthy volunteers, we analyzed the effects of oblique acquisition on magnetic susceptibility maps. We compared three tilt‐correction schemes at each step in the QSM pipeline: phase unwrapping, background field removal and susceptibility calculation, using the RMS error and QSM‐tuned structural similarity index. Results Rotation of wrapped phase images gave severe artifacts. Background field removal with projection onto dipole fields gave the most accurate susceptibilities when the field map was first rotated into alignment with B^0. Laplacian boundary value and variable‐kernel sophisticated harmonic artifact reduction for phase data background field removal methods gave accurate results without tilt correction. For susceptibility calculation, thresholded k‐space division, iterative Tikhonov regularization, and weighted linear total variation regularization, all performed most accurately when local field maps were rotated into alignment with B^0 before susceptibility calculation. Conclusion For accurate QSM, oblique acquisition must be taken into account. Rotation of images into alignment with B^0 should be carried out after phase unwrapping and before background‐field removal. We provide open‐source tilt‐correction code to incorporate easily into existing pipelines: https://github.com/o‐snow/QSM_TiltCorrection.git.


INTRODUCTION
The acquisition of oblique image slices, or an oblique slab or volume in 3D MRI, is common in clinical practice to facilitate radiological viewing of brain MRI.For example, axial slices are often aligned along the subcallosal line for longitudinal studies that require consistent repositioning of acquired images. 1 Alternatively, slices may be aligned perpendicular to the principal axis of the hippocampus for accurate hippocampal volume measurements and sharper hippocampal boundary delineation. 2Oblique slices are also acquired to reduce image artifacts from, for example, eye motion, resulting in localized blurring around the eyes and ghosting along the phase-encode direction. 3Note that acquiring oblique slices does not require the subject to rotate their head, as only the acquisition volume is tilted.Quantitative susceptibility mapping (QSM) [4][5][6][7] uses the information in the (conventionally discarded) phase component, (r), of the complex MRI signal from a gradient-echo sequence to calculate the tissue magnetic susceptibility, .A typical QSM pipeline includes three key steps: (1) phase unwrapping of wraps present due to (r) being constrained to the [−π, π) interval; (2) background field removal separating the local field perturbations due to internal  sources inside the volume of interest (eg, the brain), ΔB int (r), from unwanted background field perturbations due to external sources, ΔB ext (r); and (3) a local field-to-(r) calculation to solve an ill-posed inverse problem: where ℱ is the Fourier transform; ℱ −1 is its inverse; B 0 is the magnetic field strength in Tesla; and d z (k) is the z-component of the magnetic dipole in k-space d z (k) = z k 2 (see Equation 5).Calculation of d z (k) requires knowledge of the "z" direction of the main magnetic field, B 0 , with respect to the image volume acquired.Therefore, oblique acquisition must be taken into account within the QSM pipeline; otherwise, incorrect  estimates arise, as suggested by a preliminary study 8 and our preliminary data. 9With the increase in clinical applications of QSM, 10,11 accuracy in  estimates for oblique acquisition, typical in clinical protocols, is of paramount importance in ensuring smooth translation of QSM into clinical practice.However, accurate QSM accounting for oblique acquisition is nontrivial, and there are a number of techniques proposed to account for oblique acquisition in QSM, [12][13][14][15] including the most common methods of rotating the k-space dipole or the image volume into alignment with B 0 .The effect of these proposed tilt-correction techniques on susceptibility values has not been evaluated, and it is not known at which point in the QSM pipeline these techniques should be applied.Furthermore, it is not clear what is the optimal method for taking oblique acquisition into account in the QSM pipeline: Simply defining the dipole at an angle (see DipK or DipIm subsequently) has been shown to be nonoptimal. 8Therefore, the research presented here is the first quantitative and comparative evaluation of correction methods for oblique acquisition in QSM.We used a numerical phantom to carry out a comprehensive analysis of the effect of oblique acquisition on each step of the QSM pipeline, and propose three tilt correction schemes, analyzing their effects on susceptibility values when applied at different points in the QSM pipeline.We also acquired several images, in 5 healthy volunteers, with volumes tilted at different angles and performed the same analysis of the effects of tilting and correction schemes in vivo.We provide open-source tilt-correction code at https://github.com/o-snow/QSM_TiltCorrection.git, which uses the header information from NIfTI 16 format images to correctly orient image volumes and account for tilted acquisition for accurate QSM.

THEORY
To accurately model the magnetic dipole kernel required for the field-to- calculation and, in some cases, for background field removal, it is necessary to know where the magnetic field B 0 lies in the acquired MR images.Defining the two coordinate systems of interest 8 as the acquired image frame (û, v, ŵ) and the scanner frame (x, ŷ, ẑ), the main magnetic field can be written as B 0,im and B 0,sc in the image and scanner frames, respectively, as follows: In the case of nonoblique acquisition, the coordinate systems are aligned and B 0,im = B 0 ŵ in the image frame (Figure 1, left).
For the local field, ΔB int (r), to (r) calculation (Equation 2), the magnetic dipole kernel must be calculated.Throughout this paper, references will be made to the dimensionless k-space dipole, d z (k) (Figure 1, middle row), and the dimensionless "image-space dipole" defined in image space and Fourier-transformed into k-space, Nonoblique and oblique acquisition about the x-axis (u-axis) of axial slices (top row) with corresponding k-space dipoles (middle row) and image-space dipoles (bottom row).The image axes (u, v, w) and scanner axes (x, y, z) are shown in red and black, respectively.Note that the rotation axis is at the center of the image d z,im (k) (Figure 1, bottom row), kernels defined as follows 17,18 : where k = k u û + k v v + k w ŵ is the unit vector of k; B0 is the unit vector of B 0 ; V is the voxel volume; θ is the angle between B0 and r, the unit vector of r in image space, where r = √ u 2 + v 2 + w 2 ; and ℱ is the Fourier transform.The periodicity of the discrete Fourier transform constrains the boundaries of k-space, resulting in the dipole pattern becoming fixed along those boundaries.This causes a rotated image-space dipole to appear twisted, sheared, or distorted (Figure 1, bottom row).It is possible to obtain the direction of B0 relative to the (tilted) image axes from the image headers (eg, DICOM or Nifti format), and therefore to correctly calculate the magnetic dipole kernel using either Equation 5 or 6.

METHODS
To determine the optimal method for taking oblique acquisition into account in the QSM pipeline, we investigated three proposed tilt-correction schemes, and, for comparison, an uncorrected analysis pipeline (Figure 2): 1 RotPrior: Rotation of the oblique image into alignment with the scanner frame before phase unwrapping, background field removal, or the susceptibility calculation method.In this method, the dipole is defined in k-space in the scanner frame (using Equation 5); 2 DipK: The image is left unaligned to the scanner frame, and the dipole used is defined in k-space in the oblique image frame (using Equation 5).This is the default tilt-correction method implemented in popular QSM toolboxes. 15,19,20However, this method often requires the user to input the corrected B 0 direction, which is optional in many of these toolboxes; 3 DipIm: The image is left unaligned to the scanner frame, and the dipole used is defined in image space in the oblique image frame (using Equation 6); and 4 NoRot: The oblique image is left unaligned to the scanner frame, and the k-space dipole is mistakenly defined in the scanner frame (by wrongly assuming B 0 = B 0,sc = B 0,im in Equation 5) and is thereby misaligned to the true magnetic field direction, B 0 .This is the uncorrected method, which can easily result from users failing to input (the correct) B 0 direction.
These schemes are the general names of the methods that we applied at different points in the pipeline (ie, before phase unwrapping, background field removal, or susceptibility calculation) and for different methods or algorithms.For example, as no dipole kernel is necessary for phase unwrapping (and we substitute the dipole operations illustrated in Figure 2 with B0 orientation-independent unwrapping operations), we have called the only two schemes appropriate before unwrapping RotPrior and NoRot, where the image volume is rotated before unwrapping and after, respectively.
All rotations were carried out about the x-axis (u-axis) to simulate single oblique acquisition, the y-axis (v-axis) for confirmation, and about the y = x-axis (v = u-axis) to simulate double oblique acquisition.Rotations were undertaken using FSL FLIRT 21 with trilinear interpolation.To facilitate comparisons, all images left in the image frame after correction (DipK, DipIm, and NoRot) were rotated back into alignment with the scanner axes (see black arrow in Figure 2).Unless stated otherwise, all All tilt-correction schemes including the reference, nonoblique acquisition for rotations about the x-axis.The native (oblique) image space (û, v, ŵ) was transformed to ( û′ , v′ , ŵ′ ), aligned with the scanner frame.The black arrow denotes rotation into the scanner frame of reference.DipK, DipIm, and NoRot were rotated back into the reference (scanner) frame after correction to facilitate comparisons.RotPrior and NoRot still apply when no dipole is used processing and analysis operations were carried out using MATLAB (MathWorks, Natick, MA, USA).

Numerical phantom investigations
Multi-echo (TE = 4, 12, 20, and 28 ms) magnitude and phase images, from a numerical phantom, 22 with (originally) no phase wraps or background fields present, were used to independently investigate the effect of the three tilt-correction methods (described previously), and no correction, on each step in the QSM pipeline.We carried out these investigations with two image volumes: one unpadded with the original matrix size 164 × 205 × 205, and a second volume padded to 357 × 357 × 357.The padded matrix size was chosen as the long diagonal of the initial volume (padded to a cube: 205 × 205 × 205) and rounded up to the nearest odd integer.This was to ensure that none of the original frequency coefficients of the unit dipole field were cut off due to rotations about any of the three axes.An odd matrix size meant that there was a true center of rotation correctly located within a single central voxel.

Numerical phantom: susceptibility calculation
Local field maps, ΔB int (r), were calculated from a non-linear fit 23 over all echo times (for the most accurate field estimates 24 ) of the complex data set created by combining the magnitude and background field-free phase images.These local field maps, obtained from the supplied raw numerical phantom data, were free of any (synthetic) background fields or phase wraps, and therefore allowed investigation of the effect of oblique acquisition on  calculation alone.To simulate oblique acquisition, local field maps were rotated between ±45 • in 5 • increments.All tilt correction methods described (and no correction) were compared for three  calculation methods chosen to cover the two main approaches: direct non-iterative solutions (in k-space) and iterative solutions (in image-space).
The first method tested was direct, thresholded k-space division (TKD) 25,26 (from open-source software) where a modified dipole kernel was generated in k-space with values below a threshold, δ = 2∕3, replaced by the signed threshold value: The dipole was originally defined according to DipK and DipIm and then always thresholded in k-space.Susceptibility underestimation was corrected by multiplication with a correction factor, c  (δ), calculated according to. 27he second and third  calculation methods aim to iteratively solve for  through the minimisation of where M is a binary mask, W is a weighting term and R() is the data regularization term that reflects some prior information about .0][31] It was applied with R()= |||| 2 2 , a regularization parameter  = 0.003 (chosen through an L-Curve analysis 32 ), and W reflecting the spatially varying noise, and was also corrected for  underestimation. 27eighted linear total variation regularization (from the FANSI toolbox 19,33 ) with R() = || 1 ,  = 6.31 × 10 −5 (chosen through an L-Curve analysis) and W the magnitude of the complex data 33 was also tested.This method was chosen as total variation based iterative approaches were shown to produce the most accurate susceptibility maps in the 2019 QSM challenge 2.0. 34ean  values were calculated in five deep gray matter regions of interest (ROIs): the caudate nucleus, globus pallidus, putamen, thalamus and red nucleus.All susceptibility maps were compared using the root mean square error (RMSE) and QSM-tuned structural similarity index (XSIM) 35 metrics relative to the ground-truth susceptibility map at 0 • .

Numerical phantom: background field removal
For the background field removal step, local field maps from the numerical phantom required the addition of synthetic background fields, which were then removed following the three different tilt-correction methods (and no correction).After background field removal, the susceptibility maps were calculated from the resulting field maps using the  calculation method found to be optimal in the assessment described previously.
To investigate the effect of tilt-correction schemes on the background field removal step, synthetic background fields, ΔB ext (r) (Figure 3, bottom left), were added to the local field maps used in section 3.1.1.The background fields were calculated using the forward model, which is through a convolution, formulated as a multiplication in Fourier space, between the unit magnetic dipole field and a head-shaped susceptibility map 17,18 as follows: where χ head is the head-shaped susceptibility map, with soft tissue (−9.4 ppm) and bone (−11.4 ppm) 8,36 regions obtained by thresholding the magnitude (sum of squares over all echoes) and a pseudo-CT [37][38][39] image, respectively (Figure 3).The magnitude and pseudo-CT images were padded from their original matrix size of 164 × 205 × 205 to 512 × 512 × 512 to ensure that the edge effects from the periodic Fourier transform were minimized around the volume of interest.These synthetic background fields were then cropped back to their original matrix size and added to the local field maps obtained previously simulating a total field map, ΔB(r) = ΔB int (r) + ΔB ext (r).To simulate oblique acquisition, total field maps were rotated between ±45 • in 5 • increments.To remove the synthetic background fields, ΔB ext (r), from the tilted total field maps, three different state-of-the-art background field removal methods 40 were used, based on their widespread Method for calculating the synthetic background field from a head-shaped susceptibility map obtained by thresholding the numerical phantom magnitude image and a pseudo-CT image to delineate soft tissue and bone, respectively.The thresholded magnitude and pseudo-CT images were filtered for smoothness using a 3 × 3 × 3 box filter use, robustness, and accuracy. 40,41Projection onto dipole fields (PDF) 42 from the MEDI Toolbox 15 was used following tilt correction with all three correction schemes and no correction, because PDF is orientation-dependent (ie, it uses the dipole field d z (k) (Equation 5)).Laplacian boundary value (LBV) 43 from the MEDI Toolbox, 15 and variable-kernel sophisticated harmonic artifact reduction for phase data (V-SHARP) 44 from STI Suite, 20 were tested with RotPrior and NoRot only, as LBV and V-SHARP are orientation-independent methods (ie, they do not use the dipole field).Following rotation back into the reference frame (equivalent to RotPrior for the susceptibility calculation step), susceptibility maps were calculated from all local field maps using the iterative Tikhonov regularization (regularization parameter  = 0.003), as this was found to be optimal.Susceptibility maps were compared using RMSE and XSIM 35 relative to the ground-truth susceptibility map at 0 • .

Numerical phantom: phase unwrapping
To investigate the effect of tilt correction on phase unwrapping, the synthetic background fields added in section 3.1.2induced phase wraps when the phase was constrained to the [−π, π) interval, which were then unwrapped.Susceptibility maps were then calculated from these unwrapped field maps using background field removal and susceptibility calculation algorithms found to be optimal in the experiments described in sections 3.1.1and 3.1.2.
To investigate the effect of tilt correction on phase unwrapping, phase wraps were introduced into the wrap-free numerical phantom images through the additional synthetic background field described previously.From each total field map at each angle, ΔB(r) = ΔB int (r) + ΔB ext (r), multi-echo unwrapped phase images were simulated by scaling the tilted total field maps at each TE according to (r, TE) =  ⋅ TE ⋅ ΔB(r).At every tilt angle, a complex data set (S) was made from the multi-echo magnitude images (M) and simulated phase images () using S = M(r, TE)e i(r,TE) , which constrained the phase to the range [−π, π), resulting in phase wraps.A wrapped total field map was calculated via a nonlinear fit over all TEs, 23 which then underwent phase unwrapping using the commonly used Laplacian, 15,45 SEGUE 46 (https://xip.uclb.com/product/SEGUE),and ROMEO 47 (https://github.com/korbinian90/ROMEO)techniques with the NoRot and RotPrior tilt-correction methods.After rotating all of the unwrapped images back into the reference frame (equivalent to RotPrior for the background field removal step and susceptibility calculation step), susceptibility maps were then calculated with PDF 42 background field removal and susceptibility calculation using iterative Tikhonov regularization (regularization parameter  = 0.003), as we found these to provide optimal results.

3.2
Investigations in vivo

In vivo: MRI acquisition
All acquisitions were performed having obtained informed consent, and with approval by the local ethics committee.The 3D gradient-echo brain images of 5 healthy volunteers were acquired on a 3T Siemens Prisma-Fit MR system (National Hospital for Neurology and Neurosurgery, London, United Kingdom) using a 64-channel head coil across a range of image volume orientations.Note that the volunteers did not tilt their head but remained in the same position throughout the experiment.The image volume was tilted about the x-axis, as this is the most common in clinical practice, from −45 to +45 • in 10 • increments, with the reference image at 0 • representing a nonoblique acquisition (subject 5 angles were between ±15 • ) .Each image volume was acquired in 3 min 23 s with TR = 30 ms, TEs

3.2.2
In vivo: phase unwrapping For all angles/volumes, a total field map and a noise map were obtained using a nonlinear fit of the complex data 23 from the MEDI toolbox. 15A brain mask was created using the brain extraction tool 49 with default settings applied to the final echo magnitude image of the reference 0 • volume for a conservative brain mask estimate.This brain mask was then registered to all oblique-acquired volumes to maintain consistency.As with the numerical phantom, both the RotPrior and NoRot correction schemes were applied.Residual phase wraps were then removed using Laplacian unwrapping, 45 SEGUE, 46 and ROMEO. 47To investigate the effect of the correction schemes on this step in the pipeline, unwrapped total field maps were rotated back into the reference frame (equivalent to RotPrior for the background field removal step and susceptibility calculation step), and susceptibility maps were created using PDF background field removal and susceptibility calculation with iterative Tikhonov regularization ( = 0.017 chosen through an L-curve analysis).
As in the numerical phantom, and also due to very slight unavoidable changes in subject position between scans, the unwrapped field maps and susceptibility maps were registered into the reference image space to facilitate comparisons of results in vivo.To carry out this registration, the magnitude image (added in quadrature over all echoes) for each angle was rigidly registered to the 0 • magnitude using NiftyReg, 50 resulting in a transformation matrix per angle/volume, which was applied to bring all angles/volumes into the same common reference space.

3.2.3
In vivo: background field removal The ROMEO unwrapped field maps described in section 3.2.2 for volumes at all angles, and before any registrations or rotations, were used to investigate the effect of oblique acquisition on background field removal.The ROMEO technique was chosen because it has been shown to outperform 47 PRELUDE 51 and BEST-PATH. 52As for the numerical phantom, for each field map, at each angle, background fields were removed using PDF 42 with all tilt-correction schemes and no correction, and using LBV 43 and V-SHARP 44 with only RotPrior and NoRot.For all three background field removal methods, the brain mask was eroded by four outer voxels. 53RotPrior was performed twice: with mask erosion either before or after the rotation, to compare the effects of interpolation, particularly along the boundaries of the field map on PDF and V-SHARP, as it is known that boundary effects arise in these background field removal methods. 40or comparison purposes, after rotation and registration of the local field maps back into the reference frame (equivalent to RotPrior for the susceptibility calculation step), susceptibility maps were calculated from the local field maps using iterative Tikhonov regularization Effect of tilt correction before phase unwrapping in the numerical phantom.Phase-unwrapped field maps and the resulting susceptibility maps at 15 • for the NoRot (column b) and RotPrior (column c) tilt-correction methods relative to the reference (column a).Rotation of the wrapped field maps before phase unwrapping with Laplacian, SEGUE, and ROMEO techniques results in errors along phase wraps and incorrect unwrapping, leading to prominent artifacts in the final susceptibility maps ( = 0.017, chosen with an L-Curve).Local field maps and susceptibility maps were compared with RMSE and XSIM metrics (XSIM only for the susceptibility maps) averaged across all subjects relative to the 0 • reference image.

In vivo: susceptibility calculation
To investigate the effect of oblique acquisition on the  calculation step in the pipeline, we used the local field maps following ROMEO unwrapping and LBV background field removal (described in section 3.2.3),before any registrations or rotations.The LBV method was chosen because it is orientation-independent, thereby allowing our analysis to focus on the effect of oblique acquisition and the different correction schemes on the  calculation step alone.The three tilt correction schemes (and no correction) were compared using the same three  calculation methods as for the numerical phantom: TKD, iterative Tikhonov regularization with a regularization parameter  = 0.017 from an L-curve analysis, and weighted linear TV with a regularization parameter  = 6.31 × 10 −5 found also from an L-curve.The resulting susceptibility maps were transformed into the reference space as described in section 3.2.2.The same regions of interest (ROIs) as in the numerical phantom were investigated and obtained by registering the EVE 54 magnitude image with the reference magnitude image (at the first TE) and applying the resulting transformation to the EVE ROIs.Mean  values were calculated in these ROIs for all tilt angles and all correction schemes.The RMSE and XSIM measures averaged across all subjects were also used to compare the susceptibility maps.

Numerical phantom
All numerical phantom results shown here are for rotations of image volumes about the x-axis with unpadded matrices.Note that acquisitions tilted about the y-axis Effect of different tilt-correction schemes on QSM with three background field removal methods in a numerical phantom.Susceptibility maps for QSM-tuned structural similarity index (XSIM) comparisons were calculated with iterative Tikhonov regularization.For projection onto dipole fields (PDF) (A), the XSIM metric shows that RotPrior gives the most accurate susceptibilities, with DipIm performing the worst.When using PDF with DipK, striping artifacts (C, red ellipse) arise in the local field maps for tilted acquisitions.Laplacian boundary value (LBV) and variable-kernel sophisticated harmonic artifact reduction for phase data (V-SHARP) (C,D) are shown to be largely unaffected by oblique acquisition with differences arising primarily from rotation interpolations and y = x-axis, as well as images with padded matrices, all gave similar results (see Supporting Information Figures S1-S3).We chose to display these results as the padded matrix size leads to increased computation time, which is not recommended in a practical setting, and the x-axis is the most common axis of rotation for oblique acquisition.
When wrapped phase images are rotated before phase unwrapping with the correction scheme RotPrior, artifacts arise for Laplacian, SEGUE, and ROMEO unwrapping (Figure 4).The SEGUE and ROMEO techniques appear to fail with RotPrior, with SEGUE removing a portion of the brain mask and ROMEO leaving residual phase wraps (Figure 4C).
When using PDF for background field removal, Rot-Prior is the most accurate method, and the largest errors arise from DipIm and NoRot (Figure 5).Striping artifacts are present in the local field map from the DipK method (Figure 5C).The LBV and V-SHARP are shown to be largely unaffected by oblique acquisition.
Figure 6 summarizes the mean susceptibility in the caudate nucleus and thalamus, alongside XSIM measurements across all angles for all three  calculation methods (RMSE results are similar and are shown in Supporting Information Figure S4).The TKD and iterative Tikhonov methods are most accurate with RotPrior, and least accurate when the dipole is misaligned to the main magnetic field (NoRot).Weighted linear TV is relatively robust to oblique acquisition with RotPrior and DipK performing similarly.However, DipK shows more variability in  at the ROI level than RotPrior.Weighted linear TV with DipIm fails at nonzero angles, and NoRot results in the largest errors.Example susceptibility maps are shown in Figure 7, highlighting the widespread  errors that arise when the magnetic dipole is defined incorrectly (NoRot).

In vivo
In vivo, Laplacian, SEGUE, and ROMEO phase unwrapping with RotPrior have the same image artifacts as in the numerical phantom (not shown here) compared with NoRot, with incorrect identification of phase wraps when wrapped field maps are rotated before phase unwrapping.
Figure 8 shows that PDF background field removal is most accurate with RotPrior and least accurate with DipIm followed by NoRot, confirming the results obtained in the Mean susceptibilities in the caudate and thalamus (top rows), and XSIM (bottom row) across all tilt angles for all tilt-correction schemes, plus all three  calculation methods in the numerical phantom.The RMS error (RMSE) measurements shown in Supporting Information Figure S4 agree with the XSIM findings.NoRot performs worst across all angles.RotPrior is the most accurate tilt-correction scheme.For weighted linear total variation (TV), DipK and RotPrior have similar XSIM values but the mean thalamus  varies more over angles with DipK.Note that DipIm is not shown for weighted linear TV, as this method fails numerical phantom (Figure 5).The average XSIM differences between tilt-correction schemes in vivo (Figures 8  and 9) are smaller than in the numerical phantom (Figures 5 and 6), most likely due to issues inherent to in vivo acquisition, including motion and greater noise.Striping artifacts are present in the DipK method for PDF in the local field maps before re-orientation for comparison purposes (Figure 8C).Rotation interpolation obscures these artifacts in the in vivo images.The LBV and V-SHARP methods are shown to be largely orientation-independent in the in vivo case, as expected.
When RotPrior was performed with mask erosion before rotating the total field map, artifacts arose along the boundaries of the PDF local field map (Supporting χ maps and difference images illustrating the effects of all tilt correction schemes in the numerical phantom.An axial and a coronal slice are shown for a volume tilted at 25 • and reference 0 • volume with all χ maps calculated using the iterative Tikhonov method.The ROIs analysed are also shown (bottom left).RotPrior performs the best while NoRot results in substantial χ errors across the whole brain.The results from TKD and weighted linear TV (not shown) are very similar Information Figure S5).The PDF method performs more robustly if mask erosion is carried out after rotation, whereas V-SHARP appears to perform equally well in both scenarios.
Figure 9 shows the effect of all tilt-correction schemes on susceptibility calculation in vivo and confirms that NoRot results in the largest susceptibility errors and that RotPrior is consistently the most robust tilt-correction method compared with other methods.Both RotPrior and DipK perform better than DipIm, in agreement with results in the numerical phantom (Figure 6).Difference images (Figure 10) also confirm those obtained in the numerical phantom (Figure 7).Subtle effects found in several of the numerical phantom ROIs (Figure 6) were not apparent in vivo.

DISCUSSION
We have shown that oblique acquisition must be accounted for in the QSM pipeline to ensure accurate susceptibility estimates throughout the brain.For all background field removal and susceptibility calculation methods tested, if the magnetic dipole kernel is left misaligned to the B0 direction (NoRot), which can arise from user error in popular QSM toolboxes, then significant susceptibility errors result.Through the analysis of the effect of tilted acquisition on a numerical phantom and 5 healthy volunteers in vivo, we have shown that any rotations that are applied to a wrapped field map before phase unwrapping will result in incorrect unwrapping, using Laplacian, SEGUE and ROMEO unwrapping techniques, and subsequent artifacts in the resulting QSMs.Results indicate that, for PDF background field removal, rotating the image into the scanner frame and using a k-space dipole defined in the scanner frame (RotPrior correction method) provides the most accurate susceptibility maps.If no image rotations are desired, due to unwanted interpolation effects, LBV or V-SHARP are recommended, as they are largely unaffected by oblique acquisition.Both TKD and iterative Tikhonov susceptibility calculation methods provide the most accurate results when local field maps are rotated into alignment with the scanner axes and a k-space dipole, defined in the scanner frame, is used (Rot-Prior).The same conclusion holds for weighted linear TV, but susceptibility calculation can be carried out in the oblique image frame without any rotations, provided the Effect of different tilt-correction schemes on background field removal in vivo.The average XSIM measurements across all subjects were used to compare the  maps calculated with iterative Tikhonov regularization after background field removal to the nonoblique (0 • ) reference map.The PDF method (A) has the highest XSIM with RotPrior and the lowest XSIM with DipIm, followed by NoRot, confirming the results in the numerical phantom (Figure 5).Striping artifacts are found in local field maps when using DipK and PDF (C, red ellipses) but are obscured after rotation and registration back into the reference 0 • space due to interpolation.The LBV (B) and V-SHARP (D) methods are shown to be unaffected by oblique acquisition in vivo as well as in the numerical phantom (Figure 5B and 5D).Error bars represent the SD of the mean XSIM across subjects correct B 0 direction is used in defining the k-space dipole (DipK correction method).We therefore recommend rotating the total field map into alignment with the scanner frame after phase unwrapping but before background field removal.
Both the numerical phantom and in vivo results indicate that when wrapped phase images are rotated before phase unwrapping (with the correction scheme RotPrior), artifacts arise for the Laplacian, SEGUE, and ROMEO unwrapping methods.This is probably due to interpolation errors along phase wraps (Figure 4).Therefore, any phase unwrapping must be carried out in images left in the same orientation as acquired, with rotations only being applied afterward to avoid artifacts.
When using PDF for background field removal, numerical phantom and in vivo results show that RotPrior consistently provides the most accurate susceptibility maps, while NoRot performs the worst (Figures 5 and 8).Striping artifacts arise in local field maps in both the numerical phantom and in vivo when using PDF with the DipK method.First identified by Dixon, 8 these striping artifacts are present due to the violations in circular continuity Average XSIM plots over all angles for all tilt-correction schemes and all three  calculation methods across all subjects in vivo.The RMSE measurements in Supporting Information Figure S4 agree with these XSIM findings.These results are similar to those in the numerical phantom (Figure 3), with RotPrior consistently reporting higher XSIM measures than other methods and NoRot performing worst across all methods.At nonzero tilt angles, XSIM has a respectively high/low baseline level arising from rotation and registration interpolations.DipIm fails for weighted linear TV, and therefore is omitted from the plots in the last column.Error bars represent the SD of the mean XSIM across subjects when defining the tilted dipole in k-space and using the inverse discrete Fourier transform to transform the susceptibility maps into image space (which enforces periodicity; see Supporting Information Figure S7).Striping artifacts arise from regions of high susceptibility changes, such as on the brain boundaries (Figures 5C and 8C).DipIm also resulted in poor background field removal, most likely due to fitting to the incorrect twisted or sheared unit dipole field (bottom row of Figure 1).To avoid artifacts and robustly achieve background field removal with PDF, total field maps must be rotated into alignment with the scanner frame before PDF background field removal, as it is then possible to use the nonoblique dipole, which does not violate circular continuity (Supporting Information Figure S8).We showed that LBV and V-SHARP were mostly unaffected by oblique acquisition, with the differences between zero and nonzero tilt angles arising solely from rotation interpolation effects.
Given that RotPrior is the most accurate method for PDF, the typically necessary mask erosion must be carried out after rotation into the reference space.Artifacts that arise along the boundaries of the local field map if erosion is carried out before rotation (Supporting Information Figure S5A,B) probably arise from distortion of the dipolar background fields due to interpolation at the edges.In contrast, V-SHARP does not show substantial differences with mask erosion before v. after rotation, which suggests that the interpolation may not substantially affect the harmonic nature of the background fields on which this method relies.
We found that TKD and iterative Tikhonov regularization are affected by oblique image orientation and most accurate with RotPrior.We showed weighted linear TV to be relatively robust to oblique acquisition; however, Rot-Prior is still maintained to be the most consistently robust method (Figure 9).For all susceptibility calculation methods tested, a unit dipole field misaligned to the main magnetic field (NoRot) leads to artifacts and substantial errors in susceptibility maps.The subtle differences between correction methods found in the numerical phantom ROIs (Figure 6) were not apparent in vivo, probably due to noise, motion, and the expected variability in susceptibility maps over repeated acquisitions. 28,55t nonzero tilt angles, XSIM (and RMSE in Supporting Information Figure S4) values have a respectively high (or low) baseline level arising from rotation (no matter how small the angle) and registration interpolations, and imperfections inherent to in vivo acquisition.Additional discrepancies in similarity measures between rotated and unrotated QSM may also have occurred due to slight differences in repeated acquisitions, as evidenced by the slightly larger (0-5 • ) XSIM discrepancy in vivo than in the numerical phantom (Figures 8 and 9 compared with Figures 5 and 6).The effect of these discrepancies has been minimized by averaging across all 5 healthy volunteers.
Our results also indicate that at larger tilt angles in the numerical phantom, DipK is less accurate than Rot-Prior; therefore, for certain imaging applications including cardiac imaging 56 and pelvic imaging, 57,58 where large tilt angles up to and exceeding 45 • are often required, tilt correction is likely to be essential for accurate susceptibility mapping.
Therefore, we recommend accounting for oblique acquisition by using the RotPrior tilt-correction method

F I G U R E 10
The  maps and difference images illustrating the effects of all tilt-correction schemes on susceptibility calculation in vivo.An axial and a coronal slice are shown for a volume tilted at 45 • and a reference (0 • ) volume with all  maps calculated using the iterative Tikhonov method (top) and weighted linear TV (bottom).Weighted linear TV with DipIm fails at nonzero angles and is therefore omitted from the figure.NoRot leads to the largest differences and image artifacts throughout the brain for iterative Tikhonov and weighted linear TV methods.The EVE ROIs used are shown (bottom left).Results from TKD (Supporting Information Figure S6) are very similar before background field removal, as this method gave the most accurate susceptibility maps in both the numerical phantom and in vivo.If desired, the susceptibility map can be rotated back into the original orientation after susceptibility calculation to facilitate comparison with other (processed) images.In the future, it would be interesting to carry out a similar investigation of tilt-correction methods for total field inversion, 59,60 in which background field removal and susceptibility calculation are combined into a single step.Due to the reliance of total field inversion on the correct definition of the magnetic dipole, we would expect to see similar results, with RotPrior being more robust to oblique acquisition.However, this would need to be confirmed with further work.It is possible to build an alternative rotation-free pipeline of methods relatively unaffected by oblique acquisition (such as LBV and weighted linear TV), but those methods must be checked to ensure true independence of image orientation.However, such an approach limits the choice of methods for the steps in the QSM pipeline, which could lead to suboptimal susceptibility maps.For example, LBV's highly specific boundary approximations can be easily violated, making it easier to simply rotate the field maps in some cases.These aspects must be considered carefully when designing a QSM pipeline.

CONCLUSIONS
Oblique acquisition must be accounted for in the QSM pipeline to avoid artifacts and erroneous susceptibility estimates.We recommend rotating the total field map into alignment with the scanner frame after phase unwrapping but before background field removal (and then rotating the final susceptibility map back into the original orientation).Alternatively, a QSM pipeline relatively robust to oblique acquisition can be built from a more limited number of image orientation-independent methods (eg, LBV or V-SHARP for background field removal and weighted linear TV for susceptibility calculation).However, care must be taken in weighing up the minimal effects of image interpolation (from tilt-correction rotations) versus choosing from a smaller range of methods that are orientation-independent and may not be as robust to oblique acquisition, as they may not be as accurate nor optimal for a given data set.It would also be vital to ensure a chosen method is independent of slice orientation, which may require further investigation.Our recommended correction scheme ensures that all methods developed for each stage of the QSM pipeline can be used and optimized.

Figure S2 .Figure S8 .
Figure S8.Nonoblique k-space magnetic dipole kernels laid side by side to illustrate circular continuity.When there is no oblique acquisition, there are no violations in circular continuity (ie, identical values and no discontinuities at the boundaries of the k-space dipoles; white square)