1.

Introduction

Model observers have been widely used to perform task-based assessment of medical image quality. Two important categories of model observers are ideal and anthropomorphic observers. The ideal observer (IO) outperforms all other observers and sets an upper limit on task performance measured by figures of merit such as the area under the receiver operating characteristic (ROC) curve (AUC).¹ The IO makes optimal use of all the information in the raw data. However, it requires full knowledge of the raw data statistics. Therefore, it is typically applied in the projection domain, where the quantum noise statistics are known and easier to analyze. Moreover, IO performance is not improved by invertible operations on the raw data such as linear filtering or (invertible) reconstruction algorithms, and it thus allows optimization of instrumentation or reconstruction parameters in the projection domain and provides an alternative to image-domain optimization.

Anthropomorphic model observers, such as the channelized Hotelling observer (CHO) using anthropomorphic channels,² are designed to predict how humans would perform. They have been extensively used in the evaluation and optimization of acquisition, instrumentation, reconstruction, and compensation methods in the context of myocardial perfusion, bone, and hepatic single-photon emission computed tomography (SPECT).^3–6 The CHO has shown good agreement with human observers in a variety of clinical applications and tasks.^3–5,7–13 For SPECT, the CHO is typically applied to reconstructed images. Since task performance for both humans and the CHO depends on details of compensation and regularization methods and parameters, comparing performances in this domain should be done with images that are optimized in terms of relevant reconstruction and regularization parameters such as iteration number and the cutoff frequency of any postreconstruction smoothing filter. This optimization requires significant computational resources, and it is thus simpler and more efficient to optimize imaging systems in the projection domain using the IO rather than in the reconstructed image domain using the CHO.

Despite the potential advantages of using the IO, it has rarely been rigorously applied to real imaging tasks due to the technical difficulty of estimating the performance for realistic background variability arising from variations in patient anatomy and uptake. Its application has been limited to cases where the background and signal models could be expressed analytically, which often do not capture the variability observed in clinical studies.^14–16 To overcome this limitation, methods based on Markov Chain Monte Carlo (MCMC) techniques have been developed and applied to estimate the IO test statistic, i.e., the likelihood ratio (LR) for realistic and general backgrounds and signal models for binary defect detection tasks.^17–19

While the development of MCMC techniques has allowed the IO to be applied to more realistic objects and ensembles of objects, it is unclear if the acquisition and instrumentation parameters optimized using the IO in the projection domain will also be optimal when the reconstructed images are used for visual interpretation by human observers. It is clear that the optimization using the IO results in projection data that contain the maximum possible information about the detection task. However, for some sets of reconstruction algorithms and compensation methods, there will undoubtedly be information loss and thus the instrumentation and reconstruction parameters optimal for the IO may not be optimal for the reconstruction and compensation methods of interest.

In this work, we address the less general reconstruction optimization problem of finding optimal acquisition parameters for a particular set of reconstruction and compensation methods. In this case, we restrict ourselves to iterative reconstruction using various scatter compensation methods (including no scatter compensation) and postreconstruction low-pass filtering regularization for noise control.

Based on the above discussion, two ways to approach this problem are as follows. In sequential optimization, the acquisition and instrumentation parameters are first optimized using the IO. These parameters are then fixed and reconstruction parameters, such as number of iterations and low-frequency cutoff parameters, are optimized using an anthropomorphic observer. An alternative is joint optimization where the acquisition, instrumentation, and reconstruction parameters are optimized at the same time. As we have argued above, and present data to support below, joint optimization is very time consuming, and the sequential optimization process is appealing.

However, one potential difficulty with the sequential optimization procedure is that there can be some information loss in the reconstruction process. In this work, we investigate one particular source of such information loss: differences in the physical model used or assumed by the reconstruction algorithm compared with the true image formation process. In this work, we propose introducing information loss into the IO formalism in a way that will be described more fully below. We then investigate whether sequential optimization using this observer results in acquisition parameters, i.e., the energy windows, that are more similar to those from joint reconstruction. We investigate this in the context of reconstruction methods that include models of the image formation process with varying degrees of fidelity to the true process.

To understand how we propose to implement this information loss in the IO, recall that, since the IO has full knowledge of all the statistical information about the background and signal, it implicitly has a perfect model of the image formation process including all physical processes and factors that degrade image quality (e.g., finite energy resolution, scatter, and attenuation). In other words, the image formation models embedded in the IO match perfectly the models used to form the projection images, $p_t()$ (Fig. 1). In some sense, the IO can be viewed as providing ideal compensation for the various image degrading effects, and the task performance it achieves reflects the best achievable with perfect compensation. However, real reconstruction and compensation algorithms do not perfectly model the image formation process and, in fact, there are often physical effects for which no compensation is provided.

Fig. 1

A schematic diagram of the image formation process.

An important class of reconstruction methods to consider for the optimization task is statistical iterative reconstruction methods such as ordered subsets-expectation maximization (OS-EM). One advantage of these methods is that they allow incorporating complete (and incomplete) models of the underlying physics of the imaging process. This is accomplished by modeling the physics in the projection operator, $p_m()$, used in the algorithms to compute projections from the current activity estimate. The transpose of the projection operator is used in the algorithm to update the image estimate after comparing the computed and measured projections. In SPECT, the projection formation process is continuous-to-discrete, spatially varying, and patient dependent. Closed form expressions are not available. For realistic object shapes and Monte Carlo simulations that accurately model the nuclear medicine image formation process, it is currently not possible to model the projection without digitizing the object. As a result, the true projection operator $p_t()$ was approximated as discrete-to-discrete, but was applied to a phantom with a voxel size half that of the projection bin size. Projection operators used in the reconstruction algorithms, $p_m()$, are also approximations of the true operator. They are discrete-to-discrete and model the imaging physics to varying degrees of accuracy. This operator was also applied to the phantom having the same voxel size used for the true projection operator. The difference between the true projection model, $p_t()$, and that used in reconstruction, $p_m()$ is referred to as model mismatch (MM). We hypothesize that MM could lead to differences in optimal acquisition and instrumentation parameters when using sequential and joint optimizations.²⁰

We have previously introduced the concept of including MM into the IO framework.^21,22 We called this observer the IO with MM (IO-MM). The IO-MM still uses the LR as the test statistic, but in the estimation of the LR it uses an approximate model of the image formation process [$p_m()$] instead of the true model [$p_t()$]. In Ref. 22, we compared the performance of the IO and IO-MM with that of the CHO for the task of optimizing the energy window for myocardial perfusion SPECT (MPS) imaging. The results for the IO-MM showed good agreement with that of the CHO. We have also used the IO and IO-MM to optimize the energy window in ${}^{90}\mathrm{Y}$ bremsstrahlung SPECT for a detection task.²³ The optimal energy window was narrower when taking into account MM and was similar to that obtained previously when optimized for an estimation task.²⁴

In MPS imaging using Tc-99m, due to the interactions of the Tc-99m photons in the body and collimator-detector system and the limited energy resolution of scintillation cameras, the recorded energy spectrum is degraded. Detection of scattered photons degrades the final image quality and thus has a significant impact in limiting the diagnostic accuracy of readers interpreting the images. Therefore, it is desirable to reduce the effects of scatter on the final image quality by either reducing the number of detected scattered photons or compensating for them.

One way to reduce the number of detected scattered photons is to use narrower acquisition energy windows. However, there is a tradeoff between the width of the acquisition energy window and image noise: reducing the width of the energy window decreases the number of detected primary and scattered photons concurrently. This would result in an increase of the noise in the projection data, which could degrade the performance of human observers operating on the final images. In this work, we seek the optimal acquisition energy window width that provides the best tradeoff between reducing scatter and increasing noise, as measured by performance on a perfusion defect detection task.

Scatter compensation can result in improved quantification, lesion detection, and contrast-to-noise ratios.^4,25–28 Current methods for scatter compensation can be divided into two general categories based on how the scatter component is estimated: energy and spatial estimation methods. In the former category, scatter in the photopeak image is estimated using information acquired in one or more additional energy windows.²⁹ The dual and triple energy window methods are examples of this approach and have been widely adopted due to their simplicity and effectiveness.³⁰ In the dual energy window (DEW) method, an appropriately scaled projection image acquired in a relatively narrow energy window immediately below the photopeak window is used as an estimate of the scattered component of the photopeak window. For the triple energy window (TEW) method, images used for scatter estimation are acquired in energy windows below and above the photopeak window. A disadvantage of window-based methods is the noise in the scatter estimate that results from acquiring data in these narrow windows. Increasing the width of the scatter window reduces the noise in the scatter estimate, but provides, in general, a more biased estimate of scatter in the photopeak window.³¹

The alternative to energy-based scatter estimates is spatial modeling. Scatter modeling methods estimate the scatter component of the projection data based on an estimate derived from the reconstructed image.^32,33 They can naturally be applied as part of iterative reconstruction algorithms. An example of model-based scatter estimation methods is the effective source scatter estimation (ESSE) method.³³ It uses a set of scatter kernels that are precalculated using Monte Carlo simulation methods to estimate an effective scatter source whose attenuated projection gives the scatter contribution in the projection data. It has previously been applied in a number of SPECT imaging applications for scatter, downscatter, and crosstalk estimation and provided accurate compensation.^34,35

In this work, we extended the work introduced in Ref. 36 to compare different scatter estimation methods, including the DEW, TEW, and ESSE methods in the context of MPS, and to find the optimal acquisition energy window width that provides the best performance on a binary defect detection task using the IO and IO-MM. We also compared optimal energy windows obtained with the IO and IO-MM to that of an anthropomorphic observer (CHO) applied to images reconstructed using the OS-EM algorithm. In other words, we investigated the agreement in optimal energy windows between sequential optimization using the IO-MM or IO and joint optimization using the CHO. In addition to providing information about optimal energy windows and scatter copensation strategies for MPS, this study investigates the use of the IO and IO-MM in the context of projection-domain optimization of acquisition parameters in comparison with reconstructed image-domain optimization using an anthropmorphic CHO in the presence of varying degrees of MM.

2.

Methods

In this section, we introduce the key elements of task-based image quality assessment including the task, object, imaging system models, observer, and figure of merits.¹ We revisit the concept of MM and give a brief overview of the IO, the use of MCMC to estimate IO performance, and the modifications needed in the IO to take MM into account. Finally, we give a brief description of the CHO that was used as a surrogate for human observer performance.

2.2.

Object and Imaging System Models

It has been previously reported that the variability in patients’ anatomies and activity uptakes are important factors that limit task performance. Thus, it is important that the object model be realistic and represent the range of variability seen in clinical populations.^5,37–40 To this end, we have previously designed and developed a digital phantom population based on the three-dimensional (3-D) eXtended CArdiac Torso (XCAT) phantom and generated the corresponding Tc-99m projection data.⁴¹ The population consists of 54 phantoms including anatomical models for both genders and three variations (small, medium, and large) in body habitus, subcutaneous adipose tissue thickness, and heart size. The anatomical parameters of the phantom population were obtained from the anatomical distributions in the Emory Cardiac Database (Barclay, Emory University),⁴² which includes the anatomical measurements from 166 patients. The phantoms were digitized into 0.221 cm cubic voxels. We simulated anterolateral and inferior perfusion defects with a 10% extent and 25% severity as shown in Fig. 2. The extent and severity are the fraction of the left ventricular volume and reduction of uptake in the perfusion defect relative to the normal myocardium, respectively. The extent and severity of the defect were chosen to be challenging and clinically relevant.

Fig. 2

Sample short-axis slice of the heart showing a perfusion defect at the anterolateral and inferior locations of the myocardium wall. For illustrative purposes, defects shown have 100% severity.

Projections of the phantom population were simulated using the SimSET Monte Carlo code⁴³ and angular response functions (ARFs)⁴⁴ to model interactions in the body and the collimator detector system, respectively. We modeled a GE Infinia dual-detector SPECT system with a GE low-energy high-resolution collimator and a 9.5-mm thick NaI(Tl) crystal. Projections were simulated using noncircular phantom-specific body-contouring orbits and acquired at 60 equispaced angles over a 180 deg acquisition arc extending from 45 deg right anterior oblique to 45 deg left posterior oblique. Projection images were binned in a $128\times 114$ matrix with a pixel size of 0.442 cm. We used a smaller voxel size in the phantom and Monte Carlo simulation to realistically simulate the continuous spatial nature of the anatomy in patients.

To model uptake variability, we simulated low-noise Tc-99m projection data for the heart, liver, and body (including all other organs) for each phantom. This enabled us to separately scale the individual projection images of the different organs using random scale factors obtained from distributions that were based on data from a set of 34 patients who underwent MPS. We could then sum the scaled projections before simulating Poisson noise to generate an uptake realization. In this study, we modeled an injected activity of 10 mCi of Tc-99m to model a low-dose protocol.

Optimizing energy windows requires projection images acquired in various acquisition energy windows. The flexible and efficient design of the phantom population and the projection database in Ref. 41 permitted us to compute projections in any arbitrary energy window. The original projection data were simulated in 1 keV wide energy bins with no measurement-related energy blurring. We computed projections in different energy windows by summing the data from these narrow bins. The summing process took into account the energy-dependent energy resolution function. We investigated 11 acquisition energy windows, labeled ${\mathrm{W}}_1$ to ${\mathrm{W}}_{11}$, with the parameters shown in Table 1. Figure 3 shows the sample low-noise and noisy projection images acquired in the different energy windows. Note the tradeoff between increasing scatter, which reduces the contrast of the myocardium, and reduced noise as the energy window width increases. The goal was to seek the window among those listed below that provides the best task performance.

Table 1

Investigated Tc-99m energy window settings.

	W1	W2	W3	W4	W5	W6	W7	W8	W9	W10	W11
Lower threshold (keV)	140	136	132	130	128	126	124	122	120	112	102
Upper threshold (keV)	141	145	149	151	153	155	157	159	161	169	179
Width (keV)	1	9	17	21	25	29	33	37	41	57	77
Width (%)	0.7	6.4	12.1	15	17.8	20.6	23.5	26.3	29.2	40.6	54.8

Fig. 3

Sample low noise (top) and noisy (bottom) Monte Carlo simulated anterior projection images of a medium-sized male phantom with a medium-sized heart and fat level acquired in energy windows ${\mathrm{W}}_1$ to ${\mathrm{W}}_{11}$ (from left to right). The images are displayed using a logarithmic grayscale to better show the lower uptake organs.

In order to evaluate different scatter estimates, and thus different compensation methods, we generated scatter estimates for the DEW, TEW, and ESSE methods. For the DEW method, scatter estimates were generated in a 10-keV wide window below and adjacent to the photopeak window. Scatter estimates for the TEW method were acquired in two 5 keV wide windows immediately below and above the photopeak energy window. For both the DEW and TEW methods, we calculated the scatter estimates using methods reported in Ref. 30. For ESSE, the projections were calculated directly for the photopeak window using the appropriate scatter kernels. Figure 4 shows the true and estimated scatter projections obtained with the various methods.

Fig. 4

Low-noise scatter projection images generated using SimSET and ARFs (true), dual energy window (DEW), triple energy window (TEW), and effective scatter source estimation (ESSE) methods (from top to bottom) for a medium sized male phantom with medium size heart and fat level acquired from anterior projection view in energy windows ${\mathrm{W}}_1$ to ${\mathrm{W}}_{11}$ (from left to right). The images are displayed using a logarithmic grayscale to better show the lower uptake organs.

2.3.

Observers Models

2.3.1.

Application of the ideal observer and the ideal observer with model mismatch

In SPECT imaging, given a continuous object, $f$, the projection data, $\mathbf{g}$, are the result of the image formation process and can be represented by

Eq. (1)

$$\mathbf{g}=p_t(\mathbf{f})+\mathbf{n},$$

where $p_t()$ is a projection operator that maps the object, $\mathbf{f}$, to the projection space, and $\mathbf{n}$ is the measurement noise.

In this work, we consider the task of detecting a fixed signal, ${\mathbf{f}}_s$, such as a perfusion defect, in a random background, ${\mathbf{f}}_b$. We define the background and signal projection images as

Eq. (2)

$${\mathbf{b}}_t=p_t({\mathbf{f}}_b),$$

Eq. (3)

$${\mathbf{s}}_t=p_t({\mathbf{f}}_s)\mathrm{.}$$

For a binary detection task, the two hypotheses to be tested can be written as

Eq. (4)

$$H_0\text{:}\mathbf{g}={\mathbf{b}}_t+\mathbf{n},$$

Eq. (5)

$$H_1\text{:}\mathbf{g}={\mathbf{b}}_t+{\mathbf{s}}_t+\mathbf{n},$$

where $H_0$ and $H_1$ mean signal absent and signal present, respectively.

The IO uses the LR, $\mathrm{\Lambda }(\mathbf{g})$, of defect-present versus defect-absent as the test statistic. For the case where the background and signal are known exactly (SKE/BKE) and the only source of randomness is Poisson measurement noise, $\mathbf{n}$, the expression for the LR is given by

Eq. (6)

$${\mathrm{\Lambda }}_{\mathrm{BKE}}(\mathbf{g},{\mathbf{b}}_{\mathrm{t}})=\frac{\mathrm{pr}(\mathbf{g}|H_1)}{\mathrm{pr}(\mathbf{g}|H_0)}=\prod _i{(1+\frac{{\mathrm{s}}_{\mathrm{t}}(i)}{{\mathrm{b}}_{\mathrm{t}}(i)})}^{\mathrm{g}(\mathrm{i})}\exp [-{\mathrm{s}}_{\mathrm{t}}(\mathrm{i})].$$

For a signal known exactly embedded in a background known statistically (SKE/BKS), the expression for the LR is given by

Eq. (7)

$${\mathrm{\Lambda }}_{\mathrm{BKS}}(\mathbf{g})=\int {\mathrm{\Lambda }}_{\mathrm{BKE}}(\mathbf{g},{\mathbf{b}}_t)\mathrm{pr}({\mathbf{b}}_t|\mathbf{g},H_0)\mathrm{d}b.$$

The idea of Eq. (7) is to average over the backgrounds sampled randomly from the posterior density $\mathrm{pr}({\mathbf{b}}_t|\mathbf{g},H_0)$. However, the density $\mathrm{pr}({\mathbf{b}}_t|\mathbf{g},H_0)$ is usually not known analytically and is hard to sample using standard Monte Carlo integration methods. Accordingly, methods based on MCMC techniques have been developed to estimate the integral in Eq. (7).^17–19

We applied a previously developed method to efficiently estimate the LR.¹⁹ We parameterized the object with a parameter vector, $\stackrel{\rightharpoonup }{\theta }$, which was a combination of discrete anatomical parameters, $n$, $p,q$, and $l$, that specify the object gender, body size, heart size, and fat level, respectively, and continuous organ activity parameters, $A_{\text{heart}}$, $A_{\text{liver}}$, and $A_{\text{body}}$, corresponding to the uptake activities of the heart, liver, and body, respectively. Thus, the background ${\mathbf{b}}_t(\stackrel{\rightharpoonup }{\theta })$ parameterized by the vector $\stackrel{\rightharpoonup }{\theta }$ was represented as

Eq. (8)

$${\mathbf{b}}_t(\stackrel{\rightharpoonup }{\theta })=A_{\text{heart}}{\mathrm{Prj}}_{n,p,q,l}^{\text{heart}}+A_{\text{liver}}{\mathrm{Prj}}_{n,p,q,l}^{\text{liver}}+A_{\text{body}}{\mathrm{Prj}}_{n,p,q,l}^{\text{body}},$$

where ${\mathrm{Prj}}_{n,p,q,l}^{\text{heart}}$, ${\mathrm{Prj}}_{n,p,q,l}^{\text{liver}}$ and ${\mathrm{Prj}}_{n,p,q,l}^{\text{body}}$ were the projections of the heart, liver, and the body, respectively, filled with unit activity for a given anatomy. The advantage of this approach is that the projections can be precalculated so that the background can be estimated very rapidly simply by scaling and summing the set of three projection images.

It has been shown in Ref. 16 that, if ${\mathbf{b}}_t(\stackrel{\rightharpoonup }{\theta })$ has a one-to-one relation with the parameter vector $\stackrel{\rightharpoonup }{\theta }$, the integral in Eq. (7) can be expressed as

Eq. (9)

$${\mathrm{\Lambda }}_{\mathrm{BKS}}(\mathbf{g})=\int {\mathrm{\Lambda }}_{\mathrm{BKE}}(\mathbf{g},{\mathbf{b}}_t(\stackrel{\rightharpoonup }{\theta }))\mathrm{pr}(\stackrel{\rightharpoonup }{\theta }|\mathbf{g},H_0)\mathrm{d}\stackrel{\rightharpoonup }{\theta }.$$

Adopting a Metropolis–Hastings approach to draw samples from the posterior distribution $\mathrm{pr}(\stackrel{\rightharpoonup }{\theta }|\mathbf{g},H_0)$, we could then compute the ergodic average ${\tilde{\mathrm{\Lambda }}}_{\mathrm{BKS}}(\mathbf{g})$, which gives an estimate of the integral in Eq. (9).

It is relatively straightforward to use this formulation to take into account the “mismatch” between the models used in reconstruction and the real physical image formation process. For the IO-MM and a binary decision task, the two hypotheses to be tested were represented mathematically by

Eq. (10)

$$H_0^{*}\text{:}\mathbf{g}={\mathbf{b}}_m+\mathbf{n},$$

Eq. (11)

$$H_1^{*}\text{:}\mathbf{g}={\mathbf{b}}_m+{\mathbf{s}}_m+\mathbf{n},$$

where $\mathbf{g}$ was the input projection image generated using Eq. (1), and the background and signal projection images, ${\mathbf{b}}_m$ and ${\mathbf{s}}_m$, respectively, were defined as

Eq. (12)

$${\mathbf{b}}_m=p_m({\mathbf{f}}_b),$$

Eq. (13)

$${\mathbf{s}}_m=p_m({\mathbf{f}}_s).$$

For the IO-MM, the expression for the LR for the SKE/BKE case is given by

Eq. (14)

$${\mathrm{\Lambda }}_{\mathrm{BKE},\mathrm{m}}(\mathbf{g},{\mathbf{b}}_{\mathrm{m}})=\frac{\mathrm{pr}(\mathbf{g}|H_1^{*})}{\mathrm{pr}(\mathbf{g}|H_0^{*})}=\prod _i{(1+\frac{{\mathrm{s}}_{\mathrm{m}}(\mathrm{i})}{{\mathrm{b}}_{\mathrm{m}}(\mathrm{i})})}^{\mathrm{g}(\mathrm{i})}\exp [-{\mathrm{s}}_{\mathrm{m}}(i)].$$

When background variability is included, the equation for the LR becomes

Eq. (15)

$${\mathrm{\Lambda }}_{\mathrm{BKS},\mathrm{m}}(\mathbf{g})=\int {\mathrm{\Lambda }}_{\mathrm{BKE},\mathrm{m}}(\mathbf{g},{\mathbf{b}}_{\mathrm{m}})\mathrm{pr}({\mathbf{b}}_{\mathrm{m}}|\mathbf{g},H_0)\mathrm{d}b.$$

To compute the LR for the IO-MM, we followed the same methodology as for the case of the IO. We parameterized the object with the parameter vector, $\stackrel{\rightharpoonup }{\theta }$, applied the MCMC method to sample from the posterior distribution, $\mathrm{pr}(\stackrel{\rightharpoonup }{\theta }|\mathbf{g},H_0^{*})$, and computed the ergodic average ${\tilde{\mathrm{\Lambda }}}_{\mathrm{BKS},\mathrm{m}}(\mathbf{g})$. This was implemented using the same computer code as for the IO by providing a database of organ projections that included the desired MM instead of the true projections.

In the standard IO, we tested the hypothesis that a projection image, $\mathbf{g}$, contains or does not contain a fixed defect in a random background, where both the signal and the background images were generated using the true model of the image formation process. Thus, the IO had a perfectly matched model of all physical image-degrading processes. For the IO-MM, we still tested the same hypothesis. However, the signal and the background images provided to the observer were generated using an approximate model, $p_m()$, of the image formation process. In both cases, the input projection data were generated using the true model.

For each energy window, we computed the IO and IO-MM (in the cases where the observer had an approximate model of scatter based on the ESSE, DEW, and TEW methods) test statistics using the MCMC method for an ensemble of 2160 pairs of defect-present and defect-absent projection images modeling different anatomical and uptake variations. The projection data contained a $64\times 24\times 60$ pixel region of interest centered over the centroid of the heart. This ensured that the whole heart was included inside the region-of-interest. For each energy window, the IO test statistics were used as the input to the ROCkit code,⁴⁵which fits a binormal ROC curve to the input set of LRs to estimate the AUC.

2.3.2.

Channelized Hotelling Observer Study

We used an anthropomorphic CHO as a surrogate for human observer performance and compared optimal energy window settings and the different scatter compensation methods to those obtained from the IO and the IO-MM.

Image reconstruction

We performed CHO studies on images reconstructed using the OS-EM algorithm with four combinations of compensations: attenuation only (A); attenuation, spatially varying geometric response, and scatter using ESSE scatter modeling (AG-ESSE); attenuation, spatially varying geometric response, and scatter using the DEW method (AG-DEW); and attenuation, spatially varying geometric response, and scatter using the TEW method (AG-TEW). Figure 5 shows sample transaxial slices from images reconstructed using the different compensation methods. For each energy window and compensation method, we used 10 pairs of projection datasets per defect position per phantom, leading to a total of 1080 pairs of defect-present and -absent images. We used 12 subsets with five projections per subset in the OS-EM algorithm. We considered iterations 1, 2, 3, 5, 7, 10, 12, 15, 18, and 20 in choosing the optimal one.

Fig. 5

Sample transaxial images located at the heart centroid for different phantoms and the corresponding attenuation maps (rows 1 and 2). Rows 3 to 6 show the reconstructed slices obtained using A, AG-DEW, AG-TEW, and AG-ESSE compensation methods after the third iteration with 12 subsets per iteration.

Post-reconstruction processing

After reconstruction, images from each of the iterations identified above were filtered using a low-pass Butterworth filter with order 8 and cutoff frequencies of 0.08, 0.1, 0.12, 0.14, 0.16, 0.2, and $0.24{\text{pixels}}^{-1}$. We reoriented the filtered transaxial images to the short-axis view (standard view of MPS images; orthogonal to the long axis of the left ventricle). For each filtered and reoriented image, we extracted the short-axis slice containing the center of the defect (or the same slice in the defect-absent short-axis image). Following the methodologies adopted in Refs. 4 and 46, the pixel values in the filtered, reoriented images were windowed to 256 gray levels by scaling the image so the maximum value in the heart was 255, and negative values were set to zero.

Application of the channelized Hotelling observer

In this study, we used four nonoverlapping difference-of-mesa frequency channels with successively doubling widths and a starting frequency of $0.5{\text{pixels}}^{-1}$. We estimated CHO performance using a leave-one-out technique in which the CHO was trained using all but one image and tested using the left-out image, producing one test statistic value. This process was repeated with each of the feature vectors being left out, resulting in a set of test statistics with as many members as images. These test statistics were used as inputs to the ROCkit code to estimate the AUC. For each energy window and compensation method, the iteration and cutoff frequencies giving the highest AUC were selected as optimal. A schematic diagram showing the steps of the CHO study is shown in Fig. 6.

Fig. 6

Schematic diagram showing the steps used in the channelized Hotelling observer (CHO) study.

3.

Results

3.1.

IO and IO-MM Studies

Figure 7 shows a plot of the AUC values of the IO, i.e., the case where the observer used a perfect scatter model for the detection task as a function of acquisition window width. The standard deviations, estimated by ROCkit, were very small ($\sim 0.003$) and thus are not plotted. The performance of the IO operating on the images acquired in energy window, ${\mathrm{W}}_1$, a 1-keV wide energy window, was poor. This was a case where quantum noise was dominant due to fewer detected photons, and significantly limited the observer performance. We observed that the IO performance in the range of energy window widths 9 to 41 keV (${\mathrm{W}}_2$ to ${\mathrm{W}}_9$) did not change substantially, indicating that the observer could handle moderate amounts of scatter. Despite the fact that the IO had a perfect scatter model, the performance was slightly worse for energy windows of widths more than 41 keV (${\mathrm{W}}_{10}$ and ${\mathrm{W}}_{11}$). This is likely because the amount of scatter was very high and resulted from very large angle and multiple scatter, thus adding little information relevant to the task. These results suggest that there is no benefit in increasing the acquisition energy window width to more than $\sim 41\mathrm{keV}$ (${\mathrm{W}}_9$), even when incorporating a perfect scatter model.

Fig. 7

The ideal observer (IO) performance measured in terms of the area under the receiver operating characteristic (ROC) curve (AUC) value for the different acquisition energy windows.

AUC values of the IO and IO-MM for the cases when scatter was modeled using the ESSE, DEW, and TEW methods are shown in Table 2. The standard deviations of the AUC values for all methods were, again, very small ($\sim 0.004$). We note that the maximum AUC value for the IO and IO-MM for the different cases was for an energy window width of 29 keV (${\mathrm{W}}_6$). For each scatter modeling method (perfect, ESSE, DEW, and TEW), we tested the statistical significance of the difference between the AUC value for a width of 29 keV (${\mathrm{W}}_6$) and the rest of the energy windows under investigation. The corresponding $p$-values are shown in Table 3. Italicized AUC and corresponding p-values indicate the cases where the differences were not statistically significant at the $p=0.05$ level. The results in Tables 2 and 3 indicated that, for the IO, acquisition energy windows of widths between 17 keV (${\mathrm{W}}_3$) and 41 keV (${\mathrm{W}}_9$) were near-optimal in the sense that we could not accept the hypothesis that the AUC values were significantly less than that for the 29-keV wide window. When the IO-MM incorporated a scatter model based on the DEW, TEW, and ESSE methods, the near-optimal range of energy window widths was between 21 and 33 keV (${\mathrm{W}}_4–{\mathrm{W}}_7$), 25 and 33 keV (${\mathrm{W}}_5–{\mathrm{W}}_7$), and 25 and 37 keV (${\mathrm{W}}_5–{\mathrm{W}}_8$), respectively.

Table 2

AUC values for the ideal observer (IO) and IO with model mismatch (IO-MM).

Energy window width	IO (Perfect)	IO-MM (DEW)	IO-MM (TEW)	IO-MM (ESSE)
1	0.794	0.792	0.723	0.789
9	0.946	0.943	0.855	0.941
17	0.949	0.949	0.938	0.946
21	0.949	0.949	0.954	0.948
25	0.951	0.948	0.957	0.949
29	0.956	0.955	0.960	0.958
33	0.954	0.951	0.959	0.958
37	0.954	0.947	0.954	0.956
41	0.947	0.948	0.950	0.953
57	0.949	0.928	0.931	0.940
77	0.937	0.926	0.924	0.935

Note: DEW, dual energy window; TEW, triple energy window; ESSE, effective scatter source estimation.

Table 3

p-values for hypothesis that the area under the receiver operating characteristic (ROC) curve (AUC) for the indicated energy window was different than that for the 29 keV width energy window (W6).

Energy window width	IO (Perfect)		IO-MM (DEW)		IO-MM (TEW)		IO-MM (ESSE)
	ΔAUC	p-value	ΔAUC	p-value	ΔAUC	p-value	ΔAUC	p-value
1	0.162	0	0.163	0	0.237	0	0.169	0
9	0.010	0.004	0.012	0	0.105	0	0.017	0.023
17	0.006	0.062	0.006	0.012	0.022	0	0.012	0.048
21	0.006	0.633	0.006	0.472	0.006	0.014	0.010	0.013
25	0.004	0.065	0.007	0.095	0.003	0.661	0.009	0.07
29	0	1	0	1	0	1	0	1
33	0.002	0.243	0.004	0.063	0.001	0.371	0	0.572
37	0.002	0.278	0.008	0.002	0.006	0.012	0.002	0.578
41	0.009	0.479	0.007	0.006	0.010	0	0.005	0.026
57	0.013	0.022	0.027	0	0.029	0	0.018	0
77	0.019	0.001	0.029	0	0.036	0	0.023	0

Note: p-values less than 10−5 are shown as 0.

We compared the performance of the IO and IO-MM when incorporating each of the different scatter models for each energy window. Table 4 shows the $p$-values for the hypothesis that the AUCs were the same for the various energy windows. We see that the IO-MM that used the ESSE scatter model was statistically indistinguishable (had a $p$-value $>0.05$) from the IO, which implicitly had a perfect scatter model, for all the energy windows. This indicates that the information about the scatter provided by ESSE was very similar to that in the true scatter in terms of its effect on defect detection. The performances of the IO-MMs that used the DEW and TEW scatter models were statistically indistinguishable from that of the IO for energy windows ${\mathrm{W}}_6$ to ${\mathrm{W}}_9$ and ${\mathrm{W}}_4$ to ${\mathrm{W}}_9$, respectively. However, the performance was significantly worse than the IO for energy window of widths more than 41 keV (${\mathrm{W}}_{10}$ and ${\mathrm{W}}_{11}$). These data suggest that using an accurate scatter estimate, such as that from ESSE, would allow for the use of a wider energy window without degrading task performance. The DEW or TEW estimates resulted in a performance similar to that achieved by the IO when the appropriate acquisition energy windows were used.

Table 4

Statistical significance of the differences between the AUC values of the IO and IO-MM for each scatter model.

Energy window width	Perfect versus DEW		Perfect versus TEW		Perfect versus ESSE
	ΔAUC	p-value	ΔAUC	p-value	ΔAUC	p-value
1	0.002	0.002	0.071	0	0.005	0.457
9	0.003	0.007	0.090	0	0.004	0.185
17	0	0.001	0.012	0.020	0.003	0.327
21	0.001	0.001	$\mathit{-}\mathit{0.007}$	0.113	0.001	0.531
25	0.003	0.004	$\mathit{-}\mathit{0.006}$	0.068	0.002	0.258
29	0	0.785	$\mathbf{-}\mathbf{0.005}$	0.063	$\mathbf{-}\mathbf{0.002}$	0.268
33	0.002	0.618	$\mathit{-}\mathit{0.005}$	0.066	$\mathit{-}\mathit{0.004}$	0.341
37	0.005	0.313	$\mathit{-}\mathit{0.001}$	0.073	$\mathit{-}\mathit{0.002}$	0.108
41	$\mathit{-}\mathit{0.001}$	0.392	$\mathit{-}\mathit{0.003}$	0.388	$\mathit{-}\mathit{0.006}$	0.819
57	0.005	0.008	0.001	0.001	0.003	0.532
77	0.011	0	0.003	0.006	0.002	0.099

Note: p-values less than 10−5 are shown as 0.

3.2.

CHO Study

3.2.1.

Optimization of Reconstruction Parameters

Figure 8 shows the two-dimensional contour plots of the AUC values as a function of the iteration number and the Butterworth postreconstruction filter cutoff frequency for the different scatter compensation methods. The plots correspond to energy window widths of 21 keV (${\mathrm{W}}_4$), 25 keV (${\mathrm{W}}_5$), 25 keV (${\mathrm{W}}_5$), and 33 keV (${\mathrm{W}}_7$) for the scatter compensation methods A, AG-DEW, AG-TEW, and AG-ESSE, respectively. We also observed similar trends for the rest of the energy windows and compensation methods. Table 5 shows the optimal cutoff frequency and iteration number for each energy window and compensation method. Figure 9 shows the short-axis images reconstructed using the different compensation methods and optimal reconstruction parameters. For all methods and energy windows, we observed that the CHO performance was best for two to three iterations. Beyond that number of iterations, the AUC values changed slowly with the iteration number. The optimal cutoff frequency was between 0.12 and $0.16{\mathrm{pixels}}^{-1}$. This was very similar to what was reported in previous studies.^4,5

Fig. 8

Two-dimensional contour plots of the AUC values as a function of the iteration number and the Butterworth filter cutoff frequency for the different compensation methods. In the above, A indicates the attenuation alone and AG-X indicates the attenuation and geometric detector response and scatter compensations using the scatter estimated with method X, where X is either DEW, TEW, or ESSE.

Table 5

Optimal reconstruction parameters for each compensation method.

Energy window width	A		AG-DEW		AG-TEW		AG-ESSE
	Iterationa	Cutoff	Iteration	Cutoff	Iteration	Cutoff	Iteration	Cutoff
1	1	0.1–0.12	1–3	0.1–0.12	2–7	0.08–0.1	1–5	0.1–0.12
9	1–2	0.12–0.16	2–3	0.12–0.14	3–7	0.1–0.12	2–3	0.12–0.14
17	2–5	0.14	2–5	0.12–0.14	2–5	0.12	2–3	0.12–0.14
21	2–5	0.14–0.6	2–5	0.12–0.16	2–5	0.12	2–3	0.12–0.14
25	2–5	0.14–0.16	2–5	0.12–0.14	2–3	0.12–0.14	2–3	0.12–0.16
29	2–5	0.14–0.6	2–5	0.12–0.14	2–5	0.12–0.14	2–3	0.12–0.16
33	2–3	0.14–0.16	3–5	0.12–0.14	2–5	0.12–0.14	2–3	0.12–0.16
37	2–3	0.14–0.16	2–5	0.12–0.14	2–5	0.12–0.14	2–3	0.12–0.16
41	2–5	0.14–0.16	2–5	0.12–0.14	2–5	0.12–0.14	2–3	0.12–0.16
57	2–7	0.14–0.16	1–5	0.12–0.14	1–7	0.12	2–3	0.12–0.16
77	2–7	0.14–0.16	1–3	0.12	1–5	0.12	2–3	0.12–0.16

^aNumber of subsets=12.

Fig. 9

Sample short-axis image corresponding to the optimal reconstruction parameters for energy windows ${\mathrm{W}}_1$ to ${\mathrm{W}}_{11}$ (from left to right) reconstructed using compensation methods A, AG-DEW, AG-TEW, and AG-ESSE (from top to bottom).

3.2.2.

Comparison of Compensation Methods

A plot of the AUC as a function of energy window width for the different compensation methods is shown in Fig. 10. Again, as was the case with the IO, the standard deviations were very small ($\sim 0.004$) and so error bars are not plotted. From these data, we see that scatter compensation using AG-DEW, AG-TEW, or AG-ESSE provided better performance than when only attenuation compensation was included. This agreed with previous studies as reported in Refs. 4 and 5. These results also indicate that ESSE scatter modeling provided better observer performance than the DEW and TEW methods. This agreed with results of previous human observer studies.⁴⁷ The observer performance, when using ESSE scatter compensation, was less sensitive to changes in the energy window width. This agreed with the results from the IO-MM when using the ESSE scatter model.

Fig. 10

Plot of AUC values for the CHO for the different energy windows and scatter compensation methods using optimal reconstruction parameters.

For each compensation method, we tested the statistical significance of the differences between the AUC value of the window that had the highest AUC value and that of the rest of energy windows. We computed the $p$-values, given in Table 6, for the hypothesis that the AUC values were different. Italicized AUC and p-value pairs indicate that the differences were not statistically significant at the level $p=0.05$. Acquisition energy window widths of 17 to 25 keV (${\mathrm{W}}_3–{\mathrm{W}}_5$) were near-optimal for attenuation compensation alone. However, when scatter compensation was incorporated into the reconstruction algorithm, the optimal energy window width was larger than for attenuation compensation only. Thus, scatter compensation enabled the use of wider energy windows. We also observed that the range of the near-optimal energy window widths changed with the compensation method. These data suggest that the acquisition energy windows of widths 21 to 29 keV (${\mathrm{W}}_4–{\mathrm{W}}_6$), 25 to 33 keV (${\mathrm{W}}_5–{\mathrm{W}}_7$), and 29 to 37 keV (${\mathrm{W}}_6–{\mathrm{W}}_8$) were near optimal for the AG-DEW, AG-TEW, and AG-ESSE compensation methods, respectively.

Table 6

Results of testing hypothesis that the highest window is optimal for various compensation methods using the CHO.

Energy window width	A	AG-DEW	AG-TEW	AG-ESSE
1	0	0	0	0
9	0.023	0	0	0
17	1	0.03	0	0
21	0.239	0.427	0.036	0
25	0.367	1	1	0
29	0	0.099	0.567	0.635
33	0	0	0.196	1
37	0.012	0	0.008	0.071
41	0	0	0.031	0.048
57	0	0	0	0
77	0	0	0	0

Note: (1) In each column, the p-values are for the difference with respect to the window with the highest AUC, indicated by the cell in bold values. (2) p-values less than 10−5 are shown as 0.

3.3.

Comparison between IO-MM and CHO

We compared the near-optimal energy window settings and the ranking of the different compensation methods obtained when using the IO-MM and the CHO, as shown in Fig. 11. We see that the IO-MM had a performance similar to the three scatter estimation methods for energy windows in the near-optimal range. However, the CHO favored the ESSE method over the DEW and TEW methods and had a similar performance when using the DEW and TEW methods. The range of near-optimal energy windows was larger for the IO-MM than the CHO when using the DEW or ESSE scatter models.

Fig. 11

Comparison of the IO-MM and CHO performances for the different compensation methods. AUC values are the averaged over the optimal range of energy window widths.

4.

Discussion

This paper provides information on the use of projection-domain observers, the IO and IO-MM, as an alternative to optimizing instrumentation and acquisition parameters in the reconstructed image domain using an anthropomorphic model observer. The IO is an observer that makes optimal use of all the information in the image data about the task. Thus, it is reasonable and has been proposed by other authors to optimize acquisition and instrumentation in terms of IO performance.^14–16 This is to ensure that there is maximum information in the data about the task. In this concept, the role of reconstruction is to put that information in a form that a human observer can best interpret. Previous work by a number of authors has consistently demonstrated that using reconstruction algorithms that model the image formation process provides equal or better performance than using reconstruction without full modeling.⁴ In other words, using imperfect models of the image formation process in the reconstruction process results in a loss of information. In cases where full modeling is used, we hypothesized that optimizing instrumentation and acquisition parameters in the projection domain and then optimizing the reconstruction algorithm (sequential optimization) would yield similar results compared with jointly optimizing the acquisition parameters and the reconstruction algorithm. However, for a variety of reasons including both practical reasons such as computation time or lack of commercial availability of reconstruction methods with full modeling of the image formation process, reconstruction with MM is often used clinically. Thus, an underlying hypothesis of this work was that, in this case, the sequential and joint optimizations yield different results. We also hypothesized that, if there is known MM, then better correlation between the sequential and joint optimization processes will be realized if the IO “discards” information in a similar way as in the reconstruction.

In the case of MPS, the diagnosis is made by human observers interpreting reconstructed images. Clearly, it is desirable that images be optimized in this domain. However, optimizing reconstructed images is computationally intensive due to the computational cost of image reconstruction and the need to optimize regularization parameters. For the calculation of the IO and IO-MM test statistics, about 6 min were required to calculate the LR of one input image. Thus, the calculations presented here required a total of 18 CPU days for each acquisition energy window and scatter modeling method. The CPU times are for a single core of a 2.33 GHz Intel Xeon E5410 quad core processor. For the CHO study, it required about 23 min to reconstruct, filter, and reorient one input image for each acquisition energy window and scatter compensation method. This required a total time of 69 CPU days. This did not include the additional time required to apply, train, and test the CHO. Thus, projection-domain optimization using the IO or IO-MM leads to reduction in the computational cost by 74% compared with image-based CHO studies.

The data in this paper indicate that projection-domain optimization using the IO can give suboptimal parameter values when MM is present. The IO suggested a broad optimal range of energy window widths between 17 and 41 keV, with the maximum AUC obtained at a width of 29 keV. On the other hand, the CHO indicated that the optimal energy window width depended on whether scatter was modeled in the reconstruction and on the scatter model used. For example, when scatter compensation was not used (i.e., using attenuation compensation only), the optimal energy window width suggested by the CHO was 17 keV. In addition, when scatter compensation was included, the difference between the optimal window for the IO and CHO was smaller. This was particularly true for ESSE, the scatter estimation method with the most accurate scatter estimates and thus the least MM. The IO-MM was partially effective in accounting for MM during the optimization: the ranges of the optimal energy window widths for the various scatter modeling methods were narrower than with the IO and very similar to those obtained using the CHO, but the IO-MM and CHO did not predict the same rankings of the scatter estimation methods.

An important question about the clinical applicability of the results of this work is with respect to the realism of the object and object variability modeled. We used a previously developed adult digital phantom population based on the 3-D XCAT phantom⁴¹ with organs’ parameters sampled from the Emory Cardiac Database.⁴² The phantom itself is based on human data and provides high levels of anatomical realism. The population consists of 54 digital phantoms including variabilities in gender, body size, heart size, and fat level with no modeling of organ textural variability included. We chose to vary anatomical parameters rather than to model other variations such as the variable uptake in organs. Due to the resolution of SPECT imaging and the nature of the defect detection task, these parameters are likely to be more important limiters of task performance than organ texture variations. In addition, data on texture variations in MPS were not available, and modeling such variations would make the fast MCMC method more difficult to use. Further, based on the relative homogeneity of myocardial tissue and the uptake mechanism of myocardial perfusion tracers, we expect the tracer distribution in normal organs, especially the heart, to be uniform.

Another limitation of this study with respect to clinical realism is that we investigated the case of detecting a known signal in a randomly varying background, whereas clinically, the signal itself is variable. The case where there is randomness in the signal, i.e., the signal is known statistically, is relatively easy to implement with the IO and IO-MM.¹⁸ However, the CHO can become a very poor observer in this situation and may not model human observer performance. We thus did not include this effect in this investigation.

5.

Conclusions

In this work, we applied three different model observers for the optimization of Tc-99m acquisition energy window width and the evaluation of different scatter estimation methods in the context of MPS defect detection. Performance was evaluated in terms of the AUC. The near-optimal range of energy window widths was broader for the IO than for the other observers. This is consistent with the fact that the IO, which has the highest performance of any observer, implicitly has perfect knowledge of the image formation process and all physical image-degrading factors. The IO-MM is an observer that allows the incorporation of nonideal forward models in the optimization process and evaluates them in terms of IO performance. We used this ability of the IO-MM to compare the DEW, TEW, and ESSE scatter estimation methods. The IO-MM had poorer performance than the IO, as expected. The results showed that the ESSE scatter estimation method provided a very similar performance to the perfect scatter model implicit in the IO. We also applied the CHO, an established surrogate for a human observer, and compared the optimal energy window settings and the ranking performance of the different scatter compensation methods to those obtained from the IO and the IO-MM. The IO-MM showed good agreement with the CHO for MPS energy window optimization and the evaluation of scatter modeling methods. The results of this study also demonstrated that the disagreement between the optimal energy window between the IO and the CHO was smallest when the MM was smallest. This suggests that the IO may be appropriate for projection-domain optimization in cases where the model of the image formation process used in image reconstruction is a good representative of the true image formation process. Together, these data indicate that appropriately taking MM into account is important in order to allow projection-domain optimization of acquisition parameters.