Application of autoregressive modelling in magnetic resonance imaging to remove noise and truncation artifacts

Magnetic Resonance Imaging. Vol. 4 pp. 251-261. Printed in the USA. All rights reserved.

1986 Copyright

0730-725X/86 503.00 + .OO 8 I986 Pcrgamon Journals Ltd.

l Technical Note APPLICATION OF AUTOREGRESSIVE MODELLING IN MAGNETIC RESONANCE IMAGING TO REMOVE NOISE AND TRUNCATION ARTIFACTS M.

R.

SMITH,*

S.

T. NICHOLS,*

R.

M.

HENKELMAN~

AND

M. L. WOOD?

*Department of Electrical Engineering, University of Calgary, Calgary, Alberta, Canada, T2N lN4; TDepartment of Medical Biophysics, University of Toronto, 500 Sherbourne Street Toronto, Ontario, Canada M4X 1K9 Magnetic resonance imaging data is conventionally reconstructed using two dimensional discrete Fourier transforms. However, there is growing interest in other types of spectral estimation which minimize noise and artifacts due to truncated data. This note presents preliminary results-showing the improvement obtainable using a modified autoregressive model, the Transient Error method. Keywords: Magnetic resonance imaging, Maximum entropy methods, Autoregressive

lobes, which give rise to artifacts in the final image and may provide spurious anatomical detail. Wood and Henkelman have shown that applying amplitude and phase windows to the MR data before the DFT reduces the sidelobes but with a loss of spatial resolution.2 The trade off between sidelobes and resolution for various window types used in DFT applications is discussed by Harris.4 Harris explains that a physical insight into the appearance of the sidelobes can be obtained by realizing that the DFT places a cyclic nature on the data. If this cyclic data shows no discontinuities, then it can be described by the orthogonal basis set of the DFT (sines and cosines) and the transform of the finite data set is equivalent to the transform of the infinite set. However, discontinuities indicate that frequencies other than the basis set are present in the data and the Fourier transform will be spread across the entire frequency band. This is known as Spectral Leakage and is responsible for the sidelobes described by Wood and Henkelman. Autoregressive modelling is an alternative method for spectral estimation. The technique attempts to determine a parametric model of the infinite data set based on a short record. Provided the model fits, this technique will provide excellent spectral information with good noise rejection from short data records with no windowing artifacts. Further, since short data records can be fitted, reconstructions may be made with a small number of phase-encoded signals. This is impor-

The Kumar, Welti and Ernst technique for magnetic

resonance (MR) imaging involves sampling a spin echo in the presence of a linear magnetic field gradient. Information within the image plane, but perpendicular to the field gradient, is incorporated into the spin echo as phase encoding.’ The image is reconstructed by applying a two dimensional Fourier transform to a series of spin echoes with differing phase encoding. Wood and Henkelman have discussed the problems associated with reconstruction using the Discrete Fourier Transform (DFT) on a finite or truncated segment of data.2 These result in artifacts in the final image that introduce uncertainty in the discrimination of anatomical detail, particularly in the phase-encoding direction. These artifacts can be reduced by low-pass filtering or windowing but with an associated loss of resolution. This note gives the preliminary results of a new technique of image reconstruction which inherently reduces both the noise and truncation effects. It is based on a modified autoregressive (AR) modelling algorithm called the Transient Error method which estimates the Fourier transform of the signal and not its power spectral density function (PSDF).3 When a reconstruction of MR data is performed using the DFT on a finite segment of data, the resulting image is the required image convolved with the Fourier transform of the window function applied to the MR data. Since the transform is typically a sine function, the final image will show evidence of sideRECEIVED~/~~/~~;ACCEPTED

modelling.

11/11/85. 257

258

Magnetic

Resonance

Imaging

PREDICTION ERROR SEQUENCE

l Volume 4, Number

3.1986

MAGNETIC RESONANCE DATA

En

’ H,R(Z)=

I---

I A(Z)

AR MODEL FILTER

TRANSIENT ERROR SEQUENCE H,(Z)=A(Z)

En

X.

4

----I

Fig. 1. Schematic representation of the spectral decomposition of magnetic resonance data using the Transient Error method. tant in anisotropic three dimensional volume imaging where as few as 16 slices may be imaged with a high degree of slice to slice interference due to truncation. A small number of phase encoding gradients will also reduce the overall measurement time in two dimensional imaging. Since the MR data is not expected to be well parameterized by a strict AR model, the modified AR Transient Error method is used.3 Sibiski et af5 have reported using modelling of conventional NMR high resolution spectroscopy using a Maximum Entropy criterion. Their results show considerable noise suppression on highly spiked data. Similar signal to noise enhancements have been reported by Barkhuijsen et UP who modelled “P FID data using an AR model. The following is a brief description of the theory behind the Transient Error method. Consider the samples x, . . . x, of a single echo as a subset of the infinite data sequence x,, x2. . . . The autoregressive technique makes a parametric fit to these data points using a p th order AR model P

xn= -

x

UiX,_i

+

E,

(1)

i-l

where ai are the AR coefficients and E, is the prediction error t, = x, - x,

f:

H(z)

=- 1

’

=

1+ f

&z-i

(4)

A(z)

i-l

This system is shown schematically by the first part of Fig. 1. There are a number of different algorithms, for example Burg’s,’ for determining the AR coefficients from a finite data record.‘*‘,’ These can be used to determine the power spectral density function (PSDF) of the original data sequence.’ However, the PSDF approach removes all phase information and in addition, since the AR model is an ALL-POLE model, can lead to highly spiked spectra. The Transient Error method goes one step further than the conventional AR techniques and overcomes this problem. It is essentially a systematic way of solving a finite difference equation with initial values. It is a modified AR method that estimates the Fourier transform and not the PSDF. It is computationally more efficient than the well known extended PRONY method,’ but is less stable in high noise environments. Since the spin echo data can be modelled as the output of a linear filter excited by the prediction error noise en, we can write in the frequency domain

(2)

The predicted value of x, is given by

x, = -

function of this filter is

aix,i

i-l

This is basically the same as considering that the spin echo data x,, x2 . . . has been produced by feeding the prediction error signal t, into a filter described by the difference equation (1). In the z domain, the transfer

FT(%) FT( x,) = FT(aJ

(5)

where FT denotes the Fourier transform and a,, the AR coefficients determined by one of the AR algorithms, with a, set to 1. The unknown numerator E, can be obtained by applying an inverse filter to the data sequence as shown by the second part of Fig. 1. This approach is called the Transient Error method because if the data sequence being modelled is a strict p th

Application of autoregressive modelling l M.R.

SMITH

259

ETAL.

P = IO

P=l5

P=20

AR

Model

I

Discrete

Fourier

Transform

Fig. 2. This shows the reconstruction of the image of a water filled phantom containing tubes of doped gel using a Transient Error model of orders p - 10, 15, 20. The modelled left hand side shows a reduction of the artifacts in certain areas of the object when compared to the normal FFT reconstructed right hand side. As the order of the model increases, the artifacts and noise are reduced. However, an occasional high intensity spot may be seen because of overfitting in the AR portion of the algorithm.

260

Magnetic Resonance Imaging 0 Volume 4, Number 3, 1986

order AR mode, then the output of the prediction error filter will vanish after the p th output. Once the coefficients a, and e,,have been obtained, it is a simple matter to determine the spectrum of the original data sequence using Eq. (5). It should be noted that this removes the effect of the finite window on the original data since it models an infinite data sequence x,, x2 . . . based on the analysis of only a finite portion X, . . . x,. If the AR characteristics of the original data is insignificant then the transient sequence E, will approximate the original data and the distortions introduced by calculating its Fourier transform will be no worse than the original distortion. In the limit of a zeroth order AR model, the two sequences will be identical. Figure 2 shows the reconstruction of a series of doped water and gel phantom images using this modelling technique (left-hand side) and the normal FFT algorithms (right-hand side). The phase-encoded direction runs vertically. The sidelobes are very evident in the FFT reconstructed portion of the image. The AR coefficients were determined by a complex least square algorithm modified to take account of the shape of the phase-encoded signal. Images obtained with model order p = 1..5 showed no significant changes from the unmodelled images. As the model order increased (8 to 15), the sidelobes decrease and the resolution begins to improve. The sidelobes in the more complex parts of the image require a higher order before they become reduced. Very high model orders show virtually no sidelobes and excellent resolution. However such images are subject to overfitting, resulting in zeros in the AR model which produces the occasional high intensity spikes in the final image. The current images were generated with a fixed model order with no allowance made for variation in the AR fit. Only the phase-encoded direction was modelled. The final

images indicate that modelling of the original spin echo as well would remove some additional artifacts. No valid time comparisons with the normal FFT reconstruction can be made as the current series of programs were simply developed to check the feasibility of the technique. However, for low model orders and using the SIFFT (Special Interpolating FFT algorithm),” it is theoretically expected that the Transient Error method will be comparable in computational efficiency to the normal FFT reconstruction technique, especially for small model orders. Timing considerations and methods to determine the optimum model order will be discussed in a later paper. In this note, we have reported the application of the Transient Error method of AR modelling to reconstruct magnetic resonance images. This provides an estimate of the Fourier transform of the data and not the normal power spectral density function. Although considerable work is required before the full inherent advantages of modelling can be obtained, the results indicate that the technique provides a method that will reduce the artifacts and loss of resolution produced by noise and truncation errors. The selection of the correct model order required to avoid overfitting of the image and obtain the maximum noise reduction is not a trivial problem. Procedures and selection criteria are under investigation. A paper giving details of the modified complex least square algorithm for determining the AR coefficienties a, and the transient error coefficients e,, of the phase encoded signal is under preparation.” The authors wish to acknowledge the financial support of the University of Calgary, Natural Science and Engineering Research Council of Canada, the Ministry of Health of the Province of Ontario, the National Cancer Institute of Canada and the Ontario Cancer Treatment and Research Foundation.

REFERENCES Kumar, A.; Welti, D.; Ernst, R. NMR Fourier Zeug matography. Journal of Magnetic Resonance 18:69-83; 1975. Wood, M.L.; Henkelman, R.M. Truncation artifacts in magnetic resonance imaging. Journal of Magnetic Resonance in Medicine 2:6; 1985. Nichols, S.T.; Smith, M.R.; Salami, M.J.E. Application of ARMA modeling to multi component signals. 7th and system IFAC Symposium on “identification 1473-1478, July 1985, York, parameter estimation,” U.K. Harris, F.J. On the use of windows for harmonic analysis with the discrete Fourier transform. Proceedings IEEE 665 I-83; 1978.

5. Sibiski,

Staunton, practical 1984.

S.; Skilling, J.; Brereton, R.G.; Laue, E.D.; J. maximum entropy signal processing in NMR

spectroscopy.

Nature

311:446-447;

6. Barkhuijsen, H.; DeBeer, R.; Bovee, W.M.M.J.; Creyghton, J.H.N.; vanormondt, D. Application of linear prediction and singular valve decomposition (LPSVD) to determine NMR frequencies and intensities from the FID. Journal of Magnetic Resonance in Medicine 2:86-89; 1985. 7. Burg, J.D. Maximum entropy spectral analysis. Ph.D. dissertation, Dept. Geophysics, Stanford CA, 1987. 8. Barrowdale, 1.; Erickson, R.E. Algorithm for least

Application of autoregressive modelling 0 M.R. SMITHETAL.

square linear prediction and maximum entropy analysis: Part I Theory. Geophysics 45(3):420-432; 1980. 9. Kay, S.M.; Marple Jr., S.L. Spectrum analysis-a modern perspective. Proceedings IEEE 69: 1380-I 419; 1981. 10. Nichols, S.T.; Campbell, K.; Smith, M.R. Two special-

261

ized interpolating DFT algorithms. Submitted lo Proceedings IEEE; 1985. 11. Smith, M.R.; Nichols, ST.; Henkelman, R.M.; Wood, M.L. Application of parametric modelling in magnetic resonance imaging. Submitted to IEEE Transactions in Medical Imaging; 1986.