A Penalized-Likelihood Image Reconstruction Algorithm for Positron Emission Tomography Exploiting Root Image Size

Iterative image reconstruction methods are considered better as compared to the analytical reconstruction methods in terms of their noise characteristics and quantification ability. Penalized-Likelihood Expectation Maximization (PLEM) image reconstruction methods are able to incorporate prior information about the object being imaged and have flexibility to include various prior functions which are based on different image descriptions. Median Root Priors intrinsically take into account the salient image features, such as edges, which becomes smooth owing to quadratic priors. Generally, a 3*3 pixels neighborhood support or root image size is used to evaluate the median. We evaluate different root image sizes to observe their effect on the final reconstructed image. Our results show that at higher parameter values, root image size has pronounced effect on different image quality parameters evaluated such as reconstructed image bias as compared to the phantom image, contrast and resolution in the reconstructed object. Our results show that for the small-sized objects, small root image is better whereas for objects of diameter more than two to three times of the resolution of the reconstructed object, larger root image size is preferable in terms of reconstruction speed and image quality. Keywords—Penalized-Likelihood expectation maximization; median root priors; maximum-likelihood expectation maximization; full-width-at-half-maximum


INTRODUCTION
Positron emission tomography (PET) and Single photon emission tomography (SPECT) are used to image human body functions non-invasively.Data obtained from these scanners, is reconstructed by analytical or iterative image reconstruction methods to estimate emission object's activity distribution.Analytical image reconstruction methods use line integral model and simply ignore any underlying noise distribution [1], [2].It is also not so simple to incorporate emission or detection physics model into the reconstruction problem and these methods do not take into account the non-negativity condition due to Fourier transforms employed.On the other hand, statistical iterative image reconstruction methods can include modeling for emission and detection processes into the system matrix and automatically fulfill non-negativity constraint [3].Owing to this ability, statistical iterative image reconstruction methods are claimed to be superior to analytical image reconstruction methods with respect to noise characteristics and quantification ability [3], [4].
A very popular and basic iterative image reconstruction method is known as maximum-likelihood expectation maximization (MLEM) method [5].Various techniques, such as post reconstruction smoothing or method of Sieves, have been used to reduce this resultant image noise [7]- [9].We adopt the same definition of the prior derivative in this paper.The equally popular class of priors based on median root priors (MRPs) is used.MRPs are designed on the basic image description that images are locally monotonic.In MRP's algorithm, neighborhood support of the prior function is termed as root image and most frequently a 3*3 pixels neighborhood is used as the root image.Median value evaluated on this neighborhood is used to penalize image pixels.We evaluate the impact of variation of the root image size on the final reconstructed image by proposing a penalized-likelihood image reconstruction algorithm which exploits root prior knowledge.
To the best of our knowledge, there is not much work done to evaluate the impact of root image size on final reconstructed image characteristics [15].It is observed that the correlation between image pixels decreases as the distance between two pixels is increased.However, the value of the median evaluated on different root image sizes may vary.This work focuses on analyzing the impact of this variation on the properties of the final reconstructed image.
The rest of this paper is organized as follows.In Section II, we describe the methodological background.In Section III, we outline the proposed algorithm.Experimental results are shown in Section IV.In Section V, the discussion on the experimental results is presented.Finally, the paper is concluded in Section VI www.ijacsa.thesai.orgII.METHODOLOGICAL BACKGROUND MLEM is an iterative optimization method used to reconstruct images from scanner data known as sinogram.Mathematically, this method attempts to maximize an objective function based on the statistical average of the logarithm of the likelihood function as follows:

̂
(1) where the objective function L(X,Y) is defined as the logarithm of the likelihood function of the emission data Y given object X and the statistical expectation of Y is considered to be the mean of the independent Poisson emission model as given below.A ij is the system matrix element and characterizes the probability of an event being detected at the i-th bin and emitted from the j-th pixel point.

̅̅̅
However, with increasing iteration number this method is known to produce noisier images as image reconstruction is an ill-posed inverse problem and this technique attempts to fit the solution image to the data and does not consider any priori information about the object being imaged [6].
Various techniques such as post reconstruction smoothing or method of Sieves, have been used to reduce this resultant image noise [7]- [9].However, penalized-likelihood image reconstruction methods are more favorable due to their flexibility to incorporate various penalizing schemes.A penalty function is added to the likelihood function based on some priori image description or knowledge and logarithm of this modified objective function is maximized, instead, as described in Eq. (II).These functions induce additional constraints to reduce the solution image set [10], [11].
In relationship given below, R(X) is the prior term and depends on some priori image information or image description and is the hyper parameter which controls prior influence on the final image.Generally, quadratic prior functions are used due to their implementation simplicity [3], [12].However, these priors produce overly smoothed salient image features due to their penalizing scheme to modify pixel values on the basis of the differences in their values [3], [7], [11], [13].These priors are designed with the basic concept of image description that images are locally smooth.Hence, quadratic priors attract all image pixels in a local neighborhood which are close to each other in values [1].Following relationship describes mathematical form of the quadratic priors: In the above relationship, are the weights assigned to pixels in a smaller neighborhood and V(x) is energy function as defined in Gibbs distribution [14].Another equally popular class of priors based on median root priors (MRPs).MRPs are designed on the basic image description that images are locally monotonic.Interestingly, this includes the definition of local smoothness.Hence, we may observe some level of smoothing in the resultant image.MRPs also have their ability to automatically preserve edges inside the image without any need of an additional tuning parameter which is a disadvantage of other non-quadratic edge-preserving priors [1], [15], [16].MRPs can be described as follows:


Where M j is the median in local neighborhood of jth pixel and N b is the number of pixels in the neighborhood of the relevant disk.Pixel neighborhoods are generally considered as Markov Random Fields (MRFs) and described according to Gibbs distribution [14], [17].These prior functions are incorporated into the reconstruction algorithm according to Greens very popular One-Step-Late (OSL) algorithm.According to this algorithm, we need first derivative of the objective function, including log-likelihood and logarithm of the prior function, to obtain the current image estimate.However, while the current image is being estimated, it is not possible to evaluate first derivative at the current image estimate.Derivative of the prior function is evaluated for image estimated at the previous iteration, hence the name OSL algorithm [18].
OSL algorithm needs first derivative of the prior function and, unfortunately, it is not possible to directly evaluate the gradient of the median function as dependence of median on local neighborhood is nonlinear.We have to resort to some empirical derivation of the gradient of median based priors.MRP has been defined as Gaussian distribution of the prior function, in the form of Gibbs priors [15], [17].For the sake of the evaluation of the derivative of MRP, authors have assumed an empirical dependence of the local median on the image values as presented in Fig. 1.They have assumed local median to be constant within the functioning limits of the MRP function.Under this assumption, it becomes much easier to find an empirical derivation of the prior function.We adopt

Median
Pixel value f j < M f j > M www.ijacsa.thesai.org the same definition of the prior derivative in this paper.In MRP algorithm, neighborhood support of the prior function is termed as root image and most frequently, a 3*3 pixels neighborhood is used as the root image.

III. PROPOSED ALGORITHM
We use a 128*128 pixels image grid with a circular emission object comprised of a main background disk, having four small disks inside it placed at the same distance from center of the background disk as shown in Fig. 2 (left).The simulated phantom image was based on the simulation scheme shown in Fig. 2 (bottom).This figure shows a square grid encircling reconstruction circle inside its boundaries.A random point was generated inside the grid along with a randomly generated angle to simulate a random event.This simulated event, based on a random point and angle pair, was assumed to describe path of a randomly emitted photon pair and detected in a particular line-of-response (LOR).Index of this LOR was calculated by the following (6) and describes perpendicular distance of the assumed path of emitted pair of photons.(6) In this expression, t represents perpendicular distance from the origin to the LOR, defined by the line of flight of the pair of photons emitted at point (x; y) and traveling at an angle with the x-axis.To reconstruct images for our analysis, we have used One-Step-Late (OSL) algorithm by P. Green [18].That algorithm can be given by the following relationship.
Where C ML k is the maximum-likelihood correction factor and C PL k is the priors correction factor both evaluated at kth iteration.As mentioned by Green [18], prior function needs to be a continuous function having continuous derivatives, unluckily, median priors do not have its derivatives defined.However, Alenius [15] has developed some heuristic approach to visualize dependence of median on the local neighborhood as shown in Fig. 1.According to this approach, median can be considered constant within the operational range of the MRP.
Our work deals with the effect of root image size on the images reconstructed by the PLEM image reconstruction methods including MRPs, hence we have used varying size of root images.Generally, a 3*3 pixels or 5*5 pixels root image is used to reconstruct images with PLEM including MRPs.We used 3*3, 5*5, 7*7 and 9*9 pixels root images to analyze characteristics of the reconstructed images.It is important to note that at the boundaries of the image grid, we need zero padding in order to apply MRPs.We have zero-padded the images as shown in Fig. 3.
Initially, for our analysis, to see if size of the root image does produce any change in the reconstructed image, we evaluate reconstructed image bias (RIB) relative to the true phantom image based on the following equation: ∑ Where X k and Ph k are the relative pixels in the reconstructed image and the phantom image respectively and N is the total number of pixels in the image.We analyze the behavior of MRPs relative to the size of their root image in terms of image contrast in a local neighborhood.We select 4 4 pixels ROIs inside each small disk, at higher activity as compared to the background, to obtain positive contrast www.ijacsa.thesai.orgvalues.Contrast is evaluated using the following formula: (11) In this relation, X e is the estimated mean value in an ROI in the hot disk and X b is the estimated mean value in the background disk of the same size ROI.Finally, we study the effect of root image size on the resolution of the reconstructed image.We generate two realizations of the sinogram to evaluate reconstructed resolution, one with an impulse added to the centre pixel of each disk and the other one without an added impulse.These sinograms are then reconstructed to obtain two images f(X +X) and f(X).Resolution in Full-Width-at-Half-Maximum (FWHM) is calculated from the differential image evaluated using the following formula: (12) Where X is a very small quantity and logically tends to zero for the differential image.A Matlab Function fwhm2.m,developed by J. A. Fessler is used to calculate resolution in pixels at the centre of each small disk.The proposed algorithm is outlined in Algorithm 1.

IV. EXPERIMENTAL RESULTS
We evaluate reconstructed image bias relative to the true phantom image, to observe any effect of root image size on reconstructed images, and presented our results in Fig. 4. We use (10) to calculate image bias and this figure displays bias values for different beta values and root image sizes.We observe that image bias reduces with decreasing value of beta parameter value, reconstructed image moves more and more towards its MLEM solution, hence reducing effect of the prior and, consequently, of the root image size.It is clear that at low beta values, there is not much difference in bias values for different root image sizes.However, as beta value increases, and gets near 1, image bias increases and clearly is affected by different root image sizes in penalized-likelihood image reconstruction method using MRPs.We also analyze the effect of root-image-size in MRPs on the reconstructed contrast values inside each small disk, using (4), and present our results in Fig. 5.This image clarifies that the reconstructed contrast varies with different root image size at higher beta values for smaller objects.This contrast values becomes almost same for different root image sizes for the largest disk size having highest activity level.However, for lower beta values root image size does not have much effect on the reconstructed contrast in different objects and for different root image sizes.Fig. 6 presents results for reconstructed resolution, in FWHM (Pixels), at the centre of each small disk for two different beta parameter values.With higher beta value, reconstructed resolution is much different for different root image sizes.However, for lower beta value where priors influence is reduced, root image size does not seem to have strong influence.This effect is similar to that of different root image sizes over the image bias or reconstructed image contrast values.
Average reconstructed resolution also seems to be higher (lower FWHM values) in case of lower beta value as compared to the one with higher beta values which may be considered as higher smoothing at high beta values.This may not, strictly, be the same kind of smoothing as produced by quadratic priors [2].It should also be noted that DISK2 and DISK4 are on the vertical axis whereas DISK1 and DISK3 are on the horizontal axis.This indicates an asymmetrical reconstructed resolution response.

V. DISCUSSION
We have analyzed effects of root image size on the final reconstructed image using penalized-likelihood image reconstruction methods including median root based priors.We have evaluated this effect in terms of reconstructed image bias as compared to the true phantom image, reconstructed contrast in high activity regions and resolution in FWHM (in pixels) at certain locations inside the image.
Our results show that root image size has pronounced effect on the bias values for higher beta parameter values.Higher beta value, in case of MRPs, means higher priors influence on the final reconstructed image which means that it is attracted more towards local medians.At higher beta values, shape of the prior is more sharply peaked and locally nonmonotonic areas are dragged towards local median with stronger force.Moreover, with higher beta values and bigger root image size, response of MRPs is farthest from that of any local smoothing prior (quadratic priors), hence reduced smoothing will be observed.That could be the reason that we observe higher bias values with bigger root image size at higher beta values.This can also be attributed to less correlation between distant pixels.
Additionally for contrast values, no significant difference is observed at lower beta values with different root image size.However, at higher beta values, contrast is much reduced for bigger root image size for small object size containing low activity value as compared to the smaller neighborhood in the same object.This may also mark presence of partial volume error because this effect is reversed in the bigger size disks containing higher activity values where contrast is higher for bigger root image size.For higher values of beta parameter, resolution in FWHM (Pixels) varies markedly for different size of root image and even the trend varies in horizontal and vertical directions.DISK1 and DISK 3 are on the horizontal axis where the reconstructed resolution is worse (or higher in FWHM) as compared to the DISK 2 and DISK 4 that are lying on the vertical axis.This effect is totally strange for the smallest root image size and could represent a mixture of root image partial volume error combined due to the disk size smaller than 2 to 3 times of the resolution (FWHM in Pixels).However, at lower beta values (near to zero) difference between different root image size disappeared, though the effect of object size and activity level still persists in the form of varying resolution at the center of these disks.

Fig. 1 .
Fig. 1.An empirical dependence of median in a local neighborhood against those pixel values in the same neighborhood.

Fig. 2 .
Fig. 2. Data simulation scheme (bottom) simulated Phantom image (left) and simulated sinogram data (right) which was used to reconstruct images for our analysis evaluated on this neighborhood is used to penalize image pixels.

Fig. 3 .
Fig. 3. Zero padding scheme for the application of median.

Fig. 4 .
Fig. 4. Reconstructed image bias relative to the true phantom image using MRPs for different root-image sizes and beta parameter values.

Fig. 5 .
Fig. 5. Contrast recovery results with MRPs for different size of root images.DISK1 to DISK4 present recovered contrast for different activity values and different size of the object for two different values of beta hyper parameter values.At higher beta values, for the smallest DISK1 with least activity, contrast varies for different root-image-sizes, whereas for largest DISK4 with highest activity value, contrast is almost same.

Fig. 6 .
Fig. 6.Resolution in FWHM (Pixels) with MRPs for different size of root images at the center of each small disk.DISK1 to DISK4 present reconstructed resolution for different activity values and different size of the object for two different values of beta hyper parameter values.At higher beta values, for the smallest DISK1 with least activity, contrast varies for different root-image-sizes, whereas for largest DISK4 with highest activity value, contrast is almost same.