FULL PAPER Magnetic Resonance in Medicine 75:169–180 (2016)

POCS-Enhanced Inherent Correction of Motion-Induced Phase Errors (POCS-ICE) for High-Resolution Multishot Diffusion MRI Hua Guo,1* Xiaodong Ma,1 Zhe Zhang,1 Bida Zhang,2 Chun Yuan,1,3 and Feng Huang4 Purpose: For multishot diffusion weighted imaging (DWI), one of the challenges is to remove phase variations induced by physiological motion among different shots. In this study, a new method is proposed to iteratively solve the phase errors and DWI images simultaneously, for navigator-free acquisitions. Theory and Methods: Instead of solving phase errors and the image sequentially in the two-step parallel imaging, the proposed method, named POCS-enhanced Inherent Correction of motion-induced phase Errors (POCS-ICE), treats both the phase and DWI image as unknowns and solves them simultaneously. Multishot DWI with constant density spiral trajectory served as a specific example. Simulation and in vivo experiments were performed to evaluate the proposed method. Results: POCS-ICE shows improved image quality in terms of higher SNR and fewer artifacts than the compared method, SENSEþCG. The improvement becomes more conspicuous as the number of shots increases. The convergence behavior of POCS-ICE was also shown to be more stable. Conclusion: POCS-ICE can inherently and reliably correct motion-induced phase errors in navigator-free multishot DWI, and it is easier to determine the stopping criterion without manual interventions. The improved spatial resolution and image resolvability are beneficial to study of brain microstructures and physiological features for neuroscience. Magn C 2015 Wiley Periodicals, Reson Med 75:169–180, 2016. V Inc. Key words: diffusion weighted imaging; multishot diffusion imaging; navigator-free; SENSE1CG; POCS-ICE

INTRODUCTION Diffusion weighted MR imaging (DWI) in vivo can not only probe microstructures (1,2) but also yield functional information of tissues (3,4), thus it has been widely 1 Center for Biomedical Imaging Research, Department of Biomedical Engineering, Tsinghua University, Beijing, China. 2 Healthcare Department, Philips Research China, Shanghai, China. 3 Department of Radiology, University of Washington, Seattle, Washington, USA. 4 Philips Healthcare (Suzhou) Co., Ltd., Suzhou, China Grant sponsor: National Natural Science Foundation of China; Grant number: 61271132; Grant sponsor: Beijing Natural Science Foundation; Grant number: 7142091. *Correspondence to: Hua Guo, Ph.D., Center for Biomedical Imaging Research, Department of Biomedical Engineering, Tsinghua University, Beijing, China. E-mail: [email protected]

Received 20 September 2014; revised 4 November 2014; accepted 7 December 2014 DOI 10.1002/mrm.25594 Published online 3 February 2015 in Wiley Online Library (wileyonlinelibrary. com). C 2015 Wiley Periodicals, Inc. V

adopted in clinical diagnosis and neuroscience research. Single-shot EPI is conventionally used for DWI acquisition due to its high speed and motion insensitivity. However, because single-shot echo-planar imaging (EPI) has a limited readout window, the resultant images usually have limited spatial resolution and are distorted by magnetic field inhomogeneity and eddy currents induced by diffusion gradients (5,6). Partially parallel imaging (PPI) methods (7–9) have been routinely used to either improve the resolution or reduce geometric distortion in single-shot EPI DWI through undersampling the k-space (10,11). However, PPI methods in clinical practice are limited to a small acceleration number. Consequently, it is difficult to use single-shot EPI to achieve the high spatial resolution DWI images desired for small lesion detections (12) and intracortical fiber tracking (13). Because multishot techniques can shorten the duration of the readout window and hence provide higher spatial resolution than single-shot techniques, they have become the common choice for high resolution DWI. However, due to the application of diffusion gradients, DWI is very sensitive to any kind of motion. This includes involuntary motion such as cardiac pulsation (14), which can induce linear (15) and nonlinear phase changes (16) from shot to shot. This, in turn, causes k-space lines to displace from their original positions and thus results in severe image artifacts, including aliasing and signal cancellation. Therefore, the major difficulty in multishot DWI is how to correct motion-induced phase errors. The methods for phase error correction in multishot techniques can be categorized into three classes: extra-navigated methods, self-navigated methods, and navigator-free methods. Extra-navigated methods acquire additional navigator signals as well as the imaging signals. It is assumed that the navigator signals experience the same motion as the to-be corrected imaging signals. Hence the phase information can be derived from the navigator and used to correct the phase errors in imaging signals. This method was first implemented on interleaved EPI (17–19) and then extended to several schemes including read-out segmented EPI (20,21), interleaved EPI with SENSE (22), and other multishot EPI methods (23,24). They have been widely explored and studied. The major issue of extranavigated methods is the reduced acquisition efficiency. Self-navigated methods use the imaging data itself to generate the phase information. Some imaging trajectories inherently contain a fully sampled k-space region for each shot, such as PROPELLER (25) and multishot variable density spiral (26), providing low resolution images for phase correction. These approaches have successfully been used for high resolution DWI. In the

169

170

Guo et al.

meantime, the acquisition efficiency is also compromised because the readout window is expanded. Several multishot DWI methods without use of a navigator have been brought forward. Skare et al (27) proposed to use PPI to individually reconstruct an image of each shot and then average their magnitudes after removing the phase variation. This method is easier to implement than methods to be reviewed below. However, it suffers from a low signal-to-noise ratio (SNR) induced by PPI and is, therefore, limited to a small shot number. Truong et al (28) proposed an iterative method to directly estimate the point spread function for multishot spiral diffusion imaging, but this method is limited to a linear phase model and small shot numbers, and is time consuming in computation. Recently, a two-step PPI strategy was proposed for navigator-free multishot DWI. In the first step, PPI was used to initially reconstruct an image of each shot. In the second step, an integrated equation incorporating all shots with phase errors obtained from the initial images was again solved by PPI to reconstruct the final image. This strategy has been implemented on STEAM DWI, also known as Inverse Reconstruction (29), as well as EPI DWI (a.k.a. MUSE) (30) and spiral DWI [a.k.a. SENSEþCG (conjugate gradient)] (31). The key advantage of these methods is that they use the constraint that the images of all shots share the same magnitude. Hence, the integrated equation has more sub-equations for the same amount of unknowns, which improves matrix inversion conditioning and results in a lower g-factor. The drawback is that the initial reconstruction by PPI may not provide accurate phase information if the undersampling rate (shot number) is high, so these methods usually require that the shot number is smaller than the number of radio frequency (RF) coil elements. To solve the problems in the existing methods, an iterative reconstruction method for multishot navigator-free DWI is proposed in this study, entitled Projections Onto Convex Sets (POCS) enhanced Inherent Correction of motion-induced phase Errors (POCS-ICE). As described in our recent work (32), it borrows the assumption from the two-step PPI strategy that the image of each shot shares the same magnitude, but improves the performance of phase correction through iterative reconstruction based on POCS. A multishot spiral DWI method is investigated as a demonstration example in this study. Numerical simulation and in vivo experiments show that POCS-ICE can provide high resolution DWI images with reduced artifacts, even with larger shot number than traditionally used values, and show more stable convergence behavior than the compared method, SENSEþCG. THEORY Reformulation of Multishot DWI Reconstruction In the multishot DWI reconstruction, the energy function we want to minimize is X   E ðI Þ ¼ jjGi FT I  fi  Sj  di;j jj22 [1] i;j

where I is the complex image to be reconstructed, fi is the exponential of the phase error for shot i, Sj is the coil sensi-

tivity map of channel j, FT is the Fourier transform which transforms the image into k-space, Gi is a resampling operation transferring the k-space data from a regular Cartesian grid to the original sampling locations of shot i, and di;j is the acquired k-space data of shot i and channel j. The previous methods, for example, SENSEþCG (31) in spiral DWI, use SENSE to individually reconstruct the image of each shot and then calculate fi to extract its phase information in the first step. In the second step, the calculated fi is used to minimize E and thus to reconstruct the final image I, the only unknown. Although this idea is straightforward and can be solved using the traditional PPI algorithm, the separately calculated fi may be inaccurate due to a high reduction factor of each shot and the final reconstructed image may be degraded. What we propose is that both fi and I are treated as unknowns and solved simultaneously by minimizing E. Hence, the precision of fi can be iteratively improved along with I. A specific example to implement this idea is given below. First, we define the energy function of each individual shot: X   Ei ðI; fi Þ ¼ jjGi FT I  fi  Sj  di;j jj22 : [2] j

Then E ðI; fi Þ ¼ ¼

X X i X

  jjGi FT I  fi  Sj  di;j jj22



j

ðEi ðI; fi ÞÞ:

[3]

i

Based on this formulation, we can split the minimization into two parts, minimizing Ei of each individual shot based on Eq. [2] as an intermediate step and then minimizing the total energy function E based on Eq. [3]. The two-part minimization can be solved together using a POCS-based MRI reconstruction for Sensitivity Encoding (POCSENSE) (33) algorithm. The first part is the data projection in k-space; and the second serves as an image domain constraint. The details are explained further in the following subsections. Projections Onto Convex Sets (POCS) The POCS reconstruction method, based on the Gerchberg-Saxton iterative algorithm (34), is to solve inverse problems by iteratively projecting an initial estimate onto the predefined convex sets. This method has been modified by introducing a priori constraints for improved image restoration from incomplete or inconsistent data (35). POCSENSE shows superiority in highly undersampled situations when compared with traditional SENSE reconstructions, and can also be used to identify and discard k-space data corrupted by motion to generate artifact-free images (36). In this study, the projection operator in POCSENSE is modified for use in the reconstruction of multishot spiral diffusion imaging. Modified Data Projection Operator In the original POCSENSE, the data projection operator P projects the initial image onto V in each iteration, where

High Resolution Multishot DWI without Navigator

V is defined as a set including all the images whose kspace values are equal to originally acquired data d at sampled k-space locations (33). In nondiffusion weighted imaging, the images of all shots are assumed to be identical, so that data of all shots are reconstructed together by one single operator. In multishot DWI, however, because images of different shots have different phases, they should be handled separately during the data projection. Specifically, for multishot spiral DWI with multiple channels, the data projection operator Pi;j for shot i and channel j is defined as Pi;j Ii ðrÞ ¼ Ii ðrÞ  Sj ðrÞ n  o [4] þ IFTGI di;j ðkÞ  GFðk2Ki Þ FT Ii ðrÞ  Sj ðrÞ where r and k denote image space and k-space coordinates, respectively; Ii is the initial image of shot i, which is equal to the desired image IðrÞ multiplied by fi ðrÞ, i.e., Ii ðrÞ ¼ IðrÞ  fi ðrÞ; Sj is the sensitivity map of channel j; di;j is the acquired data at k-space locations Ki , which is the k-space trajectory of shot i; FT and IFT are the forward and inverse Fourier transforms; GF and GI are forward and inverse resampling operators which transfer data from Cartesian grids to an arbitrary trajectory and vice versa, which can be implemented together with Fourier transforms using either the Gridding algorithm (37) or the nonuniform Fourier transform (NUFFT) (38). Similar to the original POCSENSE, for multishot DWI the convex set Vi;j for each shot that Pi;j projects onto contains all the images whose k-space values are equal to the originally acquired data di;j at sampled k-space locations. In other words, the data projection operator individually updates the image of each shot by modulating it with the coil sensitivity, calculating the k-space error, transferring the error into image space and adding it to the initial estimation. The Same Magnitude as a Constraint After data projection, the same magnitude shared by different shots is taken as a constraint in the following POCS-ICE reconstruction. Specifically, images of all shots are averaged after the data projection in each iteration. Meanwhile, the phase of each shot includes its own phase error, which should be removed before averaging. In other words, a shot averaging strategy combined with phase correction is adopted in this study, and is regarded as a projector for this constraint. This is the second part of POCS-ICE which minimizes the cost function E defined in Eq. [3]. Further details can be found in the “Algorithm” section. Algorithm The proposed reconstruction algorithm, POCS-ICE, consists of five steps in each iteration: data projection, channel combination, shot averaging, image update, and ðnÞ phase recovery. Given the image estimation Ii ði ¼ 1; ::: NÞ from the nth iteration, the current iteration of number nþ1 is carried out as follows. Step 1, data projection:

171 ðnÞ

gi;j ¼ Pi;j Ii ; i ¼ 1; :::; N; j ¼ 1; :::; Nc:

[5]

Previous image estimation is projected to get the image of each individual channel gi;j . N is the total shot number, and Nc is the total channel number. Step 2, channel combination: hi ¼

Nc X j¼1

Sj  gi;j ; i ¼ 1; :::; N XNc S  Sj j¼1 j

[6]

gi;j of different channels are combined through an SNR optimal reconstruction, assuming the noise correlation is an identity matrix (39). Sj is the sensitivity map of coil j, and the notation * denotes complex conjugate. Step 3, shot averaging:    ^ ¼ IDFT Wtri  DFT h ; i ¼ 1; :::; N h i i N ^ 1X h Iavg ¼ h  i ; i ¼ 1; :::; N: ^ j N i¼1 i jh i

[7] [8]

Here, DFT is the discrete Fourier transform and IDFT the inverse discrete Fourier transform. The notation jj denotes the magnitude value pixel-by-pixel. In this step, images of the different shots hi obtained in step 2 are averaged after removing their low resolution phases, based on the fact that phase errors are relatively smooth unless severe bulk motion occurs and very high b values are used (40). The high frequency phases need to be kept to maintain the essential information for off-resonance correction; and they can be used for complex average of different shots to improve SNR (41). First, hi are Fourier transformed into k-space, windowed by a 2D triangular window Wtri and then inverse Fourier transformed ^ , as shown in Eq. (25,27), to get low-resolution images h i [7]. Second, hi are averaged after removing the low reso^ ), to get I , as shown lution phases (i.e., the phases of h i avg in Eq. [8]. Step 4, image update:  Iuðnþ1Þ ¼ IuðnÞ þ l Iavg  IuðnÞ :

[9]

The image estimation is updated by adding the differðnÞ ence between Iavg and the previous estimation Iu , weighted by a relaxation factor l which can affect the convergence speed (33). Step 5, phase recovery: ðnþ1Þ

Ii

¼ Iuðnþ1Þ 

^ h i ; i ¼ 1; :::; N: ^ j jh

[10]

i

For each shot, the low resolution phase is restored to ðnþ1Þ ðnþ1Þ Iu to get a new image of each shot Ii with complete phase information. This is then used as the initialization for the next iteration. The schematic diagram of the POCS-ICE algorithm is shown in Figure 1. The five steps mentioned above will be repeated until the maximum iteration number is reached or the iteration converges, which means the residual error dðnþ1Þ is smaller than a predefined tolerance t.

172

Guo et al.

FIG. 1. Schematic diagram of POCS-ICE. The main diagram on the left shows the algorithm, where the five steps in each iteration (enclosed by the dashed rectangles) are included. The right diagram shows the implementation details of the projection operator.

ðnþ1Þ

ðnÞ

ðnÞ

dðnþ1Þ is defined as: dðnþ1Þ ¼ jjIu  Iu jj2F =jjIu jj2F , where jjjjF denotes the Frobenius norm, i.e., the root of sum-ofsquare of all the elements in one matrix. METHODS Numerical Simulation To validate the POCS-ICE algorithm quantitatively, a numerical simulation was performed using a noise-free T1 weighted digital phantom with matrix size of 256  256, downloaded from BrainWeb (http://www.bic.mni. mcgill.ca/brainweb/). The eight-channel complex coil sensitivity maps were computed using the Biot-Savart law analytically. Based on the phantom and sensitivity maps, multishot spiral diffusion weighted data of multiple channels were generated by means of the following procedures: (i) multiplication of the phantom with a spatially second-order phase map, which is random for different shots, and with the range comparable to the phase errors estimated from in vivo data; (ii) modulation using sensitivity maps, (iii) transferral to k-space, and (iv) addition of spatially uncorrelated Gaussian noise. To test the robustness of the algorithm, shot numbers of 5 to 8 and SNR levels ranging from 5 to 20 dB were tried separately. The normalized root mean square error (nRMSE) was used to evaluate the quality of reconstructed images, which ^ F =jjmref jjF , where mref is calculated as nRMSE ¼ jjmref  mjj ^ is the reconis the ground-truth phantom image and m structed image. Each combination of shot number and SNR level was repeated 10 times, and the nRMSE value each time was recorded.

In Vivo Experiments Eight healthy volunteers were recruited and provided informed written consent for in vivo studies. Brain diffusion tensor imaging (DTI) data were acquired on a Philips 3 Tesla scanner (Philips Healthcare, Best, The Netherlands) using an eight-channel RF coil. A spinecho diffusion weighted spiral sequence without navigator was used. Three experiments were carried out using the following imaging parameters. In the first experiment, relatively benign imaging conditions were tested, with the number of spiral shots (NSS) ¼ 6, spiral readout duration for each shot ¼ 30 ms, and b ¼ 800 s/mm2. In the second experiment, NSS ¼ 8, and spiral readout duration ¼ 24 ms. To assess different SNR levels, two different b-values, 800 s/mm2 and 1600 s/mm2 were applied separately. The first and second experiments acquired data at the same imaging locations and shared the following imaging parameters: FOV ¼ 212  212 mm2, slice thickness ¼ 4 mm, TR/TE ¼ 2680/54 ms, in-plane image resolution ¼ 0.95  0.95 mm2, the number of diffusion directions ¼ 15 for each b value with additional nondiffusion weighted data (b ¼ 0), and Number of Signals Averaged (NSA) ¼ 2. It is worth noting that two NSAs were only acquired for quantitative analysis, which includes the noise map calculation that is required for SNR analysis and scatter plots of FA quantification. However, the reconstruction results presented below are of one NSA. Besides spiral DWI, single-shot EPI DTI data were also acquired as a reference using a clinical sequence, and the corresponding imaging parameters were:

High Resolution Multishot DWI without Navigator

in-plane image resolution ¼ 2  2 mm2, TR/TE ¼ 2680/102 ms, NSA ¼ 4, SENSE reduction factor ¼ 1.8. Other imaging parameters were the same as those of the spiral DWI. In the third experiment, NSS was increased to 10, which is larger than the number of RF coil elements. The imaging parameters were TR/TE ¼ 2568/54 ms, FOV ¼ 220  220 mm2, b ¼ 1000 s/mm2, the number of diffusion directions ¼ 15, in-plane image resolution ¼ 0.9  0.9 mm2, slice thickness ¼ 4 mm, NSA ¼ 1, spiral readout duration ¼ 21 ms. In this experiment, a larger shot number was used, resulting in a shorter readout duration and, therefore, a reduced number of off-resonance blurring artifacts. Like the first two experiments, single-shot EPI DTI data were acquired as a reference: in-plane image resolution ¼ 2  2 mm2, TR/TE ¼ 2500/69 ms, NSA ¼ 3, SENSE reduction factor ¼ 2. Other imaging parameters were consistent with those of the spiral DWI.

Image Reconstruction For comparison, the proposed POCS-ICE and the twostep PPI method, SENSEþCG, were used to reconstruct images from either simulated or in vivo data. The NUFFT was used to transform data from image space to k-space. For in vivo data, the coil sensitivity maps were calculated from b ¼ 0 images. In POCS-ICE, the initial images f ð0Þ were set to zero. To assure convergence, the maximum iteration number was set to 200, and the tolerance t was empirically set to 108 for simulation and 106 for in vivo data. The relaxation factor l was set to 1 because a larger relaxation parameter would not improve the quality of reconstructed image using such a large iteration number (42). The width of triangular window Wtri was set to half of the matrix size, that is, 118 for the first two experiments and 122 for the third experiment. In SENSEþCG, the algorithm and parameters used here were slightly different from those used in the original study (31), but were tested to have better performance in the reconstruction. The residual norm was measured as the stopping criterion (43) and the tolerance t was set to the same value as that in POCS-ICE. For the first step, the whole k-space data of each shot rather than the low resolution data were used to calculate the phase error. The maximum iteration number in the first and second step was empirically set to 12 and 10, respectively, for the simulation, and 12 and 8, respectively, for the in vivo data. Total variation denoising (44) (instead of median filtering from the original study) was used to smooth the initial reconstruction for better phase error estimation, and the denoising factor, which was used to balance the total variation term and data fidelity, was set to 0.4. All the reconstruction procedures were implemented using Matlab (Mathworks Inc., Natick, MA) on a PC computer with 3.10 GHz CPU and 32 GB RAM. The convergence iteration number of POCS-ICE was around 50 for 6-shot, 100 for 8-shot, and 200 for 10-shot. The computation time to reconstruct a single image was approximately 0.3, 0.5, and 1 min for three different shot

173

numbers, respectively, using SENSEþCG, and was approximately 5 times longer using POCS-ICE. To suppress blurring artifacts, automatic off-resonance correction without field maps (45,46) was applied to the images reconstructed by SENSEþCG or POCS-ICE. In the end, fractional anisotropy (FA) maps were calculated using DTI studio (47). RESULTS Numerical Simulation A representative set of simulation results using eightshot spiral acquisitions and SNR level of 10 dB is shown in Figure 2. Because of artificially introduced phase which is different for each shot, direct reconstruction without any correction results in severe aliasing artifacts and signal cancellation (Fig. 2d). Using SENSEþCG, those artifacts are reduced but residual aliasing still exists (Figs. 2b and e). On the contrary, POCS-ICE generates images with fewer artifacts and improved SNR (Figs. 2c and f). The improvements in POCS-ICE results are considered to be from more accurate phase estimation than SENSEþCG, in which SENSE with a large acceleration factor (R ¼ 8) was used for reconstruction of each shot in the first step. The evidence can be found in Figure 3, wherein the phase errors estimated from two methods are compared and it can be seen that POCS-ICE provides phase errors closer to the ground truth than SENSEþCG. Quantitative analysis of numerical simulation is shown in Figure 4. The nRMSE plots show that for different shots and different SNR levels, POCS-ICE provides more precise reconstruction results with less variation (Figs. 4a and b). When the shot number is increased to eight, the image quality is maintained well compared with five-shot (8% 6 0.83% versus 5% 6 0.64%) and is not disturbed much by the SNR (Fig. 4b). The convergence plots (Fig. 4c) demonstrate that POCS-ICE has more stable convergence behavior than SENSEþCG. Although the nRMSE of SENSEþCG decreases fast at the beginning, it goes up when the iteration increases further. This was caused by the intrinsic noise-amplification feature of CG iteration in the ill-conditioning problem, which was referred to as semiconvergence (43). The nRMSE of POCS-ICE, however, decreases monotonically during the whole reconstruction process. The persistent convergence suggests that it is easier to determine the stopping criterion for POCS-ICE than SENSEþCG, without manual interventions. In Vivo Experiments Figure 5 shows the in vivo results of six-shot acquisition using the eight-channel RF coil from the first experiment. Even though the DWI images and FA maps of the two methods are comparable with each other in terms of artifact levels, POCS-ICE shows higher SNR and anatomical fidelity, which can be visually observed in the diffusion weighed images (top row in Figure 5). The improved structural detail in POCS-ICE can be validated by the zoomed-in inset from the b ¼ 0 image at the

174

Guo et al.

FIG. 2. Reconstructed images from simulated eight-shot data with SNR ¼ 10 dB. a: Reference image without phase errors. b: SENSEþCG. c: POCS-ICE. d: direct reconstruction without correction. e,f: Error maps of b and c, respectively, amplified by a factor of 2.

bottom right corner of the POCS-ICE image. In the FA maps, POCS-ICE shows improved directionality delineation and lower noise level, highlighted by the white arrows in the bottom row of Figure 5.

DWI images and zoomed-in FA maps of the eight-shot spiral acquisitions in the second experiment are shown in Figure 6. Compared with SENSEþCG, POCS-ICE shows reduced artifacts and thus better delineation of

FIG. 3. Comparison of phase error estimation between SENSEþCG and POCS-ICE. Three representative shots are shown from top to bottom. Note: error map was the subtraction of the ground truth from the estimated phase errors after phase wrapping from the ground truth.

High Resolution Multishot DWI without Navigator

175

FIG. 4. Quantitative comparison between SENSEþCG and POCS-ICE. a: nRMSE of different shot number when SNR ¼ 10 dB. b: nRMSE of different SNR level when shot ¼ 8. c: Convergence plot of POCS-ICE and the second step of SENSEþCG for 6 shots and 8 shots, which shows nRMSE of SENSEþCG drops quickly in the first few iterations but then increases when the iteration number increases to around 8, while nRMSE of POCS-ICE decreases monotonically. For display, only the beginning 50 iterations are shown in c).

the structure for both b values tested, which is consistent with the conclusion from the simulation results. Due to the high acceleration factor, the results from SENSEþCG show increased noise, especially in regions around the brain ventricle, and strong aliasing artifacts (yellow arrowheads in the second column). The improved image quality is further demonstrated by the FA maps in the third and fourth rows, where POCS-ICE shows better directionality and higher SNR than SENSEþCG (white

FIG. 5. DWI images (upper row) of one representative direction and corresponding FA maps (bottom row) from the six-shot acquisition in the first experiment. The columns from left to right are results of SENSEþCG and POCS-ICE, respectively. While both methods provide images with the a similar artifact level, POCSICE shows improved SNR and more reliable anatomic information, e.g., the vessel indicated by the yellow arrows in the top row and the posterior callosum pointed out by white arrows in the bottom row. The zoomed-in inset at the bottom-right corner of the POCSICE image is a reference image from b ¼ 0.

arrowheads). Compared with the traditional single shot EPI DWI with low resolution, higher resolution imaging provides higher resolvability, especially in the cerebral cortex. In addition to the final reconstructed images shown above, the estimated phase errors by two methods are shown in Figure 7. As we can see, increased b values led to an increase in motion sensitivity (wrapped phase in the rightmost column) (40). Even though the true phase errors can never be known, the estimation by POCS-ICE is considered more precise because it can produce DWI images with reduced artifacts. For quantitative assessment of the proposed algorithm on in vivo data, the SNR values of the DWI images reconstructed by the proposed algorithm and SENSEþCG were compared. The noise map was obtained by subtracting the image of one NSA pffiffiffifrom another, and the SNR was  (48), where S  is mean intendetermined by SNR ¼ 2S=d sity of one NSA, and d is the standard deviation of the noise map. For the first experiment (six-shot and b ¼ 800 s/mm2), the SNR from POCS-ICE was on average 11.36% 6 8.41% (from 15 directions) higher than that of SENSEþCG. The effect became more prominent when the number of shots and b values increased. For the eightshot acquisition, the average increase percentage of SNR amounts to 17.04% 6 9.55% and 22.31% 6 8.03% for b ¼ 800 and 1600 s/mm2, respectively. It should be noted that because the images of different NSA may contain different artifacts, the noise map obtained through subtraction may be contaminated. Considering the random feature of artifacts across diffusion directions, the SNR shown here is still meaningful in some ways. Using the eight-shot acquisition with b ¼ 800 s/mm2 as an example, Figure 8 shows the relation between the single shot EPI and both POCS-ICE and SENSEþCG for FA quantification. These values were calculated from the regions of interest (ROIs) chosen on the posterior callosum (rectangle in Figure 6). The consistency for FA assessment between SENSEþCG and the reference is low regardless of SNR levels used. In contrast, POCS-ICE demonstrated better consistency, especially when SNR was increased. Figure 9 shows the results from 10-shot acquisition from the third experiment. Because the NSS is larger than the coil elements, the DWI image reconstructed by SENSEþCG has residual aliasing artifacts (yellow

176

Guo et al.

FIG. 6. DWI images of one representative direction and zoomed-in FA maps from eightshot spiral acquisitions, with two different b values (b ¼ 800 and 1600 s/mm2). In the first column, low resolution images with the same diffusion encoding using single-shot EPI are shown as the reference images. In the other columns, high resolution images reconstructed by SENSEþCG and POCS-ICE are shown from left to right. The arrows highlight the regions where POCS-ICE has fewer artifacts than SENSEþCG. The rectangular ROI in the single-shot EPI FA map (the third row) is for quantitative analysis.

arrows). The images by POCS-ICE have reduced artifacts, and the FA map shows detailed structure better. For instance, the distribution of white matter in the corona radiata (hollow arrows) is more consistent with low resolution image by single-shot EPI DTI, although its reliability is still under investigation.

DISCUSSION AND CONCLUSIONS In this study, a new reconstruction method named POCS-ICE is proposed for multishot navigator-free DWI, which can iteratively correct motion-induced phase errors based on the POCS algorithm. Numerical simulation and in vivo experiments present superior image

FIG. 7. Phase errors estimated by SENSEþCG and POCS-ICE from in vivo data with different shot numbers and b values. Only one representative shot is shown.

FIG. 8. Scatter plots show the FA measurement by SENSEþCG and POCS-ICE compared with the reference acquisitions. The FA values are from the ROIs drawn on the posterior callosum in Figure (6) on the leftmost column. Two different NSA values were evaluated for the eight-shot acquisition with b ¼ 800 s/mm2.

178

Guo et al.

FIG. 9. DWI images (upper row) of one representative direction and the corresponding FA maps (bottom row) from the 10-shot acquisition using an 8-channel coil. From left to right, results from single shot EPI, SENSEþCG and POCS-ICE are shown. The solid arrows show the region where artifacts are reduced in POCS-ICE compared with SENSEþCG. The hollow arrows in FA maps show the corona radiata, where the white matter distribution by POCS-ICE is more consistent with single-shot EPI.

quality using POCS-ICE compared with SENSEþCG, especially when the shot number is relatively high. Although multishot spiral DWI was used as an example, the POCS-ICE algorithm can be easily extended to any arbitrary trajectories, e.g., interleaved EPI. Furthermore, the phase estimation method proposed here can be extended in other applications, such as multishot fMRI (49) The robustness of the proposed method benefits from the novel iterative strategy, which has the following two advantages. First, it incorporates the constraint that the image of each shot shares the same magnitude, then an integrated equation from different shots can be set up to improve matrix inversion conditioning. Like two-step PPI methods, the g-factor effect is suppressed and SNR is well maintained. Second, as iteration proceeds, the accuracy of the image and phase error estimation is gradually and simultaneously improved by iterative calibration using the sampled data. This was confirmed in both simulation and in vivo experiments. It is worth noting that although SENSEþCG shares the first advantage, it calculates the phase errors from undersampled data of each shot using PPI, thus can be affected by large gfactor noise amplification from the high reduction factor. For spiral sequences, a larger shot number is desired to reduce off-resonance blurring if high spatial resolution is pursued. In a previous study (31), the SENSEþCG

method was used to achieve high resolution DTI successfully using six-shot acquisitions and a 32-channel coil. Using our proposed POCS-ICE method, the in vivo results acquired with NSS from 6 to 10 and an eightchannel coil can provide images with consistently higher SNR and fewer artifacts than SENSEþCG. Thus a larger shot number can be used. When the number of shot increases further, however, the images may not be solvable with either SENSEþCG or POCS-ICE. The upper limitation of shot number can be increased if an RF coil with more channels is used. For example, we once successfully used a 32-channel RF coil to solve 12-shot DWI data (not shown here), and similar attempts using an 8channel RF coil were not successful. From another point of view, a larger shot number will require a longer acquisition time and will result in extra bulk motion which induces misregistration among different directions. The convergence behavior of the two algorithms was analyzed. POCS-ICE is found to be more stable than SENSEþCG. The divergence of SENSEþCG results from the semiconvergence of CG iteration in the illconditioning problem. Because the reconstruction error decreases at the beginning and then rises as the iteration proceeds, it is difficult to determine the optimal stopping point automatically. Here the stopping iteration number of SENSEþCG was selected manually around the minimum point of the convergence curve. At the

High Resolution Multishot DWI without Navigator

stopping point, the predefined tolerance may not be reached. The semiconvergence problem can be alleviated using regularization algorithms (50). However, the degree of regularization should be carefully chosen to balance SNR and levels of artifacts. Because POCS-ICE shows monotonic decreasing of image error, either a maximum iteration number or a tolerance can be used without manual interventions to stop the iteration. Although it takes a longer reconstruction time than SENSEþCG (approximately 5 times as long), this may not be a problem considering that computation power will continue to increase in the future. Thus POCS-ICE provides a robust solution which can be applied easily in clinical application. Several strategies can be used in the future to accelerate the convergence. Because POCS has the flexibility to incorporate different types of a priori knowledge, the slow convergence can be accelerated by including appropriate constraints, such as region support (33). Due to the point spread function characteristics of spiral imaging, there still exist residual signals outside the true image region, so regional masks cannot be applied directly. A carefully tuned mask which preserves both the region of interest and the residual signals around corners has been tested and shown to be effective. However, it is only able to reduce the reconstruction time by around 10%. Furthermore, an adaptive relaxation factor may also accelerate the convergence (33). Like all the other reconstruction methods, the performance of POCS-ICE should be affected by SNR of acquired signals, which is related to b value, voxel size, etc. In this study, different SNR levels were tested on both simulated and in vivo data using fixed in-plane resolution and slice thickness, and POCS-ICE shows better results compared with SENSEþCG in all conditions. Specifically, the levels of artifacts in the DWI images reconstructed by POCS-ICE (the third column of Figure 6) are very similar, although they have different SNR levels. For higher resolution or thinner slices, the combination of large diffusion gradients and small imaging voxels can result in images with low SNR, which can be very challenging for both reconstruction methods. In such scenarios, more repetitions or longer acquisition times are necessary to compensate for the SNR. A recently introduced compressed sensing method using sharable information among different diffusion directions may help to improve the SNR (51). In conclusion, in this study we proposed a new and robust method, POCS-ICE, to reconstruct multishot diffusion images without use of navigators. The motion induced phase errors and DWI image can be simultaneously solved through the POCS algorithm. Numerical simulations and in vivo experiments demonstrate the superior image quality in terms of SNR and level of artifacts compared with SENSEþCG, especially when the shot number is high relative to the number of coil elements. ACKNOWLEDGMENTS The authors thank the editor and referees for their valuable comments, Drs. Allen W. Song and Trong-Kha

179

Truong at Duke University for their helpful discussions, and Nancy Payne at Oxford University for her assistance in proofreading this manuscript. REFERENCES 1. Basser PJ, Mattiello J, LeBihan D. MR diffusion tensor spectroscopy and imaging. Biophys J 1994;66:259–267. 2. Mori S, Crain BJ, Chacko VP, van Zijl PC. Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Ann Neurol 1999;45:265–269. 3. Muller MF, Prasad PV, Bimmler D, Kaiser A, Edelman RR. Functional imaging of the kidney by means of measurement of the apparent diffusion coefficient. Radiology 1994;193:711–715. 4. Kwee TC, Takahara T, Ochiai R, Nievelstein RA, Luijten PR. Diffusion-weighted whole-body imaging with background body signal suppression (DWIBS): features and potential applications in oncology. Eur Radiol 2008;18:1937–1952. 5. Jezzard P, Balaban RS. Correction for geometric distortion in echo planar images from B0 field variations. Magn Reson Med 1995;34:65–73. 6. Chen B, Guo H, Song AW. Correction for direction-dependent distortions in diffusion tensor imaging using matched magnetic field maps. Neuroimage 2006;30:121–129. 7. Sodickson DK, Manning WJ. Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays. Magn Reson Med 1997;38:591–603. 8. Pruessmann KP, Weiger M, Scheidegger MB, Boesiger P. SENSE: sensitivity encoding for fast MRI. Magn Reson Med 1999;42:952–962. 9. Griswold MA, Jakob PM, Heidemann RM, Nittka M, Jellus V, Wang J, Kiefer B, Haase A. Generalized autocalibrating partially parallel acquisitions (GRAPPA). Magn Reson Med 2002;47:1202–1210. 10. Bammer R, Keeling SL, Augustin M, Pruessmann KP, Wolf R, Stollberger R, Hartung HP, Fazekas F. Improved diffusion-weighted single-shot echo-planar imaging (EPI) in stroke using sensitivity encoding (SENSE). Magn Reson Med 2001;46:548–554. 11. Bammer R, Auer M, Keeling SL, Augustin M, Stables LA, Prokesch RW, Stollberger R, Moseley ME, Fazekas F. Diffusion tensor imaging using single-shot SENSE-EPI. Magn Reson Med 2002;48:128–136. 12. Bammer R, Holdsworth SJ, Veldhuis WB, Skare ST. New methods in diffusion-weighted and diffusion tensor imaging. Magn Reson Imaging Clin N Am 2009;17:175–204. 13. Leuze CW, Anwander A, Bazin PL, Dhital B, Stuber C, Reimann K, Geyer S, Turner R. Layer-specific intracortical connectivity revealed with diffusion MRI. Cereb Cortex 2014;24:328–339. 14. O’Halloran RL, Holdsworth S, Aksoy M, Bammer R. Model for the correction of motion-induced phase errors in multishot diffusionweighted-MRI of the head: are cardiac-motion-induced phase errors reproducible from beat-to-beat? Magn Reson Med 2012;68:430–440. 15. Anderson AW, Gore JC. Analysis and correction of motion artifacts in diffusion weighted imaging. Magn Reson Med 1994;32:379–387. 16. Miller KL, Pauly JM. Nonlinear phase correction for navigated diffusion imaging. Magn Reson Med 2003;50:343–353. 17. Butts K, de Crespigny A, Pauly JM, Moseley M. Diffusion-weighted interleaved echo-planar imaging with a pair of orthogonal navigator echoes. Magn Reson Med 1996;35:763–770. 18. Butts K, Pauly J, de Crespigny A, Moseley M. Isotropic diffusionweighted and spiral-navigated interleaved EPI for routine imaging of acute stroke. Magn Reson Med 1997;38:741–749. 19. Bammer R, Stollberger R, Augustin M, et al. Diffusion-weighted imaging with navigated interleaved echo-planar imaging and a conventional gradient system. Radiology 1999;211:799–806. 20. Porter D, Mueller E. Multi-shot diffusion-weighted EPI with readout mosaic segmentation and 2D navigator correction. In Proceedings of the 12th Annual Meeting of ISMRM, Kyoto, Japan, 2004. Abstract 442. 21. Holdsworth SJ, Skare S, Newbould RD, Guzmann R, Blevins NH, Bammer R. Readout-segmented EPI for rapid high resolution diffusion imaging at 3 T. Eur Radiol 2008;65:36–46. 22. Jeong HK, Gore JC, Anderson AW. High-resolution human diffusion tensor imaging using 2-D navigated multishot SENSE EPI at 7 T. Magn Reson Med 2013;69:793–802. 23. Atkinson D, Counsell S, Hajnal JV, Batchelor PG, Hill DL, Larkman DJ. Nonlinear phase correction of navigated multi-coil diffusion images. Magn Reson Med 2006;56:1135–1139.

180 24. Atkinson D, Porter DA, Hill DL, Calamante F, Connelly A. Sampling and reconstruction effects due to motion in diffusion-weighted interleaved echo planar imaging. Magn Reson Med 2000;44:101–109. 25. Pipe JG, Farthing VG, Forbes KP. Multishot diffusion-weighted FSE using PROPELLER MRI. Magn Reson Med 2002;47:42–52. 26. Liu C, Bammer R, Kim Dh, Moseley ME. Self-navigated interleaved spiral (SNAILS): application to high-resolution diffusion tensor imaging. Magn Reson Med 2004;52:1388–1396. 27. Skare S, Newbould RD, Clayton DB, Albers GW, Nagle S, Bammer R. Clinical multishot DW-EPI through parallel imaging with considerations of susceptibility, motion, and noise. Magn Reson Med 2007;57: 881–890. 28. Truong TK, Chen NK, Song AW. Inherent correction of motioninduced phase errors in multishot spiral diffusion-weighted imaging. Magn Reson Med 2012;68:1255–1261. 29. Uecker M, Karaus A, Frahm J. Inverse reconstruction method for segmented multishot diffusion-weighted MRI with multiple coils. Magn Reson Med 2009;62:1342–1348. 30. Chen NK, Guidon A, Chang HC, Song AW. A robust multi-shot scan strategy for high-resolution diffusion weighted MRI enabled by multiplexed sensitivity-encoding (MUSE). Neuroimage 2013;72:41–47. 31. Truong TK, Guidon A. High-resolution multishot spiral diffusion tensor imaging with inherent correction of motion-induced phase errors. Magn Reson Med 2014;71:790–796. 32. Ma X, Huang F, Zhang Z, Zhang B, Guo H. POCS-ICE: POCS based inherent correction of phase errors for multi-shot spiral DWI. In Proceedings of the 22nd Annual Meeting of ISMRM, Milan, Italy, 2014. Abstract 4444. 33. Samsonov AA, Kholmovski EG, Parker DL, Johnson CR. POCSENSE: POCS-based reconstruction for sensitivity encoded magnetic resonance imaging. Magn Reson Med 2004;52:1397–1406. 34. Gerchberg RW, Saxton WO. A practical algorithm for the determination of phase from image and diffraction plane pictures. Optik 1972; 35:237–246. 35. Youla DC, Webb H. Image restoration by the method of convex projections: part 1 theory. IEEE Trans Med Imaging 1982;1:81–94. 36. Samsonov AA, Velikina J, Jung Y, Kholmovski EG, Johnson CR, Block WF. POCS-enhanced correction of motion artifacts in parallel MRI. Magn Reson Med 2010;63:1104–1110. 37. Jackson JI, Meyer CH, Nishimura DG, Macovski A. Selection of a convolution function for Fourier inversion using gridding. IEEE Trans Med Imaging 1991;10:473–478.

Guo et al. 38. Fessler JA, Sutton BP. Nonuniform fast Fourier transforms using minmax interpolation. IEEE Trans Signal Process 2003;51:560–574. 39. Roemer PB, Edelstein WA, Hayes CE, Souza SP, Mueller OM. The NMR phased array. Magn Reson Med 1990;16:192–225. 40. Jones DK. Diffusion MRI: theory, methods, and applications. Oxford: Oxford University Press; 2010. 41. Newbould RD, Skare S, Bammer R. On the utility of complexaveraged diffusion-weighted images. In Proceedings of the 16th Annual Meeting of ISMRM, Toronto, Canada, 2008. Abstract 1810. 42. Hedley M, Yan H, Rosenfeld D. Motion artifact correction in MRI using generalized projections. IEEE Trans Med Imaging 1991;10:40– 46. 43. Qu P, Zhong K, Zhang B, Wang J, Shen GX. Convergence behavior of iterative SENSE reconstruction with non-Cartesian trajectories. Magn Reson Med 2005;54:1040–1045. 44. Wang YL, Yang JF, Yin WT, Zhang Y. A new alternating minimization algorithm for total variation image reconstruction. Siam J Imaging Sci 2008;1:248–272. 45. Noll DC, Pauly JM, Meyer CH, Nishimura DG, Macovski A. Deblurring for non-2D Fourier transform magnetic resonance imaging. Magn Reson Med 1992;25:319–333. 46. Man LC, Pauly JM, Macovski A. Improved automatic off-resonance correction without a field map in spiral imaging. Magn Reson Med 1997;37:906–913. 47. Jiang H, van Zijl PC, Kim J, Pearlson GD, Mori S. DtiStudio: resource program for diffusion tensor computation and fiber bundle tracking. Comput Methods Programs Biomed 2006;81:106–116. 48. Association NEM. Determination of signal-to-noise ratio (SNR) in diagnostic magnetic resonance imaging. NEMA Standards Publication MS1-2008 2008. 49. Chang H-C, Guhaniyougi S, Chou Y-H, Chen N-K. Interleaved EPI based fMRI improved by integration of multiplexed sensitivity encoding (MUSE) and simultaneous multi-band imaging. In Proceedings of the 22th Annual Meeting of ISMRM, Milan, Italy, 2014. Abstract 2991. 50. Ying L, Liu B, Steckner MC, Wu G, Wu M, Li SJ. A statistical approach to SENSE regularization with arbitrary k-space trajectories. Magn Reson Med 2008;60:414–421. 51. Shi X, Ma X, Wu W, Huang F, Yuan C, Guo H. Parallel imaging and compressed sensing combined framework for accelerating highresolution diffusion tensor imaging using inter-image correlation. Magn Reson Med 2015;73:1775–1785.