JOURNAL OF MAGNETIC RESONANCE, ARTICLE NO.
Series A 121, 60–64 (1996)
0136
Optimization of Electron Paramagnetic Resonance Image Reconstruction Using Filtered Back-Projection Followed by Two-Dimensional Deconvolution G. PLACIDI,* M. ALECCI, G. GUALTIERI,
AND
A. SOTGIU
INFM, c/o Dipartimento di Scienze e Tecnologie Biomediche, Universita’ dell’Aquila, Coppito, 67100 L’Aquila, Italy Received December 26, 1995; revised March 18, 1996
In continuous-wave electron-paramagnetic-resonance imaging, the spectra are recorded by sweeping the magnetic field H0 in the presence of a stationary field gradient G (1, 2). Each spectrum bu(r) represents the convolution of the projection pu(r) of the sample spin density along the gradient direction, at angle u, with the lineshape h(r) of the observed paramagnetic species, bu(r) Å
*
`
h(s)pu(s 0 r)ds Å pu(r) # h(r),
by the following two steps: one-dimensional deconvolution of the acquired spectra, followed by image reconstruction from projections (DR sequence). Theoretically, in the absence of noise, the order of application of the two steps, deconvolution–reconstruction, is irrelevant. In this Note, we present a new EPR method of image reconstruction: direct image reconstruction from the acquired spectra, followed by two-dimensional deconvolution (RD sequence). We show experimentally that the RD sequence generates better images in terms of spatial resolution, noise amplification, and maintenance of proportions when compared to the images obtained using the DR sequence. The method has been developed for signals, including one component of the hyperfine triplet only, in view of in vivo applications. For this reason, the experiments have been made avoiding the crossing of hyperfine peaks by limiting the strength of the field gradient. One-dimensional deconvolution techniques offer an attractive possibility for increasing spatial resolution and several methods have been developed (6–8). Fourier deconvolution is a simple and easy method to implement. If Bu(k), H(k), and Pu(k) are the 1D Fourier transforms of bu(r), h(r), and pu(r), respectively, the projection of the spin distribution is obtained evaluating the 1D deconvolution equation,
[1]
0`
where the symbol # indicates convolution. This dependence of the measured spectra or projections on instrumental parameters or on the lineshape of the measured paramagnetic species is common to other experimental situations and imaging techniques (3, 4). However, a characteristic feature of EPR is that the lineshape function h(r) is broad and a reconstruction using bu(r) would be very poor. This problem is more pronounced in in vivo applications of EPR where the probe concentration is excessively low, and high values of field gradients cannot be used (5). Under these conditions, for a given sample dimension, the maximum gradient strength is limited and the hyperfine peaks of the nitroxide probe do not overlap. Deconvolution is often used to reduce the effect of the linewidth of the measured projections (6–8). The use of deconvolution implies that the linewidth function h(r) is space-invariant along the sample. In in vivo EPR imaging, this seems very often to be true because the free radicals used are water soluble to avoid loss of sensitivity. On the other hand, in the case of an EPR signal composed of multiline spectra with variable linewidths along the sample, other techniques like spectral–spatial imaging (9, 10) are more appropriate. In EPR, it is customary to perform image reconstruction
pu(r) Å FT 01 1D {[Bu(k)/H(k)]G(k)},
where FT 01 1D is the one-dimensional reverse Fourier transform (6). The function G(k) is a low-pass filter, generally a Hamming function (11), that is characterized by two parameters: shape a, representing the smoothness of the filter, which as it goes from 0 to 1 becomes sharper, and bandwidth b, representing the maximum number of nonzero frequencies. It is used in the frequency domain to eliminate the highfrequency components where the noise dominates the signal and is much more amplified by the ratio in Eq. [2]. To perform image reconstruction, the projections pu(r), collected along a set of field-gradient orientations in polar
* To whom correspondence should be addressed.
1064-1858/96 $18.00 Copyright q 1996 by Academic Press, Inc. All rights of reproduction in any form reserved.
6j0e$$0893
06-05-96 02:24:12
[2]
60
magas
AP: Mag Res
61
NOTES
coordinates, are used to obtain the sample spin density f ( x, y) by filtered back projection (FBP) (12, 13), f ( x, y)
* p * (r)du Å * F * P (k)ÉkÉ e p
Å
u
0
`
p
02p ikr
u
0`
0
G
dk du,
[3]
where r is taken on the x–y plane such that r Å x cos u / y sin u, and p * u (r) is the projection pu(r) filtered according to the expression inside the square brackets. Deconvolution and reconstruction are both inverse problems which are performed in sequence. Neither of them is able to reconstruct exactly the theoretical function, f ( x, y), but each produces errors which are amplified by the application of the other. From this point of view, the deconvolution algorithm must implement a division (see Eq. [2]): the errors that it introduces could be very large because of the possibility of divisions by zero, or by values near zero. The reason for this is related to the effect on the reconstruction of the spurious oscillations and the distortions produced by deconvolution. The presence of noise in the experimental data enhances this effect which is more evident at low signalto-noise ratios (S/N). On the other hand, FBP, like any other reconstruction algorithm, is performed as a summation (see Eq. [3]) and it produces smaller errors than deconvolution. In this case, the errors are principally due to the limited number of measured projections. For this reason, we have developed the RD sequence, in which deconvolution is performed as the last step, because it is less subject to amplification of the errors. Moreover, the RD sequence performs the deconvolution over the whole image, while DR performs it over each single projection. This fact implies that the distortions, the amplification of spurious oscillations, and any other types of error are generated over the whole image in RD, maintaining all the relative properties of the image. In DR it is more difficult, impossible in some cases, to assemble the data correctly for the reconstruction algorithm. To obtain the sample spin density f ( x, y) by the RD method we first use the FBP equation, fb (x, y)
* b * (r)du Å * F * B (k)ÉkÉ e p
Å
u
0
p
`
u
0
0`
02p ikr
G
dk du,
[4]
where r is r Å x cos u / y sin u, b* u (r) is the measured
6j0e$$0893
06-05-96 02:24:12
magas
FIG. 1. Spectra obtained at X band of a phantom composed of four capillaries containing a 2 mM potassium nitrosodisulfonate water solution: (A) the zero-gradient spectrum; (B) the measured 1D projection at S/N Å 400; and (C) the calculated projection at S/N Å 40, obtained by adding to (B) white Gaussian noise with a zero mean. The positions and the dimensions of the four capillaries are given in Fig. 2.
projection bu(r) filtered according to the expression included in the square brackets, and fb (x, y) is the 2D image reconstructed using the measured projections. It is worth noting that the image fb (x, y) differs from the sample spin density f ( x, y) (obtained by means of Eq. [3]) because it contains the broadening effect due to the lineshape of the paramagnetic probe. Let Fb (kx , ky ) be the 2D Fourier transform of the image fb (x, y) obtained from the measured projections by using Eq. [4]. Let H(kx , ky ) be the 2D Fourier transform of the zero-gradient 2D lineshape function h(x, y) of the probe. The lineshape function h(x, y) is calculated by rotating the integral of the zero-gradient EPR spectrum h(r) acquired in derivative form (see Fig. 1A). In the present work, only the low-field EPR transition of the paramagnetic probe was acquired. Because the function h(x, y) is calculated from the experimental EPR spectra, it may contain additional distortion introduced by the spectrometer. The required image f ( x, y) is given by the 2D Fourier deconvolution f ( x, y) Å FT 01 2D {[Fb (kx , ky )/H(kx , ky )]G(kx , ky )}, [5]
AP: Mag Res
62
NOTES
where FT 01 2D is the 2D reverse Fourier transform ( 6), and G(kx , ky ) is the 2D version of the Hamming filter, obtained by rotating the 1D low-pass Hamming filter (11) around its maximum. Under these conditions, the parameters a and b have the same meaning both for the 1D and the 2D filter. For this reason, in the following, the parameters a and b of the filter, for both the 1D and the 2D case, will be the same for the two sequences when applied to the same data. It is easy to demonstrate that the sequences DR and RD are theoretically equivalent. In fact, if F(kx , ky ) is the Fourier transform of the theoretical image f ( x, y) and k Å kx cos u / ky sin u, from the projection-slice theorem (12), we have F(kx , ky ) Å Pu(k) Fb (kx , ky ) Å Bu(k).
[6]
ferent values of S/N would have been that of lowering the sample concentration or the modulation amplitude. However, the first requires the removal of the sample from the cavity, changing the experimental conditions, and the second produces only a moderate decrease of the S/N, making the comparison difficult to make over a wide span of S/N. The projections were used to reconstruct 2D images using the DR and RD sequences and an example is shown in Fig. 2. A quantitative comparison between the images has been made evaluating three parameters. The first, q, accounting for resolution, measures the degree of overlapping between the lines representing contiguous objects. It is best calculated using a simplified model of two samples along the direction of the field gradient at a given distance between them. The two lines will be coincidental or more or less overlapping, depending on the linewidth and gradient intensity. In this case, q can be defined as
Then, using Eq. [2], Eq. [5], and Eqs. [6], we can write Fb (kx , ky ) Bu(k) Å Å Pu(k). F(kx , ky ) Å H(kx , ky ) H(k)
[7]
This equation shows the theoretical equivalence between the DR and RD sequences. In fact, the relations that it contains are all equalities. Then, starting from Bu(k), it is possible to obtain f ( x, y), in one of the following two sequences: Bu(k) cEq.2 Pu(k) cEq.3 f ( x, y)
(DR)
Bu(k) cEq.4 Fb (kx , ky ) cEq.5 f ( x, y).
(RD)
The different behavior of the DR and RD sequences, when applied to data at different S/N, was tested with experimental EPR data obtained at X band (9.5 GHz) (14). The sample consisted of four contiguous capillaries containing a 2 mM potassium nitrosodisulfonate water solution. The potassium nitrosodisulfonate solution was freshly prepared following the method described in (15), and the stability of the solution was controlled by the EPR method at X band. A set of 64 projections sampled on 128 points was acquired in a 2D plane perpendicular to the axes of the capillaries. Figure 1B shows an example of 1D projections acquired with the following parameters: field span 0.7 mT; sweep time 8 s; time constant 30 ms; gradient value 150 mT/m; modulation amplitude 0.02 mT; modulation frequency 100 kHz; microwave power 2 mW. Three other sets of 1D projections at different S/N (S/N Å 100, 40, and 20) were obtained from the original data (see Fig. 1C), adding random noise normally distributed at different amplitudes. The adoption of this method is justified by the fact that the experimental noise data obtained with the X-band spectrometer follows a frequency distribution similar to that of white noise. An alternative approach to collecting the 1D projections at dif-
6j0e$$0893
06-05-96 02:24:12
qÅ10
magas
l1 , l2
[8]
where l1 is the height of the overlapping and l2 is the height of one of the two lines. Maximum resolution corresponds to q Å 1. The second parameter, h, measures the noise amplification of the deconvolution technique. It is evaluated by h Å 1 0 RMS[Nfb (x, y) 0 NfV (x, y)],
[9]
where Nf b (x, y) Å fb (x, y)/[max( fb (x, y)] is the normalized signal intensity, NfV ( x, y) Å fV ( x, y)/[max( fV ( x,y)] is the normalized DR or RD reconstructed image, and the RMS value has been evaluated on approximately 20% of the points of the image away from the useful signal. Without noise amplification the two terms are equal and h Å 1. The third parameter, j, gives an indication of how deconvolution affects the intensity ratios between different regions of the image and can be written as
j Å 1 0 maxj
SZ
*xj12,yj1 2 fb (x, y)dxdy xj ,yj
*x ,y fb (x, y)dxdy
*xj12,yj1 2 fV (x, y)dxdy xj ,yj
0
*x ,y fV (x, y)dxdy
ZD
[10]
for each region j in which the image is divided. A detailed discussion about the preceding parameters, for the one-dimensional case, is presented elsewhere (8). These three parameters, which range from zero to one, indicate a better image, relative to what they represent, as their values increase. The
AP: Mag Res
63
NOTES
FIG. 2. Two-dimensional images of the cross section of the phantom containing a 2 m M potassium nitrosodisulfonate water solution. The phantom is composed of four contiguous capillaries whose dimensions are: (1) 0.3 mm i.d., 1.2 mm o.d.; (2) 0.3 mm i.d., 1.2 mm o.d.; (3) 0.5 mm i.d., 1.5 mm o.d.; (4) 0.2 mm i.d., 1.4 mm o.d. The 2D images, obtained from the projections at S/N Å 100, were reconstructed by the sequences DR (A) and RD (B).
values of the three parameters, relative to the final images of the X-band experiments, are reported in Table 1. Additional experimental data were acquired with a radiofrequency EPR spectrometer (280 MHz) (16). The choice of a low-frequency experiment is strictly connected with the EPRI necessity of using deconvolution when the S/N ratio is low. This low S/N ratio requires the use of small field gradients and, consequently, loss of resolution for large samples. The sample consisted of two small glass flasks, one round-bottomed and the other point-bottomed, filled with a 1 mM potassium nitrosodisulfonate water solution. The flasks were about 6.5 cm in length and 2.6 cm in width. Both necks
TABLE 1 Values of a, b, q, h, and j at Different S/N Ratios for both DR and RD Sequences, Applied to X-Band Data Deconvolution parameters
h
q
j
S/N ratio
a
b
DR
RD
DR
RD
DR
RD
400 100 40 20
0.5 0.5 0.5 0.5
60 40 32 28
0.99 0.90 0.80 0.76
0.99 0.99 0.99 0.99
0.97 0.93 0.88 0.78
0.98 0.98 0.97 0.90
0.68 0.66 0.64 0.64
1.00 1.00 1.00 1.00
6j0e$$0893
06-05-96 02:24:12
magas
had two narrower parts and an average width of 1.6 cm. To obtain a 2D image along the axes of the flasks, a total of 64 projections, sampled on 128 points, were collected. The acquisition parameters were gradient intensity 10 mT/m; field span 1.3 mT; field center 9.5 mT; sweep time 50 s; time constant 300 ms; modulation amplitude 0.04 mT; modulation frequency 8.6 kHz. Figure 3 shows the images obtained with the DR and RD sequences. From the values of the parameters reported in Table 1 and in Fig. 2, we note that DR and RD give substantially equivalent results if the S/N ratio is high ( §400). On the other hand, using lower S/N, the image produced by the sequence DR has a lower resolution, even if the same parameter values of deconvolution have been used. This fact is due to the blurring effect on the image when the reconstruction is performed on data distorted by deconvolution, which is more pronounced as the S/N ratio becomes lower than 100. This fact is especially evident in the results from the low-frequency experiment because of the small distance between the two flasks. In fact, the narrow portions of the neck of one of the flasks are well evident in Fig. 3B but not in Fig. 3A. For the same reason, the intensity proportions between different zones of the image are more altered in DR. One important consideration is that it is difficult to increase the values of the three parameters q, h, and j simultaneously. This is because, for example, improving the resolution (q
AP: Mag Res
64
NOTES
FIG. 3. Two-dimensional images of the longitudinal section of a phantom composed of two small flasks obtained at radiofrequency (280 MHz). (A) Obtained by the DR sequence, and the calculated values of q, h, and j are 0.75, 0.70, and 0.72, respectively. (B) Obtained by the RD sequence, and the corresponding values of q, h, and j are 0.96, 0.98, and 0.93, respectively. The deconvolution parameters were a Å 0.5 and b Å 26 for both sequences. The S/N ratio of the projection data was 15.
going to one) is accompanied by an increase of spurious oscillations, i.e., a decrease of the value of h, and by a decrease of the proportion stability, i.e., a decrease of the value of j. Whenever the preceding fact is true, RD shows better resolution than DR and does not amplify the noise as DR does (see Table 1 and Fig. 3). The results presented in this Note demonstrate that, simply by changing the order of deconvolution and reconstruction, it is possible to obtain an improvement in the quality of the final result, especially for data having a low S/N ratio. One difference is that in the DR sequence the deconvolution is one-dimensional. On the other hand, in the RD sequence a bidimensional deconvolution is required. This fact implies that, in this case, some accurate deconvolution techniques usable in 1D problems only, such as spline deconvolution (8), must be rejected. Moreover, it is worth underlining that with the values of the field gradient used, for both the Xband and the low-frequency experiments, the nitroxide hyperfine lines do not overlap. We have not treated the case of line overlapping and it is not certain that the RD sequence is also better than the DR sequence in this case. REFERENCES 1. M. J. R. Hoch and A. R. Day, Solid State Commun. 30, 211 (1979). 2. G. R. Eaton and S. S. Eaton, Concepts Magn. Reson. 7(1), 49 (1995).
6j0e$$0893
06-05-96 02:24:12
magas
3. B. F. Carley and R. W. Joyner, J. Electron. Spectrosc. Relat. Phenom. 16, 1 (1979). 4. S. M. Riad, Proc. IEEE 74(1), 82 (1986). 5. G. R. Eaton, S. S. Eaton, and K. Ohno (Eds.), ‘‘EPR Imaging and in Vivo EPR,’’ CRC, Boca Raton, Florida, 1991. 6. R. H. T. Bates and M. J. McDonnell, ‘‘Image Restoration and Reconstruction,’’ Vol. 16, Oxford Engineering Science, Oxford, 1986. 7. F. Momo, S. Colacicchi, and A. Sotgiu, Meas. Sci. Technol. 4, 60 (1993). 8. G. Placidi, M. Alecci, and A. Sotgiu, Rev. Sci. Instrum. 65, 58 (1994). 9. M. Maltempo, S. S. Eaton, and G. R. Eaton, J. Magn. Reson. 72, 449 (1987). 10. M. Maltempo, S. S. Eaton, and G. R. Eaton, J. Magn. Reson. 77, 75 (1988). 11. A. V. Oppenheim and R. W. Schafer, ‘‘Digital Signal Processing,’’ Prentice–Hall, Englewood Cliffs, New Jersey, 1975. 12. R. M. Mersereau and A. V. Oppenheimer, Proc. IEEE 62, 1319 (1984). 13. R. A. Brooks and G. Di Chiro, Phys. Med. Biol. 21, 689 (1976). 14. G. Gualtieri, G. Placidi, A. Sotgiu, and S. Del Monaco, Rev. Sci. Instrum. 66, 3715 (1995). 15. H. M. Swartz, J. R. Bolton, and D. C. Borg, ‘‘Biological Applications of Electron Spin Resonance,’’ Wiley Interscience, New York, 1972. 16. M. Alecci, S. Della Penna, A. Sotgiu, L. Testa, and I. Vannucci, Rev. Sci. Instrum. 63, 4263 (1992).
AP: Mag Res