Anatomy assisted PET image reconstruction incorporating multi-resolution joint entropy

In previous work (Tang and Rahmim 2009a), we developed a one-step-late (OSL) MAP PET image reconstruction algorithm with the JE between anato-functional image intensities as the prior. The MAP estimate of the functional image X from the emission data g is given by

$$\begin{eqnarray}&&\mathbf{\hat{X}}=\underset{\mathbf{X}\geqslant 0}{{\text{arg max}}} \frac{P\left(\mathbf{g}\mid \mathbf{X}\right)P\left(\mathbf{X}\right)}{P\left(\mathbf{g}\right)},\end{eqnarray} \tag{ 1 }$$

where P () indicates the probability of the variable enclosed in parentheses. The prior on the image is assumed to follow a Gibbs distribution

$$\begin{eqnarray}&&P\left(\mathbf{X}\right)=\frac{1}{W} \text{exp} \left(-\beta V\left(\mathbf{X}\right)\right),\end{eqnarray} \tag{ 2 }$$

where W is a normalization factor and β V (X) is the weighted potential function, which has been conventionally used to penalize the inter-voxel intensity variations in the reconstructed image. In the JE-MAP algorithm, the JE H (X, Y) between the PET image X and MR image Y was applied to serve as the potential function. To introduce the WJE-MAP algorithm, let us denote the kth (k = 1, ..., K) wavelet analysis generated feature image pair of the PET image and the MR image by {f^k (X), f^k (Y)}. Instead of the intensity-only JE term, the potential function to apply is the weighted sum of the JE of feature image pairs represented by $\underset{k=1}{\overset{K}{\sum}} {{\beta}_{k}}H\left({{f}^{k}}\left(\mathbf{X}\right),{{f}^{k}}\left(\mathbf{Y}\right)\right)$, where β_k denotes the weighting hyperparameter for a given kth wavelet subband image pair.

For the WJE-MAP algorithm, the image reconstruction process is to seek for the image estimation $\mathbf{\hat{X}}$ that maximizes the log-posterior probability

$$\begin{eqnarray} &&\text{log} P\left(\mathbf{g}\mid \mathbf{X}\right)-\underset{k=1}{\overset{K}{\sum}} {{\beta}_{k}}H\left({{f}^{k}}\left(\mathbf{X}\right),{{f}^{k}}\left(\mathbf{Y}\right)\right)+s,\end{eqnarray} \tag{ 3 }$$

where the constant s represents the contribution of the W term in (2) and the denominator in (1). Modeling Poisson statistics for the first term in (3) and utilizing the OSL EM approach (Green 1990), we arrive at the iterative procedure

$$\begin{eqnarray}&&x_i^{{\rm{new}}} = \frac{{x_i^{{\rm{old}}}}}{{{{\left. {\sum\nolimits_j {{c_{ij}}} \, + \sum\nolimits_{k = 1}^K {{\beta _k}\frac{{\partial H({f^k}({\bf{X}}),{f^k}({\bf{Y}}))}}{{\partial {x_i}}}} } \right|}_{{x_i} = x_i^{{\rm{old}}}}}}}\sum\limits_j {\frac{{{c_{ij}}{g_j}}}{{\sum\nolimits_i {{c_{ij}}x_i^{{\rm{old}}}} }}} ,\end{eqnarray} \tag{ 4 }$$

where the new estimate of voxel i in image X is updated from the old estimate. The bin j of emission data g is represented by g_j and c_ij denotes an element of the system matrix modeling the contribution of voxel i to projection bin j.

We propose to use the subbands generated by wavelet analysis of PET and MR image pairs as feature images from which JEs are calculated. A three-dimensional (3D) image can be decomposed into low frequency (L) and high frequency (H) components, via the approximation and detail wavelet filters in a particular direction. This approach can then be repeated along a second dimension, generating LL, LH, HL, HH and finally performed along the third dimension to generate 8 image subbands (each 1/8 of total size, i.e. 1/2 size in each direction), as depicted in figure 1. In this work, we use the nearly symmetric orthogonal wavelet bases in Abdelnour and Selesnick (2001) to decompose the PET and MR image pairs.

Zoom In Zoom Out Reset image size

The totality of the 8 subbands for each 3D image will have the same total size as the original image. It is then possible to utilize the entire information provided by some or all of the wavelet subbands for the JE computations. Similar to what was used in Xu and Chen (2007), we first considered the JE of approximation (LLL) and first order detail (HLL, LHL, LLH) wavelet subbands, the former integrating the similarity of image intensity features and the latter group having features comparable to gradients (depicted in figure 2). After analyzing the results from the 4-subband WJE algorithm, we moved onto using the JE from all the 8 subbands as the potential function.

Zoom In Zoom Out Reset image size

Here we note that the proposed MAP EM algorithm (4) involves calculation of the derivatives of the subband JE terms with respect to intensity of individual voxels in the original full PET image ∂H(f^k (X), f^k (Y))/∂x_i . We are able to show that this calculation can be performed from the derivatives of the subband JE terms with respect to the corresponding subband voxels similar to how wavelet synthesis is achieved. For simplicity, let us consider wavelet analysis along a particular dimension and denote the wavelet subband image k by ${{{\mathbf{\tilde{X}}}}^{k}}\triangleq {{f}^{k}}\left(\mathbf{X}\right)$, which is given by the convolution-downsampling operation

$$\begin{eqnarray}&&{{\mathbf{\tilde{X}}}^{k}}={{\Downarrow}_{2}}\left\{\mathbf{X} \otimes \text{}\text{}{{\mathbf{W}}^{k}}\right\}\text{or}\tilde{x}{{_{j}^{k}}_{{}}}=\underset{u}{\sum} {{x}_{u}}w_{2j-u}^{k},\end{eqnarray} \tag{ 5 }$$

where W^k is the wavelet analysis filter and the downsampling ⇓₂ is reflected in the 2j index in the wavelet filter. Then we invoke the chain rule

$$\begin{eqnarray}&&\frac{\partial H\left({{{\mathbf{\tilde{X}}}}^{k}},{{{\mathbf{\tilde{Y}}}}^{k}}\right)}{\partial {{x}_{i}}}=\underset{l}{\sum} \frac{\partial H\left({{{\mathbf{\tilde{X}}}}^{k}},{{{\mathbf{\tilde{Y}}}}^{k}}\right)}{\partial \tilde{x}_{l}^{k}}\frac{\partial \tilde{x}_{l}^{k}}{\partial {{x}_{i}}}.\end{eqnarray} \tag{ 6 }$$

The first term on the right hand side, which is the derivate of subband image JE with respect to voxels in the PET subband image, may now be calculated exactly the same way as proposed in Tang and Rahmim (2009a). In brief, the joint probability density for JE calculation is approximated by Parzen density estimate through superposition of Gaussian densities (Duda et al 2001).

Let's refer the derivate to the (subband) PET image voxel from which the (subband) JE is generated as the SDJE (standard derivative of JE). From the definition of ${{{\tilde{x}}}_{j}}$ in (5), it follows that

$$\begin{eqnarray}&&\frac{\partial H\left({{{\mathbf{\tilde{X}}}}^{k}},{{{\mathbf{\tilde{Y}}}}^{k}}\right)}{\partial {{x}_{i}}}=\underset{l}{\sum} \text{SDJ}{{\text{E}}_{l}}\left({{\mathbf{\tilde{X}}}^{k}},{{\mathbf{\tilde{Y}}}^{k}}\right)w_{2l-i}^{k}.\end{eqnarray} \tag{ 7 }$$

Thus once the SDJE images are calculated, the overall derivative can be derived from a 'convolution-like' calculation. Here one may notice an interesting parallel with the wavelet synthesis process: for a synthesis filter ${{\mathbf{\bar{W}}}^{k}}$ defined such that ${{\bar{w}}_{m}}\triangleq {{w}_{-m}}$ and using the upsampling operation ⇑₂ (inserting zeros between values at odd positions), it can be shown that

$$\begin{eqnarray}&&\frac{\partial H\left({{{\mathbf{\tilde{X}}}}^{k}},{{{\mathbf{\tilde{Y}}}}^{k}}\right)}{\partial {{x}_{i}}}=\underset{l}{\sum} {{\left[{{\Uparrow}_{2}}\left\{\text{SDJ}{{\text{E}}_{l}}\left({{\mathbf{\tilde{X}}}^{k}},{{\mathbf{\tilde{Y}}}^{k}}\right)\right\}\right]}_{l}}\bar{w}_{i-l}^{k}.\end{eqnarray} \tag{ 8 }$$

This now looks similar to how wavelet synthesis process is performed (i.e. upsampling each subband, followed by convolution with the synthesis filter). The only difference is that we first perform SDJE calculations on the wavelet-generated PET-MR subband pairs. Figure 3 shows the procedure explained above for calculating the derivatives of the subband JEs with respect to the original PET image voxel intensities.

Zoom In Zoom Out Reset image size

We performed simulation study using the BrainWeb phantom (Collins et al 1998). An anatomical model from the BrainWeb database was used to create PET activity and attenuation images. The model consists of a set of 'fuzzy' tissue membership classes including gray matter (GM), white matter (WM), cerebrospinal fluid (CSF) and other 8 membership volumes. A voxel value in a volume reflects the proportion of that tissue presented in that voxel. For the anatomical model, the BrainWeb provides simulated T1-weighted MR image for specified SPLASH parameters with (1 mm)³ voxels. The simulator uses first-principles modeling based on the Bloch equations to implement a discrete-event simulation of nuclear MR signal production and realistically models noise and partial volume effects of the image production process.

We created the PET distribution by assigning a clinically realistic activity of 12 500 Bq ml⁻¹ in GM, 3250 Bq ml⁻¹ in WM, 0 Bq ml⁻¹ in air, CSF and bone and 1000 Bq ml⁻¹ in all other tissues including skin, muscle, connective tissue and fat (Vunckx et al 2012). Attenuation coefficients were assigned as 0 cm⁻¹ for air, 0.146 cm⁻¹ for bone and 0.096 cm⁻¹ for other tissues to build the attenuation map. The PET (and MR) image was interpolated to create voxel size of (1.22 mm)³ and image dimension of 256 × 256 × 207 for simulation in the geometry of the high resolution research tomograph (HRRT) (Sossi et al 2005). The simulated PET activity, attenuation and MR images (central transaxial slice) are shown in figure 4. Emission data corresponding to three noise levels were simulated, corresponding to a low-noise (10 min), medium-noise (5 min) and a high-noise (3 min) FDG PET measurement, with decay, normalization and attenuation effects taken into account. The three noise levels are ∼1.8 times apart in terms of total counts. Image reconstruction was performed using the maximum likelihood (ML) algorithm (implemented using ordered-subset (OS) EM with 16 subsets in each iteration), the intensity-only JE-MAP algorithm and the WJE-MAP algorithm.

Zoom In Zoom Out Reset image size

To quantitatively evaluate the images reconstructed from different algorithms, we used the tradeoffs between noise and bias on the GM and WM regions as the metrics. The GM (or WM) region was defined to include the voxels that more than 95% belong to the GM (or WM) based on the fuzzy tissue class model. For each region of interest (ROI), the normalized mean squared error (NMSE) was calculated to measure bias

$$\begin{eqnarray}&&\text{NMSE}={{\left(\frac{\bar{X}-{{{\bar{X}}}_{\text{true}}}}{{{{\bar{X}}}_{\text{true}}}}\right)}^{2}},\end{eqnarray} \tag{ 9 }$$

where $\bar{X}=\frac{1}{N}\underset{i=1}{\overset{N}{\sum}} {{X}_{i}}$ ; X_i denotes the ith voxel reconstructed activity value and N is the number of voxels in the ROI. Similarly, ${{\bar{X}}_{\text{true}}}$ is the regional mean value of the phantom activity image. For each ROI, the noise was measured using the normalized standard deviation (NSD) as calculated by

$$\begin{eqnarray}&&\text{NSD}=\frac{\sqrt{\frac{1}{N-1}{\sum}_{i=1}^{N}{{\left({{X}_{i}}-\bar{X}\right)}^{2}}}}{{\bar{X}}},\end{eqnarray} \tag{ 10 }$$

where X_i and $\bar{X}$ are defined as those in (9).

We applied the WJE and intensity-only JE MAP techniques to two PET imaging datasets (and corresponding MR images) with tracers for different diagnostic purposes. One set of PET data was acquired using the tracer [11C] DPA-173, which is a radioligand that binds to the translocator protein (TSPO) and serves as a marker of the neuroinflammatory burden (Endres et al 2009). The other tracer was [18F] Florbetapir, which is an FDA approved amyloid-beta imaging agent for clinical evaluation of late-life cognitive impairment in Alzheimer's disease (Wong et al 2010). PET scans were performed on the second-generation Siemens HRRT (Sossi et al 2005), an LSO-based, 2.5 mm-resolution, dedicated brain PET scanner. The data was first reconstructed using the OS EM algorithm (10 iterations and 16 subsets), with correction for attenuation, normalization, deadtime, scatter and randoms (Rahmim et al 2005). The PET image space consisted of cubic voxels with sizes (1.22 mm)³, spanning dimensions of (31 cm)² transaxially and 25 cm axially. We studied the DPA-173 data with accumulation from 0 to 30 min post-injection while for the Florbetapir data we worked on data with two noise levels. The low noise level has data from 0 to 70 min and the higher noise level has duration time spanning 50–70 min. The structural MR images were created for the corresponding subjects using the T1-weighted MPRAGE sequence. Before MR assisted PET reconstruction was performed, the images were aligned to initially reconstructed PET images using an in-house software implementing rigid mutual-information registration.

To evaluate the performance of the proposed technique and compare it with other methods, we calculated the noise over the ROIs using (10). The ROIs were drawn using an in-house software on the MR image and applied to the registered reconstructed PET images. As the true intensity values of patient experiments are not available, we plotted the regional mean value versus the noise along with iterations on the specified ROIs of the reconstructed PET images.

As stated in section 2.2, we first applied the JEs of the approximation (LLL) and the first order detail (HLL, LHL, LLH) wavelet subbands in the reconstruction (figure 2). We were able to see improvement in the reconstructed image components made of the subbands involved in the prior formation (Tang and Rahmim 2009b). Therefore we moved on to use all the 8 subbands to form the prior. We show the results from using the all-subband prior here and leave the analysis of results from different combinations of subbands in the Discussion section.

Besides the formation of the prior, the performance of a MAP reconstruction algorithm heavily depends on the weighting factor. We chose the weighting factors that result in similar end-iteration bias before comparing the NSD versus NMSE tradeoff from the MAP algorithms with the ML algorithm. Equal weights were used for the WJE subbands in this study. For the three cases at low-, medium- and high-noise levels, the NSD versus the NMSE along with iteration number from the ML, JE-MAP and WJE-MAP algorithms are plotted for GM and WM in figure 5, respectively. For the low-noise level simulation, both the JE-MAP and the WJE-MAP algorithms result in better noise versus bias performance than the ML algorithm, in both the GM (figure 5(a1)) and the WM (figure 5(a2)) regions. The results from the two MAP algorithms are close to comparable. In the medium-noise level case, the WJE-MAP algorithm and JE-MAP algorithm perform similarly in the WM region (figure 5(b2)). The WJE-MAP algorithm demonstrates its advantage over the JE-MAP algorithm in the GM region (figure 5(b1)) reaching lower noise level with comparable bias. In the high-noise level simulation, the performance of WJE-MAP algorithm clearly surpasses the JE-MAP algorithm in both the GM (figure 5(c1)) and WM (figure 5(c2)) regions in terms of noise versus bias tradeoff. The reconstructed images (center transaxial slice, 5th iteration) of the three noise-level simulations from different reconstruction algorithms are plotted in figure 6. The maximum gray level of the reconstructed images is set as twice of the maximum value in the phantom image, so that it is easier to compare the images from different reconstruction algorithms and from different noise levels.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.1. DPA-173 results.

In figure 7, we plot the regional noise versus mean value along with iteration number for different ROIs in the DPA-173 patient images reconstructed from using the ML, JE-MAP and WJE-MAP algorithms. In most of the regions, the WJE-MAP algorithm suppresses noise while arriving at mean regional values close to what ML algorithm achieves. Its improvement on noise versus bias tradeoff over the ML and also the JE-MAP algorithms is well demonstrated. The relative performance of the reconstruction algorithms in this patient study is similar to that in the high-noise simulation study. For visual comparison, figure 8 shows the transaxial center slice of the reconstructed images (10th iteration) using the different algorithms (together with the MR image used for the MAP reconstructions and the ROIs). The reconstructed images are shown with the same gray level limit for easy comparison.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2.2. Florbetapir results.

We plot the regional noise versus mean value along with iteration number for different ROIs reconstructed from the 0–70 min Florbetapir data using the ML, JE-MAP and WJE-MAP algorithms in figure 9. This result is comparable to the result from the low-noise simulation case, that is, the performance of the WJE-MAP and the JE-MAP algorithms are comparable (note the scale of the mean value is very small and the difference between the mean values from the two MAP algorithms is about ∼2% .) The regional noise versus mean value along with iteration number plots for images reconstructed from the 50–70 min data are shown in figure 10. Similar to the results from high-noise simulation case, the WJE-algorithm demonstrates improvement over the ML and the JE-MAP algorithms on the noise versus bias tradeoff in this higher noise Florbetapir patient study. Figures 11 and 12 show the transaxial center slice of the reconstructed images (10th iteration) using the different algorithms (together with the MR image used for the MAP reconstructions and the ROIs). In each figure, the reconstructed images are shown with the same gray level limit for straightforward comparison.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size