- Email: [email protected]
A novel algorithm for the reduction of undersampling artefacts in magnetic resonance images
Magnetic Resonance Imaging 22 (2004) 1279 – 1287
A novel algorithm for the reduction of undersampling artefacts in magnetic resonance images Giuseppe Placidi*, Antonello Sotgiu Centro Interdipartimentale di Risonanza Magnetica Universita` de L’Aquila, Via Vetoio 10, 67010 Coppito L’Aquila, Italy Received 4 February 2004; accepted 22 September 2004
Abstract An innovative algorithm is presented which is effective in reducing the truncation artefacts occurring in magnetic resonance images due to missing k-space samples. The algorithm works first by filling the incomplete matrix of coefficients with zeroes and then adjusting, by an iterative process, the missing coefficients by performing a reduction of the undersampling artefacts. Then, this set of coefficients is used as a basis for a superresolution algorithm that estimates the missing coefficients by modeling the data as a linear combination of increasing and decreasing exponential functions using Prony’s method. In fact, the Prony’s method consists of the interpolation of a given data set with a sum of exponential functions: the MRI signals can be well represented as a sum of exponential functions and the missing data can be extrapolated by this representation. The algorithm has been proven to perform better than either a simple algorithm, which detects and then reduces the undersampling artefacts, or an algorithm that models the measured data with approximation functions. The presented algorithm is quite simple and is applicable both to missing rows (phase-frequency acquisitions) and to radial-missing angle (acquisition from projections) undersampling. Experimental results are reported; comparisons, made between the results obtained using the presented algorithm and the alternative methods described above, clearly demonstrate the superiority of the algorithm. D 2004 Elsevier Inc. All rights reserved. Keywords: MRI; Undersampling artefacts; Threshold; Prony’s method; Exponential fitting
1. Introduction In magnetic resonance imaging (MRI), there is often the possibility that images are reconstructed using a number of k-space samples lower than that required to obtain an optimal image. This is due to several causes, the more frequent are a reduction of the acquisition time, sample motion, the effects due to the presence of metallic objects and the rapidity of the experimental conditions [1–5]. There are different ways in which undersampling can influence the resulting image because there are different ways in which it can present itself. Frequently, it occurs because several rows in the k-space grid are missing (principally in the phase encoding direction, see Fig. 1A); sometimes it occurs because several radial directions in the k-space are missing (acquisition from projections, see Fig. 1B). The nature of the truncation artefacts, depending as it does on the nature of the under* Corresponding author. Tel.: +39 862 433493; fax: +39 862 433433. E-mail address: [email protected] (G. Placidi). 0730-725X/$ – see front matter D 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.mri.2004.09.010
sampling, influences the image in different ways and, consequently, the required restoration algorithm used for treating each situation [6,7]. In this paper, both cases of missing data are treated, namely, when some rows (or columns) of the k-space samples are missing (phasefrequency imaging) along the phase encoding direction or when some projections (radial directions) are missing. In both cases, the resulting image is affected by artefacts and blurring. Several methods have been proposed to reduce undersampling effects [8–15] but the results have been frequently very poor, particularly in the presence of noise, although highly complex and computationally expensive methods have been used. Some algorithms perform better than others because they allow the extrapolation of missing coefficients by using complicated models of the measured signals and a priori knowledge of the experiments [10,15]. They often depend on many parameters, with the result that implementation of these techniques is very difficult, often unstable, and computationally expensive [15]. Moreover, when estima-
1280
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
Fig. 1. In MRI, undersampling can be caused by missing rows (indicated by dashed lines in A) in the case of phase-frequency acquisition, or by radial-missing directions (dashed lines in B) in the case of acquisition from projections.
tions have to be made utilizing inadequate data from measured samples, namely, when the number of missing k-space samples is high (long missing signal tails), these methods tend to perform poorly. Other methods, based on threshold, are very simple [16,17] but only a limited improvement in image quality is realizable because some truncation artefacts (having the form of ringing around edges) will remain. It is impossible for thresholding to eliminate truncation effects in image regions where both useful signal and artefacts are present. A novel correction algorithm is presented which, while retaining the advantages of the above-cited method, reduces truncation artefacts to a considerable degree. The idea behind this method is that the application of a particular threshold algorithm can be successfully used for a first step reduction of the truncation artefacts. In this way, it is possible to create support for a subsequent exponential approximation, based on Prony’s method [12,13,18,19], avoiding the instability associated with this method of fixing parameters, its computational requirements and its poor results when treating signals with very long missing tails. This results in a further reduction of truncation artefacts. It can also be applicable in radial undersampling. In what follows, a description of the novel algorithm will be presented, preceded by a description of both of the algorithms that it uses in cascade. The algorithm was applied to experimental MR images and some comparisons made with the previously cited techniques.
2. Details of the proposed algorithm 2.1. General notation The proposal is to reconstruct an MR image using an undersampled set of its Fourier coefficients, both in the case of linear scanning (phase-frequency acquisition) and in the case of radial scanning (acquisition from projections). To this end, let I be the theoretical image reconstructed by
using a complete set S of data coefficients. Consideration is given to the problem of having an incomplete data set S˜ =Sd U where U is the undersampling matrix filled with one in those positions where the coefficients are present ¯ represents the logical not of U). and 0 elsewhere (U Moreover, I˜ indicates the image reconstructed by using S˜ , I˜ t indicates the image I˜ after the substitution with 0 of those pixels whose amplitude falls below a given threshold t and by the subtraction of t from the other pixels. The reconstruction of the image is performed by using a twodimensional Fourier transform (2DFT) of the coefficients or, equivalently, by performing a one-dimensional Fourier transform (1DFT) along the rows of the matrix coefficients, followed by a 1DFT along the columns of the resulting transformed matrix. F r and F c indicate the 1DFT operators along the rows and the columns of a matrix given as argument, respectively. Following this notation, S˜ indicates the 1DFT of S˜ along its rows, S˜ =F r(S˜ ), and a reconstruction I˜ can be obtained by performing a 1DFT of the resulting matrix, I˜ = |F c(S˜ )|. Note that the absolute value was used in order to eliminate the phase differences due to experimental errors (magnetic field inhomogeneities, chemical shift, etc.). The inverse process is also possible, namely, to reconstruct the coefficient matrix by means of a two-dimensional inverse Fourier transform (2DIFT) of the reconstructed thresholded image I˜ t ; this process is indicated by S˜ t = 2DIFT(I˜ t ). 2.2. The threshold algorithm The threshold method is based on a postprocessing strategy which is applied to spin density images, namely, absolute value images, in order to suppress phase errors and to correct for spurious noisy negative pixels. It consists of the following steps: Algorithm Threshold. (1) (2)
Input: incomplete grid of measured, incomplete, data, S˜ ; If row-missing application,
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
(a) Fill the missing coefficients of S˜ with 0, else (angular directions are missing) (b) Interpolate the missing coefficients of S˜ using nearest neighbor; (3) Calculate I˜ = |2DFT(S˜ )|; (4) Calculate the threshold t in I˜ ; (5) Calculate the mean value e and the standard deviation r of the noise level in I˜ ; (6) If t V e+2r go to step 11; (7) Calculate I˜ t ; (8) Calculate S˜ t = 2DIFT(I˜ t ); ¯; (9) Assign S˜ = SU+S˜t U (10) Go to step 3; (11) Return S˜ .
P The threshold is assigned as t ¼ 14 4j¼1 Mj where M j = max(a j ) and a j is one of the four square submatrices of I˜ positioned at its corners (the dimensions of each submatrix are tbð1 p1ffiffi2 Þ N2 c bð1 p1ffiffi2 Þ N2 c, N N are the whole image dimensions). This value of t measures the maximum oscillation in four square regions situated at the corners of the image. The dimensions of these squares are calculated by taking into account that they have to be external to the circle inscribed in the image square, that is, external to the image region. In this way, it is possible to ensure that the value of t will depend upon undersampling artefacts or on uncorrelated noise, but not on any image features. 1 values of e and PThe P4 r are calculated as follows: e ¼ 4 4 1 j¼1 mj and r ¼ 4 j¼1 rj where m j =mean(a j ) and r j = STD(a j ). The calculations were performed in the same regions as t. The value of t is compared with e+2r in order to indicate the presence of undersampling artefacts: if the values are substantially different (t Ne+2r), then some artefacts will still occur and the process has to be repeated; otherwise, the noise on the image is independent of the artefacts and the process is stopped. Some points have to be made about the optimal threshold value and the termination criteria. The threshold value is strictly dependent on the number of missing kspace samples and on the image structure: the greater the number of missing k-space samples, the more enhanced are the truncation artefacts. In fact, having fixed the number of missing k-space samples, the truncation artefacts of a simple-shaped image are lower than for a complex-shaped image. A functional choice would be to calculate the threshold as the average of the maximum values of the peripheral image regions where the sample structures are known to be absent (note that these regions are always present on the image borders). The calculation of the parameters t, e and r is made at each repetition of the algorithm to take into account the modifications to them by artefact reduction. This last point involves the termination criteria: it would be better to stop the algorithm when the threshold value falls below e+2r. This would allow direct control of the structure of the
1281
artefacts and avoidance of the necessity to calculate other, more difficult, heuristic functions such as the image entropy [17]. These improvements make the proposed threshold method robust. Nevertheless, the application of this method to an MR image does not completely remove the effects of truncation artefacts. This fact will be clarified in the Results section. 2.3. The Prony’s approximation algorithm The presented algorithm is a modification of an existing method reported elsewhere [15]. The original method consisted of the extrapolation of missing data with Prony’s method and uses a linear combination of exponential functions to represent MRI signals. The conditions for the existence of a solution to this problem depend on the image signal-to-noise ratio (S/N), which defines the image quality, and image complexity (depending on the number of small image details that define its intrinsic resolution). The original algorithm proposed by Barone and Sebastiani [15] suffers from several problems: (a) (b) (c)
(d) (e) (f )
it depends upon the variable number of exponential functions used to approximate each signal; it depends upon the number of samples used to perform the approximation; it performs poorly when the tails of the missing coefficients are long (it has to extrapolate away from the given nodes); it is unstable and computationally expensive; it is not applicable to radial-missing angle problems (acquisition from projections); it is strongly dependent on S/N ratio data.
An improvement to this method would consist of modifying the starting image coefficients. It would be better to use the algorithm on the image modified by applying the previous threshold method instead of using the original measured coefficients matrix. Thus, the following advantages are obtained: (a) (b) (c) (d) (e)
(f )
it reduces the time spent to model each direction of coefficients; it makes the method more stable; it applies the algorithm both in the case of missing rows and in the case of missing angles; it avoids a cumbersome iterative process; it reduces the time spent choosing the best number of exponential functions to model a direction of coefficients; it fixes the number of exponential functions used to represent the model.
In the extrapolation process, some prior information is introduced: fitting with exponential functions of a particular shape to take into account the nature of the MRI data, which are recorded as echoes of oscillating exponentially
1282
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
decaying signals [15]. The part at the left of the echo is increasing and the part at its right is decreasing. The names left and right are considered with respect to the center of the echo. For this reason, the left part of each signal can be approximated with a linear combination of increasing complex exponential functions and the right part with a linear combination of decreasing complex exponential functions. It should be noted that although the requirement is to model each signal to extrapolate to its extremities, it is also important to use mainly the coefficients situated at the spectrum center where the S/ N ratio is higher in order to reduce the ill-conditioning due to noise (this represents another substantial difference from the method reported in [15]). S˜ d(l) indicates the left end and S˜ d(r) indicates the right end of the data measured along the direction d. Each direction d is composed of 2n measured elements, counted from n to n1. Given two integers n V v l b0 and 0 V v r V n1 defining the internal indices of the left and right ends of each measured direction, Prony’s method will consist of the approximation of the tails of each measured signal with a linear sum of p exponential functions: kþn p X ðlÞ ðlÞ ðlÞ S˜ d ðk Þ ¼ bj z i k ¼ n; N N ; vl 1 ð1Þ i¼1
ðrÞ S˜ d ðk Þ ¼
p X
ðr Þ
bi
ðr Þ
zi
kvr k ¼ vr ; N N ; n 1
ð2Þ
i¼1 (r) with the constraints |z(l) i |N 1 and |z i | b1, i = l,. . .,p. The constraints ensure that the left part of the signal increases and that the right part of the signal decreases when collected along the direction d. The maximum number of interpolation functions, p, was made equal to half of the number of measured points used for the model. The effective number of approximation exponential functions can be lower than p because those violating the constraints are discarded. The Prony’s method [18,19] to be used consists of the calculation of the parameters b i and z i starting from the measured coefficients S˜ d (k) (by the application of Eqs. (1) and (2) and by considering b i and z i as unknown), both for the left and the right signal tails. When the optimal values of parameters b i and z i have been estimated, the missing coefficients are calculated by applying the Eqs. (1) and (2) directly: in this case, the values of b i and z i are known and S˜ d (k) is calculated (k representing the position of a missing coefficient). It is worth noting that only those estimated missing coefficients are retained whose magnitude does not differ from the 50% of the magnitude of the same coefficients as estimated by the threshold method. The estimated coefficients are substituted for those obtained by the threshold method and the image is reconstructed. This choice is necessary to avoid bad conditioning and approximations when the method operates away from the measured data. The value 50% was chosen by experimental trial and error: a lower value leads
to a simple threshold, while a greater value requires a simple Prony’s method. The application of this approximation method will depend upon the type of undersampling problem: row-missing directions or radialmissing directions. For the first case, the algorithm is the following: Algorithm Prony_ Row. Input: S˜ (sampled data+estimated coefficients by thresholding); (2) Calculate S˜ = F r(S˜ ); (3) Assign S˜ = S˜ (4) For each column S˜ d of S˜ : (a) Estimate b i and z i for each tail of S˜ d (k) by applying Eqs. (1) and (2); (b) Calculate, using Eqs. (1) and (2), the missing (1)
coefficients S˜ d (k) from the estimated b i and z i ; For each estimated coefficient c of S˜ : If c b 0.5c or c N 1.5c, then assign c = c; (6) Calculate I˜ = |F c(S˜ )|; (7) Return I˜ .
(5)
In step 5 of the previous algorithm, c is the coefficient of S˜ corresponding to c in S˜ . S˜ is the row Fourier transform of the matrix output of Threshold. When the method is applied to radial-missing directions, some additional considerations are required. Reference should be made to Fig. 2 to illustrate the application of the method. First of all, the method is applied directly to the measured data and not to the row-transformed data: this is due to the nature of the measured data (projections). Secondly, the data are always estimated in rows and columns, although the measured data are distributed in radial directions. This implies that the sampling positions between the data along a given row are not regular (see Fig. 2 where the black thick curves cross the other curves). Thirdly, the method is
Fig. 2. The grid represents the steps followed by Prony_Radial to estimate unmeasured coefficients when radial-missing sampling occurs. In particular, black thick lines represent the measured directions; thin black lines represent the approximation of the central lines and columns; grey lines represent the approximation of the peripheral image regions (also the estimated dark grey coefficients are used for this approximation).
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
1283
to the completely sampled cross-region to cover the corners of the regular k-space image region (grey in Fig. 2). During this second application, the measured coefficients situated in the corner regions are used as control test values for the approximation (see Fig. 2 where the thick black rows cross the grey regions on the image corners). The complete description of the algorithm for this case of undersampling is the following: Algorithm Prony_ Radial.
Fig. 3. An MR image of the brain of a healthy volunteer reconstructed by using a complete set of (256256) coefficients collected by a GE Signa MRI System at 1.5 T and used as the test image.
applied twice: the first time it is applied to the central group of rows and columns (to obtain a cross of completely sampled region; thin black lines in Fig. 2), then it is applied
(1) Input: S˜ (sampled data+estimated coefficients by thresholding); (2) Assign S˜ = S˜ ; (3) For each of the central columns of S˜ : apply steps a and b of step 4 of Prony_Row; (4) For each of the central rows of S˜ : apply steps a and b of step 4 of Prony_Row; (5) Apply step 5 of Prony_Row; (6) For each corner of S˜ (I) For each peripheral column of S˜ : apply steps a and b of step 4 of Prony_Row; (II) For each peripheral row of S˜ : apply steps a and b of step 4 of Prony_ Row;
Fig. 4. The same image of Fig. 3 reconstructed by using just 50 of the 256 collected rows (from the 121st to the 170th): (A) unfiltered image; (B) image filtered by using just Threshold; (C) image filtered by using just Prony_ Row; (D) image filtered by using Restore.
1284
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
(III)
As estimated values, take the mean values between those estimated by the columns and those estimated by the rows. (IV) Apply step 5 of Prony_Row; (7) Calculate I˜ = |2DFT(S˜ )|; (8) Return I˜ .
The method used for the manipulation of the central columns or rows requires additional explanations: it is necessary to calculate the radius r of the maximum circle whose circumference is divided into arcs by the measured projections; the amplitude of these arcs must be smaller than, or equal to, the distance between two consecutive pixels. By indicating with R the radius of the circle inscribed into the image square, with N as the number of measured projections and n p as the number of points per projection (that is the same as the number of pixels per side of the image), the following inequality can be written: 2pr 2R V 2N np
ð3Þ
and the value of the radius r can be easily derived, r V (2NR/pn p).
The central columns or central rows are those whose distance from the center is lower than r. Obviously, r depends on the number of measured projections, N. 2.4. Main algorithm description The proposed algorithm can be described as follows: Algorithm Restore. (1) (2) (3) (4)
Input: matrix of the measured data, S˜ ; S˜ = Threshold(S˜ ); If row undersampling, I˜ = Prony_Row(S˜ ); else, I˜ = Prony_Radial(S˜ ); Return I˜ .
3. Results and discussion In order to demonstrate the capability of the Restore algorithm in reducing undersampling artefacts, it was applied to MRI experimental data collected at 1.5-T with a commercial GE Signa MRI System (General Electric Healthcare, Milwaukee, WI). The algorithm has been applied both in the case of row undersampling and in the
Fig. 5. Difference images obtained by subtracting from Fig. 3 the images shown in Fig. 4A–D. The difference are shown from A to D.
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
1285
Table 1 Mean square error values calculated for each of the images presented in Fig. 4 with respect to that of Fig. 3 (Row data) and for each of the images presented in Fig. 6B–D with respect to that of Fig. 6A (Radial data) Row data Radial data
Fig. 4 MSE values Fig. 6 MSE values
Zero filling
Threshold
Prony_
Restore
A 0.345 B 0.123
B 0.098 C 0.058
C 0.175 – –
D 0.035 D 0.021
case of radial-missing directions. The data referring to radial reconstruction were collected by numerical reprojection along different angular directions of the image [20] reconstructed using a complete set of spin-echo complex signals acquired in row mode (phase-frequency). In this way, a complete set of radial k-space directions was obtained. In what follows, we will show images in lowcontrast mode in order to highlight the artefacts and their reduction. In the case of row-missing directions, the image shown in Fig. 3 is the test image reconstructed by using the complete regular set of coefficients directions. The images shown in Fig. 4 are the images reconstructed by using just 50 of the 256 rows of coefficients, chosen asymmetrically between the 121st and the 170th rows. In particular. Fig. 4A
shows the reconstruction obtained by zero filling of the incomplete coefficients grid, Fig. 4B shows the result of filtering just using the Threshold algorithm, Fig. 4C shows the image as filtered just using Prony_ Row, and Fig. 4D shows the result of the application of the Restore algorithm. The images shown in Fig. 5 are the differences between the test image of Fig. 3 and the images contained in Fig. 4. The image of Fig. 4A is strongly affected by undersampling artefacts having the form of ringing around the horizontal edges of the image. These effects are reduced in Fig. 4B–D, although some unwanted new artefacts appear in the image Fig. 4C. The nature of these artefacts is due to the application of the exponential summing algorithm alone which can produce instability and poor approximations
Fig. 6. (A) The MR image obtained by reconstruction from projections using a complete set of 256 projections calculated by numerical reprojection from the image reported in Fig. 3; this image was used as the test image for comparison. The others represent the images reconstructed by using 64 regular projections: (B) without filtering, (C) applying only Threshold, and (D) applying Restore.
1286
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
Fig. 7. Difference images obtained by subtracting from Fig. 6A the images Fig. 6B–D. These differences are shown from A to C.
when operating on long tail-missing approximations. But, when Prony_Row is applied after the Threshold algorithm, as in Restore, the unwanted effects are completely absent and the image improves in quality. This is also demonstrated by the mean square error (MSE) values (Table 1, Row data), calculated for each of the images presented in Fig. 4 with respect to that shown in Fig. 3 by considering the whole rectangular image field. The comparison in Fig. 4 and in Table 1 (Row data) shows that Restore performs better than the singular application of Threshold or Prony_ Row and that Threshold gives better results than Prony_ Row. The good resemblance of the image shown in Fig. 4D to the test image clearly demonstrates the efficacy of the Restore algorithm; the undersampling artefacts are greatly reduced (although some of them at a very high frequency will always remain, as shown also by Fig. 5D). Fig. 6 shows the images reconstructed by using the radial projections numerically collected from the image shown in Fig. 3. In particular, Fig. 6A shows the image reconstructed by using the complete set of 256 projections; this is the reference image that was used for comparison purposes. Fig. 6B–D shows the images reconstructed using a set of 64 regular angular projections extracted from the complete set: without approximation (B), just using the Threshold algorithm for approximation (C), and by using the Restore algorithm (D). Note that it was not possible to reconstruct the image using only Prony_Radial because it is necessary, for its application, to have a starting regular grid that can only be obtained by using the Threshold algorithm. The images shown in Fig. 7 represent the differences between the test image Fig. 6A and each of the Fig. 6B–D. In this case as well, the images shown in Figs. 6 and 7 clearly demonstrate the improvement in image quality. Also in this case, Restore performs better than simple application of Threshold. A reduction of the radial undersampling artefacts, appearing in the form of star artefacts around the image, is also evident. Moreover, this has also been confirmed by analysis of Fig. 7 and by the values of the corresponding MSE (see Table 1, Radial data). Both in rows
and in radial undersampling, the Restore algorithm performed better than unfiltered or singularly filtered methods, producing good filtered images. 4. Conclusions An innovative filtering algorithm for MRI, capable of reducing the undersampling artefacts due to missing k-space samples, has been presented. The algorithm consists of the application of a threshold iterative method to the original kspace samples, followed by an exponential approximation method, based on Prony’s technique, on the resulting grid. The method has been proven to retain the advantages in accuracy of cumbersome exponential approximation algorithms based on Prony’s method and the simplicity of threshold algorithms. The effectiveness of the presented algorithm is, in fact, its simplicity and its dependence on a small number of approximation parameters that also make it suitable for different acquisition strategies, such as phasefrequency or projection imaging. Moreover, it always converges to a filtered image significantly improved on the initial unfiltered one. This is a very important advantage; other more complicated techniques often fail to converge (especially in conditions of low S/N ratio situations) or require long computation times. References [1] Tacke J, Adam G, Classen H, Muhler A, Prescher A, Gunther RW. Dynamic MRI of a hypovascularized liver tumor model. J Magn Reson Imaging 1997;7(July–Aug):678 – 82. [2] Wallis F, Gilbert FJ. Magnetic resonance imaging in oncology: an overview. J R Coll Surg Edinb 1999;44(Apr):117 – 25. [3] Jolesz FA, Blumenfeld SM. Interventional use of magnetic resonance imaging. Magn Reson Q 1994;10(June):85 – 96. [4] Rasche V, deBoer RW, Holz D, Proska R. Continuous radial dataacquisition for dynamic MRI. Magn Reson Med 1995;34:754 – 61. [5] Yang PC, Kerr AB, Liu AC, Liang DH, Hardy C, Meyer CH, et al. New real-time interactive cardiac magnetic resonance imaging system complements echocardiography. J Am Coll Cardiol 1998;32(Dec): 2049 – 56.
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287 [6] Wood ML, Henkelman RM. Truncation artifacts in magnetic resonance imaging. J Magn Reson Med 1985;2(December):517 – 26. [7] Lufkin RB, Pusey E, Stark DD, Brown R, Leikind B, Hanafee WN. Boundary artifacts due to truncation errors in MR imaging. AJR Am J Roentgenol 1986;147(Dec.):1283 – 7. [8] Wajer FTAW, Lethmate R, de Jong RAJ, Martinez LT, GraveronDemilly D, Fuderer M, et al. MR image reconstruction from sparse and corrupted k-space data. Proceedings of the 10th workshop on Circuits, Systems and Signal Processing ([IEEE] / [ProRISC99]). Mierlo (NL)7 STW Technology Foundation, Utrecht; 1999. p. 577 – 84. [9] Flusser J, Suk T. Degraded image analysis: an invariant approach. IEEE Trans Pattern Anal Mach Intell 1998;20(6):1 – 14. [10] Broersen PMT. Autoregressive model orders for Durbin’s MA and ARMA estimators. IEEE Trans Signal Process 2000;48(August): 2454 – 7. [11] Broersen PMT. Automatic spectral analysis with time series models. IEEE Trans Instrum Meas 2002;51(2):211 – 6. [12] Smith GK, Hawkins DM. Robust frequency estimation using elemental sets. J Comput Graph Stat 2000;9:196 – 214. [13] Osborne MR, Smyth GK. A modified Prony algorithm for exponential function fitting. SIAM J Sci Statist Comput 1995;16: 119 – 38.
1287
[14] Sebastiani G, Barone P. Truncation artifact reduction in magnetic resonance imaging by Markov random field methods. IEEE Trans Med Imaging 1995;14(September):434 – 41. [15] Barone P, Sebastiani G. A new method of magnetic resonance image reconstruction with short acquisition time and truncation artifact reduction. IEEE Trans Med Imaging 1992;11(June):250 – 9. [16] Fain SB, Block W, Barger A, Mistretta CA. Correction for artifacts in 3D angularly undersampled MR projection reconstruction. Proceedings of the International Society for Magn Reson Med. Scotland (UK)7 Glasgow; 2001. p. 759. [17] Lisboa JC, Guarini M, Irarrazaval P. A correction algorithm for undersamples images using dynamic segmentation and entropy based focus criterion. Magn Reson Imaging 2002;20:659 – 66. [18] Prony R. Essai experimental et analytique. J Ec Polytech 1795;2:24 – 76. [19] Marple SL. Digital spectral analysis with applications. Englewood Cliffs7 NJ Prentice-Hall; 1987. p. 303 – 26. [20] Placidi G, Alecci M, Colacicchi S, Sotgiu A. Fourier reconstruction as a valid alternative to filtered back projection in iterative applications: implementation of Fourier spectral spatial EPR imaging. J Magn Reson 1998;134:280 – 6.
A novel algorithm for the reduction of undersampling artefacts in magnetic resonance images Giuseppe Placidi*, Antonello Sotgiu Centro Interdipartimentale di Risonanza Magnetica Universita` de L’Aquila, Via Vetoio 10, 67010 Coppito L’Aquila, Italy Received 4 February 2004; accepted 22 September 2004
Abstract An innovative algorithm is presented which is effective in reducing the truncation artefacts occurring in magnetic resonance images due to missing k-space samples. The algorithm works first by filling the incomplete matrix of coefficients with zeroes and then adjusting, by an iterative process, the missing coefficients by performing a reduction of the undersampling artefacts. Then, this set of coefficients is used as a basis for a superresolution algorithm that estimates the missing coefficients by modeling the data as a linear combination of increasing and decreasing exponential functions using Prony’s method. In fact, the Prony’s method consists of the interpolation of a given data set with a sum of exponential functions: the MRI signals can be well represented as a sum of exponential functions and the missing data can be extrapolated by this representation. The algorithm has been proven to perform better than either a simple algorithm, which detects and then reduces the undersampling artefacts, or an algorithm that models the measured data with approximation functions. The presented algorithm is quite simple and is applicable both to missing rows (phase-frequency acquisitions) and to radial-missing angle (acquisition from projections) undersampling. Experimental results are reported; comparisons, made between the results obtained using the presented algorithm and the alternative methods described above, clearly demonstrate the superiority of the algorithm. D 2004 Elsevier Inc. All rights reserved. Keywords: MRI; Undersampling artefacts; Threshold; Prony’s method; Exponential fitting
1. Introduction In magnetic resonance imaging (MRI), there is often the possibility that images are reconstructed using a number of k-space samples lower than that required to obtain an optimal image. This is due to several causes, the more frequent are a reduction of the acquisition time, sample motion, the effects due to the presence of metallic objects and the rapidity of the experimental conditions [1–5]. There are different ways in which undersampling can influence the resulting image because there are different ways in which it can present itself. Frequently, it occurs because several rows in the k-space grid are missing (principally in the phase encoding direction, see Fig. 1A); sometimes it occurs because several radial directions in the k-space are missing (acquisition from projections, see Fig. 1B). The nature of the truncation artefacts, depending as it does on the nature of the under* Corresponding author. Tel.: +39 862 433493; fax: +39 862 433433. E-mail address: [email protected] (G. Placidi). 0730-725X/$ – see front matter D 2004 Elsevier Inc. All rights reserved. doi:10.1016/j.mri.2004.09.010
sampling, influences the image in different ways and, consequently, the required restoration algorithm used for treating each situation [6,7]. In this paper, both cases of missing data are treated, namely, when some rows (or columns) of the k-space samples are missing (phasefrequency imaging) along the phase encoding direction or when some projections (radial directions) are missing. In both cases, the resulting image is affected by artefacts and blurring. Several methods have been proposed to reduce undersampling effects [8–15] but the results have been frequently very poor, particularly in the presence of noise, although highly complex and computationally expensive methods have been used. Some algorithms perform better than others because they allow the extrapolation of missing coefficients by using complicated models of the measured signals and a priori knowledge of the experiments [10,15]. They often depend on many parameters, with the result that implementation of these techniques is very difficult, often unstable, and computationally expensive [15]. Moreover, when estima-
1280
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
Fig. 1. In MRI, undersampling can be caused by missing rows (indicated by dashed lines in A) in the case of phase-frequency acquisition, or by radial-missing directions (dashed lines in B) in the case of acquisition from projections.
tions have to be made utilizing inadequate data from measured samples, namely, when the number of missing k-space samples is high (long missing signal tails), these methods tend to perform poorly. Other methods, based on threshold, are very simple [16,17] but only a limited improvement in image quality is realizable because some truncation artefacts (having the form of ringing around edges) will remain. It is impossible for thresholding to eliminate truncation effects in image regions where both useful signal and artefacts are present. A novel correction algorithm is presented which, while retaining the advantages of the above-cited method, reduces truncation artefacts to a considerable degree. The idea behind this method is that the application of a particular threshold algorithm can be successfully used for a first step reduction of the truncation artefacts. In this way, it is possible to create support for a subsequent exponential approximation, based on Prony’s method [12,13,18,19], avoiding the instability associated with this method of fixing parameters, its computational requirements and its poor results when treating signals with very long missing tails. This results in a further reduction of truncation artefacts. It can also be applicable in radial undersampling. In what follows, a description of the novel algorithm will be presented, preceded by a description of both of the algorithms that it uses in cascade. The algorithm was applied to experimental MR images and some comparisons made with the previously cited techniques.
2. Details of the proposed algorithm 2.1. General notation The proposal is to reconstruct an MR image using an undersampled set of its Fourier coefficients, both in the case of linear scanning (phase-frequency acquisition) and in the case of radial scanning (acquisition from projections). To this end, let I be the theoretical image reconstructed by
using a complete set S of data coefficients. Consideration is given to the problem of having an incomplete data set S˜ =Sd U where U is the undersampling matrix filled with one in those positions where the coefficients are present ¯ represents the logical not of U). and 0 elsewhere (U Moreover, I˜ indicates the image reconstructed by using S˜ , I˜ t indicates the image I˜ after the substitution with 0 of those pixels whose amplitude falls below a given threshold t and by the subtraction of t from the other pixels. The reconstruction of the image is performed by using a twodimensional Fourier transform (2DFT) of the coefficients or, equivalently, by performing a one-dimensional Fourier transform (1DFT) along the rows of the matrix coefficients, followed by a 1DFT along the columns of the resulting transformed matrix. F r and F c indicate the 1DFT operators along the rows and the columns of a matrix given as argument, respectively. Following this notation, S˜ indicates the 1DFT of S˜ along its rows, S˜ =F r(S˜ ), and a reconstruction I˜ can be obtained by performing a 1DFT of the resulting matrix, I˜ = |F c(S˜ )|. Note that the absolute value was used in order to eliminate the phase differences due to experimental errors (magnetic field inhomogeneities, chemical shift, etc.). The inverse process is also possible, namely, to reconstruct the coefficient matrix by means of a two-dimensional inverse Fourier transform (2DIFT) of the reconstructed thresholded image I˜ t ; this process is indicated by S˜ t = 2DIFT(I˜ t ). 2.2. The threshold algorithm The threshold method is based on a postprocessing strategy which is applied to spin density images, namely, absolute value images, in order to suppress phase errors and to correct for spurious noisy negative pixels. It consists of the following steps: Algorithm Threshold. (1) (2)
Input: incomplete grid of measured, incomplete, data, S˜ ; If row-missing application,
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
(a) Fill the missing coefficients of S˜ with 0, else (angular directions are missing) (b) Interpolate the missing coefficients of S˜ using nearest neighbor; (3) Calculate I˜ = |2DFT(S˜ )|; (4) Calculate the threshold t in I˜ ; (5) Calculate the mean value e and the standard deviation r of the noise level in I˜ ; (6) If t V e+2r go to step 11; (7) Calculate I˜ t ; (8) Calculate S˜ t = 2DIFT(I˜ t ); ¯; (9) Assign S˜ = SU+S˜t U (10) Go to step 3; (11) Return S˜ .
P The threshold is assigned as t ¼ 14 4j¼1 Mj where M j = max(a j ) and a j is one of the four square submatrices of I˜ positioned at its corners (the dimensions of each submatrix are tbð1 p1ffiffi2 Þ N2 c bð1 p1ffiffi2 Þ N2 c, N N are the whole image dimensions). This value of t measures the maximum oscillation in four square regions situated at the corners of the image. The dimensions of these squares are calculated by taking into account that they have to be external to the circle inscribed in the image square, that is, external to the image region. In this way, it is possible to ensure that the value of t will depend upon undersampling artefacts or on uncorrelated noise, but not on any image features. 1 values of e and PThe P4 r are calculated as follows: e ¼ 4 4 1 j¼1 mj and r ¼ 4 j¼1 rj where m j =mean(a j ) and r j = STD(a j ). The calculations were performed in the same regions as t. The value of t is compared with e+2r in order to indicate the presence of undersampling artefacts: if the values are substantially different (t Ne+2r), then some artefacts will still occur and the process has to be repeated; otherwise, the noise on the image is independent of the artefacts and the process is stopped. Some points have to be made about the optimal threshold value and the termination criteria. The threshold value is strictly dependent on the number of missing kspace samples and on the image structure: the greater the number of missing k-space samples, the more enhanced are the truncation artefacts. In fact, having fixed the number of missing k-space samples, the truncation artefacts of a simple-shaped image are lower than for a complex-shaped image. A functional choice would be to calculate the threshold as the average of the maximum values of the peripheral image regions where the sample structures are known to be absent (note that these regions are always present on the image borders). The calculation of the parameters t, e and r is made at each repetition of the algorithm to take into account the modifications to them by artefact reduction. This last point involves the termination criteria: it would be better to stop the algorithm when the threshold value falls below e+2r. This would allow direct control of the structure of the
1281
artefacts and avoidance of the necessity to calculate other, more difficult, heuristic functions such as the image entropy [17]. These improvements make the proposed threshold method robust. Nevertheless, the application of this method to an MR image does not completely remove the effects of truncation artefacts. This fact will be clarified in the Results section. 2.3. The Prony’s approximation algorithm The presented algorithm is a modification of an existing method reported elsewhere [15]. The original method consisted of the extrapolation of missing data with Prony’s method and uses a linear combination of exponential functions to represent MRI signals. The conditions for the existence of a solution to this problem depend on the image signal-to-noise ratio (S/N), which defines the image quality, and image complexity (depending on the number of small image details that define its intrinsic resolution). The original algorithm proposed by Barone and Sebastiani [15] suffers from several problems: (a) (b) (c)
(d) (e) (f )
it depends upon the variable number of exponential functions used to approximate each signal; it depends upon the number of samples used to perform the approximation; it performs poorly when the tails of the missing coefficients are long (it has to extrapolate away from the given nodes); it is unstable and computationally expensive; it is not applicable to radial-missing angle problems (acquisition from projections); it is strongly dependent on S/N ratio data.
An improvement to this method would consist of modifying the starting image coefficients. It would be better to use the algorithm on the image modified by applying the previous threshold method instead of using the original measured coefficients matrix. Thus, the following advantages are obtained: (a) (b) (c) (d) (e)
(f )
it reduces the time spent to model each direction of coefficients; it makes the method more stable; it applies the algorithm both in the case of missing rows and in the case of missing angles; it avoids a cumbersome iterative process; it reduces the time spent choosing the best number of exponential functions to model a direction of coefficients; it fixes the number of exponential functions used to represent the model.
In the extrapolation process, some prior information is introduced: fitting with exponential functions of a particular shape to take into account the nature of the MRI data, which are recorded as echoes of oscillating exponentially
1282
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
decaying signals [15]. The part at the left of the echo is increasing and the part at its right is decreasing. The names left and right are considered with respect to the center of the echo. For this reason, the left part of each signal can be approximated with a linear combination of increasing complex exponential functions and the right part with a linear combination of decreasing complex exponential functions. It should be noted that although the requirement is to model each signal to extrapolate to its extremities, it is also important to use mainly the coefficients situated at the spectrum center where the S/ N ratio is higher in order to reduce the ill-conditioning due to noise (this represents another substantial difference from the method reported in [15]). S˜ d(l) indicates the left end and S˜ d(r) indicates the right end of the data measured along the direction d. Each direction d is composed of 2n measured elements, counted from n to n1. Given two integers n V v l b0 and 0 V v r V n1 defining the internal indices of the left and right ends of each measured direction, Prony’s method will consist of the approximation of the tails of each measured signal with a linear sum of p exponential functions: kþn p X ðlÞ ðlÞ ðlÞ S˜ d ðk Þ ¼ bj z i k ¼ n; N N ; vl 1 ð1Þ i¼1
ðrÞ S˜ d ðk Þ ¼
p X
ðr Þ
bi
ðr Þ
zi
kvr k ¼ vr ; N N ; n 1
ð2Þ
i¼1 (r) with the constraints |z(l) i |N 1 and |z i | b1, i = l,. . .,p. The constraints ensure that the left part of the signal increases and that the right part of the signal decreases when collected along the direction d. The maximum number of interpolation functions, p, was made equal to half of the number of measured points used for the model. The effective number of approximation exponential functions can be lower than p because those violating the constraints are discarded. The Prony’s method [18,19] to be used consists of the calculation of the parameters b i and z i starting from the measured coefficients S˜ d (k) (by the application of Eqs. (1) and (2) and by considering b i and z i as unknown), both for the left and the right signal tails. When the optimal values of parameters b i and z i have been estimated, the missing coefficients are calculated by applying the Eqs. (1) and (2) directly: in this case, the values of b i and z i are known and S˜ d (k) is calculated (k representing the position of a missing coefficient). It is worth noting that only those estimated missing coefficients are retained whose magnitude does not differ from the 50% of the magnitude of the same coefficients as estimated by the threshold method. The estimated coefficients are substituted for those obtained by the threshold method and the image is reconstructed. This choice is necessary to avoid bad conditioning and approximations when the method operates away from the measured data. The value 50% was chosen by experimental trial and error: a lower value leads
to a simple threshold, while a greater value requires a simple Prony’s method. The application of this approximation method will depend upon the type of undersampling problem: row-missing directions or radialmissing directions. For the first case, the algorithm is the following: Algorithm Prony_ Row. Input: S˜ (sampled data+estimated coefficients by thresholding); (2) Calculate S˜ = F r(S˜ ); (3) Assign S˜ = S˜ (4) For each column S˜ d of S˜ : (a) Estimate b i and z i for each tail of S˜ d (k) by applying Eqs. (1) and (2); (b) Calculate, using Eqs. (1) and (2), the missing (1)
coefficients S˜ d (k) from the estimated b i and z i ; For each estimated coefficient c of S˜ : If c b 0.5c or c N 1.5c, then assign c = c; (6) Calculate I˜ = |F c(S˜ )|; (7) Return I˜ .
(5)
In step 5 of the previous algorithm, c is the coefficient of S˜ corresponding to c in S˜ . S˜ is the row Fourier transform of the matrix output of Threshold. When the method is applied to radial-missing directions, some additional considerations are required. Reference should be made to Fig. 2 to illustrate the application of the method. First of all, the method is applied directly to the measured data and not to the row-transformed data: this is due to the nature of the measured data (projections). Secondly, the data are always estimated in rows and columns, although the measured data are distributed in radial directions. This implies that the sampling positions between the data along a given row are not regular (see Fig. 2 where the black thick curves cross the other curves). Thirdly, the method is
Fig. 2. The grid represents the steps followed by Prony_Radial to estimate unmeasured coefficients when radial-missing sampling occurs. In particular, black thick lines represent the measured directions; thin black lines represent the approximation of the central lines and columns; grey lines represent the approximation of the peripheral image regions (also the estimated dark grey coefficients are used for this approximation).
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
1283
to the completely sampled cross-region to cover the corners of the regular k-space image region (grey in Fig. 2). During this second application, the measured coefficients situated in the corner regions are used as control test values for the approximation (see Fig. 2 where the thick black rows cross the grey regions on the image corners). The complete description of the algorithm for this case of undersampling is the following: Algorithm Prony_ Radial.
Fig. 3. An MR image of the brain of a healthy volunteer reconstructed by using a complete set of (256256) coefficients collected by a GE Signa MRI System at 1.5 T and used as the test image.
applied twice: the first time it is applied to the central group of rows and columns (to obtain a cross of completely sampled region; thin black lines in Fig. 2), then it is applied
(1) Input: S˜ (sampled data+estimated coefficients by thresholding); (2) Assign S˜ = S˜ ; (3) For each of the central columns of S˜ : apply steps a and b of step 4 of Prony_Row; (4) For each of the central rows of S˜ : apply steps a and b of step 4 of Prony_Row; (5) Apply step 5 of Prony_Row; (6) For each corner of S˜ (I) For each peripheral column of S˜ : apply steps a and b of step 4 of Prony_Row; (II) For each peripheral row of S˜ : apply steps a and b of step 4 of Prony_ Row;
Fig. 4. The same image of Fig. 3 reconstructed by using just 50 of the 256 collected rows (from the 121st to the 170th): (A) unfiltered image; (B) image filtered by using just Threshold; (C) image filtered by using just Prony_ Row; (D) image filtered by using Restore.
1284
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
(III)
As estimated values, take the mean values between those estimated by the columns and those estimated by the rows. (IV) Apply step 5 of Prony_Row; (7) Calculate I˜ = |2DFT(S˜ )|; (8) Return I˜ .
The method used for the manipulation of the central columns or rows requires additional explanations: it is necessary to calculate the radius r of the maximum circle whose circumference is divided into arcs by the measured projections; the amplitude of these arcs must be smaller than, or equal to, the distance between two consecutive pixels. By indicating with R the radius of the circle inscribed into the image square, with N as the number of measured projections and n p as the number of points per projection (that is the same as the number of pixels per side of the image), the following inequality can be written: 2pr 2R V 2N np
ð3Þ
and the value of the radius r can be easily derived, r V (2NR/pn p).
The central columns or central rows are those whose distance from the center is lower than r. Obviously, r depends on the number of measured projections, N. 2.4. Main algorithm description The proposed algorithm can be described as follows: Algorithm Restore. (1) (2) (3) (4)
Input: matrix of the measured data, S˜ ; S˜ = Threshold(S˜ ); If row undersampling, I˜ = Prony_Row(S˜ ); else, I˜ = Prony_Radial(S˜ ); Return I˜ .
3. Results and discussion In order to demonstrate the capability of the Restore algorithm in reducing undersampling artefacts, it was applied to MRI experimental data collected at 1.5-T with a commercial GE Signa MRI System (General Electric Healthcare, Milwaukee, WI). The algorithm has been applied both in the case of row undersampling and in the
Fig. 5. Difference images obtained by subtracting from Fig. 3 the images shown in Fig. 4A–D. The difference are shown from A to D.
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
1285
Table 1 Mean square error values calculated for each of the images presented in Fig. 4 with respect to that of Fig. 3 (Row data) and for each of the images presented in Fig. 6B–D with respect to that of Fig. 6A (Radial data) Row data Radial data
Fig. 4 MSE values Fig. 6 MSE values
Zero filling
Threshold
Prony_
Restore
A 0.345 B 0.123
B 0.098 C 0.058
C 0.175 – –
D 0.035 D 0.021
case of radial-missing directions. The data referring to radial reconstruction were collected by numerical reprojection along different angular directions of the image [20] reconstructed using a complete set of spin-echo complex signals acquired in row mode (phase-frequency). In this way, a complete set of radial k-space directions was obtained. In what follows, we will show images in lowcontrast mode in order to highlight the artefacts and their reduction. In the case of row-missing directions, the image shown in Fig. 3 is the test image reconstructed by using the complete regular set of coefficients directions. The images shown in Fig. 4 are the images reconstructed by using just 50 of the 256 rows of coefficients, chosen asymmetrically between the 121st and the 170th rows. In particular. Fig. 4A
shows the reconstruction obtained by zero filling of the incomplete coefficients grid, Fig. 4B shows the result of filtering just using the Threshold algorithm, Fig. 4C shows the image as filtered just using Prony_ Row, and Fig. 4D shows the result of the application of the Restore algorithm. The images shown in Fig. 5 are the differences between the test image of Fig. 3 and the images contained in Fig. 4. The image of Fig. 4A is strongly affected by undersampling artefacts having the form of ringing around the horizontal edges of the image. These effects are reduced in Fig. 4B–D, although some unwanted new artefacts appear in the image Fig. 4C. The nature of these artefacts is due to the application of the exponential summing algorithm alone which can produce instability and poor approximations
Fig. 6. (A) The MR image obtained by reconstruction from projections using a complete set of 256 projections calculated by numerical reprojection from the image reported in Fig. 3; this image was used as the test image for comparison. The others represent the images reconstructed by using 64 regular projections: (B) without filtering, (C) applying only Threshold, and (D) applying Restore.
1286
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287
Fig. 7. Difference images obtained by subtracting from Fig. 6A the images Fig. 6B–D. These differences are shown from A to C.
when operating on long tail-missing approximations. But, when Prony_Row is applied after the Threshold algorithm, as in Restore, the unwanted effects are completely absent and the image improves in quality. This is also demonstrated by the mean square error (MSE) values (Table 1, Row data), calculated for each of the images presented in Fig. 4 with respect to that shown in Fig. 3 by considering the whole rectangular image field. The comparison in Fig. 4 and in Table 1 (Row data) shows that Restore performs better than the singular application of Threshold or Prony_ Row and that Threshold gives better results than Prony_ Row. The good resemblance of the image shown in Fig. 4D to the test image clearly demonstrates the efficacy of the Restore algorithm; the undersampling artefacts are greatly reduced (although some of them at a very high frequency will always remain, as shown also by Fig. 5D). Fig. 6 shows the images reconstructed by using the radial projections numerically collected from the image shown in Fig. 3. In particular, Fig. 6A shows the image reconstructed by using the complete set of 256 projections; this is the reference image that was used for comparison purposes. Fig. 6B–D shows the images reconstructed using a set of 64 regular angular projections extracted from the complete set: without approximation (B), just using the Threshold algorithm for approximation (C), and by using the Restore algorithm (D). Note that it was not possible to reconstruct the image using only Prony_Radial because it is necessary, for its application, to have a starting regular grid that can only be obtained by using the Threshold algorithm. The images shown in Fig. 7 represent the differences between the test image Fig. 6A and each of the Fig. 6B–D. In this case as well, the images shown in Figs. 6 and 7 clearly demonstrate the improvement in image quality. Also in this case, Restore performs better than simple application of Threshold. A reduction of the radial undersampling artefacts, appearing in the form of star artefacts around the image, is also evident. Moreover, this has also been confirmed by analysis of Fig. 7 and by the values of the corresponding MSE (see Table 1, Radial data). Both in rows
and in radial undersampling, the Restore algorithm performed better than unfiltered or singularly filtered methods, producing good filtered images. 4. Conclusions An innovative filtering algorithm for MRI, capable of reducing the undersampling artefacts due to missing k-space samples, has been presented. The algorithm consists of the application of a threshold iterative method to the original kspace samples, followed by an exponential approximation method, based on Prony’s technique, on the resulting grid. The method has been proven to retain the advantages in accuracy of cumbersome exponential approximation algorithms based on Prony’s method and the simplicity of threshold algorithms. The effectiveness of the presented algorithm is, in fact, its simplicity and its dependence on a small number of approximation parameters that also make it suitable for different acquisition strategies, such as phasefrequency or projection imaging. Moreover, it always converges to a filtered image significantly improved on the initial unfiltered one. This is a very important advantage; other more complicated techniques often fail to converge (especially in conditions of low S/N ratio situations) or require long computation times. References [1] Tacke J, Adam G, Classen H, Muhler A, Prescher A, Gunther RW. Dynamic MRI of a hypovascularized liver tumor model. J Magn Reson Imaging 1997;7(July–Aug):678 – 82. [2] Wallis F, Gilbert FJ. Magnetic resonance imaging in oncology: an overview. J R Coll Surg Edinb 1999;44(Apr):117 – 25. [3] Jolesz FA, Blumenfeld SM. Interventional use of magnetic resonance imaging. Magn Reson Q 1994;10(June):85 – 96. [4] Rasche V, deBoer RW, Holz D, Proska R. Continuous radial dataacquisition for dynamic MRI. Magn Reson Med 1995;34:754 – 61. [5] Yang PC, Kerr AB, Liu AC, Liang DH, Hardy C, Meyer CH, et al. New real-time interactive cardiac magnetic resonance imaging system complements echocardiography. J Am Coll Cardiol 1998;32(Dec): 2049 – 56.
G. Placidi, A. Sotgiu / Magnetic Resonance Imaging 22 (2004) 1279 –1287 [6] Wood ML, Henkelman RM. Truncation artifacts in magnetic resonance imaging. J Magn Reson Med 1985;2(December):517 – 26. [7] Lufkin RB, Pusey E, Stark DD, Brown R, Leikind B, Hanafee WN. Boundary artifacts due to truncation errors in MR imaging. AJR Am J Roentgenol 1986;147(Dec.):1283 – 7. [8] Wajer FTAW, Lethmate R, de Jong RAJ, Martinez LT, GraveronDemilly D, Fuderer M, et al. MR image reconstruction from sparse and corrupted k-space data. Proceedings of the 10th workshop on Circuits, Systems and Signal Processing ([IEEE] / [ProRISC99]). Mierlo (NL)7 STW Technology Foundation, Utrecht; 1999. p. 577 – 84. [9] Flusser J, Suk T. Degraded image analysis: an invariant approach. IEEE Trans Pattern Anal Mach Intell 1998;20(6):1 – 14. [10] Broersen PMT. Autoregressive model orders for Durbin’s MA and ARMA estimators. IEEE Trans Signal Process 2000;48(August): 2454 – 7. [11] Broersen PMT. Automatic spectral analysis with time series models. IEEE Trans Instrum Meas 2002;51(2):211 – 6. [12] Smith GK, Hawkins DM. Robust frequency estimation using elemental sets. J Comput Graph Stat 2000;9:196 – 214. [13] Osborne MR, Smyth GK. A modified Prony algorithm for exponential function fitting. SIAM J Sci Statist Comput 1995;16: 119 – 38.
1287
[14] Sebastiani G, Barone P. Truncation artifact reduction in magnetic resonance imaging by Markov random field methods. IEEE Trans Med Imaging 1995;14(September):434 – 41. [15] Barone P, Sebastiani G. A new method of magnetic resonance image reconstruction with short acquisition time and truncation artifact reduction. IEEE Trans Med Imaging 1992;11(June):250 – 9. [16] Fain SB, Block W, Barger A, Mistretta CA. Correction for artifacts in 3D angularly undersampled MR projection reconstruction. Proceedings of the International Society for Magn Reson Med. Scotland (UK)7 Glasgow; 2001. p. 759. [17] Lisboa JC, Guarini M, Irarrazaval P. A correction algorithm for undersamples images using dynamic segmentation and entropy based focus criterion. Magn Reson Imaging 2002;20:659 – 66. [18] Prony R. Essai experimental et analytique. J Ec Polytech 1795;2:24 – 76. [19] Marple SL. Digital spectral analysis with applications. Englewood Cliffs7 NJ Prentice-Hall; 1987. p. 303 – 26. [20] Placidi G, Alecci M, Colacicchi S, Sotgiu A. Fourier reconstruction as a valid alternative to filtered back projection in iterative applications: implementation of Fourier spectral spatial EPR imaging. J Magn Reson 1998;134:280 – 6.