2. Methodologies
2.1 Reconstruction with Random Correction for MicroPET Images
A new approach is proposed to reconstruct microPET images with random correction by the joint Poisson model of prompt and delay sinograms. We will assume that the data in prompt and delay sinograms follow two independent Poisson distributions that are labeled as (1) and (2).
)),
target image with size B = 128x128 in which b=1,2,…,B, and d indicates the dth LOR with total numbers D = 96x84 in which d=1,2,…,D. The term n*p(d) is the number of coincidental events in the prompt sinogram at the dth projection line of response (LOR), which is formed by two detectors with the Poisson mean,λ*(d); n*d(d) is the
number of random coincidental events in the delay sinogram with the Poisson
meanλ*r(d); P(b,d) is the system probability matrix from the bth pixel to the dth LOR.
Parameters of λt(b) and λ*r(d)are unknown and must be estimated. Parameter )
t(b
λ represents the intensities of true coincidental events. Appendix (A3) lists the log-likelihood of observed data in the prompt and delay sinograms. Since the maximum likelihood estimate (MLE) is difficult to find by maximizing Eqs. (1) and (2) numerically, the EM algorithm [30-31] can be utilized (as the details in Appendix).
Equations (3) and (4) are the ith iteration steps of the PDEM.
where i=1,2,……I, is the number of iterations.
The MLEM algorithm of joint the prompt and delay sinograms is described as follows and such a scheme is called PDEM reconstruction.
Algorithm for PDEM reconstruction:
Step 1: Set the initial parameters (i=0) using the FBP, the method of moments estimate (MME) or alternative approaches.
Step 2: The iteration number is increased by 1 (that is, i is replaced by i+1). Update
the parameters by applying Eqs. (3) and (4).
Step 3: If (that is, the difference of
log-likelihoods among new and old estimates is smaller than for a tolerance level), then the iteration is terminated; otherwise, go to Step 2.
tolerance d
b l
d b
lin(λit( ),λr*i( ))− in(λit−1( ),λ*ri−1( ))<
This method preserves Poisson properties and update estimates iteratively. In this study, P(b,d) was computed from LORs and the locations of pixels based on the geometric characteristics of the microPET R4, including number of detectors, image size, field of view (FOV), ring diameter, and number of angular views. The matrix size of one slice is 96×84. There are 96 angular views and 84 LORs for each angular view during image scanning. Furthermore, each P(b,d) can be identified from its detector pairs of LOR and image pixel location. Therefore, the PDEM reconstructs the sinogram after the step that sinograms are rebinned by the FORE approach in the microPET system.
2.2 Segmentation of 3D MicroPET Images
This section briefly introduces the segmentation algorithm of 3D microPET image.
Several pre-processing steps are used before the GMM algorithm is applied to segment 3D microPET image. The number of clusters to be used in GMM and the initialized values of each cluster are needed in applying the GMM method. Hence, the
nonparametric method of KDE is employed to estimate the density curve of image intensities. The numbers of cluster (k) can be determined by searching the number of local maximum points from the estimated density curve. In addition, the initialized values of parameters for the normal distributions in k clusters can be determined from the data consequently.
2.2.1 Determination of the Cluster Number for GMM
The kernel density estimation (KDE) [41] will be applied to determine the number of clusters used in GMM. From the estimated density curve, those high and low peaks could be used to estimate the means and standard deviations of Gaussian distributions in clusters. It is based on the empirical rule that the range of μ ±3σ covers most of observations from a normal distribution. Hence, the initialized values of GMM are determined by applying this empirical rule (as illustrated in Fig. 2.2.1-1).
This approach can automatically decide the cluster number and starting values of parameters in the EM algorithm of GMM. This is a simple and computation efficiently method. The details are reported in the following algorithm.
Fig. 2.2.1-1. An example shows the procedure to decide the number of clusters and initialized values in the EM algorithm of GMM. There are two peaks from the estimated densities of . Hence, the number of clusters in GMM will be two and initialized values will be computed from the above equations.
) ˆ y( f
Algorithm of KDE for Determining the Number of Clusters
A Gaussian kernel function is used to estimate the density of data as in Eq. (5), , kernel to estimate a probability density function, n is the sample size, and j = 1, 2, …, 128.
Step 1: Input dataX ={x1,x1,L,xn}.
Step 2: Find 128 grid points using equally spaced as Eq. (6):
. Step 3: Calculate the data-driven bandwidth for KDE as Eq. (7):
Step 5: Determine the number of clusters by counting the number of local maximum.
The high and low peaks of estimated density can be used to estimate the means and
standard deviations of Gaussian distributions in clusters like the illustration in Fig.
2.2.1-1.
2.2.2 Gaussian Mixture Model (GMM)
For the study of a set of 3D images from microPET, Ixyz represents the intensity at the xyzth voxel. The range of the position xy in one slice has the size of 128x128 and the number of slices z has the size of 10 in this study. Since the variance scale in different segments of 3D images may be different, we consider Gaussian mixture model with different parameters of means and variances in various clusters [28-29].
We will suppose that the image data Ixyz follow a mixture of K distributions with the
mixing probability πk such that
That is, the probability density function of Ixyz is , parameter θk. All parameters are collected to form a parameter vector
).
Φ The MLE of every parameter is usually difficult to obtain directly by numerical methods. Alternatively, the EM algorithm can be applied to find the MLE iteratively [28-29]. Firstly, we introduce an index function as follows,
⎪⎩ the complete data for applying the EM algorithm. Given Cxyz, the conditional density of Ixyz becomes The EM algorithm finds MLE of parameters with the constraint of Eq. (8) by the following iterative procedure:
The Gaussian distribution for the kth cluster is listed below,
.
The resulting EM algorithm for GMM is described as follows.
EM Algorithm for GMM
Step 1: Set the initial parameters Φ (old).
Step 2: Update the parameters by using Eqs. (12) - (15).
Step 3: If lin(Φ(new) ) –lin (Φ(old)) < tolerance (that is, the difference of log-likelihoods among new and old estimates is smaller than for a tolerance level), then the iteration stops. Otherwise, go to Step 2 with the old values of parameters replaced by the new values.
2.2.3 Slice Normalization
The intensity heterogeneity in slices of microPET images comes from the differences in emission activities and reconstruction processes. The heterogeneity between slices will affect the accuracy of the segmentation of 3D images. Hence, slice normalization is needed before 3D segmentation to adjust image levels between slices.
We will use the slice normalization defined in Eq. (17),
z z
z new xyz
xyz
Max Min
Min a I
I −
= −
, (17)where is the new intensity after slice normalization, is the original intensity at xyz
new
Ixyz Ixyz
th voxel, Minz is the minimum value of the zth slice, z is the
maximum value of the z
Max
th slice, and a is a nonzero constant for the range of adjusted
image level.
2.3 Segmentation of Microarray Images
This section discusses three segmentation algorithms for spotted microarrays.
Empirical evaluation of the performance for segmentation methods will use different criteria based on target ratios of spike genes, duplicate spots in a microarray image, and swapped microarray images Cy3 and Cy5 dyes. In this study, the aim is to segment pixels into foreground or background. Hence, there are two clusters used for the segmentation of the pixel intensities in a spot. Besides the GMM with EM algorithm, we will propose the KDE and GMM incorporated the KDE (GKDE) to segment spotted microarray image in this study.
2.3.1 Kernel Density Estimation (KDE)
The KDE with automatic bandwidth selection [41] is used to estimate the density function of pixel intensities for each spot. The Gaussian kernel function and 128 grid points are used for the KDE for each spot as Eq. (5) in section 2.2.1. The details are reported in the following algorithm.
Algorithm for Segmenting One Spot by KDE
Step 1: Input dataX ={x1,x1,L,xn}.
Step 2: Find 128 grid points that are equally spaced as Eq. (6).
Step 3: Calculate the data-driven bandwidth for KDE as Eq. (7).
Step 4: Calculate the KDE using Eq. (5).
Step 5: Search the cut-off point (CP) that is the first local minimum of the KDE at yj*
and let CP = yj*.
Step 6: Segment the pixelx into foreground ifi xi >CP, else into background.
2.3.2 Gaussian Mixture Model (GMM)
The GMM assumes that the distribution of foreground intensities is a Gaussian
distribution f1(μ1,σ12) with meanμ and variance1 ; while the distribution of background intensities is another Gaussian distribution with mean
2
and variance . Hence, the distribution of the intensity at every pixel in a spot can be modeled as a mixture of two Gaussian distributions as Eq. (18).
2
πm is the mixing (or prior) probability for the foreground and the background constrained by 0≤πm≤ and 1 π1+π2 =1. The foreground intensities typically include those of the signals and noise. Therefore, the mean foreground intensity usually exceeds the mean background intensity. Accordingly, the condition μ1 ≥μ2 is considered in this study, which is also commonly used in the literature of normal mixtures [28]. The log-likelihood of the observed data in the model of two mixtures is
Eq. (19). To estimate the above parameters, the EM algorithm can be applied [26]. The segmentation algorithm of one spot using the GMM is listed below.
Algorithm for Segmenting One Spot by GMM
Step 1: Input initial parameters: k = 0 and . In this study, the initial parameters are set as follows. Initial μ
( )k { m( )k , m( )k , m2( )k ,m 1, 2.
φ = π μ σ =
1 and μ2 are set to the first andthird quartiles of pixel intensities in one spot. Initial σ1 and σ2 are the standard deviations of the pixel intensities below thefirst quartile and above thanthe third quartile, respectively. Initial π1 and π2 values are set to 0.5.
Step 3: Calculate new estimates of
.
here, then the iteration is terminated. Otherwise, k Å k+1, , and the iteration goes to Step 2.
( )k (k 1) { m(k 1), m(k 1), m2(k 1),m 1, φ ←φ + = π + μ + σ + = 2}
Step 5: Segment the pixel x into foreground or background according to the i maximum of posterior probabilities with the final values of the parameters,
(k 1)=
2.3.3 GMM Incorporated with KDE (GKDE)
We can combine the methods of GMM and KDE, which will be abbreviated as GKDE. The GMM can provide the initial segmentation and the KDE can further improve the segmentation by relaxing the assumption of normality in the GMM. Once the foreground and background are found using GMM, the KDE can be applied to find their estimated densities. Then, a cut-off point for segmentation of pixels in a spot into two clusters is determinate by the cut-off point that has the near equality of two estimated densities. The details are reported below.
Algorithm for Segmenting One Spot by GKDE
Step 1: Segment a spot initially using the GMM algorithm in Section 2.3.2.
Step 2: Estimate the kernel densities for foreground ( ) and background ( ) following the steps in Eqs. (5) – (7).
fˆf fˆg
Step 3: Find a cut-off point (CP) that has the near equality of fˆf and fˆg. Step 4: Segment a spot as follows.
, ;
2.3.4 Empirical Evaluation of the Performance for Segmentation Methods
Spike genes (or spots) with known contents on microarrays are used in the empirical studies. The target ratios of spike genes thus represent the gold standard for evaluating the accuracy of segmentation methods investigated in this study. The sum of squared relative error (SSRE) and the sum of squared error (SSE) are used to evaluate accuracy according to Eqs. (20) and (21):
∑∑
= = ⎪⎭where is the feature estimated from the ratio of means between Cy3 and Cy5 arrays for the j
b
Tˆj,
th spike gene in the bth block, and Tj is target ratio of the jth spike gene.
The number of blocks is B = 32 and the number of spike genes is M = 8. The smallness of SSRE and SSE indicate the closeness to target ratios. Table 2.3.4-1 and 2.3.4-2 present those two types of spike genes with various target contents and ratios used in this study.
Table 2.3.4-1: One set of eight spike genes with target contents and ratios in real microarrays.
Spike Gene 1 2 3 4 5 6 7 8
Target Content
(Cy5: Cy3) 50:500 10:100 50:250 20:100 200:500 40:100 200:200 20:20 Target Ratio
1:10 1:10 1:5 1:5 1:2.5 1:2.5 1:1 1:1 (Cy5: Cy3)
Table 2.3.4-2: Another set of eight spike genes with target contents and ratios in real microarrays.
For those tw sets of microarrays are produced and each set of microarrays consists of one pair of two dye-swapped microarrays. Therefore, each type of spike gene is associated with eight microarrays, of which a total of 16 microarrays are tested herein.
The Pearson and concordance correlation coefficient [46] of two random variables Y1 and Y2 is shown as Eq. (22) and (23):
o types of spike genes, four
),
easure the reproducibility among the ing the log ratios of means in Cy5 to Cy3 dyes from one microarray image, which is expected to be close to 1.
These correlations are used in this study to m
expression of every gene and that of its duplicated spot us
These correlation coefficients of the swapped microarrays are also considered to evaluate the performance with reference to selected features with high log ratios of means in Cy5 to Cy3 dyes. As the dyes of Cy3 and Cy5 in the swapped arrays are exchanged, the negative correlation coefficient is obtained from the features of the swapped arrays and is expected to be close to -1.
3. Applications on MicroPET Images
This chapter explains the reconstruction and segmentation methods applied to microPET images. The FBP and OSEM are available methods of the built-in software of microPET manager V1.6.4 associated with the scanner to reconstruct microPET image. From our evaluation, the PDEM will demonstrate more accurate and less noisy reconstruction images than the FBP and OSEM. Moreover, the microPET manager V1.6.4 does not have any segmentation method for 3D images. Section 3.3 will present the GMM algorithm to segment 3D microPET images obtained from PDEM.