Joint B 0 and image estimation integrated with model based reconstruction for field map update and distortion correction in prostate diffusion MRI

In prostate Diffusion Weighted MRI, differences in susceptibility values exist at the interface between the prostate and rectal-air. This can result in off-resonance magnetic field leading to geometric distortions including signal stretching and signal pile-up in the reconstructed images. Using a set of EPI data acquired with blip-up and blip-down phase encoding gradient directions, model based reconstruction has recently been proposed that can correct these distortions by using a B 0 field estimated from a separate B 0 scan. However, change in the size of the rectal air region across time can occur that can result in a mismatch of the B 0 field to the EPI scan. Also, the measured B 0 field itself can be erroneous in regions of low Signal to Noise ratio around the prostate rectal air interface. In this work, using a set of single shot EPI data acquired with blip-up and blip-down phase encoding gradient directions, a novel joint model based reconstruction is proposed that can account for changes in the off resonance effects between the B 0 and EPI scans. For ten prostate patients, using a measured B 0 field as an initial B 0 estimate, on a 5-point scale (1–5) image quality scores evaluated by an experienced radiologist, the proposed framework achieved scores of 3.50 ± 0.85 and 3.40 ± 0.51 for b-values of 0 and 500s/mm 2 , respectively compared to 3.40 ± 0.70 and 3.30 ± 0.67 for model based reconstruction. The proposed framework is also capable of estimating a distortion corrected EPI image even without an initial B 0 field estimate in situations where a separate B 0 scan cannot be obtained due to time constraint.


Introduction
Prostate diffusion MRI scans are used for tumour detection within a multi-parametric MRI protocol. However, in some patients the diffusion scans are of limited use due to geometric distortions including signal pile-up and signal drop out at the interface between prostate and rectalair. These distortions occur due to susceptibility differences at the interface creating an additional off-resonance field (called B 0 field) so that the position-frequency relationship is changed. This results in signal stretching in areas within the image where the gradient of the B 0 field has the same polarity as the phase encoding gradient. Conversely, signal compression or pile-up occurs in regions where the B 0 field gradient direction is opposite to that of the phase encoding gradient.
Several techniques have been proposed to correct for both signal pile-up and stretching distortions in diffusion EPI images. These techniques acquire the same slice twice with opposite phase encoding gradient directions, resulting in two data sets, blip-up and blip-down [1][2][3]. Image registration based distortion correction methods [1,2] attempt to correct the distortions by finding the B 0 field as a symmetrical displacement field that will result in identical corrected images from blip-up and blip-down data sets. In areas of low SNR such as the region around the prostate rectal-air interface, these techniques may fail due to registration problems [3][4][5]. Model based reconstruction [6,7] provides an alternative that formulates the distortion correction problem as a model linking the corrupted k-space data to the corrected image. Using a previous estimate of the B 0 field obtained from a separate B 0 scan, this technique solves a full inverse problem using iterative reconstruction and can achieve better reconstruction results in terms of resolution and distortion correction than image registration based methods [7]. However, the model based reconstruction framework assumes the B 0 field to be static except for an offset resulting from scanner frequency drifts. In cases of bowel movements and/or rectal gas changes in the rectal area adjacent to the prostate region between acquisitions of the B 0 and the EPI scans, the estimated B 0 field may be mismatched to the EPI scans. Also, there can be errors in the estimation of B 0 field itself in the low SNR regions [8]. Thus, the model based distortion correction techniques may benefit from an update of the B 0 field estimation as part of the reconstruction problem.
There has been some works on dynamic B 0 field estimation in functional MRI as part of reconstruction process [9][10][11][12]. Sutton at al [10] and Matakos et al. [11] proposed joint image and dynamic B 0 field estimation framework for brain functional MRI that used a set of interleaved gradient echo EPI data sets acquired with different echo times to facilitate the correction of B 0 field errors due to magnetic field drift and respiratory induced phase oscillations. However, as the EPI data sets were acquired with only one phase encoding gradient direction, they did not provide complementary information from opposite phase encoding direction to be able to correct for signal pile-up artifacts [1,13].
In this work, we propose a model based joint image and B 0 field estimation that combines the strength of model based reconstruction with dynamic B 0 field estimation. Using a set of single shot EPI data acquired with blip-up and blip-down phase encoding gradient directions, the proposed framework can account for changes in the off resonance effects between B 0 and EPI scans and can simultaneously compensate for geometric distortions including both signal pile-up and drop out in the reconstructed EPI images. Results with and without an initial measured B 0 field are presented, as well as using the model based reconstruction alone.

Material and methods
The proposed framework acquires two EPI data sets (blip-up and blip-down) with opposite phase encoding gradient directions. For data acquired with a b-value of 0 s/mm 2 , starting from an initial B 0 field, a joint estimation framework is proposed that can estimate both the corrected EPI image and the corrected B 0 field. The corrected B 0 field is then used in subsequent steps of phase correction and model based reconstruction for diffusion weighted data (b-value > 0 s/mm 2 ).
The proposed reconstruction framework has the following three steps: Step 1) Dynamic B 0 field estimation: Starting from an initial B 0 field, this step estimates the corrected B 0 field by iteratively solving a joint image and B 0 field estimation problem. The initial B 0 field can be estimated from a separate dual echo gradient echo scan. In case, a separate scan is not available, the initial B 0 field is set to 0.
With B 0 field inhomogeneities, the corrupted k-space Y j from the j th coil corresponding to trajectory k(t) at time t can be related to the undistorted image x via the following model [14].
where x(r n ) and ΔB 0 (r n ) are image and B 0 field at position r n , r n denotes a position vector corresponding to the n th pixel in the image, C j (j = 1,2, …,J) is the j th coil sensitivity estimated from a separate low resolution SENSE reference scan. N is the number of pixels in the image, k(t). r n denotes the inner product.
For M number of data points in k-space acquired at sampling time points t 1 , t 2 , …, t M , we can expand the above equation in a matrixvector notation as: We can summarize the above equation as: where Y j has dimensions Mx1 and is the acquired distorted k-space from the j th coil, E j (ΔB 0 ) is the measurement matrix with dimensions MxN, x is the undistorted image with dimensions Nx1. The above expression implicitly takes into account any undersampling that is used for accelerated acquisition such as SENSE, Partial Fourier etc. Let k up and k down be the k-space trajectory for blip-up and blip-down scans, the acquired k-space data Y up j and Y down j at b-value of 0 s/mm 2 corresponding to blip-up and blip-down EPI scans for the j th coil are given as: By stacking data and encoding operators from all J coils, we can write: The data from both phase encoding directions can be combined into a single formulation by setting The proposed joint image and B 0 field estimation framework can be summarized as a minimization problem with optimization over both image and B 0 field [11,15]: and R(ΔB 0 ) are quadratic regularization terms ‖Dx‖ 2 and ‖DΔB 0 ‖ 2 respectively, D is the first order finite difference operator that computes derivatives in all in-plane spatial dimensions; β 1 , β 2 are regularization weights that provide a balance between resolution and noise  in the reconstructed image and B 0 field, respectively; ||.|| 2 denotes the l 2 norm. The above joint image and B 0 field problem can be solved iteratively using a two stage alternating minimization scheme [11,15,16] that splits the joint problem into two sub-problems. In the first stage within each iteration, the image update is estimated using a previous field estimate. In the second stage, we minimize over the B 0 field using the estimated image from the first stage. Mathematically, the image update x k in the k th iteration (k = 1,2, …,K) is estimated using a previous field estimate B k 0 1 as: Using estimate x k , the updated B 0 field in the k th iteration B k 0 is estimated as: By using the assumption that the dynamic changes between the current and updated estimates of the B 0 field are small, Eq. (10) can be linearized using first order Taylor series expansion [11] to formulate it in a similar form as Eq. (9). The details of linearization procedure are given in Appendix I. The standard Conjugate Gradient (CG) method can then be used to solve both sub problems in Eqs. (9) and (10).
The convergence of the CG iterations is achieved when either the maximum number of iterations is reached or the normalized residual (||.|| 2 denotes the l 2 norm) in the current iteration becomes smaller than ϵ (ϵ being a small number). The B 0 field in the final K th iteration (ΔB 0,corr ) is then used in the subsequent phase correction and model based reconstruction steps for data acquired at b-value > 0.
Step 2) Phase correction: For diffusion weighted data (b-value > 0 s/ mm 2 ), a phase correction has to be performed because even small physiological motion occurring during diffusion sensitization gradients can cause phase changes that may lead to signal cancellation when data from opposite phase encoding gradient directions is combined, leading to signal drop out in the final reconstructed images.
Using the corrected B 0 field (ΔB 0,corr ) found in Step 1 of the proposed framework, the model based reconstructions x up and x down of the separate blip-up and blip-down phase encoding gradient direction scans are given as: With k-space data Y down as reference, Y up is corrected with the phase correction ΔФ obtained by taking the Hermitian inner product between the two reconstructions x up and x down [7].
Step 3) Model based reconstruction Using the corrected ΔB 0,corr , the phase corrected Diffusion weighted data Y up and Y down from blip-up and blip-down phase encoding gradient directions can be combined into a single formulation Y=E(ΔB 0,corr )x by setting up down up down 0,corr 0,corr 0,corr T T Model based reconstruction using the corrected ΔB 0,corr is performed on k-space data Y by solving the following minimization problem [7]: The above reconstruction from diffusion weighted data is performed for each b-value and diffusion direction.

Experiments
Ten male patients (median weight 83 (range: 68-98) kg and age 73 (57-94) years old were recruited from the clinical prostate imaging pathway and were consented for additional image acquisitions. The study was approved by the local Ethics committee and written signed consent was obtained from all patients for the research scans. Patients were placed feet first into the scanner. No antispasmodic agent was administered.
Scanning was performed on a 3T scanner (Achieva, Philips Healthcare) equipped with 16 anterior +16 posterior channel receive coil array. Single shot EPI data in both blip-up and blip-down phase encoding directions were acquired with a SENSE factor of 2. The EPI scans had the following parameters: resolution = 2 × 2 × 4 mm 3 , FOV = 180-220 × 180-220 × 55-90 mm 3  . For reference, axial T2 weighted images were acquired using a turbo spin echo scan with the following parameters: resolution = 2x2x4 mm 3 , FOV = 180-220 × 180-220 × 55-90 mm 3 , SENSE acceleration factor = 2, TE/TR = 100 ms/4700 ms, scan time = 40 s. For calculation of the B 0 field, a separate 3D dual echo gradient echo scan was acquired at the end of scanning session with the following parameters: resolution = 2 × 2 × 2 mm 3 , FOV = 200-250 × 200-250 × 70-120 mm 3 , flip angle = 6 0 , right to left phase encoding direction, SENSE acceleration factor = 1, echo time difference ΔTE = 2.3 ms, TE1/TE2/TR = 4.6 ms/6.9 ms/8.7 ms, scan time = 1 min. Volume shimming was performed in all scans to cover the whole prostate and surrounding areas. To keep the same shimming from one scan to the next, SAMEPREP option was selected in the scanner software that forces to use the same preparation data for all the scans. In our scans, the patients were asked to lie still and no breathing management was used.

Data post processing
To save the raw data together with the relevant information needed for the reconstruction framework, a software patch was implemented using ReconFrame software (Gyrotools Zurich, CH). EPI phase correction was performed using the ReconFrame tool to correct for ghosts originating from the opposing directions of alternate readouts. Subsequent post processing was implemented in MATLAB. The coil sensitivities were calculated using ReconFrame tool by dividing individual coil images by the body coil images followed by smoothing and extrapolation as proposed in SENSE paper [17]. The measured B 0 field was processed using the quantitative susceptibility mapping toolbox [18] that estimates a B 0 field map by a weighted least squares fit of temporally unwrapped phases in each voxel over echo time. A robust spline based smoothing [19] was applied to the B 0 field map in the image domain to smooth out noisy components. The volumes from B 0 field map and T2-weighted images were resampled to match the EPI scan resolution and the associated field of view. The resampled B 0 field map was further refined by compensating for any frequency offset that may exist due to scanner frequency drift or coordinate shift in scanner software [7].

Joint image and B 0 field estimation
For joint image and B 0 field estimation, we alternated 30 times between updating the image and then updating the field map. In each update of image or B 0 field, 50 sub iterations of CG method were performed. For the CG method, the matrix vector multiplications used for computing the gradients in CG were performed efficiently with time segmentation and FFT algorithms available in Fessler's Image Reconstruction Toolbox [20]. To balance between resolution and noise amplification, the regularization parameters β 1 and β 2 in Eqs. (9) and (10) were chosen to achieve a specific target spatial resolution [21,22], such that the Full Width at Half Maximum (FWHM) of Point Spread Function (PSF) was 1.01 pixels for image estimation and 1.2 pixels for B 0 field estimation, respectively. The computation of PSF was done by using quadratically penalized shift invariant least squares method as implemented in function qpwls_psf.m of Fessler's Image Reconstruction Toolbox [20].
To investigate the dependence of the proposed framework on the initial B 0 field estimate, the proposed joint B 0 and image reconstructions were performed both without an initial estimate (i.e. initial B 0 field as zero) and with an initial estimate set to measured B 0 field. The joint reconstructions were compared against the uncorrected reconstruction without B 0 field (setting ΔB 0 to 0 in Eq. (9) and model based distortion corrected reconstruction using a measured single B 0 field map. The same value for the regularization parameter β 1 was used both for joint and uncoupled reconstructions.
For qualitative assessment, the images reconstructed with different methods were presented in a random order and blinded to the reconstruction method to avoid subjective assessment. A radiologist with 25 years of experience scored each image in terms of overall image quality on a 5-point scale (1:poor, 2: below average, 3: average, 4:above average, 5:excellent) [23].

Results
For patient 1, results for the proposed joint reconstruction without an initial B 0 field estimate are shown in Fig. 1 for data acquired at bvalues of 0 and 500 s/mm 2 . Without distortion correction, the variation in B 0 at the prostate rectal air interface resulted in significant artifacts (caption on next page) and blurring in the reconstructions (Fig. 1b). The proposed framework was able to correct for most of the distortion artifacts in the final reconstructed images even in case of starting with a zero initial B 0 field estimate (Fig. 1c).
For ten patients, the reconstruction results are shown for b-values of 0 and 500 s/mm 2 in Figs. 2 and 3, respectively. For patient 2 and patient 5, significant errors in the measured B 0 field map resulted in imperfect distortion corrections in the model based reconstruction. With erroneous B 0 field map as an initial estimate, the proposed joint reconstruction framework was able to correct for those imperfections. There was no obvious bulk motion observed in any of the patient prostate scans.
A comparison of measured B 0 field maps and the B 0 field maps estimated using the proposed framework is shown in Fig. 4 for selected patients to include both cases a) when measured B 0 is accurate and results in reasonable distortion correction (patient 8), and b) when measured B 0 has significant errors that result in imperfect distortion correction (patient 2 and patient 5) and joint reconstruction helps to correct these errors.
A summary of image quality scores for different reconstruction methods corresponding to b-value of 0 s/mm 2 and b-value of 500 s/ mm 2 are given in Fig. 5a and b, respectively. In the uncorrected reconstruction, the image quality scores were 1.90 ± 0.99 and 1.60 ± 0.51. When no B 0 field map was available, the proposed joint reconstruction was able to achieve significant improvement in image quality scores (3.00 ± 0.67 and 2.90 ± 0.57). When a B 0 field map is available, the joint reconstruction still leads to an improvement in the average image quality score (3.50 ± 0.85 and 3.40 ± 0.51 versus 3.40 ± 0.70 and 3.30 ± 0.67 for model based reconstruction).

Discussion
A novel joint reconstruction framework is proposed that combines the strength of model based reconstruction [7] with dynamic B 0 field estimation [11]. The power of complementary encoding information available from both blip-up and blip-down directions enables the correction of both signal pile-up and signal stretching artifacts within the diffusion EPI images. The joint B 0 field and image estimation method used in our proposed framework allows for compensation of the inaccuracies in the B 0 field due to motion or physiological changes near the prostate-rectal interface between the EPI and B 0 scans. Furthermore, in the case of not having an initial B 0 field estimate, by alternating between the stages of B 0 field update and image update, the proposed method is able to correct the distortions in most of the patients. This is the main advantage of the proposed method as it eliminates the need of a separate B 0 scan. In cases where a B 0 scan is available, starting with the measured B 0 field map, the proposed method performs better than all the other reconstructions. Thus, the proposed framework can be beneficial either when a B 0 scan is performed or in cases where a B 0 scan cannot be performed due to time constraints. Our proposed method estimates the changes in the B 0 field as a single update that occurs between B 0 scan and one set (blip-up and blip-down) of EPI scans. In case of having multiple sets of bvalue = 0 s/mm 2 EPI data acquired in an interleaved manner (for example, in DWI applications requiring long scans due to many b-values, diffusion directions and number of averages, Diffusion Tensor Imaging, functional MRI etc.), the proposed joint reconstruction could be used to estimate corrected B 0 fields and corrected EPI images corresponding to each EPI data set.
The proposed method may be further improved by having slightly different echo times for the two EPI scans such that centre of k-space in the two EPI scans is sampled at two different echo times [10,11]. For gradient echo scans, this time shift has been shown to facilitate better initial B 0 field estimation that can help for better convergence of joint estimation. For the spin echo EPI sequences used in our framework, a time shift may be achieved by shifting the timing of the refocussing pulse using pulse sequence development tools. In our proposed framework, the joint B 0 and image reconstruction was done only for EPI data acquired at b-value of 0 s/mm 2 . Reconstructing images jointly for all bvalues and one B 0 map may result in better reconstruction performance at the expense of increased computational complexity.
The proposed framework used computationally efficient quadratic regularization with first order differences. This results in slightly blurred reconstructions in both image and B 0 field updates. Better regularizers such as l 1 norm based minimization [24] may achieve better results. In our framework, standard CG based optimization was used to solve the sub problems in the joint reconstruction. In future, the CG based optimization might be replaced with alternatives such as quasi Newton methods that ensure the global cost function decreases.
Our proposed reconstruction method solves the full inverse problem constrained to the acquired original raw k-space data rather than the scanner reconstructed images. Image based distortion correction methods (such as Topup method [1]) are designed for practical cases where the imaging community do not have access to k-space data and only the scanner reconstructed distorted images are available. The Topup method reconstructs undistorted images based on image registration based optimization. This may result in an erroneous calculation of the B 0 field attributed to the lack of a unique solution between the corresponding locations in the blip-up and blip-down images especially in regions with severe signal pileup, leading to artifacts or blurring in the final images [3]. Moreover, the Topup reconstructed images may also contain the bias from noise when magnitude images are combined [25][26][27][28]. Our proposed reconstruction avoids this bias due to complex summation/averaging done implicitly via the encoding operators. This is likely to be beneficial for high b-value DWI images with low SNR.
The proposed framework can correct for changes or errors in the B 0 field due to any potential origin such as scanner frequency drift or motion etc. However, it does not correct for physiological motion effects that may occur between the EPI blip-up and blip-down scans. In some cases, modification to the B 0 field might compensate for the motion effects, for example, a rigid translational shift in the PE direction can be compensated by an offset correction in the B 0 field. Integrating the framework with motion field estimation and correction [29,30] may further make the method robust against both distortion and motion corruption effects.
Distortion correction in prostate diffusion MRI is challenging in some patients due to erroneous measured B 0 field that is used to compensate for off-resonance effects. The proposed framework improves the overall image quality by combining the strength of model based reconstruction with dynamic B 0 field estimation. Validation of proposed framework was performed successfully in ten clinical patients. The proposed framework offers a potential to improve the diagnostic value of prostate images for tumour detection in diffusion weighted imaginga technique that is now commonly used for detecting prostate cancer.

Fig. 2.
In-vivo reconstruction results for 10 patients (P1 to P10) for data acquired at b-value of 0 s/mm 2 . The reference T2-weighted image (left column), uncorrected reconstruction (no B 0 map), distortion corrected reconstruction with measured B 0 field map (measured B 0 map), proposed model based joint reconstruction results without initial B 0 estimate (initial B 0 zero) and with initial B 0 estimate set to measured B 0 field map (initial B 0 from map) are shown. Fig. 3. In-vivo reconstruction results for 10 patients (P1 to P10) for Diffusion Weighted data acquired at b-value of 500 s/mm 2 . The reference T2-weighted image (left column), uncorrected reconstruction (no B 0 map), distortion corrected reconstruction with measured B 0 field map (measured B 0 map), proposed model based joint reconstruction results without initial B 0 estimate (initial B 0 zero) and with initial B 0 estimate set to measured B 0 field map (initial B 0 from map) are shown. Fig. 4. Example B 0 fields. The reference T2W image and measured B 0 field map are shown in the first and second columns. The B 0 fields obtained from the joint reconstruction framework without initial B 0 estimate (initial B 0 zero) and with initial B 0 estimate set to measured B 0 field map (initial B 0 from map) are shown in third and fourth columns, respectively. In Patients 2 and 5, the jointly estimated B 0 fields yield improved reconstructions (see Figs. 2 and 3 for details). In Patient 8, the joint reconstruction without initial B 0 estimate (initial B 0 zero) was of inferior quality compared to the other reconstructions.  5. Image quality assessment of proposed framework: Bar plots showing average expert scores for overall image quality (1: poor to 5: excellent) for different reconstruction methods (uncorrected reconstruction (No B 0 map), Model based distortion corrected reconstruction using measured B 0 field map (measured B 0 ), joint reconstruction with initial B 0 field set to zero (initial B 0 zero) and joint reconstruction with initial B 0 field set to measured B 0 map (initial B 0 from map) are shown in (a) and (b) for b-value of 0 s/mm 2 and b-value of 500 s/mm 2 , respectively. The associated standard deviations are also indicated. The proposed joint reconstruction with measured B 0 field map as an initial estimate performed better on average than all other reconstructions.
Similar to Eq. (8), the equivalent joint minimization problem can be expressed as: The above problem can be solved iteratively using two stage alternating minimization scheme.
In the first stage, the image estimate x k in the k th iteration (k = 1,2,..,K) is found by using previous field map estimate = =ω k−1 . Thus, Eq. (A6) reduces to: In the second stage, an updated field map estimate ω k in the k th iteration is found by setting x = x k in Eq. (A6) The above expression is summarized as: Both Eqs. (A7) and (A9) are of similar form and this linear approximation means that both can be solved efficiently using standard Conjugate Gradient methods.