|Home | About | Journals | Submit | Contact Us | Français|
The problem we address is the optimization and comparison of window-based scatter correction (SC) methods in SPECT for maximum a posteriori reconstructions. While sophisticated reconstruction-based SC methods are available, the commonly used window-based SC methods are fast, easy to use, and perform reasonably well. Rather than subtracting a scatter estimate from the measured sinogram and then reconstructing, we use an ensemble approach and model the mean scatter sinogram in the likelihood function. This mean scatter sinogram estimate, computed from satellite window data, is itself inexact (noisy). Therefore two sources of noise, that due to Poisson noise of unscattered photons and that due to the model error in the scatter estimate, are propagated into the reconstruction. The optimization and comparison is driven by a figure of merit, the area under the LROC curve (ALROC) that gauges performance in a signal detection plus localization task. We use model observers to perform the task. This usually entails laborious generation of many sample reconstructions, but in this work, we instead develop a theoretical approach that allows one to rapidly compute ALROC given known information about the imaging system and the scatter correction scheme. A critical step in the theory approach is to predict additional (above that due to to the propagated Poisson noise of the primary photons) contributions to the reconstructed image covariance due to scatter (model error) noise. Simulations show that our theory method yields, for a range of search tolerances, LROC curves and ALROC values in close agreement to that obtained using model observer responses obtained from sample reconstruction methods. This opens the door to rapid comparison of different window-based SC methods and to optimizing the parameters (including window placement and size, scatter sinogram smoothing kernel) of the SC method.
There has been extensive work in the development of model observers , , , ,  to optimize or compare, via a scalar figure of merit (FOM), imaging systems and image reconstruction strategies in emission tomography. Such methods allow one to prototype time consuming human observer studies. The use of model observers usually entails sample methods involving many reconstructions. Since sample methods themselves can be time consuming, we shall use theoretical methods to compute a scalar FOM. Such theory methods can be run orders of magnitude faster than sample methods, and this is one motivation for the present work.
The problem we address is the optimization and comparison of window-based scatter correction (SC) methods in SPECT for MAP (maximum a posteriori) reconstruction. It turns out that this is equivalent to a problem of optimizing and comparing reconstruction strategies. We use a realistic task in which a model observer must detect as well as localize a signal (lesion).
The current work builds on, but is a significant extension of, two of our previous works , . In , we introduced theory methods to rapidly evaluate performance of an observer for the task of detecting and localizing a signal in a SPECT MAP reconstruction. This had nothing to do with SC. However, in , we addressed the problem of formulating theoretical expressions for the additional bias and covariance (above that due to photon noise) in a SPECT MAP reconstruction that included window-based SC. Here, we extend  and combine it with  to address a problem of practical import.
We will use certain notational conventions. Bold lower-case letters refer to lexicographically ordered vectors and calligraphic and non-calligraphic capital letters correspond to matrices. The corresponding notation indicating components of these quantities follow obviously. Thus x is a vector and matrix, with xn and mn are elements of these. Also, diag(x) is a diagonal matrix with diagonal elements given by xn. If vector x is random, then is its estimate, and its mean.
Here, we give a basic background on window-based SC methods. In SPECT, scatter can account for a significant fraction of the counts detected within the photopeak energy window. The corresponding sinogram thus contains scatter counts, and when reconstructed, this scatter can lead to a degradation of image quality in addition to that due to other sources such as photon noise. Despite the development of more sophisticated reconstruction-based SC methods , , , , , , , , energy-window SC is a simple and oft-used method , , , , , ,  in which counts collected in “satellite” energy windows are used to estimate the scatter contribution in the photopeak window (See Fig. 1). One could simply subtract such an estimated scatter sinogram (or a smoothed version of it) and reconstruct the remaining counts. An improvement  involves eschewing the subtraction step in favor of modeling scatter directly in reconstruction as discussed below. It is this latter window-based SC method that we use.
Figure 1 shows an energy spectrum with two peaks. The central peak of width Wm is the photopeak. This photopeak is bounded by two satellite windows of width Ws. The satellite windows are shown overlapping the photopeak window but they needn’t overlap. Generally, Ws Wm. The counts outside of the photopeak are mostly due to scatter and some counts in the photopeak are due to scatter. It is the scatter contribution to the photopeak that we would like to estimate. In the TEW method one uses an interpolation assumption tantamount to counting all counts below the trapezoidal region (below the dotted line) as scatter counts. Let gl and gr are measured sinogram counts in the left and right satellite windows. Define s as the true scatter contribution to the photopeak window, and ŝ as its estimate. Then the TEW “trapezoid strategy” yields a scatter estimate ŝ given by
Here B is a smoothing operator used to smooth the counts in the noisy satellite windows sinogram. For a single peak radiopharmaceuticals such as Tc99m, the same strategy applies but only a left window is needed.
Our task will be to detect and localize a signal that will be somewhere on a search grid. Figure 3 illustrates. Figure 3(a) shows one slice of a simple phantom (64 × 64 × 16) with a hot and cold cylinder. The light rectangle circumscribing the hot-spot is a 20 × 20 region comprising the search region. The signal can appear on any of those 400 locations. The goal is to correctly detect if a signal is present, and if so, to correctly localize it within a search tolerance. Figures 3(b) and 3(c) are reconstructions of the phantom including compensation for detector and collimator bur, attenuation and TEW SC. In this case the signal is just a 1 × 1 hot-spot. Figure 3(b) shows such a signal correctly detected and localized within the search grid of the phantom. Figure 3(c) shows a situation where a noise blob has been incorrectly identified as the signal. Note that in our analysis, at most one signal can appear.
The LROC curve figure 2 is useful in analyzing detection and localization performance. The observer views the image and reports a scalar response λ. (Note that the object and its reconstruction is 3D but detection and localization is typically carried out within a 2D slice.) The signal is deemed present/absent if λ τ, where τ is a decision threshold. Define PFP(τ) and PTP(τ) as the probabilities of a false positive report and true positive report with τ a threshold, respectively. The ROC curve is a plot of PFP(τ) vs. PTP(τ). In a LROC study, the observer reports a response λ and a location. PCL(τ) is the probability of deciding that the signal is present when the signal is actually present and also correctly localized. The LROC curve is a plot of PCL(τ) vs. PFP(τ). The area under the LROC curve, ALROC, is a figure of merit (FOM) useful in summarizing an LROC study. Note 0 ≤ ALROC ≤ 1, with ALROC = 1 indicating perfect performance.
Define f as the N-dim lexicographically ordered vector representing the object and g as the M -dim lexicographically ordered vector representing the sinogram data in the photo peak. Then, for SPECT, the imaging model is given by g ~ Poisson() where indicates the mean sinogram and g is distributed as an independent Poisson process at each bin. Importantly the mean data is given by
Here HP is the sparse system matrix accounting for the primary photons in the photopeak and HS is the dense system matrix accounting for the scatter photons in the photopeak. So here [HP]mn is proportional to the probability of a photon emanating from voxel n and being detected at the detector bin m. A similar definition for scatter photons obtains for [HS]mn. To accurately compute HSf one would need to resort to a photon tracking Monte Carlo (MC) simulator such as SIMSET or SIMIND . To compute HPf an analytical projector suffices. Inevitably in a reconstruction, forward projection and backprojection operations involving these system matrices are needed. It is completely impractical to use HS in projection and backprojection operations. There have been many studies in which approximate versions of HS are used ,  in reconstructions. However, we consider window-based methods that are different and much less computationally demanding. In fact, there is considerable evidence , , , , ,  that window-based SC can often perform as well as so called reconstruction based SC (RBSC) and indeed, window-based methods are in frequent use.
In window-based scatter one replaces HSf with an scatter estimate ŝ and treats ŝ as a constant independent of f. As mentioned above is computed directly form sinogram data. The imaging model effective replacing (2) becomes the affine model
Note that g is a noisy quantity by virtue of photon noise. However ŝ is itself derived from noisy quantities viz. (1). To the extent that ŝ is modeled as independent of f and to the extent that it differs from HPf, the scatter estimate ŝ is a form of “model error”, i.e. an error in the system model. A reconstruction that uses the imaging model (3) is thus a function of g and ŝ and hence is driven by two sources of noise: the Poisson noise in g and the noise in ŝ.
Define as the reconstruction of f, and the operator as the reconstruction operator. Then by the above argument = (g, ŝ; f). So far we have not committed to a particular reconstruction algorithm, but in this study we shall use MAP (maximum a posteriori) reconstruction. There are several reasons for this: (i) Unlike analytical methods MAP methods incorporate the imaging model in (3) (ii) MAP can model the Poisson noise in g (iii) In many cases the statistics of can be explicitly written for MAP reconstructions. We note that any iterative method, e.g OSEM, that incorporates an imaging model and an approximate noise model will give a reconstruction with image quality roughly equal to MAP, but MAP is much easier to analyze. Also the speed of MAP reconstructions is no longer an insurmountable barrier.
The MAP estimate is given by the following optimization
A reconstruction is then an optimization over f for a given instance of noisy ŝ and noisy g. Here Φ is the MAP objective, ΦL and ΦP are the log likelihood and log prior, respectively. The matrix is a digital analytic projector that approximates the “true” primary photon projection matrix HP. Note is usually carried out as a specialized function call. The last term in the (4) is a quadratic smoothing prior that penalizes voxel differences within a small neighbourhood (n) about n. Here wnn′ are weights proportional to the inverse distance between n and n′. The scalar β > 0 is weights the prior and controls the noise resolution trade-off in the reconstruction. Note that the prior can be written as the quadratic form
In (4) the scatter estimate ŝ has been incorporated into Φ. As an alternative one could have computed ŝ and simply subtracted it from the photopeak sinogram before reconstruction and indeed this is commonly done. However as shown in , , ,  one obtains improved results by directly modeling scatter in to to the reconstruction. In our work we use this strategy of modeling scatter directly into the reconstruction. We shall call this type of reconstruction an “SC-MAP reconstruction”.
Here, we consider how the noise in g and ŝ eventually propagated through and manifests itself as a mean and covariance of the SC-MAP reconstruction. We will need these first and second-order statistics of as components of our theoretical method for predicting ALROC of SC-MAP reconstruction.
As a first stage in the noise propagation chain, we consider the first and second order statistics of s. The mean scatter estimate and the covariance using (1) is given by
where l and r are the mean sinogram counts in the left and right satellite windows, respectively.
Given these moments, we may propagate them through to obtain the mean and covariance of the SC-MAP reconstruction. The MAP reconstruction is a function of g and ŝ, the noise in g and ŝ propagates into the reconstruction. By linearizing the reconstruction operator using a Taylor series approximation about the mean and the mean of , we derived  theoretical expressions for mean and covariance of the reconstruction. The mean of the MAP reconstruction given by is , that is, the mean SC-MAP reconstruction is well approximated as a reconstruction of the mean data.
The covariance of the SC-MAP reconstruction comprises terms involving the covariance of g, the covariance of and terms representing the cross-covariance of g, ŝ. In the case of overlapping windows let the counts in the part of the left satellite window that overlap the main window be gli and correspondingly gri be for the right overlap region. In the case of non-overlapping windows, g and ŝ will be uncorrelated. We derived the covariance of the MAP reconstruction in . It is given by
where is the Fisher information matrix. For the case of non-overlapping energy windows, the third and fourth terms of (10) are zero.
For a nonlinear shift-variant MAP reconstruction, resolution is characterized by a local point spread function that is object and position dependent. The lpsf expression we derived in  is given by
where ej is an unit impulse at location j in the object and Hl and Hr are the system matrices for photon propagation into the left and right satellite windows.
We have already described the LROC curve, but have given no information as to what the model observer is or how the observer response is computed. In this section, we describe an observer relevant to our task.
In detection + localization tasks using the LROC curve, one considers signal-present and -absent images. Let f− = b be the background (signal-absent object) and f+,j = b + sigj be the signal-present phantom with signal at j. Note that sigj is a known compact signal centered around its unknown location j. The signal sigj (if present) can appear somewhere on a signal grid with elements. Using our notation convention − and +,j are the reconstructions, and and are the mean reconstructions of f− = b and f+,j = b + sigj, respectively. Note that and . Consider a location loc from a set Ω of NΩ locations where one searches for the signal. The search region Ω is a dense grid comprising consecutive pixels as shown by the light rectangular region in Fig. 3(a). The signal grid is a sparser grid within the search grid on which the signal, if present, will appear. Figure 4 shows a 3 × 3 signal grid. In our ensuing discussion, we will need the concept of a mean reconstructed signal. It is given by
Note that from (12) the mean reconstructed signal is obtained by a difference of mean signal- present and absent reconstructions. Equation (12) shows the object-dependent resolution  property of MAP reconstruction.
After defining the constituent quantities required for the LROC task, we discuss the LROC model observer itself. For every location loc within the search grid Ω the model observer delivers a local observer response λloc. The observer then uses a maximization operator to localize the signal. It delivers a global observer response
The local observer response is given by
where wloc is the observer template applied to the background subtracted reconstruction to get the local observer response at loc for every location in the search grid Ω. The background subtraction becomes intuitively necessary when we consider a signal present in a cold region of the object. The max detector will yield a much higher response in a hot region of the object (where there is no signal) unless some form of background subtraction is included.
Our model observer is based on LROC model observer given by Gifford et al in , . Their observer tries to emulate humans and they have demonstrated good quantitative agreement with human observers. Psychophysical studies  indicate that the human visual system processes the retinal image by a feature extraction step in which the image is operated on by a series of bandpass channels to form NC features. Define is an N × NC matrix of channels centered at loc. Gifford  proposed a human emulating observer, the scanning Channelized Non-prewhitening observer (CNPW), that is given by
where is a grand average of the mean reconstructed signal averaged over the signal grid and centered at loc. The CNPW observer is similar to the well-known Channelized Hotelling Observer (CHO) ,  but without the prewhitening operation in the channel space.
In (15) we have defined the observer, but a few more constraints are needed to aid in explaining our theory calculations for ALROC. Given that − and +,j are the signal present and absent reconstructions with j being an element of the signal grid , using (15) we get the local observer responses λloc|− and λloc|(+, j) for every location loc in the search grid Ω. The local response λloc|(+, j) indicates that it is the response yielded by the observer at location loc when the reconstructed image has a signal present at j. The global observer response of the signal-absent reconstructions is obtained for (13) and is given by λ|− = maxlocΩλloc|−. The signal-present global observer response when the signal is present at j a location in the signal grid is given by λ|(+, j) = maxlocΩλloc|(+, j). Note that for the signal-present case when the signal is present at j the max can occur at any location and does not necessarily have to occur at j. For convenience we make further definitions so that we may write observer responses in a compact form. Denote the vector of all local observer responses λloc, loc = 1, …, NΩ a vector Λ. Define Λ|− as the response vector when there is no signal present and Λ|(+, j) as the response vector when the signal is present at j. From our previous definitions of λ|− and λ|(+, j), we may write λ|− = max Λ|− and λ|(+, j) = max Λ|(+, loc). We define W as a matrix whose NΩ columns comprise of wloc. Then from the definitions of Λ, Λ|−, Λ|(+, j), we get
Whether through theory methods or sample reconstruction methods, our net output will be a collection of global observer responses λ|− for signal-absent cases and λ|(+, j) for signal-present case with a signal at j S. To compute the ALROC - our final FOM - we do the following.
For a given tolerance distance one takes each of the λ|(+, j)’s and computes the distance from j to the location where the signal was deemed to be present. If this distance exceeds the tolerance distance then this observer response is eliminated from the set of signal-present observer responses. We form two histograms, one for the signal-absent responses and one for the remaining signal-present responses. We normalize the signal-absent histogram by the total number of signal-absent sample so that the area under the signal-absent histogram is unity. We normalize the signal-present histogram by the total number of signal-present samples (including those which were mislocalized i.e. exceeded the tolerance distance). This leads to a histogram area that is generally less than unity.
We integrate the area under each histogram using integration limits ranging from the decision threshold τ to +∞. The two areas correspond to PFP(τ) and PCL(τ) for the signal-absent and present cases, respectively. By sweeping τ we generate the LROC curve. We then use trapezoidal integration to calculate the ALROC.
The purpose of our paper is to show a theoretical procedure for computing ALROC that is much faster than that using a sample reconstruction method. However, we use the sample reconstruction method to validate the theory method. The accuracy of the of the sample reconstruction method is limited only by sample number, i.e. the number of reconstructions, while the accuracy of the theory method is limited by various theoretical approximations.
We generate Nsamp− signal-absent realization of g and ŝ, and from these obtain Nsamp− signal-absent reconstructions. For a given signal in the signal grid we generate Nsamp+ realizations of signal-present g and ŝ and from these obtain Nsamp+ signal-present reconstructions. Since there are signal locations we thus obtain a total of Nsamp+ signal-present reconstructions. We note that the this method has a implicit assumption that a signal can appear with equal probability at any location in , but it is trivial to modify this method if the prior probability Pr(j) of signal appearing at j is nonuniform.
For each signal-absent reconstruction we scan the image with the observer, take the max, and thus obtain the global observer response λ|−. For each signal-present reconstruction we do the same operation thus obtaining global signal-present observer responses λ|(+, j). Armed with the set of these signal-present and absent observer responses we can use the methods mentioned above to generate ALROC.
In this section we give a theoretical method to calculate the ALROC which supplants the slower sample reconstruction method explained earlier.
Our basic assumption is that Λ has a multi-variate Gaussian (Normal) distribution, i.e. Λ ~ (, ) under both the signal absent and present hypothesis. Here is the mean Λ and is its covariance. The assumption will hold if itself is normal and may hold approximately even if is not normal . If one could develop a theoretical expression for and then one need only sample this normal distribution to generate observer responses. This sampling is far faster than obtaining observer responses via sample reconstructions.
We now present  the mean and covariance of Λ for signal-absent and present cases.
We note that the physical effects of attenuation, blur and scatter are modeled in .
While we have expressions for (10) this is matrix is too large to deal with directly. Instead we use a locally stationary approximation ,  that applies to matrix vector products about the point of interest i. In this case the i, jth element of can be given by using the approximation
The cost in obtaining the local approximation is approximately one full backprojection per location i . Such an approximation is accurate and is quite efficient if there are few locations in the search grid. A SPECT-specific alternate method to calculate was presented in . This method, while complex, leads to local covariance approximations at speeds orders of magnitude faster than one backprojection per location.
The other theoretical expression we need is . This can be obtained using (12) applied to sigj rather than ej. One could instead obtain via the two-reconstruction method in , however, the theoretical method is far faster.
Once we have and we can use standard methods in ,  or using the MATLAB function mvnrnd to generate normal samples of signal-present and absent Λ’s. After a max operation and a tolerance test applied to each sample Λ we obtain samples of λ|− and λ|(+, j). With these quantities we follow the procedures in Sec. III-D to generate LROC.
Here we use the sample reconstruction method to validate the theory computation of ALROC.
As seen in figure 5 we used a 3D 64 × 64 ×16 object comprising 3 cylinders: small hot and cold cylinders in a medium intensity cylindrical background. The cylindrical background had a diameter of 37cm. We used a cubic voxel size of 0.625cm. Each camera face comprised 96×32 square detector bins of size 0.625cm. We simulated a circular-orbit parallel-beam geometry with 65 equispaced projection angles. We used SIMSET to simulate noiseless image formation through a parallel hole collimator with a depth dependent collimator response. The collimator bore radius and length was 0.1 cm and 3.0 cm, respectively, and the radius of rotation (measured from the center of the object to the collimator face) was 30 cm. These parameters were consistent with a depth-dependent resolution having a Gaussian psf whose σ was modeled as σ (cm) = 0.026 d (cm) + 0.0392 cm where d is depth (measured from the collimator face). We used a 140 Kev Tc99m source and H2O attenuator. The entire object had uniform attenuation. Counts were collected into a photopeak window between 128–152 Kev and the non-overlapping low energy satellite window was 124–128 Kev. The energy resolution of the detector was Gaussian with a FWHM of 10% for Tc99m. On average the total primary counts in the photopeak window was 4.4120 × 106, scatter counts in the photopeak window was 1.1680 × 106, primary counts in the satellite window was 8.4019 × 104 and total scatter counts in the satellite window was 3.2979 × 105.
We performed 3D reconstructions as described below, but of the 16 transaxial slices comprising the object, we used the 9th reconstructed slice to perform the detection and localization task. In the 9th slice we used the 20 × 20 search grid and 3 × 3 signal grid as shown in figure 4. The signal itself was a single voxel cubic signal. We performed validation runs at two signal contrasts, 10:1 and 18:1.
We performed MAP reconstructions using the COSEM MAP reconstruction  run to high iterations to ensure convergence. The initial reconstructed estimate was simply a noiseless reconstruction. To accomplish sample reconstructions, we use a sinogram g obtained by adding Poisson noise , and a TEW estimate ŝ obtained by using gl which is itself a Poisson noise realization of l. To form the estimate ŝ we did not implement a smoothing filter B. We used a smoothing parameter β = 0.005 for a reasonable noise resolution tradeoff.
We computed Nsamp−= 2000 signal-absent reconstructions. We computed 2700 signal-present reconstructions using Nsamp+ = 300 at each of the 9 locations in . This sample number leads to a very small sample error in ALROC.
Figure 6 summarizes our validation results. Sample quantities were computed as described in Sec. III-E and theory quantities were computed as described in Sec. III-F. Figure 6(a) shows LROC curves for the 10:1 signal contrast at three tolerances. Both theory and sample curves are shown but are hard to distinguish since they overlap almost exactly. In figure 6(b) the same situation obtains though with the higher contrast the LROC curves are higher. Figure 6(c) shows ALROC vs tolerance radius for theory and sample for the 10:1 signal contrast. Figure 6(d) shows the ALROC versus tolerance radius for the higher signal contrast of 18:1. In both cases the theory and sample results agree closely. Thus our theory method is validated.
We have validated our fast theory method for computing ALROC in TEW scatter-corrected SPECT MAP reconstruction. The theory method can be generalized to any window-based scatter correction method. One technical contribution of this paper has been in the developments of expressions for and when model error, attenuation and depth-dependent blur are included. This contribution extends the much simpler work in . As mentioned earlier the main computational bottle-neck for the theory method is in the computation of the local approximation to the covariance matrix . Using techniques in  to compute , the entire theory calculation can be done 3–4 orders of magnitude faster than sample based reconstruction based methods (assuming a number of samples adequate to obtain reasonable sample error). Armed with accurate ALROC theory estimates, we plan to compare a variety of window-based SC methods and to optimize them by, for example, optimizing sinogram smoothing (), reconstruction smoothing (β) or window parameters.
This work was supported by NIH NIBIB 02629.
Santosh Kulkarni, Electrical & Computer Engineering Department, Stony Brook University, Stony Brook NY.
Parmeshwar Khurd, Philips Corp. Cleveland OH.
Lili Zhou, Department of Radiology, Stony Brook University, Stony Brook NY, USA.
Gene Gindi, Department of Radiology, Stony Brook University, Stony Brook NY, USA.