Model‐based reconstruction framework for correction of signal pile‐up and geometric distortions in prostate diffusion MRI

Purpose Prostate diffusion‐weighted MRI scans can suffer from geometric distortions, signal pileup, and signal dropout attributed to differences in tissue susceptibility values at the interface between the prostate and rectal air. The aim of this work is to present and validate a novel model based reconstruction method that can correct for these distortions. Methods In regions of severe signal pileup, standard techniques for distortion correction have difficulty recovering the underlying true signal. Furthermore, because of drifts and inaccuracies in the determination of center frequency, echo planar imaging (EPI) scans can be shifted in the phase‐encoding direction. In this work, using a B0 field map and a set of EPI data acquired with blip‐up and blip‐down phase encoding gradients, we model the distortion correction problem linking the distortion‐free image to the acquired raw corrupted k‐space data and solve it in a manner analogous to the sensitivity encoding method. Both a quantitative and qualitative assessment of the proposed method is performed in vivo in 10 patients. Results Without distortion correction, mean Dice similarity scores between a reference T2W and the uncorrected EPI images were 0.64 and 0.60 for b‐values of 0 and 500 s/mm2, respectively. Compared to the Topup (distortion correction method commonly used for neuro imaging), the proposed method achieved Dice scores (0.87 and 0.85 versus 0.82 and 0.80) and better qualitative results in patients where signal pileup was present because of high rectal gas residue. Conclusion Model‐based reconstruction can be used for distortion correction in prostate diffusion MRI.


| INTRODUCTION
Prostate cancer is the most commonly diagnosed cancer in males and the second cause of cancer-related deaths in men. 1 Early detection of prostate cancer can help in better management and treatment of disease when the cancer is still localized. Multiparametric prostate magnetic resonance imaging (mpMRI) is now becoming a common tool for early detection and staging of prostate cancer. Typically, mpMRI is done as a combination of T2-weighted (T2W), diffusion-weighted MRI (DWI), and dynamic contrast-enhanced MRI (DCE) scans, which determine the likelihood of clinically significant cancer at a particular location within prostate. 2,3 DWI is the main sequence for cancer detection in the peripheral zone of the prostate, 4 where 75% of tumors usually occur. 5 Prostate cancers normally show abnormal diffusion restrictions and high signal in DWI images.
In DWI, single-shot echo planar imaging (SS-EPI) has been the preferred k-space read out technique because of its fast speed and robustness against motion artefacts. However, SS-EPI has a disadvantage of low bandwidth along the phaseencoding direction, making it sensitive to inhomogeneities in the magnetic field. In MRI, spatial encoding is achieved by using magnetic field gradients, and the measured signal represents the Fourier transform of the object. In the presence of magnetic field inhomogeneities, an additional off-resonance field (called B 0 field) will add an additional component to the linear magnetic field gradient so that the position-frequency relationship is changed. This results in signal stretching in areas within the image where the gradient of B 0 field has the same polarity as the phase-encoding gradient. Conversely, signal compression or pileup occurs in regions where the B 0 field gradient direction is opposite to that of the phaseencoding gradient. This results in multiple pixels along the phase-encoding direction being merged into a single pixel. The main source of B 0 field inhomogeneities is the susceptibility differences that arise at the interface between tissues of different susceptibilities. Larger local susceptibility differences will create a stronger off-resonance field that will result in more-severe signal pileup or stretching within the image. For DWI prostate images, susceptibility differences result in severe distortions and poor depiction of zonal anatomy. 6 This can result in nondiagnostic images, errors in image interpretation, or a requirement for biopsy. Additionally, DWI images from prostate cancer patients with metallic hip prostheses 7 may also suffer from susceptibility-induced distortions across the whole prostate region. In addition to susceptibilityinduced distortions, a translational shift may occur along the phase-encoding direction 8 in the reconstructed EPI images attributed to drifts in center frequency of the magnet.
Several methods have been proposed to address the B 0 distortion problem in EPI images. Some of these techniques correct distortions directly in the image space using a B 0 field calculated from a separate dual-echo gradient echo scan. [9][10][11][12] Having an accurate B 0 field, these methods can generally correct for signal-stretching artefacts in the images. However, using EPI data acquired with only 1 phase-encoding direction may not provide sufficient information to correct for signal pile-up artefacts. To address this issue, image-registrationbased reverse-phase-encoded gradient methods have been proposed for brain functional MRI that acquire the same slice twice with opposite phase-encoding gradient directions, resulting in 2 data sets, blip-up and blip-down. 8,13 Using the assumption that the signal pile-up effects with blip-up data set will correspond to the signal-stretching effects in the same region with blip-down data set and vice versa, these methods attempt to find the B 0 field as a symmetrical displacement field by using image registration that will result in identical corrected images in the 2 phase-encoding gradient directions. However, instead of solving the full inverse problem constrained to the acquired raw k-space data, these methods perform image-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), 14 leading to artifacts or blurring in the final images. Furthermore, in image-based methods, the registration becomes more difficult when the image signal-to-noise ratio (SNR) is poor 15,16 and registration errors may lead to failure in distortion correction. By defining a forward model linking the distortion free image to the acquired corrupted raw k-space data, an exact solution of the full inverse problem may be achieved by using the forward and conjugate transpose operators with a conjugate gradient scheme. 17 Distortions attributed to B 0 field inhomogeneities may be corrected directly during the reconstruction process using this full inverse problem solution, although a previous estimation of the B 0 field is needed. The complex averaging used in the reconstruction avoids the bias from noise that can occur when magnitude images are combined. [18][19][20][21] This is likely to be especially beneficial for high b-value DWI images with low SNR. 18,20,21 In this work, we propose a model-based reconstruction framework to solve the full inverse problem of prostate DWI distortion correction. By using the EPI raw k-space data from acquisitions with blip-up and blip-down phase-encoding gradients, we model the distortion correction problem linking the corrupted k-space data to the corrected image and solve it in a manner analogous to sensitivity encoding (SENSE), 22 using the conjugate gradient (CG) iterative method. To solve the issue of translational shifts attributed to the uncertainties in center frequency, we include the centre frequency offset correction as an initial unknown in our framework. In addition, a phase correction is also incorporated into our framework to avoid any phase cancellation issues that may arise because of small motion of the tissue within the diffusion-encoding gradients. Results using the proposed method are compared in 10 patients against a neuroimaging distortion correction method based on reverse-phase-encoding gradient, known as Topup. 8

| Acquisition
The proposed framework acquires 2 EPI data sets (blip-up and blip-down) with opposite phase-encoding gradient directions. A B 0 scan is acquired as a separate dual-echo gradient echo scan for estimation of the B 0 field.

| Reconstruction
With B 0 field inhomogeneities, the corrupted k-space Y j corresponding to j th coil (j = 1, 2,…,J; J being total number of coils) can be related to the undistorted image x by the following model (Equation 1): where m, n are image coordinate indices, M, N are image dimensions, k, l are k-space coordinate indices, t(k, l) is the sample time for location (k, l) in k-space, ΔB 0 is the B 0 field in Hz that is estimated from a separate dual echo gradient echo scan, C j is the j th coil sensitivity, and j = 1, 2,..,J In the presence of a center frequency offset (∆f 0 ), the model in Equation 1 is modified as Equation 2: In the above expression, the center frequency offset ∆f 0 is assumed to be a constant in space. For a given ∆f 0 , translational shifts in image space are equal and opposite for opposite phase-encoding gradient directions (e.g., positive shift for blip-up and negative for blip-down scans or vice versa).

| Proposed Framework
The proposed reconstruction framework has the following 3 steps. The optimal center frequency offset ∆f 0,opt is estimated as a parameter that maximizes the mutual information (MI) similarity measure 23 between preliminary model-based conjugate phase reconstructions of the blip-up and blipdown EPI data (x cp1 and x cp2 ) acquired with a b-value of 0 s/mm 2 by solving the following unconstrained problem (Equation 6): The B 0 field is corrected by updating it with estimated optimal center frequency offset ∆f 0,opt . The definition of Mutual Information function used in our method is described in Appendix I.

| Step 2. Phase correction
For diffusion-weighted data (b-value > 0 s/mm 2 ), the phase of the blip-up and blip-down data in the image space might be different because of small physiological motion that may occur during the diffusion sensitization gradients. If data are combined from the 2 phase-encoding directions without phase correction, this may result in phase cancellation leading to signal dropout in the final reconstructed image. To cater for this issue, we propose to calculate a phase correction ∆Ф obtained by taking the Hermitian inner product between the preliminary model-based conjugate phase reconstructions (x cp1 and x cp2 ) of blip-up and blip-down EPI data using corrected B 0 field (Equation 7): where <.> indicates the inner product, (.)* denotes the complex conjugate, and m, n are image coordinate indices. In the above expression, the phase correction ΔФ varies in space because of the spatial variance of phases in x cp1 and x cp2 . The raw k-space data are corrected in phase by applying the phase correction ΔФ in the image space as follows.
With k-space Y j 2 selected as a reference, phase correction was applied to Y j 1 by (1) first applying the inverse Fourier transform (IFT) to it, (2) multiplying by an exponential term containing the phase correction ΔФ to get the corrected image ỹ j 1 , and (3) transforming back to k-space by the Fourier transform (FT) to get a final phase-corrected k-space Ỹ j 1 .

| Step 3. Model-based reconstruction
The phase-corrected multicoil k-space data Ỹ 1 and Ỹ 2 and the encoding operators E 1 and E 2 for blip-up and blip-down scans can be expressed mathematically by stacking data and encoding operators from all J coils (Equation 9): Using the center frequency offset corrected ΔB 0 obtained from step 1, for each b-value and diffusion direction, the phase-corrected EPI data Ỹ 1 and Ỹ 2 can be combined into a single formulation by setting Ỹ = [Ỹ 1Ỹ2 ] T and E = [E 1 E 2 ] T in Equation 2. Model-based reconstruction is done ( Figure 1a) from combined k-space data Ỹ in a manner analogous to an iterative SENSE reconstruction 22 based on the conjugate gradient method by solving the following minimization problem (Equation 10): The above minimization problem is strictly convex because the encoding operator E is a linear function in x. Large negative local B 0 field gradients in the phase-encoding direction can make the encoding operator E 1 or E 2 singular or badly ill-conditioned. 12 The convergence of the CG iterations is achieved when the normalized residual r = ||E H Ex−E HỸ || 2 ||E HỸ || 2 (|| ⋅ || 2 denotes the l 2 norm) in the current iteration becomes smaller than ϵ (ϵ being a small number) giving the final distortion corrected output image ̇x. Figure 1b and Figure 1c shows visual illustrations for details of implementation of forward and adjoint encoding operators E i and E H i , respectively; i = 1 and i = 2 refer to blipup and blip-down scans, respectively. The forward encoding operator E i involves multiplication of input image x by individual coil sensitivities C j that is followed by a modified Fourier transform operator G that maps the product of image and coil sensitivity from image to Fourier space, taking into account the susceptibility effects determined by B 0 field and phase-encoding direction dependent k-space sampling times t i . The adjoint encoding operator E H i involves transforming each individual coil k-space data to image space by adjoint of modified Fourier transform operator (G H ) followed by multiplication by the corresponding complex conjugated coil sensitivity (C j ) * and the final summation over all the coils.

| Experiments
Ten male patients (median weight, 84 [range, 68-98] kg and age 68 [57-79] years old) were recruited from the clinical prostate imaging pathway and were consented for additional image acquisitions. No antispasmodic agent was administered. Patients were placed feet first into the scanner and imaging was carried out during free breathing for all patients. The study was approved by the ethics committee, and written signed consent was obtained from all patients for the research scans.

| Data Postprocessing
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, Zurich, Switzerland). EPI phase correction was performed using the ReconFrame tool to correct for ghosts originating from alternating offsets of phase encode lines in k-space. Subsequent postprocessing was implemented in MATLAB (The MathWorks, Inc., Natick, MA).
The B 0 field was calculated using the quantitative susceptibility mapping toolbox 24 that estimates the field by a weighted least squares fit of temporally unwrapped phases in each voxel over echo time. A robust spline-based smoothing 25 was applied to the B 0 field in image domain to smooth out noisy components. The 3D B 0 field and T2W images were resampled to match the EPI scan resolution and the associated FOV. The iterative optimization in Equation 6 to estimate the optimal center frequency offset Δf 0,opt was carried out using the fminsearch function in MATLAB that uses the Nelder-Mead simplex algorithm. 26 The number of iterations in the algorithm was set to 10, and the optimal frequency offset Δf 0,opt was set to the value of Δf 0 found in the last iteration.
This was followed by phase correction (Equations 7 and 8) and model-based reconstruction in Equation 10. The threshold ∊ for convergence of CG iterations in Equation 10 was set to 0.0025 in all our experiments. The model-based corrected reconstruction was compared against uncorrected blip-up and blip-down reconstructions and the neuroimaging distortion correction Topup method.

| Image Analysis
For quantitative evaluation, reconstructed images were converted to DICOM format and exported to Horos software, 27 which is an open-source medical image viewer. A region of interest (ROI) was manually drawn around the boundary of whole prostate in each T2W image and each reconstructed image by a radiologist with 12 years' experience in reporting prostate MRI. Dice similarity coefficient scores 28 were computed from the ROI overlap of each reconstruction and the reference T2W image. For qualitative assessment, the images reconstructed with different methods were put in a random order and blinded to method to avoid subjective assessment. The radiologist scored each image in terms of resolution, distortion, demarcation, and zonal anatomy. 29 The "resolution" was defined as the ability to recognize detailed anatomical structures within the prostate, and it was assessed on a 5point scale (1 = poor; 2 = below average; 3 = average; 4 = above average; and 5 = excellent). The "distortion" was defined as the presence of artifacts, including signal pile-up, signal dropout, warping, ghosting, and blurring. It was assessed on a 5-point scale (1 = severe influence; 2 = significant influence; 3 = moderate influence; 4 = low influence; and 5 = no influence). The "demarcation" was defined as the ability to depict the prostatic capsule in a continuous fashion around the prostate. It was assessed on a 5-point scale (1 = poor; 2 = below average; 3 = average; 4 = above average; and 5 = excellent). "Zonal anatomy" was defined as the ability to distinguish the transitional zone of the prostate from the peripheral zone and was assessed on a 5-point scale (1 = poor; 2 = below average; 3 = average; 4 = above average; and 5 = excellent). F I G U R E 1 Implementation of iterative model-based reconstruction (step 3 of proposed framework) using data from both blip-up and blipdown EPI scans. (a) Given the input k-space data Ỹ = [Ỹ 1Ỹ2 ] T that have been phase corrected in step 2, in each iteration of the conjugate gradient (CG) algorithm, we iterate back and forth between the k-space and image x by encoding operator E = [E 1 E 2 ] T and its adjoint E H that include the center frequency offset corrected B 0 field. The convergence is achieved when the residual r = E H Ex−E HỸ 2 E HỸ 2 in the current iteration becomes smaller than ϵ (ϵ being a small number) giving the final distortion corrected output image ̇x. (b) Details of the forward encoding operator E i that takes the input single image x to the multicoil output k-space Y i for phase-encoding direction i, i = 1 corresponds to the blip-up and i = 2 corresponds to the blip-down scans, respectively. The image x is first multiplied by the coil sensitivity C j (j = 1, 2,..,J). This is followed by a modified Fourier transform operator G that maps the product of image and coil sensitivity from image to Fourier space, taking into account the susceptibility effects determined by B 0 field and k-space sampling times t i , resulting in k-space data Y j i (j = 1, 2,..,J). (c) Details of the adjoint encoding operator E H i that takes the multicoil k-space data Y j i (j = 1, 2,..,J) to single image x for phase-encoding direction i, i = 1 corresponds to the blip-up and i = 2 corresponds to the blip-down scans, respectively. Each individual coil k-space data Y j i is transformed to image space by the adjoint of a modified Fourier transform operator (G H ). The resulting images are individually multiplied by complex conjugated coil sensitivities (C j i ) * (j = 1, 2,..,J) and summed to give the final image x | 1985 M. USMAN et al.

| Phase correction procedures
The phase of DWI data (b-value, >0 s/mm 2 ) can change dynamically because of small physiological motion occurring during the diffusion sensitization gradients. To evaluate the performance of the phase correction procedure used in our proposed framework, a mid-prostate single slice was acquired for b-value of 500 s/mm 2 on 2 patients for 40 dynamics with blip-up phase-encoding gradient. Model-based reconstruction was done with and without phase correction for each pair of dynamics (dynamic 1 and dynamic 2, dynamic 1 and dynamic 3,….) in each diffusion direction with dynamic 1 selected as a reference. The goal is to check whether the proposed phase correction can work in cases where the dynamics have significantly different phases from one another in the conjugate phase reconstructions.

| Estimation of center frequency offset (∆f 0 )
For all patients, the values of optimal center frequency offset (in Hz) for a mid-prostate slice are shown in Figure 2a, together with a convergence plot of objective function (negative of mutual information function) as a function of iteration number in Figure 2b. The curve in Figure 2b shows that the objective function minimum (or equivalently mutual information function maximum) was achieved in 5 to 6 iterations. Mean frequency offset across all patients was 47.15 Hz with a standard deviation of 5.2 Hz.

| Model-based reconstruction
In Figure 3, for patient 1, the reference T2W image, the B 0 field, uncorrected blip-up and blip-down reconstructions, correction using blip-up data only and blip-down data only, and proposed model-based reconstructions are shown for bvalues of 0 and 500 s/mm 2 , respectively. Because of susceptibility differences at the interface between different tissues in the prostate region, the variation in B 0 was around 120 Hz that corresponds to a shift of 7 to 8 pixels at a bandwidth/ pixel of ~15.80 Hz in phase-encoding direction. This resulted in pileup in regions where the B 0 field gradient magnitude was high and its direction was opposite to that of phaseencoding gradient (see regions pointed out by red arrows in Figure 3). The reconstruction results for patient 2 are shown in Figure 4. The proposed method corrected for signal pileup artefacts that cannot be corrected using data from only 1 phase-encoding gradient direction. For a mid-prostate axial slice in a selective group of 3 patients (patient 4, patient 6, and patient 9), the proposed and Topup reconstructions for b-values of 0 and 500 s/mm 2 are shown in Figures 5 and 6, respectively. The complete set of results for patient 3 to patient 10 for b-value of 0 and 500 s/mm 2 can be seen in Supporting Information Figures S1 and S2, respectively. Supporting Information Figure S3 shows an example plot of normalized residual error r as a function of CG iteration number.  Tables 1 and 2, respectively. Mean percentage of improvement in qualitative scores using the proposed method compared to other reconstructions is shown in Supporting Information Figure S4. The proposed method performed better, on average, than all other reconstructions for each qualitative measure.

| DISCUSSION
A novel model-based reconstruction framework is proposed that can correct for geometric distortions, signal pileup, and signal dropout in diffusion-weighted prostate images. By using the power of complimentary encoding information within data from both blip-up and blip-down directions, the model-based framework was able to correct most of the pileup in regions of severe distortions and F I G U R E 4 In vivo patient 2 reconstruction results for b-value of 0 and b-value of 500 s/mm 2 are shown in (a) and (b), respectively. The reference T2W image and estimated B 0 field (in Hz) are shown in (a). For EPI scans, the AP axis was selected as the phase-encoding direction with fat shift in the direction "P" for blip-up and "A" for blip-down scans, respectively. Whole prostate (red) was delineated on reference T2W and overlaid on reconstructions without distortion correction (uncorrected blip-up [UC P] and uncorrected blip-down [UC A]), reconstruction with distortion correction using data from 1 direction only (corrected blip-up [C P] and corrected blip-down [C A]), model-based reconstruction using both blip-up and blip-down data (C AP), and Topup method (Topup). The proposed method performed better than all other reconstructions in terms of distortion correction performed better than the Topup method (commonly used for neuroimaging) or model-based reconstructions using data from 1 phase-encoding gradient direction only. The proposed mutual information maximization used for correction of center frequency offset in the B 0 field may also be used as a method for overcoming uncertainties in coordinates that may happen in the process of reconstruction from raw data.
The proposed method assumes the B 0 field to be static except for center frequency offset attributed to the frequency drift between B 0 and EPI scans. In case of motion or changes in the rectal area adjacent to the prostate region between the EPI and B 0 scans, the B 0 field estimated from B 0 scan may be mismatched, resulting in inaccurate distortion correction. Some of the remaining distortions in the proposed method reconstructions can be attributed to dynamic changes in B 0 F I G U R E 5 In vivo reconstruction results for selected patients (P4, P6, and P9) for data acquired at b-value of 0 s/mm 2 . For EPI scans, the AP axis was selected as the phase-encoding direction with fat shift in the direction "P" for blip-up and "A" for blip-down scans, respectively. Whole prostate (red) was delineated on a reference T2W image (left column) and overlaid on uncorrected blip-up (UC P), uncorrected blip-down (UC A), Topup, and proposed distortion-corrected (C AP) reconstructions F I G U R E 6 In vivo reconstruction results for selected patients (P4, P6, and P9) for data acquired at b-value of 500 s/mm 2 . For EPI scans, the AP axis was selected as the phase-encoding direction with fat shift in the direction "P" for blip-up and "A" for blip-down scans, respectively. Whole prostate (red) was delineated on a reference T2W image (left column) and overlaid on uncorrected blip-up (UC P), uncorrected blip-down (UC A), Topup, and proposed distortion-corrected (C AP) reconstructions field because no antispasmodic agent drug was administered in our scans that would suppress bowel movements and/or rectal gas. Changes in B 0 might be addressed by a joint estimation of both B 0 and the corrected EPI images, starting with the initial B 0 fields estimated from B 0 scan. 30 For higher b-values (b-value > = 1000 s/mm 2 ), the low SNR in the reconstructed images might affect the inverse problem and the technique may benefit from preconditioning or additional regularization.
Our proposed method uses a B 0 field estimated from a dual-echo gradient echo scan. In image regions where the individual echoes have low SNR or missing signal (especially in the rectal-air region), the measurement of B 0 field in that area might not be always possible, leading to noisy and unwrapped phases in the B 0 field. To address this issue, B 0 field might be modeled in this region using a susceptibility map distribution 31 that could be obtained from a segmentation of air and tissue areas on a reference T2W image. Alternatively, a projection onto dipole fields (PDF) method 32 might also be used for estimating the B 0 field inside the rectal-air region. The PDF method calculates the B 0 field inside an ROI by projection of known dipole fields from outside the ROI.
Last, the proposed method does not include any physiological motion effects that may occur between the EPI blipup and blip-down scans. In case of motion between blip-up and blip-down scans, both motion parameters and the image would need to be estimated simultaneously and the optimization in Equation 10 becomes nonconvex. The optimization may be simplified by estimating motion in a preceding step similar to the framework in a previous work 33 by registering model-based reconstructions from blip-up and blip-down scans with either blip-up or blip-down reconstruction being set as the reference. The estimated motion fields can then be incorporated into Equation 10 by motion matrix transformation 17,33 that will relate the acquired corrupted k-space  data to corrected images through a combined model that is a concatenation of B 0 distortion and motion corruption effects.
Prostate diffusion MRI is recognized as a potential biomarker for tumor detection, but, currently, it is unusable in some patients because of significant distortions. We proposed a novel model-based reconstruction framework that can correct these distortions by using data from opposite phase-encoding gradient directions. The proposed method was applied successfully in 10 clinical patients, despite no antispasmodic drug being administered. The proposed technique may offer potential to radiologists and clinicians by increasing the diagnostic value of prostate images for tumor detection, thus making prostate MRI a more reliable and reproducible biomarker in future.

ACKNOWLEDGMENTS
We acknowledge research support provided by Philips Healthcare.

David Atkinson
http://orcid.org/0000-0003-1124-6666 In vivo reconstruction results for patients 3 (P3) to patient 10 (P10) for data acquired at b-value of 500 s/mm 2 . For EPI scans, the AP axis was selected as the phase-encoding direction with fat shift in the direction "P" for blip-up and "A" for blip-down scans, respectively. Whole prostate (red) was delineated on reference T2W image (left column) and overlaid on uncorrected blip-up (UC P), uncorrected blipdown (UC A), Topup, and proposed distortion-corrected (C AP) reconstructions