Spatial Resolution and Noise Tradeoffs in Pinhole Imaging System Design: A Density Estimation Approach

This paper analyzes the tradeo between spatial resolution and noise for simple pinhole imaging systems with photon-counting detectors. We consider image recovery algorithms based on density estimation methods using kernels that are based on apodized inverse filters. This approach allows a continuous-object, continuous-data treatment of the problem. The analysis shows that the pinhole size that minimizes the estimate variance for a specied reconstructed spatial resolution is directly proportional to that spatial resolution. For a Gaussian pinhole, the variance-minimizing full-width half maximum (FWHM) of the pinhole equals the desired object spatial resolution divided by p2. Simulation results confirm this conclusion empirically. The general approach is a potentially useful addition to the collection of tools available for imaging system design.


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 systemdependent 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 increase1 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.

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 1 , X 2 , . ..} drawn from a Poisson spatial point process [25].In particular, the joint probability that the first n 0 emissions originate in any (measurable) regions B n ⊆ Ω is given by2 : 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 (density3 ) that that emission originated from spatial location x is given by where r = λ(x)s(x) dx is the total rate of detected events (with units detected counts per unit time4 ).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 X n = 1 m V n .Given the recorded positions of many such photons V 1 , . . ., 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: 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: 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.
This is reasonable assumption, except possibly at high count rates when deadtime factors and pulse pile-up effects are significant, e.g.[29].

The estimation problem
Suppose the imaging system records a total of N events during a prespecified period of time t 0 .By assumption, N is a Poisson random variable with mean 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 1 , . .., V N that each have the following marginal pdf 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 n=1 .We assume that t 0 , 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 Thus a natural estimator for the emission-rate density λ(x) is simply 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).

Kernel-based indirect density estimator
In this paper, we consider the following class of kernel-based indirect density estimators5 : 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: 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 6 , 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 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 λ(x).

Mean function
The mean function for the estimator λ(x) is derived as follows 7 : Thus we have the following linear relationship between the estimator mean and the true emission density: where 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).

Second-moment functions
Before computing the autocorrelation function of λ, we first note that since N is Poisson, Then from (7), the autocorrelation function for λ(x) is derived as follows: Therefore the autocovariance function for λ is To simplify, note that so the autocovariance function is In particular, the variance function is (This equation is a space-variant generalization of (10.31) in [32].)Note that the variance depends inversely on the scan time t 0 , 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.

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 system8 .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 0 for some positive constant s 0 .Then the above expressions simplify as follows.
The mean expression (9) becomes: where x = v − x and * denotes d-dimensional convolution.Thus we have the following convolution relationship (cf (10.11) of [32]): 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 with corresponding frequency response or overall transfer function where F (u) = f(x)e −i2πu•x dx denotes the d-dimensional Fourier transform of f(x).
Similarly, the inner variance term in (12) becomes: 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 Therefore in the shift-invariant case it is straightforward to compute variances (approximately) using FFTs to calculate the convolutions.

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.
where we define λ = h * λ (cf (10.41) of [32]).For small β or for spatially smooth objects this approximation should be fairly accurate 9 .For pinhole imaging, g β is a real function, so by combining the above approximation with Parseval's theorem: This is a very tractable approximation to the estimator variance.

Resolution-noise tradeoffs
In general, both the sensitivity s 0 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).

Apodized inverse filter
Consider a general pinhole with transmissivity function t(x) ≥ 0, which we assume is normalized so that t(x) dx = 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(x) = 1 w d t(x/w) for which H(u) = T (wu).Define the apodized inverse filter 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 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 0 = c 0 w p for some constant c 0 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 2 .From (17) the variance is approximately: where z = wu and c 1 = λ(x)/t 0 c 0 .To find the pinhole width w that minimizes the variance, we zero the partial derivative of the variance σ 2 with respect to the width w: or, by defining α = β/w: The above equality depends only on the ratio α = β/w.So if there is a root α 0 > 0 that corresponds to a global minimizer of the variance σ 2 , then the variance-minimizing w is proportional to the reconstructed spatial resolution β through the relationship w min = α −1 0 β.

Relationship to sieves
The above apodized inverse filter is closely related to the method of sieves for density estimation [34], in the sense that λ(x) is an unbiased estimate of β −d a(x/β) * λ(x).Profile through an approximate Gaussian pinhole.
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 However, if µl is sufficiently large, then we can approximate this transmissivity by where w is the FWHM of the pinhole response (i.e.τ w (w/2) = 1/2) and The sensitivity of this pinhole is therefore which is proportional to w d as expected.The normalized transmissivity (for unit pinhole width w = 1) is with corresponding frequency response We choose a Gaussian apodizing function A(u) = e −π(ρ/κ) 2 so that the PSF corresponding to A(βρ) has FWHM β.The corresponding recovery filter is thus with corresponding space-domain recovery kernel Substituting A(•) into (19), the variance function is approximately  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: 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.Profile through an approximate Laplacian pinhole.
A somewhat more conventional pinhole is the Laplacian10 pinhole shown in Fig. 3.The exact transmissivity of such a pinhole is τ w (r) = e −µlr/rb , r ≤ r b e −µl , r≥ r b .
If µl is sufficiently large, we can approximate this transmissivity by τ w (r) = e −γr/w , where γ = 2 log 2 and w is the FWHM of the pinhole response.The sensitivity of this pinhole is which is also proportional to w 2 .The normalized transmissivity is We again choose a Gaussian apodizing function A(ρ) = e −π(ρ/κ) 2 where κ was defined in (20), so the PSF again has FWHM β.Substituting into (19), the estimate variance is approximately After some tedious integration, we arrive at 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.

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 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 λ(x) was computed for each realization using the apodized inverse filter (18), which in this case corresponds to a Gaussian filter with FWHM β 2 − w 2 , 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.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.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.

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 varianceminimizing 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.
Figure 1.Profile through an approximate Gaussian pinhole.

Figure 2 .
Figure 2.Standard deviation of estimate versus Gaussian pinhole width.
Figure 3.Profile through an approximate Laplacian pinhole.
which has corresponding frequency response[35]

Figure 4 .
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.