Activity and attenuation recovery from activity dataonly in emission computed tomography

Activity and attenuation recovery from activity dataonly in emission computed tomography

Alvaro R. de Pierro^*; Fabiana Crepaldi^**

State University of Campinas, Department of Applied Mathematics, CP. 6065, 13081-970 Campinas, SP, Brazil, E-mail: [e-mail]

ABSTRACT

We describe the continuous and discrete mathematical models for Emission Computed Tomography (ECT) and the need for attenuation correction. Then we analyse the problem of retrieving the attenuation directly from the emission data, nonuniqueness and its consequences. We present the existing and new approaches for solving the problem. Methods are compared and illustrated by numerical simulations. The presentation focuses on Positron Emission Tomography, but we also discuss the extensions to Single Photon Emission Computed Tomography.

Mathematical subject classification: 44A12, 65R30, 65R32.

Key words: Attenuated Radon Transform, maximum likelihood, emission tomography.

1 Introduction

Emission Computed Tomography (ECT), in its two modalities, SPECT and PET (Single Photon and Positron Emission), is the main tool to understand physiological processes inside the body at a molecular level. The pure mathematical model (no attenuation, no scattering, high statistics, etc) for ECT is the Radon transform, the same as for X-ray Computed Tomography (CT). In this article we describe the main issues concerning the introduction of attenuation in the ECT model and the different approaches to the direct computation of the activity and the attenuation maps from emission data only. Our emphasis will be on iterative methods and Maximum Likelihood models for PET, although, it is worth mentioning that each method for PET generates a similar approach for SPECT and viceversa. We discuss the nonuniqueness issue for PET and present some new methods for solving the problem as well as numerical simulations. Open problems and research directions are discussed in Section 5.

1.1 What is Emission Computed Tomography (ECT)

In ECT we are trying to understand a physiological process, like glucoseconsumption, blood flow, oxidative metabolism, neurotransmission, and many others (see [17] for a long list for PET). How this is done? There is a compound that is metabolized during that physiological process, and that could be tagged using a radioative isotope. In the case of PET, the isotope is artificial and decays emitting positrons. After a very short path, the positron meets an electron, generating two photons in almost opposite directions, that are detected by pairs of detectors. Each pair of detectors determines a line and the data is the total number of coincidences counted along all possible lines (Figure 1 shows a PET scanner with its ring of detectors).

We aim at reconstructing the emission density that is correlated with the process being studied. The most common example in PET is glucose, that is tagged with F₁₈, giving rise to FDG (fluorodeoxyglucose). This measures glucose consumption, an indirect description of many processes for clinical diagnosis (malignancy, heart failure, etc) as well as medical reserach (brain function, blood flow, etc). In the case of SPECT [5], the isotope is natural and decays by generating a single photon (see Figure 2), so, the line is determined by using a collimator.

Many photons are absorbed by the collimator before detection, thus obtaining a very low statistics, poorer than PET.

1.2 A continuous 'pure' mathematical model

In a continuous setting, we aim at reconstructing a function f(x) representing the mean emission density at the point x given the integrals along lines L, that is,

where d_L is the mean number of coincidences detected by the pair of detectors (or by a single detector and its collimator for SPECT) that determine L. Equation (1.1) means that the 'pure' problem reduces to invert the well known Radon Transform [13]

1.2.1 The difference between ECT and Computed Tomography (CT)

We observe that the inversion above is the basic model also for Transmission Computed Tomography (CT), where d_L, instead of number of coincidences stands for log , being in this case and the mean number of photons emitted by a source and the detected ones respectively.

Inversion formulas for (1.2) can be easily obtained using the Fourier Projection Theorem (see [13]). However, the higher levels of noise in ECT, as compared to CT, made models that take into account noise characteristics and the use of iterative methods more reliable, producing better images [25]. This leads to to the necessity of discretization, in order to develop Maximum Likelihood (ML) models.

1.3 A discrete model

We consider now the discretized two-dimensional PET model. In what follows, x_j denotes the expected number of coincidences per unit area in the jth pixel (j = 1,... , n) (we are assuming here that the ''tubes'' are of constant cross-section and that dividing the radioactivity concentration per unit volume by the area of this constant cross section gives us the ''number of coincidences per unit length''); x represents the n-dimensional vector whose jth component is x_j and will be referred to as the image vector (see Figure 3).

Suppose we count coincidences along m lines. We denote by g the m-dimensional vector whose ith element, g_i (g_i > 0), is the number of coincidences which are counted for the ith line (i = 1,..., m) during the data collection period; we shall refer to g as the measurement vector. If a_ij (a_ij > 0) denotes the probability that a positron emitted from pixel j results in a coincidence at the i-th line of response (in this work, as usual, this probability is approximated by the weighted length of the intersection), then g_i is a sample of a Poisson distribution whose expected value is

where á, ñ is the standard inner product and aⁱ is the ith column of the transpose A^T of the m × n projection matrix A = (a_ij) [28].

The probability of obtaining the measurement vector g if the image vectoris x (i.e., the likelihood function) is

The reconstruction problem is to estimate the image vector x given the data measurements g.

As mentioned before one approach to this problem, first suggested in [24], is the ML method which estimates the x that maximizes P_L(g|x) subject to nonnegative constraints on x, or, equivalently, finds the x > 0 which maximizes

In [26] and [16], the Expectation Maximization (EM) algorithm was proposed for the maximization of (1.5). This results in an iterative algorithm, which, starting with a strictly positive vector x⁽⁰⁾, successively estimates the image vector by

The EM algorithm's convergence rate is extremely slow and a significant amount of computation is often needed to achieve an acceptable image, (essentially, applying filtered backprojection costs one iteration of the EM, and, in general, at least fifty iterations are required). In the 90's faster methods, like ordered subsets [14] and RAMLA [4], made possible the adoption of iterative methods by commercial scanners.

2 The need for attenuation correction: attenuated transforms, problems

In practice, a more accurate model for ECT should be considered because of the attenuation effect, that is, the fact that a substantial number of photons is absorbed by the tissue, and the formulas in (1.1) should be substituted by(see [21])

for PET and

for SPECT, where µ is now the attenuation density and d_L as before. Considering that the previous formulas are now strongly nonlinear, the usual way to compensate for atenuation is to estimate it by a previous (or simultaneous) CT scan for estimating µ. An important improvement could be achieved if it is possible to retrieve the attenuation map directly from the activity data (meaning less time, less patient motion, automatic registration, no error amplification,etc). Recently, PET scanners coupled with CT scanners, called PET-CT, have been developed, overcoming the problem (at a very high cost). However, still many of the existing scanners are not PET-CT, and there is not such a capability in most of the gamma-camera based systems [22] with coincidence detection (SPECTscanners adapted for coincidence detection), and in small high resolution tomographs for research purposes [19].

Methods trying to perform this task, without CT corrections, tend to retrieve solutions that generate an undesired cross talk between the activity and attenuation images, unless a considerable amount of information about attenuation values is assumed as in [20]. That is, visually, the main drawback is that the shadow of the attenuation appears in the emission image and vice versa. This is a consequence of nonlinearity and nonuniqueness of the inversion of the operators just described. In the following we analyse the problem and we describe the existing and new approaches for its solution, now restricted to PET ((2.7)).

3 Recovering the attenuation from activity data only: nonuniqueness and different approaches

It is very easy to show that, for data generated by circularly symmetric activity and attenuation images, the inversion in (2.7) has an infinite number of solutions. Also, examples of the influence of 'approximate' nonuniqueness could be found in [2]. In any case, methods based on consistency of the data, first proposed by Natterer [20] have been the most effective ones, producing the best images. In the following µ and l will stand for the images corresponding to the attenuation and the activity values respectively. The results of the numerical experiments corresponding to the methods described in 3.2 and 3.3 are describedin Section 5.

3.1 Natterer

For the sake of completeness, and because it was the first to introduce in the problem the necessity of data consistency, we briefly describe in this section Natterer's approach [20]. As we said before, this approach relies heavily on the knowledge of additional information on the attenuation, but, on the other hand, it is the point of departure of the other successful approaches based on consistency. From the Helgason-Ludwig conditions for PET data [11], [18], if we take the Fourier Transform of the m-th moment, all coefficients greater than m should be zero, that is

and T(s,q) is the Radon Transform of µ. g is the projection data as before. (3.9) determines a set of equations to be satisfied by the attenuation. However, these equations are not enough to completely determine µ. So, the idea is to use additional information in order to reduce the indetermination. To do this, Natterer et al. [29] consider that the sought attenuation is a shifted version of a given one. That is, they consider µ(x) = µ₀ (Ax+a) where A is a 2x2 matrix and a a 2-vector. If the attenuation is not too different from µ₀ sometimes it works fine.

3.2 Bronnikov

In this Section we present a method for the attenuation recovery from PET emission data, that follows closely Bronnikov's approach for SPECT. Our implementation uses known properties of pseudoinverses, allowing a relatively low cost reconstruction, as compared to Bronnikov's implementation for SPECT.We show experimentally that, even without regularization, the obtained images are stable. Following [7, 6], the discretization of the attenuated Radon transform for PET (2.7) can be described by the equation:

where A_µÎ ^m×n represents the projection matrix, m Î ⁿ, the attenuation vector, l Î ⁿ, the activity, and g Î ^m the emission data vector as before. As in [7], µ appears as a matrix parameter. The elements of the matrix A_µ = (a_ij) are now defined by:

where s_ij and s_ik stand for the intersection length of the ith ray with the jth and kth pixels of activity and attenuation, respectively. In our approach, for the sake of simplicity, we consider the same discretization for both, activity and attenuation.

The number of equations m in (3.10) should be greater than or equal to the number of variables, that is, the number of attenuation voxels (l) plus the number of activity voxels (n): m > n + l. In the particular case where the images have the same resolution, (l = n), Bronnikov observed experimentally that the choice m = 3n is the more efficient one. For a general Hilbert space H, let be a closed subspace of H and ^{^} its orthogonal complement. If we assume that g Î H, it can be written as

where Î and g_^Î ^{^}. In our case is ^m and is the range of the operator defined by (3.10). So, let R(A_µ) be the range of the operator A_µ and R(A_µ)^{^}. Its orthogonal complement. then ^m = R(A_µ) Å R(A_µ)^{^} and we can write g = +g_^, where Î R(A_µ) and g_^ Î R(A_µ)^{^} and g Î ^m. Then the consistency condition for (3.10) is that the data should not have components other than those in the operator's range, that is,

If Î ^n×m stands for the Moore-Penrose pseudoinverse of the matrix A_µ [3], then:

where

are the orthogonal projections onto R(A_µ) and R(A_µ)^{^}, respectively.

With the previous notation, the consistency condition assumes the form:

In the case of PET, A_µ is the product of two matrices: a diagonal matrix, that depends on the attenuation, and a projection matrix, S, that does not, that is:

with

Considering that, in this particular case (D_µ- diagonal and positive), (D_µ-S)⁺ = S⁺(D_µ)^-1, and replacing (3.17) in (3.16), we obtain

where

Therefore the consistency condition is equivalent to:

In system (3.20), the matrix D_µ- can be eliminated, by multiplying by the inverse. We have, then, a 'linear system' for x = D_µg, nonlinear for µ, that can be written in a more convenient form as:

with

M = (I - SS⁺)x_i + u_ig_i, x_i > g_iu_i =

In the next section we describe two different approaches for solving (3.21)

3.2.1 Using the Newton-Raphson method

System (3.21) can be solved using the Newton-Raphson method for non-linear systems with nonnegative constraints, defined by the iteration:

where

(x is given by (3.21)) and J_µ is the Jacobian matrix with components:

J_m or:

In matrix form,

where B_µ is a matrix with components defined by

Using the QR decomposition and (3.21), the matrices S and M can be written as:

with Q₁Î ^m×n and Q₂Î ^m×(m-n). As in the original Bronnikov's approach, the Jacobian matrix is singular with rank m - n. An alternative, equivalent, but more expensive way to obtain M is the SVD decomposition. In [9], details could be found on the different ways to deal with the matrix M.

3.3 The ML Approach. Crepaldi and De Pierro

According to our analysis in (1.3), the ML approach, that takes into account the statistical properties of the noise, tends to retrieve better solutions than nonstatistical models. So, back to the discretized model, given the activity and the attenuation maps for PET, where the attenuation is independent of the position along the projection line, the expected number of counts is given by

with

and

with s_ij being the probability that an annihilation occurred at j is detected in tube i and s_ij the intersection of tube i with box k:

Now, if g_i is the number of counts measured at the i-th tube of detectors, as in Section 1.3 the g_i's, are independent Poisson variables with mean r_i. So, after excluding terms not containing l or µ, the new log-likelihood function that we want to maximize is

or equivalently,

It can be proven that, for fixed µ, L is a concave function of l and vice versa [23]. However, L is not jointly concave for (l, µ). Following [23], we develop an alternate maximization method, where each maximization step is based on the minorizing functions approach proposed in [10]. The difference between the algorithm in [23] and ours is essentially the guaranteed monotonicity for the latter. If we assume first that the attenuation map is constant, using the logarithm's concavity, and, taking into account that

where l^k is the current emission at iteration k and, for any given vectors u = (u)_j and v = (v)_jáu, vñ = S_ju_jv_j stands for the scalar product, we have that:

is greater equal than

On the other hand, using the same procedure to separate the µ variables for the (negative) exponential term. That exponential part is greater or equal than

Using (3.36) and (3.37), we obtain, for each k:

with

Now, the algorithm's global iteration is defined in two steps: for a given activity-attenuation vector pair (l^k, µ^k), the first step consists of solvingthe problem

and the second step iterates in µ solving the problem

The solution of the first problem (3.40) is nothing but the EM algorithm applied to L with fixed µ and this gives the updating

with

The solution of the problem defined by the second step (3.41) requires solving a set of nonlinear optimization problems in one variable that can be very easily solved using the Newton-Raphson method with a positivity constraint. In our implementation, enough iterations of Newton-Raphson's method were performed in such a way that the maximum in (4.41) is guaranteed (in practice one iteration of Newton Raphson is enough to increase the function value) as well as the monotonicity of the algorithm (4.40-41). In [15], l and µ are updated sequentially by applying Newton's method to the first order necessary conditions for the maximum. Trial and error is needed to choose the ascent parameters, in order to preserve monotonicity, as in [23].

Unfortunately, the application of the algorithm above in its pure form, as well as Nuyts', still generates images with the undesired 'cross-talk', so, following Natterer's original idea, consistency should be imposed. The way we found to do that is by using Iterative Data Refinement.

3.3.1 Iterative Data Refinement (IDR)

In this Section we briefly describe the general methodology known as Iterative Data Refinement (IDR). A more detailed description together with some historical background and applications in tomography could be found in [8].

Let x be a mathematical representation of a physical object. It might represent, in our problem, the activity density. The idealized measurer is defined as the operator that provides the ideal measurement vector, , when applied to the vector x:

x =

The matematical operator that recovers the input from the output is called the recovery operator and denoted by :

= .

In general is a good approximation for x.

However, in practice we do not have an exact physical implementation of the idealized measurer . We have an actual measurer which, given x, produces the actual measurement vector :

x = .

The problem arises because, in general, = is not a good approximation for x.

This difficulty is the motivation for Iterative Data Refinement (IDR). The main idea is: if y^k is a reasonable approximation for , then

y^k - h(k)

y^k » - h(k)

where h(k) is a relaxation parameter. From this equation, it is natural to propose the iterative process:

Suppose that y^k is a good approximation for . Then

y^k »

y^k.

and (3.44) becomes

In the next Section we describe how we apply IDR to the activity attenuation problem.

3.3.2 IDR and the optimization approach

The main numerical problem that we have to solve for the optimization model is the fact that iterative methods like (3.40)-(3.41) tend to choose the 'wrong' solution (fixed point). Solutions showing the undesired 'crosstalk' look more equilibrated and highly likely than the ones we are looking for. This rationale plus the observation of the results of many experiments led us to the introduction of a multiplicative factor in (3.42), that, essentially ''breaks'' the equilibrium, but still converging to an ML solution. That is, we substituted (3.42) by

where a_l is as before and

is the multiplicative factor that retains the asymptotic convergence properties of the EM part of the algorithm, x is a constant parameter to be chosen. f_x(k) tends to one, so that the monotonicity property of the likelihood remains essentially valid when updating the activity variables. An important observation is that the multiplicative factor enhances the data, preparing a better starting point for the next step, IDR.

In order to apply IDR to our problem, now x represents the vector whose elements are the activity and attenuation densities, that is, x = (l,µ), is the operator that we want to 'invert', that is, the discrete Radon transform with the corresponding attenuation, given by the left hand side of (2.7), and the recovery method of simultaneous reconstruction, previously described (the alternated maximization using the multiplicative factor for the EM step). In the appropriate notation, we can write (3.45) as

Observe that, in (3.48), if there is convergence to a fixed point, say , the ideal measurement (and that is the case), that means that the data has been modified in such a way that our activity-attenuation solution is one that the forward operator (the attenuated Radon Transform) applied to it is consistent with the original data (actual measurement).

Summarizing, our method consists of two iterative processes: the IDR iterations and the alternated iterations from section 2 that recover both images simultaneously. This can be described by the following steps; for a given starting point, and the data , and, in general, for a point (l^k, µ^k) and the modified data g^k:

Step 1: Apply (4.41-42) until convergence using g^k as data, that is, compute g^k.

Step 2: Update the data in (4.45) using the forward operator .

Step 3: Go back to Step 1. Stop if a suitable criterion is satisfied.

In the next sections we present our simulations and describe how the various parameters have been chosen, including the value of the constant x in the multiplicative factor and the number of iterations used for each method, the alternated maximization and IDR.

4 Simulations

In this Section we present a visual comparison of the methods described in Sections 4.2 and 4.3, and the results for pure ML approaches without IDR. That is, we consider only methods not adding any information about the solution, through the use of approximate. For the simulations we used a simple model of a 2D thorax cross section phantom. The activity source was modeled as a ring, representing the heart, inside an ellipse, which represents the background activity in the human body, approximately 10 times less than the ring's activity. The attenuation map was an ellipse representing the body, two ellipses representing the lungs and a small circle for the spine bone (see Figure 4). For all the 50×50 reconstructions it is also shown the profile of the the 35th line and 16th line for 32×32 reconstructions.

For all the experiments, the starting images were uniform using an average value; there were no modifications in the results when varying these values in an acceptable range. Images were reconstructed using a 50×50 grid, except for Bronnikov's, that was 32×32, because of the huge amount of computations.

Figures 5 and 6 show the reconstructions for Nuyts et al's method and pure Crepaldi & De Pierro's (without IDR) respectively. The cross-talk between images is clearly seen.

Figure 7 shows the result of applying Bronnikov's approach (that computes only the attenuation) and Figure 8 the result for Crepaldi & De Pierro's method with IDR. In both cases the cross-talk almost dissapeared.

Finally Figure 9 shows the result obtained by Bronnikov's approach withnoise. Noise was generated by using a Poisson generator with mean given by the data for each projection. Our results with noise for IDR did not show substantial changes.

5 Some conclusions and questions still to be answered

The methods presented in Sections 4.2 and 4.3 are a first step to obtain quality images in ECT from emission data only, without any additional information. However, there is a relatively long list of open problems still to be addressed.It follows.

1. One advantage of iterative methods with IDR, as compared to Bronnikov's approach is the possibility of considering ML models for the noise, but the other one is the computational burden. Iterative methods 'could be' computationally much less time consuming if faster methods like RAMLA [4] are used for the inner part in (3.40)-(3.41). The overall number of iterations (including IDR) should be also reduced.

2. Reasonable stopping criteria are needed for the outer (IDR) and inner iterations.

3. Stability: together with the number of iterations, several other parameters are involved. Also, there is some dependence of these parameters with respect to the discretization.

4. Some kind of regularization should be considered. That could help to improve the images in the presence of noise, but could also help for the choice of the 'right' solution.

5. Local convergence results could be obtained along the lines of [12].

6. We are currently working on simulations for more complex objects, for example, nonconvex. Preliminary results show that our approach with IDR does not produce a convexification as in [23].

7. We are now implementing the methods for SPECT simulations, and with real data, in a joint research with the Heart Institute of the University of São Paulo.

It is also worth mentioning that these methods for problems with two 'separated' unknowns and many possible solutions depending on the data, could be a reference in other important inverse poblems like optical tomography and phase retrieval.

[17] http://www.crump.ucla.edu/software/lpp/lpphome.html

Received: 19/XII/06. Accepted: 19/XII/06.

#698/06.

• [1] M. Avriel, Nonlinear Programming: Analysis and Methods, Prentice Hall, Englewood Cliffs, NJ, 1976.
• [2] C. Bai, P. Kinahan, D. Brasse, C. Comtat, D.W. Townsend, C.C. Metzer, V. Villemagne, M. Charron and M. Defrise, An analytic study of the effects of attenuation on tumor detection in whole-body PET oncology imaging, J. Nucl. Med., 44 (2003), 1855-1861.
• [3] A. Ben Israel and T. Greville, Generalized Inverses: Theory and Applications, Springer, 2 edition, 2003.
• [4] J.A. Browne and A.R De Pierro, A row-action alternative to the EM algorithm for maximizing likelihoods in emission tomography, IEEE Transactions on Medical Imaging, 15(4) (1996), 687-699, October.
• [5] T.F. Budinger, G.T. Gullberg and R.H. Huesman, Emission computed tomography, in Image Reconstruction from Projections: Implementation and Applications, G.T. Herman (Ed), Springer Verlag, Berlin, Heidelberg. 5 (1979), 147-246.
• [6] A.V. Bronnikov, Numerical solution of the identification problem for the attenuated Radon transform, Inverse Problems, 15 (1999), 1315-1324.
• [7] A.V. Bronnikov, Reconstruction of attenuation map using discrete consistency conditions, IEEE Trans. Medical Imaging, 19(5) (2000), 451-462.
• [8] Y. Censor, T. Elfving and G.T. Herman, A method of iterative data refinement and its applications, Math. Meth. Appl. Sc., 7 (1985), 108-123.
• [9] F. Crepaldi, On the computation of the attenuation and the activity in emission tomography from activity data, PHd Dissertation, Institute of Physics, State University of Campinas, 2004.
• [10] A.R. De Pierro, A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography, IEEE Trans. Med. Imaging, 14(1) (1995), 132-137.
• [11] S. Helgasson, The Radon Transform, Boston, Birkhauser, 1980.
• [12] E.S. Helou and A.R. De Pierro, Convergence results for scaled gradient algorithms in positron emission tomography, Inverse Problems, 21(6) (2005), 1905-1914, October.
• [13] G.T. Herman, Image Reconstruction from Projections: The Fundamentals of Computerized Tomography, Academic Press, New York, NY, 1980.
• [14] H.M. Hudson and R.S. Larkin, Accelerated image reconstruction using ordered subsets of projection data, IEEE Trans. Med. Imaging, 13(4) (1994), 601-609.
• [15] M. Landmann, S.N. Reske and G. Glatting, Simultaneous iterative reconstruction of emission and attenuation images in positron emission tomography from emission data only, Med. Phys., 29 (2002), 1962-1967.
• [16] K. Lange and R. Carson, EM reconstruction algorithms for emission and transmission tomography, J. Comput. Assisted Tomog., 8 (1987), 306-316.
• [18] D. Ludwig, The Radon Transform on euclidean space, Comm. Pure Appl. Math. XIX (1966), 49-81.
• [19] R.H. Huesman and T.F. Budinger, Design of a high resolution, high sensitivity PET camera for human brains and small animals, IEEE Trans. Nucl. Sci. 44 (1997), 1487-1491.
• [20] F. Natterer, Determination of tissue attenuation in Emission Tomography of optically dense media, Inverse Problems, 9 (1993), 731-736.
• [21] F. Natterer, The Mathematics of Computerized Tomography, SIAM, 2001.
• [22] P. Nelleman, H. Hines, W. Braymer, G. Muehllehner and M. Geagan, Performance characteristics of a dual head SPECT scanner with PET capability, IEEE Nuclear Science Symposium and Medical Imaging Conference, Conference record, (1995), 1751-1755.
• [23] J. Nuyts, P. Dupont, S. Stroobants, R. Benninck, L. Mortelmans and P. Suetens, Simultaneous Maximum A Posteriori reconstruction of attenuation and activity distributions from emission sinograms, IEEE Trans. Med. Imaging, 18(5) (1999), 393-403.
• [24] A.J. Rockmore and A. Macovski, A maximum likelihood approach to emission image reconstruction from projections, IEEE Trans. Nucl. Sci., 23 (1976), 1428-1432.
• [25] M.S. Rosenthal, J. Cullom, W. Hawkins, S.C. Moore, B.M.W. Tsui and M. Yester, Quantitative SPECT imaging: a review and recommendations by the focus committee of the Society of Nuclear Medicine Computer and Instrumentation Council, J. Nucl. Med., 36 (1995), 1489-1513.
• [26] L.A. Shepp and Y. Vardi, Maximum likelihood reconstruction for emission tomography, IEEE Trans. Med. Imaging, 1 (1982), 113-121.
• [27] M.M. Ter-Pogossian, M. Raichle and B.E. Sobel, Positron emission tomography, Scientific Amer., 243(4) (1980), 170-181.
• [28] Y. Vardi, L.A. Shepp and L. Kaufman, A statistical model for positron emission tomography, J. Amer. Statist. Assoc., 80 (1985), 8-20.
• [29] A. Welch, C. Campbell, R. Clackdoyle, F. Natterer, M. Hudson, A. Bromiley, P. Mikecz, F. Chillcot, M. Dodd, P. Hopwood, S. Craib, G.T. Gullberg and P. Sharp, Attenuation correction in PET using consistency information, IEEE Trans. Nuclear Science, 45(6) (1988).

*

Work by this author was partially supported by CNPq Grants No. 300969/2003-1 and 476825/2004-0 and FAPESP Grant No. 2002/07153-2, Brazil

**

Work by this author was partially supported by FAPESP Grant No. 04/07238-3R Brazil

Publication Dates

• Publication in this collection
  19 Mar 2007
• Date of issue
  2006

History

• Accepted
  19 Dec 2006
• Received
  19 Dec 2006