# Blind Deblurring Reconstruction Technique with Applications in PET Imaging

Authors: Heng Li; Osama R. Mawlawi; Ronald X. Zhu; Yibin Zheng
Journal: International Journal of Biomedical Imaging
 (2009)

---

## Abstract

We developed an empirical PET model taking into account system blurring and a blind iterative reconstruction scheme that estimates both the actual image and the point spread function of the system. Reconstruction images of high quality can be acquired by using the proposed reconstruction technique for both synthetic and experimental data. In the synthetic data study, the algorithm reduces image blurring and preserves the edges without introducing extra artifacts. The localized measurement shows that the performance of the reconstruction image improved by up to 100%. In experimental data studies, the contrast and quality of reconstruction is substantially improved. The proposed method shows promise in tumor localization and quantification.


---

## Body

## 1. Introduction

In the past decade, positron emission tomography (PET) has rapidly become a popular diagnostic imaging modality for tumor detection and cancer staging. PET generates three-dimensional (3D) tomographic images of the distribution of positron-emitting radiotracers within an object, which provides quantitative functional information in vivo. However, the physics of the photon emission and detection process limits the resolution of PET images. The sensitivity and resolution of the PET system is a complex function of many factors, such as positron range, scattering, medium attenuation, photon noncollinearity, detector size, and intercrystal crosstalk. Iterative methods are often used in PET reconstruction because of the ability to handle the imaging system more precisely and to incorporate a prior information in the process.Figure1 illustrates a PET imaging system [1]. A positron is emitted following a nuclear transmutation. After traveling a short distance (positron range), the positron annihilates with an electron and produces two photons of 0.511 MeV almost 180° to each other. The photons are then detected by the detector ring. Finally, by recognizing the coincidence photons, a line of response (LOR) can be established to determine the origin of the positron emission. It should be noted that LOR is not a real line due to the finite size of the detector elements. As shown in Figure 1, photons may be scattered to travel along a totally different direction upon reaching the detector ring, and the origin of the photon pair cannot easily be described by the LOR. Therefore, reconstructions obtained from iterative algorithms that determine the system matrix solely based on LOR or extended LOR are often blurry and contaminated by artifacts.Figure 1
Diagram of PET data collection.We recently developed a blind deblurring reconstruction technique to reduce blurring in pinhole SPECT imaging [2]. It is natural to apply a similar technique to PET imaging, because SPECT and PET share a similar underlying physic. Our approach is to model the discrepancy of LOR, or the measurement, and the true emission by an independent random variable E, as demonstrated in Figure 1, and to derive the maximum likelihood (ML) solution of the new system incorporating measurement error. In this report, we describe the development of the reconstruction algorithm.


## 2. Blind Deblurring Reconstruction

Equation (1) illustrates the widely used Expectation-Maximization (EM) algorithm for PET reconstruction [3]:(1)λ̂in+1=λ̂in∑jpij∑jnjpij∑i'λ̂i'npij,
where pij, element (i,j) of system matrix P, denotes the probability of detecting an emission from voxel i,i=1,…,S, at detector pair j,j=1,…,T. The emitted photon pairs Y, radiotracer concentration  Λ, and detected photon pairs N are as follows:(2)N=PY,nj=∑ipijyi,yi~Poisson(λi).
The system matrix P can be factorized and correction factors can be applied as follows [4]:(3)P=Pdet.sensPdet.blurPattnPgeomPpositron,
where Pgeom is the geometric projection matrix with each element (i,j) representing the probability of a photon pair produced in voxel i reaching detector pair j, ignoring  attenuation, and assuming perfect photon-pair collinearity. The values in the matrix are defined by the geometries of each LOR. Pattn and Pdet.sens are the attenuation correction and detector sensitivity normalization factors, respectively.Ppositron, the positron range factor, is usually omitted in 18F studies in which the positron range is submillimeter, and Pdet.blur is the detector blurring factor used to model photon-pair noncollinearity, intercrystal scatter, and penetration, which cause the same photon detected by several adjusting detector elements. In general, the correction can be viewed as weighting and broadening the LOR on top of Pgeom.In this work, instead of adding correction factors to the system matrixP, we model the uncertainty caused by positron range, photon emission angle, scatter, and detector response as a blurring factor E that introduces measurement error, which is a random variable independent from emission Y, and modify the system model accordingly. The measurement error introduces blur in the reconstruction image, and a blurred PET reconstruction can be viewed as the convolution of a low-pass point spread function (PSF) with the actual image, where both the PSF and the actual image are unknown in practice. Inspired by the blind deconvolution algorithm introduced by Holmes [5] and along the same line of research we have applied to SPECT imaging [2], we formulated a blind deblurring reconstruction algorithm. The algorithm is a modified EM algorithm that includes a convolution kernel to model the blurring factor E. The algorithm consists of two iterative updates, instead of one, to reconstruct both the object and the PSF.In a PET imaging system, supposesk,k=1⋯Λ,sk∈{1⋯S}, where Λ is the total number of emitted photons, denotes the index of location from which the kth pair of photons are emitted, and tk,tk∈{1⋯T} denotes the location where the kth pair of photons are detected, as shown in Figure 1. We call these emission locations “true emission points.” A finite number of these points form an inhomogeneous Poisson random-point process having the intensity function λi. With the presence of measurement error, the positional measurement of each emission point is corrupted by a random translation. Let ek,ek∈{1⋯S} denote this error vector, then the measured data for detector pixel j is related to sk and ek by(4)P(tk=j∣sk+ek=i)=pij,nj=∑k=1Λp(sk+ek)j.
Here, ek is statistically independent of all sk’s, and they are assumed to be all statistically independent of each other for all photons emitted and identically distributed with a probability density gi, which is also the PSF of the PET system. In addition, the set of error vectors E={e1⋯eΛ} also constitutes an inhomogeneous Poisson random point process with intensity function γi=Λgi [6].Now letyi be the actual number of photon pairs emitted from voxel i, and let bi be the corresponding error vectors within voxel i; they then follow the Poisson distribution with mean λi and γi, respectively, from the results above. We then use the EM algorithm to find the maximum likelihood solution of the system. In our application, nj and  Λ are known measured data, whereas yi and bi are unknown data. Here we note the set of true emission vectors Y and the set of error vectors B as(5)Y={y1,y2,y3,…},B={b1,b2,b3,…}.
Then the log likelihood of I can be expressed as(6)L(I∣λ)=-∑iλi+∑iln(λi)yi,
where λ is the vector notation for all λi. Similarly, the log likelihood of B can be written as(7)L(Y∣g,Λ)=-∑iγi+∑iln(γi)bi,
where g is the vector notation for all gi. The log likelihood of the complete data then can be equivalently expressed in two ways, assuming Y or B being known, that is,(8)L(I,Y∣λ,g,Λ)=L(I∣λ,g,Λ)+l(Y∣λ,g,Λ)=-∑iλi+∑iln(λi)yi+∑k=1Λln(g(sk))=-∑iγi+∑iln(γi)bi+∑k=1Λln(λ(ek)Λ).
By maximizing (8) using a derivation similar to those of Holmes [5] and Li [2], the following iteration can be shown to converge to the maximum likelihood estimate of λi and gi:(9)λ̂in+1=λ̂in∑jnj∑kgi-i''pi''j∑i'(λ̂i'n*gi')pi'j,(10)ĝin+1=ĝinΛ∑jnj∑i''λi-i''pi''j∑i'(λi'*ĝi'n)pi''j,
where * denotes convolution. The initial λi0 is an image of all 1’s, and gi0 is the same image normalized to 1. Equations (9) and (10) are then evaluated to acquire a new pair of estimates of λ and  g. The PSF of the system is assumed to be real, nonnegative, band limited, and limited in extend. Letting Fz be the frequency components of the PSF that are known to be zero and Fr be the space components of the PSF that are known to be zero, the band-limited and limited-extend constraints are incorporated by executing the following steps in each iteration.(1)
The Fourier transform ofĝn+1 is taken, and any frequency components that lie within Fz are set to zero.(2)
The inverse Fourier transform in step 1 above is taken, and any negative or complex values or values withinFr in the spatial domain are set to zero.The first step of the process ensures the band-limited constraint, and the second step ensures the reality, nonnegativity, and limited-extend of the PSF. Realness and nonnegativity are implicitly applied toλ. Equations (9) and (10) and steps (1) and (2) are then iterated until convergence occurs.The blind deblurring reconstruction algorithm estimates both the spatial radioactivity distribution and the system PSF from the set of blurred projection images. The iteration for reconstruction can be understood as replacing the forward projector in the original EM (denominator of  (9) with the new projector using the convolved radioactivity map, and the iteration for solving the PSF can be understood as blind deblurring. This iteration differs from the general image-blind deconvolution in the sense that the kernel is partly known: pij, the system matrix, is in fact part of the blurring kernel. The more precise the model pij is, the closer the remainder of blurring kernel, or g, is to a true delta function. In addition, instead of deconvolving an image where both the input and output are 2D images, the input of blind deblurring reconstruction is a series of projection images, and the output is a 3D-image array. This property gives us much more knowledge of the noise distribution within the object, because instead of a single-shot image, we now essentially have multiple samples for each point in the 3D array (although mixed with other points). Both simulation and experimental data were used to validate and evaluate the performance of the blind deblurring EM (BDEM) technique.


## 3. Methods

### 3.1. Monte Carlo Simulation

PET simulations were performed utilizing the NCAT phantom covering the chest region [7, 8]. The PET emission data were generated using the Monte Carlo method, with 1 million counts per slice using a geometry corresponding to the design of the GE discovery DSTE PET/CT scanner (GE Healthcare, Milwaukee, WI, USA). The images were reconstructed using the standard EM and the BDEM.   The reconstruction image sets were then evaluated using visual inspections and line profiles. The EM reconstruction image was postsmoothed using a Gaussian kernel of 1 pixel full width at half medium, and no post processing was done for BDEM.


### 3.2. Convergence Study

One major concern with regard to the type of blind deblurring technique used is the effectiveness and convergence of the algorithm. Two hundred iterations of EM and BDEM reconstructions of the NCAT phantom were performed to evaluate the convergence property of the BDEM algorithm. In each iteration, the measurement log-likelihoodL of the reconstructed image λ for EM was calculated as follows [6]:(11)L(λ∣N)=∑j(njlog∑ipijλi-∑ipijλi-lognj!)
whereas for BDEM, the log-likelihood L of the reconstruction image was calculated as(12)L(λ∣N)=∑j(njlog∑ipij(λi*gi)-∑ipij(λi*gi)-lognj!).
Two different settings of limited-extend constraint were used for BDEM reconstruction, and the difference in log-likelihood was compared.


### 3.3. Comparison with Image Deconvolution

The blurred reconstruction image could always be modeled as a convolution of the true radioactivity with a blurring kernelhi:(13)fi=λi*hi,
where fi is the standard EM reconstruction image and * denotes 2D linear convolution. Assuming hi can be estimated, one could first compute fi using standard EM iterations  (with λ̂i being replaced by f̂i) and then deconvolve f̂i with hi. However, the λ̂i so obtained is not the maximum likelihood estimate of λi given the measurement error; also, the kernel hi is generally a complex unknown function and is hard to measure. We used the Wiener filter and image blind deconvolution algorithm to denoise and deconvolve the EM reconstruction image, and compared the results with the BDEM reconstruction.


### 3.4. Physical Phantom Study

A physical phantom study with a water-tank phantom containing 4 spheres of diameters 1.5 cm, 2.0 cm, 3.0 cm, and 3.0 cm was conducted to assess the performance of a BDEM algorithm for a real scanner. Both the water tank and the spheres were filled withFDG18, with activity concentrations of 0.20μCi/cc and 1.01 0.20μCi/cc, respectively. The tank was imaged with a GE Discovery DSTE PET/CT scanner (GE Healthcare) in 2D list-mode; the sinogram was extracted and reconstructed using the EM and BDEM techniques.


## 3.1. Monte Carlo Simulation

PET simulations were performed utilizing the NCAT phantom covering the chest region [7, 8]. The PET emission data were generated using the Monte Carlo method, with 1 million counts per slice using a geometry corresponding to the design of the GE discovery DSTE PET/CT scanner (GE Healthcare, Milwaukee, WI, USA). The images were reconstructed using the standard EM and the BDEM.   The reconstruction image sets were then evaluated using visual inspections and line profiles. The EM reconstruction image was postsmoothed using a Gaussian kernel of 1 pixel full width at half medium, and no post processing was done for BDEM.


## 3.2. Convergence Study

One major concern with regard to the type of blind deblurring technique used is the effectiveness and convergence of the algorithm. Two hundred iterations of EM and BDEM reconstructions of the NCAT phantom were performed to evaluate the convergence property of the BDEM algorithm. In each iteration, the measurement log-likelihoodL of the reconstructed image λ for EM was calculated as follows [6]:(11)L(λ∣N)=∑j(njlog∑ipijλi-∑ipijλi-lognj!)
whereas for BDEM, the log-likelihood L of the reconstruction image was calculated as(12)L(λ∣N)=∑j(njlog∑ipij(λi*gi)-∑ipij(λi*gi)-lognj!).
Two different settings of limited-extend constraint were used for BDEM reconstruction, and the difference in log-likelihood was compared.


## 3.3. Comparison with Image Deconvolution

The blurred reconstruction image could always be modeled as a convolution of the true radioactivity with a blurring kernelhi:(13)fi=λi*hi,
where fi is the standard EM reconstruction image and * denotes 2D linear convolution. Assuming hi can be estimated, one could first compute fi using standard EM iterations  (with λ̂i being replaced by f̂i) and then deconvolve f̂i with hi. However, the λ̂i so obtained is not the maximum likelihood estimate of λi given the measurement error; also, the kernel hi is generally a complex unknown function and is hard to measure. We used the Wiener filter and image blind deconvolution algorithm to denoise and deconvolve the EM reconstruction image, and compared the results with the BDEM reconstruction.


## 3.4. Physical Phantom Study

A physical phantom study with a water-tank phantom containing 4 spheres of diameters 1.5 cm, 2.0 cm, 3.0 cm, and 3.0 cm was conducted to assess the performance of a BDEM algorithm for a real scanner. Both the water tank and the spheres were filled withFDG18, with activity concentrations of 0.20μCi/cc and 1.01 0.20μCi/cc, respectively. The tank was imaged with a GE Discovery DSTE PET/CT scanner (GE Healthcare) in 2D list-mode; the sinogram was extracted and reconstructed using the EM and BDEM techniques.


## 4. Results and Discussion

### 4.1. Monte Carlo Simulation

Figure2 compares the image slices through a 2-cm-diameter tumor in the right lung produced by the EM and BDEM reconstruction techniques. It is clear that the image quality of the BDEM reconstruction is superior to that of the EM reconstruction. The contrast of the lesion of interest is greatly improved, edges are preserved, and artifacts are suppressed. The transverse-view line profile across the tumor also confirms the improvement. The peak contrast of the tumor in the BDEM reconstruction is double that of the EM reconstruction. Although no postsmoothing was performed for the BDEM reconstruction, the noise level in the background of the BDEM reconstruction is lower than that of the EM reconstruction, with postsmoothing by a 1-pixel full width at half medium Gaussian kernel. The mean and standard deviation of a small region of interest in the background region are 36.20 and 25.12 (resp.) for the Gaussian-smoothed EM and 34.83 and 21.28 (resp.) for BDEM, with no postprocessing.Reconstruction of NCAT utilizing EM and BDEM algorithms. From top to bottom are the coronal, sagittal, and transverse views.
(a)
BDEM(b)
EM(c)
Profile through a 2-cm-diameter tumor as displayed on the transverse viewFigure 3
Log-likelihood of EM and BDEM algorithms as a function of iteration number. The two BDEM plots have different PSF extend constraints; small PSF extend limits the PSF to10×10, whereas large extend limits the PSF to 20×20. The image size is 128×128 pixels.


### 4.2. Convergence Study

Figure2 shows the log-likelihood as a function of iteration numbers as calculated in  (11) and (12). The measurement log-likelihood of BDEM is a monotonic increasing function of iteration number, which confirms the convergence of the algorithm. It also varies with different PSF constraints and lies within the envelope defined by the EM algorithm, as demonstrated in the figure. Therefore, setting a proper constraint on PSF is important. If the initial PSF is a delta function and no other constraint is set, the BDEM would regress to the regular EM solution. Figure 2 also shows the log-likelihood of BDEM with a 10×10- and a 20×20-PSF spatial constraint being applied for a 128×128 reconstruction. The small discrepancy in the convergence of two different constraints indicates that BDEM is not a strong function of PSF constraint: as long as the constraint is reasonable, BDEM would converge to very similar solutions. The figure also indicates that BDEM converges a little faster than conventional EM in terms of number of iterations. It should be noted, however, that the amount of computation for each BDEM iteration is about twice that of a computation of an EM iteration. We have not tested this premise, because the algorithm converges reasonable fast even without order-subset, but the principle of order-subset EM (OSEM) should be applied if necessary.


### 4.3. Comparison with Image Deconvolution

The results obtained by deblurring the EM reconstruction image with Wiener filter and the image-blind deblurring technique were displayed and compared with the results of the EM and BDEM reconstructions (Figure4). These image deconvolution results are clearly inferior compared to BDEM reconstruction. The images are still noisy, and the visual improvement of deconvolution is limited. Our explanation is that the BDEM algorithm utilizes the statistical information of both emission and measurement noise in the projection data, which is loss in the reconstruction image. The noise in the projection data commonly appears as streaks or other local or global artifacts in the reconstruction image, which makes it much harder to either identify or clear just from the reconstruction image.EM and BDEM compared to image deblurring techniques.
(a)
BDEM(b)
EM(c)
Wiener filter(d)
Image blind deconvolutionWater phantom reconstruction. From top to bottom are the coronal, sagittal, and transverse views.
(a)(b)(c)


### 4.4. Physical Phantom Study

Figure4 shows one slice of water phantom with spheres. As with the Monte Carlo study, the image contrast of the BDEM has been improved over the conventional EM, and the mass and edges are well preserved in the reconstruction. The image blurring is reduced, and the contrast is greatly improved (up to 50%), as observed from the line profile from the transverse slice.  The shape of the lesions is not distorted, indicating that the algorithm preserves the general shape of objects.


## 4.1. Monte Carlo Simulation

Figure2 compares the image slices through a 2-cm-diameter tumor in the right lung produced by the EM and BDEM reconstruction techniques. It is clear that the image quality of the BDEM reconstruction is superior to that of the EM reconstruction. The contrast of the lesion of interest is greatly improved, edges are preserved, and artifacts are suppressed. The transverse-view line profile across the tumor also confirms the improvement. The peak contrast of the tumor in the BDEM reconstruction is double that of the EM reconstruction. Although no postsmoothing was performed for the BDEM reconstruction, the noise level in the background of the BDEM reconstruction is lower than that of the EM reconstruction, with postsmoothing by a 1-pixel full width at half medium Gaussian kernel. The mean and standard deviation of a small region of interest in the background region are 36.20 and 25.12 (resp.) for the Gaussian-smoothed EM and 34.83 and 21.28 (resp.) for BDEM, with no postprocessing.Reconstruction of NCAT utilizing EM and BDEM algorithms. From top to bottom are the coronal, sagittal, and transverse views.
(a)
BDEM(b)
EM(c)
Profile through a 2-cm-diameter tumor as displayed on the transverse viewFigure 3
Log-likelihood of EM and BDEM algorithms as a function of iteration number. The two BDEM plots have different PSF extend constraints; small PSF extend limits the PSF to10×10, whereas large extend limits the PSF to 20×20. The image size is 128×128 pixels.


## 4.2. Convergence Study

Figure2 shows the log-likelihood as a function of iteration numbers as calculated in  (11) and (12). The measurement log-likelihood of BDEM is a monotonic increasing function of iteration number, which confirms the convergence of the algorithm. It also varies with different PSF constraints and lies within the envelope defined by the EM algorithm, as demonstrated in the figure. Therefore, setting a proper constraint on PSF is important. If the initial PSF is a delta function and no other constraint is set, the BDEM would regress to the regular EM solution. Figure 2 also shows the log-likelihood of BDEM with a 10×10- and a 20×20-PSF spatial constraint being applied for a 128×128 reconstruction. The small discrepancy in the convergence of two different constraints indicates that BDEM is not a strong function of PSF constraint: as long as the constraint is reasonable, BDEM would converge to very similar solutions. The figure also indicates that BDEM converges a little faster than conventional EM in terms of number of iterations. It should be noted, however, that the amount of computation for each BDEM iteration is about twice that of a computation of an EM iteration. We have not tested this premise, because the algorithm converges reasonable fast even without order-subset, but the principle of order-subset EM (OSEM) should be applied if necessary.


## 4.3. Comparison with Image Deconvolution

The results obtained by deblurring the EM reconstruction image with Wiener filter and the image-blind deblurring technique were displayed and compared with the results of the EM and BDEM reconstructions (Figure4). These image deconvolution results are clearly inferior compared to BDEM reconstruction. The images are still noisy, and the visual improvement of deconvolution is limited. Our explanation is that the BDEM algorithm utilizes the statistical information of both emission and measurement noise in the projection data, which is loss in the reconstruction image. The noise in the projection data commonly appears as streaks or other local or global artifacts in the reconstruction image, which makes it much harder to either identify or clear just from the reconstruction image.EM and BDEM compared to image deblurring techniques.
(a)
BDEM(b)
EM(c)
Wiener filter(d)
Image blind deconvolutionWater phantom reconstruction. From top to bottom are the coronal, sagittal, and transverse views.
(a)(b)(c)


## 4.4. Physical Phantom Study

Figure4 shows one slice of water phantom with spheres. As with the Monte Carlo study, the image contrast of the BDEM has been improved over the conventional EM, and the mass and edges are well preserved in the reconstruction. The image blurring is reduced, and the contrast is greatly improved (up to 50%), as observed from the line profile from the transverse slice.  The shape of the lesions is not distorted, indicating that the algorithm preserves the general shape of objects.


## 5. Conclusion

In this work, we demonstrated highly desired reconstruction results with no complex assumption about the imaging system or the object. The blind deblurring reconstruction technique can significantly improve the quality and contrast of the reconstruction as demonstrated in both simulation and experimental scans. This algorithm does not only reconstruct the radiotracer map, but also determines the PSF of the system. The masses and edges are well preserved in the reconstruction image, which can be extremely useful when doctors need to localize, segment, or tally the activities in the possible tumor. Future studies would involve the application of this system in patient imaging and quantitative studies.



---

*Source: 718157-2009-06-11.xml*