Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2892877

Formats

Article sections

Authors

Related links

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.4436989PMCID: PMC2892877

NIHMSID: NIHMS210592

Santosh Kulkarni, Electrical & Computer Engineering Department, Stony Brook University, Stony Brook NY;

Contact: G. Gindi (telephone: 631-444-2539, Email: ude.bsynus.lim@idnig)

See other articles in PMC that cite the published article.

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 [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.

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 *x _{n}* and []

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.

Figure 1 shows an energy spectrum with two peaks. The central peak of width *W _{m}* is the photopeak. This photopeak is bounded by two satellite windows of width

$$\widehat{\mathbf{s}}=\frac{{W}_{m}}{2{W}_{s}}\mathbf{B}({\mathbf{g}}_{l}+{\mathbf{g}}_{r})$$

(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 Tc^{99}* ^{m}*, 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].

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.

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 *λ* *τ*, where *τ* is a decision threshold. Define *P _{FP}*(

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

$$\overline{\mathbf{g}}={\mathbf{H}}^{P}\mathbf{f}+{\mathbf{H}}^{S}\mathbf{f}$$

(2)

Here **H*** ^{P}* is the sparse system matrix accounting for the primary photons in the photopeak and

In window-based scatter one replaces **H**^{S}**f** 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

$$\overline{\mathbf{g}}={\mathbf{H}}^{P}\mathbf{f}+\widehat{\mathbf{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 **H**^{P}**f**, 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

$$\begin{array}{l}\widehat{\mathbf{f}}=\mathcal{O}\{\mathbf{g},\widehat{\mathbf{s}};\mathbf{f}\}=arg\underset{\mathbf{f}\ge \mathbf{0}}{max}\mathrm{\Phi}(\mathbf{f},\mathbf{g},\widehat{\mathbf{s}})\\ =arg\underset{\mathbf{f}\ge \mathbf{0}}{max}\phantom{\rule{0.16667em}{0ex}}({\mathrm{\Phi}}_{L}(\mathbf{g},\widehat{\mathbf{s}};\mathbf{f})+\beta {\mathrm{\Phi}}_{P}(\mathbf{f}))\\ =arg\underset{\mathbf{f}\ge \mathbf{0}}{max}\sum _{m}{g}_{m}log(\sum _{n}{\mathcal{H}}_{mn}{f}_{n}+{\widehat{s}}_{m})-\end{array}$$

(4)

$$\sum _{m}(\sum _{n}{\mathcal{H}}_{mn}{f}_{n}+{\widehat{s}}_{m})$$

(5)

$$-\frac{\beta}{2}\sum _{n}\sum _{{n}^{\prime}\in \mathcal{N}(n)}{w}_{n{n}^{\prime}}{({f}_{n}-{f}_{{n}^{\prime}})}^{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 Φ

$${\mathrm{\Phi}}_{P}(\mathbf{f})=-\frac{\beta}{2}{\mathbf{f}}^{T}\mathcal{R}\mathbf{f}$$

(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”.

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
$\overline{\widehat{\mathbf{s}}}$ and the covariance using (1) is given by

$$\overline{\widehat{\mathbf{s}}}=\frac{{W}_{m}}{2{W}_{s}}\mathbf{B}({\overline{\mathbf{g}}}_{l}+{\overline{\mathbf{g}}}_{r})$$

(8)

$${\mathcal{K}}_{\widehat{\mathbf{s}}}=\frac{{W}_{m}^{2}}{4{W}_{s}^{2}}\mathbf{B}\text{diag}\left(({\overline{\mathbf{g}}}_{l}+{\overline{\mathbf{g}}}_{r})\right){\mathbf{B}}^{T}$$

(9)

where * _{l}* and

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
$\overline{\widehat{\mathbf{s}}}$, we derived [7] theoretical expressions for mean and covariance of the reconstruction. The mean of the MAP reconstruction given by
$\overline{\widehat{\mathbf{f}}}$ is
$\overline{\widehat{\mathbf{f}}}\approx \mathcal{O}(\overline{\mathbf{g}},\overline{\widehat{\mathbf{s}}};\mathbf{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 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 **g*** _{li}* and correspondingly

$${\mathcal{K}}_{\widehat{\mathbf{f}}}={[\mathcal{F}+\beta \mathcal{R}]}^{-1}\mathcal{F}{[\mathcal{F}+\beta \mathcal{R}]}^{-1}+\frac{{W}_{m}^{2}}{4{W}_{s}^{2}}{[\mathcal{F}+\beta \mathcal{R}]}^{-1}{\mathcal{H}}^{T}\text{diag}\left(\frac{\overline{\mathbf{g}}}{{(\mathcal{H}\overline{\widehat{\mathbf{f}}}+\overline{\widehat{\mathbf{s}}})}^{2}}\right)\mathbf{B}\times \text{diag}({\overline{\mathbf{g}}}_{l}+{\overline{\mathbf{g}}}_{r}){\mathbf{B}}^{T}\text{diag}\left(\frac{\overline{\mathbf{g}}}{{(\mathcal{H}\overline{\widehat{\mathbf{f}}}+\overline{\widehat{\mathbf{s}}})}^{2}}\right)\mathcal{H}{[\mathcal{F}+\beta \mathcal{R}]}^{-1}-\frac{{W}_{m}}{2{W}_{s}}{[\mathcal{F}+\beta \mathcal{R}]}^{-1}{\mathcal{H}}^{T}\text{diag}\left(\frac{1}{\mathcal{H}\overline{\widehat{\mathbf{f}}}+\overline{\widehat{\mathbf{s}}}}\right)\text{diag}({\overline{\mathbf{g}}}_{li}+{\overline{\mathbf{g}}}_{ri}){\mathbf{B}}^{T}\text{diag}\left(\frac{\overline{\mathbf{g}}}{{(\mathcal{H}\overline{\widehat{\mathbf{f}}}+\overline{\widehat{\mathbf{s}}})}^{2}}\right)\mathcal{H}{[\mathcal{F}+\beta \mathcal{R}]}^{-1}-\frac{{W}_{m}}{2{W}_{s}}{[\mathcal{F}+\beta \mathcal{R}]}^{-1}{\mathcal{H}}^{T}\text{diag}\left(\frac{\overline{\mathbf{g}}}{{(\mathcal{H}\overline{\widehat{\mathbf{f}}}+\overline{\widehat{\mathbf{s}}})}^{2}}\right)\mathbf{B}\text{diag}({\overline{\mathbf{g}}}_{li}+{\overline{\mathbf{g}}}_{ri})\text{diag}\left(\frac{\mathbf{1}}{\mathcal{H}\overline{\widehat{\mathbf{f}}}+\overline{\widehat{\mathbf{s}}}}\right)\mathcal{H}{[\mathcal{F}+\beta \mathcal{R}]}^{-1}$$

(10)

where $\mathcal{F}\equiv {\mathcal{H}}^{T}\text{diag}\left({\scriptstyle \frac{\overline{\mathbf{g}}}{{(\mathcal{H}\overline{\widehat{\mathbf{f}}}+\overline{\widehat{\mathbf{s}}})}^{2}}}\right)\mathcal{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

$${\overline{\widehat{\mathbf{e}}}}_{j}\approx {(\mathcal{F}+\beta \mathcal{R})}^{-1}{\mathcal{H}}^{T}\text{diag}\left(\frac{1}{\mathcal{H}\overline{\widehat{\mathbf{f}}}+\overline{\widehat{\mathbf{s}}}}\right)\phantom{\rule{0.16667em}{0ex}}({\text{H}}^{P}+{\text{H}}^{S}){\mathbf{e}}_{j}-\frac{{W}_{m}}{2{W}_{s}}{(\mathcal{F}+\beta \mathcal{R})}^{-1}{\mathcal{H}}^{T}\text{diag}\left(\frac{\overline{\mathbf{g}}}{{(\mathcal{H}\overline{\widehat{\mathbf{f}}}+\overline{\widehat{\mathbf{s}}})}^{2}}\right)\times \mathbf{B}({\text{H}}^{l}+{\text{H}}^{r}){\mathbf{e}}_{j}$$

(11)

where **e*** _{j}* is an unit impulse at location

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}* =

$$\overline{\widehat{{\mathbf{sig}}_{j}}}={\overline{\widehat{\mathbf{f}}}}_{+,j}-\overline{\widehat{\mathbf{b}}}$$

(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

$$\lambda =\underset{\mathit{loc}\in \mathrm{\Omega}}{max}{\lambda}_{\mathit{loc}}$$

(13)

The local observer response is given by

$${\lambda}_{\mathit{loc}}={\mathbf{w}}_{\mathit{loc}}(\widehat{\mathbf{f}}-\overline{\widehat{\mathbf{b}}})$$

(14)

where **w*** _{loc}* is the observer template applied to the background subtracted reconstruction to get the local observer response at

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 *N _{C}* features. Define is an

$${\mathbf{w}}_{\mathit{loc}}={\mathcal{T}}_{\mathit{loc}}{\mathcal{T}}_{\mathit{loc}}{<\overline{\widehat{\mathbf{sig}}}>}_{\mathit{loc}}$$

(15)

where
${<\overline{\widehat{\mathbf{sig}}}>}_{\mathit{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 _{−} and _{+}* _{,j}* are the signal present and absent reconstructions with

$$\mathbf{\Lambda}={\mathbf{W}}^{T}(\widehat{\mathbf{f}}-\overline{\widehat{\mathbf{b}}})$$

(16)

$$\mathbf{\Lambda}\mid -={\mathbf{W}}^{T}({\widehat{\mathbf{f}}}_{-}-\overline{\widehat{\mathbf{b}}})$$

(17)

$$\mathbf{\Lambda}\mid (+,j)={\mathbf{W}}^{T}({\widehat{\mathbf{f}}}_{+,j}-\overline{\widehat{\mathbf{b}}})$$

(18)

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 *P _{FP}*(

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 *N _{samp}*

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 [6]. 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 [6] the mean and covariance of **Λ** for signal-absent and present cases.

$$\begin{array}{c}\overline{\mathbf{\Lambda}}\mid -=\mathbf{0}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\overline{\mathbf{\Lambda}}\mid (+,j)={\mathbf{W}}^{T}\overline{\widehat{{\mathbf{sig}}_{j}}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\forall j\in \mathcal{S}\\ {\mathcal{K}}_{\mathbf{\Lambda}\mid -}={\mathcal{K}}_{\mathbf{\Lambda}\mid (+,j)}={\mathbf{W}}^{T}{\mathcal{K}}_{\widehat{\mathbf{f}}}\mathbf{W}\end{array}$$

(19)

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
${\mathcal{K}}_{\widehat{\mathbf{f}}}^{i}$ [33], [34] that applies to matrix vector products
${\mathcal{K}}_{\widehat{\mathbf{f}}}^{i}\mathbf{x}$ about the point of interest *i*. In this case the *i, j ^{th}* element of can be given by using the approximation

$${[{\mathcal{K}}_{\mathbf{\Lambda}}]}_{ij}=\frac{1}{2}{\mathbf{w}}_{i}^{T}{\mathcal{K}}_{\widehat{\mathbf{f}}}^{i}{\mathbf{w}}_{j}+\frac{1}{2}{\mathbf{w}}_{i}^{T}{\mathcal{K}}_{\widehat{\mathbf{f}}}^{j}{\mathbf{w}}_{j}$$

(20)

The cost in obtaining the local approximation
${\mathcal{K}}_{\widehat{\mathbf{f}}}^{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
${\mathcal{K}}_{\widehat{\mathbf{f}}}^{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
$\overline{\widehat{{\mathbf{sig}}_{j}}}$. This can be obtained using (12) applied to **sig*** _{j}* rather than

Once we have and 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.

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 37*cm*. We used a cubic voxel size of 0.625*cm*. Each camera face comprised 96×32 square detector bins of size 0.625*cm*. 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 Tc^{99}* ^{m}* source and H

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 9* ^{th}* reconstructed slice to perform the detection and localization task. In the 9

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 , and a TEW estimate **ŝ** obtained by using **g*** _{l}* which is itself a Poisson noise realization of

We computed *N _{samp}*

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
$\overline{\widehat{{\mathbf{sig}}_{\mathit{loc}}}}$ 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
${\mathcal{K}}_{\widehat{\mathbf{f}}}^{i}$. Using techniques in [34] to compute
${\mathcal{K}}_{\widehat{\mathbf{f}}}^{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 (), 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.

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 ^{99}^{m}Tc 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.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |