Spatial resolution and noise tradeoffs in pinhole imaging system design: a density estimation approach

Jeffrey A. Fessler

doi:10.1364/OE.2.000237

1. Introduction

The design of imaging systems and image recovery algorithms generally involves tradeoffs between spatial resolution and noise. For example, in a simple pinhole imaging system, a larger pinhole allows more photons to pass through, which reduces the relative uncertainty of the measurements, but at a price of degraded spatial resolution. The problem of specifying system parameters such as pinhole size is therefore frequently encountered in the system design process. This paper considers the image recovery problem as an indirect density estimation problem, and considers the following design criterion: minimize the variance of the object estimate subject to a prespecified object spatial resolution. We show analytically that the variance-minimizing spatial resolution of the imaging system is proportional to the desired spatial resolution of the object estimate, when the kernel of the density estimator is based on an apodized inverse filter. This is an intuitive relationship, but one that has not been previously established theoretically to our knowledge.

There are a variety of methods that have been proposed for “optimal” choice of system parameters in imaging system design. Each such design method has its own merits and limitations, and it is unlikely that any single design method will be universally accepted as the canonical choice. Since imaging systems are often built to serve multiple purposes, system designers can benefit from exploring multiple design criteria. We make no pretense that the criterion analyzed in this paper is always preferable over alternatives, but we believe that it is a potentially useful addition to the collection of tools available for imaging system design.

One very principled approach to imaging system design is to optimize the system for the performance of a certain task or collection of tasks, e.g. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. In the context of detecting a known Gaussian signal in a stationary nonuniform background, Myers et al.[3] found that the optimum aperture size was fairly close to the Gaussian signal width. One can evaluate and optimize task performance with respect to system parameters using either human observers or machine observers. When imaging systems are designed for specific tasks, such as detecting myocardial perfusion defects [11], task performance is a natural metric for design. Often imaging systems must serve multiple purposes, so more generic measures of performance, such as spatial resolution and noise, are useful to complement task-specific performance measures. Manufacturers of medical imaging instruments typically report only spatial resolution and sensitivity, despite their indirect relationships to task performance.

Another approach to analyzing system performance is the Cramer-Rao (CR) bound. The ordinary CR bound is applicable only to unbiased estimators, which limits its utility in imaging problems where bias is typically inevitable. The uniform CR bound [12] is a recent extension of the CR bound that allows for biased estimators. The uniform CR bound provides the minimum achievable variance of an estimator whose bias-gradient length is below a specified threshold. Although the bias gradient is related to spatial resolution in some cases [13], in general it is currently fairly challenging to interpret the tradeoff between variance and bias gradient length. In particular, we have observed some counter-intuitive results concerning optimal collimator resolution as a function of target image resolution, perhaps in part due to nonlinear and system-dependent relationships between bias gradient length and spatial resolution [14, 15]. Our intuition is that as one reduces the required reconstructed spatial resolution, the variance-minimizing collimator size should increase¹ This intuitive relationship has not always been apparent in our uniform CR bound experiments, which motivated the work described in this paper.

One appeal of the uniform CR bound is that it is estimator independent. However, data from any system must eventually be reconstructed by some estimator, and the class of reasonable estimators is arguably fairly small. So as a part of exploration of system performance, it is sensible to also investigate resolution/noise tradeoffs for broad classes of estimators, albeit without the full generality of the uniform CR bound.

Another difficulty with CR bounds is that they (apparently) require an inherently discrete formulation both for the detector space (which is often but not always natural) and for the image space (which is somewhat unnatural since emission distributions are continuous entities). The discrete formulation leads to large matrix inversion problems, and, as shown in [15], challenges in interpretation due to differences in performance even for neighboring pixels depending on “small” discretization effects.

In this paper, we adopt a completely continuous formulation. The only discretization is at the final step (numerical integration), which is fundamentally different from initially formulating a discrete problem. The treatment is closely related to “indirect” density estimation [16, 17]. A good reference on direct density estimation is [18]. Other publications that are relevant to density estimation an image reconstruction include [19, 20, 21, 22, 23, 24].

Section 2 describes the problem generally. Section 3 describes the density estimation approach and analyzes it statistically. Section 4 focuses on the shift-invariant case, considers a specific kernel for the density estimator based on an apodized inverse filter, and derives analytic results for specific pinhole shapes. Section 5 reports numerical simulations that confirm the analysis.

2. Problem

Consider an emitting object with emission-rate density λ(x̲) having units emissions per unit time per unit volume. The emission-rate density λ(x̲) is defined over a subset Ω of IR^d, where typically d = 2 for planar imaging and d = 3 for volumetric imaging. We assume that the time-ordered sequence of emissions originate from statistically independent random spatial locations {X̲ ₁, X̲ ₀,…} drawn from a Poisson spatial point process [25]. In particular, the joint probability that the first n ₀ emissions originate in any (measurable) regions B_n ⊆ Ω is given by²:

P [\cap_{n = 1}^{n_{0}} {{\underline{X}}_{n} \in B_{n}}] = \prod_{n = 1}^{n_{0}} P [{\underline{X}}_{n} \in B_{n}] = \prod_{n = 1}^{n_{0}} \frac{\int_{B_{n}} λ (\underline{x}) d \underline{x}}{\int_{Ω} λ (\underline{x}) d \underline{x}} .

Spatial locations x̲ ∈ Ω over which the emission-rate density λ(x̲) has relatively greater values are the “hot regions” of the object.

For any emission imaging system, not all emitted photons are detected. Let s(x̲) denote the sensitivity function of the emission system, i.e., s(x̲) is the probability that a photon emitted from location x̲ is detected (somewhere) by the system. Then when the system detects an emission, the probability (density³) that that emission originated from spatial location x̲ is given by

f (\underline{x}) = \frac{λ (\underline{x}) s (\underline{x})}{\int λ (\underline{x}') s (\underline{x}') d \underline{x}'} = \frac{λ (\underline{x}) s (\underline{x})}{r},

where

r ≜ \int λ (\underline{x}) s (\underline{x}) d \underline{x}

is the total rate of detected events (with units detected counts per unit time⁴).

Unfortunately, emission imaging systems never observe the emission locations {X̲_n } directly. Instead, the nth emitted photon is detected by a position-sensitive measurement device, which records a position V̲_n . (The detector may also record other event attributes such as energy, and our formulation allows for this generality, but for simplicity one can think of V̲_n as position.)

For a planar emitting object imaged by an ideal planar detector through an ideal pinhole located at the center of the transverse coordinates, the recorded spatial locations would be related to the locations of the emissions through the simple relationship V̲_n = mX̲_n , n = 1,2,…, where m is the (negative) source magnification factor [26]. In this hypothetical case one could exactly recover the emission location from the measured event positions via the simple relationship ${\underset{̲}{X}}_{n} = \frac{1}{m} {\underset{̲}{V}}_{n}$ . Given the recorded positions of many such photons V̲ ₁,…, V̲_N , and therefore the positions of many emissions, one could estimate the density f(x̲) by a variety of well-known methods for density estimation, such as a simple kernel estimate:

\hat{f} (\underline{x}) = \frac{1}{N} \sum_{n = 1}^{N} \frac{1}{β} k (\frac{\underline{x} - \frac{{\underline{V}}_{n}}{m}}{β}) = \frac{1}{N} \sum_{n = 1}^{N} \frac{1}{β} k (\frac{\underline{x} - {\underline{X}}_{n}}{β}),

where k is a nonnegative 2D kernel function (e.g. a Gaussian kernel) that integrates to unity [18]. This type of problem is called “direct” density estimation, since the measurements {X̲_n } are drawn directly from the density f(x) that we wish to estimate. The “bandwidth” parameter β controls the tradeoff between spatial resolution and “noise” (variance of f̂). When more detected events are available, one typically uses narrower kernels [18]. The problem of choosing the bandwidth for kernel density estimation for direct observations is well studied, and data-driven methods that are efficient in terms of mean squared-error are available [27]. However, in the context of image recovery, the mean squared-error metric, which equally weights bias and variance, may not be the appropriate loss function.

An ideal pinhole does not allow very many photons to pass, so in practice one must use a finite-sized pinhole. Furthermore, the position-sensitive detector does not record the exact position of the incident photon, but rather a noisy version thereof. In general the recorded positions {V̲_n } are only indirectly related to the emitted positions {X̲_n } through some conditional pdf:

f (\underline{v} | \underline{x})

which describes the “infinitesimal probability” that an emission at location x̱ that is detected will be recorded at detector position v̲. The pdf f(v̲|x̲) includes both the pinhole collimator response function as well as the detector response function. A maximum smoothed likelihood approach to the indirect density estimation problem of estimating f(x̲) from the measurements {V̲_n } is considered in [28], although without analyzing the variance or spatial resolution of the estimator.

Note that since f(v̲|x̲) is a conditional pdf, it integrates to unity over v̲. We assume that this conditional pdf holds for all detected events, i.e.

f ({\underline{v}}_{1}, {\underline{v}}_{2}, \dots ∣ {\underline{x}}_{1}, {\underline{x}}_{2}, \dots) = \underset{n}{Π} f ({\underline{v}}_{n} ∣ {\underline{x}}_{n}) .

This is reasonable assumption, except possibly at high count rates when deadtime factors and pulse pile-up effects are significant, e.g. [29].

2.1 The estimation problem

Suppose the imaging system records a total of N events during a prespecified period of time t ₀. By assumption, N is a Poisson random variable with mean

E [N] = t_{0} \int λ (\underline{x}) s (\underline{x}) d \underline{x} = t_{0} r .

Note that for a pinhole system the sensitivities s(x̲) of the system increase with pinhole size, and thus so does the expected number of recorded events. For each of these N events the system records independent and identically distributed position attributes V̲ ₁,…,V̲_N that each have the following marginal pdf

f_{V} (\underline{v}) = \int f (\underline{v} ∣ \underline{x}) f (\underline{x}) d \underline{x}

by total probability. This is a list-mode formulation [30, 31].

We would like to estimate λ(x̲) from the observed random variables N and {V̲_n $}_{n = 1}^{N}$ . We assume that t ₀, s(x̲), and f(v̲|x̲) are known, i.e., previously determined by some combination of modeling and system measurements. If we can find an estimate f̂(x̲) of f(x̲), then we can also easily estimate λ(x̲). Combining (1) and (3) we see that

λ (\underline{x}) = \frac{f (\underline{x})}{s (\underline{x})} r = \frac{f (\underline{x})}{s (\underline{x})} \frac{E [N]}{t_{0}} .

Thus a natural estimator for the emission-rate density λ(x̲) is simply

\hat{λ} (\underline{x}) = \frac{\hat{f} (\underline{x})}{s (\underline{x})} \frac{N}{t_{0}} .

We now turn to the problem of finding a suitable estimator f̂(x̲). This is called an indirect density estimation problem, since the observed measurements {V̲_n } are only indirectly related to the density f(x̲) of interest through (4).

3. Kernel-based indirect density estimator

In this paper, we consider the following class of kernel-based indirect density estimators⁵:

\hat{f} (\underline{x}) = \frac{1}{N} \sum_{n = 1}^{N} g_{β} (\underline{x}, {\underline{V}}_{n}),

where g_β (x̂, v̲) is a user-defined function that typically partially inverts the blurring caused by the system response f(v̲|x̲). A concrete example of g_β is given in (22) below. The function g_β depends on a user-selected parameter β that determines the spatial resolution of f̂(x̲). For f̂(x̲) to be a valid pdf, it must integrate to unity. Therefore the function g_β (x̲, v̲) should integrate to unity over x̲:

1 = \int g_{β} (\underline{x}, \underline{v}) d \underline{x} .

In direct density estimation, one usually chooses kernel functions k(∙) in (2) that are nonnegative, since f(x̲) must be nonnegative, although it is possible to reduce bias by allowing kernels with negative values [18, p. 66]. In the context of indirect density estimation, the function g_β generally must contain negative values in order to partially deconvolve the blur in f(v̲|x̲). In the context of image reconstruction, one can think of the estimator (6) as an event-by-event backprojector⁶, where the backprojector includes the ramp filter, apodizer, etc. Note that estimators in this class are probably suboptimal since they treat all photons equally. Nevertheless it is a useful class for examining resolution/noise tradeoffs.

Combining (6) with (5), the corresponding estimator for the emission-rate density is

\hat{λ} (\underline{x}) = \frac{\hat{f} (\underline{x})}{s (\underline{x})} \frac{N}{t_{0}} = \frac{1}{t_{0} s (\underline{x})} \sum_{n = 1}^{N} g_{β} (\underline{x}, {\underline{V}}_{n}) .

This emission-rate density estimate is a function defined for all x̲. There are no pixels or voxels involved, which simplifies the analysis. The effective number of “degrees of freedom” is determined by N and by β. In the following we examine the statistical properties of the above estimator $\hat{λ}$ (x̲).

3.1 Mean function

The mean function for the estimator $\hat{λ}$ (x̲) is derived as follows⁷:

μ (\underline{x}) = E [\hat{λ} (\underline{x})]

Thus we have the following linear relationship between the estimator mean and the true emission density:

μ (\underline{x}) = \int psf (\underline{x}, \underline{x}') λ (\underline{x}') d \underline{x}'

where

psf (\underline{x}, \underline{x}') ≜ \frac{s (\underline{x}')}{s (\underline{x})} \int g_{β} (\underline{x}, \underline{v}) f (\underline{v} | \underline{x}') d \underline{v}

is effectively the overall point-spread function (PSF) for the combined image acquisition / reconstruction process. This PSF depends on the system response, which is contained in f(v̲|x̲), as well as the regularization in the reconstruction algorithm, which is contained in g_β . Equations (9) and (10) are space-varying generalizations of equation (10.32) in Barrett and Swindell’s text [32] for the mean of a filtered Poisson point process. If one uses a g_β function that has negative values, then the PSF may also have negative values. The reconstructed spatial resolution is controlled by the PSF (10), so for good spatial resolution, g_β must partially “deconvolve” any blur caused by f(v̲|x̲).

3.2 Second-moment functions

Before computing the autocorrelation function of $\hat{λ}$ , we first note that since N is Poisson,

E [N^{2}] = Var {N} + {(E [N])}^{2} = E [N] + {(E [N])}^{2} = t_{0} r + {(t_{0} r)}^{2}

so E[N ² - N] = (t ₀ r)². Then from (7), the autocorrelation function for $\hat{λ}$ (x̲) is derived as follows:

R_{\hat{λ}} ({\underline{x}}_{1}, {\underline{x}}_{2}) = E [\hat{λ} ({\underline{x}}_{1}) \hat{λ} ({\underline{x}}_{2})]

Therefore the autocovariance function for $\hat{λ}$ is

K_{\hat{λ}} ({\underline{x}}_{1}, {\underline{x}}_{2}) = E [\hat{λ} ({\underline{x}}_{1}) \hat{λ} ({\underline{x}}_{2})] - μ ({\underline{x}}_{1}) μ ({\underline{x}}_{2})

To simplify, note that

E [g_{β} ({\underline{x}}_{1}, \underline{V}) g_{β} ({\underline{x}}_{2}, \underline{V})] = \int g_{β} ({\underline{x}}_{1}, \underline{v}) g_{β} ({\underline{x}}_{2}, \underline{v}) f_{V} (\underline{v}) d \underline{v}

so the autocovariance function is

K_{\hat{λ}} ({\underline{x}}_{1}, {\underline{x}}_{2}) = \frac{r}{t_{0} s ({\underline{x}}_{1}) s ({\underline{x}}_{2})} \int [\int g_{β} ({\underline{x}}_{1}, \underline{v}) g_{β} ({\underline{x}}_{2}, \underline{v}) f (\underline{v} | \underline{x}') d \underline{v}] s (\underline{x}') λ (\underline{x}') d \underline{x}' .

In particular, the variance function is

σ^{2} (\underline{x}) ≜ Var {\hat{λ} (\underline{x})} = K_{\hat{λ}} (\underline{x}, \underline{x}) = \frac{1}{t_{0} s^{2} (\underline{x})} \int [\int g_{β}^{2} (\underline{x}, \underline{v}) f (\underline{v} | \underline{x}') d \underline{v}] s (\underline{x}') λ (\underline{x}') d \underline{x}' .

(This equation is a space-variant generalization of (10.31) in [32].) Note that the variance depends inversely on the scan time t ₀, which is expected.

For a specific imaging system f(v̲|x̲), object λ(x̲), and reconstruction method g_β of interest, one could compute (9) and (12) for a range of β values or pinhole sizes to investigate the resolution/noise tradeoff. The computational tractability of such evaluations will depend on the complexity of f(v̲|x̲) and g_β . To obtain insight into the tradeoffs, we consider the simpler shift-invariant case in the remainder of this paper.

4. Shift-invariant case

Suppose the system is shift-invariant, i.e. f(v̲|x̲) = h(v̲ - x̲), where for example h is a normalized pinhole response function. A pinhole that is mechanically scanned over the emitting object is an example of a shift-invariant system⁸. Such systems have been used for many years [26] and continue to find specialized applications, e.g. [33]. Note that since f(v̲|x̲) is a pdf, it must integrate to unity over v̲, so we must also have h integrate to unity. Suppose also that the reconstruction algorithm is shift invariant, i.e. g_β (x̲,v̲ ) = g_β (x̲ - v̲) (with a slight notation abuse/reuse). Finally, assume that the sensitivity is also space-invariant, i.e. s(x̲) = s ₀ for some positive constant s ₀. Then the above expressions simplify as follows.

The mean expression (9) becomes:

μ (\underline{x}) = \int [\int g_{β} (\underline{x}, \underline{v}) f (\underline{v} ∣ \underline{x}') d \underline{v}] λ (\underline{x}') d \underline{x}'

where x̲′′ = v̲ - x̲′ and ∗ denotes d-dimensional convolution. Thus we have the following convolution relationship (cf (10.11) of [32]):

μ = g_{β} * h * λ,

i.e., the estimator mean is the convolution of the underlying emission-rate density with the system PSF h(∙) and the recovery kernel g_β (∙). Therefore the spatial resolution is controlled by

psf (\underline{x}) ≜ (g_{β} * h) (\underline{x}),

with corresponding frequency response or overall transfer function

PSF (\underline{u}) ≜ G_{β} (\underline{u}) H (\underline{u}),

where F(u̲) ≜ ∫ f(x̲)e ^{-i2πu̲∙x̲} d x̲ denotes the d-dimensional Fourier transform of f(x̲).

Similarly, the inner variance term in (12) becomes:

\int g_{β}^{2} (\underline{x}, \underline{v}) f (\underline{v} ∣ \underline{x}') d \underline{v} = \int g_{β}^{2} (\underline{x} - \underline{v}) h (\underline{v} - \underline{x}') d \underline{v}

which is equivalent to the “noise kernel” of (10.35) of [32]. Thus, in the shift-invariant case the variance function (cf (10.10) of [32]) simplifies to

σ^{2} (\underline{x}) = \frac{1}{t_{0} s_{0}} \int (g_{β}^{2} * h) (\underline{x} - \underline{x}') λ (\underline{x}') d \underline{x}'

Therefore in the shift-invariant case it is straightforward to compute variances (approximately) using FFTs to calculate the convolutions.

4.1 Spatially smooth objects λ(x̲)

If the object is spatially smooth, i.e. the scale of the spatial variations of (h ∗ λ)(x̲) is large relative to the support of $g_{β}^{2}$ (x̲), then the variance expression simplifies as follows.

t_{0} s_{0} σ^{2} (\underline{x}) = (g_{β}^{2} * h * λ) (\underline{x})

where we define $\tilde{λ}$ ≜ h ∗ λ (cf (10.41) of [32]). For small β or for spatially smooth objects this approximation should be fairly accurate⁹. For pinhole imaging, gβ is a real function, so by combining the above approximation with Parseval’s theorem:

σ^{2} (\underline{x}) \approx \frac{\tilde{λ} (\underline{x})}{t_{0} s_{0}} \int g_{β}^{2} (\underline{x'}) d \underline{x'} = \frac{\tilde{λ} (\underline{x})}{t_{0} s_{0}} \int {∣ G_{β} (\underline{u}) ∣}^{2} d \underline{u} .

This is a very tractable approximation to the estimator variance.

4.2 Resolution-noise tradeoffs

In general, both the sensitivity s₀ and the overall transfer function PSF(u̲) = G_β (u̲)H(u̲) depend on the pinhole size. Therefore the expression (17) does not immediately provide the optimal choice for the pinhole size. In the following we consider a specific class of choices for g_β (∙), and show that the variance-minimizing pinhole size is proportional to the specified reconstructed spatial resolution.

The relationships (15) and (17) epitomize the resolution-noise tradeoff. For good spatial resolution in (15), we would like G_β (u̲) ≈ 1/H(u̲), but if H(u̲) is small, then such a G_β (u̲) is large, which amplifies the variance term in (17).

4.3 Apodized inverse filter

Consider a general pinhole with transmissivity function t(x̲) ≥ 0, which we assume is normalized so that ∫t(x̲)d x̲= 1. Let T(u̲) be the d-dimensional Fourier transform of t(x̲). The design problem is to choose the pinhole size w, where the normalized pinhole response function is defined by $h (\underset{̲}{x}) = \frac{1}{w^{d}} t (\frac{\underset{̲}{x}}{w})$ for which H(u̲) = T(wu̲). Define the apodized inverse filter

G_{β} (\underline{u}) ≜ \frac{A (β \underline{u})}{H (\underline{u})} = \frac{A (β \underline{u})}{T (w \underline{u})},

where A(βu̲) is a user-chosen apodizing function which we assume to be real and symmetric. Without loss of generality, we assume A(u̲) has been defined so that the FWHM of a(x̲) has unit length. From (15), the overall transfer function of this system is

PSF (\underline{u}) = G_{β} (\underline{u}) H (\underline{u}) = A (β \underline{u}),

so the overall PSF is simply psf(x̲) = β ^-d a(x̲|β). Therefore the FWHM of the overall PSF is precisely β for this estimator for any pinhole size w. We now show that the variance-minimizing choice for the pinhole width w is directly proportional to β.

We assume s ₀ = c ₀ W^p for some constant c ₀ independent of w and for some power p > 0. Typically p = d; for example, the sensitivity of a circular pinhole is proportional to its area, which is proportional to w ². From (17) the variance is approximately:

σ^{2} (\underline{x}) \approx \frac{\tilde{λ} (\underline{x})}{t_{0} c_{0} w^{p}} \int {∣ G_{β} (\underline{u}) ∣}^{2} d \underline{u}

where z̲ ≜ wu̲ and c ₁ ≜ $\tilde{λ}$ (x̲)/t ₀ c ₀. To find the pinhole width w that minimizes the variance, we zero the partial derivative of the variance σ ² with respect to the width w:

0 = \frac{- (p + d) c_{1}}{w^{p + d + 1}} \int {∣ T (\underline{z}) ∣}^{- 2} A^{2} (\frac{\underline{z} β}{w}) d \underline{z}

or, by defining α = β/w:

0 = \int \frac{\frac{p + d}{2} A^{2} (α \underline{z}) + A (α \underline{z}) \dot{A} (α \underline{z}) α \underline{z}}{{∣ T (\underline{z}) ∣}^{2}} d \underline{z} .

The above equality depends only on the ratio α = β/w. So if there is a root α ₀ > 0 that corresponds to a global minimizer of the variance σ ², then the variance-minimizing w is proportional to the reconstructed spatial resolution β through the relationship W _min = $α_{0}^{- 1}$ β.

4.4 Relationship to sieves

The above apodized inverse filter is closely related to the method of sieves for density estimation [34], in the sense that $\hat{λ}$ (x̲) is an unbiased estimate of β^-da(x̲| $\underline{β}$ ) ∗ λ(x̲).

4.5 Gaussian pinhole example

Figure 1. Profile through an approximate Gaussian pinhole.
Download Full Size | PDF

As a concrete example, consider the Gaussian pinhole illustrated in Fig. 1 for d = 2 dimensional imaging. To simplify notation, define r = ∥x̲∥ and ρ = ∥u∥ for circularly symmetric 2D imaging.

The exact transmissivity of this aperture is

τ_{w} (r) = {\begin{matrix} e^{- μl {(\frac{r}{r_{b}})}^{2}}, r \leq r_{b} \\ e^{- μl}, r \geq r_{b} \end{matrix} .

However, if μl is sufficiently large, then we can approximate this transmissivity by

τ_{w} (r) = \exp (- π {(\frac{k}{w} r)}^{2})

where w is the FWHM of the pinhole response (i.e. τ_w (w/2) = 1/2) and

κ ≜ 2 \sqrt{\frac{\ln 2}{π}} .

The sensitivity of this pinhole is therefore

s_{w} = \int τ_{w} (∣ ∣ \underline{x} ∣ ∣) d \underline{x} = {(\frac{w}{κ})}^{d},

which is proportional to w^d as expected. The normalized transmissivity (for unit pinhole width w = 1) is

t (r) = \frac{τ_{1} (r)}{s_{1}} = κ^{d} e^{- π {(kr)}^{2},}

with corresponding frequency response

T (ρ) = e^{- π {(\frac{ρ}{κ})}^{2}} .

We choose a Gaussian apodizing function A(u̲) = e ^-π(ρ/κ)² so that the PSF corresponding to A(βρ) has FWHM β. The corresponding recovery filter is thus

G_{β} (\underline{u}) = \frac{A (βρ)}{T (wρ)} = \frac{e^{- π {(\frac{βρ}{κ})}^{2}}}{e^{- π {(\frac{wρ}{κ})}^{2}}} = \exp (- π {(\frac{ρ \sqrt{β^{2} - w^{2}}}{κ})}^{2}),

with corresponding space-domain recovery kernel

g_{β} (\underline{x}) = \frac{κ^{2}}{β^{2} - w^{2}} \exp (- π {(\frac{rκ}{\sqrt{β^{2} - w^{2}}})}^{2}) .

Substituting A(∙) into (19), the variance function is approximately

σ^{2} (\underline{x}) \approx \tilde{λ} (\underline{x}) \frac{κ^{d}}{t_{0} w^{2 d}} \int \int e^{π 2 {(\frac{ρ}{κ})}^{2}} \exp (- π 2 {(\frac{ρβ}{wκ})}^{2}) d \underline{u}

Figure 2. Standard deviation of estimate versus Gaussian pinhole width.
Download Full Size | PDF

where we have applied Parseval’s theorem in conjunction with the Hankel transform to evaluate the integral. Note that we must have w < β for the integral to be finite, i.e. the pinhole width must be no larger than the desired spatial resolution. Figure 2 plots the variance versus pinhole width w. Differentiating the variance with respect to w and zeroing yields the following relationship:

w_{min} = \frac{β}{\sqrt{2}} .

Taking the second derivative confirms that this is the variance-minimizing choice. A plot of the variance as a function of pinhole width is shown in Figure 2.

Therefore, we have shown that for a Gaussian pinhole imaging system and an apodized inverse filter reconstruction method, the variance-minimizing pinhole width is proportional to the desired reconstructed spatial resolution.

5. Laplacian pinhole example

Figure 3. Profile through an approximate Laplacian pinhole.
Download Full Size | PDF

A somewhat more conventional pinhole is the Laplacian¹⁰ pinhole shown in Fig. 3. The exact transmissivity of such a pinhole is

τ_{ω} (r) = {\begin{matrix} e^{\frac{- μ l r}{r_{b}}}, r \leq r_{b} \\ e^{- μ l}, r \geq r_{b} . \end{matrix}

If μl is sufficiently large, we can approximate this transmissivity by

τ_{w} (r) = e^{\frac{- γr}{w}},

where γ = 2 log 2 and w is the FWHM of the pinhole response. The sensitivity of this pinhole is

s_{w} = \int τ_{w} (∣ ∣ \underline{x} ∣ ∣) d \underline{x} = 2 π {(\frac{w}{γ})}^{2},

which is also proportional to w ². The normalized transmissivity is

t (r) = \frac{τ_{1} (r)}{s_{1}} = \frac{γ^{2}}{2 π} e^{- γr}

which has corresponding frequency response [35]

T (ρ) = \frac{γ^{3}}{{[{(2 πρ)}^{2} + γ^{2}]}^{\frac{3}{2}}} .

We again choose a Gaussian apodizing function A(ρ) = e ^-π(ρ/κ)² where κ was defined in (20), so the PSF again has FWHM β. Substituting into (19), the estimate variance is approximately

σ^{2} (\underline{x}) \approx \frac{\tilde{λ} (\underline{x})}{\frac{t_{0} 2 π w^{2}}{γ^{2}}} \int_{0}^{2 π} \int_{0}^{\infty} \frac{1}{γ^{6}} {[{(2 πρ)}^{2} + γ^{2}]}^{3} \exp (- π 2 {(\frac{ρβ}{wκ})}^{2}) ρ dρ dθ

After some tedious integration, we arrive at

σ^{2} (\underline{x}) \approx \frac{\tilde{λ} (\underline{x})}{t_{0}} \frac{1}{2 {(βκ)}^{4}} [6 y^{2} + 6 y + 3 + y^{- 1}] where y = 2 π {(\frac{w}{κλβ})}^{2} .

The variance-minimizing pinhole size can be found numerically to be w _min ≈ 0.3β.

The variance for the Laplacian pinhole is also plotted in Fig. 2. The minimum standard deviation for the Gaussian pinhole is about 2.4 times lower than that of the Laplacian pinhole, presumably because the Gaussian pinhole is better matched to the Gaussian-apodized inverse filter.

6. Simulation results

Since the analytical development for the variance-minimizing pinhole width involved approximations, we performed Monte Carlo simulations to evaluate the empirical variance of the estimators as a function of pinhole size.

We used a 1D object of the form

λ (\underline{x}) \propto 9 δ (x - 146) + rect (\frac{(x - 208)}{64} + 2 Λ \frac{(x - 64)}{44}),

where Λ(x) = (1 - |x|) rect(x/2) is the unit triangular function. The desired reconstructed spatial resolution was arbitrarily chosen to be β = 3 mm. The system response f(v̲|x) was a 1D Gaussian pinhole whose FWHM w varied from 0.9 to 2.9 mm. For each pinhole size, we performed 4000 realizations, where the mean number of photons per realization was 100w, i.e. the sensitivity increased linearly with pinhole size (see (21)). An estimate $\hat{λ}$ (x̲) was computed for each realization using the apodized inverse filter (18), which in this case corresponds to a Gaussian filter with FWHM √β ² - w ², as shown in (22).

Figure 4 shows the sample means of the 4000 realizations for each of the 21 pinhole sizes considered, ranging from 0.9 to 2.9 mm FWHM in 0.1 mm increments. The 21 curves are indistinguishable because we are fixing the reconstructed spatial resolution to β = 3mm FWHM, as confirmed by Figure 4. We also computed the sample standard deviations from the 4000 realizations for each pinhole size. Three of the 21 curves are shown in Fig. 5. Note that the variance-minimizing pinhole size is 3/√2 ≈ 2.1mm FWHM, which has the lowest empirical variance of the three curves shown. (The plot would be difficult to interpret if all 21 curves were shown.) To verify that the theoretically predicted variance-minimizing pinhole size indeed yields the lowest variance, Fig. 6 shows the relative empirical standard deviations for each of the 119 spatial positions x for which λ(x̲) > 0 as a function of pinhole size w. All of the curves have a minimum near the predicted value of 2.1 mm FWHM.

Figure 4. Sample means of 4000 realizations of the estimates of λ(x), for 21 Gaussian pinhole sizes ranging from 0.9 to 2.9 mm FWHM. The 21 mean curves are virtually indistinguishable since the reconstructed spatial resolution has been held fixed at 3mm FWHM.
Download Full Size | PDF

Figure 5. Sample standard deviations of 4000 realizations of the estimates of λ(x), for 3 of the 21 Gaussian pinhole sizes: 1, 2.1, and 2.9 mm FWHM.
Download Full Size | PDF

Figure 6. Normalized sample standard deviations of 4000 realizations of the estimates of λ(x), versus the FWHM of the Gaussian pinhole sizes. There are 119 plots, one for each of the x positions for which λ(x) > 0. The minimum is consistently near the theoretical prediction of 3/√2 ≈ 2.1.
Download Full Size | PDF

7. Discussion

We have analyzed the performance of a kernel-based indirect density estimation method for image recovery from list-mode measurements. We showed, under a few simplifying assumptions, that the variance-minimizing pinhole width is proportional to the desired reconstructed spatial resolution. The simplifying assumptions include consideration of a shift-invariant imaging system, a spatially smooth emitting object, and a particular kernel based on an apodized inverse filter. Empirical results confirmed that the predicted variance-minimizing pinhole size yielded the lowest variability estimates, even for an object that was far from “spatially smooth.” We conjecture that there should be a monotonic relationship between desired reconstructed spatial resolution and variance-minimizing pinhole width even for broader classes of image recovery methods and more general imaging systems. Exploring this conjecture will be the subject of future work.

Although we have focused on pinhole size in this paper, a more general question would be ‘what is the optimal pinhole transmissivity function for a given target reconstructed spatial resolution?’. We conjecture that the density estimation approach described in this paper may be useful for exploring this question.

Acknowledgement

This work was supported in part by NIH grants CA-60711 and CA-54362 and by the Whitaker Foundation. The author gratefully acknowledges K. M. Brown and W. L. Rogers for helpful discussions.

Footnotes

¹	This intuition is somewhat consistent with the findings of Myers et al.[3] on the relationship between optimum aperture size and signal size in the context of signal detection, although low image variance need not necessarily be associated with high detection SNR since correlation properties are also important.
²	All integrals over d x̲ are d-dimensional.
³	Strictly speaking this is a conditional pdf, conditioned on the event that the emission is detected. We consider only the detected emissions, so for simplicity we omit the notation for conditioning on detection. With such notation, (1) is just a form of Bayes rule.
⁴	Without loss of generality, one can rescale the time axis exponentially to account for radioactive decay.
⁵	Silverman [18, p. 27] refers to such methods as general weight function estimates in the context of direct density estimation.
⁶	Our purpose here is to analyze such estimators for the goal of system design, not to argue the merits of such estimators over alternatives.
⁷	Equation (8) is closely related to equation (3.6) on on p. 36 of [18] for direct density estimation; the remainder of the derivation is distinct to indirect density estimation.
⁸	Neglecting edge effects at the boundaries of the field-of-view, and assuming that any magnification factor has already been accounted for in the V̲_n ’s [26].
⁹	As a further approximation, one can assume $\hat{λ}$ ≈ λ if the scale of the spatial variations in λ is large relative to the FWHM of h.
¹⁰	The transmissivity of the 1D version of this pinhole has the form of the Laplacian pdf ½e ^-\|x\|, hence the name—for lack of a better name.

References

1. B M Tsui, C E Metz, F B Atkins, S J Starr, and R N Beck, “A comparison of optimum detector spatial resolution in nuclear imaging based on statistical theory and on observer performance,” Phys. Med. Biol. 23654–676 (1978). [CrossRef] [PubMed]

2. H H Barrett, J N Aarsvold, H B Barber, E B Cargill, R D Fiete, T S Hickernell, T D Milster, K J Myers, D D Patton, R K Rowe, R H Seacat, W E Smith, and J M Woolfenden, “Applications of statistical decision theory in nuclear medicine,” In C N de Graaf and M A Viergever, editors, Proc. Tenth Intl. Conf. on Information Processing in Medical Im. , (Plenum Press, New York, 1987) pp. 151–166.

3. K J Myers, J P Rolland, H H Barrett, and R F Wagner, “Aperture optimization for emission imaging: effect of a spatially varying background,” J. Opt. Soc. Am. A 71279–1293 (1990). [CrossRef] [PubMed]

4. H H Barrett, “Objective assessment of image quality: effects of quantum noise and object variability,” J. Opt. Soc. Am. A 71266–1278 (1990). [CrossRef] [PubMed]

5. J P Rolland, H H Barrett, and G W Seeley, “Ideal versus human observer for long-tailed point spread functions: does deconvolution help?” Phys. Med. Biol. 361091–1109 (1991). [CrossRef] [PubMed]

6. H H Barrett, J Yao, J P Rolland, and K J Myers, “Model observers for assessment of image quality,” Proc. Natl. Acad. Sci. 909758–65 (1993). [CrossRef] [PubMed]

7. C K Abbey and H H Barrett, “Linear iterative reconstruction algorithms: study of observer performance,” In Information Processing in Medical Imaging, (Kluwer, Dordrect, 1995) pp 65–76.

8. H H Barrett, J L Denny, R F Wagner, and K J Myers, “Objective assessment of image quality. II. Fisher information, Fourier crosstalk, and figures of merit for task performance,” J. Opt. Soc. Am. A 12834–852 (1995). [CrossRef]

9. N E Hartsough, H H Barrett, H B Barber, and J M Woolfenden, “Intraoperative tumor detection: Relative performance of single-element, dual-element, and imaging probes with various collimators,” IEEE Trans. Med. Imaging 14259–265 (1995). [CrossRef] [PubMed]

10. T Kanungo, M Y Jaisimha, J Palmer, and R M Haralick, “A methodology for quantitative performance evaluation of detection algorithms,” IEEE Trans. Image Process. 41667–1674 (1995). [CrossRef] [PubMed]

11. E P Ficaro, J A Fessler, P D Shreve, J N Kritzman, P A Rose, and J R Corbett, “Simultaneous transmission/emission myocardial perfusion tomography: Diagnostic accuracy of attenuation-corrected 99m-Tc-Sestamibi SPECT,” Circulation 93463–473 (1996). http://www.eecs.umich.edu/?fessler [CrossRef] [PubMed]

12. A O Hero, J A Fessler, and M Usman, “Exploring estimator bias-variance tradeoffs using the uniform CR bound,” IEEE Trans. Signal Process. 442026–2041 (1996). http://www.eecs.umich.edu/?fessler [CrossRef]

13. J A Fessler and A O Hero, “Cramer-Rao lower bounds for biased image reconstruction,” In Proc. Midwest Symposium on Circuits and Systems, Vol. 1, (IEEE, New York, 1993) pp 253–256. http://www.eecs.umich.edu/~fessler

14. Mohammad Usman, “Biased and unbiased Cramer-Rao bounds: computational issues and applications,” PhD thesis, Univ. of Michigan, Ann Arbor, MI, 48109-2122, Ann Arbor, MI., 1994.

15. Chor-Yi Ng, “Preliminary studies on the feasibility of addition of vertex view to conventional brain SPECT imaging,” PhD thesis, Univ. of Michigan, Ann Arbor, MI, 48109-2122, Ann Arbor, MI., January 1997.

16. F O’Sullivan and Y Pawitan, “Bandwidth selection for indirect density estimation based on corrupted histogram data,” J. Am. Stat. Assoc. , 91(434):610–26, June (1996). [CrossRef]

17. P P B Eggermont and V N LaRiccia, “Nonlinearly smoothed EM density estimation with automated smoothing parameter selection for nonparametric deconvolution problems,” J. Am. Stat. Assoc. , 92(440):1451–1458, December (1997). [CrossRef]

18. B W Silverman, Density estimation for statistics and data analysis, (Chapman and Hall, New York, 1986).

19. I M Johnstone, “On singular value decompositions for the Radon Transform and smoothness classes of functions,” Technical Report 310, Dept. of Statistics, Stanford Univ., January 1989.

20. I M Johnstone and B W Silverman, “Discretization effects in statistical inverse problems,” Technical Report 310, Dept of Statistics, Stanford Univ., August 1990.

21. I M Johnstone and B W Silverman, “Speed of estimation in positron emission tomography,” Ann. Stat. 18251–280 (1990). [CrossRef]

22. P J Bickel and Y Ritov, “Estimating linear functionals of a PET image,” IEEE Trans. Med. Imaging 1481–87 (1995). [CrossRef] [PubMed]

23. B W Silverman, “Kernel density estimation using the fast Fourier transform,” Appl. Stat. 3193–99 (1982). [CrossRef]

24. B W Silverman, “On the estimation of a probability density function by the maximum penalized likelihood method,” Ann. Stat. 10795–810 (1982). [CrossRef]

25. D L Snyder and M I Miller, Random point processes in time and space, (Springer Verlag, New York, 1991). [CrossRef]

26. A Macovski, Medical imaging systems, (Prentice-Hall, New Jersey, 1983).

27. M C Jones, J S Marron, and S J Sheather, “A brief survey of bandwidth selection for density estimation,” J. Am. Stat. Assoc. , 91(433):401–407, March (1996). [CrossRef]

28. P P B Eggermont and V N LaRiccia, “Maximum smoothed likelihood density estimation for inverse problems,” Ann. Stat. 23199–220 (1995). [CrossRef]

29. Y-C Tai, A Chatziioannou, M Dahlbom, and E J Hoffman, “Investigation on deadtime characteristics for simultaneous emission-transmission data acquisition in PET,” In Proc. IEEE Nuc. Sci. Symp. Med. Im. Conf., (IEEE, New York, 1997).

30. D L Snyder and D G Politte, “Image reconstruction from list-mode data in emission tomography system having time-of-flight measurements,” IEEE Trans. Nucl. Sci. 201843–1849 (1983). [CrossRef]

31. H H Barrett, Timothy White, and Lucas C Parra, “List-mode likelihood,” J. Opt. Soc. Am. A 142914–2923 (1997). [CrossRef]

32. H H Barrett and W Swindell, Radiological imaging: the theory of image formation, detection, and processing, (Academic, New York, 1981).

33. V Ochoa, R Mastrippolito, Y Charon, P Laniece, L Pinot, and L Valentin, “TOHR: Prototype design and characterization of an original small animal tomograph,” In Proc. IEEE Nuc. Sci. Symp. Med. Im. Conf., (IEEE, New York, 1997).

34. S Geman and C R Hwang, “Nonparametric maximum likelihood estimation by the method of sieves,” Ann. Stat. 10401–414 (1982). [CrossRef]

35. R Bracewell, The Fourier transform and its applications, (McGraw-Hill, New York, 1978).