1. Introduction

Among the most important characteristics of many imaging, communication and measurement systems are their capacity to encode or transmit large quantities of information, resolve sufficiently fine spatial details and provide output with high signal-to-noise ratio (SNR) at lowest possible radiation doses or input energy levels [1, 2]. For mainly historical reasons (abundance of photons in typical visible light imaging applications), imaging performance characteristics related to the incident particle fluence have not been extensively studied in classical optics. However, such characteristics, and the interplay between them, have attained additional relevance in recent years in the context of biomedical imaging, where the samples are often sensitive to the radiation doses [3], in certain astronomical methods where the detectable photon flux can be extremely low [4], as well as in some problems related to foundations of optics and quantum physics [5, 6]. In X-ray medical imaging, in particular, it is critically important to minimize the radiation dose delivered to the patient, while still being able to obtain 2D or 3D images with sufficient spatial resolution and SNR in order to detect the features of interest, such as small tumours [7, 8]. In this context, an imaging system (e.g. a CT scanner) must be able to maximize the amount of relevant information that can be extracted from the collected images, while keeping sufficiently low the number of X-ray photons impinging on the patient. The present paper addresses some mathematical properties of generic imaging systems, that are likely to be important in the context of designing medical imaging instruments, and may also have relevance to some fundamental aspects of quantum physics and information theory.

We have recently introduced a notion of intrinsic imaging quality, which incorporates both the spatial resolution and the noise propagation properties of an imaging system [6–10]. Note that in our treatment we consider all the components, beginning from the source of imaging particles/radiation and ending with the optional computer-aided post-processing of the images, to be parts of a complete “imaging system”, with all these components potentially affecting the overall imaging quality of the system. Compared to previous publications on this and related topics, Section 2 of this paper extends the notion of intrinsic imaging quality from “flat-field” (featureless) images to the situation were the signal and noise levels can vary across the images, but the rate of variation is slow compared to that of the point-spread function (PSF) of the imaging system. In Section 3, we consider the notion of information capacity of imaging systems, leading to results which are consistent with the ones previously reported on the basis of the classical approach originated by C. Shannon in the context of communication systems [11] and later extended by others to imaging systems [12–17]. These results are also related to the theory of communication channels and modes [18–22], as discussed in the Appendix. In Section 4, we apply the obtained results for the intrinsic imaging quality and information capacity of linear shift-invariant (LSI) systems to analysis of a simple model of quantitative imaging of weakly scattering samples. In particular, we derive an estimate for the maximum information about the sample that can be extracted in the relevant scattering experiments. We show that such a limit for the extractable information depends only on the total number of imaging particles used (i.e. on the radiation dose), the number of spatial resolution units and the “scattering power” of the sample, which is defined as a product of the average scattering coefficient and the linear size of the sample. The main results of the present paper are summarized in Section 5.

2. Noise and spatial resolution

Consider a simple model for formation of images in an LSI system which involves a stochastic flux of imaging particles incident on a position-sensitive detector, and a linear filter function describing the deterministic post-detection action of the imaging system. Here the fluence of particles (i.e. the particle density) in an image registered by the detector will typically result from interaction of particles emitted by a source with a scattering object (sample) that is being imaged, but we leave the analysis of effects of such interaction until Section 4 below.

Different realizations (instances) of images registered by the detector under the same conditions correspond to different sample functions of the stochastic detected particle fluence distribution $S_{in}(x)$, $x\in {ℝ}^n$, where n is the dimensionality of the image. Note that here we generalize the concept of a “detection system” to any positive integer dimension n, which can be convenient e.g. for description of 3D imaging. We will call the expectation value (ensemble average) ${\overline{S}}_{in}(x)$ the “detected signal” (which serves as an “input” for the post-detection signal processing system). The difference ${\tilde{S}}_{in}(x)=S_{in}(x)-{\overline{S}}_{in}(x)$ is the corresponding “noise process”.

The second-order correlation properties of the particle fluence registered by the detector are described by the autocovariance function $\Gamma (x,y)=\overline{{\tilde{S}}_{in}(x){\tilde{S}}_{in}(y)}$. We assume that the autocovariance function has the following form:

As the noise power spectrum $W_x(\xi )$ is a non-negative and even function of ξ, the same is true for $\widehat{H}(\xi )$. Therefore, there exists a real-valued even function $P(x)$, such that $\widehat{H}(\xi )=|\widehat{P}(\xi ){|}^2$, and hence $H(x)={\displaystyle \int P(y)P(x+y)dy}$. The function $P(x)$ represents the “native” point-spread function (PSF) of the imaging system, which can reflect both the correlation properties of the incident particle fluence (i.e. the properties of the source and the imaging set-up) and the PSF of the detector. We assume here that $P(x)$ is non-negative. Then, taking into account that $\widehat{H}(0)=\left|\right|H|{|}_1=|\left|P\right|{|}_1^2$ and $H(0)=\left|\right|P|{|}_2^2$, where $\left|\right|f|{|}_p=({\displaystyle \int |f(x){|}^{p}dx{)}^{1/p}}$ is the usual L_p norm, we obtain in particular that $W_x(0)={\displaystyle \int {\Gamma }_x(h)\text{d}h}={\sigma }_{in}^2(x)\widehat{H}(0)/H(0)={\sigma }_{in}^2(x)\left|\right|P|{|}_1^2/|\left|P\right|{|}_2^2.$According to its definition and the properties of the covariance function assumed above, the native PSF $P(x)$ is much narrower than both the detected signal ${\overline{S}}_{in}(x)$ and the noise variance ${\sigma }_{in}^2(x)$. Note that such a behavior is “natural” for a well-designed imaging system: it ensures, in particular, that the detector PSF does not significantly smear typical incident signals. Model examples of the native PSF are represented by $P_{Gauss}(x)={(2\pi )}^{-n/2}{\sigma }^{-n}\exp [-|x{|}^2/(2{\sigma }^2)]$ and $P_{rect}(x)={(2\sigma )}^{-n}{\chi }_{{[-\sigma ,\sigma ]}^n}(x)$, where ${\chi }_{{[-\sigma ,\sigma ]}^n}(x)$ is the function that is equal to one in the cubic area around the origin of coordinates with side length equal to 2σ, and is equal to zero outside this area. An autocovariance function of the type described above is obtained, for example, in the case of a filtered Poisson point process [1] when the width of the filter is small compared to the typical variation length of the mean photon fluence.

We want to find out how the output of the imaging system changes after the action of an LSI transformation with a non-negative PSF (filter function) $T(x)$,

Having introduced the notion of SNR in the LSI system defined by Eq. (3), we now turn to the notion of spatial resolution. The intrinsic spatial resolution Δx of a system described by Eq. (3) can be defined via the variance of its PSF as follows,

Table 1. Spatial resolution associated with different PSFs. All PSFs are normalized such that $\left|\right|T(x;\sigma )|{|}_1=1$.

View Table

Following the ideas discussed in [2, 23] and elsewhere, it is also possible to define the spatial resolution of the LSI system by an alternative expression:

Combining Eqs. (5) and (7), we can obtain:

If the detector area Aⁿ is fixed, the spatial resolution “area” and the total number $M_2$ of resolvable effective “pixels” (resolution units) are inversely proportional to each other: $M_2=A^n/{({\Delta }_{2}x)}^n$. Then Eq. (8) implies that the product of the number of resolvable spatial units (pixels) and the square of the SNR is also invariant under the action of LSI systems with suitable PSFs:

Note that, while the invariant defined by Eq. (9) is dimensionless, the invariant defined by Eq. (8) has the dimensionality of inverse length taken to the power of n. One way to “naturally” normalize this function and make it dimensionless is to divide it by ${\overline{S}}_{in}(x)$. Indeed, in the case of a filtered Poisson point process with the filter function $T(x)$ which is much narrower than the variation length of detected signal ${\overline{S}}_{in}(x)$ and the noise variance ${\sigma }_{in}^2(x)$, one has [1, Sect.11.3.9]: $\overline{S}(x)=(T\ast {\overline{S}}_{in})(x)\cong \left|\right|T|{|}_1{\overline{S}}_{in}(x)$ and ${\sigma }^2(x)=(T^2\ast {\overline{S}}_{in})(x)\cong \left|\right|T|{|}_2^2{\overline{S}}_{in}(x)$. Then using Eq. (7) we obtain:

We have previously introduced and studied [6–10] the following dimensionless quantity, which is close in form to Eq. (10) and incorporates both the noise propagation and the spatial resolution properties of LSI system:

We will also need an estimate for the average value of $SN{R}^2(x)$ over all spatial resolution units in the image. In the case of incident fluences with Poisson statistics, we obtain from Eq. (10): $SN{R}^2(x_m)={(\Delta x_2)}^n{\overline{S}}_{in}(x_m)=\overline{n}(x_m)$, which can be identified with the average number of imaging particles registered in the “pixel” (spatial resolution unit) with the center at $x=x_m$. Let us sum up this quantity over all points $x_m,m=1,2,\mathrm{...},M_2$, corresponding to the centers of the output spatial resolution units: $M_2〈SN{R}^2〉={\displaystyle \sum _{m=1}^{M_2}SN{R}^2(x_m)}={\displaystyle \sum _{m=1}^{M_2}{(\Delta x_2)}^n{\overline{S}}_{in}(x_m)}={\displaystyle \sum _{m=1}^{M_2}\overline{n}(x_m)}=N$, i.e.

A useful expression for the intrinsic imaging quality characteristic can be obtained by first noting that $Q_S^2={({\Delta }_{2}x)}^n/{(\Delta x)}^n=M/M_2$, where $M=A^n/{(\Delta x)}^n$, and then substituting the expression $M_2=N/〈SN{R}^2〉$ obtained from Eq. (13), with the result given by:

3. Information capacity of imaging systems

Consider a simple “imaging” experiment that involves N incident particles (“imaging quanta”, e.g. photons) and a 2D detector with M pixels (where both N and M are non-negative integer numbers). We would like to calculate first how many distinct “images” can be produced in this situation, in other words, how many distinct sequences $(n_1,n_2,\mathrm{...},n_M)$ exist, with ${\displaystyle \sum _{m=1}^M{n}_m=N}$ (all n_m are non-negative integers). Two sequences $(n_1,n_2,\mathrm{...},n_M)$ and $({n^{\prime }}_1,{n^{\prime }}_2,\mathrm{...},{n^{\prime }}_M)$ are considered distinct if $n_m\ne {n^{\prime }}_m$ for at least one index m, $1\le m\le M$. This problem is equivalent to the one about the number of different ways for distributing N indistinguishable (identical) balls into M distinguishable bins. The corresponding total number of different combinations, that we denote by $X(M,N)$, can be found by utilizing the following recursive relationship:

Now let us consider the possibility that some of the particles incident on the object can be absorbed or strongly deflected by the object, etc., and as a result will not be registered by the detector. This means that some of the resultant “images” in this situation can be created with a different number $N^{\prime }$ of registered quanta, where $N^{\prime }$ can be any integer satisfying $0\le N^{\prime }\le N$. Let us calculate the number, $Y(M,N)$, of all distinct images that can be formed this way. Obviously,

A simpler approximation for the expression $Y(M,N)$ in Eq. (19) can be derived for large values of M and N. According to the refined version of Stirling's formula for the factorial [26]:

Let us introduce a notion of “information capacity” $D(M,N)={\log }_{2}Y(M,N)$ by analogy with the original information capacity definition introduced by Shannon [11] for communication systems and later extended to imaging systems [12–17]. We will also consider information capacity per single particle: $D_S(M,N)=D(M,N)/N=(1/N){\log }_{2}Y(M,N)$. From Eqs. (19)-(20) we obtain:

Usually in practice the number of imaging particles (e.g. photons) in an image is much larger than the number of detector pixels, which corresponds to the case $N>>M$ in the above equations. In that case Eq. (21) leads to the following expressions

Let us now postulate that two images are distinguishable if and only if the difference between their signal values in at least one pixel exceeds the sum of standard deviations of noise in that pixel in the two images (this can be easily modified to include a distinguishability factor SNR_min which can e.g. be equal to 5, if the Rose criterion is applied). In the “noise-free” case considered above, the distinguishability was equivalent to any difference in detected particle numbers in the two images by at least one particle in at least one pixel, which formally corresponds e.g. to the standard deviation of noise being equal to 1/2 or less in all pixels of both images.

Consider the case where the detector is not perfect and contains a fixed amount of additive noise in each pixel. Assuming that the level of noise is σ, two signals, corresponding to n₁ and n₂ registered particles in a pixel, respectively, can be considered distinguishable e.g. if $|n_1-n_2|\ge 2\sigma$. This is equivalent to $|n_1/(2\sigma )-n_2/(2\sigma )|\ge 1$. In other words, the number of distinguishable images can now be calculated in terms of “bunches” of $2\sigma$ particles, instead of the individual particles (note that the fractions of the bunches, that correspond to signal increments below the level of distinguishability, may lead to a reduction of the total number of full bunches). Therefore, the result obtained above in Eqs. (22)-(23) should remain valid in the case of additive detector noise, with the replacement of N by $N/(2\sigma )$. Note also that because ${\log }_2(N/(2\sigma M))={\log }_2(N/M)-{\log }_2(2\sigma )={\log }_2(N/M)+O(1),$ and hence the form of Eqs. (22-23) can remain unchanged in this case.

Finally consider the case of a photon-counting detector, which corresponds to the Poisson statistics with the expectation and variance both equal to the mean number $\overline{n}$ of photons per pixel, $\sigma ={\overline{n}}^{1/2}$. Let us define two signal levels as distinguishable if they differ at least by the sum of their standard deviations. If the number of photons in a pixel is equal or smaller than $\overline{n}$, then there are ${\overline{n}}^{1/2}$ distinguishable possible signal levels. For example, for $\overline{n}=25$, the ${\overline{n}}^{1/2}=5$ distinguishable levels (bands) correspond to 25 ± 5, 16 ± 4, 9 ± 3, 4 ± 2, and 1 ± 1 photons. Therefore, the result stated in Eqs. (22) and (23) should remain valid in the case of Poisson statistics, with the replacement of N by $N^{1/2}$. Note also that $M{\log }_2(\sqrt{N}/M)+O\left(M\right)=0.5M{\log }_2(N/M)+O\left(M{\log }_{2}M\right).$Although the above considerations regarding the extension of Eqs. (22) and (23) to the case of noisy images are not rigorous, it turns out that they lead to correct formulae for the information capacity in the case of Poisson (shot) and Gaussian (additive) noise. A rigorous derivation for these cases is given in the Appendix where more accurate formulae are derived from the results in [18, 20, 22]. Specifically, an estimate for the information capacity for the case of Poisson statistics when $N/M$ is large, is given by

Identifying the individual “pixels” with the isotropic spatial resolution units having the side length given by Eq. (7) and substituting Eq. (13) into Eq. (24), we obtain that in the case of Poisson noise:

4. Effect of the imaged sample

The above analysis has been carried out without a reference to an imaged sample. In the present section we consider one simple generic model which takes into account the interaction of the incident particle fluence with a sample and its consequences for the information about the sample that can be obtained in an imaging experiment.

Consider an experiment involving interaction of a flux of incident imaging particles with a weakly scattering sample, after which the scattered particles are registered by a position-sensitive detector. We assume for simplicity that the particle fluence over the “incident surface” of the sample can be described by its spatially uniform expectation value ${\overline{S}}_0^{(n-1)}$. Note that the latter quantity is expressed in the number of particles per unit “area”, i.e. it has the dimensionality of inverse length in the power of n - 1, while the detected signal, ${\overline{S}}_{in}(x)$, introduced in Section 2 was expressed in the number of particles per unit “volume”, i.e. it has the dimensionality of inverse length in the power of n. We assume here that the detected scattered fluence $S_{in}(x)$ has a Poisson distribution with the mean ${\overline{S}}_{in}(x)={\overline{S}}_0^{(n-1)}\gamma (x)$. This corresponds to the case of so-called “dark-field imaging”, when the primary unscattered beam does not reach the detector. Modification to the case of bright-field imaging is quite straightforward and the corresponding result is given below. Here $\gamma (x)$ is a deterministic non-negative function (“scattering coefficient”) that describes the interaction of the incident fluence with the sample. This interaction is assumed to be linear with respect to the interaction length, i.e. $\gamma (x)$ has the dimensionality of inverse length. As a model for $\gamma (x)$ one may consider, for example, the linear attenuation, $\mu (x;\lambda )=(4\pi /\lambda )\beta (x;\lambda )$, or refraction, $(2\pi /\lambda )\delta (x;\lambda )$, coefficients in X-ray imaging, where $n=1-\delta +i\beta$ is the complex refractive index of the sample and λ is the X-ray wavelength (we omit the argument λ in the notation for $\gamma (x)$ for brevity). It can sometimes be helpful to express $\gamma (x)=\sigma \rho (x)$, where σ is the scattering cross-section and $\rho (x)$ is the number density of scattering centers (usually, atoms) in the sample. The interaction is also assumed to be weak, so that $\gamma (x)A<<1$, where A is the linear size of the sample. In particular, this allows us to ignore the weakening of the incident flux in the course of its propagation through the sample. Finally, the interaction has a “single-scattering” nature, i.e. each incident particle interacts at most with a single scattering center the sample, as is usually assumed in the first Born approximation or in the so-called “kinematical approximation” in X-ray and electron diffraction [27]. This model is somewhat similar to the one used in diffraction tomography [28], where the first Born approximation is invoked for the description of interaction of the radiation with the sample. However, unlike diffraction tomography, our simple model ignores the effect of propagation of the scattered radiation to the detector, assuming only that there is a one-to-one correspondence between the resolution volumes in the sample and the effective “detector pixels”. The totality of detector pixels may possibly be aggregated over multiple 2D exposures of the sample. For example, in conventional or diffraction CT imaging, such one-to-one correspondence may be achieved after a reconstruction procedure, e.g. after the application of the inverse Radon transform to the data collected in real detector pixels at different rotational positions of the sample. In the case of weakly absorbing or weakly scattering samples, such reconstruction procedure is usually linear and can be expressed by a convolution with the corresponding filter function [29], and hence these cases fit well into the LSI approach considered above.

Consider now a post-detection filtering of the Poisson process with a narrow filter function $T(x)$ as in Section 2 above (i.e. the filter function is required to be narrow with respect to the characteristic variation length of the scattering properties across the sample). We then obtain: $\overline{S}(x)=(T\ast {\overline{S}}_{in})(x)\cong \left|\right|T|{|}_1{\overline{S}}_0^{(n-1)}\gamma (x)$ and ${\sigma }^2(x)=(T^2\ast {\overline{S}}_{in})(x)\cong \left|\right|T|{|}_2^2{\overline{S}}_0^{(n-1)}\gamma (x)$. Accordingly,

Using the expression for $SN{R}^2(x)$ from Eq. (28) and arguing similarly to the derivation of Eq. (13) above, i.e. averaging over all spatial resolution units into which the sample is split, we can obtain

Let us consider for comparison a more specific example of conventional computed tomography (CT) imaging. It was shown in [30] that in the case of filtered back-projection reconstruction (FBP) with the nearest-neighbor interpolation, one has:

Using Eq. (29) in combination with Eq. (26) for the information capacity, we obtain:

Substituting the expression for the SNR² from Eq. (28), together with the expression $S_{in}(x)={\overline{S}}_0^{(n-1)}\gamma (x)$ for the incident fluence, into the definition of $Q_S^2$ in Eq. (11), we find that the imaging quality here has exactly the same form as in the case with no sample considered in Section 2. In particular, Eq. (12) and the estimate $Q_S^2\le 1/C_n$ remain unchanged. This happens because, according to Eq. (12), the imaging quality can be expressed in terms of the two definitions of the system's spatial resolution alone, and these resolutions do not depend on the effect of the scattering power of the sample under the conditions of the simple model considered here.

5. Conclusions

We have considered the relationship between the spatial resolution, signal-to-noise ratio and information capacity of linear shift-invariant imaging systems. Convenient imaging performance characteristics of such systems include the imaging quality $Q_S$ and the information capacity per single imaging particle (e.g. photon), D_S. From the results of this paper it is clear that, apart from relatively insignificant factors, the information capacity per single particle and the intrinsic imaging quality have similar mathematical expressions, namely $D_S(M,N)~(M/N){\log }_2<SN{R}^2>$ and $Q_S(M,N)~(M/N)<SN{R}^2>$, respectively. However, the behavior of these two quantities under linear transformations of the imaging system is different due to the presence of the logarithm in the first expression. It is usually possible to increase the signal-to-noise (SNR) ratio at the expense of spatial resolution (or equivalently, the total number of spatial resolution units, M) and vice versa using operations like image filtering, denoising or deconvolution. While the intrinsic imaging quality remains essentially invariant under such (linear) transformations, the information capacity can change significantly: for example, a system with a single spatial resolution unit has quite low information capacity per particle, ~$({\log }_{2}N)/N$. On the other hand, when the number of particles is fixed, the information capacity per particle can grow without a limit when the number of resolution units grows: consider the fact that an unlimited number of new images can be created, in principle, by directing a photon or a fixed number of photons into new pixels (one pixel at a time).

In imaging problems that involve quantitative analysis of samples on the basis of transmitted or, more generally, scattered radiation, the notion of “useful signal” is modified. We have studied one example of this type of problems in Section 4 of the paper. We have found that the information capacity is proportional to the logarithm of the product of the total number of imaging particles and the scattering power of the sample, the latter depending only on the average scattering coefficient and the size of the sample (or alternatively, on the total number of scattering centers, e.g. atoms, in the sample, and the average relative scattering cross-section of each such center). On the other hand, the imaging quality per single particle is not affected by the scattering properties of the sample, and remains the same as in the sample-free case. It means that, in accordance with its name and the intended role, the intrinsic imaging quality characteristic depends only on the properties of the imaging system itself, and does not depend on the imaged sample.