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 Trans Nucl Sci. Author manuscript; available in PMC 2010 May 3.
Published in final edited form as:
IEEE Trans Nucl Sci. 2003; 4(MP-3): 2516–2520.
doi:  10.1109/TNS.2005.851458
PMCID: PMC2862501
NIHMSID: NIHMS195647

Rapid Computation of LROC Figures of Merit Using Numerical Observers (for SPECT/PET Reconstruction)

Parmeshwar Khurd, Student Member, IEEE and Gene Gindi, Member, IEEE*

Abstract

The assessment of PET and SPECT image reconstructions by image quality metrics is typically time consuming, even if methods employing model observers and samples of reconstructions are used to replace human testing. We consider a detection task where the background is known exactly and the signal is known except for location. We develop theoretical formulae to rapidly evaluate two relevant figures of merit, the area under the LROC curve and the probability of correct localization. The formulae can accommodate different forms of model observer. The theory hinges on the fact that we are able to rapidly compute the mean and covariance of the reconstruction. For four forms of model observer, the theoretical expressions are validated by Monte Carlo studies for the case of MAP (maximum a posteriori) reconstruction. The theory method affords a 102 – 103 speedup relative to methods in which model observers are applied to sample reconstructions.

Index Terms: Image Quality, SPECT, PET, MAP reconstruction, model observer, LROC

I. Introduction

Image Quality in emission computed tomography (ECT) reconstruction comprising PET or SPECT, is measured by how well the reconstruction supports a task, such as lesion detection by human or model observers. A figure of merit (FOM) summarizing detection task performance could be used, for instance, to optimize the parameters of a reconstruction algorithm. Such FOM’s can be obtained via human observer studies, but these quickly become impractical if many parameters are to be considered. One can substitute a model observer chosen to emulate human performance [1], and apply this observer to a test set of sample reconstructions. Even in this case, the computational requirements in generating a sufficient number of sample reconstructions can be daunting. A third alternative is to develop theoretical formulae that rely only on computable first and second-order reconstruction statistics to rapidly compute the desired FOM’s. In this paper, we focus on the third alternative.

The key to using theoretical methods is to have available rapidly computable expressions for the reconstruction covariance and mean. For linear reconstructions such as FBP, these are easy to compute, and expressions have also been developed for EM-ML [2] and MAP reconstructions [3]–[9]. Previous work has used these theoretical tools to calculate FOM’s for the case of SKE/BKE (signal known exactly/background known exactly) detection tasks. Such SKE/BKE FOM’s have been developed for the Hotelling and non-prewhitening observers [10] as applied to FBP and MAP reconstructions, and for the CHO (Channelized Hotelling) observer as applied to MAP reconstructions [11], [8], [12].

The SKE/BKE case does not include a search task where the signal location is unknown. In [13], [14], the model was extended to include forms of signal and background variability, but not signal location variability. Human performance is affected by a search task, thus it is desirable to include signal location variability in our theoretical computations of an appropriate FOM. While the addition of a search task can affect absolute human performance, it often leaves relative performance unaffected, e.g. the rank order of a performance FOM vs. some imaging parameter is preserved when search is added [15]. Our own interest is in the effect of the regularization parameter (i.e. the amount of smoothing) in an ECT MAP reconstruction on lesion detectability. Previous work with human [16] and model [8] observers for SPECT has shown an SKE/BKE FOM to be rather insensitive to regularization. It was suggested in [12] that addition of a search mechanism may affect this insensitivity to regularization for ECT.

In our work, we are given a reconstruction method, a model observer, a fixed background and a signal known except for location. Here, the model observer is one which presumably emulates human performance for the task at hand. Under this context, we show how to rapidly compute relevant FOM’s. The problem of LROC FOM computation via model observers using sample methods was addressed by Gifford et al. in [17] and [18]. They also studied correlation with human performance. The fast theoretical methods proposed here were inspired by the work in [17] and [18]. To explain our methods, we first review the ROC and LROC curves.

II. ROC and LROC curves

We consider model observers that compute an observer response λ when applied to the reconstruction. The response λ is compared to a decision threshold τ and we decide signal present if λτ and signal absent if λ < τ. Let “+” or “−” designate a signal present or absent in the reconstruction. We define λ|+ as the response for signal present and λ|− for signal absent. The ROC curve is generated as a plot of the true-positive (TP) rate vs. the false positive rate (FP) rate, where these are defined as PTP (τ) = Pr(λ|+ > τ) and PFP (τ) = Pr(λ|− > τ). The area under the resulting curve, AROC, is a widely accepted measure of detectability. We shall make use of the common acronyms pdf and cdf to designate “probability density function” and “cumulative distribution function”, respectively. If we denote Pλ|− (τ) = 1 − Pr(λ|− > τ) as the cdf of λ|− and pλ|+ (τ), Pλ|+ (τ) as the pdf and cdf of λ|+, then AROC may be computed [19] as

AROC=pλ+(τ)Pλ(τ)dτ.
(1)

The LROC curve is more subtle. The reconstructed images presented to the observer may either contain no signal or may have a signal present at any one of Nloc locations. As before, the observer computes a response λ and determines the presence or absence of a signal by comparison with a threshold τ. In addition to this, the observer also identifies a location where he deems the signal to be present. The abscissa of the LROC curve, the false-positive rate PFP (τ) = Pr(λ|− > τ), is defined as before. The ordinate of the LROC curve is given by PCL (τ), the probability of detecting a signal when one is present and also of correctly localizing the signal. In Sec. III, we give expressions for PCL (τ).

Fig. 1 shows examples of both the ROC curve and the LROC curve. The quantity PCL = PCL (τ = −∞) is the probability of correct localization given that the observer always decides “target-present” when the image indeed contains a target. Just as AROC is a FOM obtained from the ROC curve, the area under the LROC curve ALROC, and PCL are relevant FOM’s for the detection task with signal location uncertainty.

Fig. 1
The ROC and LROC curves

III. Theory for Rapid FOM Computation

A. Closed Form Expressions for LROC FOM’s

We will now obtain closed-form expressions for the LROC FOM’s for the case when the observer’s detection and localization strategy takes the form of a “max detector” explained below. In the remainder of the paper, we shall use, when dealing with a continuous random variable X, the notations pX (x) and PX (x) for the pdf and the cdf respectively, where x is the point at which p or P is evaluated. Note that PCL (τ), PTP (τ) and PFP (τ) are exceptions to this notational rule since they are not cdf’s. Consider a location loc from a set Ω of Nloc locations, where the signal may appear with probability Pr (loc). Also, define Ωloc as the set of all locations in Ω excluding loc. Define observer response λloc to be the local observer response at location loc, and take the “global” observer response as λ = maxloc λloc. The location reported by the observer is the location where the maximum occurs. (This form of a max detector for λ is commonly used for detection tasks with location uncertainty. One could consider other forms, but the max detector formulation naturally fits in with our work.)

Define λloc|− as the observer response at loc for an image that contains no signal. Define λloc|(+, loc′) as the local observer response at loc when the signal is located at loc′. When loc = loc′, the quantity λloc|(+, loc) can be interpreted as the local observer response at the signal location. Also, we define λ|− as the global observer response for an image that contains no signal, and λ|(+, loc) as the global observer response when there is a signal present and located at loc. It follows then that λ|− = maxloc[set membership]Ω λloc|− and that λ|(+, loc) = maxloc[set membership]Ω λloc |(+, loc). Also, we will define [lambda with tilde]|(+, loc) = maxloc[set membership]Ωloc λloc |(+, loc), i.e. [lambda with tilde]|(+, loc) is the maximum of the local observer responses at all non-signal locations when there is a signal present at loc. Then λ|(+,loc) = max {λloc|(+, loc), [lambda with tilde]|(+, loc)}.

We will make the following assumptions:

  1. We assume that for a given loc [set membership] Ω and ∀ loc[set membership] Ωloc, λloc|− and λloc|(+, loc′) have the same distributions. In other words, the distribution of a local observer response at a non-signal location depends upon the location itself, but not upon whether there is a signal present elsewhere or whether there is no signal present.
  2. The Nloc random variables λloc|−, loc [set membership] Ω are assumed to be jointly independent. The Nloc random variables λloc |(+, loc), loc[set membership] Ω are also assumed to be jointly independent ∀ loc [set membership] Ω. It then follows that the random variables λloc |(+, loc) and [lambda with tilde]|(+, loc) are independent ∀ loc [set membership] Ω.
  3. We also assume that Nloc is large. This is discussed further below.

Define λs as the signal response, i.e. the response at any signal location. The signal response is not the same as λloc |(+, loc). The latter term refers to a signal response at a specific location loc. The former term may correspond to any loc [set membership] Ω and hence it has the mixture pdf

pλs(λ)=locΩpλloc(+,loc)(λ)Pr(loc).
(2)

That is, if one took an ensemble of signal-present images with signals located according to Pr(loc), and computed the histogram of local observer responses at the signal location, then the histogram (suitably normalized) would be given by pλs (λ).

By Assumption 2 and using the fact that the cdf of the maximum of independent random variables is the product of the cdf’s of each individual random variable, we get

Pλ(λ)=locΩPλloc(λ).
(3)

By Assumptions 2 and 1,

Pλ(+,loc)(λ)=locΩlocPλloc(+,loc)(λ)=locΩlocPλloc(λ)locΩ.
(4)

Now, we define PCL (τ)|(+, loc) as the probability of correct detection and localization when the signal is present at location loc and the decision threshold τ is used. Note that to correctly detect and localize a signal, the response at the signal location must be greater than τ, and also greater than the maximum of the responses at all other locations. Therefore,

PCL(τ)(+,loc)=Pr(λloc(+,loc)>τ,λloc(+,loc)>λ(+,loc)).
(5)

Furthermore, since the r.v.’s λloc |(+, loc) and [lambda with tilde]|(+, loc) are independent by Assumption 2, it follows that

Pr(λloc(+,loc)>τ,λloc(+,loc)>λ(+,loc))=(τpλloc(+,loc)(λ)Pλ(+,loc)(λ)dλ).
(6)

Then, by averaging over location, PCL (τ) becomes

PCL(τ)=locΩPCL(τ)(+,loc)Pr(loc)=locΩ(τpλloc(+,loc)(λ)Pλ(+,loc)(λ)dλ)Pr(loc).
(7)

To simplify this integral, we will make the approximation

Pλ(+,loc)(λ)locΩPλloc(λ)=Pλ(λ).
(8)

We will assume that Nloc is large enough so that this approximation is fairly accurate when used to evaluate PCL (τ). (The test of this and other assumptions is supported by our final validation results in Sec. V.) Then from (7), (8) and (2)

PCL(τ)locΩ(τpλloc(+,loc)(λ)Pλ(λ)dλ)Pr(loc)=τPλ(λ)(locΩpλloc(+,loc)(λ)Pr(loc))dλ=τpλs(λ)Pλ(λ)dλ.
(9)

By substituting τ = −∞ in (9), we obtain

PCL=pλs(λ)Pλ(λ)dλ.
(10)

Using (9), we may obtain an expression for ALROC

ALROC=0PCL(1PFP(τ))dPCL(τ)=Pλ(τ)(pλs(τ)Pλ(τ))dτ=pλs(τ)(Pλ(τ))2dτ
(11)

which is the LROC analog to (1).

Our observer model was inspired by work in [20], [21], though it differs in many respects. In [20], it is assumed that all non-signal local observer responses were independent and identically distributed, whereas our Assumption 2 assumes independence of all local observer responses, but does not assume that they are identically distributed. Our Assumption 1 does not assume that all responses at non-signal locations are identically distributed, and allows these non-signal observer responses to be location dependent. However, the model in [21] is more general than ours and the conditions in [21] are satisfied with our assumptions. Equations (9)(11) may be directly obtained using the results in [21].

B. Using Reconstruction Statistics in LROC FOM Computation

Equations (10)(11) are the FOM’s we seek to rapidly evaluate. In what follows, we use additional assumptions along with our knowledge of reconstruction mean and covariance to derive efficient expressions for the terms pλs (λ) and Pλ|− (λ) that appear in (10)–(11).

Consider an N-dim lexicographically ordered object vector f, a (fixed) N-dim background vector b and an N-dim signal vector s. Note that s is known except for location and we use the symbol sloc to denote the signal centered at location loc. We model image formation by a simple Poisson model g ~ Poisson(Hf), where H is the system matrix whose element Hmn is proportional to the probability of receiving a count in detector bin m that emanates from pixel n, and g is the data vector (sinogram) whose elements are gm, m = 1, ···, M. For a signal-present object, we may append subscripts to write f+,loc = b + sloc and for a signal-absent object, f = b. We append a carat superscript to quantities to denote their reconstructions so that f = b is the reconstruction of the background and f^+,loc=b+sloc^ is the reconstruction of a signal-present object. We shall use the symbol f as a notational convenience to denote either b, or b+sloc^ as needed. We denote the reconstruction covariance by the N × N matrix Kf and note that Kf is approximately the same for both the signal-present and signal-absent cases. This approximation will be fairly accurate for a small signal and is commonly made in SKE/BKE studies [10], [8]. We also append a bar superscript to indicate the relevant reconstruction averaged over data noise, so that e.g. f^¯ is the mean signal-absent reconstruction.

We shall make repeated use of the approximation that f^¯ is well estimated by a noiseless reconstruction [3]. We define

s^¯locb+sloc^¯b^¯.
(12)

Note that in the above definition for s^¯loc, we have deviated from our usual policy of using a carat to denote a reconstructed quantity. Due to the fact that the reconstruction operator may be nonlinear, s^¯loc no longer represents the mean reconstruction of sloc itself, but rather the residual of signal-present and signal-absent mean reconstructions.

By assumption, we take our local observer to be the affine model observer λloc=wlocTf^+vloc, where wloc is an N-dim observer template and vloc is a scalar, and as before, λ = maxlocλloc. The scalar vloc is added to the familiar inner product wlocTf^ for reasons that will later become clear. In the ensuing discussion, we will use the convenient definitions λn|loc [equivalent] λloc| − and λs|loc [equivalent] λloc|(+, loc). From our definitions in Sec. II, we get

λ=maxloc[wlocTf^+vloc]=maxlocλnloc
(13)
λsloc=wlocTf^+,loc+vloc.
(14)

For the case of these affine local observers, we present some conditions under which the joint independence assumption (Assumption 2) in the previous section may hold. For all the local observer responses to be jointly independent, they need to be pairwise independent. If the local covariance of f (about location loc) and the template wloc have low amplitude tails, we may argue that λn|loc1 and λn|loc2 are independent if loc1 and loc2 are further apart than the signal width.

We make a second assumption that λn|loc is normally distributed as N(μnloc,σnloc2). This normality assumption holds exactly or approximately in many cases of interest. The local observer λloc has a mathematical form similar to observers for SKE/BKE detection tasks, which are commonly found to be normally distributed [19]. With these assumptions, (3) becomes

Pλ(τ)=locΩΦ(τμnlocσnloc).
(15)

Here, Φ denotes the cdf of the standard normal random variable. Using (13), we obtain

μnloc=wlocTb^¯+vlocσnloc2=wlocTKf^wloc.
(16)

We may claim that λs|loc is also normally distributed as N(μsloc,σsloc2). Using (14) and (12), we obtain

μsloc=wlocT(b^¯+s^¯loc)+vlocσsloc2=wlocTKf^wloc.
(17)

Using (2) and our normality assumption, we may write

pλs(τ)=locΩ1σslocφ(τμslocσsloc)Pr(loc)
(18)

where [var phi] denotes the standard normal pdf.

Thus, our goal, to evaluate equations (10)(11) by numerical 1-D integration, is attainable if we can compute the quantities μn|loc, σn|loc, μs|loc, and σs|loc These four quantities depend on the choice of observer (wloc, vloc).

We now consider a few common local observers. Each of the template expressions in (19)–(22) makes use of the definition in (12). The local non-prewhitening (NPW) observer corresponds to

wloc,NPW=(b+sloc^¯)b^¯=s^¯loc.
(19)

The local Hotelling (HO) observer corresponds to

wloc,HO=Kf^1s^¯loc.
(20)

Psychophysical studies [22] indicate that the human visual system processes the retinal image by a feature reduction step in which the image (here f) is operated on by a series of bandpass channels tk to form NC features. Define Tloc as a N × Nc matrix of the channels. The sub-script loc indicates that the channels have been centered at loc. Channelized versions of the local NPW and HO observers are given by

wloc,CNPW=TlocTlocTs^¯loc
(21)

wloc,CHO=Tloc(TlocTKf^Tloc)1TlocTs^¯loc
(22)

respectively.

We have defined the templates wloc, but not the offset term vloc. Here, we assume that vloc accounts for background subtraction, as recommended in [17], and is given by vloc=wlocTb^¯. In this case, λloc=wlocT(f^b^), showing the local observer to operate on a background subtracted image. The need for background subtraction becomes evident if 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.

The means and variances in (16)–(17) may be re-stated using the particular forms for our four observers given in (19)–(22). From (16), our form of background subtraction implies that μn|loc = 0 for all four local observers. For the NPW observer, we combine (19) and (16)–(17) to get

μsloc,NPW=s^¯locTs^¯locσnloc,NPW2=σsloc,NPW2=s^¯locTKf^s^¯loc.
(23)

For the HO observer, we combine (20) and (16)–(17) to get

μsloc,HO=σnloc,HO2=σsloc,HO2=s^¯locTKf^1s^¯loc.
(24)

For the CNPW observer, we combine (21) and (16)–(17) to get

μsloc,CNPW=s^¯locTTlocTlocTs^¯locσnloc,CNPW2=σsloc,CNPW2=s^¯locTTlocKf^TlocTs^¯loc.
(25)

For the CHO observer, we combine (22) and (16)–(17) to get

μsloc,CHO=σnloc,CHO2=σsloc,CHO2=s^¯locTTloc[TlocTKf^Tloc]1TlocTs^¯loc.
(26)

Note that in all the above expressions, b^¯ drops out, and only s^¯loc and Kf need be considered.

IV. Fast LROC FOM’s for MAP reconstruction

Our theory so far will apply to any scenario in which s^¯loc and Kf can be rapidly calculated. We still need to specify a particular reconstruction algorithm. We are interested in MAP reconstructions given by f = argmaxf0 ΨL(g|f) + βΨP(f), where ΨL is a Poisson log-likelihood and ΨP is a smoothing regularizer with weight β. We will assume a quadratic regularizer, given by: ΨP(f)=nnN(n)wnn(fnfn)2=12fTRf. The term N(n) is a local neighborhood about n. The weights are symmetrical (wnn = wnn) and positive. For a 2-D problem, a typical neighborhood N(n) comprises the eight nearest neighbors of n. The neighborhood structure is typically invariant with position (except for edge effects), so that R is an (approximate) doubly block circulant matrix [4].

For this class of reconstruction algorithms, it turns out [3] that Kf is approximated by

Kf^(F+βR)1F(F+βR)1
(27)

where F denotes the Fisher information matrix, given by F=HTdiag(1g¯)H. The other quantity we need is s^¯loc, which may be approximated using an expression for the local impulse response [23]. The local impulse response e^¯i due to ei, a point source at location i, is given by

e^¯i(F+βR)1Fei.
(28)

Then for the small, compact signal sloc, our approximation is given by

s^¯loc(F+βR)1Fsloc.
(29)

The N × N matrices F and R are formidably large and calculation of the inverses in (27) and (29) appears intractable. However, in some cases [5], we can use the fact that doubly block circulant matrices may be diagonalized by a 2D DFT (Discrete Fourier Transform) to obtain efficient expressions for (27) and (29). The matrix R is doubly block circulant as mentioned earlier. For a space-invariant imaging system, the quantity HTHei, equivalent to a projection followed by backprojection of a point source, yields a point response independent of i. For such a system, HTH is also doubly block circulant. In this case, it can be shown [5] that (29) can be approximated as

s^¯locQ1diag(ζn(ζn+βκloc2ηn))Qsloc.
(30)

Here, An external file that holds a picture, illustration, etc.
Object name is nihms195647ig1.jpg denotes the 2-D DFT. The vector ζ comprises the eigenvalues of HTH and the vector η comprises the eigenvalues of R. The quantity κ is calculated according to the formula: κn2=mHmn2/gm¯mHmn2. It turns out [5] that (27) may be similarly approximated as

Kf^κloc2Q1diag(ζn(ζn+βκloc2ηn)2)Q.
(31)

For a detailed discussion of the quantities κ, ζ and η, please see [5], [6], [8]. The approximation in (31) is one that applies locally at location loc, that is, the approximation is accurate only when computing matrix-vector products of the form Kfaloc, where aloc is a compact vector centered at loc.

Equations (30) and (31) can be substituted in the expressions for the observer templates (19)–(22), and the local means and variances, (23)–(26). The resulting expressions (32)–(35), crucial to our calculation of ALROC and PCL, can be evaluated extremely rapidly. In particular, the vectors κ, ζ and η may be pre-computed for a given object and system geometry. The dominant computation is that of κ which takes one backprojection for a new object if the system geometry remains fixed. In fact, we need only compute κ at Nloc locations, so that even a full backprojection is not needed. The remainder of the calculations for the means and variances take on the order of a 2D-FFT

For the NPW observer, we obtain

wloc,NPW=Q1diag(ζn(ζn+βκloc2ηn))Qslocμsloc,NPW=slocTQ1diag(ζn2(ζn+βκloc2ηn)2)Qslocσnloc,NPW2=σsloc,NPW2=κloc2slocTQ1diag(ζn3(ζn+βκloc2ηn)4)Qsloc
(32)

and for the HO observer, we obtain

wloc,HO=κloc2Q1diag(ζnβκloc2ηn)Qslocσnloc,HO2=σsloc,HO2=μsloc,HO=κloc2slocTQ1diag(ζn)Qsloc.
(33)

We can also develop similar expressions for the channelized observers in (21) and (22). To do this, we denote the 2D DFT of the channel matrix Tloc by Tloc and we use the superscript H to denote the Hermitian transpose. Then, for the channelized NPW observer, we obtain

wloc,CNPW=Q1TlocTlocHdiag(ζn(ζn+βκloc2ηn))Qslocμsloc,CNPW=slocTQ1diag(ζn(ζn+βκloc2ηn))TlocTlocHdiag(ζn(ζn+βκloc2ηn))Qslocσnloc,CNPW2=σsloc,CNPW2=κloc2slocTQ1diag(ζn(ζn+βκloc2ηn))Tlocdiag(ζn(ζn+βκloc2ηn)2)TlocHdiag(ζn(ζn+βκloc2ηn))Qsloc.
(34)

For the channelized Hotelling observer, we obtain

wloc,CHO=κloc2Q1Tloc[TlocHdiag(ζn(ζn+βκloc2ηn)2)Tloc]1TlocHdiag(ζn(ζn+βκloc2ηn))Qslocσnloc,CHO2=σsloc,CHO2=μsloc,CHO=κloc2slocTQ1diag(ζn(ζn+βκloc2ηn))Tloc[TlocHdiag(ζn(ζn+βκloc2ηn)2)Tloc]1TlocHdiag(ζn(ζn+βκloc2ηn))Qsloc.
(35)

This completes the theory needed to evaluate (10)–(11) for an affine observer. For example, for the NPW observer, in order to compute ALROC and PCL, we would first compute the terms in (32) and substitute them into the expressions (15) and (18), which may in turn be substituted into (10)–(11). Finally, (10)–(11) are evaluated by 1-D numerical integration.

V. Monte Carlo Validation of Theory Expressions

We validate, for the four observers, our theoretical expressions (10)–(11) using Monte Carlo (MC) methods. Our fixed background b (Fig. 2a, 2c) was a large 64 × 64 Gaussian lump, scaled to produce 1M counts over 64 angles. We chose this Gaussian lump in order to get some degree of background variation. The signal, 5 pixels in the shape of a cross, was placed on one of 25 possible locations. The array of locations constituted a square grid with a 6-pixel separation. Figures 2b, 2d show the signal placed on the background. We simulated an idealized ECT system ignoring degrading factors such as attenuation, randoms, scatter or depth-dependent blur. Activity was summed along parallel strips and Poisson noise added. Thus, there was no distinction of PET or SPECT in this simple simulation. Noisy MAP reconstructions were done for β = 0.1 and signal contrasts of (0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7). The value of β was chosen to yield signal-present reconstructions that were neither too noisy nor too smooth. Signal contrasts were defined as (amplitude of s)/(peak of b). The phrase “peak of b” refers to the unique overall peak of the unimodal background. Figure 3 shows two anecdotal signal-present reconstructions at contrasts 0.2 and 0.6. We used a positivity-constrained regularized ordered-subset convergent MAP algorithm [24] for the reconstruction.

Fig. 2
(a) The 64×64 background (b) Background+sign al at location (33, 27) for signal contrast 0.4 (c) Background profile (d) Back ground+signal profile. Note (a) and (b) use different grey scales for clarity.
Fig. 3
(a) Anecdotal reconstruction of background+signal for signal contrast 0.2 (b) Anecdotal reconstruction of background+signal at same location for signal contrast 0.6. Arrows show signal location.

500 signal-absent reconstructions were done. For each contrast, 20 signal-present reconstructions were done at each of the 25 locations. (Pr(loc) was uniform.) Thus, for each of the 7 signal contrasts, 500 reconstructions were done. Observer templates for the NPW, HO, CNPW and CHO observers were computed (32)–(35). (We note that we used theoretical formulae to compute the observer templates for MC experiments so that both MC and theory computations with (10)–(11) would utilize exactly the same template.) For the channelized observers, we used 4 square dyadic channels with a starting cutoff frequency of 0.041 cycles per pixel.

Applications of each observer resulted in a set of 500 values for λ|−. Also, for each observer, 500 values of λ|+ were obtained at each signal contrast. If a value of λ|+ corresponded to an incorrect location, it was not counted. In order to plot the LROC curve for a given signal contrast, we determined a set of 200 thresholds from the range of values of λ|+ and λ|−. For each threshold τ, the probability of correct localization PCL(τ) and the false positive rate PFP(τ) were determined and the LROC curve was plotted. Then ALROC was obtained by simple trapezoidal integration of the LROC curve, while PCL could be read off the curve.

Figure 4 shows the MC vs. theory results for ALROC. (Results for PCL are not shown.) The theoretical ALROC values for each signal contrast and observer were computed using the formulae in Secs. III and IV. Error bars on the MC estimates of ALROC for each observer and each signal contrast were obtained by the following procedure [25]: Using the bootstrap, 1000 re-sampled sets of 500 values for λ|− and λ|+ were obtained. For each set, the values of ALROC were obtained using the MC procedure described earlier. We then determined 95% confidence intervals from histograms of these 1000 values.

Fig. 4
(a) ALROC vs. signal contrast for NPW observer (b) ALROC vs. signal contrast for Hotelling observer (c) ALROC vs. signal contrast for CHO observer (d) ALROC vs. signal contrast for CHO observer. Error bars indicate 95% confidence intervals.

VI. Discussion

We have proposed an efficient means to calculate, via theory, LROC figures of merit for the case of a BKE and a variable signal location. We argue that our theory method results, for MAP reconstruction, in a speedup of the order of 102–103 relative to methods employing samples of reconstructions. To see this, we first examine the number of reconstructions needed by sample methods. This number comprises training-set reconstructions and test-set reconstructions.

For training set reconstructions, one must calculate enough samples to learn the observer templates. For the NPW and CNPW observers, this involves one noiseless signal-present reconstruction per grid location and one signal absent reconstruction, since the template involves only s^¯loc. That is, the quantity s^¯loc=b+sloc^¯b^¯ and b+sloc^¯ needs one reconstruction per grid location, while b^¯ requires a single reconstruction. Thus, Nloc + 1 reconstructions are needed for the CNPW and NPW training phase. For the CHO, one must additionally estimate a channel covariance matrix inverse at each location. To obtain a stable estimate of the sample channel covariance matrix, one requires a number of reconstructions 10 – 100 times the number of channels [8], [26]. We presume that if one is using a sample method, then this requirement is followed. This typically requires a few hundred signal-absent reconstructions [8]. Thus, of the order of 102, denoted by An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102), reconstructions are needed for the CHO training phase. For the HO observer, the size of the inverse covariance matrix often [27] precludes learning from samples. Note that in our MC validations, we intentionally bypassed this training phase by analytically computing the observer templates. This assured that the MC and analytical FOM calculations used the same templates and that any MC-theory discrepancies were not affected by template differences.

For the test-set reconstructions, the number of reconstructions is roughly independent of the observers and is given by the number of reconstructions needed to get enough observer responses to obtain an accurate LROC curve. We found that to obtain error bars of the magnitude seen in Fig. 4 we required 500 signal present and 500 signal absent reconstructions. This number is of the order of a few tens per location for the signal present case plus a few hundred for the signal absent case, for a total of An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(10)Nloc + An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102).

Adding the training and test requirements, we get for the NPW and CNPW observers, a total of Nloc + 1 An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(10)Nloc + An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102) ≈ An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(10)Nloc + An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102) reconstructions. For our own case of Nloc = 25, this number is ≈ An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102) or higher. For the CHO, the total is An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102) + An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(10)Nloc + An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102) ≈ An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(10)Nloc + An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102), which in our case would be An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102). For the HO, the computation is intractable.

The theory calculation is dominated by the computation of the term κ which takes about one backprojection. Since the sample methods (as we have prescribed them) take An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102) reconstructions and each iterative reconstruction An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(10) iterations (one iteration takes approximately one projection plus one backprojection), then the theory method results in a reduction of An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg (103) in effort, relative to sample methods.

One might attempt a fast sample method in lieu of our theory approach where, for example, for the CHO, only An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(10) reconstructions were used for both the training and test phases. It remains to be seen whether the error in the FOM due to propagating this considerable sample error is of the same magnitude as that of the theory method. Even for the fast sample methods, however, our theory method still offers an An external file that holds a picture, illustration, etc.
Object name is nihms195647ig2.jpg(102) speedup.

The MC validations for the four observers showed reasonable agreement with theory. The agreement was not perfect, especially for the Hotelling observer. We also note that the agreement would not have been as good if we had chosen a higher value of β that yielded overly smooth reconstructions. We are looking to improve agreement of theory vs. MC by re-examining some of the assumptions in our theory. Two crucial assumptions were independence and normality of the local observer responses. We are currently working to extend our theory to include correlation between observer responses so as to handle problems due to local covariance and template widths that are wide with respect to the distance between signal locations. We also note that our theory apparently rested on the space-invariant nature of the imaging system (cf. equations (30) and (31)). However, we may extend our results to space-variant systems using results from [6], [8] and [9].

While our theory methods are fast, sample methods are more flexible. They are not dependent on independence assumptions and allow precise placement of signals on finely spaced location grids. The fast Fourier based expressions in our theory become less accurate for non-quadratic edge-preserving regularizers, a problem which does not affect sample-based methods.

No attempt was made to see if observer performance correlated with that of humans. This correlation ultimately depends on observer choice. Gifford et al. [17], [18] have conducted human vs. model observer LROC tests (using sample methods) with some of the observers listed here. In [18], the observer responses were converted to rating data and curve-fitting software [21] used to obtain ALROC. For the task in [18], it appears that channelized forms of the observer correlated well with human performance. Our work may allow rapid computational prototyping for the types of sample-based testing reported in [17], or for human LROC testing.

Acknowledgments

This work is being supported by the grant NIH-NIBIB EB02629. We thank Prof. Petar Djuric for useful discussions.

References

1. Barrett HH, Yao J, Rolland JP, Myers KJ. Model Observers for Assessment of Image Quality. Proc Natl Acad Sci USA. 1993 November;90:9758–9765. [PubMed]
2. Barrett HH, Wilson DW, Tsui BMW. Noise Properties of the EM Algorithm. I. Theory. Phys Med Biol. 1994 May;39(5):833–846. [PubMed]
3. Fessler J. Mean and Variance of Implicitly Defined Biased Estimators (such as Penalized Maximum Likelihood): Applications to Tomography. IEEE Trans Image Processing. 1996 March;5(3):493–506. [PubMed]
4. Wang W, Gindi G. Noise Analysis of MAP-EM Reconstruction Algorithms in Emission Tomography. Phys Med Biol. 1997;42(11):2215–2232. [PubMed]
5. Qi J, Leahy R. A Theoretical Study of the Contrast Recovery and Variance of MAP Reconstructions from PET Data. IEEE Trans Med Imaging. 1999 April;18(4):293–305. [PubMed]
6. Qi J, Leahy R. Resolution and Noise Properties of MAP Reconstruction for Fully 3-D PET. IEEE Trans Med Imaging. 2000 May;19(5):493–506. [PubMed]
7. Xing Y, Hsiao IT, Gindi G. Efficient Calculation of Resolution and Variance in 2D Circular-Orbit SPECT. Conf. Rec. IEEE Nuc. Sci. Symp. Med. Imaging Conf; Nov. 2001; IEEE; pp. 1298–1302.
8. Xing Y, Hsiao IT, Gindi G. Rapid Calculation of Detectability in Bayesian SPECT. Phys Med Biol. 2003 Nov;48(22):3755–3773. [PubMed]
9. Stayman W, Fessler J. Fast Methods for Approximation of Resolution and Covariance in SPECT. Conf Rec. IEEE Nuc. Sci. Symp. Med. Imaging Conf; Nov. 2002; IEEE; pp. 786–788.
10. Qi J, Huesman R. Theoretical Study of Lesion Detectability of MAP Reconstruction Using Computer Observers. IEEE Trans Med Imaging. 2001 August;20(8):815–822. [PubMed]
11. Bonetto P, Qi J, Leahy R. Covariance Approximation for Fast and Accurate Computation of Channelized Hotelling Observer Statistics. IEEE Trans Nuclear Science. 2000 August;47(4):1567–1572.
12. Fessler J, Yendiki A. Channelized Hotelling Observer Performance for Penalized-Likelihood Image Reconstruction. Conf. Rec. IEEE Nuc. Sci. Symp. Med. Imaging Conf; Nov. 2002; IEEE; pp. 1040–1044.
13. Qi J. Theoretical Evaluation of the Detectability of Random Lesions in Bayesian Emission Reconstruction. Conf. Rec. Info. Processing In Med. Imag; July 2003; pp. 354–365. [PubMed]
14. Qi J. Analysis of Lesion Detectability in Bayesian Emission Reconstruction With Nonstationary Object Variability. IEEE Trans Med Imaging. 2004 March;23(3):321–329. [PubMed]
15. Gifford H, Wells R, King M. A Comparison of Human Observer LROC and Numerical Observer ROC for Tumor Detection in SPECT. IEEE Trans Nuclear Science. 1999 Aug;46(4):1032–1037.
16. Xing Y, Oldan J, Kulkarni S, Khurd P, Gindi G. Effects of Image Quantization and Regularization on Lesion Detectability in SPECT MAP Reconstruction. Conf. Rec. IEEE Nuc. Sci. Symp. Med. Imaging Conf; Oct. 2003; IEEE;
17. Gifford HC, Farncombe TH, King MA. Ga-67 Tumor Detection Using Penalized-EM with Nonanatomical Regularizers. Conf. Rec. IEEE Nuc. Sci. Symp. Med. Imaging Conf; Nov. 2002; pp. 1397–1401.
18. Gifford HC, Pretorius PH, King MA. Comparison of Human-and Model-observer LROC Studies. Conf. Rec. SPIE Med. Imaging; Feb. 2003; pp. 112–122.
19. Abbey C. PhD thesis, Graduate Interdisciplinary Program in Applied Mathematics. University of Arizona; Tuscon, Arizona, USA: Apr, 1998. Assessment of Reconstructed Images.
20. Swensson R, Judy P. Detection of Noisy Visual Targets: Models for the Effects of Spatial Uncertainty and Signal-to-noise Ratio. Perception and Psychophysics. 1981;29(3):521–534. [PubMed]
21. Swensson R. Unified Measurement of Observer Performance in Detecting and Localizing Target Objects on Images. Medical Physics. 1996 Oct;23(10):1709–1725. [PubMed]
22. Myers KJ, Barrett HH. Addition of a Channel Mechanism to the Ideal-observer model. J Optical Soc America A. 1987 Dec;4(12):2447–2457. [PubMed]
23. Fessler JA, Rogers WL. Spatial Resolution Properties of Penalized-likelihood Image Reconstruction: Space-invariant Tomographs. IEEE Trans Image Processing. 1996 Sept;5(9):1346–1358. [PubMed]
24. 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.
25. Zoubir AM, Boashash B. The Bootstrap and its Application in Signal Processing. IEEE Signal Processing Magazine. 1998 Jan;15(1):56–76.
26. Sain J, Barrett HH. Performance Evaluation of a Modular Gamma Camera Using a Detectability Index. J Nuc Med. 2003;44:58–66. [PubMed]
27. Barrett H, Myers K, Gallas B, Clarkson E, Zhang H. Megalopinakophobia: Its Symptoms and Cures. Conf. Rec. SPIE Med. Imaging; SPIE; 2001. pp. 299–307.