Abstract
Using the boxcar representation in the spatial domain and a signal-space representation of its frequency-weighted -space, an iterative prediction method is developed to derive an improved low-resolution phase approximation for phase correction. Compared to the homodyne filter, the proposed predictor is found to be more efficient due to its capability of exhibiting an equivalent degree of performance using a lower number of fractional lines. The phase correction performance is illustrated using partially acquired susceptibility weighted images (SWI). An extension of the predictor into higher frequency regions of phase-encodes in conjunction with a signal-space projection in the frequency-weighted partial k-space is shown to provide restoration of fine structural details of sparse magnitude images. The application of subspace projection filtering is demonstrated using magnetic resonance angiogram (MRA).
1. Introduction
Phase errors in MRI can result from off resonance effects due to imperfections of the static magnetic field or significant changes in the susceptibility within the imaged field-of-view (FOV) [1–4]. The latter can often result in image distortion, particularly along natural and artificial air-filled cavities. When phase errors are present, the complete -space is no longer conjugate symmetric. Coverage of complete -space requires acquisition of all the phase-encode lines and frequency-encode points for a given resolution. Consequently, constrained image reconstruction using either half the set of phase-encode lines or the complete set of partial echoes excluding the dephasing lobe and covering the entire -space, may result in loss of anatomical details. Thus, application of phase correction methods necessitates acquisition of fractional lines in the phase-encode direction, or a set of -space samples spanning the acquired phase-encode locations at a predetermined set of low-frequency locations along the -space frequency-encode direction. The low-frequency locations are such that each echo should contain an equal number of acquired samples on either side of its peak. The set of -space samples at the low-frequency locations on the truncated side of -space form a set of fractional columns. A 2D partial -space is truncated in both directions and contains a predetermined set of both fractional lines and columns. The number of fractional lines and columns is determined by the nature of echo energy distribution in the central portion of the -space. In the energy-concentric regions, the echo magnitudes show a prominent peak at zero frequency and decay rapidly at higher frequencies. The number of fractional lines or columns is chosen such that it lies in the energy-concentric regions.
For real data, phase errors are shown to arise from main field inhomogeneities which is caused by bone-air-tissue interface (also known as susceptibility effects). Apart from phase error, imposition of temporally constrained acquisition, or acquisition of incomplete set of phase-encodes can result in truncation effects [5]. Alternatively, the effects can also be produced from under-sampling in either phase or frequency-encode directions. Often, these artifacts referred to as “Gibbs Ringing” are manifested as striations of sharper image edges in the form of repeated bands parallel to the edge [6]. Alleviation of these artifacts by windowing in the partial -space can lead to blurring and loss of spatial resolution. A second type of artifact arises from performing extrapolation of missing information in -space derived from high pass filtered low signal-to-noise ratio (SNR) images [7]. This type of artifact, often manifested as streak lines, results mainly from noise amplification due to high pass filtering. The solutions for alleviation of phase errors, and simultaneous measures for artifact suppression have been broadly addressed using three distinct approaches. These consist of (1) methods based on low-resolution symmetric data, (2) parametric models, and (3) statistical estimation. Statistical methods use information obtained from the phase of first and higher order autocorrelation, calculated from complex valued phase-distorted image [8]. While the statistical methods are applicable only in situations where the phase variation due to field inhomogeneity is minimal, the usage of higher order autocorrelation results in phase wrapping due to large pixel shifts. The proposed method comprises of a combination of the first two approaches. This results in a more efficient phase correction, since the same performance can be achieved using a lesser number of fractional lines.
In the first group of methods, a narrow strip of data having symmetric phase-encode coverage is selected to provide a low-resolution approximation of the spatial phase variation. With the acquired data consisting of the entire set of positive phase-encodes, the region of this symmetric phase-encode coverage includes the set of fractional lines in the negative phase-encode region. The phase of this low-resolution data is then used for phase correction. In the homodyne method [9], the partial -space data is multiplied with a weighting function called merging filter. The weighted partial -space is then phase corrected using the low-resolution symmetric data. In the iterated version of this approach known as POCS, the missing -space lines in each iteration are replaced from the -space of the phase corrected image [10]. It has been possible to reduce the blurring effects associated with conventional homodyne by applying phase correction to the high pass filtered image [11]. Nevertheless, attempts on prior high pass filtering can result in noise amplification. The limitations of low-resolution approximation for phase correction, applied to in vivo susceptibility weighted image (SWI) are discussed in Section 2.
2. Limitations of Homodyne Phase Correction Using Low-Resolution Approximation
One reason that partial Fourier reconstruction algorithms are not as successful in practice as in theory is associated with distortions introduced during the phase correction. Multiplying an image with a phase correction in the spatial domain is equivalent to circularly convolving the spatial frequency domain with FT of the correction factor. Many of the points in the corrected -space are invalid since they were not actually acquired. This results in areas of increased/decreased signal magnitude in which phase estimation using the low-resolution image breaks down due to rapid spatial frequency variations [12, 13]. Further, low-resolution phase correction methods require phase unwrapping with subsequent high pass filtering, which introduces artifacts in areas with very steep transitions, such as areas near interfaces between parenchyma and bone or air [14].
The effect of phase correction using homodyne method is illustrated using magnitude and phase images of Susceptibility Weighted Image (SWI) shown in Figure 1(a). Magnitude image reconstructed from the offline truncated 2D partial -space with fractional lines in the negative phase-encode direction is shown in the top panel (b). The bottom panel shows the image reconstructed from the same partial -space after application of homodyne phase correction. For ease of comparison, we have chosen regions-of-interest (ROIs) in the frontal, brain stem, and occipital areas. The ROIs are enclosed within colored bounding boxes in (b) to highlight areas of structural distortion after homodyne phase correction. Panels (c) and (d) represent the individual image sections within respective ROIs, prior to and following application of phase correction. Though phase correction serves to improve the venous contrast as evident from the individual ROIs, it introduces loss of structural information in all the three regions.
In this paper, we address the problem of improving the efficiency of low-resolution phase correction by combining the same with a model-based approach. Using this model, signal samples in the unacquired higher phase-encode regions are predicted without introducing artifacts in the reconstructed image. Prior to introducing the details of the proposed method in Section 4.1, we begin with a brief description of the existing model-based approaches in MR literature in Section 3, and the need for adapting the models for phase correction and artifact-free reconstruction from 2D partial -space. Subsequent sections are devoted to the description of theory and mathematical formulation, applicable to the different variants of model-based filtering approaches discussed in this paper.
3. Model-Based Methods for Partial -Space Filling
In model-based methods, extrapolation of missing -space data is accomplished using parametric models. A form of ARMA known as the transient error reconstruction method has been used to model MRI data successfully [15]. Unfortunately, the autoregressive portion of the ARMA model can result in poles in the transfer function, which produce high intensity spikes in the image. It would be much easier to model MRI signals if the corresponding image-domain data contained only sharp features. It is possible to transform raw MRI data such that Fourier transformation produces an image, which consists of predominantly sharp lines [16]. Such a transformed image can be obtained by applying a linear high pass filter. Equivalently, representation of the low-resolution image as a series of boxcar functions [17], leads to the Fourier Transform (FT) of the spatial derivatives modeled as a summation of complex sinusoids in noise. This framework entails the application of linear prediction for extrapolation of partial -space samples weighted using the appropriate frequency terms [18, 19], referred to as the frequency-weighted -space. A possible solution for increasing the efficiency of phase correction is to use a predictor in the frequency-weighted -space along phase-encode direction, followed by low-resolution phase correction. However, the predictor may result in unstable filter weights, leading to spurious spikes or high intensity streaks in the final reconstruction. In view of these difficulties, we propose -space filling methods based on signal-space representation of -space data.
The following sections describe the mathematical background required for understanding the algorithms used for partial -space filling and image reconstruction.
4. Formulation of -Space as a Signal-Space Model
Consider a tissue having proton density distribution , located within a rectangular field-of-view (FOV) extending from cm to cm and cm to cm. For convenience, the FOV is divided into pixels. For two-dimensional imaging, the -direction is frequency-encoded using a gradient-field pulse of height T/cm, and the -direction is phase-encoded with a short gradient-field pulse of duration samples and height T/cm. The MRI signal originating from a single pixel location at is then given by where is the FID signal originating at location . Inverse FT of the local spin-spectral distribution is weighted by , is the sampling interval, echo-time , , , , , and denotes the gyromagnetic ratio. Assuming that the FOV is such that , yields . The composite signal from all locations within the FOV is given by where . For ease of representation, the region of positive and negative phase-encodes is denoted using for , and for .
Using the boxcar representation [17], the spatial derivative of the image in the -direction can be represented using the weighted summation of discrete impulses as a function of . Hence, the 2D-FT of the spatial derivative can be modelled as a summation of finite number of sinusoids, which is linear-predictable to a certain order [18, 19]. From (2), the FT of the spatial derivative can be expressed as
Using the differentiation property of FT, the Transform of the spatial derivative is the same as the Fourier original -space weighted by . The given set of fractional lines of the weighted -space in the negative phase-encoding direction can then be used as input to predict the missing data using where denotes the FIR filter coefficients, and represents the prediction error coefficients. In further discussions, we omit the subscripts and for brevity, and consider only , , and . Partial -space filling is accomplished by scaling each row of the predicted -space by . At each prediction step, the unbiased autocorrelation of the filtered -space is computed and input to the Levinson algorithm [19] for estimating the filter coefficients. In the equivalent signal space model, (4) takes the form where represents th component of the signal vector. denotes the residual component consisting of the sum of noise and the nonregular component of prediction error. Parameters and denote the damping coefficient and spatial frequency, respectively. The relationship between signal model using linear prediction in (4) and signal space representation in (5) is well described in [18, 19]. This enables the signal equation in (5) to be represented in terms of the prediction polynomial roots as where denotes roots of the prediction polynomial
The prediction coefficients are computed using the Levinson algorithm. Both phase correction and magnitude reconstruction only require the knowledge of .
Using the filter coefficients, the next sample in the higher phase-encode line is predicted using (4). In the iterative approach, this estimated sample is reused to form the new input to the Levinson algorithm. In the first prediction step, the error output of the predictor is realized as a regular process of the Wold’s decomposition [20]. However, by inclusion of the predicted sample into the dataset for the succeeding prediction step, the error generated by the predictor would now consist of an additional component that cannot be modelled as a regular process. The steps are summarized in the block schematic shown in Figure 2.
The next section describes how the Levinson algorithm applied to frequency-weighted -space may be used for phase compensation.
4.1. Phase Compensation Using Iterated Prediction
As the algorithm is iterated through a larger number of steps, the nonregular component of prediction error will introduce artifacts in both the magnitude and phase of the reconstructed image. Increase in the number of steps implies prediction in the higher frequency region of the -space, accompanied by a larger nonregular component of the prediction error. Increase in noise power, resulting from the large component of prediction error, introduces artifacts in the magnitude image reconstructed from the predicted -space. Hence, the number of prediction steps has to be limited to the low-frequency region. The phase information in the low-resolution image reconstructed from the central symmetric portion of this predicted -space will contain more information compared to the central symmetric portion of the acquired partial -space. The procedure for phase error compensation using iterated prediction is depicted in Figure 3. The advantages of this approach in comparison to homodyne phase compensation are discussed in Section 5.1. An extension of the iterated prediction, for improved reconstruction of sparse magnitude images from 2D partial -space, is described in the following section.
4.2. Image Reconstruction of Sparse MRI from 2D Partial -Space
As discussed in the previous section, a low-resolution approximation of the predicted -space is used for phase error compensation. The low-resolution approximation is used to avoid prediction of -space data into the high-frequency regions. If prediction is carried out into the higher-frequency region, the magnitude image of the predicted -space is observed to be distorted. This is due to accumulation of nonregular component of the prediction error. For the negative phase-encode half, the high-frequency region includes lines distal from the -space center. The distortion will be reduced if the angle between predicted and true signal vectors is small. The subspace projection method is used to reduce this distortion. In this method, the subspaces are constructed using an equal number of -space samples at any given location along the frequency-encoding axis, in the positive and negative phase-encode directions, respectively. It is assumed that all signal samples in the positive phase-encode half are completely acquired, and hence known apriori. For the projection to work successfully, the angle between true signal vectors that form elements of the positive and negative phase-encode subspaces should be small. Under this condition, it is shown in the succeeding discussion that error compensation using subspace projection will be effective in eliminating the magnitude distortion subject to satisfaction of a condition involving the norm of the projected prediction vector. The condition is such that the norm of the difference between predicted () and true signal vectors () in the negative phase-encode direction is more than the difference between projection of the predicted signal vector () onto the subspace of the positive phase-encoding half and the true signal vector ().
The subspace projection approach is found to work well in sparse MR images. Sparsity means that there are only few significant pixels with nonzero values [21]. In the signal space representation, this is equivalent to a compact representation of the signal vector. Since the prediction errors are reduced due to compact representation, the small angle between predicted and true signal vectors enables subspace projection to reconstruct the magnitude image directly using the predicted -space, with less distortion and restoration of fine structural information. Because of the compactness, the signal vector can be represented using a fewer number of basis functions. Hence, the prediction order can take a value less than . With increase in step size, the -space samples can be predicted with a much lower filter order at each step, without leading to increase in prediction error high enough to cause distortion in the magnitude image. In Section 5.2, the application of subspace projection is illustrated using a sparse GE DQA phantom and MR angiogram.
The dataset consisting of the positive phase-encode steps at a given location along frequency encode axis is first distributed into an dimension subspace. This is achieved by casting the signal samples in the form of a prediction matrix [22] where . By application of Gram-Schmidt orthogonalization procedure to dimensional subspace , we obtain where represents the orthogonal basis vectors spanning the dimension subspace with , and is an upper triangular matrix with unit diagonal elements. In the current approach, the closest approximation of the iteratively predicted data from the acquired set of fractional lines in the negative phase-encoding direction to the subspace in (9) is calculated using the projection theorem [23]. The projected data is then lowpass filtered by multiplying with to estimate the missing -space data. Thus, the compensation method will be effective only in cases where the norm of the predicted signal vector is less than that of the original signal vector.
5. Results
This section includes illustrations using in vivo examples to demonstrate salient features of the FIR filtering approaches discussed in Sections 4.1 and 4.2. Phase compensation is illustrated using SWI, subspace projection method using MR angiogram. All images are acquired on 1.5 T clinical MR scanner (Magnetom-Avanto, Siemens, Erlangen, Germany) with a 12-channel head coil. The in-plane resolution of the images used is 0.25 mm for SWI and 0.7 mm for MR angiogram. The out-of-plane resolution is 5 mm for SWI and 1.4 mm for MR angiogram. Each subsection highlights application of the relevant filtering method as applied to SWI and MR angiogram, respectively.
5.1. Application of Phase Compensation Using Iterated Prediction to SWI
SWI provides a high resolution mapping of the brain’s venous vasculature. It combines the phase and magnitude data to generate an image whose contrast is sensitive to venous blood and iron content. When all data samples from full -space is available, the susceptibility induced contrast variations can be enhanced using a phase mask [24]. However, when the image reconstruction is performed using a partially filled -space, there will be distortion in the magnitude image. The nature of this distortion and its correction using a fixed number of fractional lines was discussed in Section 2. Though homodyne phase correction improves the venous contrast, it was shown to introduce loss of structural information. At the onset of iterated prediction applied to SWI data, we use the same number of fractional lines as in the bottom panel of Figure 1(b). Also, the three regions showing structural distortion after homodyne correction remain the same as in Figure 1(d1)–(d3). In Figure 4, we reproduce these images in the top panel (a) and the ROIs in (b1)–(b3). The bottom panels show the image and ROIs after application of iterated prediction and low-resolution phase compensation. The iterations are limited such that the number of predicted samples is less than , for . For the given resolution of , this corresponds to a low-frequency region in the phase-encode direction. In this example, the low-resolution phase data is reconstructed from the zero-filled centro-symmetric partial-kspace filled using samples within a phase-encode range indexed from to . This phase data is then used for phase compensation by subtracting it from the original phase data reconstructed from the unpredicted -space. The real part of the complex image obtained by combining the magnitude image reconstructed from the original unpredicted -space and the compensated phase data is shown in bottom panel (a) as the compensated image after application of iterated prediction. The ROIs displayed in (c1)–(c3) represent regions within colored bounding boxes of the resultant phase compensated image in bottom panel Figure 4(a). The image restoration achieved using iterated prediction is clearly seen by comparing the respective panels of (b) and (c). It is to be particularly noted that the magnitude image of predicted -space has not been used in this procedure due to the reasons mentioned in Section 4.1.
5.2. Application of Subspace Projection to Sparse MRI
Figures 5 and 6 show sparse MR images reconstructed using iterated prediction, with and without subspace projection. Figure 5 shows the image reconstructed using GE DQA (Daily Quality Assurance) phantom scanned in Axial Plane with 18 cm FOV, 3 mm thick slices with TR of 500 ms, ms with a matrix acquisition, and ±15.63 kHz readout. The raw data is first truncated offline prior to image reconstruction. Panels (a1)–(a3) show a graphical sketch of the 2D partial -space with fractional lines, magnitude of the image reconstructed from zero-filled 2D partial -space, and an ROI to highlight the artifacts generated due to truncation effects. Panels (b) and (c) represent images reconstructed after iterated prediction, prior to and after application of subspace projection, respectively. Since the negative phase-encode half consists of 256 lines, the number of iterations is chosen to exceed 150–200, so as to extend the prediction steps into the higher frequency region. Even though the truncation artifacts are eliminated in the process, it is seen that image reconstructed without application of subspace projection contains artifacts due to large prediction errors in the higher phase-encode locations. The artifacts are clearly visible in the cutout portion of the image shown in (b3). The image reconstructed after subspace projection in (c3) is observed to be artifact-free. Figure 6 illustrates similar effects observed using MR angiogram.
6. Discussion and Summary
This work presents FIR filters for reduction of artifacts and phase errors resulting from truncation of conjugate asymmetric -space. Whereas, truncation artifacts are significant in low-resolution images, phase errors can manifest in both high and low-resolution images. The latter mainly occurs due to magnetic field inhomogeneities. The FIR filters presented in this paper are useful for elimination of both types of artifacts. Recent methods for phase correction use low-resolution approximation of the phase image. The low-resolution approximation is obtained from the symmetric portion of -space, inclusive of a set of fractional lines in the negative phase-encode direction. The drawbacks of such phase correction algorithms include generation of high intensity spikes and intensity distortion leading to loss of fine structural information.
When the number of fractional lines acquired is less compared to the number of lines available in the energy-concentric regions, the resulting intensity distortion can smear intensity variations in areas with steep phase transitions. In this paper, we have developed a linear prediction based -space filling approach to address this problem. In earlier literature relating to application of linear prediction to frequency-weighted -space, prediction of a set of -space samples using a single FIR filter has been shown to result in unstable filter coefficients accompanied by a form of streak artifacts in the reconstructed image. This is true even for the case where predicted lines are used only for estimation of low-resolution phase approximation. That is, when the magnitude reconstructed from the original unpredicted -space is combined with its phase, compensated using the low-resolution phase approximation derived from the symmetric portion of the predicted -space. To address the instability issue, we resort to a form of one-step predictor applied iteratively to fill a predetermined set of -space samples, using a successively updated set of filter coefficients. In this paper, the Levinson algorithm is used for updating the filter coefficients. When the number of predicted -space samples are limited in number, just sufficient to compensate for intensity distortions, the low-resolution phase approximation derived from the symmetric portion of predicted -space has been shown to successfully restore the lost portions of a homodyne phase corrected SWI image. However, if the number of predicted samples extends into higher frequency regions of -space, the nonregular component of prediction errors will begin to introduce distortion in both the magnitude and phase images reconstructed from predicted -space. In an equivalent signal-space representation of frequency-weighted -space, it is shown that this distortion will be reduced only under special conditions. The conditions require (1) the angle between predicted and true signal vectors to be small (2) the angle between true signal vector and projection of predicted signal vector onto a signal subspace representation of the positive phase-encode half to be small. Though both conditions are sensitive to the presence of additive noise, the satisfaction of condition-1 by sparse MR images is construed from a compact nature of their signal-space representation. For sparse MR images, it is shown that magnitude reconstruction using subspace projection of the predicted -space is able to restore the fine structural information, otherwise lost due to -space prediction extended into the higher frequency regions. Though the prediction-based approaches discussed in this paper are extremely sensitive to noise, their main advantages arise from the improved phase correction capability and ability to restore fine structural information.