2.1.

Definitions and notation

As noted in the introduction, radiance is a function of spatial position and ray or flux direction. The position of a point in space is specified by a 3D vector R from some arbitrary origin, and the flux direction is specified by a unit vector ŝ. Thus radiance is denoted L(R, ŝ). Because ŝ is a unit vector, it has only two independent components, say s_x and s_y; the third component is given by , so that .

To define the radiance for the given point and flux direction, we construct an arbitrary plane P passing through the point. We denote the normal to this plane by the unit vector . The equation of the plane is · R = constant, so the position of the observation point within the plane is determined by the two components of R perpendicular to . We denote this 2D vector as r; if we choose coordinates so that the reference plane is perpendicular to the z axis (i.e. parallel to ẑ), then R = (x, y, z) and r = (x, y)

With this geometry, the radiance is defined as

where Φ is radiant flux (W), dω is the differential solid angle of a cone of flux directions from a point on the plane, and dA is the differential area on the reference plane P (see Fig. 1). The cos θ in (2.1) is defined as the scalar product between the flux direction ŝ as and the unit vector normal to the reference plane:

Figure 1.

Geometry for specifying radiance on a reference plane. The 3D vector R specifies a point on the plane; the 3D unit vector is perpendicular to the plane, and the 3D unit vector ŝ specifies the flux direction.

The cosine can be associated with either of the partial derivatives in (2.1), leading to two alternative definitions of radiance:

where dA_proj is the projected differential area and dΩ is the projected differential solid angle. If we again choose coordinates so that the reference plane is perpendicular to the z axis, these projections are given by

The two alternative definitions in (2.3) show that radiance can be interpreted either as the flux per unit projected area per unit solid angle or as the flux per unit projected solid angle per unit actual area (on the reference plane). In both cases the word projected refers to projection onto a plane perpendicular to the flux direction through multiplication with cos θ. With either interpretation, radiance is invariant as the orientation of the reference plane is changed.

Since radiant flux is measured in watts (W), the SI units of radiance are watts per square meter per steradian, which we indicate as [L] = W/(m² · sr).

2.2.

Other radiometric quantities derived from radiance

Irradiance, radiant intensity and radiant flux can all be found by appropriate integration if the radiance is known on a reference plane.

Irradiance (power per unit area) incident on the plane P is given by

The 3D observation point R is fully specified by the plane P indicated in the superscript and the 2D vector r in that plane. For parallel to the z axis, for example, E^P(r) is interpreted as the power per unit area on a plane of constant z, and r = (x, y). So long as the plane is specified, r and R determine the same point. The notation E^P(r) indicates that the irradiance depends on (x, y) and on the location of the plane, z. Moreover, tilting the plane tilts the z axis and changes the irradiance; unlike radiance, the irradiance on a plane does depend on the orientation of the plane.

The two integrals in (2.6) run over hemispheres centered on the normal to the plane, so (2.5) applies, but ∫_hemi dω = 2π while ∫_hemi dΩ = π.

Radiant exitance, denoted M^P(r), has the same definition as irradiance, (2.6), but it refers to flux reflected or emitted from plane P rather that incident on it. For a planar Lambertian source or reflector, where the radiance at the plane is independent of direction ŝ, it follows from (2.6) that M = πL.

Radiant intensity, used to express the directional dependence of the flux from some area A on plane P, is given by

The radiant flux (power) through some area A on plane P is given by

2.3.

Spectral radiance

If we wish to specify the spectral dependence of the radiance (per unit wavelength), we define [cf. (2.3)]

We can also define the spectral radiance in terms of the frequency of the radiation, ν = c/λ:

or in terms of the photon energy ε = hν:

2.4.

Flux and radiance in photon units

So far radiant flux has been in SI units of watts or joules per second, but for monochromatic radiation or for the spectral radiances defined above, we can express the flux in photons per second just by dividing the flux in joules/s by the photon energy in joules, denoted as ε_J. Thus, for example, the spectral radiance per unit wavelength from (2.9) becomes a spectral photon radiance, given by

The photon energy in joules is equal to the photon energy in electron volts multiplied by the charge of an electron in coulombs: .

The word photon in this usage has no quantum-mechanical implications; it is just a way of changing the units of flux and spectral radiance.

3.1.

Physical optics

Physical optics deals with coherence and diffraction. The basic difficulty in integrating radiometry with coherence theory and scalar diffraction theory is that scalar fields are complex, while radiometric quantities are real and nonnegative.

For quasimonochromatic light of frequency ν (wavelength λ =c/ν), Adriaan Walther^2–4 defined the generalized spectral radiance (actually a spectral radiance per unit frequency) at a point R on plane P by

where θ is the usual angle between ŝ and the normal to the plane, k = 2π/λ = 2πν/c_m, e(R) is a scalar electric field (e.g., one component of the electric field vector), the angle brackets denote a statistical average over an ensemble of scalar fields, and R and R′ are 3D position vectors specifying points in the plane.

The expectation in (3.1) is a familiar quantity from coherence theory, the ‘mutual intensity, and (3.1) states that (within constants) the generalized spectral radiance is the 2D Fourier transform of the mutual intensity. Thus Walther’s radiance definition makes an important connection between radiometry and coherence.

Another connection, noted by Walther in his 1968 paper, is that the integral in (3.1) is a Wigner Distribution Function (WDF), which originated in quantum mechanics but is now an important tool in phase-space optics] see the book by Testorf et al.⁵ and the review article by Alonso.⁶

By construction, the WDF is real, but it can take on negative values, which of course the radiance defined by (2.1) cannot. Nevertheless, when one uses the generalized radiance to compute irradiance, radiant intensity or flux with the integral formulas from Sec. 2, the results are nonnegative.

3.2.

Quantum optics

A quantum-mechanical definition of radiance can be given by a straightforward generalization of (3.1)¹:

where ê⁺ and ê⁻ are electric field operators and the angle brackets now denote a quantum-mechanical expectation. The other radiometric quantities can then be computed just as in the classical case.

An example of a quantum imaging systems that can be described radiometrically by (3.2) is shown in Fig. 2. In the optical or infrared realm, the source of entangled photons in this figure can be a process called spontaneous parametric down-conversion (SPDC), but in medical imaging with positron-emission tomography (PET) it can be the two entangled gamma-ray photons produced by electron-positron annihilation.

Figure 2.

Quantum imaging with entangled photons.

4.1.

Boltzmann transport equation

The Boltzmann equation is an equation for the time derivative of the spectral photon radiance L_p,ε (R, ŝ, ε, t), or L_P,ε for short. This derivative has contributions from the physical processes of absorption, emission, scattering and propagation radiation, so we have

where the subscripts have the obvious meanings. The four terms in this equation are derived in Sec. 10.3 of Ref. 1, and the resulting transport equation is

where c_m is the speed of light in the medium, μ_tot is the total attenuation coefficient (with contributions from both absorption and scattering processes), is an integral operator describing the angular and energy dependence of the scattering, and Ξ_p,ε describes the photon source. All of these quantities can depend in general on R, ŝ, ε and t.

For a self-luiminous volumetric source, as in fluorescent microscopy or nuclear medicine, the units of the source term are [Ξ_p,ε] = (photons/s) / (m³ · sr · eV). If the source emits isotropically, we can write

where f (R, ε, t) is the volumetric source we wish to image. Note that the 4π in (4.3) carries the units of sr.

For planar Lambertian sources (a common model in optics), the source is confined to a plane (say, z = 0), and the source term is written as

Now f (r, ε, t) is the planar Lambertian source we wish to image. Note that [δ(z)] = 1/m.

If the source is independent of time, or at least slowly varying, we can derive a steady-state form of the Boltzmann equation by setting dL_p,ε/dt = 0; the result is

If we use either of the two source models, (4.3) or (4.4), then (4.5) can be written abstractly as

where f is the object to be imaged and (within constants) is the operator defined by the left side of (4.5).

If we assume there are no sources outside some volume V, then has a left inverse, which we denote by , so that

If we know the source distribution and the characteristics of the propagation medium and/or imaging system, we can calculate the spectral photon radiance at any point outside V.

4.2.

Conservation of radiance along a ray in lossless media

A simple but important example of the Boltzmann formalism is to consider a volumetric, monochromatic, time-independent source f(R) inside a medium with constant speed of light c_m and no absorption, scattering or other radiation sources. For points outside the source region, (4.5) becomes

To solve this equation,¹ we choose a Cartesian coordinate system such that the z axis is parallel to ŝ. Then ŝ · , and (4.8) is an ordinary differential equation in z, which integrates to

If we define ℓ = z – z′ and go back to a general vector notation, we can also write this result as

As illustrated in Fig. 3, the interpretation of this equation is that L_p,ε(R, ŝ, ε) is obtained by integrating the source distribution along a line parallel to ŝ and passing through the point R. The operator that performs this integration is a special case of the operator defined in (4.7), and [f](R, ŝ, ε) is a spectral photon radiance with the indicated arguments.

Figure 3.

Diagram showing that radiance in a homogeneous medium is a line integral of a volumetric source that does not absorb or scatter its own radiation. The cross indicates the origin of coordinates for the three 3D position vectors.

It follows immediately from (4.10) that

where R₁, R₂ and R₃ are three different points along the same line through the source, but all three lying outside the source region (see Fig. 3).

Since this expression holds for all photon energies and all units for the flux, it is valid for the total radiance as well as the spectral photon radiance:

Noting that the three radiances in (4.12) have the same angular argument ŝ, we can call this line a ray, and (4.12) shows that radiance is constant along a ray in a homogeneous, lossless medium. Alternatively, we can define a ray as the path along which the radiance is constant, a path found by solving the Boltzmann transport equation.

Moreover, this same conclusion holds if we consider a lossless optical system, i.e., one with no absorption, scattering or reflection losses, and trace any ray through it as in Fig. 4. Again, radiance is constant along the ray, but now the ray direction depends on which portion of the ray we examine, so we have to write

Figure 4.

Diagram showing the invariance of radiance along a ray in a lossless optical system, as expressed by (4.13). The three points are vectors from a common origin of coordinates, as in Fig. 3, but the vectors themselves are omitted for simplicity in this figure. Note that all three points are in media of the same refractive index; otherwise L/n² would be the conserved quantity.

so long as all three points are in media of the same refractive index (e.g., air). If the points are in homogeneous media of different refractive indices, but again there are no losses, it can be shown that L/n² is the conserved quantity:

4.3.

Conservation of radiance in a lossless paraxial imaging system

ABCD matrices are an elementary tool for analyzing paraxial optical systems. If we consider two arbitrary planes perpendicular to the optical axis in a system with rotational symmetry, the ray height x and the angle of the ray with respect to the optical axis, α, in the two planes are well known to be related by

If the system is paraxial but not rotationally symmetric, then we wirite

where s_⊥ is the projection of the flux direction ŝ onto the plane (i.e., perpendicular to ), and A, B, C and D are all 2 × 2 matrices, so M is 4 × 4. We refer to matrices in this form as ABCD matrices rather than ABCD matrices.

The principle of conservation of radiance in a lossless paraxial system is most easily stated if we write the radiance as

With this notation, the radiance on an output plane is related to the radiance on a parallel input plane by

4.4.

Transport equations for physical optics

A lesser-known application of ABCD matrices is to paraxial diffraction problems. The field propagation from input to output plane is described in the Fresnel approximation by⁷

Friberg⁸ studied the propagation of Walther’s generalized radiance through systems described by ABCD matrices. With notation similar to (4.18), he showed that

Thus the physical-optics radiance is constant along the ray defined by geometrical optics. The same conclusion was reached by Walther (1978) by a stationary-phase approximation, valid asymptotically in the limit of zero wavelength. Friberg’s derivation does not require this limit (though it does use the Fresnel approximation), and it is valid for any state of coherence of the fields.^{1, 6}

A fully quantum mechanical version of Fresnel diffraction theory has been given by Shapiro.⁹ Combining this theory with the quantum-mechanical form of the Walther generalized radiance in (3.1) leads to a transport equation for quantum radiance, at least for paraxial propagation.

5.1.

Integrating detectors

A CCD or CMOS camera produces a pixellated image related to the optically induced electrical charge accumulated at each pixel. A first approximation is that the mean charge (averaged over many images of the same object) is proportional to mean number of incident photons with energy above the bandgap energy during the exposure time. In fact, however, the mean charge produced is influenced also by the spectral, spatial, angular and temporal distribution of the flux. If the spectral photon radiance is known in the detector plane, and if the sensor array responds linearly to the radiation incident on it, the most general expression for the mean output at pixel m has the form¹

where η_m(R, ŝ, ε) specifies the quantum efficiency of the detector, τ is the exposure time and R is a general 3D vector, but the spatial integral is over the two variables needed to specify points in the detector plane.

Other locations for P may also be chosen. Various optical elements such as lenslets, pinholes, spectral filters, choppers or other optical modulators may be placed in the optical path between P and the detector array. The general form of the mean CCD/CMOS output signal in this case is

The function d_m(R, ŝ, ε, t) describes the propagation of radiance through the optical elements between plane P and the detector as well as the response of the detector to the the radiance incident on it. One way to summarize (5.2) is to say each mean pixel output is a linear functional of the radiance in plane P, where the optical elements between the reference plane and the detector give us great freedom in choosing the weighting functions d_m(R, ŝ, ε, t). Cell-phone cameras and DSLRs are now available with 50 MP sensors, which means that we can measure 50 million independent linear functionals of the radiance function in an arbitrary plane.

One application of (5.2), of considerable current interest, is the light-field camera, pioneered by Ren Ng¹⁰ and illustrated in Fig. 5. The microlens array makes the photosensor array sensitive to the angular dependence of the radiance and not just the irradiance, and that allows the image to be refocused on different planes after it is recorded. Plane P in (5.2) is tangent to the left face of the microlens array in Fig. 5, so the light-field camera essentially captures samples of .

Figure 5.

Schematic of a light-field camera (adapted from Ng¹⁰). The object is imaged onto the microlens array, and each microlens images the pupil of the main imaging system onto the photosensor array.

5.2.

Pixellated photon-counting detectors

So-called photon counting detectors have circuitry at each pixel to detect the presence of a photoelectric interaction in the pixel and record the number of interactions (“counts”) during some exposure time for later readout. What is counted, therefore, is photoelectric interactions rather than photons; even a classical electromagnetic field induces photoelectric interactions. Nevertheless, we will continue to use common parlance and refer to a photoelectric interaction as absorption of a photon.

As a practical matter, it is easier to distinguish radiation-induced photoelectric interactions from thermally induced electrons (dark current) if some form of optical gain is used. Intensified CCDs, for example, use an optical image intensifier ahead of the CCD sensor so that a single optical photon absorbed in the photocathode of the intensifier produces a burst of light and hence multiple photons absorbed in the CCD element.

Another way of incorporating gain is to consider an array of avalanche photodiodes (APDs) as the image detector. In this case an absorbed photon produces an electron-hole pair, and each charge carrier is accelerated by an applied electric field to energies greater than the bandgap energy, so it can produce additional electron-hole pairs.

For x rays or gamma rays, the gain can be provided by scintillation materials for which each absorbed high-energy photon produces a burst of low-energy optical photons that can be detected by optical sensors.

Perhaps surprisingly, (5.1) and (5.2) are still valid descriptions of the mean output at pixel m for arrays of photon-counting detector. The actual outputs, g_m without the bar, are integers for photon-counting arrays, but the means are still linear functionals of the radiance so long as the event rate is low enough that there is no loss of counts from temporally overlapping events (pile-up) and no overflow of the registers that store the number of counts at each pixel.

5.3.

Position-sensitive photon-counting detectors

It is also possible to build photon-counting detectors without pixels. A common configuration is shown in Fig. 6, which shows a large-area silicon photodiode with strip electrodes on the four edges. Light incident in a spot, as shown, creates signals in all four electrodes, and a simple centroid calculation provides a measure of the x-y position of the spot.

Figure 6.

Position-sensitive silicon detector.

The primary application of position-sensitive detectors is to track a beam or localized spot of light, but the addition of an image intensifier makes it possible for them to respond to individual photoelectric interactions. As in Sec. 5.2, a single interaction in the photocathode of the intensifier produces a burst of light and hence measurable signals at each of the four output electrodes.

Thus the x-y coordinates of the individual interactions can be computed and stored in real time, as the events occur; no integration time is needed, and there are no detector pixels. Instead, the data consist of a list of the coordinates, possibly tagged with the time at which the event occurred. This form of data storage is referred to as list mode, in contrast to the pixellated or binned mode implicit in (5.1) and (5.2). Listmode data storage has received considerable attention in the radiology literature,^11–13 and it is used routinely in PET.

5.4.

Photon-processing detectors

A photon-absorption event in a detector can be characterized by several attributes, including the position R of the interaction, the energy ε deposited in the detector by the event, the direction ŝ of the radiant flux that produced the event and possibly the time t at which the event occurs. For the j^th event (j = 1…, J), We can organize these attributes into an attribute vector with 5-7 components (2 or 3 components for R_j, depending on the thickness of the detector; one component for energy ε_j; two components for the unit vector ŝ, and possibly time t_j). We denote this attribute vector for the j^th event by A_j.

Of course we cannot directly measure the components of an attribute vector, but we can estimate them from suitable data sets. For example, if the position-sensitive detector of Fig. 6 is operated with an image intensifier, the pulse signals from the four electrodes can be used to estimate x_j and y_j, which are two of the 5-7 components of A_j. To get sufficient information to estimate additional components, we need to add optical elements such as lenslets or color filters in front of the detector to alter its sensitivity to the spectral photon radiance [see (5.2) and Fig. 5).

The goal of a photon-processing detector is estimate the components of each attribute vector as accurately as possible. We denote the estimated attributes with tildes. Thus, if A_j = (R_j, ŝ_j, ε_j), then .

Of course these estimates will not constitute perfect measurements of the true photon attributes, but statistical estimation theory^{1, 14, 15} provides us with the tools to analyze and minimize the bias and variance of the estimates with respect to the true attributes.

An efficient estimator is one that is unbiased and has the smallest possible variance, as specified by the Cramèr-Rao lower bound. It is well known from statistical estimation theory that the maximum-likelihood (ML) estimator is efficient in any problem where an efficient estimator exists, and in any case the ML estimator is asymptotically efficient as more or better data are acquired about each photon-absorption event.

For a review of ML estimation as applied to photon-processing detectors, see Barrett et al.¹⁶

To summarize what we know about how to implement optimum estimation of photon attributes, an ideal photon-processing detector must:

Figure 7 shows a conceptual diagram of a photon-processing detector capable of operating in the visible or NIR portions of the spectrum. A practical realization of a photon-processing detector for x rays, gamma rays and charged particles is the intensified Quantum Imaging Detector (iQID) developed by Brian Miller.^{17, 18} Computer systems for performing the MLE step in real time are discussed by Hesterman at al.¹⁹

Figure 7.

Generic photon-processing detector for optical wavelengths.

6.1.

List-mode data and point processes

In practice, the estimated photon attributes will be stored in listmode, where each list entry is a K-dimensional vector, with K being the number of estimated attributes per event; K = 6 in our example where 3D position, two components of a unit vector and an energy are estimated.

A mathematically equivalent representation of the same information is a random point process, a sum of delta functions defined by

For a discussion of angular delta functions, , see Sec. 2.4.7 in Barrett and Myers,¹ and for a thorough treatment of random point processes, see Sec. 11.3.

Because of the independence of the photon-absorption events, the number of events J is a Poisson random variable, and u(R, ŝ, ε) is called a Poisson random process.

The mean of this point process is obtained by averaging over the statistics of each of the J vectors A_j and over J itself for a given object f that produced the data. The result is a function of A which depends on f, so we can denote it as ū(A|f), or more explicitly as ū(R, ŝ, ε|f). Technically, this function can be regarded as a vector ū(f) in the Hilbert space 𝕃₂(ℝ^K), which is the space of square-integrable functions of K continuous variables.

6.2.

Linear imaging of emissive objects

For self-luminous objects as in astronomy, nuclear medicine or fluorescence microscopy, Caucci and Barrett¹³ have shown that

where is a linear operator that maps 𝕃₂(ℝ³) to 𝕃₂(ℝ^K). We emphasize that is linear even when the relation between raw data (e.g., sensor signals) and estimated attributes is nonlinear. Linearity depends only on the independence of the individual photon interactions.

For comparison, all linear photon-counting imaging systems that record only the number of counts in each of M bins are described by

where now maps 𝕃₂(ℝ³) to an M-dimensional Euclidean space 𝔼_M. In the terminology of Barrett and Myers,¹ is a continuous-to-discrete (CD) operator but is a continuous-to-continuous (CC) operator. The crucial difference between CC and CD operators is that the latter must necessarily have an infinite-dimensional null space; the system maps a vector in an infinite-dimensional Hilbert space to a finite set of numbers, so there is an infinite set of objects f_null that yield no data at all: . Because real objects almost always have null components for any CD system, the object itself cannot be recovered from the discrete data, and even simple linear functionals such as integrals of the object over voxels cannot be estimated without a possibly infinite bias. A CC operator like may have a null space, but if K ≥ 3 the null space is not demanded by dimensionality considerations, so it is possible that many of the indeterminacy problems inherent in the world of conventional digital imaging may disappear with photon processing rather than mere photon counting.

6.3.

Nonlinear imaging of transmissive objects

In many situations in optics and radiology, an image is formed by illuminating a transmissive object with a known radiation source. In these cases, the object is described by the spatial and spectral distributions of its attenuation and scattering coefficients. In terms of the Boltzmann equation (4.5), the object properties are embedded in the attenuation term c_mμ_totL_p,ε and the scattering term rather than the source term c_mΞ_p,ε. The solution of the Boltzmann equation shows that the radiance is exponentially attenuated as it passes through the object, so the radiance on the image detector depends nonlinearly on the object.

We can still use (6.3) to describe the mean image data with a photon-counting or integrating detector, and (6.2) for a photon-processing detector, but now the operators and are nonlinear, and it is much more difficult to say anything about null functions.

6.4.

Probability density functions for photon-processing data

The full probability law for a univariate Poisson random variable is determined by its mean. The same thing is true for a Poisson random process but now the mean is a function, ū(R, ŝ, ε) and the probability law is a probability density function or PDF.

The 6D PDF for the attribute estimates is obtained by normalizing the mean:

Because the events are statistically independent (conditional on f), the overall PDF for the entire data set can be expressed as a product of 6D PDFs like those in (6.4).^{11, 13}

6.5.

Computation of the conditional mean of the Poisson Point Process

According to (6.4), the mean of the Poisson point process tells you everything about the statistics of the process, conditional on a a specific object f, but how do we compute that mean?

The first approach is to recognize that, within constants, the mean ū(R, ŝ, ε|f) is just a blurred version of the spectral photon radiance in plane P, , with the blurring arising from the inaccuracies in estimating the photon attributes from the sensor data. Ways of accounting mathematically for this estimation blur have been discussed by Caucci and Barrett,¹³ so finding the mean of photon-processing data, and hence the full PDF by (6.4), boils down to computing the spectral photon radiance in plane P.

As we saw in Sec. 4.1, the spectral photon radiance is the solution of the Boltzmann transport equation, and often analytical solutions can be found; a simple example is (4.10). In many other situations in a plane of interest P can be computed from the radiance in another plane, such as a source plane, by simple invocation of the principle of conservation of radiance along a ray, as discussed in Secs. 4.2 and 4.3. As we know from Sec. 4.4, these conservation and propagation laws apply to physical-optics and quantum-optics versions of radiance also.

For emissive objects, it may be possible to approximate the linear operator with a very large matrix and find its elements by simulation of he imaging system. For transmissive objects, is nonlinear, but nevertheless the Boltzmann can be solved to get the spectral photon radiance on the detector and from there get the mean data for any of the four detection modes discussed in Sec. 5. For an example of this approach as applied to diffuse optical imaging, see Jha et al.^20–22

When analytical solutions are not readily available, we can use Monte-Carlo simulation to estimate the spectral photon radiance. There is even a version of Monte-Carlo simulation, called adjoint Monte Carlo,¹ that allows direct estimation of the spectral photon radiance blurred by the estimation errors.

7.1.

Task-based assessment

In much of the medical-imaging community, and increasingly in other areas of imaging, image quality is defined in terms of the performance of some observer on some clinically or scientifically relevant task. Generically, the task can be either classification or estimation. An estimation (or quantitation) task is to determine a numerical value for one or more parameters associated with the object being imaged. Examples include cardiac ejection fraction or the concentration of a dye or contrast agent in some region of interest. In a classification task, the goal is to assign the object to one of two or more classes. A simple example is tumor detection, where the classes are tumor-present and tumor-absent. The observer for a classification task can be either a human or a mathematical algorithm. The performance of a human observer can be measured by well-established psychophysical methods and analysis of the resulting ROC (receiver operating characteristic) curves. Psychophysical studies are complicated and time-consuming, however, so there is considerable interest in developing mathematical model observers that will predict human performance. Another form of mathematical observer is the ideal observer, which makes optimal use of the available data in such a way as to maximize an ROC-based figure of merit. If we know the performance of the ideal observer, we can compare it to the performance of a human and determine how well the human can extract the information contained in the data. Often the performance of the true ideal observer is difficult to compute, so we use the optimal linear observer, which we call the Hotelling observer. It is optimal among all observers constrained to perform only linear observations on the data. An important variant on the Hotelling observer is the channelized Hotelling observer (CHO), which is the Hotelling observer applied not to an image but to a reduced data set consisting of the sampled outputs of a small number of linear filters, called channels. An important example is spatial-frequency-selective channels such as those known to be operative in the human visual system.

7.2.

Object variability

For all of these tasks and observers, task performance is limited by the randomness in the data. So far we have discussed randomness arising from photon noise in different kinds of detectors, but there is also an inevitable statistical variability in the objects being imaged,²³ and in some cases the imaging system itself introduces a form of noise.²⁴

To account for the object variability in a classification task, we have to understand the statistics of the data conditional on a class C of objects. For any data set g, including the listmode data acquired by photon-processing detectors, it follows from elementary laws of probability that

where pr(g|f) is the PDF of the data given a particular object f, and the expectation is over all objects drawn from the class.

To visualize this expectation and its impact on classification performance, we can imagine using an accurate simulation program to generate N sample objects from class C. Denoting the n^th sample object as f_n (n = 1, …N), we can approximate pr(g|C) by a sample average as

For any form of photon noise, which is the subject of this paper, the PDF pr(g|f_n) is fully determined by its conditional mean ḡ(f_n). For a photon-counting detector with M pixels, therefore, (6.3) and (7.2) show that

where g is a vector in an M-D space with components (g₁, g₂, g₃, …,g_M). These components are statistically independent conditional on an object f, so pr(g|) is a product of M univariate Poissons, but the average over objects in (7.1) introduces statistical dependence. The accuracy of the approximation in (7.3) increases as the number of sample objects, N, increases.

To display this M-D PDF, we consider two components at a time. For example, we can plot the bivariate PDF pr(g_m,g_m′|C) as in Fig. 8 for different total exposure times τ. The axes on these plots represent g_m and g_m′ (not m and m′), and the gray level is the PDF. The leftmost panel is for the largest τ, and the rightmost panel is for the smallest τ.

Figure 8.

Illustration of the bivariate class PDF pr(g_m, g_m′|C) for a class of objects C and a photon-counting pixel array. Distributions are shown for three different values of the exposure time τ (increasing right to left). See text for explanation.

The points shown by x on all three panels correspond to the values of for the N sample objects used in (7.3). For this illustration, N = 20, and the bivariate PDFs for different objects are well separated on the leftmost panel where τ is large so the standard deviation of the counts per pixel is small compared to its mean. For smaller τ, however, the standard deviation of the counts on each pixel increases relative to its mean, and the individual PDFs, pr(g_m, g_m′|) merge into one big blob, which for large N will be a good approximation to the class PDF, pr(g_m, g_m′|C). The overall variance and covariance for the class PDF depends on both the exposure time (mean number of counts per pixel for a given object) and the randomness of the objects themselves.

A similar illustration holds for photon-processing detectors, but there the discrete index m is replaced by continuous indices like (R, ŝ, ε). The exposure time τ controls the total number of recorded events J in the data list for a particular object, and again the overall variance and covariance for the class PDF depend on both the exposure time and the randomness of the objects. For mathematical details, see Caucci.¹³

7.3.

Information content of a photon

Consider an imaging problem where the task is to decide if the object that produced a particular image g was from class C₁ or class C₂. To perform this task, we need some knowledge the two class PDFs, pr(g|C₁) and pr(g|C₂). These densities can again be approximated by expressions like (7.3) if we can accurately simulate objects from the two classes.

A depiction of the data statistics like that in Fig. 8, but for two classes, is given in Fig. 9. It is visually evident that the classification task can be performed with high accuracy at large exposure time (leftmost image in the figure) but not for shorter exposures where the photon noise dominates.

Figure 9.

Illustration of the bivariate PDFs, pr(g_m, g_m′|C₁) and pr(g_m, g_m′|C₂), for two classes and a photon-counting pixel array. Distributions are shown for three different values of the exposure time τ (increasing right to left). See text for explanation.

It is beyond the scope of this paper, but it is possible to make these statements much more precise by specifying the observer that will perform the tasks (e.g., human or Hotelling observer) and by adopting some measure of the task performance, such as area under an ROC curve. If one then plots this performance metric as a function of the mean number of absorbed photons, then one can define the information content per photon for this particular imaging system, task and observer. For a full mathematical exposition of this approach, with application to imaging with ionizing radiation, the reader can consult a forthcoming tutorial review article in Physics in Medicine and Biology.²⁵