PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Nucl Sci Symp Conf Rec (1997). Author manuscript; available in PMC 2010 June 28.
Published in final edited form as:
IEEE Nucl Sci Symp Conf Rec (1997). 2007; 5(4436989): 3986–3993.
doi:  10.1109/NSSMIC.2007.4436989
PMCID: PMC2892877
NIHMSID: NIHMS210592

Rapid Optimization of SPECT Scatter Correction Using Model LROC Observers

Abstract

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.

I. Introduction

There has been extensive work in the development of model observers [1], [2], [3], [4], [5] 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 [6], [7]. In [6], 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 [7], 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 [6] and combine it with [7] to address a problem of practical import.

II. Background

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 [An external file that holds a picture, illustration, etc.
Object name is nihms210592ig1.jpg]mn are elements of these. Also, diag(x) is a diagonal matrix with diagonal elements given by xn. If vector x is random, then x̂ is its estimate, and [x with macron] its mean.

A. Window Based Scatter Correction

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 [9], [10], [11], [12], [10], [13], [14], [15], energy-window SC is a simple and oft-used method [16], [17], [8], [18], [19], [20], [21] 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 [2] 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.

Fig. 1
Energy histogram of sinogram data illustrates window-based SC. Counts in the narrow (Ws) satellite windows are used to estimate the scatter contribution in the photopeak window (width Wm) [8]

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 [double less-than sign] 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

s^=Wm2WsB(gl+gr)
(1)

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.

The TEW is but one of a family of window-based SC methods. Our analysis will apply to any such window-based method. These include methods found in [16], [22], [18], [20], [17], [23], [24] and [15].

B. Detection and Localization using the LROC Curve

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.

Fig. 3
SPECT Reconstruction with SC (a) Shows a phantom with the search grid (light rectangle about hot spot). (b) Anecdotal signal-present recon with correct localization (shown by arrow). (c) Signal-present as in (b) but incorrectly localized (shown by arrow). ...

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 λ [greater, less] τ, 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.

Fig. 2
(a) The LROC and ROC curves. (See text for discussion)

III. Methods

A. Incorporating Scatter Correction into the Reconstruction

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(g) where g indicates the mean sinogram and g is distributed as an independent Poisson process at each bin. Importantly the mean data is given by

g¯=HPf+HSf
(2)

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 [25]. 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 [26], [14] in reconstructions. However, we consider window-based methods that are different and much less computationally demanding. In fact, there is considerable evidence [13], [23], [24], [15], [9], [27] 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 f is computed directly form sinogram data. The imaging model effective replacing (2) becomes the affine model

g¯=HPf+s^
(3)

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 f as the reconstruction of f, and the operator An external file that holds a picture, illustration, etc.
Object name is nihms210592ig2.jpg as the reconstruction operator. Then by the above argument f = An external file that holds a picture, illustration, etc.
Object name is nihms210592ig2.jpg (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 f 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

f^=O{g,s^;f}=argmaxf0Φ(f,g,s^)=argmaxf0(ΦL(g,s^;f)+βΦP(f))=argmaxf0mgmlog(nHmnfn+s^m)
(4)
m(nHmnfn+s^m)
(5)
β2nnN(n)wnn(fnfn)2
(6)

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 An external file that holds a picture, illustration, etc.
Object name is nihms210592ig3.jpg is a digital analytic projector that approximates the “true” primary photon projection matrix HP. Note An external file that holds a picture, illustration, etc.
Object name is nihms210592ig3.jpg 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 An external file that holds a picture, illustration, etc.
Object name is nihms210592ig4.jpg(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

ΦP(f)=β2fTRf
(7)

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 [28], [29], [2], [30] 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”.

B. Covariance and LPSF of SC-MAP Reconstructions

Here, we consider how the noise in g and ŝ eventually propagated through An external file that holds a picture, illustration, etc.
Object name is nihms210592ig2.jpg and manifests itself as a mean and covariance of the SC-MAP reconstruction. We will need these first and second-order statistics of f 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 s^¯ and the covariance An external file that holds a picture, illustration, etc.
Object name is nihms210592ig5.jpg using (1) is given by

s^¯=Wm2WsB(g¯l+g¯r)
(8)

Ks^=Wm24Ws2Bdiag((g¯l+g¯r))BT
(9)

where gl and gr are the mean sinogram counts in the left and right satellite windows, respectively.

Given these moments, we may propagate them through An external file that holds a picture, illustration, etc.
Object name is nihms210592ig2.jpg 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 g and the mean of s^¯, we derived [7] theoretical expressions for mean and covariance of the reconstruction. The mean of the MAP reconstruction f given by f^¯ is f^¯O(g¯,s^¯;f), 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 f 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 [7]. It is given by

Kf^=[F+βR]1F[F+βR]1+Wm24Ws2[F+βR]1HTdiag(g¯(Hf^¯+s^¯)2)B×diag(g¯l+g¯r)BTdiag(g¯(Hf^¯+s^¯)2)H[F+βR]1Wm2Ws[F+βR]1HTdiag(1Hf^¯+s^¯)diag(g¯li+g¯ri)BTdiag(g¯(Hf^¯+s^¯)2)H[F+βR]1Wm2Ws[F+βR]1HTdiag(g¯(Hf^¯+s^¯)2)Bdiag(g¯li+g¯ri)diag(1Hf^¯+s^¯)H[F+βR]1
(10)

where FHTdiag(g¯(Hf^¯+s^¯)2)H 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 [7] is given by

e^¯j(F+βR)1HTdiag(1Hf^¯+s^¯)(HP+HS)ejWm2Ws(F+βR)1HTdiag(g¯(Hf^¯+s^¯)2)×B(Hl+Hr)ej
(11)

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.

C. LROC Model Observer

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 An external file that holds a picture, illustration, etc.
Object name is nihms210592ig6.jpg with An external file that holds a picture, illustration, etc.
Object name is nihms210592ig7.jpg elements. Using our notation convention f [equivalent] b and f+,j are the reconstructions, and f^¯b^¯ and f^¯+,j are the mean reconstructions of f = b and f+,j = b + sigj, respectively. Note that f^+,jb+sigj^ and f^¯+,jb+sigj^¯. 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

Fig. 4
Phantom with rectangular search grid (light region over hot-spot) and the 3 × 3 signal grid.
sigj^¯=f^¯+,jb^¯
(12)

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 [31] 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

λ=maxlocΩλloc
(13)

The local observer response is given by

λloc=wloc(f^b^¯)
(14)

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 [3], [4]. Their observer tries to emulate humans and they have demonstrated good quantitative agreement with human observers. Psychophysical studies [32] 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 An external file that holds a picture, illustration, etc.
Object name is nihms210592ig8.jpg is an N × NC matrix of channels centered at loc. Gifford [4] proposed a human emulating observer, the scanning Channelized Non-prewhitening observer (CNPW), that is given by

wloc=TlocTloc<sig^¯>loc
(15)

where <sig^¯>loc 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) [1], [32] 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 f [equivalent] b and f+,j are the signal present and absent reconstructions with j being an element of the signal grid An external file that holds a picture, illustration, etc.
Object name is nihms210592ig6.jpg, 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[set membership]Ωλ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[set membership]Ωλ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

Λ=WT(f^b^¯)
(16)
Λ=WT(f^b^¯)
(17)
Λ(+,j)=WT(f^+,jb^¯)
(18)

D. ALROC Computation from Observer Responses

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 [set membership] 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.

E. ALROC with Sample Reconstruction Method

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 An external file that holds a picture, illustration, etc.
Object name is nihms210592ig6.jpg we generate Nsamp+ realizations of signal-present g and ŝ and from these obtain Nsamp+ signal-present reconstructions. Since there are An external file that holds a picture, illustration, etc.
Object name is nihms210592ig7.jpg signal locations we thus obtain a total of An external file that holds a picture, illustration, etc.
Object name is nihms210592ig7.jpg 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 An external file that holds a picture, illustration, etc.
Object name is nihms210592ig6.jpg, 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.

F. ALROC with Theory Computation

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. Λ ~ An external file that holds a picture, illustration, etc.
Object name is nihms210592ig4.jpg([Lambda with macron], An external file that holds a picture, illustration, etc.
Object name is nihms210592ig9.jpg) under both the signal absent and present hypothesis. Here [Lambda with macron] is the mean Λ and An external file that holds a picture, illustration, etc.
Object name is nihms210592ig9.jpg is its covariance. The assumption will hold if f itself is normal and may hold approximately even if f is not normal [6]. If one could develop a theoretical expression for [Lambda with macron] and An external file that holds a picture, illustration, etc.
Object name is nihms210592ig9.jpg 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 [6] the mean and covariance of Λ for signal-absent and present cases.

Λ¯=0Λ¯(+,j)=WTsigj^¯jSKΛ=KΛ(+,j)=WTKf^W
(19)

We note that the physical effects of attenuation, blur and scatter are modeled in An external file that holds a picture, illustration, etc.
Object name is nihms210592ig10.jpg.

While we have expressions for An external file that holds a picture, illustration, etc.
Object name is nihms210592ig10.jpg (10) this is matrix is too large to deal with directly. Instead we use a locally stationary approximation Kf^i [33], [34] that applies to matrix vector products Kf^ix about the point of interest i. In this case the i, jth element of An external file that holds a picture, illustration, etc.
Object name is nihms210592ig10.jpg can be given by using the approximation

[KΛ]ij=12wiTKf^iwj+12wiTKf^jwj
(20)

The cost in obtaining the local approximation Kf^i is approximately one full backprojection per location i [33]. 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 Kf^i was presented in [34]. 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 sigj^¯. This can be obtained using (12) applied to sigj rather than ej. One could instead obtain sigj^¯ via the two-reconstruction method in [7], however, the theoretical method is far faster.

Once we have [Lambda with macron] and An external file that holds a picture, illustration, etc.
Object name is nihms210592ig9.jpg we can use standard methods in [35], [6] 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.

IV. Validation

Here we use the sample reconstruction method to validate the theory computation of ALROC.

A. Experimental Setup

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.

Fig. 5
Simulation geometry with camera face in one position (not to scale). Object with background: hot: cold= 2:4:1.

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 [36] 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 g, and a TEW estimate ŝ obtained by using gl which is itself a Poisson noise realization of gl. 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 An external file that holds a picture, illustration, etc.
Object name is nihms210592ig6.jpg. This sample number leads to a very small sample error in ALROC.

B. Results

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.

Fig. 6
PCL(τ) vs PFP (τ) for signal contrast of (a) 10:1 and (b)18:1. ALROC vs search tolerance for signal contrast of (c)10:1 and (d) 18:1

V. Discussion and Summary

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 An external file that holds a picture, illustration, etc.
Object name is nihms210592ig10.jpg and sigloc^¯ when model error, attenuation and depth-dependent blur are included. This contribution extends the much simpler work in [6]. As mentioned earlier the main computational bottle-neck for the theory method is in the computation of the local approximation to the covariance matrix Kf^i. Using techniques in [34] to compute Kf^i, 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 (An external file that holds a picture, illustration, etc.
Object name is nihms210592ig11.jpg), reconstruction smoothing (β) or window parameters.

Acknowledgments

This work was supported by NIH NIBIB 02629.

Contributor Information

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.

References

1. Barrett HH, Myers KJ. Foundations of Image Science. John Wiley & Sons; 2004.
2. Farncombe T, Gifford HC, Narayanan M, Pretorius H, Frey E, King MA. Assessment of scatter compensation strategies for Ga-67 SPECT using numerical observers and human LROC studies. J Nuc Med. 2004 May;45(5):802–812. [PubMed]
3. Gifford H, Kinahan P, Lartizen C, King M. Evaluation of multiclass model observers in PET LROC studies. IEEE Trans Nuclear Science. 2007 Feb;54(1):116–123. [PMC free article] [PubMed]
4. Gifford HC, King MA, Pretorius PH, Wells RG. A comparison of human and model observers in multislice LROC studies. IEEE Trans Med Imaging. 2005 Feb;24(2):160–169. [PubMed]
5. Gilland SD, Tsui BMW, Qi Y, Gullberg GT. Comparision of channelized Hotelling and human observers in determining optimum OS-EM reconstruction parameters for myocardial SPECT. IEEE Transactions on Nuclear Science. 2006 June;53(3):1200–1204.
6. Khurd P, Gindi G. Fast LROC analysis of Bayesian reconstructed emission tomographic images using model observers. Phys Med Biol. 2005 March;50(11):1519–1532. [PMC free article] [PubMed]
7. Khurd P, Kulkarni S, Gindi G. In: Noo F, Kudo H, Zend L, editors. Noise propagation from scatter correction in SPECT MAP reconstruction; Proc 8th Intl Meeting on Fully Three-Dimensional Image Reconstruction in Radiology and Nuclear Medicine; Jul, 2005. pp. 97–100.
8. Ogawa K, Harata Y, Ichihara T, Kubo A, Hashimoto S. A practical method for position-dependent Compton-scatter correction in single photon emission CT. IEEE Trans Med Imaging. 1991 Sep;10(3):408–412. [PubMed]
9. Xiao J, de Wit TC, Staelens SG, Beekman FJ. Evaluation of 3D Monte Carlo-Based scatter correction for 99mTc cardiac perfusion SPECT. J Nuc Med. 2006;47:1662–1669. [PubMed]
10. Frey E, Tsui B. A new method for modeling the spatially-variant, object-dependent scatter response function in SPECT. Conf. Rec. IEEE Nuc. Sci. Sym. Med. Imaging Conf; Nov. 1996; pp. 1082–1086.
11. Kadrmas DJ, Frey EC, Karimi SS, Tsui BMW. Fast implementations of reconstruction-based scatter compensation in fully 3D SPECT image reconstruction. Phys Med Biol. 1998;43:857–873. [PMC free article] [PubMed]
12. Kadrmas DJ, Frey EC, Tsui BMW. Application of reconstruction-based scatter compensation to Thallium-201 SPECT: Implementations for reduced reconstruction noise. IEEE Trans Med Imaging. 1998 June;17(3):325–333. [PMC free article] [PubMed]
13. Dewaraja Y, Ljungberg M, Fessler J. 3-D Monte Carlo-based scatter compensation in quantitative I-131 SPECT reconstruction. IEEE Trans Nuclear Science. 2006 Feb;53:181–188. [PMC free article] [PubMed]
14. Beekman F, de Jong H, van Geloven S. Efficient fully 3-D iterative SPECT reconstruction with Monte Carlo-based scatter compensation. IEEE Trans Med Imaging. 2002 Aug;21(8):867–877. [PubMed]
15. Wernick MN, Aarsvold JN. Emission Tomography: The Fundamentals of PET and SPECT. Elsevier Academic Press; 2004.
16. Jaszczak RJ, Greer KL, Floyd CE., Jr improved spect quantification using compensation for scattered photons. J Nuc Med. 1984;25:893–900. [PubMed]
17. Wang X, Koral K. A regularized deconvolution-fitting method for Compton-scatter in SPECT. IEEE Trans Med Imaging. 1992 Sept;11(3):351–360. [PubMed]
18. Haynor DR, Kaplan MS, Miyaoka RS, Lewellen Thomas K. Multiwindow scatter correction techniques in single-photon imaging. Medical Physics. 1995 December;22(12):2015–2024. [PubMed]
19. King MA, Hademenos GJ, Glick SJ. a dual-photopeak window method for scatter correction. J Nuc Med. 1992;33:605–612. [PubMed]
20. Moore S, Kijewski M, Mueller S, Rybicki F, Zimmerman R. Evaluation of scatter compensation methods by their effects on parameter estimation from SPECT projections. Medical Physics. 2001 February;28(2):278–287. [PubMed]
21. Gagnon D, Todd-Pokropek A, Arsenaualt A, Dupras G. Introduction to holospectral imaging in nuclear medicine for scatter subtraction. IEEE Trans Med Imaging. 1989 Sept;8(3):245–250. [PubMed]
22. Ljungberg M, King MA, Hademenos GJ, Strand SE. Comparison of four scatter correction methods using Monte Carlo simulated source distributions. J Nuc Med. 1994;35:143–151. [PubMed]
23. Zaidi Habib, Koral Kenneth. Scatter modelling and compensation in emission tomography. Eur J Nucl Med Mol Imaging. 2004 May;31(5):761–781. [PubMed]
24. Zaidi Habib. Quantitative Analysis in Nuclear Medicine Imaging. Springer; 2006.
25. Ljungber M, Strand S, King M. Monte Carlo Calculations in Nuclear Medicine Applications in Diagnostic Imaging. Institute of Physics Publishing; Bristol, UK: 1998.
26. Zeng GL, Bai C, Gullberg GT. A projector/backprojector with slice-to-slice blurring for efficient three-dimensional scatter modeling. IEEE Trans Med Imaging. 1999;18:722–732. [PubMed]
27. Narayanan MV, Pretorius PH, Dahlberg ST, Leppo JA, Botkin N, Krasnow J, Berndt W, Frey E, King M. Evaluation of scatter compesation strategies and their impact on human detection performance Tc-99m myocardial perfusion imaging. IEEE Trans Nuclear Science. 2003;50(5):1522–1527.
28. King MA, deVries DJ, Pan TS, Pretorius P, Case J. An investigation of filtering of TEW scatter estimates used to compensate for scatter with ordered subset reconstructions. IEEE Trans Nuclear Science. 1997 June;44(3):1140–1145023.
29. Farncombe T, Gifford HC, Narayanan M, Frey E, King MA. A comparison of triple energy window scatter compensation methods for Ga-67 tumor imaging. J Nuc Med. 2002;43 Abstract 54P.
30. Beekman FJ, Kamphuis C, Frey EC. Scatter compensation methods in 3D iterative SPECT reconstruction: A simulation study. Phys Med Biol. 1997;42:1619–1632. [PubMed]
31. Fessler JA, Rogers WL. Spatial resolution properties of penalized-likelihood image reconstruction: Space-invariant tomographs. IEEE Trans Image Proc. 1996 Sep;5(9):1346–1358. [PubMed]
32. Myers KJ, Barrett HH. The addition of a channel mechanism to the ideal-observer model. J Opt Soc Am A. 1987;4:2447–2457. [PubMed]
33. Qi J, Leahy RM. Resolution and noise properties of MAP reconstruction for fully 3-D PET. IEEE Trans Med Imag. 2000 May;19(5):493–505. [PubMed]
34. Stayman JW, Fessler J. Efficient calculation of resolution and covariance for penalized-likelihood reconstruction in fully 3-D SPECT. IEEE Trans Med Imaging. 2004 Dec;23(12):1543–1556. [PubMed]
35. Golub GH, Van Loan CF. Matrix Computation. John Hopkins Univ. Press; 1989.
36. Hsiao IT, Rangarajan A, Gindi G. A new convergent MAP reconstruction algorithm for emission tomography using ordered subsets and separable surrogates. Conf. Rec. IEEE Int. Symp. Biomed. Imaging; July 2002; IEEE; pp. 409–412.