- Email: [email protected]
Determining detector requirements for medical imaging applications1
Nuclear Instruments and Methods in Physics Research A 409 (1998) 501—507
Determining detector requirements for medical imaging applications1 Neal H. Clinthorne!,*, Chor-yi Ng", Joerg Strobel#,2, Chia-ho Hua", James W. LeBlanc$, Scott J. Wilderman$, W. Leslie Rogers! ! Division of Nuclear Medicine, The University of Michigan Medical Center, 3480 Kresge III Box 0552, 204 Zina Pitcher Place, Ann Arbor, MI 48109-0552, USA " Department of Biomedical Engineering, The University of Michigan, Ann Arbor, MI 48109-0552, USA # Technical University of Munich, Munich, Germany $ Department of Nuclear Engineering and Radiological Sciences, The University of Michigan, Ann Arbor, MI 48109-0552, USA
Abstract The Crame´r—Rao lower bound on estimator variance is used as a tool for assessing radiation detector requirements for nuclear medicine imaging. The classical bound for unbiased estimators is reviewed and its relationship to the well-known propagation-of-error formula is demonstrated. The uniform Crame´r—Rao bound is described for situations where unbiased estimation is impractical such as in tomographic imaging. The bounds are used to evaluate the necessary detector requirements for a CsI(Tl)/Si photodiode scintillation camera, and to compare the performance of a Comptonscatter camera to that of a parallel-hole collimator/Anger camera system. ( 1998 Published by Elsevier Science B.V. All rights reserved. Keywords: Noise; Noise-propagation; Crame´r—Rao lower bound
1. Introduction Determining the necessary detector requirements for imaging systems that directly measure a source distribution (e.g., a CCD camera) is both straightforward and intuitive. Here, the concepts of detector spatial and energy resolution have a direct * Corresponding author. Tel.: #1 734 764 4289; fax: #1 734 764 0288; e-mail: [email protected]. 1 This work was supported by the US DHHS under NIH Grants R01 CA32846 and R01 CA65637. 2 Presently at: The University of Michigan, Ann Arbor, MI 48109-0552, USA.
relationship to imaging performance. Increasingly, however, in medicine and as well as in other fields, detectors are only one component of a complex imaging system that makes indirect measurements of the source. An obvious example is tomographic imaging where an image is reconstructed from projections of the object. And while detector properties remain important, it is not so clear how they influence overall system performance because the source distribution must be decoded or reconstructed from these indirect measurements. Perhaps the question of greatest significance is: given a desired system performance what are the detector characteristics necessary to achieve it?
0168-9002/97/$19.00 ( 1998 Published by Elsevier Science B.V. All rights reserved PII S 0 1 6 8 - 9 0 0 2 ( 9 7 ) 0 1 3 0 3 - X
VI. MEDICAL APPLICATIONS
502
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
A number of methods exist for analyzing requirements; unfortunately, most include the reconstruction (or more generally the estimator) as an integral element in the analysis. While this must of course be done eventually, it is undesirable from the viewpoint of the initial system design. Many well-known estimators (e.g., filtered backprojection reconstruction) are approximations that may have some computational advantage but often have significant disadvantages such as poor noise handling properties. Comparisons among imaging systems, some of which may be poorly chosen estimators, may not lead to valid conclusions. In contrast, the Crame´r—Rao (CR) lower bound on estimator variance is independent of the estimator, can quantify the intrinsic quality of the measurements, and can be used to determine optimum system design parameters. In the next sections, various forms of the CR bound are reviewed and the relationship with the more standard propagation-of-error method is demonstrated. Subsequently, the bounds are applied to analyze the performance of two imaging systems currently of interest in nuclear medicine.
the Fisher information matrix F [1] $%& E[+2h log f (yD h)]~1 KhK 5F~1 "
(1)
and lower bounds the estimator covariance matrix KhK in the sense that KhK !F~1 is a non-negativedefinite matrix. Importantly, the limiting covariance specified by the CR bound can be asymptotically achieved using the maximum-likelihood (ML) estimator. For conditionally Poisson distributed measurements as well as for Gaussian measurements in which the covariance is not a function of the parameters, the Fisher information matrix assumes the following form: F " (+yM (h))TK~1 y +yN (h),
(2)
where K is the covariance matrix of the measurey ments and yN (h) is a vector representing the mean measurements as a function of the parameters. For Poisson measurements, the data covariance matrix Ky " diagMyN (h)N. 3. Relation to the propagation of error formula
2. The Crame´ r—Rao lower bound There are two essential components for the bound calculations: (1) a finite parameterization of the source, which we define as a vector of M para$%& [h , ,h ]T; and (2) a statistical meters h " 12 M model, represented by the conditional probability density or mass function f (y D h), which relates these parameters to a vector of N measurements $%& [y ,2,y ]T. The source parameterization is y" 1 N somewhat arbitrary (a pixellated object, for example) although it should be meaningful in regards to the desired imaging tasks (one would not want pixels that are too large). In contrast, the statistical model must be accurate to realistically evaluate performance. It is convenient at this point to define an estimator of h as hK (y) where the dependence on the measurements is explicit. With the above definitions, the CR bound for unbiased estimators (those whose bias $%& E(hK )!h " 0) is defined as the inverse of bhK (h) "
Although CR bounds are gaining popularity in design problems, many investigators are more familiar with standard techniques such as the “error propagation formula” [2] (Ch. 4). This common expression can be derived by expanding the estimator hK (y) in a Taylor’s series around yN (h) and keeping the constant and linear terms: hK (y) + hK (yM )#+yhK (yN )(y!yN ). The approximate variance is easily calculated using the linearized estimator (3) K 4 + +yhK (yN )Ky(+yhK (yN ))T, h which is the well-known error propagation formula extended to a vector of parameters. The limitation of this technique is that it requires an estimator to be specified in advance, which as noted above is not always the most desirable route. To see the relation to the CR bound, consider the ML estimator, which is defined implicitly as hK " arg maxh f (y D h).
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
Using techniques similar to those detailed in Ref. [3] it is straightforward to show that for Poisson distributed measurements +y hK (yN ) " [(+yN (h))TK~1 y +yN (h)]~1(+yN (h))TK~1 y for the ML estimator. Consequently, +yhK (yN )Ky(+yhK (yN ))T , [(+yN (h))TK~1 y +yN (h)]~1.
(4)
In other words, the propagation-of-error formula and the CR bound are equivalent for the ML estimator. This relationship holds not only for Poisson measurements but also for Gaussiandistributed measurements in which the covariance is not a function of the parameter vector h.
503
pair of equations: d (h,j) " E(I#jF)~1e E, (5) p p (6) B (h,j) " j2eT(I#jF)~1F(I#jF)~1e , p p p where B (h,j) is the limiting variance for estimating p h , d (h,j) is the corresponding bias-gradient length, p p I is the identity matrix, and e is an M]1 vector of p zeros except for a one in the pth position. The bound and the corresponding bias-gradient length are calculated parametrically by sweeping j through the interval [0,R). Like the CR bound for unbiased estimators, the UCR bound is useful for system design because it is asymptotically achievable using the appropriately regularized ML estimator.
4. The uniform Crame´ r—Rao bound A pitfall to use of the classical bound in imaging system design is that unbiased estimation is often not practical. This is typically the case in image reconstruction where one must often accept some resolution loss — equivalent to adding bias — in order to sufficiently reduce the mean-squared error in reconstructed images.3 While a form of the CR bound for biased estimation exists, it is impractical for general design since it bounds performance only for those estimators constrained to having a given gradient of the bias with respect to the parameters (or bias-gradient). The uniform CR (UCR) bound removes this restriction by specifying the limiting variance of any estimator as a function of the length of the biasgradient [4] (as opposed to the entire bias-gradient vector for each parameter). As an index, bias-gradient length is similar to the concept of point-spread function width in the reconstructed image. It can — with caution — be interpreted as a measure of spatial resolution. The UCR bound divides the variance/bias-gradient length plane (or sigma—delta plane) into achievable and unachievable regions as shown in Fig. 1. For the pth element of the parameter vector h, the uniform bound is described by the following
3 The conditional mean-squared error of an estimator is equal to the sum of its squared bias and variance.
5. Applications To illustrate possible uses of the bounds described above, two applications are presented. The first evaluates the spatial resolution of a CsI(Tl)/ silicon photodiode camera as a function of detector noise. Because the estimation problem is relatively well-posed, the vector form of the Crame´r—Rao bound for unbiased estimators is used. The second application compares the performance of a conventional mechanically collimated imaging system with that of an electronically collimated Comptonscatter aperture. Because unbiased estimation is generally not practical in this setting, the uniform CR bound is applied. 5.1. CsI(¹l)/Silicon Photodiode Camera With the advent of low-noise silicon-pad detectors [5], there has been significant interest in small, solid-state scintillation cameras for nuclear medicine. The feasibility of a camera based on the Anger principle using silicon photodiodes is straightforward to ascertain using the CR bound for unbiased estimators. Referring to Fig. 2, the camera is constructed by mating a silicon photodetector array with a CsI(Tl) scintillator 5 mm thick. This thickness offers a good compromise between detection efficiency at 141 keV (99mTc) and a narrow light spread function (LSF),
VI. MEDICAL APPLICATIONS
504
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
Fig. 1. UCR bound curves separate the sigma—delta performance plane into a region that can potentially be achieved by the appropriate image reconstruction, and an unachievable region. Shown plotted are curves for two imaging systems. The lower curve corresponds to a system that gives uniformly better performance.
Fig. 2. Diagram of a solid-state scintillation camera. A siliconpad photodetector is optically coupled to a 5 mm thick CsI(Tl) scintillation crystal. From the relative photodiode outputs the position and energy of the interacting c-ray can be determined.
which is important for good spatial resolution. The LSF and its gradient at the photodetectors were estimated using optical Monte Carlo methods [6] for two scintillator surface treatments. The first was a diffuse surface, which resulted in a LSF having
a width of 6.7 mm FWHM at the photodetectors. The second surface treatment employed retroreflectors (similar to cube-corners) to narrow the distribution of reflected light. The width of the second distribution was 3.5 mm FWHM. We assumed that a 141 keV photon resulted in the generation of a mean of 1400 electron—hole pairs in the pad detector (4000 intial scintillation photons, half were absorbed by the blackened crystal edges, 70% quantum efficiency). The detector consisted of a 16 ] 16 array of 1.875]1.875 mm2 pads. The statistical model was straightforward. The measurement vector y was generated by integrating the LSF over the area of each photodetector pad. Measurements were assumed to have the following distribution y & PoissonMyN (h)#nN, where h is a vector consisting of the x, y, and z components of position as well as the photon energy, yN (h) is the integrated LSF (in units of e~), and n is the equivalent noise charge for each
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
505
Fig. 3. FWHM spatial resolution in the x and y directions as a function of photodetector noise (ENC). Results are shown for two LSF widths: 3.5 mm FWHM (retroreflectors) and 6.7 mm FWHM (ground surface).
photodiode.4 The corresponding bound is KhK 5[(+yN (h))T(diag(yN (h)#n))~1+yN (h)]~1.
(7)
Once the LSFs and gradients were estimated, the bound calculation simply required inversion of the 4]4 Fisher matrix. Limiting resolutions in the x and y directions are plotted in Fig. 3. The FWHM resolution was estimated as 2.35] the limiting standard deviation obtained from the appropriate diagonal elements of F~1. Note that the narrow light distribution (retroreflector top) gives significantly better resolution than the more conventional ground surface treatment. Requirements for the silicon detector can be easily determined from these plots. For example, to operate at 1 mm FWHM resolution requires a detector noise of approximately 20 e~ ENC for the scintillator having a ground surface, 40 e~ ENC for the scintillator with retroreflectors. The next step is 4 While the photodiode noise contains components that can be considered as both Gaussian- and Poisson-distributed, calculations are simplified greatly by assuming the diode noise is entirely Poisson.
of course to corroborate these predictions with experimental results. 5.2. Compton vs. mechanical collimation The Compton-scatter aperture has been of interest in nuclear medicine for nearly two decades [7]. Referring to Fig. 4, photons from the source undergo Compton-scattering in a first detector and are detected in time-coincidence by a second. Knowledge of the interaction locations in the two detectors as well as the recoil electron and scattered photon energies allows the emission location to be reconstructed to within a conical ambiguity. The actual emitter distribution can be reconstructed from a large number of these measurements. Since Compton-scatter apertures can have a large active solid-angle, a counting efficiency two orders of magnitude higher than conventional, mechanically collimated systems is possible; in addition, they are not subject to the inherent resolution-sensitivity tradeoff that constrains mechanical collimation. Despite the promise of Compton cameras, the actual improvement in performance over that of more
VI. MEDICAL APPLICATIONS
506
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
Fig. 4. Diagram of the ring Compton-scatter camera. Photons originating in the source Compton-scatter in the first detector and are detected in time-coincidence in a second. Origin of the detected event can be localized to the surface of a cone given the detected energies and positions.
conventional systems has been unclear because the coded Compton camera data need to be reconstructed to form an image whereas a conventionally acquired image does not.5 There are many such reconstruction methods and a poor choice would unfairly penalize the Compton camera. Again, the CR bound can be used to evaluate performance irrespective of the estimator (reconstruction algorithm). In this case the dimension of the reconstructed images relative to the number of measurements results in an inversion problem that is “practically ill-posed”; i.e., unbiased estimation will result in unusable images due to large variance. The uniform CR bound is employed to compare performance of the two systems as a function of bias-gradient length or reconstructed resolution. Models for the two systems (yM (h)) are described in detail in Ref. [8] and are only briefly summarized
here. The comparison was performed for an incident energy of 141 keV, which corresponds to the emission line of the most important radionuclide in nuclear medicine, 99mTc. The source model was taken to be a plane 30]30 cm2, located 10 cm from both the conventional and Compton camera. It was quantized into a 64]64 element array of pixels. The conventional system consisted of a parallel hole collimator and an Anger camera having a combined resolution at 10 cm of 5.8 mm FWHM and an efficiency of 1.8]10~4. The Compton camera was comprised of a cylindrical second detector with radius 25 cm and length 38 cm along with a position-sensitive silicon first detector of sufficient size to achieve a counting efficiency 25] higher than that of the conventional system. Energy resolution of the silicon detector, which is critical to performance, was taken as 1 keV FWHM. Because of the size of the problem (F is 4096]4096 for this case), a conjugate-gradient algorithm was employed to solve the following problem for the 4096]1 vector d (I#jF)d " e . p Using d, the limiting performance for a single pixel in the object, (d (h,j),B (h,j)) was calculated for p p several values of j. Fig. 5 plots performance in terms of the limiting standard deviation in estimating the intensity of the on-axis pixel in the image. Because of the significant multiplexing in the projection domain, Compton camera performance depends strongly on object size, and this is reflected in the bound. Results have been plotted for 5, 10 and 30 cm diameter uniformly emitting disk sources. Performance of the conventional system is only weakly dependent on source size. Results for the 10 cm diameter object are shown. Note that the Compton camera can provide significantly improved performance over parallel-hole collimation. However, because the projections must be decoded, gains are less than the 5] improvement one might predict based solely on increased counting efficiency. Of course, the bound quantifies this decoding penalty.
6. Discussion and conclusion 5 Although the conventional system does not require explicit reconstruction, some processing — either for resolution-recovery or noise reduction — is often desirable.
While various forms of the CR bound can are useful for assessing detector requirements, there are
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
507
Fig. 5. Uniform CR bound plots comparing the performance of the ring Compton camera shown in Fig. 4 (solid lines) and parallel-hole collimator/Anger camera (dashed line).
several caveats. First, bound predictions are only as good as the statistical model. If the model does not accurately reflect reality, inferences made from the calculations may not be useful. This is no different than for any design tool — at some point predictions must be validated with measurements. Second, when comparing the performance of two or more designs, the same parameterization of the source must be used. It hardly needs to be mentioned that this parameterization must also be meaningful in regards to the imaging tasks of interest. The judgement of the system designer must guide these modelling decisions; the CR bound cannot itself be used to evaluate the quality of either the statistical model or the parameterization. Finally, there are situations in which the CR bounds provide an overly optimistic estimate of performance. In particular, we have found this true for Poisson measurements when the photon rate is very low or the ML estimator is highly non-linear in the data. These tend to be situations where the propagationof-error method also has poor accuracy. Often alternative lower bounds can be developed [9] when this occurs.
Mindful of the above cautions, CR bounds are invaluable tools for the design of imaging systems and for the assessment of necessary detector requirements.
References [1] R. Blahut, Principles and Practice of Information Theory, Addison-Wesley, Reading, MA, 1987. [2] G.F. Knoll, Radiation Detection and Measurement, Wiley, New York, 1979. [3] J.A. Fessler, IEEE Trans. Image. Process. 5 (3) (1996) 493. [4] A.O. Hero, J.A. Fessler, M. Usman, IEEE Trans. Signal. Process. 44 (8) (1996) 2026. [5] P. Weilhammer, E. Nygard, W. Dulinski, A. Czermak, F. Djama, S. Gadomski, S. Roe, A. Rudge, F. Schopper, J. Strobel, Nucl. Instr. and Meth. A 383 (1996) 89. [6] G.F. Knoll, T.F. Knoll, T.M. Henderson, IEEE Trans. Nucl. Sci. NS-35 (1988) 872. [7] M. Singh, R. Leahy, R. Brechner, T. Hebert, IEEE Trans. Nucl. Sci. NS-35 (1) (1988) 772. [8] N.H. Clinthorne, C.-Y. Ng, C. Hua, J.W. LeBlanc, S.J. Wilderman, W.L. Rogers, Conf. Record of the 1996 IEEE Nuclear Science Symp. and Medical Imaging Conf. 3 (1997) pp. 788—792. [9] N.H. Clinthorne, A.O. Hero, N. Petrick, W.L. Rogers, Nucl. Instr. and Meth. A 299 (1990) 157.
VI. MEDICAL APPLICATIONS
Determining detector requirements for medical imaging applications1 Neal H. Clinthorne!,*, Chor-yi Ng", Joerg Strobel#,2, Chia-ho Hua", James W. LeBlanc$, Scott J. Wilderman$, W. Leslie Rogers! ! Division of Nuclear Medicine, The University of Michigan Medical Center, 3480 Kresge III Box 0552, 204 Zina Pitcher Place, Ann Arbor, MI 48109-0552, USA " Department of Biomedical Engineering, The University of Michigan, Ann Arbor, MI 48109-0552, USA # Technical University of Munich, Munich, Germany $ Department of Nuclear Engineering and Radiological Sciences, The University of Michigan, Ann Arbor, MI 48109-0552, USA
Abstract The Crame´r—Rao lower bound on estimator variance is used as a tool for assessing radiation detector requirements for nuclear medicine imaging. The classical bound for unbiased estimators is reviewed and its relationship to the well-known propagation-of-error formula is demonstrated. The uniform Crame´r—Rao bound is described for situations where unbiased estimation is impractical such as in tomographic imaging. The bounds are used to evaluate the necessary detector requirements for a CsI(Tl)/Si photodiode scintillation camera, and to compare the performance of a Comptonscatter camera to that of a parallel-hole collimator/Anger camera system. ( 1998 Published by Elsevier Science B.V. All rights reserved. Keywords: Noise; Noise-propagation; Crame´r—Rao lower bound
1. Introduction Determining the necessary detector requirements for imaging systems that directly measure a source distribution (e.g., a CCD camera) is both straightforward and intuitive. Here, the concepts of detector spatial and energy resolution have a direct * Corresponding author. Tel.: #1 734 764 4289; fax: #1 734 764 0288; e-mail: [email protected]. 1 This work was supported by the US DHHS under NIH Grants R01 CA32846 and R01 CA65637. 2 Presently at: The University of Michigan, Ann Arbor, MI 48109-0552, USA.
relationship to imaging performance. Increasingly, however, in medicine and as well as in other fields, detectors are only one component of a complex imaging system that makes indirect measurements of the source. An obvious example is tomographic imaging where an image is reconstructed from projections of the object. And while detector properties remain important, it is not so clear how they influence overall system performance because the source distribution must be decoded or reconstructed from these indirect measurements. Perhaps the question of greatest significance is: given a desired system performance what are the detector characteristics necessary to achieve it?
0168-9002/97/$19.00 ( 1998 Published by Elsevier Science B.V. All rights reserved PII S 0 1 6 8 - 9 0 0 2 ( 9 7 ) 0 1 3 0 3 - X
VI. MEDICAL APPLICATIONS
502
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
A number of methods exist for analyzing requirements; unfortunately, most include the reconstruction (or more generally the estimator) as an integral element in the analysis. While this must of course be done eventually, it is undesirable from the viewpoint of the initial system design. Many well-known estimators (e.g., filtered backprojection reconstruction) are approximations that may have some computational advantage but often have significant disadvantages such as poor noise handling properties. Comparisons among imaging systems, some of which may be poorly chosen estimators, may not lead to valid conclusions. In contrast, the Crame´r—Rao (CR) lower bound on estimator variance is independent of the estimator, can quantify the intrinsic quality of the measurements, and can be used to determine optimum system design parameters. In the next sections, various forms of the CR bound are reviewed and the relationship with the more standard propagation-of-error method is demonstrated. Subsequently, the bounds are applied to analyze the performance of two imaging systems currently of interest in nuclear medicine.
the Fisher information matrix F [1] $%& E[+2h log f (yD h)]~1 KhK 5F~1 "
(1)
and lower bounds the estimator covariance matrix KhK in the sense that KhK !F~1 is a non-negativedefinite matrix. Importantly, the limiting covariance specified by the CR bound can be asymptotically achieved using the maximum-likelihood (ML) estimator. For conditionally Poisson distributed measurements as well as for Gaussian measurements in which the covariance is not a function of the parameters, the Fisher information matrix assumes the following form: F " (+yM (h))TK~1 y +yN (h),
(2)
where K is the covariance matrix of the measurey ments and yN (h) is a vector representing the mean measurements as a function of the parameters. For Poisson measurements, the data covariance matrix Ky " diagMyN (h)N. 3. Relation to the propagation of error formula
2. The Crame´ r—Rao lower bound There are two essential components for the bound calculations: (1) a finite parameterization of the source, which we define as a vector of M para$%& [h , ,h ]T; and (2) a statistical meters h " 12 M model, represented by the conditional probability density or mass function f (y D h), which relates these parameters to a vector of N measurements $%& [y ,2,y ]T. The source parameterization is y" 1 N somewhat arbitrary (a pixellated object, for example) although it should be meaningful in regards to the desired imaging tasks (one would not want pixels that are too large). In contrast, the statistical model must be accurate to realistically evaluate performance. It is convenient at this point to define an estimator of h as hK (y) where the dependence on the measurements is explicit. With the above definitions, the CR bound for unbiased estimators (those whose bias $%& E(hK )!h " 0) is defined as the inverse of bhK (h) "
Although CR bounds are gaining popularity in design problems, many investigators are more familiar with standard techniques such as the “error propagation formula” [2] (Ch. 4). This common expression can be derived by expanding the estimator hK (y) in a Taylor’s series around yN (h) and keeping the constant and linear terms: hK (y) + hK (yM )#+yhK (yN )(y!yN ). The approximate variance is easily calculated using the linearized estimator (3) K 4 + +yhK (yN )Ky(+yhK (yN ))T, h which is the well-known error propagation formula extended to a vector of parameters. The limitation of this technique is that it requires an estimator to be specified in advance, which as noted above is not always the most desirable route. To see the relation to the CR bound, consider the ML estimator, which is defined implicitly as hK " arg maxh f (y D h).
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
Using techniques similar to those detailed in Ref. [3] it is straightforward to show that for Poisson distributed measurements +y hK (yN ) " [(+yN (h))TK~1 y +yN (h)]~1(+yN (h))TK~1 y for the ML estimator. Consequently, +yhK (yN )Ky(+yhK (yN ))T , [(+yN (h))TK~1 y +yN (h)]~1.
(4)
In other words, the propagation-of-error formula and the CR bound are equivalent for the ML estimator. This relationship holds not only for Poisson measurements but also for Gaussiandistributed measurements in which the covariance is not a function of the parameter vector h.
503
pair of equations: d (h,j) " E(I#jF)~1e E, (5) p p (6) B (h,j) " j2eT(I#jF)~1F(I#jF)~1e , p p p where B (h,j) is the limiting variance for estimating p h , d (h,j) is the corresponding bias-gradient length, p p I is the identity matrix, and e is an M]1 vector of p zeros except for a one in the pth position. The bound and the corresponding bias-gradient length are calculated parametrically by sweeping j through the interval [0,R). Like the CR bound for unbiased estimators, the UCR bound is useful for system design because it is asymptotically achievable using the appropriately regularized ML estimator.
4. The uniform Crame´ r—Rao bound A pitfall to use of the classical bound in imaging system design is that unbiased estimation is often not practical. This is typically the case in image reconstruction where one must often accept some resolution loss — equivalent to adding bias — in order to sufficiently reduce the mean-squared error in reconstructed images.3 While a form of the CR bound for biased estimation exists, it is impractical for general design since it bounds performance only for those estimators constrained to having a given gradient of the bias with respect to the parameters (or bias-gradient). The uniform CR (UCR) bound removes this restriction by specifying the limiting variance of any estimator as a function of the length of the biasgradient [4] (as opposed to the entire bias-gradient vector for each parameter). As an index, bias-gradient length is similar to the concept of point-spread function width in the reconstructed image. It can — with caution — be interpreted as a measure of spatial resolution. The UCR bound divides the variance/bias-gradient length plane (or sigma—delta plane) into achievable and unachievable regions as shown in Fig. 1. For the pth element of the parameter vector h, the uniform bound is described by the following
3 The conditional mean-squared error of an estimator is equal to the sum of its squared bias and variance.
5. Applications To illustrate possible uses of the bounds described above, two applications are presented. The first evaluates the spatial resolution of a CsI(Tl)/ silicon photodiode camera as a function of detector noise. Because the estimation problem is relatively well-posed, the vector form of the Crame´r—Rao bound for unbiased estimators is used. The second application compares the performance of a conventional mechanically collimated imaging system with that of an electronically collimated Comptonscatter aperture. Because unbiased estimation is generally not practical in this setting, the uniform CR bound is applied. 5.1. CsI(¹l)/Silicon Photodiode Camera With the advent of low-noise silicon-pad detectors [5], there has been significant interest in small, solid-state scintillation cameras for nuclear medicine. The feasibility of a camera based on the Anger principle using silicon photodiodes is straightforward to ascertain using the CR bound for unbiased estimators. Referring to Fig. 2, the camera is constructed by mating a silicon photodetector array with a CsI(Tl) scintillator 5 mm thick. This thickness offers a good compromise between detection efficiency at 141 keV (99mTc) and a narrow light spread function (LSF),
VI. MEDICAL APPLICATIONS
504
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
Fig. 1. UCR bound curves separate the sigma—delta performance plane into a region that can potentially be achieved by the appropriate image reconstruction, and an unachievable region. Shown plotted are curves for two imaging systems. The lower curve corresponds to a system that gives uniformly better performance.
Fig. 2. Diagram of a solid-state scintillation camera. A siliconpad photodetector is optically coupled to a 5 mm thick CsI(Tl) scintillation crystal. From the relative photodiode outputs the position and energy of the interacting c-ray can be determined.
which is important for good spatial resolution. The LSF and its gradient at the photodetectors were estimated using optical Monte Carlo methods [6] for two scintillator surface treatments. The first was a diffuse surface, which resulted in a LSF having
a width of 6.7 mm FWHM at the photodetectors. The second surface treatment employed retroreflectors (similar to cube-corners) to narrow the distribution of reflected light. The width of the second distribution was 3.5 mm FWHM. We assumed that a 141 keV photon resulted in the generation of a mean of 1400 electron—hole pairs in the pad detector (4000 intial scintillation photons, half were absorbed by the blackened crystal edges, 70% quantum efficiency). The detector consisted of a 16 ] 16 array of 1.875]1.875 mm2 pads. The statistical model was straightforward. The measurement vector y was generated by integrating the LSF over the area of each photodetector pad. Measurements were assumed to have the following distribution y & PoissonMyN (h)#nN, where h is a vector consisting of the x, y, and z components of position as well as the photon energy, yN (h) is the integrated LSF (in units of e~), and n is the equivalent noise charge for each
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
505
Fig. 3. FWHM spatial resolution in the x and y directions as a function of photodetector noise (ENC). Results are shown for two LSF widths: 3.5 mm FWHM (retroreflectors) and 6.7 mm FWHM (ground surface).
photodiode.4 The corresponding bound is KhK 5[(+yN (h))T(diag(yN (h)#n))~1+yN (h)]~1.
(7)
Once the LSFs and gradients were estimated, the bound calculation simply required inversion of the 4]4 Fisher matrix. Limiting resolutions in the x and y directions are plotted in Fig. 3. The FWHM resolution was estimated as 2.35] the limiting standard deviation obtained from the appropriate diagonal elements of F~1. Note that the narrow light distribution (retroreflector top) gives significantly better resolution than the more conventional ground surface treatment. Requirements for the silicon detector can be easily determined from these plots. For example, to operate at 1 mm FWHM resolution requires a detector noise of approximately 20 e~ ENC for the scintillator having a ground surface, 40 e~ ENC for the scintillator with retroreflectors. The next step is 4 While the photodiode noise contains components that can be considered as both Gaussian- and Poisson-distributed, calculations are simplified greatly by assuming the diode noise is entirely Poisson.
of course to corroborate these predictions with experimental results. 5.2. Compton vs. mechanical collimation The Compton-scatter aperture has been of interest in nuclear medicine for nearly two decades [7]. Referring to Fig. 4, photons from the source undergo Compton-scattering in a first detector and are detected in time-coincidence by a second. Knowledge of the interaction locations in the two detectors as well as the recoil electron and scattered photon energies allows the emission location to be reconstructed to within a conical ambiguity. The actual emitter distribution can be reconstructed from a large number of these measurements. Since Compton-scatter apertures can have a large active solid-angle, a counting efficiency two orders of magnitude higher than conventional, mechanically collimated systems is possible; in addition, they are not subject to the inherent resolution-sensitivity tradeoff that constrains mechanical collimation. Despite the promise of Compton cameras, the actual improvement in performance over that of more
VI. MEDICAL APPLICATIONS
506
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
Fig. 4. Diagram of the ring Compton-scatter camera. Photons originating in the source Compton-scatter in the first detector and are detected in time-coincidence in a second. Origin of the detected event can be localized to the surface of a cone given the detected energies and positions.
conventional systems has been unclear because the coded Compton camera data need to be reconstructed to form an image whereas a conventionally acquired image does not.5 There are many such reconstruction methods and a poor choice would unfairly penalize the Compton camera. Again, the CR bound can be used to evaluate performance irrespective of the estimator (reconstruction algorithm). In this case the dimension of the reconstructed images relative to the number of measurements results in an inversion problem that is “practically ill-posed”; i.e., unbiased estimation will result in unusable images due to large variance. The uniform CR bound is employed to compare performance of the two systems as a function of bias-gradient length or reconstructed resolution. Models for the two systems (yM (h)) are described in detail in Ref. [8] and are only briefly summarized
here. The comparison was performed for an incident energy of 141 keV, which corresponds to the emission line of the most important radionuclide in nuclear medicine, 99mTc. The source model was taken to be a plane 30]30 cm2, located 10 cm from both the conventional and Compton camera. It was quantized into a 64]64 element array of pixels. The conventional system consisted of a parallel hole collimator and an Anger camera having a combined resolution at 10 cm of 5.8 mm FWHM and an efficiency of 1.8]10~4. The Compton camera was comprised of a cylindrical second detector with radius 25 cm and length 38 cm along with a position-sensitive silicon first detector of sufficient size to achieve a counting efficiency 25] higher than that of the conventional system. Energy resolution of the silicon detector, which is critical to performance, was taken as 1 keV FWHM. Because of the size of the problem (F is 4096]4096 for this case), a conjugate-gradient algorithm was employed to solve the following problem for the 4096]1 vector d (I#jF)d " e . p Using d, the limiting performance for a single pixel in the object, (d (h,j),B (h,j)) was calculated for p p several values of j. Fig. 5 plots performance in terms of the limiting standard deviation in estimating the intensity of the on-axis pixel in the image. Because of the significant multiplexing in the projection domain, Compton camera performance depends strongly on object size, and this is reflected in the bound. Results have been plotted for 5, 10 and 30 cm diameter uniformly emitting disk sources. Performance of the conventional system is only weakly dependent on source size. Results for the 10 cm diameter object are shown. Note that the Compton camera can provide significantly improved performance over parallel-hole collimation. However, because the projections must be decoded, gains are less than the 5] improvement one might predict based solely on increased counting efficiency. Of course, the bound quantifies this decoding penalty.
6. Discussion and conclusion 5 Although the conventional system does not require explicit reconstruction, some processing — either for resolution-recovery or noise reduction — is often desirable.
While various forms of the CR bound can are useful for assessing detector requirements, there are
N.H. Clinthrone et al./Nucl. Instr. and Meth. in Phys. Res. A 409 (1998) 501—507
507
Fig. 5. Uniform CR bound plots comparing the performance of the ring Compton camera shown in Fig. 4 (solid lines) and parallel-hole collimator/Anger camera (dashed line).
several caveats. First, bound predictions are only as good as the statistical model. If the model does not accurately reflect reality, inferences made from the calculations may not be useful. This is no different than for any design tool — at some point predictions must be validated with measurements. Second, when comparing the performance of two or more designs, the same parameterization of the source must be used. It hardly needs to be mentioned that this parameterization must also be meaningful in regards to the imaging tasks of interest. The judgement of the system designer must guide these modelling decisions; the CR bound cannot itself be used to evaluate the quality of either the statistical model or the parameterization. Finally, there are situations in which the CR bounds provide an overly optimistic estimate of performance. In particular, we have found this true for Poisson measurements when the photon rate is very low or the ML estimator is highly non-linear in the data. These tend to be situations where the propagationof-error method also has poor accuracy. Often alternative lower bounds can be developed [9] when this occurs.
Mindful of the above cautions, CR bounds are invaluable tools for the design of imaging systems and for the assessment of necessary detector requirements.
References [1] R. Blahut, Principles and Practice of Information Theory, Addison-Wesley, Reading, MA, 1987. [2] G.F. Knoll, Radiation Detection and Measurement, Wiley, New York, 1979. [3] J.A. Fessler, IEEE Trans. Image. Process. 5 (3) (1996) 493. [4] A.O. Hero, J.A. Fessler, M. Usman, IEEE Trans. Signal. Process. 44 (8) (1996) 2026. [5] P. Weilhammer, E. Nygard, W. Dulinski, A. Czermak, F. Djama, S. Gadomski, S. Roe, A. Rudge, F. Schopper, J. Strobel, Nucl. Instr. and Meth. A 383 (1996) 89. [6] G.F. Knoll, T.F. Knoll, T.M. Henderson, IEEE Trans. Nucl. Sci. NS-35 (1988) 872. [7] M. Singh, R. Leahy, R. Brechner, T. Hebert, IEEE Trans. Nucl. Sci. NS-35 (1) (1988) 772. [8] N.H. Clinthorne, C.-Y. Ng, C. Hua, J.W. LeBlanc, S.J. Wilderman, W.L. Rogers, Conf. Record of the 1996 IEEE Nuclear Science Symp. and Medical Imaging Conf. 3 (1997) pp. 788—792. [9] N.H. Clinthorne, A.O. Hero, N. Petrick, W.L. Rogers, Nucl. Instr. and Meth. A 299 (1990) 157.
VI. MEDICAL APPLICATIONS