Improving chemical shift encoding-based water-fat separation based on a detailed consideration of magnetic field contributions.

INTRODUCTION

Chemical shift encoding‐based water–fat imaging (WFI) has a wide range of clinical applications including fat suppression1, 2 and fat quantification.3, 4 In WFI, magnetic field distortions are modeled as a nonlinear “field map” parameter in a complex water–fat signal model.5 If the field map term is known, the linear water and fat parameters can be determined unambiguously. However, if the field map is unknown and there is no reasonable a priori estimate of the field map, solving the water–fat separation problem can result in ambiguous solutions often appearing as the infamous water–fat swaps—regions in the WFI output images, where the water and fat signal seem to be interchanged.1, 6

To resolve water–fat swaps in WFI, there are many algorithms available that rely on the assumption that the field map term resembling the magnetic field in the MR scanner during a scan is a spatially smooth function.7, 8, 9, 10, 11, 12, 13, 14 Given good enough field map estimates in the neighborhood of a voxel (seeds), the assumption of an only mildly varying field map term over voxels allows to limit the range of possible field map values in that voxel, which can mitigate the ambiguity problem of the water–fat separation.1 However, in cases where the assumption of field map smoothness over the voxel range is not met or reliable seed voxels are not available, water–fat separation remains challenging.

A recent work proposed the use of object‐based information of the magnetic field to obtain a valuable field map initialization for standard WFI algorithms.15 In this work, field distortions due to the object in the scanner are simulated starting from a rough estimate of the susceptibility map of the object in the scanner, which is a binary tissue–air mask of the FOV assigned corresponding susceptibilities. A demodulation of the so‐called object‐based fast field map estimate (OBFFME) from the complex source images served the role of the initialization of several standard water–fat separation algorithms from the ISMRM water–fat toolbox16 and was shown to resolve many water–fat swaps in challenging anatomies in which air–tissue interfaces create field map terms fast varying in space. However, the air–tissue susceptibility map obtained in the OBFFME method by thresholding only holds information of the object inside the FOV, which results in a missing residual field after the demodulation of the OBFFME. Furthermore, results were shown primarily at 1.5 T and without considering the effect of other field map contributions.17, 18

The purpose of the present work is to develop a field map demodulation method for WFI considering—in addition to the previously proposed OBFFME15—the magnetic field contribution from the missing residual field, inhomogeneities of the scanner magnet and the shimming gradients, all at 3 T.

THEORY

Assuming the widely used single‐ multi‐fat‐peak water–fat voxel signal model the complex signal at the n‐th echo is composed of the complex signal of water W and fat F that share a common transverse relaxation rate . The fat signal is modeled by an a priori known spectrum with P spectral peaks of corresponding relative amplitudes a and chemical shifts . 19 The spatially dependent field map, , accounts for the averaged magnetic field changes in a voxel affecting the effective precession frequency of each spectral peak and as such has a strong nonlinear effect on the phase of the signal . We model four different field map contributions: the scanner magnet inhomogeneities , the shim field , the object‐based field map estimate and a residual field B r. γ is the proton's gyromagnetic ratio.

Contribution of scanner magnet inhomogeneities

The volume inside a MR scanner, where the main magnetic field itself is relatively constant, is limited to a region around the iso‐center. In scans with large FOV or at large offcenter locations the main magnetic field is not constant in the whole imaging volume. Based on the magnet design of the scanner, the inhomogeneities of the main magnetic field can be described by a spherical harmonic expansion: with the spherical harmonic functions of order and degree m with corresponding coefficients up to a maximum order L. These coefficients typically do not change from scan to scan and are determined by the hardware and its calibration.

Contribution of shim field

When an object is placed in the main magnetic field, depending on its susceptibility the object gets magnetized and becomes itself a source of magnetic field.20 To counter the magnetic field distortions due to the susceptibility‐induced object‐based field , shimming is routinely used in clinical MR scans. After a 1D or 2D sequence module for measuring the object‐based magnetic field prior to the actual scan, electrical currents through special shim coils are computed to best possibly compensate for . These coils are designed to create certain spherical harmonic magnetic fields and are therefore also typically described by a spherical harmonic expansion like Equation (3) with different coefficients . While clinical systems up to 3 T often only have coils that can produce shim fields up to order L = 2, high‐field systems have shim coils producing fields up to order L = 3 or higher. Since the magnetic field in the scanner changes from scan to scan, the coefficients can be newly calculated and set before each scan. Spherical harmonics in Cartesian form are typically used for shimming (see Appendix).

Contribution of object‐based field map

The magnetic field distortions that the shim field tries to compensate for are mainly created by the object in the scanner. Given a susceptibility map χ of the whole object and assuming , the resulting susceptibility‐induced field map distortions can be approximated in the forward model with the dipole kernel which at the origin is set to a “DC offset” of 0.21, 22, 23 A susceptibility map of the whole object is typically not available and information about it can normally only be obtained from a scan with a FOV smaller than the object itself. According to Equation (4), the object's total susceptibility‐induced field is, therefore, separated in the so‐called object‐based field map estimate and a residual field B r created by the object's susceptibility distribution outside the FOV24:

METHODS

In four separate steps, the theoretical contributions to the field map (2) are estimated and demodulated from the complex source data prior to standard WFI algorithms. Figure 1 gives a flowchart overview of the proposed method, where each step is described in detail in the following sections.

Figure 1

Flowchart of the proposed method: four three‐dimensional field map contributions (second row) are estimated from the coefficients of spherical harmonic expansions and the complex source images (magnitude and phase): magnet inhomogeneities , shim field , and object‐based field map (computed through an object‐based fast field map estimation15) and a residual linear field B r (computed through a linear field estimation). Each field map contribution is demodulated (dashed arrows) from the phase before the input into a chemical shift encoding‐based water–fat imaging (WFI) algorithm to separate water and fat images

Scanner magnet inhomogeneities term

The information about the inhomogeneities of the scanner magnet was directly available in the scanner software as the coefficients of the spherical harmonic expansion in Equation (3) in the scanner x–y–z system. To be able to demodulate , it was then transformed to matrix size, voxel size, and orientation of the complex source images.

Shim field term

The shim coefficients up to order L = 2 were calculated by the scanner and available through the scanner software similar to the coefficients , but given in Cartesian scanner coordinates equal to the conventions in Ref. 25 given in the Appendix. The same transformation as for needed to be applied before demodulation.

Object‐based field map term

An object‐based field map estimate was obtained by the method described in Ref. 15. A maximum intensity projection across echoes was thresholded at 5% of its maximum value. The resulting binary tissue mask was assigned a susceptibility value of ppm, the average of the reported susceptibilities of water, ppm, and fat, ppm, inside the tissue region and 0.34 ppm for the outside air region.26, 27 From the obtained crude two‐component susceptibility map, the OBFFME was forward simulated with Equation (4) by zero‐padding the susceptibility by a factor of 2 in all dimensions, solving the convolution in Fourier‐space and removing the padding again. Lastly, the object‐based field map was centered around its mean value inside the tissue region to not alter the center resonance frequency of the scan data once the object‐based field map was demodulated.15

Residual field term

The yielded OBFFME is able to detect strong field changes near air–tissue interfaces inside the FOV. Background fields from tissue susceptibility sources outside the FOV cannot be incorporated in the OBFFME and are therefore still encoded in the phase evolution of the demodulated signal . Through computer simulations it has recently been shown that forward simulations of the object‐based field map with a FOV smaller than the object tend to result in a linear residual field.28 We therefore model the resulting residual field as a spatially linear field contribution with slope and location . Due to the Fourier shifting theorem, this linear field corresponds to a constant shift in the Fourier domain of the n‐th echo source image, . After the demodulation of the other field contributions the source images were bilaterally zero‐padded in the direction where the FOV clips the object before applying a 3D fast Fourier transformation. Figure 4 gives an overviewing flowchart of the linear field estimation step. The distance of the maximum of the k‐space magnitude with respect to the volumetric center of k‐space corresponds to a shift for each echo n. With the assumption of a linear residual field as the origin of that linear k‐space shift, the measured should then be linear across echo times t

Figure 4

Flow chart illustrating proposed residual linear field estimation. After the demodulation of all other field map contributions, the complex multi‐echo source data is first zero‐padded in the dimension where the FOV crops the imaging object. Second, each echo is Fourier‐transformed and the offset of the voxel of maximum k‐space energy and the k‐space center is measured. The fitted slope of measured offsets versus echo time yields the a in Equation 6

Due to experimental inaccuracies in the estimation of , possibly due to phase errors in the signal s, it is important to stabilize the estimated linear phase contribution across echo times instead of directly demodulating the measured . Anything but a perfect linear time dependence of would interfere with the assumed signal model (1). Therefore, we linearly fitted the measured against the echo times t as , which directly holds the slope in Equation (5). According to Equation (6) the intercept should be zero. Generally, we observed that was well described by a linear time dependence, however possible phase errors especially in the first echo, which is more prone to those errors,29 can lead to worse fittings with a considerably nonzero intercept . This was countered by replacing the first echo point by the origin in the ‐t plane and fitting only with and , which yields fit parameters closer to the assumption (6). The possible content of the intercept is discussed below.

To investigate the validity of the residual field being linear in the case of the present in vivo scans, a numerical simulation based on the Duke phantom was performed similar to Ref. 30. The whole body mesh data were converted into a three‐dimensional volume corresponding to a susceptibility map of 1.5 mm isotropic resolution. Three tissue regions together with the air region were assigned magnetic susceptibility values of −11.31 ppm for bone,26 −9.04 ppm for soft tissue/water,31 0.40 ppm for air,31 and 4.32 ppm for lung.32 The above susceptibility values were referenced to air, so the outside region had 0 ppm susceptibility. The object's susceptibility map is shown in the first column of Figure 2. The corresponding field map in the cervical region was forward simulated as described above via Equation (4) with four different model sizes in the head–feet direction (HF0–HF3). Line profiles in head–feet direction through the middle of the center slice for each cropped model size (HF1–HF3) were compared to the same profile of the whole body field map (HF0) and their differences were linearly fitted in the head–feet region. The same simulation was repeated for the right leg of the Duke phantom (see Figure 3).

Figure 2

Object‐based fast field map estimates of a numerical whole body susceptibility phantom based on Duke30: Starting from a susceptibility map χ (Column 1) with values stated at the end of Section 3, four object‐based field maps (Column 2, 3) are forward simulated via Equation 4, each with a different head–feet coverage (HF0–HF3). Visual comparisons inside a FOV equal in size to HF1 shows different results based on the head–feet coverage of the initial susceptibility map. Line plots show that the difference (blue) between the field map of the whole body simulation (yellow) and the cropped coverage (red) can be approximated by a linear function, indicated by a linear fit (black) in the FOV HF1

Figure 3

Similar to Figure 2: Simulation of the object‐based fast field map estimation in the ankle with the numerical whole body susceptibility phantom based on Duke.30 The object‐based fast field map estimate is computed for three different head–feet coverages (HF1–HF3) of only the right leg and related to the whole body simulation. Line plots show that the difference (blue) between the field map of the whole body simulation (yellow) and the cropped coverage (red). The residual linear field slope is indicated by a linear fit (black) in the FOV HF1. Note how the slope of the linear fit of the residual field decreases for larger head–feet coverages

In vivo measurements

To test the proposed method in an anatomy difficult for WFI due to concave geometries and many air–tissue interfaces, we scanned the cervical region of 10 healthy volunteers (6 female, 4 male, average age [ years]) in a 3 T scanner (Ingenia, Philips, Release 5.1.8, Best, The Netherlands); Informed written consent by all volunteers and approval by the institutional review board (Klinikum rechts der Isar, Technical University of Munich, Munich, Germany) was granted beforehand. All scans used a multi‐echo gradient echo sequence with flyback gradients (monopolar) acquiring three echoes with one shot in a single TR, with  = 1.06 ms,  = 1.6 ms, flip angle = 3°, orientation = coronal, readout direction = feet–head, FOV =  mm3, acquisition voxel size =  mm3, bandwidth/pixel = 1924.2 Hz, scan time = 1:08.7 minutes. Each cervical region per subject was scanned twice—once without any shim gradients and once with the default shim option switched on, where the shim volume was covering the whole FOV. This resulted in a total of 20 datasets to test the reproducibility of the proposed method across subjects.

Spherical harmonic expansion coefficients of the scanner magnet inhomogeneities and the shim field were automatically written to the raw data by modifications of the Gyrotools Recon Frame scanner patch33 and later extracted from the raw data with the MRecon software,33 which was also used to reconstruct the raw data and transform the field contributions to the same coordinate system.

To estimate the effect of concomitant gradients and k‐space misalignments due to gradient delays, we performed one coronal scan of the cervical region with the same sequence as described above, but with the implemented scanner option to calculate the concomitant gradients online as well as an additional measurement of first order echo misalignments. The echo misalignments were measured by enabling the repetition of the monopolar gradient readout through the k‐space center with opposite polarity and without any phase encoding gradients similar to the method described in Ref. 34. The same scan was repeated axially to cover the case where the readout dimension is anterior–posterior, different from the body axis along which the FOV clips the imaging object. For comparison the axial dataset and all estimated field map contributions where reformatted to the coronal perspective.

The proposed method was also tested in a different anatomy also typically challenging for WFI. The left ankle of one volunteer was scanned with the same sequence, also with two different frequency encoding directions, along the feet–head axis and the anterior–posterior axis.

Demodulation of field map terms

Demodulation of all field map contributions is done by which results in a complex signal modeled as where the term ψ is now small compared to f in Equation (1), so that the smoothness assumption of the WFI algorithm is supported.

Water–fat separation

Water–fat separation was performed with two previously published methods with MATLAB implementations available in the ISMRM water–fat toolbox16: a hierarchical IDEAL11 and a graph cut algorithm.6 The default settings of most algorithmic parameters as provided in the toolbox were kept. Only the range of field map values was adapted for the graph cut algorithm to allow values between [−1000 Hz, 1000 Hz] in steps of 4 Hz. For the assumed fat spectrum, a 10‐peaks fat model was used as previously measured in bone‐marrow.35

Three different post‐processing schemes were performed: (1) standard water–fat separation with the raw source signal as direct input to the WFI algorithms without any demodulations, (2) the previously published method15 with demodulation of the object‐based field map estimate as the only preprocessing step, and (3) the proposed method with all four demodulation steps.

Based on the water and fat images of the WFI results of all volunteer scans in the cervical region, datasets were visually rated and categorized in datasets with and without water–fat swaps, where regions in the heart and only a few voxels close to the object boundaries were not regarded. The overall counts of datasets with and without residual water–fat swaps served as a more global metric on how the proposed method performed in a challenging anatomy across all 20 datasets.

RESULTS

Numerical simulation results

Columns 2 and 3 of Figure 2 show the results of the Duke phantom simulation. The line profiles in Column 4 show that the difference in the field maps (red) of a whole body simulation (HF0, yellow) to a smaller size simulation (HF1–HF3, light blue) are for the most part linear in the head–feet direction, which is also confirmed by the linear fit of the difference (black) inside the FOV (HF1 in Column 4). As more coverage in head–feet direction is used for the simulation, the slope of the linear fit vanishes, which indicates that the residual field is created by the cropped objects regions outside the FOV. The same effect can be seen from the results of the numerical simulation of Duke's leg shown in Figure 3.

In vivo WFI results

The phase of all echoes after stepwise demodulation of each estimated field map contribution gradually becomes more homogeneous with less phase wraps, which improves the intermediate WFI results until swap‐free water–fat separation is achieved, as shown in Figure 5. The demodulation of the scanner magnet inhomogeneities removes the characteristic “ripple” phase wraps near the edge of the FOV (red arrows). More phase wraps are clearly removed by demodulating the shim field (blue arrow). In the concave neck regions it is easily observable how the demodulation of the object‐based field map spatially flattens the phase (orange arrows). After the OBFFME demodulation step, the number of phase wraps almost perpendicular to the head–feet direction increases with echo time (green arrow). This indicates the residual linear background field contribution, which is removed from the phase evolution by the linear field demodulation step, resulting in homogeneous phase images in space and time, which furthermore allows for swap‐free water–fat separation.

Figure 5

Example of stepwise demodulation of field map contributions. First row shows scanner magnet inhomogeneities , shim field , object‐based field map , residual linear field B r. Row 2 shows phase images of the third echo at . Row 3 shows fat images of intermediate hIDEAL11 results. Each column shows the phase and fat images after the additional demodulation of the field map contribution in row 1. Arrows indicate specific features of field map contributions that are removed from the phase images by the demodulation step. Extended version, Supporting Information Figure S1

To determine the performance of the proposed method over a larger number of acquired volunteer scans, Table 1 sums up the cervical datasets that show water–fat swaps somewhere in the imaging volume disregarding the heart and the very edge voxels of the object. Numbers in parenthesis do not count the characteristic ripples at the edge of the FOV due to the scanner magnet inhomogeneities. None of the WFI algorithms were able to reconstruct swap‐free results on their own. While the previous method15 is able to reduce only few swaps if no shimming is used, the proposed method leads to a 100% reduction of water–fat swaps for the hIDEAL algorithm and a significant reduction from 20 to 5 datasets suffering swaps for the graph cut algorithm. In these nonresolved graph cut datasets, fat suppression was still good; remaining swaps would not reduce clinical evaluation as only a few slices showed a few swapped voxels on a thin fringe line, which might be resolvable by further optimizing the parameters of the water–fat separation algorithm.

Table 1

Number of cervical datasets showing water–fat swaps (numbers in parenthesis do not count the characteristic ripples at the edge of the FOV due to the scanner magnet inhomogeneities)

	hIDEAL			Graph cut
	Total	Shim off	Shim on	Total	Shim off	Shim on
std. WFI	20 (16)	10 (6)	10 (10)	20 (10)	10 (7)	10 (3)
previous	20 (14)	10 (4)	10 (10)	20 (11)	10 (5)	10 (6)
proposed	0	0	0	5	2	3

Compared are the results of two algorithms—hIDEAL11 and graph cut6—in two scans with and without shimming, where three different processing schemes were performed: Standard water–fat imaging (std. WFI) without any demodulation steps, the previously investigated method15 with one demodulation steps (previous), and the proposed method with four incorporated demodulation steps (proposed).

Figure 6 shows the water images of two example datasets with shim off (left) and shim on (right), where only the proposed method was able to successfully resolve all water–fat swaps. Arrows point to regions, where the additional field map contributions could remove characteristic swaps. Although some datasets reconstructed by both algorithms without (“std. WFI”) and with (“previous”) the OBFFME demodulation showed good water–fat separation, all datasets suffered from high frequency water–fat swaps at the same location where the scanner magnet inhomogeneities cause ripple phase wraps. These ripple like water–fat swaps were resolved in all cases by the demodulation of the scanner magnet inhomogeneities.

Figure 6

Comparison of the proposed method to previous methods without (“std. WFI”) and with (“previous”) a single demodulation step of the object‐based field map estimate before two WFI methods.6, 11 Water images of the WFI results, shown for one case with shim off (left two columns) and one case with shim on (right two columns). Only the proposed method solved the water–fat separation correctly resulting in swap‐free water–fat separation in the shown datasets. WFI results for the two algorithms can be subject to different scaling as the images are direct output of the implementations

Figures 7 and 8 show the effect of two remaining field map contributions—concomitant gradients and gradient delay echo misalignments. In the same scheme as Figure 5, the evolving phase images and the intermediate hIDEAL WFI results are shown columnwise after demodulation of each field map contribution. To compare the order of magnitude of the effect of each field contribution on the phase, the two terms and are displayed, where is the measured k‐space shift due to gradient delays, is the concomitant gradient field, and ω 0 is the center frequency. Column 1 starts with the phase images after demodulation of the scanner magnet inhomogeneities, the shim field and the OBFFME. While Figure 7 shows the coronal scan with the frequency–phase–slice encoding coordinates indicated in the upper left, Figure 8 shows the coronal reformatting of the dataset from the same volunteer which was scanned axially and therefore has a different frequency–phase–slice encoding coordinate system. While in the former case (Figure 7) the frequency‐encoding direction is parallel to the body axis on which the object is clipped by the FOV, in the latter case (Figure 8) the two directions are perpendicular. By comparing Figures 7 and 8, one can see that the residual field that was assumed to be linear along the axis of the object region outside the FOV is not only much larger than the concomitant gradients and gradient delay contributions but also separable when the clip axis and the frequency‐encoding axis is not the same.

Figure 7

Example of stepwise demodulation of two additional field map contributions—concomitant gradients and gradient delays—in coronal scans of the cervical region. Row 1 displays the field map contributions. Columns display phase at the latest echo TE 3 and fat images from intermediate hIDEAL11 results after demodulating scanner magnet inhomogeneities , shim field , object‐based field map (Column 1) and after additionally demodulating concomitant gradients (Column 2), gradient delays (Column 3), and residual linear field B r (Column 4). Note that the residual linear field B r is an order of magnitude greater than the phase term resulting from the measured k‐space shift due to gradient delays. ω 0 is the center frequency. Here, the dimension in which the FOV crops the object coincides with the frequency dimension. The phase term is two orders of magnitude smaller than the estimated residual linear field. Extended version, Supporting Information Figure S2

Figure 8

Similar to Figure 7: Stepwise demodulation of field contributions in the dataset from the same volunteer scan with different frequency–phase–slice encoding coordinate system. The two additional field map contributions—concomitant gradients and gradient delays—remain much smaller than the estimated residual linear field B r. Here the frequency encoding dimension is perpendicular to the dimension of the cropped object indicating that the residual field B r cannot be explained by timing errors in the readout. Extended version, Supporting Information Figure S3

Figure 9 similarly shows the stepwise results of the proposed method for the shim‐off ankle scan with alternating frequency‐ and phase‐encoding directions. Due to the smaller FOV and the geometry of the anatomy, the magnet inhomogeneities are smaller and ripple‐creating field variation are only present in air regions inside the FOV, which is why the magnet inhomogeneity column is not displayed. In this scenario without shimming, the WFI algorithms cannot resolve the swap in the leg as illustrated by the graph cut results in rows 4 and 5. The presence of phase wraps along the body axis indicates a linear residual field contribution after the demodulation of the object‐based field map estimate. After the proposed demodulation of the estimated residual linear field, the graph cut algorithm yields swap‐free water–fat images. This is true in both subplots of Figure 9, where the frequency and phase encoding directions are switched, which shows that the residual linear field remains varying along the head–feet axis and is independent of encoding directions.

Figure 9

Stepwise demodulation in the ankle dataset with shim off for alternated frequency and phase encoding direction. Left: frequency encoding direction along feet–head. Right: frequency encoding direction along anterior–posterior. Due to the smaller FOV the magnet inhomogeneities are negligible in the tissue regions (not shown). Demodulation of the object‐based fast field map estimate leaves a linear residual field component. After demodulation of the estimated linear field, the graph cut algorithm is able to resolve all swaps as shown in row 3. Extended version, Supporting Information Figure S4

DISCUSSION

The focus of the present study was to develop a method to improve WFI in challenging anatomical regions by reducing water–fat swaps via a better field map initialization of existing algorithms. Compared to the previous work,15 three additional contributions to the field map were taken into account: inhomogeneities of the scanner magnet, the shim field, and a susceptibility‐induced field from object field sources outside the FOV.

As clinical scanners only have shim coils that can produce spherical harmonic fields up to a certain order, the shim field will never totally compensate for the object‐based field. In Ref. 15, it was therefore recommended to not perform shimming and only rely on the demodulation of the object‐based field map estimate. In the presented in vivo measurements without shimming, WFI with solely demodulating the OBFFME could only resolve swaps in 2 out of 10 datasets compared to WFI without any demodulation.

The characteristic ripples at the edge of the large‐FOV phase images come from high orders in the spherical harmonic expansion of the scanner magnet inhomogeneities, which have very high field map values compared to the other contributions (Figure 5, top row). Such strong scanner magnet inhomogeneities can lead to increased intra‐voxel dephasing and consequently signal loss, which was observed in some datasets at the top of the skull. Having the shim options on and later demodulating it, proved successful in avoiding some of such intravoxel dephasing‐induced signal losses (not shown).

The coarse two‐component object‐based field map estimate is targeted to simulate the large field map variations originating from the air–tissue interfaces inside the FOV and faces three main challenges. First, the OBFFME is only applicable in 3D‐tomography and not in 2D imaging techniques where there are a limited number of slices. Second, the thresholding process in the OBFFME to obtain the initial crude two‐component susceptibility map falsely assigns signal voids due to short T 2/ bone regions to also have the susceptibility of air. This leads to only very few incorrectly estimated voxels in the brain and did not significantly diminish the quality of the water–fat separation. In anatomies where no air is located inside the tissue, for example, in the extremities, it is easy to assign the susceptibility of bone to any signal voids inside the object and simulate the OBFFME starting with a three‐component susceptibility map (air, tissue, bone). Third, neglecting susceptibility sources outside the FOV proved to be of much higher importance than the mis‐assignment of susceptibility in air/bone regions.

The phase images after the third demodulation step showed a missing linear field map contribution and k‐space shifts increasing with the echo time. The forward simulation based on the Duke four‐component susceptibility map showed the presence of a residual field term originating from object regions that lie outside of the FOV. The relation of the extent of the clipped object regions to the object region inside the FOV influences the nonlinearity (and monotony) of the residual field term along the axis where the FOV clips the objects, most often the body axis. This heuristic result was also obtained in a previous work.28 Here, we approximated the residual field by only its first order, which allowed to estimate its presence by a shift in k‐space, which linearly depends on echo time.

Instead of demodulating the OBFFME and then measuring a residual linear field, one could also try to interpolate the susceptibility map from the thresholding step in the OBFFME method beyond the FOV, for example, by using a replicative padding of the edge voxels instead of zero‐padding in the forward simulation step. While this can work in scenarios where the object regions outside the FOV are fairly cylindrically shaped, in general mismatch between true and interpolated geometry lead to the introduction of more artificial field contributions.

Once the k‐space shifts are measured, a robust fit (6) of the ‐t slope needs to be performed. We found that neglecting the shift of the first echo and instead using the point where the shift should be zero resulted in the best performance of the currently tested WFI algorithms. Neglecting the first echo shift improves the fit as the phase image of the first echo is known to be more prone to phase errors.29 With the assumption of a linear background field one can also subsequently deduce that at the interpolated zero echo time, no shift in k‐space should be present. In our gradient echo sequence with flyback gradients, there are also other possible origins of linear phase shifts, such as gradient delays in the hardware or eddy currents, which can lead to the nonzero intercept of the linear –t fit.34, 36 Therefore, if the frequency‐encoding direction is parallel to the axis on which the FOV clips the object, demodulation of linear field might also remove contributions from other possible origins of linear phase shifts. The cervical and ankle scans with switched frequency–phase–slice encoding axis (Figures 7, 8, 9) show that the linear phase rolls due to concomitant gradients and gradient delays are present and separable from the estimated residual linear field from object regions outside the FOV. The residual linear field varies along the axis on which the FOV clips the objects and the orientations of concomitant gradients and linear field ramps due to gradient delays depend on the frequency‐encoding direction. However, in scans in which both axes coincide (Figures 7 and 9) the concomitant gradient and the gradient delay contributions were at least an order of magnitude smaller than the others. This was also confirmed by one dominant component in the measured in Equation (6), precisely the one corresponding to the axis where the object is clipped.

The additional demodulation of the concomitant gradients and artificial field map variations due to echo misalignments caused by gradient delays showed how the proposed method is easily extensible to also include more field map contributions. In scans that aim at measuring quantitative parameters in the water–fat signal model with more echoes often in several interleaves, echo‐wise concomitant gradients as well as echo misalignments are more important than in the here demonstrated qualitative WFI application.34, 37

While the present study shows the benefit of incorporating the described field map contributions as a priori knowledge to the water–fat separation problem, it has several limitations. First, the results shown here focused on the WFI performance with respect to fat suppression and used a three echo sequence to scan the cervical anatomy in a limited number of volunteers. However, the cervical region is one of the most challenging areas and the double scanning of each volunteer with two shim options gave 20 datasets that can be regarded as worst case scenarios of a priori unknown field distortions, in which, due to the large FOV, WFI was tested in a high number of slices. The successful applicability of the proposed method could also be shown in a total of four ankle scans—in combinations of coronal and axial scans with and without shim.

Second, only two algorithms from the ISMRM water–fat toolbox were used to compare the different demodulation/initialization schemes. Both represent different classes of algorithms. The hIDEAL algorithm solves the water–fat problem voxel‐wise at multiple scales, whereas the graph cut algorithm searches an optimal cut through a graph consisting of all voxels as the graph's nodes. The available implementations gave reasonable computation times to perform the different schemes for all acquired datasets. The graph cut needed approximately up to 30 minutes per data set, the hIDEAL algorithm only around 30 seconds per data set without significant time dependence on the demodulated field contributions. Other available implementations of different algorithms in the toolbox were not studied in the context of this work, either because they do not assume strictly the same signal model (1), they try to separate water and fat in the Fourier domain or they were much slower to be practical in comparisons of many datasets.

In experimental scenarios where phase information is of interest, for example, in susceptibility‐weighted imaging or quantitative susceptibility mapping, one is interested in small phase changes originating from differences in tissue susceptibilities much smaller than between air and tissue. Our proposed method of demodulating most of the largest field contributions might be viewed as a way of separating those large phase contributions from the fine “local” phase differences, resulting in a local tissue field map after a WFI algorithm with a nonsmoothed quantitative field map output. Total field inversion38 on such a field map might also prove interesting as the presently proposed demodulation step also resembles a possible preconditioning.

As the proposed method starts from complex source data and scanner output information independent of the MR pulse sequence, the technique can therefore also easily be translated to different sequences such as, for example, bipolar multi‐gradient echo or single‐ or multi‐shot pulse sequences, whose treatment is outside of the scope of this work. Furthermore, the information about the field map contribution estimated here is potentially valuable in MR application that use field map estimates for correction purposes such as, for example, k‐space WFI methods, echo planar imaging or deblurring in nonCartesian imaging.

Example MATLAB code to estimate and demodulate the considered magnetic field contributions is freely available for download at https://github.com/maxdiefenbach/MRI_field_contributions.git (SHA‐1 = 11d95a993fc5c70bc20ebf2b2107dbe3a3886534). Extended versions of Figures 5, 7, 8, and 9 are available online as supplementary material (Supporting Information Figures S1–S4).

CONCLUSION

We proposed a methodological framework for good field map initialization for WFI through modeling and demodulating four major field map contributions, inhomogeneities of the scanner magnet, the shim field, an object‐based field estimate and a residual field, from the complex multi‐echo signal prior to standard algorithms. The proposed method resulted in almost swap‐free WFI results and performed significantly better than the stand‐alone algorithms with and without demodulation of only the object‐based field map estimate. The complex signal after all four demodulations contains interesting phase information possibly valuable for quantitative parameter estimation methods.

CONFLICT OF INTEREST

Dr. Holger Eggers and Dr. Jakob Meineke are employees of Philips Healthcare.

Additional Supporting Information may be found online in the supporting information tab for this article.

Figure S1 Example of stepwise demodulation of field map contributions. First row shows scanner magnet inhomogeneities , shim field , object‐based field map , residual linear field B r. Rows 2–4 show phase images of –. Rows 5 and 6 show intermediate hIDEAL result water and fat images, respectively. Each column shows the phase and water–fat images after the additional demodulation of the field map contribution in row 1.

Figure S2 Example of stepwise demodulation of two additional field map contributions—concomitant gradients and gradient delays—in coronal scans of the cervical region. Each column shows the phase and water–fat images after the additional demodulation of the field map contribution in row 1. Rows 2–4 show phase images of – after demodulating scanner magnet inhomogeneities , shim field , object‐based field map (Column 1) and demodulating concomitant gradients (Column 2), gradient delays (Column 3), and residual linear field B r (Column 4). Rows 5 and 6 show intermediate hIDEAL result water and fat images, respectively.

Figure S3 Stepwise demodulation of field contributions in the dataset from the same volunteer scan as in Figure 11 but with different frequency–phase–slice encoding coordinate system. The two additional field map contributions—concomitant gradients and gradient delays—remain much smaller than the estimated residual linear field B r. Here the frequency encoding dimension is perpendicular to the dimension of the cropped object indicating that the residual field B r cannot be explained by timing errors in the readout.

Figure S4 Stepwise demodulation in the ankle dataset with shim off for alternated frequency and phase encoding direction. Left: frequency encoding direction along feet–head. Right: frequency encoding direction along anterior–posterior. Due to the smaller FOV the magnet inhomogeneities are negligible in the tissue regions (not shown). Demodulation of the object‐based fast field map estimate leaves a linear residual field component. After demodulation of the estimated linear field, the graph cut algorithm is able to resolve all swaps as shown in row 5 and 6.

Click here for additional data file.