Image reconstruction by conditional entropy maximisation for PET system P.P. Mondal and K. Rajan Abstract: The authors show that the conditional entropy maximisation algorithm is a generalised version of the maximum likelihood algorithm for positron emission tomography (PET). Promising properties of the conditional entropy maximisation algorithm are as follows: an assumption is made that the entropy of the information content of the data should be maximised; it is a consistent way of selecting an image from the very many images that fit the measurement data; this approach takes care of the positivity of the reconstructed image pixels, since entropy does not exist for negative image pixel values; and inclusion of prior distribution knowledge in the reconstruction process is possible. Simulated experiments performed on a PET system have shown that the quality of the reconstructed image using the entropy maximisation method is good. A Gibbs distribution is used to incorporate prior knowledge into the reconstruction process. The mean squared error (MSE) of the reconstructed images shows a sharp new dip, confirming improved image reconstruction. The entropy maximisation method is an alternative approach to maximum likelihood (ML) and maximum a posteriori (MAP) methodologies. 1 Introduction Entropy maximisation is a general tool for solving a large number of inverse problems where the possibility of direct solution is ruled out. Defining an entropy measure in the context of image reconstruction was originally due to Friden [1, 2]. The maximum likelihood (ML) algorithm plays a vital role in emission tomography (ET) [3, 4]. Though a few generalisations of ML estimation, such as Bayesian approaches to ET do exist, we propose to take an entropybased approach to ET. It is shown theoretically that the proposed algorithm leads to another generalisation of the ML algorithm for ET. ET is an imaging technique potentially useful in the study of human physiology and organ function. The basic aim of ET is to obtain a quantitative map of the spatial and temporal distribution of radionuclides inside the human body which represent functional activity such as pain [5], motor action [6], etc. This is achieved by realistic modelling of a physical phenomenon involved along with the event measurement count of gamma rays either in coincidence or singly, depending upon the modality used. For positron emission tomography (PET), coincidence counts are detected due to the positron –electron annihilation process, while a single count is detected in single photon emission tomography (SPECT). The maximum likelihood estimation (MLE) algorithm based on expectation maximisation (EM) produces a sequence of images that converge to the maximum likelihood estimate of the emission density. The quality of the image produced by MLE is found to be better than that generated by deterministic algorithms such as convolution back projection (CBP) [3, 4, 7], especially with low measurement count data. The ML algorithm has not been found suitable for routine clinical application owing to its excessive computational requirements and the deteriorating nature of the reconstructed images [8]. Time constraint is not a serious concern because it can be reduced considerably by parallelising the algorithm [9, 10] and by using fast digital signal processor (DSP) chips as processing elements [11]. Various algorithmic methods have been suggested to accelerate the convergence of the ML algorithm [12 – 16]. In particular, one can accelerate the image reconstruction process using an ordered subset (OS) method [17]. Many variants of the ordered subset algorithm, such as OSEM [18], row –action maximum likelihood (RAMLA) [19], rescaled bock iterative (RBI-EMML) [20] and RBI-MAP [21] have been proposed in the literature to reduce the number of iterations. It is well known that as the iterations of the ML algorithm continue, estimates become more and more noisy. This has been termed ‘dimensional instability’ and occurs due to the application of unconstrained maximum likelihood estimation of density functions based on a point process [22]. To overcome this serious problem, two approaches have been suggested: (i) The iteration has to be stopped before the quality of the reconstructed image starts deteriorating [8]. (ii) Inclusion of prior distribution knowledge of the emission densities in the reconstruction process [23, 24]. Obviously, approach (ii) is more promising, as this enables the possibility of improved reconstruction over the ML algorithm. In view of this possibility, a generalised version of the ML algorithm was developed [25, 26]. This algorithm considers conditional entropy as a cost function rather than a likelihood function. The motivation behind conditional entropy maximisation is that the resulting image has minimum configurational information so that there must be information in the measurement data for any structure to be seen [27, 28]. Maximum entropy estimation provides a consistent way of selecting a single image from the very many images that fit the measurement data. It is found that maximum entropy estimation does not introduce correlation in the image beyond that required by the data [29 – 31]. Also this addresses the important factor of the positivity of image pixels, since entropy does not exist for negative image pixels. Moreover, this method provides a handle for the possible inclusion of prior knowledge in the image reconstruction process. 2 Maximum entropy reconstruction for positron emission tomography Byrne [32, 33] has shown that maximisation of objective functions such as likelihood, posterior, entropy etc. is equivalent to minimisation of functionals of the form FðlÞ P ¼ aKLðy; QlÞ þ ð1  aÞDðlÞ; where KLðy; QlÞ ¼ n yn logðyn =½Qln Þ þ ½Qln  yn is the Kullback – Leibler (KL) distance between two non-negative vectors, i.e. the actual count y and the expected counts Ql; DðlÞ is some penalty function that forces a degree of smoothness upon l. Minimising FðlÞ by choosing DðlÞ ¼ KLðP; xÞ for l  0 and 0 < a 1 produces the iterative MAP algorithm of Lange et al. [14], where PðlÞ is a prior estimate of l: Titterington [34] proved that maximising the likelihood function derived from a Poisson model is equivalent to minimising KLðy; QlÞ: Instead, given the prior estimate PðlÞ; minimising the functional KLðl; PðlÞÞ; subject to non-negativity ðl  0Þ and data consistency ðQl ¼ yÞ is equivalent to entropy maximisation [32, 33]. The basic idea behind the conditional entropy maximisation algorithm is as follows: let Y be the observable measurement space and y 2 Y be the measurement data vector. The probability density function (or the likelihood function) for the measurement process Y is PY ðy; lÞ; defined over the space Y where l 2 L is an unknown parameter vector to be estimated. In general, PY ðy; lÞ is complicated to analyse and the estimation of l by direct maximisation of the likelihood function PY ðy; lÞ is difficult. The entropy maximisation formulation postulates the existence of an unobservable space X which is related to the space Y by a many to one mapping x ! yðxÞ taking each element of x 2 X to an element yðxÞ 2 Y: It often happens that the number of parameters needed for an adequate description is considerably greater than the number of available measurements and so there are an infinite number of possible solutions. Given this formulation, we have to pick the best in some sense. The procedure begins with the assignment of entropy measure to indicate the uncertainty in the statistical environment. Each measurement or piece of prior information reduces the randomness and restricts statistical freedom. Once a given model is adopted, the maximum entropy approach proceeds by finding that reconstruction which maximises the entropy measure and at the same time remains compatible with both the measurement data and a priori available information. In situations where the data are insufficient to eliminate all uncertainty, the principle of maximum entropy can be used to assure that the admittedly incomplete description of the situation reflects only the statistical structure imposed by the available information. We propose to use this approach for the PET system, but this will remain valid for other modalities of emission tomography such as SPECT. The measurements in PET, yj ; j ¼ 1; . . . ; M are modelled as independent PN Poisson random variables with mean parameters i¼1 li pij for j ¼ 1; . . . ; M; where li ; i ¼ 1; . . . ; N are the mean parameters of the emission process and pij is the probability that an annihilation in the ith pixel is detected in the jth detector. The conditional probability for observing Y ¼ y given the emission parameter L ¼ l is the joint probability of the individual Poisson process expressed as "  P   PN  yj # M Y exp  Ni¼1 li pij i¼1 li pij Pðy=lÞ ¼ ð1Þ yj ! j¼1 Now, the conditional entropy of the measurement process Y ¼ y given the object image L ¼ l is given by " # X X Hðy=lÞ ¼  Pðy=lÞ log Pðy=lÞ ð2Þ PðlÞ y l Conditional entropy measures the information content of y when l is known and our aim is to maximise the information content in the data set when the object image is given. So, the first derivative of the cost function with respect to li ; i ¼ 1; . . . ; N produces " # @Hðy=lÞ @PðlÞ X Pðy=lÞ log Pðy=lÞ ¼ @li @li y " !# @ X ð3Þ  PðlÞ Pðy=lÞ log Pðy=lÞ @li y Expanding and simplifying (3) we obtain X @PðlÞ @Hðy=lÞ ¼ ðPðy=lÞ log Pðy=lÞÞ @li @li y þPðlÞPðy=lÞð1 þ log Pðy=lÞÞ @ ðlog Pðy=lÞÞ @li ð4Þ The Kuhn –Tuker condition [35] for li to be the optimal solution of the nonlinear optimisation problem is li @Hðy=lÞ @li @Hðy=lÞ @li li ¼l^ i li ¼l^ i ¼ 0; if li > 0 ð5Þ 0; if li ¼ 0 ð6Þ for i ¼ 1; . . . ; N: Solution of the maximisation problem can be obtained by solving the set of nonlinear equations (5) and (6). Substituting the derivative of Hðy=lÞ from (4) in (5), we obtain X y @Pðli Þ ðPðy=lÞ log Pðy=lÞÞ þ Pðli ÞPðy=lÞ l^ i @li ð1 þ log Pðy=lÞÞ @ ðlog Pðy=lÞÞ @li li ¼l^ i ¼0 ð7Þ Since we have only one realisation of the measurement, (7) becomes @Pðli Þ ðPðy=lÞ log Pðy=lÞÞ þ Pðli ÞPðy=lÞ l^ i @li ð1 þ log Pðy=lÞÞ @ ðlog Pðy=lÞÞ @li li ¼l^ i Differentiating (1) with respect to li produces ¼0 ð8Þ " # M M X X yj pij @ pij þ ðlog Pðy=lÞÞ ¼  PN @li i¼1 li pij j¼1 j¼1 ð9Þ Substitution of (9) in (8) gives @Pðli Þ logPðy=lÞþli Pðli Þð1þlogPðy=lÞÞ @li # ! M M X X yj pij li Pðli Þð1þlogPðy=lÞÞ pij ¼0 PN i¼1 li pij j¼1 j¼1 l ¼l^ Pðy=l^ Þ li i i ð10Þ Since, Pðy=l^ Þ 6¼ 0; the second product term of (10) becomes zero. Rearranging the second term we obtain l^ i l^ i ¼ P Pðl^ i Þð1 þ log Pðy=li ÞÞ M j¼1 pij " @Pðli Þ ^ ^ log Pðy=li Þ @li li ¼li # M X y p j ij ð11Þ þ Pðl^ i Þð1 þ log Pðy=l^ i ÞÞ PN ^ i¼1 li pij j¼1 A closed form formula for solving (11) does not exist; it is necessary to devise schemes which will estimate a solution iteratively. Heuristically, a simple iterative algorithm for kþ1 solving (11) can be built by considering l^ i on the leftk hand side of (11) as the new estimate and l^ i on the righthand side as the old estimate. kþ1 l^ i ¼ k l^ i    PM P l^ i 1 þ log P y=lki j¼1 pij "    @P lki ^k li ¼l^ i log P y=li @li !# ð12Þ We call (12) the entropy maximisation algorithm for PET, which recursively produces a series of reconstructed images that converge to the maximum entropy estimate. k It is instructive to see that as Pðl^ i Þ tends to uniform distribution, the entropy maximisation algorithm becomes the ML algorithm k l^ kþ1 l^ i ¼ PM i j¼1 pij X yj pij PN ^ k i¼1 li pij ð14Þ For the rest of the paper we use the above simplified version of the proposed algorithm. It can be seen that this is slightly different from Green’s MAP algorithm [36]. We cannot give a general proof for the convergence of the above algorithm as this will depend on the form of the prior function P. 3 Modelling prior distribution The next step in the reconstruction process is the proper modelling of the prior distribution function, which incorporates pixel – pixel interaction in the reconstructed image. It is found that when the prior tends to a uniform distribution, the proposed algorithm tends to the ML algorithm. Only when the prior is chosen to be physically more accurate and meaningful, can one expect to obtain improved reconstruction over the ML approach. One such property of an image is nearest neighbour interaction between the pixels. The Markov random field (MRF) has been extensively used to incorporate such interactions in an image lattice [37, 38]. The MRF defined on a lattice is the collection of random variables lj ; corresponding to the site li such that Pðli jlj ; i 6¼ jÞ ¼ Pðli jlj ; j 2 Ni Þ  k  M  X  yj pij k k þ P l^ i ð1 þ log P y=l^ i PN ^ k j¼1 i¼1 li pij With this assumption, (12) simplifies to the following simple iterative procedure for updating the estimates: " # M X ^ ki y p l @Pðl Þ k kþ1 j ij i ^  k þP li l^ i ¼ PN ^ k k PM @li li ¼l^ i p P l^ i j¼1 i¼1 li pij j¼1 ij ð13Þ In other words, the entropy maximisation algorithm tends to the maximum likelihood algorithm when the prior distribution tends to a uniform distribution because a uniformly distributed prior gives equal probability for any estimate to be a solution. Further simplification of (12) is possible if we make the following assumption:   k k 1 þ log P y=l^ i  log P y=l^ i This assumption is justifiable on the basis that Pðy=l^ Þ is a probability distribution bounded between 0 and 1. In a PET system Pðy=l^ Þ has very small positive values, much smaller than 1. Hence, the assumption is valid whenever Pðy=l^ Þ  1: ð15Þ where Ni represents the subset of lattice sites in the neighborhood of the ith lattice site. The Hammersly – Clifford theorem [39] states that a random field over a lattice is a MRF if and only if its distribution function corresponds to a Gibbs distribution function and vice versa. So, a MRF model of image l characterised by a Gibbs distribution is given by PðlÞ ¼ 1 b1UðlÞ e Z ð16Þ where Z is the normalising constant for the distribution, b is the Gibbs hyper-parameter, and UðlÞ is the Gibbs energy. For homogeneous MRF, the energy function has the following form: XX Vðli ; lj Þ ð17Þ UðlÞ ¼ i j2Ni where Ni are the nearest neighbours of pixel i, and Vðli ; lj Þ is termed the potential at the site containing the subset Ni : Although a large number of potential functions have been suggested in the literature [37, 38, 40], we have used a simple quadratic potential for evaluation of the proposed algorithm. Consider the following potential function first suggested by Geman and Geman [41]: Vðli ; lj Þ ¼ ðli  lj Þ2 ð18Þ The potential Vðli ; lj Þ utilises the correlation between neighbouring pixels and forces them to be of similar value. This is justifiable because we expect neighbouring pixels not to vary drastically in an image field. It should be noted that the first term in (14) corresponds to the correction term due to the nearest neighbour contribution if the potential given by (18) is assumed, i.e.   I lki ¼ k l^ i @Pðli Þ  k PM @li P l^ i j¼1 pij k li ¼l^ i ð19Þ with PðlÞ given by (16). The second term of (14) is essentially the maximum likelihood contribution. For the simple potential defined by (18), the first correction term becomes k   l^ i 2 X X ^ k ^ k  li  lj ð20Þ I lki ¼ k PM b i j2N P l^ i i j¼1 pij As, the iteration increases, the numerator decreases since the value of l^ i slowly approaches the value of l^ j and in the limit k ! 1, Iðlki Þ ! 0: This is due to the smoothening   nature of the potential. The contribution of I lki is from P ^ ^ j2Ni ðli  lj Þ and b: For small values of b; the contribution is more, whereas for large b it is essentially ML. So, we need to search for an appropriate value of b: It has been found that for b > 2:5 104 ; the proposed algorithm converges. 4 Simulated experimental results Implementation of the proposed algorithm was performed on the simulated PET system shown in Fig. 1. The PET system consists of a ring detector system with 64 detectors and the object space is decomposed into 64 64 square pixels. The object space is a square region inscribed within the circle of detectors. To exploit the eight-fold symmetry we have assumed the object to be bounded by a circle inscribed in this square object space [15, 16]. An electron – positron annihilation event occurring inside a pixel is governed by a number of physical phenomena such as attenuation, scattering, absorption, and detector characteristics. All these physical processes have a bearing on the probability matrix. Each element of the probability matrix defines the probability of a photon being detected in the detector j after emanating from the object pixel i. For simplicity, we assume that the probability of emission in pixel i that is detected in a detector tube j is proportional to the angle yij seen by the centre of the pixel i into the detector tube j, i.e. pij ¼ yij =p (Fig. 2) [3, 4]. Before the reconstruction begins, the probability matrix Q ¼ ½ pij ; i ¼ 1; . . . ; N and j ¼ 1; . . . ; M is computed and stored. It is found that using the symmetry available in the system the data size of the Q matrix can be reduced considerably [16]. Fig. 1 head Typical PET system with detector ring surrounding the g-rays emanating from the pixels (shown as white squares) of the head phantom are subsequently detected Fig. 2 Pixel-detector probabilities ð p0ij Þ for simulated PET system computed from angle of view of the centre of pixel i into the tube j, ð pij ¼ yij =pÞ For simulating measurement data, a Monte Carlo procedure is used in which each emission is simulated as follows [3, 15]. First, a random pixel is chosen in the test phantom. The concentration of the radionuclei at the given pixel is assumed to be proportional to the emission density of the pixel. For each of the accepted emission points, a randomly oriented line (between 0 and p radians) is selected and the pair of detectors this line intersects is found. The random line corresponds to the direction in which the pair of annihilated photons travel and the pixel corresponds to the annihilation point in the phantom. The detectors which this line intersects are assumed to detect this annihilation event and the count in the corresponding detector tube is incremented. In this way all emissions are simulated and counted in the respective tubes, and used as the measurement data. A Shepp – Logan mathematical phantom with 100 000 emissions has been used for simulation studies. For iterative update of the estimates, the algorithm needs repeated passage from data space to image space and vice versa through the operator P. So, it is convenient to write (14) in an additive form i.e. kþ1 k k ¼ l^ þ Dðl^ Þ l^ ð21Þ Starting with a small positive initial value of the estimate, k (21) updates the old estimate l^ ; by adding the correction k kþ1 factor Dðl^ Þ to produce the new estimate l^ : The mean squared error (MSE) between the original Shepp – Logan phantom and the reconstructed phantom is considered to evaluate the performance of the proposed algorithm. The plot of MSE in Fig. 3 illustrates the general nature of the proposed algorithm for various b values. It is clear from the plot that only for restricted values of b does the algorithm produce superior quality reconstruction. We term this interval of b; the optimum interval. The image quality is found to be degraded outside this optimal interval, irrespective of the number of iterations performed. The reason for this is the smoothing behaviour of the potential. This potential increasingly penalises the difference between neighbouring pixels, thus forcing them to become the same. In terms of entropy, it is clear that with increasing b the information content decreases, thereby resulting in Fig. 3 Plot of mean square error against hyper-parameter b shown after different numbers of iterations a 40 iterations b 100 iterations saturation of the mean square error (MSE) as shown in Fig. 3. On the other hand if b is in the optimal interval, then there is an increase in the information content gained at the cost of the local correlation invoked in the reconstruction procedure via prior knowledge. In the limit b ! large value, the maximum entropy algorithm approaches the ML algorithm. Below the optimal interval, the estimate becomes unstable due to the dominating nature of the prior. This basically removes the contribution of measurement data in the reconstruction algorithm. These can be visually verified from the reconstruction shown in Fig. 6. Once b is chosen in the optimal interval, estimates will produce better reconstruction as the iterations continue. This is evident from the lower MSE of the reconstructed images as shown in Fig. 3 for 40 and 100 iterations. In Fig. 4, plots of the MSE against the number of iterations are shown with varying b values from 1:5 104 to 2:5 107 : The MSE is found to decrease with increasing iterations for b ¼ 2:5 104 to 2:5 105 (see Fig. 4). On the other hand, MSE increases with increasing iterations for b  2:5 105 (see Fig. 4k and Fig. 4l). The b value indicates the amount of prior to be introduced into the reconstruction procedure. When b is increased beyond a certain large value, then the reconstruction is effectively due to the likelihood contribution; the general nature of which is depicted in Fig. 4k and Fig. 4l. It is well known that the reconstruction improves during the earlier iterations of ML estimation, but becomes noisier as the iterations continued [8]. Fig. 4 Plots of mean square error against iterations a b ¼ 1:5 b b ¼ 2:5 c b ¼ 3:5 d b ¼ 4:5 e b ¼ 5:5 f b ¼ 6:5 g b ¼ 7:5 h b ¼ 1:0 i b ¼ 2:0 j b ¼ 2:5 k b ¼ 2:5 l b ¼ 2:5 104 104 104 104 104 104 104 105 105 105 106 107 Fig. 5 Original 64 64 Shepp – Logan mathematical phantom Fig. 7 Mean squared error against number of iterations showing the change in MSE for ML, MAP with different hyper-parameter values, and the proposed algorithm Fig. 6 Visual representation of image reconstruction after 100 iterations of entropy maximisation algorithm for varying b values a b ¼ 1:5 b b ¼ 2:5 c b ¼ 3:5 d b ¼ 4:5 e b ¼ 5:5 f b ¼ 6:5 g b ¼ 7:5 h b ¼ 1:0 i b ¼ 2:0 j b ¼ 2:5 k b ¼ 2:5 l b ¼ 2:5 104 104 104 104 104 104 104 105 105 105 106 107 It is evident from the proposed algorithm that with an increase in the b value, the effect of the prior is nullified by the data and the maximum entropy estimation becomes equivalent to the EM estimation. This can be further seen in the reconstruction shown in Figs. 5 and 6. With increasing b values, the image makes a transition from smooth reconstruction to non-smooth. This non-smoothness in the reconstruction is an indicator of noise. It is clear from the reconstruction shown in Fig. 6 that noise artefacts are much more severe in Fig. 6k and Fig. 6l than Figs. 6a– 6j, which are the reconstructions produced by b values in the optimal interval. Good reconstruction is found to occur for 2:5 104 b 2:0 105 : It is found that the entropy maximisation estimate becomes unstable for b < 104 : The reconstructed images become increasingly noisy for large b values. This can be visually verified in the reconstructions shown in Fig. 6. We also compare the results of the proposed algorithm with the well known maximum a posteriori (MAP) estimation [36, 40]. The potential used for the MAP estimation algorithm [36] and the proposed algorithm is the same (i.e. (18)) so as to compare the performance. Simulations were carried out for Gibbs hyperparameter b ranging from 103 to 107 . For brevity, we have shown MSE against iterations plots for some of the key b values in Fig. 7. Comparing MSE values for the proposed algorithm and the MAP algorithm (Fig. 7), it is clear that the proposed algorithm performs slightly better than the MAP algorithm. It should be noted that more realistic potentials are capable of producing better estimates but require additional parameter estimation [37, 38]. The second potential used for reconstruction is similar to the Huber function: VH ¼ ( 1 2m ðli  lj Þ 2 ; mjli  lj j; if jli  lj j otherwise m ð22Þ where m is an additional parameter to be chosen. Simulations are carried out using (22) as the potential and the results are shown in Fig. 8. The first row and second row correspond to m ¼ 1; 10; 25 and 50 for MAP and the proposed algorithm, respectively. It is clear that quality decreases as m increases because for small pixel differences, the penalty is quadratic while the penalty is linear when the differences are large. Figure 9 shows images reconstructed using ML, MAP and the proposed algorithm with potential V (18), after 80 iterations. For visual inspection and quality of reconstruction the original test phantom is also shown in Fig. 9a. The ML reconstructed image is shown in Fig. 9b. The first and second row show reconstructed images for different Gibbs hyper-parameter values using MAP and the proposed algorithm, respectively. It can be seen that the images reconstructed using the proposed algorithm look sharper than the images reconstructed using the MAP algorithm. Fig. 8 Reconstructions obtained using MAP algorithm and the proposed algorithm a b c d e f g h MAP, m ¼ 1 MAP, m ¼ 10 MAP, m ¼ 25 MAP, m ¼ 50 Proposed algorithm, m ¼ 1 Proposed algorithm, m ¼ 10 Proposed algorithm, m ¼ 25 Proposed algorithm, m ¼ 50 Fig. 9 Images reconstructed after 80 iterations a b c d e f g h 5 Original phantom Reconstructed by ML algorithm Reconstructed by MAP algorithm ðb ¼ 104 Þ Reconstructed by MAP algorithm ðb ¼ 105 Þ Reconstructed by MAP algorithm ðb ¼ 107 Þ Reconstructed by proposed algorithm ðb ¼ 104 Þ Reconstructed by proposed algorithm ðb ¼ 105 Þ Reconstructed by proposed algorithm ðb ¼ 107 Þ Conclusions A new algorithm is proposed for the PET system, which is capable of good reconstruction. This is based on the IEE Proc.-Vis. Image Signal Process., Vol. 151, No. 5, October 2004 maximisation of conditional entropy as a cost function rather than the likelihood or posterior function. The proposed method ensures positivity of the reconstructed image pixels since entropy is not defined for negative pixels. 351 Furthermore, the proposed methodology enables the inclusion of prior distribution knowledge of the object during the image reconstruction process. The Gibbs distribution function is used for modelling the prior. The attractive feature of a Gibbs distribution prior is its ability to incorporate local correlation between neighbouring pixels in an image lattice. Simulated experimental results have shown much improved reconstruction for a broad range of Gibbs hyper-parameter values. It is shown that entropy is a better candidate as a cost function than the likelihood or posterior function. 6 Acknowledgments The authors wish to thank the anonymous reviewers, who helped in improving the quality of this manuscript. The first author would like to thank Mr Hemanth, research scholar at the Indian Institute of Science for helping during the simulation experiment. Also, he wishes thanks to the Council of Scientific and Industrial research (CSIR), Government of India, for providing a Junior Research fellowship. He acknowledges that this work is partially supported by CSIR, New Delhi, India. 7 References 1 Frieden, B.R.: ‘Restoring with maximum likelihood and maximum entropy’, J. Opt. Soc. Am., 1972, 62, pp. 511–518 2 Frieden, B.R., and Wells, D.C.: ‘Restoring with maximum entropy: Poisson sources and backgrounds’, J. Opt. Soc. Am., 1978, 68, (1), pp. 93–103 3 Vardi, Y., Shepp, L.A., and Kaufmann, L.: ‘A statistical model for positron emission tomography’, J. Am. Stat. Assoc., 1985, 80, pp. 8– 37 4 Shepp, L.A., and Vardi, Y.: ‘Maximum likelihood estimation for emission tomography’, IEEE Trans. Med. Imaging, 1982, 1, pp. 113–121 5 Craig, A.D., Reiman, E.M., Evans, A., and Bushnell, M.C.: ‘Functional imaging of an illusion of pain’, Nature, 1996, 384, pp. 217 –218 6 Decety, J., Perani, D., Jeannerod, M., Bettinardi, V., Tadary, B., Woods, R., Mazziotta, J.C., and Fazio, F.: ‘Mapping motor representations with positron emission tomography’, Nature, 1994, 371, pp. 600–602 7 Chernoboy, E.S., Chen, C.J., Miller, M.I., Miller, T.R., and Snyder, D.L.: ‘An evaluation of maximum likelihood reconstruction for SPECT’, IEEE Trans. Med. Imaging, 1990, 9, pp. 99–110 8 Veclerov, E., and Llacer, J.: ‘Stopping rule for MLE algorithm based on statistical hypothesis testing’, IEEE Trans. Med. Imaging, 1987, 6, pp. 313–319 9 Chen, C.M., and Lee, S.Y.: ‘A parallel implementation of 3-D CT image reconstruction on hypercube multi processor’, IEEE Trans. Nucl. Sci., 1990, 37, (3), pp. 1333–1346 10 Chen, C.M., and Lee, S.Y.: ‘Parallelization of the EM algorithm for 3-D PET image reconstruction’, IEEE Trans. Med. Imaging, 1991, 10, (4), pp. 513–522 11 Rajan, K., Patnaik, L.M., and Ramakrishna, J.: ‘High speed computation of the EM algorithm for PET image reconstruction’, IEEE Trans. Nucl. Sci., 1994, 41, (5), pp. 1–5 12 Ranganath, M.V., Dhawan, A.P., and Mullani, N.: ‘A multigrid expectation maximization algorithm for positron emission tomography’, IEEE Trans. Med. Imaging, 1988, 7, pp. 273–278 13 Tanaka, E.: ‘A fast reconstruction algorithm for stationary positron emission tomography based on a modifie EM algorithm’, IEEE Trans. Med. Imaging, 1987, 3, (1), pp. 98–105 14 Lange, K., Bahn, M., and Little, R.: ‘A theoretical study of some maximum likelihood algorithm for emission and transmission tomography’, IEEE Trans. Med. Imaging, 1987, 6, (2), pp. 106– 114 15 Rajeevan, N., Rajgopal, K., and Krishna, G.: ‘Vector-extrapolated fast maximum likelihood estimation algorithms for emission tomography’, IEEE Trans. Med. Imaging, 1992, 11, (1), pp. 9 –20 16 Kaufmann, L.: ‘Implementing and accelerating the EM-algorithm for positron emission tomography’, IEEE Trans. Med. Imaging, 1987, 6, pp. 37– 51 17 Hudson, H.M., and Larkin, R.S.: ‘Accelerated image reconstruction using ordered subsets of projection data’, IEEE Trans. Med. Imaging, 1994, 13, (4), pp. 601–609 18 Hudson, H.M., Hutton, B.F., and Larkin, R.: ‘Accelerated EM reconstruction using ordered subsets’, J. Nucl. Med., 1992, 33, p. 960 19 Browne, J., and De Pierro, A.: ‘A row-action alternative to the EM algorithm for maximizing likelihoods in emission tomography’, IEEE Trans. Med. Imaging, 1996, 15, pp. 687–699 20 Byrne, C.L.: ‘Accelerating EMML algorithm and related iterative algorithm by rescaled block iterative methods’, IEEE Trans. Image Process., 1998, 7, pp. 100–109 21 Lalush, D.S., Frey, E.C., and Tsui, B.M.W.: ‘Fast maximum entropy approximation in SPECT using RBI-MAP algorithm’, IEEE Trans. Med. Imaging, 2000, 19, pp. 286–294 22 Snyder, D.L., and Miller, M.I.: ‘The use of sieves to stabilize images produced with the EM-algorithm for emission tomography’, IEEE Trans. Nucl. Sci., 1985, 32, (5), pp. 3864–3872 23 Hebert, T., and Leahy, R.: ‘A generalized EM algorithm for 3-D Bayesian reconstruction from Poisson data using Gibbs priors’, IEEE Trans. Med. Imaging, 1989, 8, (2), pp. 194– 202 24 Levitan, E., and Herman, G.T.: ‘A maximum a posteriori probability expectation maximization algorithm for image reconstruction in emission tomography’, IEEE Trans. Med. Imaging, 1987, 6, (3), pp. 185– 192 25 Mondal, P.P., and Rajan, K.: ‘Entropy maximization algorithm for positron emission tomography’. Presented at 9th Int. Conf. on Neural Information Processing, Singapore, 18–22 Nov. 2002 26 Mondal, P.P., and Rajan, K.: ‘Conditional entropy maximization algorithm for PET’. Presented at IEEE Int. Conf. on Systems, Man and Cybernetics, Washington DC, USA, 5 –8 October 2003 27 Skilling, J., Strong, A.W., and Bennett, K.: ‘Use of maximum entropy method in gamma-ray astronomy’, Mon. Not. R. Astron. Soc., 1979, 187, pp. 145–152 28 Skilling, J., and Bryan, R.K.: ‘Maximum entropy image reconstruction: general algorithm’, Mon. Not. R. Astron. Soc., 1984, 211, pp. 111–124 29 Gull, S.F., and Daniel, G.J.: ‘Image reconstruction from incomplete and noisy data’, Nature, 1978, 272, pp. 686–690 30 Minerbo, G.: ‘MENT: A maximum entropy algorithm for reconstructing a source from projection data’, Comput. Graph. Image Process., 1979, 10, pp. 48– 68 31 Willingale, R.: ‘Use of maximum entropy method in X-ray astronomy’, Mon. Not. R. Astron. Soc., 1981, 194, pp. 359–364 32 Byrne, C.L.: ‘Likelihood maximization for list-mode emission tomographic image reconstruction’, IEEE Trans. Med. Imaging, 2001, 20, (10), pp. 1084– 1092 33 Byrne, C.L.: ‘Iterative image reconstruction algorithms based on crossentropy minimization’, IEEE Trans. Image Process., 1993, 2, (1), pp. 96– 103 34 Titterington, D.: ‘On the iterative image space reconstruction algorithm for ECT’, IEEE Trans. Med. Imaging, 1987, 6, pp. 52–56 35 Zangwill, W.I.: ‘Nonlinear programming: a unified approach’ (Prentice Hall, Englewood Cliffs, NJ, USA) 36 Green, P.J.: ‘Bayesian reconstruction from emission tomography data using a modified EM algorithm’, IEEE Trans. Med. Imaging, 1990, 9, (1), pp. 84–93 37 Zhou, Z., Leahy, R.M., and Qi, J.: ‘Approximate maximum likelihood hyperparameter estimation for Gibbs prior’, IEEE Trans. Image Process., 1997, 6, (6), pp. 844–861 38 Hebert, T., and Leahy, R.: ‘Statistic based MAP image reconstruction from Poisson data using Gibbs Priors’, IEEE Trans. Signal Process., 1992, 40, (2), pp. 2290–2303 39 Besag, J.: ‘Spatial interaction and the statistical analysis of lattice systems’, J R. Stat. Soc. B, 1974, 36, pp. 192– 236 40 Nuyts, J., Bequ, D., Dupont, P., and Mortelmans, L.: ‘A concave prior penalizing relative differences for maximum-a-posteriori reconstruction in emission tomography’, IEEE Trans. Nucl. Sci., 2002, 49, (1), pp. 56– 60 41 Geman, S., and Geman, D.: ‘Stochastic relaxation, Gibbs distribution and the Bayesian restoration of images’, IEEE Trans. Pattern Anal. Mach. Intell., 1984, 6, pp. 721–741