A higher dimensional homodyne filter for phase sensitive partial Fourier reconstruction of magnetic resonance imaging

Magnetic Resonance Imaging 33 (2015) 1114–1125

Contents lists available at ScienceDirect

Magnetic Resonance Imaging journal homepage: www.mrijournal.com

A higher dimensional homodyne filter for phase sensitive partial Fourier reconstruction of magnetic resonance imaging Joseph Suresh Paul ⁎, Uma Krishna Swamy Pillai Medical Image Computing and Signal Processing Laboratory, Indian Institute of Information Technology and Management-Kerala, Trivandrum India

a r t i c l e

i n f o

Article history: Received 16 July 2014 Revised 29 May 2015 Accepted 20 June 2015 Keywords: Partial k-space Phase errors Incidental phase artifacts Signal losses Ringing artifacts

a b s t r a c t The aim of this paper is to introduce procedural steps for extension of the 1D homodyne phase correction for k-space truncation in all gradient encoding directions. Compared to the existing method applied to 2D partial k-space, signal losses introduced by the phase correction filter are observed to be minimal for the modified approach. In addition, the modified form of phase correction retains the inherent property of homodyne filtering for elimination of incidental phase artifacts due to truncation. In parallel imaging, this new form of homodyne filtering is shown to be effective for minimizing the signal losses, when each of the channel k-spaces is truncated along both phase and frequency-encode directions. This is illustrated with 2D partial k-space for flow compensated multichannel susceptibility weighted imaging. Extension of this method to 3D partial k-space shows improved reconstruction of flow information in phase contrast magnetic resonance angiography with reduced blur and enhanced background suppression. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Scan time reduction is achieved by limiting the acquired points in either phase-encode, frequency-encode or both directions, resulting in a partial k-space. A 2D partial k-space is obtained by constraining the signal acquisition in both phase and frequencyencode directions. Under the assumption that the image is real valued, a major part of the unacquired k-space can be filled by imposing the conjugate symmetry constraint [1]. The remaining parts are filled with zeros (truncated). The zero-filling leads to truncation artifacts in the resultant image. Broadly, the truncation effects appear in the form of blurring and Gibbs ringing artifacts. The latter is manifested as repeated bands appearing parallel to the sharper image edges [2]. For the same degree of truncation, the artifacts are more conspicuous in low-resolution as compared to high-resolution images due to the larger period of repetition of such artifacts in low-resolution images. In phase sensitive Fourier reconstruction, the assumption of a real image is no longer valid. The spatial phase variations can arise due to off-resonance effects [3–5], non-centering of echoes in the readout direction [6], or tissue susceptibility variations [7]. In contrast with the real-valued case, a constrained reconstruction is then only possible by replacing all the unacquired k-space points with zeros. This introduces additional errors in the reconstruction, ⁎ Corresponding author. E-mail address: [email protected] (J.S. Paul). http://dx.doi.org/10.1016/j.mri.2015.06.005 0730-725X/© 2015 Elsevier Inc. All rights reserved.

often referred to as phase errors. These include additional components of blur, signal losses, and truncation artifacts introduced by the presence of phase. The truncation artifacts caused by phase variations appear in the form of localized high frequency intensity changes and signal drop in the vicinity of rapid phase transitions [8]. This arises due to presence of phase wraps, or local non-linearities in spatial phase variation. The characteristic spectral shape of intrinsic phase components posit a narrowband structure, skewed toward higher frequencies. These locally generated non-linear components appear superimposed on the slow varying background phase and the information carrying components of the phase sensitive magnetic resonance imaging (MRI). In homodyne phase correction procedure, the ringing introduced by phase wraps, or incidental phase variations is easily suppressed by application of smooth weighting functions to the partial k-space, prior to image reconstruction. Mitigation of blur calls for compensation of low-frequency phase. If only partially compensated, this leads to additional filter induced signal losses. Application of projection on to convex sets (POCS) based approach does not introduce additional losses as in homodyne method [9,10]. However, it fails to remove the incidental phase artifacts (IPA). Hence, an improved version of the homodyne filter with reduced signal losses would be a desirable solution. In lieu of this, we present a modified form of homodyne phase correction. In the proposed method, the partial k-space is separately weighted along each gradient direction, followed by compensation using phase obtained from low resolution symmetric k-space along the respective gradient direction.

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125

2. Methods

1115

(a)

In an ideal full Fourier acquisition, the k-space data for a real object exhibit Hermitian symmetry. In this case, a partial Fourier acquisition requires knowledge of only one half of the k-space to reconstruct an image of a given resolution and field of view (FOV). However, this advantage of halving the acquisition time has a penalty of SNR reduction by √2. In practice, however, due to the presence of spatially dependent phase shifts, an imaginary part is introduced which destroys Hermitian symmetry. Hence, the image function actually obtained is

y

N

y

(b)

s

y

ð1Þ

N

N

y

(c) y

The homodyne filter corrects these phase variations by using a low resolution image to estimate the phase. In the ideal case, a 1D partial k-space is obtained by truncating the full k-space KF using a boxcar function along the phase-encoding (PE) direction. From signal theory, it is well known that a point-by-point multiplication of KF with a boxcar introduces oscillatory artifacts (Gibbs ringing) in the image. Consider KF to be acquired with N PE lines, with the zero-PE index at j = Ny. In further notations, a PE line is either referred using index j, or using spatial frequency ky. In the absence of complex conjugate symmetry, partial Fourier reconstruction calls for acquisition of an additional n number of fractional lines, leading to a total of Ny + n PE lines. Thus, a 1D partial k-space can be mathematically represented as a truncated form of KF using
 K pk ¼K F ∘ W ky ; n

ay

s

s F ðx; yÞ ¼ sðx; yÞexpð jϕðx; yÞÞ

y

2.1. Theory

ð2Þ

8 f or ky N nΔky ;
 <2 ky ; n ¼ 1 þ ky =nΔky f or −nΔky ≤ky ≤nΔky ; : 0 elsewhere: and
  1 f or ky ≥0; R W y ky ; n ¼ 0 elsewhere

N

N

where W(ky,n) can be either a boxcar function WyR(ky,n), or a boxcar function with linear weighting WyL(ky,n). These are mathematically represented using

y

Fig. 1. (a) Sample 1D phase with a phase-wrap at y = y0, (b) result of convolving a signal with unity magnitude and sample phase in (a) with the IFT of WyR(ky,n), (c) result of convolving a signal with unity magnitude and sample phase in (a) with the IFT of WyL(ky,n).

L Wy

Usage of symbol W without any superscript refers to the general case, where the weighting can be either rectangular or linear. For the definition of W, it is assumed that the PE index is negative for the bottom half of k-space. The reconstructed image using WyR in Eq. (2) will be h i 1 h
 i1 0 −jπ Ny þ n−1 y sin π Ny þ n y=N A A spk ðx; yÞ ¼ s F ðx; yÞ⊗@exp@ N sin½πy=N  0

ð3Þ When used as a kernel for convolution in Eq. (3), the size along y-direction is limited to 2N/(Ny + N). Phase wrapping and effects of field inhomogeneities can be seen as banding effects and wide-scale phase variations. Fig. 1(a) shows a one-dimensional phase ϕ(y) with a phase-wrap at y = y0. Assuming magnitude to be unity across y, panel (b) shows the result of convolving sF(y) = exp(jϕ(y)) with the inverse Fourier transform (IFT) of WyR(ky,n) function as defined in Eq. (3). It is seen that the signal magnitude |spk(y)| drops at y = y0. The ringing in the magnitude is reduced in (c) after convolution with the IFT of WyL(ky,n). The signal loss is dependent on the amount of discontinuity in phase at y0. Usage of a linear weighting in (c) reduces the signal loss at the wrap location. In the presence of a phase-wrap in Fig. 1(a), an oscillatory response is observed in the vicinity of the wrap, when reconstructed

after weighting the k-space with WyR(ky,n). Usage of WyL(ky,n) reduces the oscillations. The oscillatory artifacts observed in the vicinity of a phase-wrap can be treated as one form of IPA. In the presence of rapidly varying field inhomogenities and uneven RF excitations, a component of the high frequency phase (incidental phase) will be superimposed on the phase response as shown in Fig. 2(a). A typical spectrum of such high frequency artifact and its spatial variation are shown in panels (b)–(c). As seen, the central portion has a dip in the incidental phase spectrum, with two prominent side lobes on either side of the frequency scale. Such frequency response patterns can be modeled with non-linear phase components represented using a local polynomial expansion. In Fig. 2(a), if the incidental phase is added to a linear background phase in the vicinity of y0, the high-frequency component of phase (ϕH) can be represented as ϕH ðy0 þ β δyÞ ¼

P X

ai β

i

for

−1 ≤β ≤1

ð4Þ

i¼1

where P represents the order of nonlinearities. Combined background and incidental phase are ϕðyÞ ¼ γϕH ðyÞ þ ð1−γÞϕL ðyÞ

f or

0≤γ ≤1

ð5Þ

and ϕH ðyÞ ¼ 0;

f or

y0 −δy≤y≤y0 þ δy

ð6Þ

1116

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125

respectively. The effect of incidental phase is observed in the reconstructed image as a set of impulses, accompanied by damped oscillatory intensities in their vicinity. These oscillatory responses, when viewed in an image, can be interpreted as the other form of IPA. The oscillatory amplitudes on either side of each impulse may not be exactly symmetric. As before, weighting the k-space with WyL(ky,n) eliminates the artifacts. The homodyne filter [8,11,12] corrects for the low-frequency phase variations by using a low-resolution image to estimate the second phase term in Eq. (5). The phase corrected image is given by

3.14

φ (y)

(a)

0

1.57

0

y0

Ny

N

y ⌢ s F ðyÞ

70

|FT (φ (y))|

(b)

35

0

-NyΔ ky

NΔk

0

ky

y

y

ð8Þ

Since the symmetric phase ϕsym(y) approximates only the lowfrequency component, phase compensation will not reduce the IPA. The 2D partial k-space is obtained by truncating KF in both directions, with n fractional lines along PE, and m fractional samples along frequency encode directions. The fractional lines and sample points are counted with reference to the row and column indices (Ny,Mx), denoting the peak magnitude location. We model the truncation process using two different methods. A partial k-space is first constructed by truncating along PE. The weight matrix used for truncation can be represented as

(c) ð1Þ

0.12 0

h i n
o −1 ¼ Re F T K pk ∘exp −jϕsym ðyÞ

Mx

 X W y kx ; ky ¼ W ky ; n ⊗ δðkx −iΔkx Þ

ð9Þ

φ (y)

i¼−Mx

Taking the IFT,

-0.89

ð1Þ

0

y0

N

y

spk(y)

(d)

y

N

wy ðx; yÞ ¼ Awðy; nÞ∘δðxÞ

where A is a scaling factor resulting from the IFT operation. The image obtained after 1D truncation is ð1Þ

spky ðx; yÞ ¼ s F ðx; yÞ⊗

2.3

ð10Þ

kery ðy; nÞ

ð11Þ

where kery denotes a 1D kernel constructed from the samples within the main lobe of Aw(y,n) along y direction. In the second step, the 2D truncation is modeled using

1.8

ð1Þ

spkxy ðx; yÞ ¼ spky ðx; yÞ⊗ kerx ðx; mÞ ð1Þ Δ spkx ðx; yÞ⊗ kery ðy; nÞ ¼

1.3

0

N

y

N

y

spk(y)

(e)

ð12Þ

The signal loss can be minimized by compensating for the low-frequency phase in spkxy ðx; yÞ. Estimates of low-frequency phase require use of centrally symmetric portions of the k-space. Extraction of the symmetric k-space can be performed using a set of weighting matrices defined separately along PE and frequency-encode directions. Symmetric weighting in the PE direction is obtained as

2

1.4

ð1Þ

Mx

 X R W symy kx ; ky ¼ W sym ky ; n ⊗ δðkx −iΔkx Þ

ð13Þ

i¼−Mx

0.8

0

N

y

N

y

Fig. 2. (a) Nonlinearly generated incidental phase shown superimposed on background linear phase, (b) incidental phase spectrum with dip at zero frequency, (c) spatial variation of incidental phase, (d) result of convolving a signal with unity magnitude and sample phase in (a) with the IFT of WyR(ky,n), (e) result of convolving a signal with unity magnitude and sample phase in (a) with the IFT of WyL(ky,n).

As γ increases, the relative contribution of ϕH increases. Panels (d)–(e) show |spk(y)| reconstructed in the presence of incidental phase, after weighting the k-space with WyR(ky,n) and WyL(ky,n)

where W

R sym


  1 f or −n≤ky ≤n; ky ; n ¼ 0 elsewhere:

Taking the IFT, ð1Þ

wsymy ðx; yÞ ¼ Awsym ðy; nÞ∘δðxÞ

ð14Þ

The low-frequency image obtained after 1D truncation is ð1Þ

ssymy ðx; yÞ ¼ s F ðx; yÞ⊗ sery ðy; nÞ

ð15Þ

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125

where sery denotes a 1D kernel constructed from samples within the main lobe of Awsym(y,n) along y direction. In the second step, the low-frequency phase for compensation is modeled as the phase angle of a complex image given by

1117

Taking the IFT, ð2Þ

    −jπðMx þ mÞx sin½π ðMx þ m þ 1Þx=M wsymy ðx; yÞ ¼ Awsym ðy; nÞ∘ exp M sin½πx=M

ð22Þ

ð1Þ

ssymxy ðx; yÞ ¼ ssymy ðx; yÞ⊗ serx ðx; mÞ ð1Þ Δ ssymx ðx; yÞ⊗ sery ðy; nÞ ¼

ð16Þ

The low-resolution image is then obtained as ð2Þ

ssymy ðx; yÞ ¼ s F ðx; yÞ⊗ wsymy ðx; yÞ In the first step, the low frequency phase is estimated from a narrow portion of the k-space within the strip −nΔky to nΔky. In the second step, this phase is further low-pass filtered using the frequency-band from −mΔkx to mΔkx. In this process, a major portion of the low-frequency component is missed out, leading to errors in the compensated phase. This is interpreted as filter-induced signal loss. In the absence of high-frequency phase components, the filter-induced signal loss can be measured using absolute value of cosine of the phase corrected angle. If the phase-correction is perfect, the cosine values will be ideally unity. With signal loss, the values deviate away from unity. The final phase-corrected image using this method is ⌢ s pk ðx; yÞ

¼ realðspkxy ðx; yÞ∘expð−j∠ssymxy ðx; yÞÞÞ

ð17Þ

ð23Þ

In the proposed 2D-homodyne method, the phase corrected image is treated as



 ⌢ ⌢ ⌢ s pk ðx; yÞ ¼ real s pky exp −j∠ssymy ðx; yÞ þ real s pkx exp −j∠ssymx ðx; yÞ

ð24Þ To compare the signal loss computed using both models, we use a histogram of the absolute values of cosine of compensated phase. For a low-frequency Gaussian shaped spatial phase image, the histograms obtained using both methods are shown in Fig. 3. The steps used for filtering process are summarized in the block schematic in Fig. 4. 2.2. Performance evaluation

To minimize signal losses, we treat the 2D truncation as the sum of two phase-compensated components. The first component is obtained using weighting and compensation along PE. Limiting the span of PE weighting matrix along frequency-encode direction yields m
 X
 ð2Þ W y kx ; ky ¼ W ky ; n ⊗ δðkx −iΔkx Þ ð18Þ i¼−Mx

Taking the IFT, ð2Þ

    −jπ½Mx þ m−1x sin½πðMx þ mÞx=M  wy ðx; yÞ ¼ Awðy; nÞ∘ exp M sin½πx=M  ð19Þ

The second component is obtained by using weighting and compensation along frequency-encode direction, and limiting its span along PE direction. Each kernel is circularly symmetric with main lobes determined by the number of fractional lines (sample points) along respective directions. The kernels are derived by treating points outside the main lobes as zeros. The reconstructed image is now modeled as the sum of two convolutions:

A quantitative evaluation of 1) zero-filling, 2) POCS, 3) 1D homodyne, and 4) 2D homodyne phase correction method is performed by calculating the error in the reconstructed image using error ¼

jjSoriginal −Spartial jj2

where Soriginal and Spartial are the images reconstructed from the full k-space and partial k-space respectively. 2.3. Data acquisition To demonstrate salient features of proposed 2D filtering method, two in vivo examples are used, viz. flow compensated susceptibility weighted imaging (SWI) and low flip angle gradient echo MRI. All images were acquired on 1.5 T clinical MR scanner (Magnetom-Avanto, Siemens,

4000 3500

model 1 model 2

The signal loss in each component of Eq. (20) will be different, and can be minimized by compensating for the low-frequency phase in each component separately. Extraction of the symmetric k-spaces is performed as before, but with the exception that the limits used for summation of the delta functions are now restricted in each direction. For each component in Eq. (20), the symmetric weighting matrices can be defined as

Frequency

3000

ð20Þ

ð25Þ

jjSoriginal jj2

2500 2000 1500 1000 500 0

0

0.2

0.4

0.6

0.8

1

Cosine of phase difference (bins) ð2Þ

m

 X R W symy kx ; ky ¼ W sym ky ; n ⊗ δðkx −iΔkx Þ i¼−Mx

ð21Þ

Fig. 3. Histograms of cosine of compensated phases computed using the two models for 2D truncation. The k-space is constructed using a synthetic image having unity magnitude and Gaussian shaped spatial phase.

1118

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125

KF (2)

Wx ( k x , k y )

(2)

Wy ( k x , k y )

Fourier p reconstruction

s pk y

Symmetricdata data

Symmetricdata data

Fourier reconstruction

Fourier reconstruction

ssymy ( x , y )

ssymx ( x , y )

Re( spk y exp( − j∠ssymy ( x , y )))

Real

Fourier reconstruction

s pk x

Re( spk x exp( − j∠ssymx ( x , y )))

Real

0.5

0.5

Desired image

spk ( x , y ) Fig. 4. Block schematic of the proposed method.

Erlangen, Germany) with a 12-channel head coil. The SWI was acquired at TE = 23 ms, TR = 260 ms with 0.25 mm in-plane and 5 mm out-of-plane resolution, and the low flip angle gradient echo image (GREI) at TE = 16.4 ms, TR = 50 ms with 0.25 mm in-plane and 4 mm out-of-plane resolution. 3D homodyne reconstruction is illustrated using phase contrast magnetic resonance angiography (PC-MRA) datasets of two volunteers. Both datasets were acquired on 1.5 T clinical MR scanner (Magnetom-Avanto, Siemens, Erlangen, Germany) with a 6-channel head coil at TE = 9 ms, TR = 56.70 ms and velocity encoding of 10 cm/s. For balanced four-point acquisition, 4 partitions of data were acquired, with each k-space volume of dimension 512 × 167 × 64. 3. Results 3.1. Numerical simulation The IPA artifacts originate from the high-frequency phase components with a spectral shape as shown in Fig. 2(b). A 2D

phase having a similar spectral profile in both x and y directions is simulated using a set of concentric disks with uniform phase within each segment, and abrupt transition of phase across segments. It is to be noted that by using circular disks, the 2D spectral shape becomes circularly symmetric. However, circular symmetry is not actually need in practice. Numerical phantoms shown in Fig. 5 consist of (a1) magnitude image, (b1) low frequency phase image, (c1) combination of low and high frequency phase components. The magnitude image in Fig. 5(a1) is constructed using a circular object with radius of R0 and uniform intensity. The spatial low-frequency phase component is simulated using 2

2r ϕlow ¼ πð 2 −1Þ; R0

0b r ≤R0

ð26Þ

Panels (b1–c1) show images of low-frequency phase, and combined low and high-frequency phase. Three sample images are synthesized by using the magnitude alone, magnitude combined with low-frequency phase, and magnitude combined with high and

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125 1

(a1)

3

(b1)

2

0.8

1119 3

(c1)

2

1

1

0

0

0.6 0.4

-1

-1

0.2

-2

-2

0 1

-3

-3

(a2)

1

(b2)

1

(c2)

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

Fig. 5. Simulation using numerical Phantom. (a1) Magnitude image, (b1) low frequency phase image, (c1) phase image with combination of low and high frequency components, (a2) magnitude of image reconstructed from the k-space of magnitude image in (a1) after 2D truncation, (b2) magnitude of image reconstructed from the k-space of magnitude image in (a1) in combination with phase image of (b1) after 2D truncation, (c2) magnitude of image reconstructed from the k-space of magnitude image in (a1) in combination with phase image of (c1) after 2D truncation.

(a1)

(b1)

(c1)

(d1)

(a2)

(b2)

(c2)

(d2)

(a3)

(b3)

(c3)

(d3)

(a4)

(b4)

(c4)

(d4)

Fig. 6. Comparison of various phase correction methods applied to numerically simulated data. The magnitude image is constructed using a circular object with uniform intensity. The spatial low frequency phase is simulated using Eq. (21) and the high frequency component is simulated using a set of concentric disks with uniform phase within each segment and abrupt transition of phase across segments. The amount of incidental phase variation is numerically controlled using Eq. (5) for different values of the high frequency boost factor γ. The images shown correspond to row-wise increasing values of γ indicating higher levels of incidental phase variation. Column-wise panels illustrate different reconstruction methods. (a1)–(a4) images reconstructed using zero-filled partial k-space, (b1)–(b4) 1D homodyne, (c1)–(c4) modified homodyne, and (d1)–(d4) POCS method.

1120

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125

low-frequency phase. A k-space is first constructed for each image. 2-D partial k-spaces for each case are then constructed by truncating along both phase and frequency encoding directions. For a 256 × 256 resolution, the partial k-space is obtained by choosing n = 30 fractional lines in the phase-encode direction, and m = 60 frequency points in the frequency-encode direction. The unfilled portions are replaced with zeros. The magnitudes of partial Fourier reconstructed images are shown in Fig. 5 (a2–c2). It is seen that IPA is observed only for the case where the high frequency components are present. The remainder of the two cases exhibits only the Gibbs ringing artifacts. Varying degree of incidental phase component is introduced into the simulation, as per the combined phase model in Eq. (5). The images reconstructed using zero-filled partial k-space for different γ values are shown in panels (a1)–(a4) of Fig. 6. Presence of IPA is seen to increase with the boost factor γ. In addition to this, the impulses observed at locations of phase jumps (see Figs. 1(b–c) and 2(d–e), Section 2.1) are seen as ring like structures superimposed in the magnitude image. Panels (b)–(d) show images reconstructed using 1D homodyne, modified homodyne and POCS filters respectively. It is seen that with increasing value of γ, the reduction of IPA is more obvious in both homodyne methods. As an illustration of signal loss reduction outlined in Section 2.1, an image having a localized phase-wrap, superimposed on a lowfrequency background phase is shown in Fig. 7(a). Panels (b–c) illustrate images of cosine of phase-corrected angles computed using methods 1 and 2. Bright regions indicate locations that are fully compensated. It is clearly seen that the number of pixels which are fully compensated is more in method-2. 3.2. Illustration of signal loss and IPA in real data The acquired k-space data are first truncated offline in both phase and frequency-encode directions. For a 1024 × 1024 resolution, approximately 30% of the full k-space is retained by choosing n = 40

fractional lines in the PE direction, and m = 80 samples in the frequency-encode direction. The remaining portion of k-space is zero-filled. The images obtained using Fourier reconstruction applied to the zero-filled k-spaces of SWI and GREI data are shown in Fig. 8(a) and (b) respectively. The top part of each panel depicts a graphical representation of the k-space from which the images are reconstructed. Left panels correspond to the full k-space, and right panels correspond to the partial k-space. For the SWI example, effect of phase errors is seen in the magnitude image in the form of blurring and regions with signal losses. For ease of comparison, we have chosen regions-of-interest (ROIs) around the frontal and occipital region. The insets of ROIs highlight areas with signal losses in the magnitude image. For the GREI example, IPA is observed in the magnitude image reconstructed from partial k-space. The insets show ROIs highlighting regions affected by IPA. 3.3. Phase correction in 2D partial k-space Fig. 9 illustrates the phase corrected images of low flip angle GREI using POCS, 1D homodyne, and 2D homodyne methods. The leftmost panel shows the zero-filled image without phase correction. The yellow contours highlight regions with IPA and signal loss due to truncation. The result of application of POCS method is shown in panel (b). As discussed in Section 2, IPA is not removed using POCS method. It is also seen that POCS does not compensate for signal loss due to truncation. Compensated images obtained using 1D homodyne and 2D homodyne methods are shown in panels (c)–(d). While the homodyne method is able to remove IPA, it introduces an additional component of signal loss in the region enclosed within the yellow contour. The 2D homodyne method is able to compensate for this additional signal loss, together with the IPA. Hence the proposed method is a more suitable choice for phase compensation in the presence of IPA. The relative signal loss for each method is indicated in the bar graph. Here, the signal losses are computed by finding the sum total of absolute intensity differences between a reference image reconstructed from

(b)

1

0.8

0.6

0.4

(a)

3 2

0.2

1

0 1

0

(c)

-1

0.8

-2

0.6

-3

0.4

0.2

0 Fig. 7. Comparison of signal loss reduction using numerical simulation. (a) A low-frequency phase image with a localized phase wrap, (b) image showing cosine of phase corrected angle using method-1, (c) image showing cosine of phase corrected angle using method-2.

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125

Ny

n

Ny

Ny

m Mx

Mx

1121

n

Ny

m Mx

Mx

(1) (1) (2)

C

(2)

Fig. 8. Illustration of signal losses and IPA in images obtained using Fourier reconstruction applied to the zero-filled k-spaces of (a) SWI and (b) GREI data. The top part of each panel depicts a graphical representation of the k-space from which the images are reconstructed. Left panels correspond to the full k-space and right panels correspond to the partial k-space. Middle and bottom panels show magnitude and phase images respectively. Insets in (a) show regions with signal losses. Insets in (b) show regions with IPA.

the full k-space and the image reconstructed using each method. The relative signal loss is then computed as the ratio of this difference to the sum of intensities in the reference image. Severe signal losses are observed for the flow compensated SWI image, after 2D offline truncation. SWI images contain information in the high-frequency phase components, which are required for further processing. SWI processing is performed by deriving a mask from the high-frequency phase, which when applied to the magnitude image,

1

0.7 (c)

1

(b)

1

0.8

0.8

0.8

0.6

0.6

0.6

0.4

0.4

0.4

0.2

0.2

0.2

0

0

0

(c)

0.6

Relative signal loss

(a)

leads to improved contrast of iron carrying regions [13]. In the case of partial k-space acquisition of SWI, the task of high-pass filtering is automatically performed in the phase compensation steps described in Section 2.1. In the context of partial k-space processing of SWI, therefore, a second high-pass homodyne filtering of phase is not required. As discussed in Section 2.1, the filter induced signal losses will be minimal if the reconstruction is performed using the proposed 2D homodyne filter. It is to be noted, however, that as a consequence of

1

(d)

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.5 0.4 0.3

(d) (a)

(b)

0.2 0.1 0

Fig. 9. Application of phase correction methods to low flip angle GREI data. Top left panel displays reference image reconstructed from full k-space. (a) Image reconstructed from zero-filled k-space, (b) image reconstructed using POCS phase correction, (c) image reconstructed using 1D homodyne phase correction, (d) image reconstructed using 2D homodyne phase correction. Yellow contours show areas exhibiting signal losses. Bar graph shows relative signal losses in each method.

1122

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125

Phase encode direction

1,2,3,…………..….…109 1,2,3,……………………………….…167

the phase compensation step in SWI, a remnant high-frequency phase is retained. Therefore, the signal losses are measured from the magnitude of the phase corrected image, rather than the real part as in Eq. (24). Representative SWI image reconstructed without truncation (full k-space) is shown in full on the left side of Fig. 10. Three different ROIs from the frontal, brain stem and occipital regions are labeled as (a), (b) and (c). The array of panels on the right illustrates the three ROIs of images reconstructed from 2D partial k-space. Column-wise panels depict zoomed versions of each labeled region (a–c). Row-wise panels from top-to-bottom are arranged in the order of zero-filled, POCS, 1D homodyne and 2D homodyne respectively. For the zero filled reconstruction in top row, the effect of truncation appears in the form of blurring and signal losses. Application of 1D homodyne filtering leads to significant signal losses in these areas. As before, 2D homodyne reduces the amount signal losses.

1,2,3,……………………………....,326

1,2,3,…………………………………………………...,512 Frequency encode direction

Fig. 11. Illustration of 3D partial k-space.

3.4. Phase correction in 3D partial k-space Each PC-MRA dataset consists of 4 partitions with flow encoding first moments altered in pairs. Balanced four point method is used to reconstruct the magnitude flow images along x, y and z directions [14]. The 3D partial k-space for each partition is obtained by offline truncation of the respective k-space. For 512 × 167 × 64 array, 3D

(a)

partial k-space is synthesized by zero filling each k-space partition as indicated in Fig. 11. Each 3D partial k-space is phase corrected using 3D homodyne phase correction, with steps identical to the 2D homodyne discussed in Section 2.1. In addition to the fractional samples and PE lines of 2D

(a1)

(b1)

(c1)

(a2)

(b2)

(c2)

(a3)

(b3)

(c3)

(a4)

(b4)

(c4)

(c) (b)

Fig. 10. Application of phase correction methods to flow compensated SWI data. SWI image reconstructed without truncation (full k-space) is shown in full on the left side. Three different ROIs from the frontal, brain stem and occipital regions are labeled as (a), (b) and (c). Array of panels on the right illustrates the three ROIs of images reconstructed from 2D partial k-space. Column-wise panels depict zoomed versions of each labeled region (a–c). Row-wise panels from top-to-bottom are arranged in the order of zero-filled, POCS, 1D homodyne and 2D homodyne reconstruction respectively.

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125

homodyne, fractional slices along z-direction are also considered for volume reconstruction. Consequently, we obtain three components ⌢ s pky , ⌢ s pkx and ⌢ s pkz in Eq. (20). Finally, the phase correction is also separately performed for the three components. Balanced four-point method is now applied to these phase corrected images. The phase difference speed is obtained by maximum intensity projection of the combined x, y and z flow images. This is illustrated in Fig. 12, where row-wise panels correspond to each volunteer dataset. Column-wise panels show images reconstructed using zero-filled 3D partial k-space, 3D POCS reconstruction, and proposed 3D homodyne filter. Improved background suppression and reduced blur levels are clearly evident in the image reconstructed using the proposed method.

3.5. Application to parallel imaging In this section, we illustrate application of the 2D homodyne filter to multi-channel partial k-spaces. To decouple the effects of undersampling and aliasing artifacts [15], we show examples in which full FOV channel images are reconstructed from the individual channel partial k-spaces without undersampling. When combined using the root sum-of-squares (rSoS) method, the 2D phase-corrected channel images result in reduced blur and signal losses. Individual phasecorrected channel images, and resulting rSoS image reconstructed from 2D truncated channel k-spaces are shown in Fig. 13. The truncation is performed for a k-space fraction of 0.4. This means that only 40% of the total k-space is considered as acquired. The remaining portions are filled with zeros. The first four column-wise panels show images from the individual coils. The fifth column shows the rSoS image. Row-wise panels (a)–(d) show images reconstructed using zero-filled channel k-space, POCS, 1D homodyne and 2D homodyne methods respectively. ROI enclosed in red bounding box shows regions exhibiting IPA in the zero-filled and POCS reconstructed images. The yellow arrows point to regions exhibiting signal loss. Filter induced signal losses are

1123

more pronounced in channel images reconstructed using 1D homodyne as compared to 2D homodyne filter. 3.6. Performance comparison Fig. 14 shows the dependence of error on the percentage of acquired k-space for reconstruction using zero-filling, POCS filter, 1D homodyne filter, and the proposed 2D homodyne filter. It is seen that the error decreases with increase in the percentage of acquired k-space. The errors are minimum for the 2D homodyne filter. Similar performance can be seen for 3D reconstruction applied to phase contrast magnetic resonance angiogram in Fig. 15. It is seen that the error decreases with increase in the fraction of acquired k-space, and minimum for the 3D homodyne method. The errors are computed using maximum intensity projection (MIP) images, with an MIP image reconstructed from the full k-space as reference. 4. Discussion and summary Fourier analysis has revealed that the phase of a complex signal carries information that could potentially be used for post-processing and filtering applications. In MRI, this is best illustrated with scans using SWI and PC-MRA. In SWI, the high-frequency phase information is used to derive weights, which when applied to the magnitude, enhance contrast of venous structures. Likewise in PC-MRA, phase difference of acquisitions using bipolar gradients with alternately reversed polarities is combined with the magnitude information to derive directional flow images. Efforts in reduction of scan times often dictate acquisition of a limited region of the k-space. This is achieved by truncating the k-space in phase-encode, frequency-encode and slice-select directions. In phase-sensitive imaging, this form of truncation often leads to loss of information, and artifacts originating from spurious phase components, mutilating the information required for further post-processing. Non-linearity of field inhomgeneities can manifest in the form of narrowband frequency components in phase, mostly restricted to localized regions in the FOV. The characteristic spectral shape of such

(a1)

(a2)

(a3)

(b1)

(b2)

(b3)

Fig. 12. Application of 3D homodyne phase correction to PC-MRA volunteer data sets. Row-wise panels correspond to each volunteer dataset. Column-wise panels show images reconstructed using zero-filled 3D partial k-space, 3D POCS reconstruction, and proposed 3D homodyne filter.

1124

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125

(a1)

(a2)

(a3)

(a4)

(a)

(b1)

(b2)

(b3)

(b4)

(b)

(c1)

(c2)

(c3)

(c4)

(c)

(d1)

(d2)

(d3)

(d4)

(d)

Fig. 13. Image reconstruction using 2D truncated coil k-spaces. (a1)–(a4) Individual channel images reconstructed using zero-filled k-space, (b1)–(b4) individual channel images reconstructed using POCS, (c1)–(c4) individual channel images reconstructed using 1D homodyne, (d1)–(d4) individual channel images reconstructed using 2D homodyne method. (a) Combined SoS image using zero-fill, (b) combined SoS image using POCS, (c) combined SoS image using 1D homodyne (d) combined SoS image using 2D homodyne method.

intrinsic phase components posits a narrowband structure, skewed toward higher frequencies. These locally generated non-linear components appear superimposed on the slow varying background phase and the information carrying components distributed across the frequency spectrum. The effect of non-linearities is manifested locally, as rapid phase transitions. It is to be pointed out that such local transitions can also be generated due to the presence of phase wraps. Truncation of k-space will generate rapid intensity changes at the locations of these transitions, accompanied by oscillatory artifacts in their vicinity. This is visually interpreted as IPA. We have shown realistic example of such artifacts in a low flip angle GRE image, and multi-channel SWI. Using numerical simulation, we have shown that a

4.5

x 10

-3

Zero filling POCS Convential Homodyne Extended Homodyne

4 3.5

Error

3 2.5 2 1.5 1 0.5 0 30

35

40

45

50

55

60

65

70

% of k-space Fig. 14. Reconstruction error measured with reference to the image reconstructed from full k-space versus percentage of acquired k-space.

nonlinear polynomial model of phase in a local neighborhood can mimic the effects seen in the presence of such artifacts. At a glance, these can often be mistaken for Gibbs ringing artifacts, also a result of k-space truncation. The latter is caused by the under-sampling of high spatial frequencies at sharp boundaries in the object being imaged. However, a major difference between Gibbs ringing artifacts and IPA is that the former appears in the vicinity of physical edges of objects in the FOV, while ringing due to incidental phase occurs near to the boundaries of rapid phase transitions. In an effort to correct for these artifacts, we have explored possible solutions and compared the performance of existing phase correction methods, mainly, POCS and homodyne filters. As opposed to POCS, homodyne filtering involves prior weighting of partial k-space. Application of smooth weighting functions has been shown to effectively reduce the IPA. However, pre-weighting penalizes reconstruction performance by introduction of additional signal losses. In this paper, we have demonstrated that the modified form of 2D homodyne filtering significantly improves the reconstruction by way of reducing the extra signal losses and blur. The effect of blur is prominently seen in MIP images obtained from phase-corrected PC-MRA. The performance worsens when the truncation is performed additionally in the slice-select direction. In this situation, the modified form of 3D homodyne filter has been shown to suppress the background effects and reduce the blur considerably. Signal losses due to truncation are mainly observed at locations of phase wraps and rapid intensity changes generated by incidental phase components. Pre-weighting increases the extent of signal losses, particularly, in higher dimensional truncation. This constrains the use of existing homodyne filter for pre-processing and reconstruction of phase-sensitive MRI. In this context, the proposed filter finds potential application for pre-processing of multi-channel MRI. From a theoretical perspective, the reduction of filter induced signal losses is justified by the unique nature of phase compensation in our filtering method. In a conventional form of 2D phase compensation (method-1 in Section 2.1), the phase used for compensation is synthesized from a cascaded model for formulation of a 2D centrally-symmetric k-space. The estimated

J.S. Paul, U. Krishna Swamy Pillai / Magnetic Resonance Imaging 33 (2015) 1114–1125

1125

0.32 Zero-filling Modified homodyne POCS

0.3 0.28 0.26

Error

0.24 0.22 0.2 0.18 0.16 0.14 0.12 0.15

0.16

0.17

0.18

0.19

0.2

0.21

0.22

0.23

0.24

0.25

Fraction of k-space Fig. 15. Reconstruction error measured from maximum intensity projection of image volumes reconstructed using 3D zero-fill, POCS and homodyne methods versus fraction of acquired k-space.

phase is insufficient to fully compensate the low-frequency phase. In the modified form of our filter, the phase-compensation and pre-weighting are performed independently along phase-encode and frequency-encode directions. This form of compensation has proven to be effective in significantly reducing the additional losses due to pre-weighting. Though our filter does not completely eliminate the signal losses, it is able to decrease it to an extent comparable to that of zero-filled reconstruction and POCS. In summary, our filtering method would be an optimal choice for concomitant reduction of blur and IPA, especially when phase-sensitive MRI is acquired partially in two and three dimensions. 5. Conclusion In this study, we have shown how to extend the principle of homodyne phase correction for two and three dimensional partial k-space acquisition. Homodyne phase correction has been shown to be useful for removal of artifacts due to phase wraps and incidental phase variations. The main advantage of the modified form of 2D and 3D homodyne filtering is in reduction of filter induced signal losses as compared to the existing approach. Application of this extended form of homodyne filtering to 3D PC-MRA has shown to result in reconstruction with reduced extent of blur. The proposed method can serve as a practical solution to obtain artifact-free images without signal losses in the presence of local high-frequency phase components due to field inhomogeneities. Acknowledgement The authors wish to thank the Department of Science and Technology (DST SR/S3/EECE/0107/2012) of India for scholarship and operating funds. We also thankfully acknowledge Dr. Jaladhar Neelavalli, Assistant Professor and Uday Bhaskar Krishnamurthy, Research

Assistant, Department of Radiology, Wayne State University, Michigan for providing the in-vivo MR data. We would like to extend our thanks to Dr. Sairam Geethanath, and Ms. Rashmi reddy, Biomedical Research Centre, Dayanand Sagar Institution for providing phase contrast magnetic resonance angiogram data. References [1] Liang ZP, Boda FE, Constable RT, Haacke EM, Lauterbur PC, Smith MR. Constrained reconstruction methods in MR imaging. Rev Magn Reson Med 1992;4:67–185. [2] Wood ML, Henkelman RM. Truncation artefacts in magnetic resonance imaging. Magn Reson Med 1985;2(6):517–26. [3] Smith TB, Nayak KS. MRI artifacts and correction strategies. Imaging Med Rev 2010;2(4):445–57. [4] Dietrich O, Reiser MF, Schoenberg SO. Artifacts in 3-tesla MRI: physical background and reduction strategies. Eur J Radiol 2007;65:29–35. [5] Schomberg H. Off-resonance correction of MR images. IEEE Trans Med Imaging 1999;18(6):481–95. [6] Hou Ping, Hasan Khader M, Sitton Clark W, Wolinsky Jerry S, Narayana Ponnada A. Phase-sensitive T1 inversion recovery imaging: a time-efficient interleaved technique for improved tissue contrast in neuroimaging. AJNR Am J Neuroradiol 2005;26:1432–8. [7] Carlsson A. Susceptibility effects in MRI and H MRS. Doctoral Thesis SE, Sweden: Dept. Radiation Physics, University of Gothenburg; 2009. [8] Noll DC, Nishimura DG, Macovski A. Homodyne detection in magnetic resonance imaging. IEEE Trans Med Imaging 1991;10:154–63. [9] Pauly J. Partial k-space reconstruction”. http//users.fmrib.ox.ac.uk/~karla/reading_ group/lecture_notes/Recon_Pauly_read.pdf. [10] Chen J, Zhang L, Zhu Y, Luo J. MRI reconstruction from 2D partial k-space using POCS algorithm. Int. Conf. on Bioinformatics and Biomedical Engineering, Beijing; 2009. p. 1–4. [11] Margosian P. Faster MR imaging-imaging with half the data. Health Care Instrum 1986:195–7. [12] MacFall SR, Pelc NJ, Vaurek RM. “Correction of spatially dependent phase shifts for partial Fourier imaging”. Magn Reson Imaging 1988:143–55. [13] Haacke EM, Xu Y, Cheng YC, Reichenbach JR. Susceptibility weighted imaging (SWI). Magn Reson Med 2004;52(3):612–8. [14] PeIc NJ, Bemstein MA, Shimakawa A, Glovem GH. Encoding strategies for threedirection phase-contrast MR imaging of flow. J Magn Reson Imaging 1991;1:405–13. [15] Blaimer M, Breuer F, Mueller M, Heidemann RM, Griswold MA, Jakob PM. SMASH, SENSE, PILS, GRAPPA how to choose the optimal method. Magn Reson Imaging 2004;15:223–36.