1.

Introduction

Objective assessment of image quality through task-specific metrics has a long history in medical imaging and is regarded by many as ultimately being the most meaningful approach to medical image evaluation.^1–4 However, the application of task-based assessment to x-ray computed tomography (CT) is recent relative to its application to planar imaging modalities and nuclear medicine. One reason for this delay is that metrics based on the Hotelling observer (HO),³ such as those considered in this work, involve the image covariance matrix, and in CT, this matrix is often extremely large (well over ${10}^9$ elements), poorly conditioned, and possesses few, if any, simplifying structural properties. In order to address the challenge of large dimensionality, efficient channels have been proposed^5–7 which essentially constitute a transformation of the image into a new basis, where the number of basis functions is substantially less than the number of image pixels. Another common means of circumventing the dimensionality problem is to assume noise stationarity, so that HO metrics can be obtained with relative computational efficiency through discrete Fourier transform (DFT) operations (see Sec. 2). Meanwhile, in order to address the somewhat unpredictable structure of the image covariance, various estimation strategies have been proposed which rely on samples of noisy images in order to construct an estimate for the image covariance when an analytic formulation of image covariance is impossible or infeasible.^8,9

Building on previous work,¹⁰ we propose a formalism, which, in certain cases, may be more appropriate than the use of efficient channels or estimation techniques, while still addressing the issues of dimensionality and structural complexity of the image covariance matrix. Our approach differs from alternative methods in that the resulting metrics are nonstochastic, and we impose no restrictions on the structure of the signal or on the nature of the noise correlations. Instead, similar to Refs. 11–12.13, we only assume that the relevant noise correlations and signal extent can be restricted to a region of interest (ROI) in the image. In order to investigate the impact and validity of various assumptions that simplify the computation of HO metrics, we herein compare the results of parameter optimization using our proposed method to two alternatives: a set of efficient [Laguerre-Gauss (LG)] channels and a DFT domain approach. Each of these approaches corresponds to a different assumption regarding the signal, the image noise, or both. We then also investigate two statistical estimation strategies for HO metrics and compare these results with our proposed method.

The specific context in which our formalism is intended to be applied is detection or classification tasks where the object of interest is small (on the order of several pixels), and the reconstruction algorithm is a direct, linear algorithm, such as filtered backprojection (FBP). In this case, we hypothesize that most of the relevant information for performing the detection or classification task is likely to be contained in pixels within a small ROI surrounding the signal. For direct linear algorithms, an analytic form of the image covariance matrix for a small ROI can be constructed and its elements can be stored directly in computer memory so long as the ROI contains at most several thousand pixels. While we consider only two-dimensional (2-D) reconstruction in this work, the extension to three-dimensional (3-D) reconstruction is straightforward so long as the corresponding 3-D ROI contains only several thousand voxels.

Since the metrics evaluated in this work are task-specific, we restrict ourselves to two relevant tasks in a specific CT application, namely dedicated breast CT. The tasks considered are microcalcification detection and the Rayleigh discrimination task, which measures the ability of an observer to distinguish between a single object and two smaller objects. The specific metric we investigate is the HO efficiency, which is a summary scalar metric describing the preservation of relevant information from the projection data measurements to the reconstructed image. The mathematical construction of this metric is described in Sec. 2, along with a description of the two imaging tasks considered. Section 3 demonstrates the use of the proposed ROI HO method for the optimization of reconstruction filter width in the FBP algorithm. Section 3 also contains results for the application of a channelized HO (CHO), a DFT-domain HO, and two estimation approaches to the same parameter optimization. Finally, a discussion and brief conclusion are included in Sec. 4.

2.

Methods

2.1.

Task Modeling

Two tasks are considered in this work. The first is a Rayleigh resolution task, wherein an observer must classify an image as either an image of a single 1.4 mm bar convolved with a Gaussian with a full width at half maximum (FWHM) of 0.53 mm or two distinct Gaussians of $\mathrm{FWHM}=0.53\mathrm{mm}$, separated by a 0.8 mm trough. These two signals were modeled on a discrete grid with pixels an order of magnitude smaller than the detector pixels backprojected to the center of the field of view (FOV). The two signals are shown in Fig. 1.

Fig. 1

The two signals simulated for the Rayleigh discrimination task.

The second task is a signal detection task, where we model a microcalcification as a Gaussian with FWHM of $80\mu \mathrm{m}$. Here the two classes into which an observer classifies the image are a signal-present class and a signal-absent class. For both tasks, we operate under the signal-known-exactly-but-variable paradigm,³ where three signal locations are considered: the center of the FOV, the position (2 cm, 2 cm), and the position (4 cm, 4 cm); however, results among the three signal locations are not substantially different.

2.3.

Model Observer Method

For each of the approaches to system optimization considered, we employ the HO (Ref. 3) and its associated metric, the HO signal-to-noise ratio (SNR). The HO quantifies the degree to which an ideal observer with full knowledge of the relevant distributions can classify an image as belonging to one of two classes (e.g., signal-absent or signal-present), using only linear operations. Let us first consider the performance of the HO in the line-integral data domain (i.e., the projection data after applying the negative logarithm). Hereafter, bold font indicates a random variable, and bars over variables indicate population averages. Given a data vector $\mathit{g}$, the HO will classify an image based on the outcome of a scalar value $\mathit{t}$, computed as

Eq. (1)

$$\mathit{t}=w_g^T\mathit{g}.$$

If we take $\mathrm{\Delta }\overline{g}$ to be the mean difference between data belonging to class 1 and data belonging to class 2, then the Hotelling template, $w_g$, is defined as

Eq. (2)

$$w_g=K_g^{-1}\mathrm{\Delta }\overline{g},$$

where $K_g$ is the data covariance matrix, whose $(i,j)$’th component is given by $\mathrm{Cov}({\mathit{g}}_i,{\mathit{g}}_j)$. We take $K_g$ to be diagonal, with diagonal elements given by²¹

Eq. (3)

$${(K_g)}_{i,i}=(e^{{\overline{g}}_i}+1)/{\overline{N}}_0,$$

where ${\overline{N}}_0$ is the average total number of incident photons in a blank scan. Further, throughout this work, we consider the data noise to be Gaussian. Since the matrix $K_g$ is clearly signal dependent, in all further discussions, we consider $K_g$ to be the average covariance matrix between the two signal classes. Under the assumption of equal class prevalence, the resulting HO is still the optimal linear observer, but may not be equivalent to the ideal observer.³ We then have that the HO SNR in the data domain (${\mathrm{SNR}}_g$) is given by

Eq. (4)

$${\mathrm{SNR}}_g^2=w_g^T\mathrm{\Delta }\overline{g}.$$

In order to evaluate the optimal reconstruction filter width for the task considered, we consider classification in the reconstructed image domain rather than in the projection data domain. In other words, the reconstructed image in vector form is given by $\mathit{y}=A\mathit{g}$. The matrix $A$ denotes the FBP reconstruction operator. Since this operator is linear and discrete-to-discrete, its action is equivalent to left-multiplication with a matrix. We then denote this matrix as $A$. The HO SNR in the image domain (${\mathrm{SNR}}_y$) is then given by

Eq. (5)

$${\mathrm{SNR}}_y^2=w_y^T\mathrm{\Delta }\overline{y},$$

where $K_y{w}_y=\mathrm{\Delta }\overline{y}$, $K_y=A{K}_g{A}^T$, $w_y$ is the Hotelling template in the image domain, and $\mathrm{\Delta }\overline{y}$ is the mean difference in the reconstructed image between the two hypotheses. Finally, we define the HO efficiency as $\epsilon ={({\mathrm{SNR}}_y/{\mathrm{SNR}}_g)}^2$.

We consider $\epsilon$ a more useful parameter for algorithm optimization since for many CT applications, one would not operate at the threshold for observer performance within which SNR is a meaningful metric. Instead, $\epsilon$ is a useful metric for algorithm and system evaluation that remains independent of the difficulty of a given task, while still retaining the exactness of computation inherent in HO performance calculation through the exact covariance matrix $K_y$. A further motivation for selection of the efficiency metric over more conventional metrics, such as SNR or area under the ROC curve (AUC),²² is that we are interested in performing system component optimization rather than full system evaluation. A single value of SNR or AUC carries little information regarding how much one can hope to improve performance by optimizing a single component. Rather, we would prefer a metric like efficiency, which relates the quality of the input to the system component to the quality of the component’s output. Here, the component that we optimize is a parameter of the reconstruction algorithm, so that the efficiency values are a reflection of how well the reconstruction preserves information relevant to the classification task.

In order to allow for direct manipulation of the image covariance matrix $K_y$, we restrict the reconstruction to an ROI image. As stated previously, 50 projection views are simulated. This angular sampling is sufficient for accurate reconstruction in the immediate vicinity of the signal; however, undersampling artifacts are visible elsewhere in the full image. We, therefore, restrict the ROI size based on the radius from the signal at which these artifacts begin to appear. Specifically, we consider the signal energy in the reconstructed image within an annulus of radius $r$ and width $2\mathrm{\Delta }r$, given by

Eq. (6)

$$E_y(r)={\int }_0^{2\pi }\mathrm{d}\theta {\int }_{r-\mathrm{\Delta }r}^{r+\mathrm{\Delta }r}{[\mathrm{\Delta }\overline{y}(r^{\prime },\theta )]}^2{r}^{\prime }\mathrm{d}{r}^{\prime }.$$

The integrand of Eq. (6) is shown in Fig. 2(a) for the microcalcification signal and a single Hanning filter width. In order to select the ROI size, we approximate the above integral at evenly spaced values of $r$ and set the ROI radius to the value of $r$ that minimizes $E_y(r)$. An example of the function $E_y(r)$ is shown in Fig. 2(b), corresponding to the image shown in Fig. 2(a).

Fig. 2

(a) The integrand of Eq. (6) for the microcalcification detection task. (b) The reconstructed microcalcification signal energy as a function of radius from the signal, $E_y(r)$. In this case, the resulting region of interest (ROI) would have dimension of $14\times 14\text{pixels}$. The results shown here correspond to 50 projection views.

For the ROI sizes investigated here, ranging from $4\times 4$ to $40\times 40$, the reconstruction matrix $A$ ranges in size from $16\times 4100$ to $1600\times 4100$. Meanwhile, the covariance matrix $K_y$ ranges from $16\times 16$ to $1600\times 1600$. The method proposed remains numerically feasible for ROI sizes up to roughly $100\times 100$. Hereafter, we refer to this approach as ROI-HO.

3.

Results

The results of applying the ROI-HO for microcalcification detection and Rayleigh discrimination are shown in Fig. 3 for a range of Hanning filter widths. The ROI-HO is noticeably sensitive to the reconstruction filter width, showing a clear maximum in performance for Hanning windows in the range of $0.75{\nu }_N$ to $0.825{\nu }_N$ for each task. However, this result should not be interpreted as giving a universally optimal filter width, but rather as a demonstration of the sensitivity of HO efficiency to relevant reconstruction algorithm parameters.

Fig. 3

Efficiency values are shown for various views and Hanning window widths (relative to the Nyquist frequency on the detector) for both the detection task (a) and the Rayleigh task (b). The Nyquist frequency in this case is ${\nu }_N=(1/2\mathrm{\Delta }u)\approx 1.3{\mathrm{mm}}^{-1}$. While the efficiency values for moderate filtering to no filtering were seen to have a dependence on ROI size, the same trend pictured here was seen for ROI sizes up to $\sim 1.5\mathrm{cm}$ diameter.

Next, we consider the use of channels for estimating HO performance in microcalcification detection. As with the ROI-HO, we restrict the CHO to an ROI, however, due to the Gaussian envelope that modulates the LG channels, this only has an effect for the smallest ROIs used. For the CHO using LG channels, results correspond solely to the microcalcification task, since the Rayleigh task involves a signal that is not radially symmetric. Figures 4 and 5 demonstrate the dependence of CHO efficiency on the number of channels and the scale factor that modulates the width of the Gaussian envelope of the LG functions. For the majority of filter widths considered, the CHO efficiency estimate is completely stable above $\sim 50$ channels, while the optimum scale factor of the channels (full width at half maximum of the Gaussian) is roughly seven times the width of the microcalcification diameter. The remainder of the results presented corresponds to these CHO parameters. In our case, however, the dependence of the CHO performance estimates upon each of these parameters is weak, so that fewer channels or a slightly different scale factor could likely produce comparable results.

Fig. 4

Dependency of channelized Hotelling observer (CHO) efficiency on number of channels for a range of reconstruction filter widths from 0.125 to 0.5 (a), and 1.0 to the (unapodized) ramp filter (b), in units of the detector Nyquist frequency. For most filter widths, the CHO efficiency stabilizes with 50 channels; however, the performance estimates for this task and system are not sensitive to the number of channels used, in general. Subsequent results shown are for 50 channels.

Fig. 5

CHO efficiency using 50 channels and 50 projection views for a range of Gaussian envelope widths in the Laguerre-Gauss channels. The $x$ axis is normalized to the $80\mu \mathrm{m}$ microcalcification diameter. A scale factor of roughly seven times the microcalcification width was determined to be optimal. As in Fig. 4, (a) and (b) correspond to narrower and wider filters, respectively.

Figure 6 compares three approaches to optimization of the reconstruction filter width parameter, namely the ROI-HO approach proposed in this work, the CHO with 50 LG channels, and the DFT-domain approach (labeled FHO). Since each of the methods utilizes the same ROI images, the ROI-HO that computes the exact HO performance should be taken as the ground truth for the sake of evaluating the alternative approaches. Clearly, for microcalcification detection, the CHO is an excellent approximation of the true HO. Further, the CHO allows for greater computation efficiency, since, in general, the ROIs used contained $>50\text{pixels}$. However, the efficiency of this CHO comes with a loss of generality, as it is not applicable to general signals lacking radial symmetry, as in the Rayleigh discrimination task.

Fig. 6

A comparison of optimization of reconstruction filter width for 50 projection views. Shown are results obtained with the proposed ROI-HO, an approximate HO computed using the discrete Fourier transform (labeled FHO), and the CHO. The curves in (a) correspond to the detection task, while (b) corresponds to the Rayleigh task.

Meanwhile, the DFT-domain HO was capable of reproducing the general trend of the HO results; however, the quantitative reliability of this method varied widely depending on the specific parameter setting used. For example, while wider filter widths allowed for accurate computation of detection task performance with the DFT-domain approach, the estimates of Rayleigh task performance using this method could vary from the HO by as much as 40%, as seen in the right side of Fig. 6. While this may not impact the optimization of filter width substantially, if absolute HO performance with a fixed filter width of, say, $0.625{\nu }_N$ were of interest, the DFT-domain HO could be misleading.

Figure 7 demonstrates the use of noisy sample images to estimate HO efficiency when the mean image under each hypothesis is known. In our case, since we consider only linear image reconstruction algorithms, the mean image is produced simply by reconstructing an image without noise. The left and right plots correspond, respectively, to the use of independent noisy data sets for each filter width and the reuse of the same simulated data for each filter width. The variability seen in the left-hand figure provides intuition as to the variance of the efficiency estimates, while the figure on the right illustrates that, although subject to the same variability, reusing the same data realizations correlates the estimates for different filter widths, potentially allowing for more reliable rank ordering of parameter settings. In other words, the estimated curves on the right could undergo vertical translation due to the variance of the estimator, but are more likely to preserve their shape than the curves shown on the left.

Fig. 7

Results of estimating the HO efficiency from sample images using the method proposed in Ref. 9 with 300 and 700 noisy sample images. (a) An independent set of images is used for each filter width. (b) The same noisy data realizations are used for each filter width. Error bars corresponding to two standard deviations of the estimates for 700 sample images are shown in Fig. 8.

Fig. 8

Results from Fig 7 for 700 noisy sample images are shown with error bars corresponding to two standard deviations derived from jackknifing. The error bars illustrate statistical variation, but do not account for inherent estimator bias, which can be seen by comparing with the analytically computed ROI-HO curves. It is possible that other sources of uncertainty not investigated here, aside from statistical variations, could influence these results. As in Fig. 7, independent sets of images were used in (a), while the same set of images was used for each point in (b).

Figure 9 illustrates the use of linear classifier training and testing in order to estimate HO performance. In this case, prior knowledge of the mean images is not exploited, and sample estimates of the images are computed instead. The left-hand plot corresponds to the hold-out approach, where independent image sets are used for the training and testing phases, while the right-hand plot is generated using resubstitution. As seen in the figure, these two approaches introduce negative and positive biases in the estimates, respectively. Note that thousands of images are required in order to construct quantitatively meaningful estimates. However, the general trend of the results can be seen for comparatively few samples. This suggests that for certain tasks, such as rank-ordering of only a few options for algorithm implementation, the training and testing approach could be useful if one lacks an adequate model of the image covariance or class means.

Fig. 9

(a) HO ${\mathrm{SNR}}^2$ estimates from training and testing performed using the hold-out approach for 500, 1000, and 1500 training images, with an equal number of testing images. The prevalence of images from each class is also equal. (b) HO ${\mathrm{SNR}}^2$ estimates resulting from training and testing performed using resubstitution for 1000, 2000, and 3000 total images. The bias and variance of the estimates is worst for narrow filter widths, where the size of the ROI used is largest. Variance of the estimates is illustrated in Fig. 10 through 95% confidence intervals derived from bootstraping.

Fig. 10

Shown here are the results from Fig. 9 for the largest number of sample images with error bars denoting 95% confidence intervals derived from 1000 bootstrap samples. These errors derive from the variance of the training and testing estimator, while inherent bias in the estimator contributes additional error, especially for small filter widths where the number of ROI pixels is large. As in Fig. 8, the error bars here only convey statistical uncertainties, and bias can be inferred by comparison with the analytically computed ROI-HO curves. The results in (a) correspond to the hold-out approach, while the results in (b) correspond to resubstitution.

4.

Conclusions

We have demonstrated that for classification tasks involving small signals in CT, the HO performance computed within an ROI constitutes an efficient and objectively meaningful approach to reconstruction algorithm and system optimization. Specifically, we have employed the HO to optimize reconstruction algorithm filter width for two tasks: microcalcification detection in dedicated breast CT and Rayleigh discrimination.

We performed a comparison of this approach with several alternatives for objective assessment. In broad terms, these alternatives fall into two categories: methods based on computation of statistical estimates and methods based on approximations regarding the signal and/or image noise. While our proposed methodology assumes that information beyond a small ROI is irrelevant to the classification task, the CHO imposes assumptions of rotational symmetry on the signal and Hotelling template, and the Fourier-domain approach assumes stationary noise. Strictly speaking, none of these assumptions are valid in CT for the tasks considered here; however, each approach possesses strengths and weaknesses. The Fourier-domain HO is the most computationally efficient method for large ROIs, but is less accurate than the ROI-HO or CHO approach. Meanwhile, the CHO with LG channels is limited to cases when the signal of interest and Hotelling template possess radial symmetry.

The statistical estimation approaches investigated possess an attractive degree of flexibility, in that they do not rely on any assumptions regarding the image covariance matrix or, in the case of the training/testing approach, the signals themselves. However, this flexibility comes at the expense of computational efficiency, as this family of approaches requires excessively large numbers of noisy samples, making any extensive system or algorithm optimization prohibitively time-consuming. Finally, while we have illustrated some of the trade-offs, benefits, and shortcomings of several methods for objective assessment in CT, ultimately future work will be necessary to ensure that the proposed methodology yields metrics that correlate with improved performance of humans for the given tasks.