Regularized joint water-fat separation with B0 map estimation in image space for 2D-navigated interleaved EPI based diffusion MRI.

PURPOSE: To develop a new water-fat separation and B0 estimation algorithm to effectively suppress the multiple resonances of fat signal in EPI. This is especially relevant for DWI where fat is often a confounding factor.
METHODS: Water-fat separation based on chemical-shift encoding enables robust fat suppression in routine MRI. However, for EPI the different chemical-shift displacements of the multiple fat resonances along the phase-encoding direction can be problematic for conventional separation algorithms. This work proposes a suitable model approximation for EPI under B0 and fat off-resonance effects, providing a feasible multi-peak water-fat separation algorithm. Simulations were performed to validate the algorithm. In vivo validation was performed in 6 volunteers, acquiring spin-echo EPI images in the leg (B0 homogeneous) and head-neck (B0 inhomogeneous) regions, using a TE-shifted interleaved EPI sequence with/without diffusion sensitization. The results are numerically and statistically compared with voxel-independent water-fat separation and fat saturation techniques to demonstrate the performance of the proposed algorithm.
RESULTS: The reference separation algorithm without the proposed spatial shift correction caused water-fat ambiguities in simulations and in vivo experiments. Some spectrally selective fat saturation approaches also failed to suppress fat in regions with severe B0 inhomogeneities. The proposed algorithm was able to achieve improved fat suppression for DWI data and ADC maps in the head-neck and leg regions.
CONCLUSION: The proposed algorithm shows improved suppression of the multi-peak fat components in multi-shot interleaved EPI applications compared to the conventional fat saturation approaches and separation algorithms.

INTRODUCTION

DWI, which is used to probe the thermal motion of water molecules, plays an important role in many clinical applications. This includes the diagnosis of tumors and stroke and the exploration of structural connectivity. , Single‐shot EPI (ssh‐EPI) readout is often preferred because of its high sampling efficiency and its ability to virtually freeze motion. However, ssh‐EPI with its relatively small bandwidth along the phase‐encoding direction is prone to geometric distortions caused by strong local B0 inhomogeneities. Moreover, the long echo train may lead to substantial blurring, ultimately limiting the spatial resolution and the SNR. To mitigate the above‐mentioned problems, multi‐shot interleaved EPI (msh‐EPI) has become a commonly used technique, which allows uniform k‐space sampling by acquiring multiple interleaves for full k‐space coverage. However, msh‐EPI is particularly sensitive to subject motion when used for DWI. Motion can introduce significant phase inconsistencies between shots, thereby spoiling the overall coherence of the k‐space data. Therefore, a low‐resolution navigator can be added to probe and subsequently correct for physiological motion‐induced phase alterations.

Beyond the brain, DWI can provide valuable information to characterize lesions and determine appropriate treatment strategies in body applications. DWI was reported to be successfully applied to diagnose and monitor the treatments of head and neck squamous cell carcinoma. , , However, residual fat signals reduce the image quality and make improvements of the challenging fat suppression desirable. Although msh‐EPI has an increased bandwidth in the phase‐encoding direction compared to ssh‐EPI, the fat signals (with the dominant CH2 methylene proton peak at 1.3 ppm) still shift significantly with respect to water. When not suppressed, fat can obscure tiny structures in regions with overlapping water–fat signals. Furthermore, because of the multi‐peak nature of the fat spectrum, , , the different lipid signals corresponding to the different chemical sites within the fat molecule are shifted to different locations. Conventional spectrally selective fat saturation methods such as spectral presaturation with inversion recovery/spectral attenuated inversion recovery (SPIR/SPAIR) , have been used to suppress fat signals before sampling. Nevertheless, such spectral selective techniques tend to fail in regions of high B0 inhomogeneity, as encountered in the head‐neck region. , Such failures can even result in an undesired suppression of the water signal, while leaving the fat untouched. The imperfect suppression will affect image quality in most EPI based applications, but especially in DWI, because remaining minor fat resonances can still cause substantial ghosting at higher b‐values because of their neglectable attenuation by diffusion compared to water. , ,

As an alternative to fat saturation, water–fat separation approaches based on chemical‐shift encoding using multiple TEs , , gained growing popularity. Several approaches have been proposed to improve water–fat separation in EPI‐based DWI. Burakiewicz et al proposed an approach to acquire DW images using a chemical‐shift encoded spin‐echo ssh‐EPI sequence with an additional low‐resolution navigator, assuming a single line fat spectrum and ignoring the ssh‐EPI trajectory when disentangling water and fat. To deal with the multi‐peak nature of the fat spectrum in DWI, previous works , attempted to correct for artifacts from those multiple fat peaks by combining fat saturation and chemical‐shift encoding approaches. Recently, Hu et al presented an alternative msh‐EPI approach to remove the fat signals in DWI by using the point‐spread‐function (PSF)‐EPI approach adding an additional encoding dimension to encode chemical‐shift.

In this study, a simpler approach is proposed using a navigator‐based DW msh‐EPI approach (image reconstruction using image‐space sampling function [IRIS]) extended by chemical‐shift encoding to suppress the fat. After reconstruction of the chemical‐shift encoded multi‐shot images, a new regularized water–fat separation algorithm is applied, which was adapted to the EPI sampling process. This algorithm is working in the image domain addressing the chemical‐shift displacements of fat in the phase‐encoding direction. This was done using a series of shift matrices, applied to the common B0 field map, in the joint estimation of B0, water and all fat components. This work is mainly focused on improving image quality for the water images, to especially support diffusion applications. Therefore, chemical‐shift encoding is applied for fat suppression, not risking touching the water magnetization by any magnetization preparation. The algorithm was tested in simulations, providing a ground truth for numeric/statistic comparison, and on healthy volunteer data acquired with a TE‐shifted spin‐echo msh‐EPI with and without diffusion encoding. For comparison, a voxel‐independent iterative decomposition of water and fat with echo asymmetry and least‐squares (IDEAL) algorithm and conventional fat saturation techniques were applied. Numeric/statistical analyses were performed among the 3 different techniques.

METHODS

Sequence

To enable water–fat chemical‐shift encoding, the 2D‐navigated DW msh‐EPI approach was complemented by adding appropriate image echo shift functionality as shown in Figure 1. For a given b‐value, a series of chemical‐shift encoded (eg, N = 3) source data can be acquired, changing slightly the timing of the msh‐EPI train while keeping the timing of the low‐resolution phase navigator fixed.

FIGURE 1

Scheme of the used segmented diffusion EPI sequence allowing for water–fat chemical‐shift encoding and navigator correction. The multi‐shot interleaved (msh)‐EPI sampling window can be shifted back and forth by a time interval () relative to the time of the spin echo () to facilitate chemical‐shift encoding. Please note that the timing of the navigator echo remains fixed to sense potential phase distortions, introduced because of physiological motion during the diffusion sensitizing process. The gray colored gradients, applied on all channels simultaneously, are used for diffusion sensitization

To avoid phase discontinuities because of off‐resonance effects among the individual segments of the basic msh‐EPI acquisition, time shifting is applied using a shift increment τ (duration of 1 EPI read‐out lobe divided by the segmentation factor). This time base is also used to facilitate chemical‐shift encoding by shifting the msh‐EPI block by integer multiples of τ, for simplicity reasons.

Signal model

Because of the low bandwidth in the phase‐encoding direction in EPI, the displacement of fat in this direction can be extremely large. For a given voxel in an EPI image that contains both water and fat, the fat portion is shifted and originates from another voxel/location than the water signal. Therefore, the fat signal in a voxel will also have experienced a different B0 inhomogeneity compared to the water reconstructed within the same voxel. The multi‐peak nature of fat poses extra challenges. Each fat peak has a different spatial shift and therefore, has experienced a different B0 inhomogeneity compared to the other fat signals reconstructed in the same voxel; moreover, each fat peak will appear with a different amplitude in the image. , , , Therefore, the total signal model of 1 voxel can be written as: where  and  denote the complex‐valued water and fat components, [s] the time shift between the readout center of the ‐th TE‐shifted scan and the spin echo [s], [Hz] the B0 map at the original  and shifted location.

Furthermore, and [Hz] give the relative amplitude and chemical‐shift for each fat peak m, whereas  denotes the complex noise. Please note, this signal model is a simplification, it is addressing the fat shift effects in phase encoding direction, but neglects the small chemical‐shift effects in the (odd/even) readout direction with the aim to yield a feasible solution for EPI in image‐space.

Algorithm

In this work, a 2‐step approach is chosen, as in the voxel‐independent IDEAL algorithm, to solve Equation (1) for , and B0: (1) a water–fat separation step, followed by (2) a field map estimation step, is performed in an iterative process. Unless stated otherwise, the initial B0 field map is set to for the first iteration of (1). After the incremental error is calculated using the Gauss‐Newton approach in (2), we update and repeat steps (1) and (2) until the number of iterations reaches the threshold (often ≤ 10 iterations). Because of the low phase encoding bandwidth of EPI, established approaches that ignore spatial shifts , , are no longer suitable to estimate individual water–fat contents and the B0 inhomogeneity induced phase in every voxel. , Therefore, a joint estimation algorithm is implemented that considers the spatial shift both for the B0 map and for the individual fat peaks. These spatial shifts are implemented via shift matrices.

Water–fat separation

Equation (1) can be discretized and written as where S and X are vectorized representations of the N source images and the unknown water–fat components, respectively, that is, and where is the number of voxels. Equation (2) is written into the minimization problem

is the coefficient matrix with and the B0 field at voxel q, such that

For each fat peak , a shifted version of the field map, , can be written as where is a shift matrix that spatially shifts the whole fat image voxels in the phase‐encoding direction. Similar to the image‐domain approach in Hernando et al, the spatial displacement can only take integer numbers. The voxel shift for each fat peak can be calculated from where BW is the bandwidth in phase‐encoding direction [Hz/pixel] and [·] denotes rounding to the nearest integer.

Finally, Equation (5) is written into the linear least squares system and solved with the conjugate gradient (CG) method.

Field map estimation

The field map estimation is based on a Gauss–Newton search method. Similar to Reeder et al, the error terms of each unknown are introduced by writing , , and in Equation (1) as , , and and neglecting second‐order error terms. Using the first‐order Taylor expansion of yields:

Equation (13) can be discretized and written into the minimization problem: where and , and are errors of the unknowns , and . is the current estimation obtained from solving Equation (12). In each iteration, the error map is calculated corresponding to the estimated B0 map and updated to . is the well‐known total variation regularization operator, which can be used to enforce smoothness , , of the updated field map error. Specifically, where and are backward first derivative operators for and directions, respectively.

Furthermore, λ in Equation (14) is the regularization factor that weights the smoothness of the field map and is tuned empirically. The coefficient matrix can easily be built by using the block matrix elements directly from matrix , that is, where

The error for updating the field map can be calculated by solving with CG, using

Because the B0 values are treated as real values in the current implementation, the complex vectors and are split into real and imaginary parts according to Reeder et al as well as all the corresponding matrix elements.

B0 map extrapolation

Because of the large displacement of fat, estimation of the local B0 experienced by a fat component can be compromised at locations where fat is shifted outside the boundaries of the subject. Therefore, an additional extrapolation step was implemented after a few iterations of the algorithm to produce a more stable B0 reinitialization for the remaining iterations. The extrapolation step assumes that the inner area of the object has been shimmed well and is free from large inhomogeneities. This is similar to the region‐growing IDEAL approach, which is initialized by an automatically selected center‐of‐mass seed voxel. In the present approach, a thin‐plate spline extrapolation method similar to Liu et al is applied. An example of this approach can be found in Supporting Information Figure S1. Alternatively, the estimated field map from other methods , , , , would also provide a reasonable initialization (although without consideration of the chemical‐shift effects).

Experiments

Simulation study

Shepp‐Logan water–fat phantom data were simulated in k‐space to generate sets of 3 chemical‐shift encoded source images. One simulated B0 map and the chemical‐shift of each fat peak were both modulated for each k‐space data point with the actual sampling time, including time variations in both phase‐encoding/readout directions. The sampling time map was generated based on the chosen phase‐encoding bandwidth [Hz/pixel] and segmentation factors. Using a 3‐point Dixon, the encoding steps were uniformly distributed over the full encoding circle, matching the dephasing period between methylene peak and water, which corresponds to 2.3 ms at 3T (ΔTEs = 0.24 ms, 1.00 ms, 1.76 ms). The SNR ratio was set to be 100, calculated as the maximal water signal intensities divided by the SD of the noise. Using those data, water–fat separation was performed with the proposed algorithm using 10 iterations.

A voxel‐independent multi‐peak IDEAL algorithm implemented based on Reeder et al and Yu et al, , using a median filter to smooth the estimated B0 map, was applied for comparison. The maximum number of iterations of the IDEAL algorithm was set to 30. This basic voxel‐independent algorithm was used as one of the reference methods, because it works also in image space, and represents the basic approach for water–fat separation. For simplicity, we will refer to this voxel‐independent multi‐peak IDEAL algorithm as IDEAL in all figures.

Numeric simulations were conducted to explore the effect of spatial integer shifts using varying bandwidths from 15 to 50 Hz/pixel (step size, 0.5 Hz/pixel), at a fixed relatively homogeneous B0 field (2D Gaussian profile ranging from −110 to 110 Hz), a segmentation factor of 6, and a 7‐peak fat model.

Two regions of interest (ROIs) were defined, in a mixed region (water–fat overlap) and in a water‐only pure region (no overlap). Overlapping and non‐overlapping regions can be identified by taking the ground truth water–fat distribution, the given EPI bandwidth in the phase encoding direction and the chemical shift of each individual fat peak into account. It should be noticed that, because of the presence of B0‐introduced geometric distortions, the ground‐truth images were also modulated with the same B0 and time maps at each given bandwidth.

To reduce the potential “integer shift” discretization error, a simple zero‐filling interpolation (ZIP) algorithm was implemented. When the shift between water and the dominating fat resonance (methylene peak) corresponds to a fraction as half a pixel (half‐pixel shift), the resulting mismatch because of the “integer shift” character of the image‐domain algorithm, can be mitigated by interpolating the source images to a larger matrix. The separated water images of adopting the ZIP approach are shown in Supporting Information Figure S2. In this simple validation, the half‐pixel effects of all the minor peaks, which have relatively lower amplitude, were not considered.

The normalized root mean squared error (NRMSE) between separated water image and the known ground truth water image were calculated through , where represents the ground truth water, the separated water signal, its average amplitudes and the number of voxels in the water image.

In vivo volunteer study

Experiments were performed involving 8 healthy volunteers with informed consent obtained and approved by the local ethics committee. Interleaved multi‐slice diffusion weighted EPI imaging was performed in the leg and the head/neck regions on a 3T scanner (Philips, Best, The Netherlands), using an 8‐channel knee and a 16‐channel head‐neck array receive coil, respectively. Two basic protocols were applied comprising 3 b‐values (0, 300, 600 s/mm2) with all 3 gradients active in parallel, measuring 4 slices with a gap of 10 mm at a fixed TR of 2s. Therefore, a single diffusion direction was used for DWI. For each b‐value, 3 TE‐shifted source images were acquired. Ramp sampling was used for all measurements and already corrected in the scanner reconstruction. Conventional fat saturation techniques were performed for comparison. In the leg experiments, 6 volunteers were measured with SPIR and 2 were measured with SPAIR. In the head‐neck experiments only SPIR was applied. The main distinction between SPIR and SPAIR is the use of different RF pulses. By using an adiabatic inversion pulse, SPAIR has a better ability to cope with inhomogeneities, at the penalty of increased measuring time or reduced number of slices. The other sequence parameters are reported in Table 1.

TABLE 1

Sequence parameters for data acquisition

Anatomy	Resolution (mm3)	Matrix size	Seg. factor a	TE (ms)	ΔTE (ms)	Bandwidth b (Hz/pixel)	Fat sat.	Acquisition time (s)
Leg	1.4 × 1.5 × 4	160 × 150	6	59	–	36.2	SPIR	51
	1.2 × 1.2 × 4	168 × 162	6	64	–	29.2	SPIR	51
	1.4 × 1.5 × 4	160 × 150	6	59	0.26/1.00/1.74	36.2	–	152
	1.2 × 1.2 × 4	168 × 162	6	64	0.26/1.00/1.74	29.2	–	152
	1.4 × 1.5 × 4	160 × 150	6	59	0.26/1.00/1.74	36.2	SPAIR	152
	1.2 × 1.2 × 4	168 × 162	6	64	0.26/1.00/1.74	29.2	SPAIR	152
Head–neck	1.4 × 1.5 × 4	160 × 150	6	59	–	36.2	SPIR	51
	1.4 × 1.5 × 4	160 × 150	6	59	0.26/1.00/1.74	36.2	–	152
	1.4 × 1.5 × 4	160 × 152	8	53	–	47.5	SPIR	67
	1.4 × 1.5 × 4	160 × 152	8	53	0.16/1.00/1.84	47.5	–	200

seg.factor: the segmentation factor.

bandwidth: the EPI bandwidth in phase‐encoding direction.

DWI image reconstruction for the individual ΔTEs and b‐values was performed off‐line using a Python implementation of the IRIS algorithm. The physiological motion induced phase distortions were sensed by the low‐resolution navigator for each individual shot and were corrected by the IRIS algorithm, producing complex‐valued chemical‐shift encoded source images. No further interpolation of the 3 source images was performed (ZIP), because at the given EPI bandwidths the fraction of the methylene peak pixel shift was below ±0.15 pixels. The coil sensitivity maps were pre‐acquired using gradient echo imaging. , Therefore, they do not match exactly the distortions as encountered in the EPI and navigator images because of chemical‐shift and B0 displacements. Ignoring this mismatch during SENSE reconstruction will result in some residual unfolding and/or ghosting artifacts, even under smooth conditions. In the present work, a simple extrapolation of the coil‐sensitivity maps is used to get a rough sensitivity estimate in areas where the pre‐scan has seen no signal but the EPI did (see Supporting Information Figure S3 for more details).

To improve convergence of the water–fat separation for each slice, the B0 map obtained from b = 0 s/mm2, was chosen as an initialization for the water–fat separation for higher b‐values. This helped to reduce the number of iterations. Therefore, for the in vivo data, the number of iterations was chosen to be 10 for b = 0 s/mm2 data and 5 for higher b‐values. The regularization factor λ was empirically chosen to be 1 for simulations, and 104 for in vivo data. The threshold of the normalized residual norm for convergence of CG was set to 0.001 both for the water–fat separation step and the field map estimation step. The computational time per iteration was 1.5/9.9 s for a single‐peak/multi‐peak fat model in the simulation data (single slice, matrix size 144 × 144), and 17.4/21.6 seconds for a multi‐peak fat model and medium/higher resolution in vivo data (single slice, matrix size 160 × 150/168 × 162). All computations were performed on a Windows 10 computer with an Intel Core i7 CPU (3.0 GHz, 8 cores) and 32 GB of RAM.

Because of the lack of ground truth for the in vivo data, numeric/statistical analyses were conducted by evaluating the ADC fitting. Two ROIs were manually selected as pure water (no water–fat overlap) and mixed (overlap) regions. A total of 200 pixels of each ROI were selected within the same muscle for each slice. Paired t tests were performed between the fitted ADC values of each ROI, comparing the proposed algorithm with IDEAL and SPIR in 6 volunteers. A P‐value < .05 indicated statistical significance.

Fat spectrum self‐calibration

To achieve optimal water–fat separation, a correct fat spectrum model is essential. However, because of T2 relaxation , and potential J‐coupling , effects, Ren’s spectral fat model becomes inaccurate for the TE range used in this work. Therefore, a relative amplitude calibration was performed in a process like the 1 described in Yu et al, taking the fat resonance frequencies as known priors.

In this approach, self‐calibration for the individual fat peaks was done sequentially because of the limited number of TEs. It was assumed that the calibration for each fat peak was not affected by the other peaks. The peaks were sorted and calibrated in descending order according to their weights in Ren’s model. To do calibration for each peak separately, peak‐specific fat masks were constructed in a subcutaneous fat only region of 1 volunteer’s leg. Data for b = 0 s/mm2 were used, neglecting the small water content in fatty tissue. After its definition, this fat mask was shifted along the phase‐encoding direction by the corresponding number of pixels, determined by the peak’s resonance frequency according to Equation (11). To avoid potential biases from restricting shifts to integer steps, an additional boundary erosion step with 1 iteration was performed for each mask, assuring that all pixels inside the mask contain the fat signals of the calibrated peak. The methylene peak was set as the reference fat peak because it is supposed to be the most abundant component in the spectrum. Calibration started with the second peak, combined with the main peak to form a “temporary” 2‐peak fat model. The proposed algorithm was performed while minimizing the residual fat signal leaking into the water channel as a function of the peak’s amplitude. This process was repeated peak by peak until all peaks were included.

The final fat model obtained was = [−485.41, −434.32, −397.27, −341.07, −312.96, −246.54, 77.92] Hz with normalized weights = [0.067, 0.797, 0.000, 0.057, 0.010, 0.009, 0.059] (original Ren’s 7‐peak fat model weights: [0.085, 0.625, 0.071, 0.095, 0.066, 0.016, 0.042]) and was used for all other in vivo reconstructions. The reduced relative amplitude of the fat peaks at −341.07, −312.96, and −246.54 Hz at a given TE may be because of their shorter T2 compared to the remaining peaks (−485.41, −434.32, and 77.92 Hz) , or because of J‐coupling effects. , Zero amplitude was found for the peak at −397.27 Hz, probably because of the difficulty to distinguish it from the dominant methylene peak. The performance of the proposed algorithm with different numbers of peaks in the fat model is shown in Supporting Information Figure S4.

RESULTS

Figure 2 shows water–fat separation results of phantom data sets using IDEAL and the proposed algorithm. The proposed algorithm produces good water images and B0 estimations in the simulation. Although a benign B0 map is simulated, IDEAL shows artifacts in the water and the B0 channel mainly because of the neglection of spatial shifts of fat and the B0 map.

FIGURE 2

Interleaved EPI water–fat separation using simulation example data. One selected interleaved 6‐shots EPI source image (bandwidth of 36 Hz/pixel), simulated with a B0 map (range −110 to 110 Hz) and a 7‐peak fat model. Estimated water, fat, and B0 maps from the proposed and the iterative decomposition of water and fat with echo asymmetry and least‐squares (IDEAL) algorithms are shown. For IDEAL, the omission of spatial shifts results in inaccurate B0 map estimation and insufficient water–fat separation (marked by red arrows), whereas the proposed algorithm shows no distinct artifacts in the water image and B0 map. Furthermore, the absolute difference maps (displayed with 5× magnified scale) show improved accuracy for the proposed algorithm especially in water image and B0 estimation compared to IDEAL. However, because of the simplification of the EPI signal model, for fat the errors in the fat image are larger than those for water (marked by the blue arrow)

Figure 3 shows the NRMSE plots of the simulation. The NRMSEs of the proposed algorithm are maintained at reasonable values close to the noise level (~0.019) in both regions. However, a slight oscillation in the error train can be seen for both mixed and pure ROIs. The maximum error amplitude always appears at those bandwidths where the dominant fat (methylene) peak is shifted to a “half pixel.” To mitigate errors resulting from the half‐pixel shift effects, interpolation of the image to a larger matrix size was proposed to ensure an “integer shift” for the methylene peak; this indeed did reduce the enhanced NRMSE (see Supporting Information Figure S2).

FIGURE 3

Normalized root mean squared errors (NRMSE) for the water images as a function of the EPI bandwidth obtained in simulations. NRMSE between separated water and ground truth water image shown for mixed/pure regions of interest (ROIs), respectively. The source data is simulated with a 7‐peak fat spectrum model and varying bandwidth. NRMSE in both regions of the proposed algorithm always remains around the noise level, whereas a slight error increase can be seen at bandwidths where the methylene peak is almost shifted to a “half pixel” (eg, at 45.5 Hz/pixel). In comparison, IDEAL shows larger errors, especially in the mixed ROI

Figure 4 shows decomposed water and fat images of 1 subject’s leg and 1 subject’s head‐neck for both IDEAL and the proposed algorithm for a non‐diffusion case. Neglecting the spatial shifts in the signal model, IDEAL results in an inaccurate water–fat decomposition. Moreover, in the head‐neck slice, water–fat swap artifacts appear in regions of severe B0 inhomogeneities. The proposed algorithm produces significantly improved results in both cases.

FIGURE 4

Comparison between the iterative decomposition of water and fat with echo asymmetry and least‐squares (IDEAL) and the proposed algorithm for non‐diffusion sensitized (b = 0 s/mm2) leg and head‐neck data. (A) One source image of the reconstructed leg data and corresponding water–fat separation results from the proposed and IDEAL algorithms. IDEAL leaves rim‐like artifacts (marked by red arrows) in the fat overlapping region in the B0 map and water image. This artifact is not present when reconstructed with the proposed approach. Furthermore, some residual fat ghosting can be seen in both the source and the separated fat images of both algorithms because of the neglected chemical‐shift effects in readout direction (indicated by the blue arrow). (B) One source image from the head‐neck region and the water–fat separation results calculated through the proposed and voxel‐independent algorithms. Significant errors in B0 estimation cause artifacts in the corresponding regions of the water image for the IDEAL algorithm (marked by red arrows). With the help of the TV regularization and extrapolation used in the proposed algorithm, the B0 errors and species swaps in the water–fat images from severe B0 deviations appearing in the results of the IDEAL algorithms are reduced

Figure 5 shows 4 slices of 2 subject’s legs at b = 0 and 300 s/mm2, comparing the proposed, voxel‐independent IDEAL, and fat saturation method. The proposed algorithm provides more reliable water images in all slices for non‐diffusion and diffusion cases in contrast to IDEAL. For SPIR (volunteer 1), although B0 is rather homogeneous in the leg region, the inhomogeneity of the body coil transmission is another compromising factor, , contributing to incomplete fat suppression (marked by the blue arrow) because of wave propagation effects, which is not the case for SPAIR (volunteer 2). The image quality matches well between the proposed method and SPAIR especially regarding the methylene peak, whereas the fat signal contributed from the olefinic peak (5.31 ppm) can be further removed through the proposed approach (marked by the white arrows).

FIGURE 5

Comparison of 3 different approaches in 2 subject’s legs. Four slices of water‐only images selected from 2 volunteers at 2 b‐values are shown. For the fat saturation data, spectral presaturation with inversion recovery (SPIR) and spectral attenuated inversion recovery (SPAIR) were applied in volunteer 1 and 2, respectively. The proposed algorithm produces good water–fat decompositions both with and without diffusion gradients. The fat‐related artifacts as shown in Figure 4 appear in the iterative decomposition of water and fat with echo asymmetry and least‐squares (IDEAL) images as well (marked by red arrows). For spectral presaturation with inversion recovery (SPIR) in volunteer 1, an artifact resulting from inhomogeneity can be seen (marked by the blue arrow). In volunteer 2, spectral attenuated inversion recovery (SPAIR) avoids this artifact. Two regions of interest (ROIs) are selected (marked by the white squares) and used for all approaches, illustrated at the bottom with adjusted level/window. With the help of shift matrices, the fat components from different fat peaks can be suppressed simultaneously, whereas other methods show residual fat signals in the water‐only images (eg, olefinic peak, marked by the white arrows)

Figure 6 shows results of 1 subject’s head–neck DWI at 2 b‐values comparing SPIR, IDEAL, and the proposed method. SPIR and IDEAL suffer from large field inhomogeneities in the outer image areas because of bad B0 shimming conditions, showing remaining artifacts in the water‐only images. The proposed algorithm can resolve the ambiguities from chemical‐shifted fat and produce better water–fat separation.

FIGURE 6

Comparison of the 3 different approaches in 1 subject’s head‐neck region. Four slices with 2 b‐values are displayed. Because of the severe field inhomogeneities in this case, the fat signal can neither be completely suppressed by spectral presaturation with inversion recovery (SPIR), nor separated by iterative decomposition of water and fat with echo asymmetry and least‐squares (IDEAL) (marked by the red arrows). The proposed algorithm produces significantly improved water–fat separation, especially in regions where IDEAL generates water–fat swaps

Figure 7 shows 2 example ADC maps for the 3 fat suppression techniques both for the leg and for the head‐neck data. The proposed algorithm produces better ADC quantifications in both anatomies, compared to IDEAL and SPIR with the previously described artifacts.

FIGURE 7

ADC maps of 3 different approaches in 1 subject’s leg and 1 subject’s shoulder. Source images (non‐diffusion) with the corresponding ADC maps are shown, comparing the proposed algorithm, iterative decomposition of water and fat with echo asymmetry and least‐squares (IDEAL), and spectral presaturation with inversion recovery (SPIR). Two regions of interest (ROIs) were manually selected in the source images as mixed (water–fat overlap) and pure water (no overlap) regions. In the leg data, IDEAL and SPIR show fat‐related artifacts in the mixed region of the ADC map (marked by the blue arrow), which does not appear in the proposed approach’s results. For the head‐neck region, both IDEAL and SPIR suffer from severe B0 inhomogeneities, contributing to artifacts in the ADC measurements (marked by blue arrows). In contrast, the proposed algorithm produces smooth ADC maps in the same regions

Figure 8 shows the natural logarithm of the water signal intensities varying as a function of the b‐values. In the leg anatomies (Figure 8A), IDEAL shows inconsistent log‐signal decays between the 2 chosen ROIs. The proposed algorithm and SPIR adhere more to the expected linear decrease in the logarithmic plots although SPIR was slightly affected by inhomogeneities (marked by the blue arrow in Figure 7). In the head‐neck regions (Figure 8B), the critical B0 inhomogeneities contribute to species swap in the IDEAL water images and incomplete fat suppression for SPIR (as shown in Figure 6), breaking the log‐linearity for both ROIs.

FIGURE 8

ADC fitting assessment across different volunteers. Natural logarithm plots of average water signal in mixed/pure regions for the leg and the head‐neck regions in 6 volunteers (3 volunteers/anatomy) as a function of b‐value. For each anatomy, 12 slices (4 slices/volunteer, 3 volunteers) with data in mixed/pure regions are included in the averaging and standard deviation calculations. (A) Leg data. The logarithmic signal decays show reasonable linearity in both regions of interest (ROIs) for the proposed method. In contrast, iterative decomposition of water and fat with echo asymmetry and least‐squares (IDEAL) shows larger error bars in the mixed regions and significant signal inconsistencies between the 2 ROIs. Spectral presaturation with inversion recovery (SPIR) shows a slight difference of logarithmic signal decays between the 2 ROIs, mainly because of the strong inhomogeneity. (B) Head‐neck data. The problematic B0 inhomogeneity in the head‐neck region spoils the log‐linearity of the water signal evolution for IDEAL and SPIR in both regions (mixed/pure water). In comparison, the proposed algorithm produces good consistency of the logarithmic signal decays and relatively small deviations in both mixed/pure regions improving ADC fitting consistency

The improved performance of the proposed algorithm was also confirmed by calculating the apparent diffusion coefficients. For both anatomies and each technique, the average ADC values, and the P‐values from the t tests between pure/mixed ROIs of water ADC, are displayed in Table 2. The proposed algorithm yields accordant ADC values for the 2 ROIs with no statistically significant difference between them, whereas IDEAL and SPIR did show significant differences in water ADC between the 2 ROIs for both anatomies.

TABLE 2

Quantitative assessment of ADC (×10−3 mm2/s) among the 3 methods

	Proposed	IDEAL	SPIR
Leg
Mixed	1.62 ± 0.26	1.93 ± 0.49	1.74 ± 0.51
Pure	1.62 ± 0.22	1.61 ± 0.23	1.58 ± 0.32
P‐value a	0.498	0.001	0.001
Head–neck
Mixed	1.59 ± 0.45	1.09 ± 1.31	1.44 ± 0.75
Pure	1.58 ± 0.38	1.32 ± 1.03	1.69 ± 0.61
P‐value	0.218	0.001	0.001

P‐values are derived from group‐wise comparison using a paired t test between ADC values of mixed and pure regions for each technique.

DISCUSSION

In this work, we proposed an alternative water–fat separation algorithm for chemical‐shift encoded EPI. This approach is using a multi‐peak fat spectrum signal model and a novel regularized algorithm to jointly estimate water, fat and a single B0 map using chemical‐shift encoded EPI data. This algorithm deals with the finite bandwidth problem in the phase‐encoding direction of EPI. The proposed algorithm can remove most of the present fat signals as demonstrated for simulated and in vivo data. In vivo validation has been done for interleaved spin‐echo EPI and for more challenging interleaved DWI applications. In all measurements the proposed algorithm produced proper water images, even under challenging ΔB0 conditions. It should be noted that, based on the signal model of Equation (1), the proposed algorithm is also applicable to other EPI‐based approaches.

There are 2 difficulties for obtaining accurate water–fat separation for EPI: (1) the chemical‐shift displacements of individual fat peaks; and (2) the field map estimation in regions with severe field inhomogeneities. The proposed algorithm can resolve them and provides water–fat decomposition for both diffusion and non‐diffusion cases. Nevertheless, in cases of severe B0 inhomogeneities, the effects of these 2 features are hard to distinguish. This is one of the shortcomings of our validation approach. However, in the relatively homogeneous B0 experiments (eg, in Figure 2) and in the leg measurements, the benefit of correcting the spatial spectral shifts can be demonstrated independent of B0 inhomogeneity effects.

For data acquisition, the TE‐shifted DW spin‐echo msh‐EPI sequence was implemented to enable chemical‐shift encoding and reduce geometric distortions. It should be noted that, because of the intrinsic properties of EPI, both the acquisition and reconstruction were based on some assumptions. One limitation of the present work is that alternating chemical‐shift displacements of fat in the odd and even “echoes” of the EPI train were ignored in the signal model. This issue results in residual ghosting artifacts in the fat images (see Figures 2 and 4). However, those ghosts are chemical‐shift encoded, and hence, only appear dominantly in the final fat image. Therefore, such artifacts do not noticeably affect the quality of the final water images and the associated ADC measurements.

Another limitation of the proposed image‐domain‐based algorithm is the integer shift correction of the fat signal for each individual chemical shifted line, applied in the separation procedure. The small error (Figure 3), can either be mitigated by choosing an appropriate EPI bandwidth during acquisition or by image interpolation during post‐processing to realize close‐to‐integer shifts for the dominant methylene fat peak signal in the image data (see Supporting Information Figure S2). Such an interpolation step can potentially also be directly included in the shift matrices to avoid increasing the problem size.

An alternative option would be to solve the water–fat separation problem using a k‐space‐based approach. , , , In the approach of Brodsky et al, the actual k‐space sampling time for each point is used in the separation to correct for chemical‐shift displacement effects in both the phase‐ and frequency‐encoded directions. However, the B0 map is still estimated and demodulated in image space, assuming a smooth B0 field with neglectable differences between neighboring regions. This method and improved similar approaches show great ability in correcting the relatively minor chemical‐shift displacements of fat in non‐EPI acquisitions. , , However, in EPI, the increased chemical‐shift displacements in the phase‐encoding direction, may not be neglected as in other sequences. One weakness of the present work is that we compare the proposed approach to only 1 rather basic water–fat separation algorithm. Further comparisons to different approaches assessing the impact of the spatial signal displacement for fat in EPI can be considered as interesting future work. In the work of Honorato et al, the B0 estimation and demodulation, as well as water–fat separation, were performed in k‐space. Compared to the proposed algorithm, solving the whole problem in k‐space would theoretically correct the ghosting in the fat image originating from the alternating chemical‐shift in the readout direction. Moreover, it would at the same time correct susceptibility‐introduced geometric distortions. However, Honorato et al has only been implemented for 1D spin warp trajectories with already long computation times. A complete 2D solution is pending and would need further investigation, but is expected to be too slow for clinical applications.

The present work is focused on multi‐shot EPI acquisitions, which typically have better SNR compared to single‐shot EPI based approaches. , , Because of the increased bandwidth along the phase‐encoding direction, msh‐EPI shows strongly reduced geometric distortions. , Unlike the work of Hu et al, in which the geometric distortions and chemical‐shift displacements are corrected by estimating an appropriate k‐space kernel including the PSF‐dimension, the proposed algorithm only starts the water–fat separation after the 3 complex source images were reconstructed. This simplification also limits its scope to correct geometric distortions in advance. Nevertheless, the estimated B0 map could be used to further correct the geometric distortions by established post‐processing approaches , , as shown in Supporting Information Figure S5.

Similar B0 demodulation methods in k‐space , , , , can also be considered in combination with k‐space based water–fat separation. , The necessary B0 map can be estimated in the image space through the proposed algorithm in the first place. With the help of regularization and shift matrices, the main chemical‐shift effects in the phase encoding direction have already been corrected, and the estimated smooth B0 map is more reliable. Therefore, the geometric distortions and chemical‐shift effects in both directions can be addressed well in this full‐model‐based hybrid approach, and is a step to be explored in future work.

Apart from the acquisition and reconstruction, another potential source of inconsistency is the variability of the fat spectrum, also depending on the actual sequence parameters, such as the sequence type, TE, and TR. In the present work, spin‐echo EPI was exclusively used in a small TE range around 60 ms, which is why invariance of the self‐calibrated fat spectrum was assumed for all the sequences used. More accurate calibration of the spectrum can be used, for example, using single‐voxel stimulated echo acquisition mode (STEAM) spectroscopy‐based pre‐calibration , , or self‐calibration approaches using more echoes.

The DWI application shown in this work is based on the low‐resolution navigator phase information delivered for the IRIS reconstruction. This navigator is not chemical‐shift encoded and contains some fat signal and related SENSE unfolding errors. To mitigate the latter, an extrapolation of the coil sensitivity map has been applied (described in detail in Supporting Information Figure S3). Please note that this is only a simple approximation. GRAPPA could also be a potential candidate for mitigation but needs further investigation. The navigator phase information is also slightly compromised by the remaining fat that is spatially shifted. Approaches to entirely suppress the fat signals in the navigator, for example by using gradient reversal approaches, should be the subject of future research.

With respect to the DWI application, this study presents a feasibility study for navigated water–fat separated DWI using interleaved EPI. For clinical application, more diffusion orientations should be acquired to produce reliable trace ADC quantifications. Considering the large number of images acquired in a clinical setting, further improvements to the computational speed are necessary. To achieve higher efficiency, already estimated B0 maps from any b‐value, diffusion orientation, or neighboring slice could serve as a better initialization, which would lead to fewer iterations. Undersampling in the multi‐shot and chemical‐shift encoding space are promising approaches for future research when SNR is sufficiently high. Other EPI based applications like DTI or intravoxel incoherent motion (IVIM) imaging approaches are also considered for future investigations.

CONCLUSION

Applying chemical‐shift encoding to msh‐EPI can significantly improve the fat suppression and image quality of DWI. The proposed algorithm provides a feasible solution for water–fat separation with EPI and has been validated using synthetic data from simulations and acquired data from healthy volunteers. This algorithm can also be applied to other EPI‐based works to provide reliable fat‐free images.

CONFLICT OF INTEREST

Peter Börnert is an employee of Philips research.

FIGURE S1 The thin‐plate spline extrapolation step during the field map estimation (the total number of iterations is 10). (A) The “current” field map estimated in the proposed water–fat separation algorithm through 5 iterations after starting with 0 initialization throughout FOV. The algorithm may fall into local minima near the object edges (marked by the red arrow), contributing to inaccurate B0 values. (B) A binary mask created through a chosen signal threshold of source images and an erosion step to avoid the boundary effect of the object’s edge. The number of erosion iterations is empirically chosen to be 22 for the image matrix size of 160 × 150. (C) Resulting B0 map initialization for the residual iterations through the extrapolation of the B0 map. A function in MATLAB (The MathWorks, Natick, MA) called “tpaps” is used to perform the 2D thin‐plate spline extrapolation. The extrapolation is applied once after the 5th iteration, based on the masked inner voxels, which are assumed to have correct B0 values. A total of 2000 random control pixels are chosen within the mask and the smoothing factor p is empirically set to be 0.0002. For the remaining iterations, the field map is initialized with (1) estimated B0 values of the previous iterations inside the mask; and with the (2) extrapolated B0 values outside the mask. (D) The final B0 field map can be obtained after finishing the residual iterations

FIGURE S2 The extended simulation validation of the zero‐filling interpolation (ZIP). The phantom source images are simulated in k‐space as described in the Methods section, with matrix size 144 × 144, a B0 map ranging between −110 and 110 Hz and 7‐peak fat spectrum model. (A) Separated water images of the simulation at different bandwidths. The spatial shifts of the methylene peak for each simulation are illustrated. The images are divided into 2 groups: (1) “no half‐pixel” effect of methylene fat peak, with bandwidths = 36.5, 39.5, 42.5 Hz and the corresponding spatial shifts of 11.90, 11.00, 10.10 pixels, and (2) “with half‐pixel” effect of methylene fat peak, bandwidths = 35, 38, 45.5 Hz and the corresponding spatial shifts 12.41, 11.43, 9.54 pixels. The water images are estimated from the proposed algorithm. The absolute difference maps are calculated between the separated water and the ground truth water images. For the group “with half‐pixel” spatial shift of the methylene peak, as described in the Method section, the ZIP is performed to increase the matrix size of the 3 chemical‐shift encoded source images, to ensure that the spatial shift of the methylene peak in phase‐encoding direction is close to an integer. In this work, for a given bandwidth per pixel in the phase‐encoding direction, the source images were retrospectively interpolated to a matrix size, which ensures the methylene peak having a spatial shift close to the nearest integer (minimum fraction of less 0.1 pixels accepted). The interpolated target matrix size is chosen to be 150 × 150, 152 × 152, 150 × 150 for the given bandwidths, with corresponding spatial shift of 12.93, 12.06, 9.94 pixel of the methylene peak. To calculate the absolute difference maps, the ground truth water images are also interpolated to the corresponding matrix size. Compared to the “no half‐pixel” group, some artifacts can be seen in the water images and absolute difference maps of the “with half‐pixel” group. With the help of interpolation, artifacts can be reduced to an acceptable level, producing better water images that are consistent with ground truth water images. (B) Normalized root mean squared error (NRMSE) between separated water image and ground truth water image with/without ZIP. A fraction threshold of 0.1 pixel was empirically chosen to decide whether ZIP is applied for each simulation (eg, ZIP is performed when the fraction of the total pixel shift of methylene fat is smaller than 0.1). A slight error increase can be seen at bandwidths where the methylene peak is shifted to close to a “half pixel” (the blue curves). After interpolating the image to a larger matrix size, which ensure an “integer shift” for the methylene peak, the errors could significantly be reduced to an acceptable level (the green curves)

FIGURE S3 Validation of 2 extrapolation methods for the sensitivity map and the associated effect on the image/navigator quality for a subject’s leg (b = 300 s/mm2) data set. (A) Coil‐sensitivity map of 1 selected receive element, (B) reconstructed source image, (C) navigator magnitude image, and (D) phase image are shown with no coil sensitivity extrapolation applied. The next rows (E‐H) and (I‐L) show corresponding data with nearest‐neighbor extrapolation and thin‐plate spline extrapolation applied to the coil sensitivities. The raw coil‐sensitivity map (A) was acquired via gradient echo imaging , where chemical‐shift displacements of fat along the phase‐encoding direction are not present. Therefore, fat signals that are shifted outside the object in both EPI/navigator images are insufficiently unfolded because of missing sensitivity information. Therefore, part of the fat signal outside the subject is missing (marked by red arrows) and fat‐related artifacts (marked by blue arrows) are folding into the EPI/navigator images (C‐D). As a result, the reconstructed image quality is degraded because of imperfect SENSE unfolding and impaired navigator phase correction (B). Coil sensitivity extrapolation from inside the object to the full FOV can help to mitigate this effect to a certain extent. This is a similar task compared to the extrapolation of the B0 map (Supporting Information Figure S1). There are many different methods proposed to do the extrapolation for the coil‐sensitivity maps. , , , , In this work, we implemented a simple fix with “nearest‐neighbor extrapolation” to find voxel values outside the subject with the value of nearby translated voxel values under the Voronoi diagram principle, shown in (E‐H). Moreover, the extrapolation based on thin‐plate spline functions proposed by Liu et al was re‐implemented and produced similar results illustrated in (I‐L). Both extrapolations are helpful to mitigate most of the associated artifacts in EPI/navigator reconstructions, resulting from chemical shift and ∆B0

FIGURE S4 Investigation of the performance of the proposed algorithm with different numbers of fat peaks included. (A) One b = 300 s/mm2 DWI source image of 1 selected subject’s leg, and the associated separated water images calculated from the proposed algorithm with a 1‐peak model and 7‐peak model, respectively. A small region of interest (ROI) is selected and displayed in the adjusted level/window of each water image, to show the effect of the proposed algorithm to remove small peaks of fat signal. (B) The normalized root mean squared error (NRMSE) plot with different numbers of fat peaks included in the model. The fat peaks are arranged in descending order with their abundance in the multi‐peak spectrum model (eg, 1‐peak means a single‐peak model with the dominating methylene fat peak at 1.3 ppm). It should be noticed that in the 7‐peak spectrum, 1 peak at 1.59 ppm is detected with no composition. The water‐only ROIs for each slice are selected by the threshold‐based binary masks created from the separated water image with multi‐peak spectrum at b = 0 s/mm2. The NRMSE is calculated between the separated water images with different numbers of fat peaks included in the models. NRMSE in 4 slices of 1 volunteer’s leg data at b = 300 s/mm2 were shown. The NRMSE curves decay with an increasing number of fat peaks, indicating that the fat signals of the different fat peaks are effectively removed from the water‐only image. In principle, for simplicity and computational efficiency, a 4‐peak fat spectrum can be considered, because the 4‐peak separation yielded a small NRMSE around 0.01 compared to the 7‐peak spectrum separation

FIGURE S5 Validation of the post‐processing geometric distortion correction. Non‐diffusion water images of 1 volunteer’s leg and 1 volunteer’s head‐neck region reconstructed by the proposed algorithm are shown. With the estimated B0 map, a further distortion correction step based on Munger et al was performed for the separated water images. The comparisons are made among the original water images without distortion correction, the water images with distortion correction and the reference gradient echo mDixon image. As shown in Figure 4, the B0 field is relatively homogeneous in the leg region and is severely inhomogeneous in head‐neck region. Benefiting from the multi‐shot acquisition, the susceptibility introduced geometric distortion is relatively small in the leg. Compared to the reference image, a small geometric mismatching can be seen in the blue arrow marked region, whereas the distortion correction can get rid of it. In the head‐neck region, under certain critical B0 field conditions, the geometric distortions are relatively large, causing signal pile‐ups in the corresponding regions. With the help of estimated B0 map, this post‐processing step allows to reduce the geometric distortions, producing a water image that better matches the reference image

Click here for additional data file.