research papers
accessSignal-to-noise and spatial resolution in in-line imaging. 2. Phase-contrast tomography
aSchool of Physics, University of Melbourne, Parkville, Victoria 3010, Australia, bSchool of Physics and Astronomy, Monash University, Clayton, Victoria 3800, Australia, cSchool of Physical & Chemical Sciences, University of Canterbury, Christchurch, New Zealand, dSchool of Science and Technology, University of New England, Armidale, New South Wales, Australia, and eAustralian Synchrotron, Australian Nuclear Science and Technology Organisation, Sydney, New South Wales, Australia
*Correspondence e-mail: [email protected], [email protected]
In the first part of this paper, quantitative aspects of propagation-based phase-contrast imaging (PBI) were investigated using theoretical and numerical approaches, as well as experimental two-dimensional PBI images collected with plane monochromatic X-rays at a synchrotron beamline. In this second part, signal-to-noise ratio, spatial resolution and contrast are studied in connection with the radiation dose in three-dimensional PBI images of breast tissue samples obtained using propagation-based phase-contrast computed tomography (PB-CT) with energy-integrating and photon-counting detectors. The analysis is based on the theory of PBI and PB-CT using the homogeneous Transport of Intensity equation (Paganin's method). A biomedical X-ray imaging quality characteristic, suitable for quantitative assessment of X-ray images of biological samples, is introduced and applied. The key factors leading to high values of the biomedical X-ray imaging quality in PBI and to relatively low values of the same quality metric in CT imaging are identified and discussed in detail. This study is aimed primarily at developing tools for quantitative assessment and optimization of medical PB-CT imaging, initially at synchrotron facilities, with the prospect of subsequent transfer of the technology to medical clinics.
Keywords: X-ray imaging; computed tomography; phase contrast; spatial resolution.
1. Introduction
In the first part of this paper (Gureyev et al., 2024
), we presented theoretical results, numerical simulations and analysis of experimental data collected using propagation-based phase-contrast X-ray imaging (PBI) at a synchrotron source. Background information about PBI and its general advantages over conventional attenuation-based X-ray imaging were also discussed in the first part of this paper and therefore will not be repeated here. The theoretical approach used there, which is also used in the present second part of the paper, was largely based on the noise-resolution uncertainty (NRU) relationship (Gureyev et al., 2014
; Gureyev et al., 2016
; Gureyev et al., 2020
; de Hoog et al., 2014
). The NRU states that, for a broad class of linear imaging systems, the ratio of the squared signal-to-noise ratio (SNR) to the minimal spatially resolvable volume, Δn (where Δ is the spatial resolution and n = 2 for planar imaging systems, n = 3 in 3D imaging, etc.), depends only on the incident photon fluence. In particular, SNR2/Δn remains unchanged after detector pixel binning, linear image filtering or denoising, etc. The invariance of this ratio under linear photon-number-conserving transformations of imaging systems is based on the fact that, at least for Poissonian photon-counting statistics, the ratio SNR2/Δn corresponds to the mean number of detected photons that interacted with a minimal resolvable volume of the sample in the process of imaging. The latter quantity determines the amount of Shannon information that the system is capable of extracting about the imaged sample (Gureyev et al., 2016
). As such, SNR2/Δn is an essential intrinsic characteristic of the imaging system. It is, of course, trivial to decrease the ratio SNR2/Δn by applying an SNR-reducing operation to an image, e.g. adding artificial noise to it. Note, however, that if an image processing operation utilizes a priori information about the imaged sample, it can also potentially increase the ratio SNR2/Δn without violating the NRU. In particular, this can take place in deep machine learning (artificial intelligence) based image denoising methods, such as those employing the highly successful UNet architecture (Ronneberger et al., 2015
; Pakzad et al., 2025
).
It has been argued (see, for example, Gureyev et al., 2017
; Gureyev et al., 2024
) that the key benefit of PBI is its ability to `beneficially violate' the NRU, i.e. increase SNR2/Δn in comparison with an equivalent attenuation-based imaging setup (with zero propagation distance) without increasing the or using a priori information. This can be achieved in PBI either by improving the spatial resolution without lowering the SNR (in the forward process of PBI image formation) or by increasing the SNR without spoiling the spatial resolution [in the combined process of forward PBI imaging and subsequent `phase retrieval' in accordance with Paganin's method (Paganin et al., 2002
)]. Since the ratio SNR2/Δn can be increased this way in PBI, as demonstrated in multiple publications to date (see, for example, Gureyev et al., 2014
; Gureyev et al., 2024
; Kitchen et al., 2017
; Brombal et al., 2018
), this seems to represent a true violation of the NRU. It may appear that additional information is created `out of nothing' in PBI, since the average number of photons per unit volume does not change in the process of forward free-space propagation of the beam. In fact, this phenomenon is rooted in the imaging physics, and more specifically in the process of interaction of coherent radiation with an imaged sample. The NRU implicitly assumes that any information about the imaged sample can only be obtained by directly counting the photons that have been collected at the detector pixels. In that process, each photon carries no more than a single bit of information: essentially either the photon was transmitted through the sample and then detected (1), or it was absorbed or scattered away by the sample and hence not detected (0). However, in phase-contrast imaging, one can also detect and measure the phase shifts that the photons acquire upon coherent scattering by the atoms of the imaged sample. This information is made detectable in PBI by means of free-space propagation which allows the transmitted photons to interfere and thus reveal the phase information in the registered pattern of photon fluence. This detected signal is no longer binary, as in attenuation-based imaging, but can take any real value in the form of the measured phase shift. This means that the assumptions used in NRU no longer hold in phase-contrast imaging. The difference between the absorption and phase imaging here resembles the difference between the classical and quantum computing, with the former operating with binary bits and the latter utilizing qubits that can be in a weighted superposition of two states (Nielsen & Chuang, 2010
).
An even more important factor leading to the `beneficial violation' of NRU and the `unreasonable' effectiveness of PBI (Gureyev et al., 2017
) is the difference between the `coupling strengths' of the scattering channels corresponding to absorption and phase shifts in the interaction of incident photons with the imaged sample. In the case of imaging soft biological tissues with hard X-rays, the phase `coupling' is about three orders of magnitude stronger than the absorption `coupling' (Momose, 2020
). Therefore, a phase-contrast method, like PBI, can overcome the limits set in NRU by utilizing stronger scattering channels that are not taken into account in the theory underpinning NRU. The gain in SNR in PBI, in comparison with the equivalent absorption-based imaging at the same incident photon fluence and spatial resolution, was discussed in the first part of this paper and is further studied below in the case of three-dimensional (3D) imaging. It turns out that this gain is determined by the ratio of two dimensionless parameters, γ/NF, where γ is essentially the ratio of the strengths of the phase and absorption channels of scattering in the sample and NF is the `minimal Fresnel number' that is equal to Δ2 divided by the X-ray wavelength and by the sample-to-detector propagation distance. In biomedical PBI and propagation-based phase-contrast computed tomography (PB-CT), the ratio γ/NF can be of the order of 10–100. This effectively means that each detected photon in PB-CT can potentially carry hundreds of bits of information about the imaged sample, which appears to be an interesting and non-trivial phenomenon. It also leads to a potential decrease of the radiation dose of the order of 102 to 105 in PB-CT in comparison with conventional attenuation-based CT (Kitchen et al., 2017
).
For the subject of the present paper it is important to consider another, this time `detrimental', violation of NRU which occurs in computed tomography (CT). It is well known that CT is a mathematically `ill-posed' technique (Natterer, 2001
). While `ill-posedness' of a system of equations can indicate non-existence or non-uniqueness of solutions, CT reconstruction actually does deliver a unique solution, i.e. a 3D reconstructed distribution, for any correctly sampled input projection data. However, CT is mathematically unstable, i.e. it amplifies the noise present in the input data. This can be understood by looking at the eigenvalues of the (forward) X-ray transform (consisting of X-ray projections collected at a set of view angles over 180° or 360°), which decrease with the radial frequency of the corresponding eigenfunctions (Natterer, 2001
). As a consequence, the X-ray transform progressively suppresses higher-order spatial frequencies in the input data. Accordingly, the eigenvalues of the inverse X-ray transform, i.e. the CT reconstruction operator, increase in magnitude with the radial frequency of its eigenfunctions. Therefore, CT reconstruction amplifies higher spatial frequencies more strongly than the lower ones. This can be also observed through the presence of the `ramp' filter in the continuous CT reconstruction (see the next section for details). As a typical signal in CT decreases more rapidly with increasing spatial frequencies compared with a typical noise profile, this means that CT reconstruction usually amplifies the noise more than the signal, thus reducing the SNR. This leads to a detrimental violation of the NRU via the fact that SNR2 in 3D CT is proportional to Δ4, rather than Δ3, as it should be in NRU (Howells et al., 2009
; Nesterets & Gureyev, 2014
). Nevertheless, when CT is combined with PBI, into the method known as PB-CT, it becomes significantly more efficient in terms of the increased SNR2/Δn ratio in the reconstructed data, in comparison with conventional attenuation-based CT. This was previously demonstrated both theoretically and experimentally (Nesterets & Gureyev, 2014
; Kitchen et al., 2017
; Brombal, 2020
). Mathematically, the effectiveness of combining PBI with CT comes from the suppression of the noise-amplifying CT ramp filter by the noise-suppressing filter function of Paganin's method (Bronnikov, 1999
; Paganin et al., 2002
; Gureyev et al., 2006
). As a result, SNR2 in PB-CT becomes almost independent of the spatial resolution (Nesterets & Gureyev, 2014
). This subject is also central to the present paper.
In order to make this second part of our paper largely self-contained, we have recalled the necessary basic facts about PB-CT in Section 2
below. In Section 3
, we consider the question of quantitative assessment of performance of imaging systems in general and phase-contrast X-ray imaging of biological samples in particular. Several of our previous studies, including the first part of this paper (Gureyev et al., 2024
), utilized the `intrinsic imaging quality' characteristic, QS, the square of which is equal to the ratio of the quantity SNR2/Δn discussed above to the corresponding incident n-dimensional photon fluence, Iin, QS 2 = . The NRU is essentially equivalent to the statements that QS cannot be larger than unity and that QS is conserved in photon-number-conserving linear transformations of the imaging system. The intrinsic imaging quality characteristic corresponds to the amount of Shannon information that the imaging system is capable of extracting per single incident photon (Gureyev et al., 2016
). It is typically measured in `flat fields', i.e. in images or parts of images obtained with no sample interaction with the beam, in which case QS depends only on the properties of the imaging system alone. If QS is measured in a uniform area of the image in the `shadow' of the sample, X-ray absorption in the sample does affect QS by reducing the SNR. Although this effect is generally trivial, it cannot be ignored completely. More importantly, as in practice, especially in biomedical imaging, the image contrast generated by the sample and the radiation dose delivered to the sample in the process of imaging are both critically important, it appears useful to consider a modified form of the intrinsic imaging characteristic that would explicitly include these two parameters into an image quality metric. We introduce such a metric, which we have tentatively named the `biomedical X-ray imaging quality characteristic', QC, in Section 4
. In Section 5
, we discuss the theoretical gain that PBI and PB-CT can deliver for QS and QC, and then evaluate the gain in imaging quality of experimental images of breast tissue samples obtained using PB-CT in Section 6
. This is followed by concluding remarks in Section 7
.
A few technical derivations, the results of which are referenced in the paper, can be found in Appendices A
, B
and C
. The main text of the paper also contains a relatively large number of mathematical expressions, which are essential to the considered problems and their solutions. For each mathematical expression or equation, we discuss both their physical meaning and their role in the physical picture and practice of PB-CT imaging in general and in biomedical imaging applications in particular.
2. Basic theory of propagation-based phase-contrast tomography
Let ,
=
, be the (slowly varying envelope of the) complex amplitude of a scalar monochromatic paraxial X-ray wave (beam) with wavenumber k = 2π/λ, where λ is the wavelength,
=
are Cartesian coordinates in three-dimensional (3D) space, z is the beam propagation direction and
= (x,y) are the coordinates in transverse planes. The square modulus of the complex amplitude, I(r), will be identified with the photon fluence and expressed in photon number per unit area. We assume that the X-ray wave near the vicinity of the z-origin of coordinates is a plane wave:
=
. The process of transmission of this incident plane wave through a thin object with complex refractive index n(r; λ) = 1 − δ(r; λ) + iβ(r; λ), located in the vicinity of the origin of coordinates immediately upstream of the `object plane', z = 0, can be described by the Beer–Lambert law,
where is the transmitted complex amplitude in the `object' plane, z = 0, δ1(x) is the one-dimensional Dirac delta function and P is the X-ray transform (projection) operator (Natterer, 2001
) which integrates its argument, f, along straight lines parallel to the beam direction. The primed Cartesian coordinates, r′ = (x′, y, z′), in equation (1b)
are associated with the sample rotated by angle θ around the y axis with respect to the fixed coordinates r, which are associated with the fixed incident beam direction and the detector plane. Note that, starting from equations (1a)–(1b)
, we omit the dependence of any quantities on the wavelength for brevity. The propagation of the transmitted complex amplitude from the object plane z = 0 to the image (detector) plane z = R is described by the Fresnel diffraction integral (Mandel & Wolf, 1995
),
The assumption about the `thinness' of the object means that, for any rotational angle θ, δ(r′) = β(r′) = 0 outside the narrow slab −L < z < 0, where 0 < L << R.
A complex amplitude is called monomorphous (or homogeneous) within some region of space (Paganin et al., 2002
; Paganin et al., 2004
) if the proportionality coefficient γ between its phase and the logarithm of its modulus is independent of the position in that region. Monomorphous amplitudes arise, for example, in the object plane after transmission of an incident plane monochromatic X-ray wave through an object having the same chemical composition everywhere, possibly with a spatially varying density. Indeed, for such objects the ratio of the real decrement to the imaginary part of is the same at any point inside the sample: δ(r)/β(r) = γ = const (Paganin et al., 2004
). Therefore, the phase, =
, and the square modulus,
=
, of a wave transmitted through such an object satisfy a simple relationship,
If the object-plane intensity also varies sufficiently slowly, such that <<
, where a ≡ [γRλ/(4π)]1/2 and
=
is the transverse 2D Laplacian, the square modulus of equation (2)
can be accurately approximated by the homogeneous Transport of Intensity equation (TIE-Hom) (Paganin et al., 2002
),
Note that equation (4)
describes the propagation of the transmitted from the object plane to the image plane without any explicit references to the phase distribution, i.e. without reference to the real part of the of the imaged object (`the sample'). Of course, in the case of non-trivial propagation distances R, this can only be possible due to the `homogeneity assumption', δ(r) = γβ(r), about the sample, which allows for the phase information to be included in equation (4)
implicitly. Interestingly, equation (4)
has the structure of a finite-difference approximation to the diffusion equation, with −a2/R playing the role of a negative Under this view, the use of the assumption implies that it is now transverse intensity gradients (rather than phase gradients) that induce propagation-based transverse flow. The notion that intensity gradients drive flow is suggestive of diffusive energy transport, with the negative corresponding to sharpening rather than blur (cf. Gureyev et al., 2004
).
It has been shown (Thompson et al., 2019
; Gureyev et al., 2022
) that equation (4)
can be re-written as
where
is the `PBI contrast function' and
=
is the 3D Laplacian. In particular, equation (5)
effectively states that the TIE-Hom operator commutes with the projection operator. As a result, the `propagated' at any rotational position of the sample can be obtained by applying the 3D TIE-Hom operator first, followed by the X-ray projection at the corresponding rotational position of the sample.
If a full set of propagated transmission images at different view angles in a 180° range, , is available, equation (5)
can be inverted to reconstruct the 3D distribution of the n(r) = 1 + (i − γ)β(r) in a homogeneous sample (Gureyev et al., 2006
; Thompson et al., 2019
),
where P−1 is the inverse X-ray transform (i.e. CT reconstruction) operator (Natterer, 2001
), while and (1 − a2∇2)−1 are the inverse TIE-Hom operators in 2D and 3D, respectively. Note that at R = 0 we have a = 0 and also
=
Therefore, at R = 0, equation (6)
transforms into the conventional CT reconstruction formula, =
(Natterer, 2001
). Of course, in practice, the infinite set of projections, , is substituted by a properly sampled finite discrete set of measured fluence values,
l =
m =
n =
, collected at detector pixel positions (xl, ym) and view angles θn (Natterer, 2001
). A properly sampled set of such projection values allows one to reconstruct a unique distribution of discretely sampled values of β(rlmn) at LMN voxels using a discretized version of equation (6)
.
The inverse X-ray transform operator P−1 can be expressed either numerically, e.g. using a gridding-based implementation of the projection-slice theorem (Natterer, 2001
; Gureyev et al., 2022
), or analytically, e.g. with the help of the filtered back-projection (FBP) form of the inverse X-ray transform (Natterer, 2001
),
where =
is the 2D Fourier transform and
=
is dual to
= (x,y).
The 2D and 3D TIE-Hom retrieval operators, and (1 − a2∇2)−1, have simple analytical representations in Fourier space, where they can be expressed as multiplication by the functions
and 1/(1 + 4π2a2ρ2), respectively [here and below we use the notation
=
and
=
]. In real space, the 2D operator
can be represented as a convolution with the function
, where K0 is the zero-order modified Bessel function of the second kind (Nesterets & Gureyev, 2014
). The 3D operator (1 − a2∇2)−1 can be expressed as a convolution with the Yukawa (screened Coulomb) potential =
(Sakurai, 1967
).
Equations (5)
and (6)
provide one-to-one forward and inverse mappings between (i) the 3D distributions of the imaginary part of in a homogeneous sample and (ii) the PBI contrast function collected in an imaging experiment with the sample scanned over 180° of rotation around a fixed axis perpendicular to the X-ray beam direction and a detector located in the so-called `near-Fresnel' region, where <<
. The key question that we want to address in the present paper is that of the `gain' in the SNR-to-spatial-resolution ratio in the PB-CT reconstruction of β(r) from the propagated projection images
, compared with the conventional CT reconstruction at R = 0, as a function of propagation distance R, the X-ray wavelength, the spatial resolution of the imaging system and the complex refractive index of the imaged sample.
3. Quantitative performance measures of an imaging system
In the case of images collected with a position-sensitive detector, we shall distinguish between a stochastic registered distribution of the , and its mean value,
, which can be obtained by point-wise statistical averaging of an ensemble of equivalent images collected repeatedly under the same conditions. More generally, we consider random functions I(r) in n-dimensional space, where, for example, n = 2 in the case of conventional planar images and n = 3 in the case of CT-reconstructed volumes.
To avoid extensive repetitions, we shall refer to some definitions included in the first part of the present paper (Gureyev et al., 2024
). The SNR was defined in equation (5)
of Gureyev et al. (2024
) as the ratio of the mean of a random function I at a point r to its standard deviation, =
. It was mentioned that, when the spatial ergodicity condition is satisfied, the signal and noise variance can be evaluated via spatial integrals over a flat vicinity of a given point.
The spatial resolution can be defined in terms of two different definitions of the width of the point-spread function (PSF) of the imaging system, P(r):
(a) Δ[P], such that =
/
, defined in equation (6) of Gureyev et al. (2024
), and
(b) =
, where
, as defined in equation (10) of Gureyev et al. (2024
). It was mentioned by Gureyev et al. (2024
) that Δ[P] and produce the same or very similar values for the width of many types of PSFs.
An important concept that was already discussed by Gureyev et al. (2024
) and will be used here is that of the noise-resolution duality. Let be a convolution of a random function I(r) with a filter function O(r). It was shown by Gureyev et al. (2016
) that if O(r) is much broader than P(r), but varies much faster than and
, then the following ratio remains unchanged after the filtering,
In the case of 2D imaging with Poisson photon statistics, SNR2 is equal to the number of registered photons per detector pixel. Therefore, in this case, the `noise-resolution duality', equation (8)
, is just a restatement of a simple fact that larger `voxels' with the volume , effectively created as a result of image filtering, contain more registered photons compared with the original `voxels' with the volume
, leading to the proportionally larger SNR2.
Now we can introduce the notion of the (modified) intrinsic imaging quality characteristic, , of an imaging system, which is defined as
where SNR2[I] represents an average of SNR2[I](r) over a suitably representative flat area of an n-dimensional image, is the width of the PSF, as defined above, and
denotes the mean incident defined as the number of incident photons per n-dimensional volume. The subscript `S' in
serves as a reminder that this quantity is meant to be a characteristic of the imaging system rather than a characteristic of an image. Note that
is a positive dimensionless quantity. Compared with the definition of intrinsic imaging quality characteristic, QS, used previously (for example, by Gureyev et al., 2014
, Gureyev et al., 2016
, Gureyev et al., 2018
, Gureyev et al., 2020
, Gureyev et al., 2024
), equation (9)
utilizes the spatial resolution instead of Δ[P]. Although, as mentioned above, the two measures of the spatial resolution are identical or very similar for many popular functions P,
generally has more favourable formal mathematical properties compared with QS.
Note that =
, where
is the mean number of photons incident on a single voxel during the imaging process and
is the n-dimensional volume of the voxel. When the incident fluence is Poissonian, we have
=
, and hence
where p = is the ratio of the linear size of the voxel to the spatial resolution, and q =
= DQE(0) is the at zero frequency (Bezak et al., 2021
). It has also been shown (Gureyev et al., 2016
) that corresponds to the Shannon information that the imaging system is able to extract from each incident photon. Therefore, imaging systems with higher
allow one to obtain better images that contain more information about the samples, compared with systems with lower
.
It follows from equation (8)
that remains unchanged in any post-detection linear filtering operation preserving the integral value of I, such as, for example, detector pixel binning, low-pass filtering or `phase retrieval' using the inverse TIE-Hom operators
and (1 − a2∇2)−1. In the first part of this paper, we discussed in detail that, although the forward TIE-Hom operator in equation (4)
is also linear and accurately describes the free-space propagation of the slowly varying mean fluence distribution, , from the object plane z = 0 to the detector plane z = R, it does not correctly describe the corresponding propagation of the noise term,
(we omit the projection angle θ here for brevity). According to the conservation of in paraxial free-space propagation, the noise variance in flat areas of images is the same at z = 0 and z = R, regardless of the propagation distance R (at least within the near-Fresnel region). On the other hand, the spatial resolution in TIE-Hom imaging improves after free-space propagation, so that
in the 2D case (Gureyev et al., 2017
). The increase of =
is then enabled by the improvement in the spatial resolution:
]
, where
is the `minimal' Fresnel number and
=
. As mentioned above,
remains unchanged after the TIE-Hom retrieval. Therefore, one obtains
=
, where
is the intrinsic imaging quality after the free-space propagation followed by the TIE-Hom retrieval (Gureyev et al., 2017
). Since γ/NF can be much larger than 1 in hard X-ray imaging of biological samples, can be much larger than
, and
=
can be significantly larger than 1 (Gureyev et al., 2017
; Gureyev et al., 2024
). Note that such large values of only occur in imaging of `monomorphous' objects to which equation (4)
applies. We will investigate the corresponding behaviour of in 3D, in the case of PB-CT, in Section 5
below.
4. Characteristics of image quality dependent on the imaged sample
We have so far considered different performance metrics of imaging systems (SNR, ,
) without making reference to imaged samples. In practice, it is often required to analyse the image quality and optimize the performance of an imaging system for a particular class of samples. In this respect, one such metric that has not been considered above is the image contrast. Indeed, useful image contrast can typically appear only in the presence of an imaged sample, and the strength of contrast inevitably depends on the sample, as well as on the imaging system. The contrast is often defined as C =
, where
and
are the mean pixel values in two adjacent areas of the image. We will use a slightly modified version of this definition, Cm =
, which has similar properties to C. Indeed, both C and Cm can take values between 0 and 1; both reach the maximum value of 1 when one of the two intensities is equal to zero, and the minimum value of 0 when the two intensities are equal. Moreover, in the case where
=
, we get C = 1/3 and Cm = 1/2, with the latter result possibly being more natural on an intuitive level. We can then define the contrast to noise ratio as
where SNR = is measured in the (`background') area with higher pixel values,
=
. The CNR is an important characteristic of an image, and the goal is typically to maximize it. Alternatively, when the CNR is close to zero or is below the noise level, the image is likely to be considered useless.
One type of contrast which is often used in X-ray imaging is that between an absorbing sample and the flat-field background corresponding to the incident illumination. In this case we have I2 = Iin and I1 = , where μ = 4πβ/λ is the of the sample, the integral is along the rays passing through the sample, and the relevant mean values can be calculated over the transverse coordinates in the region of interest (ROI) in the image. When the absorption is weak, such that
<< 1 and hence
, the contrast is approximately equal to Cm
.
In the case of PBI and PB-CT, an important model case of image contrast is represented by an edge of a slab of weakly absorbing (`transparent') material with the real decrement of refractive index, δ. In the case of a parallel monochromatic X-ray beam, the corresponding phase shift is =
. The `edge contrast', representing a simple case of propagation-based phase contrast, is defined as the difference between the maximum and the minimum of pixel values in the first Fresnel fringe generated by the edge at the sample-to-detector distance R, divided by the sum of the same two extreme intensities. It has been shown (Gureyev et al., 2008
) that such phase edge contrast can be expressed as
=
. The latter expression is valid under the condition
=
>>
, in the so-called near-Fresnel region.
Another important characteristic of an image is the radiation dose that has been delivered to the sample in the process of acquiring the image. The absorbed dose, Dab, is the energy per unit mass absorbed by the sample, which is expressed in Gy (= J kg−1). The importance of the radiation dose can be due to the carcinogenic effect of ionizing radiation in medical imaging, or due to the damage by electrons used for imaging in high-resolution electron microscopy, or for other reasons relevant to a particular imaging context. The absorbed dose is closely related to the kerma, K, which is the energy transferred from photons to kinetic energy of electrons in the unit mass of matter in which the kerma is measured,
where the coefficient g is the fraction of energy lost to Bremsstrahlung and other radiative processes (Bezak et al., 2021
). Both the absorbed dose and the kerma can be factorized into the product of the incident photon fluence, the mass absorption coefficient, μ/ρ, where ρ is the sample density, and the average energy absorbed, , or transferred,
, respectively, at each interaction:
=
,
=
and K =
,
=
.
In the context of X-ray imaging of the breast, which is the most relevant to the experimental results analysed below, the standard measure of the dose is the mean glandular dose (MGD), , which is calculated from the measurements of the entrance air kerma,
=
, according to the formula
where is a dimensionless conversion factor which depends on the X-ray energy, the breast thickness and composition (such as glandularity) (Bezak et al., 2021
). This conversion factor is usually obtained by Monte Carlo simulations with breast-equivalent phantoms (Johns & Yaffe, 1985
; Nesterets et al., 2015
). The maximum permissible MGD for a mammographic image is currently around 1 mGy, with some variation in the standards between different countries (Liu et al., 2022
).
In order to include the effect of the sample on imaging quality, we first multiply the intrinsic imaging quality characteristic by the contrast, Cm, which can be calculated for a feature of interest in the image. Then, in order to take into account the detrimental effect of the radiation dose, in accordance with the `linear no-threshold model' (Bezak et al., 2021
), we also divide the new characteristic by the absorbed dose or by the MGD, depending on the context. However, the resultant quantity would depend on the measurement units, such as kg J−1 or g calorie−1. In order to avoid this cumbersome dimensionality and make the new imaging quality metric dimensionless, we can follow the general idea behind the Hounsfield units and normalize our metric accordingly. It appears that, in the case of biomedical X-ray imaging, a natural normalization can be achieved with respect to X-ray absorption in air,
or, in the special case of breast imaging,
Here we used the fact that, at X-ray energies typical for breast imaging (∼20–40 keV), the coefficient gair is much smaller than unity (see, for example, Kato, 2014
). We propose to call the metric `biomedical X-ray imaging quality'. Note that the word `intrinsic' is no longer relevant here, since, unlike
, the new metric
substantially depends both on the imaging system and on the imaged sample. We used the normalization coefficient
related to air in equations (14) and (14a)
, and not, for example, to water as in the case of Hounsfield units, because the sample is assumed to be imaged in air. In the case of breast imaging, the conversion factor for MGD is also defined with respect to the air kerma. If the sample is absent and there is only air, then equations (14)
and (14a) both convert simply to =
. The latter quantity is likely to be equal to zero, since in a well designed imaging system the contrast should be close to zero in the absence of a sample.
Let us consider first the case of 2D imaging,
or, in the case of breast imaging,
The mean absorbed dose (MAD) or MGD are typically measured with the help of ionization chambers or dosimeters during an experiment, while CNR and can be measured later in the collected images. In order to calculate the imaging quality metric
, one also needs to know the normalization constant
. In the case of monochromatic X-rays, we have
=
, where μen/ρ is the mass energy-absorption coefficient and Eph = hc/λ is the energy of a photon (Hubbell & Seltzer, 2004
). Calculating the photon energy is trivial and the mass energy-absorption coefficient of air at a given X-ray energy can be found, for example, in NIST databases (Higgins et al., 1992
; Hubbell & Seltzer, 2004
) and elsewhere. For example, at 32 keV, hc/λ ≅ 6.626 × 10−34 [kg m2 s−1] × 2.998 × 108 [m s−1]/3.875 × 10−11 [m] ≅ 5.127 × 10−15J, μen,air/ρair ≅ 0.01366 m2 kg−1 for dry air at sea level (Kato, 2014
), and hence
7.003 × 10−17 m4 s−2. In a typical medical X-ray imaging setup,
and
can be of the order of 10−3 Gy and
can be of the order of 100 µm. In such cases, the dimensionless ratios
and
will be of the order of 10−6. If the contrast is not too small, one will need to have SNR2 of the order of 106, or about 106 detected photons per pixel, for
to approach unity. Of course, any increase in SNR via an increased incident fluence will be associated with a simultaneous increase in the dose, so improving the biomedical X-ray imaging quality characteristic of an imaging setup is not a trivial task.
While equation (9)
for has the same mathematical form for any n in n-dimensional imaging, equations (15) and (15a)
do not retain the same form for n ≠ 2, due to the intrinsically 2D nature of the definition of the dose. Let us see how equation (15)
changes in a 3D-imaging case, such as CT imaging. Using equation (14)
, we obtain
In order to correctly re-write the denominator of the right-hand side of equation (16)
in terms of the 3D dose, let us first express the sample volume as V = ΩL, where Ω is the illuminated `front' surface area and L ≡ V/Ω is the effective depth of the sample. In the situation most frequently encountered in CT, the reconstructed sample volume is cylindrical: V = , where RC is the radius of the cylinder and H is its height. The air kerma K =
is typically measured with respect to the fluence
corresponding to a flat entrance surface with the area Ω = 2RCH and hence L = V/Ω = (π/2)RC in this case. Let N3D = N2DMa be the total number of photons, or, more precisely, the number of noise equivalent quanta, NEQ (Bezak et al., 2021
), collected in the detector area corresponding to Ω during a CT scan. Here N2D is the average total number of NEQs detected in each 2D projection and Ma is the number of angles at which the CT projections have been acquired. Then, assuming that the measurements are performed in a uniform area of the images, =
=
=
and hence the denominator of equation (16)
is equal to =
=
, where
=
=
is the total MAD accumulated during the scan. Equation (16)
can now be written as
Similarly, in the case of breast imaging, we obtain
Note that the appearance of the factor L in equations (17) and (17a)
is the consequence of the fact that MAD and MGD are defined with respect to the 2D rather than the 3D which appears in the 3D version of . In a limit case, when the number of projections is equal to 1, one has Dab,3D = Dab,2D and Dg,3D = Dg,2D, the voxel volume is
=
,
=
, and hence equations (17)
and (17a) morph into equations (15)
and (15a), respectively. Note also that, as can be seen in the latter example, the resolution volume in equations (17)
and (17a) may not be equal, in general, to , where
is the resolution area used in equations (15)
and (15a). Therefore, and
are just shorthand notations for
and
, where
and
are the spatial resolutions in the 2D and 3D imaging, respectively.
It is important to remember that both and
have been defined as imaging quality measures per single incident photon. In some cases, it can be useful to evaluate the quality of a particular 2D or 3D image, obtained with a given incident fluence or a given dose. Natural metrics for this purpose are represented by the product of the above metrics, defined per single photon, and the total number of the incident photons, Ntot, that was used to create the image, i.e.
and
. Note that
=
, where V is the n-dimensional volume of the imaged sample. Therefore,
=
, where Mvox is the total number of n-dimensional voxels, and hence
where is the squared `total' SNR of the image, obtained by summing all the squared SNRs measured in individual voxels. The latter quantity is equal to the total number of noise-equivalent quanta, NEQtot, in the image (Bezak et al., 2021
). Therefore, the intrinsic image quality, , provides a direct quantitative measure of the (Shannon) information content of the image. For comparison, recall that
is closely related to Shannon's channel information capacity of the imaging system (Gureyev et al., 2016
), which is the amount of information per single photon.
Similarly,
i.e. the square of the intrinsic quality of an image of a sample is equal to the total number of NEQs in the image, multiplied by the square of the contrast and by the ratio of the mass energy-absorption coefficients of air and the tissue. Note that, when using equation (19)
, the contrast should be approximately uniform in the whole image. Alternatively, an appropriate kind of average contrast across the image can be used, or the metric from equation (19)
could be used for a particular part of the image, etc.
5. Noise and spatial resolution in PB-CT
We proceed with the evaluation of the ratio in the distribution of β(r) reconstructed from PB-CT scans collected in an image plane z = R, in comparison with the same ratio obtained in conventional CT scans collected in the object plane z = 0. The (average) `gain factor' in the intrinsic imaging quality of the PB-CT imaging, compared with the conventional CT, can be defined as follows (Nesterets & Gureyev, 2014
; Gureyev et al., 2017
; Gureyev et al., 2024
),
where the subindexes `0' and `retr' correspond to quantities obtained from the data collected by conventional CT at z = 0 and by PB-CT with TIE-Hom retrieval from the data collected at z = R, respectively, with and
denoting the mean incident fluence in the object and the image planes.
If we consider the gain factor at a certain fixed level of the incident X-ray fluence, then, as discussed above, the mean fluences at z = 0 and z = R are going to be equal, due to the conservation of in free-space propagation (ignoring the absorption in air, for simplicity). Moreover, if the TIE-Hom retrieval is performed with a = [δRλ/(4πβ)]1/2, then, according to equation (6)
, the same distribution of is obtained after the reconstruction, regardless of the distance R. Indeed, the mean values of the input to equation (6)
, , are the same in flat areas of images at different R within the near-Fresnel region, due to the conservation of in free-space propagation. Furthermore, the inverse X-ray transform, P−1, does not depend on R, while the inverse TIE-Hom operators,
and (1 − a2∇2)−1, both preserve the mean value of their input for any R, since they can be represented as convolutions with filter functions
and T3,inv(r) which both have unit integrals. Therefore, both the signals and the spatial resolutions are the same in the case of conventional and PB-CT reconstructions, regardless of the value of R, i.e.
=
and
=
. In this case, equation (20)
reduces to a simpler expression, =
=
(Nesterets & Gureyev, 2014
).
Note that although equation (20)
is expressed in terms of the intrinsic imaging quality characteristic, the gain factor still depends on the type of imaged samples. Indeed, the ratio δ/β used in the TIE-Hom retrieval step is specific to a particular class of samples, and only for such samples can the gain in SNR be achieved without a loss of spatial resolution after the TIE-Hom retrieval of the images collected at z = R using the specified parameter a = [δRλ/(4πβ)]1/2. Therefore, it is logical to take other relevant sample-specific factors, such as contrast and dose, into account and consider the gain coefficient for the biomedical X-ray imaging quality,
. Since the distribution of the mean values
, obtained as a result of a PB-CT reconstruction, is invariant with respect to the propagation distance R, the contrast Cm =
in the reconstructed images, e.g. between the mean values of β in adjacent areas of adipose and glandular breast tissues, is independent of R. Therefore, the improvement in the biomedical X-ray imaging quality characteristic
upon free-space propagation followed by the TIE-Hom retrieval, compared with the conventional CT at R = 0, is the same as in equation (20)
,
In the case of breast imaging, the MADs in equation (21)
can be replaced by the corresponding MGDs.
Furthermore, as mentioned in the first part of this paper (Gureyev et al., 2024
), when imaging a class of samples with the same δ/β ratio, the gain factor is independent of the parameter a = [δRλ/(4πβ)]1/2 used in the TIE-Hom retrieval, as long as this procedure satisfies the conditions required for validity of equation (8)
. With smaller values of parameter a, the reconstructed images will have lower SNR, but higher spatial resolution, and the opposite will be true for larger values of a. The invariance of the gain factor with respect to the choice of parameter a is important, because TIE-Hom retrieval in PB-CT imaging is often deliberately performed with smaller values of a in order to retain some image-sharpening effects of coherent free-space propagation (Tavakoli Taba et al., 2019
). Another reason for using a smaller value of a in TIE-Hom retrieval is that the penumbral effect of the source contributes to the spatial resolution and already does some `phase-retrieval blur' (Beltran et al., 2018
).
While both the signal and the spatial resolutions are the same in the case of conventional and PB-CT reconstructions, regardless of the value of R, the noise variance in the reconstructed distribution may vary dramatically as a function of R in PB-CT at a fixed level of the incident fluence (Nesterets & Gureyev, 2014
; Kitchen et al., 2017
; Brombal, 2020
). Let us consider the latter effect quantitatively, on the basis of the noise-resolution duality, equation (8)
. According to the last expression in equation (6)
, conventional CT and PB-CT differ only by the free-space propagation and the application of the 3D TIE-Hom retrieval operator, (1 − a2∇2)−1. As mentioned above, the action of the operator (1 − a2∇2)−1 can be expressed as a convolution with the function T3,inv(r). This fact is used in Appendix A
to derive the following approximation for the 3D gain factor,
where f(A) = , A =
and
is the Fresnel number corresponding to the spatial resolution after the CT reconstruction. Equation (22)
was obtained under the condition that the width of T3,inv is much larger than the detector resolution, hence >>
, A >> 1 and f(A) ≃ 1, implying that G3 > 1. For completeness, we note that G1 = (4/π)1/4(γ/NF)1/4 and G2 = (γ/NF)1/2 (Nesterets & Gureyev, 2014
; Gureyev et al., 2017
; Gureyev et al., 2024
), where
and
is the spatial resolution in the projection images.
While equation (22)
was obtained using the last equality in equation (6)
, a somewhat different approximation for the 3D gain factor was obtained earlier on the basis the first part of equation (6)
(with 2D TIE-Hom retrieval) (Nesterets & Gureyev, 2014
),
A brief derivation of equation (23)
is given in Appendix B
for completeness. This form of the gain factor was previously investigated and discussed (Kitchen et al., 2017
; Brombal et al., 2018
; Brombal, 2020
). We shall see in the next section how well these theoretical estimates of the gain factor agree with our experimental results.
It follows from the analytical expression for the biomedical X-ray imaging quality of CT derived in Appendix C
that
where is the average and MR =
=
is π/2 times the number of spatial resolution units in the radius of the CT reconstruction cylinder. An analytical expression for conventional CT, i.e.
, is obtained from equation (24)
when G3 = 1. The constant factor (12/π2) in equation (24)
corresponds to a particular implementation of the CT reconstruction, namely the filtered back-projection with a nearest neighbour interpolation; it can be slightly different for other implementations of the CT reconstruction (Nesterets & Gureyev, 2014
). As could be expected, the biomedical X-ray imaging quality in equation (24)
is directly proportional to the image contrast, Cm, while being inversely proportional to the square root of the ratio of the MGD to the air kerma, =
. The term
in the numerator of equation (24)
achieves the maximum value of approximately 0.54 at = 2. The factor MR in the denominator of equation (24)
is a `proxy' for the highest radial frequency present in the reconstructed distribution of β, which determines the degree of ill-posedness of the CT reconstruction, as mentioned earlier (Natterer, 2001
). Since in practice the numerical factor MR is likely to be of the order of 102 to 104, and the other dimensionless parameters in equation (24)
are likely to be of the order of unity, the biomedical X-ray imaging quality of CT is likely to be much smaller than unity. Equation (24)
can be used for design and optimization of PB-CT imaging setups, and for comparison with experiments, as demonstrated in the next section.
6. Experimental results
We carried out PB-CT imaging experiments at the Imaging and Medical Beamline (IMBL) of the Australian Synchrotron. We used a quasi-parallel X-ray beam with an energy of 32 keV (λ = 0.3875 Å) and monochromaticity Δλ/λ ≃ 10−3. The source-to-sample and the sample-to-detector distances were R1 = 138 m and R2 = 5 m, respectively. Therefore, the geometric magnification factor was M = (R1 + R2)/R1 ≅ 1.036 and the `effective defocus' distance was R′ = R2/M ≅ 4.83 m. We will explicitly apply below the simple corrections required for this `quasi-spherical' wave geometry (see, for example, Gureyev et al., 2008
) to the results obtained in the previous sections for the parallel-beam geometry. In particular, the expression for the Fresnel number becomes =
, where
is the spatial resolution `projected' from the image plane back to the object plane.
In order to be consistent with the first part of this paper (Gureyev et al., 2024
), we use here a practical measure of the spatial resolution, Res, which is equal to twice the standard deviation of the corresponding PSF. For PSFs with a Gaussian and a rectangular shape, we get, respectively (Gureyev et al., 2016
),
A flat-panel (energy integrating) detector, Xineos 3030HR, with 99 µm × 99 µm pixels, and a photon-counting detector, Eiger2 3MW, with 75 µm × 75 µm pixels, were used in this experiment. The spatial resolutions, as measured in the collected projections, were Res[Xineos] ≅ 150.7 µm and Res[Eiger] ≅ 86.2 µm in the detector plane, which corresponded to values in the object plane:
π1/2150.7 µm/1.036
257.8 µm and
86.2 µm/1.036
83.2 µm. These calculations of
assumed a Gaussian-shape PSF for Xineos and a rectangular-shape PSF for Eiger. In the case of the Xineos detector, previous measurements produced similar values for Res values (Arhatari et al., 2021
). In the case of the Eiger detector, the width of the rectangular PSF was supposed to be close to the pixel size of 75 µm (Dectris AG, 2025
). However, that assumed that the lower energy threshold in the detector was set to a level ensuring complete rejection of split-pixel photon detection events. In the case of 32 keV X-rays, this threshold corresponds to 16 keV. However, in the present experiment we used the lower energy threshold setting of 4 keV in the Eiger detector, in order to increase the detective efficiency, even though it led to slightly worse spatial resolution. The resultant Fresnel numbers were, respectively, NF, M[Xineos] = (257.8 µm)2/(4.83 m × 0.3875 Å) ≅ 355.1 and NF, M[Eiger] = (83.2 µm)2/(4.83 m × 0.3875 Å) ≅ 37.0.
Two mastectomy samples from different donors were imaged in accordance with the Human Ethics Certificate of Approval and with written consent from each donor. The known average characteristics of adipose and glandular breast tissues for 32 keV X-rays implied that the ratio of the differences in real decrement to the differences in imaginary part of the refractive index of these two types of tissue was γ ≅ 869.4 (TS-Imaging, 2025
). Each sample was placed for imaging in a cylindrical thin-walled plastic container with a diameter of 11 cm and scanned in the PB-CT setup with the parameters listed above. Each scan contained 600 equispaced projections over 180° rotation. By measuring the air kerma in the sample plane and using Monte Carlo simulations with breast-equivalent phantoms of size corresponding to the imaged mastectomies, the incident X-ray flux and the rotation speed were adjusted to achieve the MGD of 4 mGy in each scan (Nesterets et al., 2015
). Examples of projection images of the two samples with and without the TIE-Hom retrieval in accordance with equation (4)
, =
, a2 = γR′λ/(4π), are shown in Fig. 1
. Table 1
contains the results of measurements of the gain factor in the individual projections. According to the measurements, the average experimental gain factor in 2D projections was approximately 1.9 for Xineos and 4.2 for Eiger. The theoretical gain values, G2 = (γ/NF)1/2, were respectively G2[Xineos] ≅ (869.4/355.1)1/2 ≅ 1.6 and G2[Eiger] ≅ (869.4/37.0)1/2 ≅ 4.8. The fact that these theoretical values were not exactly equal to the measured ones was due primarily to the uncertainty about the shape of the PSFs of the detectors, especially in the case of the Xineos detector. Using the NRU, equation (8)
, with the measured value
83.2 µm and the ratio of average measured SNRs for Xineos and Eiger from Table 1
, the value of in the object plane for the Xineos detector can be estimated as
83.2 µm × (85.6/33.6)
212.0 µm. The latter value is in between the values of
that are obtained for a rectangular and a Gaussian shape PSFs, given the measured value of Res[Xineos] ≅ 150.7 µm/1.036 ≅ 145.5 µm in the object plane [see equations (25)
and (26)
]. We will use the value of
212.0 µm in the calculations below, which corresponds to NF,M[Xineos] = (212 µm)2/(4.83 m × 0.3875 Å) ≅ 240.1 and G2[Xineos] ≅ (869.4/240.1)1/2 ≅ 1.9.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| Figure 1 PB-CT projections, normalized by flat-field images, of the two mastectomy samples collected at MGD ≃ 4 µGy, E = 32 keV and 5 m sample-to-detector distance, using the Eiger detector. (a) Sample 1, no TIE-Hom retrieval. (b) Sample 1, with TIE-Hom retrieval. (c) Sample 2, no TIE-Hom retrieval. (d) Sample 2, with TIE-Hom retrieval. Square inserts show magnified areas with fine details. |
The intrinsic imaging quality corresponding to the above imaging setup with the Xineos and Eiger detector can be evaluated according to equation (9)
. Using the measured air kerma, we can calculate the incident fluence in the object plane: Iin,2D = K3D,air/(MaRtr,air) ≅ 8 × 10−3Gy/(600 × 7 × 10−17Gy m2) ≅ 0.19 µm−2. Substituting the measured values of SNR and for Xineos and Eiger, we obtain
85.6/[(0.19 µm−2)1/2 × 212.0 µm]
0.93 and
33.6/[(0.19 µm−2)1/2 × 83.2 µm]
0.93. With the TIE-Hom retrieval, the intrinsic imaging quality is increased to
0.93 × 1.9
1.77 and
0.93 × 4.2
3.91. We see, in particular, that the latter values are larger than unity, which would be impossible in conventional attenuation-based imaging.
We performed the reconstruction of `coronal' slices from CT scans collected using two different detectors, as described above. Examples of the reconstructed coronal slices of the two mastectomy samples are shown in Fig. 2
. The results of measurements of the SNR in the reconstructed slices and the 3D gain factors are presented in Table 2
. The theoretical reasons for the measured values of SNR[β] in Table 2
being much smaller than the SNR[I] values in Table 1
are discussed in Appendix C
. The results presented in the last two columns of Table 2
show that the SNR gain factor, compared with absorption-based imaging at the same dose and spatial resolution, was approximately 4.3 for the Xineos detector and 10.9 for the Eiger detector. The theoretical 3D gain factors in the present experiment, calculated in accordance with equations (22)
and (23)
, were, respectively, 3.2 and 4.9 for the Xineos detector and 9.5 and 11.6 for the Eiger detector. As can be seen here, the experimentally measured values of the 3D gain factor lay in between the theoretical estimates given by equations (22)
and (23)
for each of the two detectors.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| Figure 2 PB-CT reconstructions (coronal slices) of the two mastectomy samples from a CT scan at MGD = 4 mGy, E = 32 keV and 5 m sample-to-detector distance, using the Eiger detector. (a) Sample 1, no TIE-Hom retrieval. (b) Sample 1, with TIE-Hom retrieval. (c) Sample 2, no TIE-Hom retrieval. (d) Sample 2, with TIE-Hom retrieval. These images contain the scaled distributions β(r) × 1011. |
While it was possible to perform the measurements of SNR[I0] in the flat-field projection images, where the sample did not affect the measurements, a similar sample-independent measurement was impossible in 3D in the case of SNR[β0]. Indeed, the mean reconstructed value (i.e. `the signal') in CT slices outside of the sample area is equal to zero in accordance with the formalism of CT reconstruction, since = 0 in the absence of X-ray attenuation. Therefore, it was impossible to calculate sample-independent values of SNR and
from the experimental data in this case. As measurements of SNR[β0] presented in Table 2
already depend on the X-ray absorption in the sample, it makes more sense to instead calculate the 3D biomedical X-ray imaging quality characteristic in this case. The contrast between the mean values of β in the glandular and adipose tissue in the reconstructed slices was Cm[β] ≃ (1.0 × 10−10 − 7.3 × 10−11)/1.0 × 10−10 ≃ 0.27. The biomedical X-ray imaging quality characteristic can be evaluated according to equation (17a)
,
, taking into account the MGD,
= 4 mGy, the radius of the cylinder containing the sample, R = 5.5 cm = 2L/π, and the value of
7.003 × 10 −17 m4 s−2 calculated above. Substituting all the relevant parameters into equation (17a)
, we obtained
9.7 × 10−3 and
7.2 × 10−3, without the TIE-Hom retrieval, and
4.1 × 10−2 and
7.8 × 10−2, with the TIE-Hom retrieval. We see that the use of PB-CT substantially increased the biomedical X-ray imaging quality characteristic, in proportion to the 3D gain coefficients G3 measured earlier (4.3 for the Xineos detector and 10.9 for the Eiger detector).
As explained above, the key reason for the values of being much smaller than unity is the mathematical ill-posedness of the CT reconstruction. If we substitute the relevant parameters of the present experiment, e.g. in the case of the Eiger detector, into equation (24)
with G3 = 1, we obtain
9.5 × 10−3 (see Appendix C
for details). This theoretical value is close to the value
7.2 × 10−3 measured in the experiment, confirming that in the present configuration
should indeed be expected to be much smaller than unity. As explained after equation (24)
and in Appendix C
, the small value of is determined by the degree of ill-posedness of the CT reconstruction in the present experimental configuration. The ill-posedness of the CT reconstruction is also behind the fact that, while the
in 2D projections was found to be the same (equal to 0.93) for the Xineos and Eiger detectors,
was slightly higher for the Xineos detector (9.7 × 10−3) than for the Eiger detector (7.2 × 10−3). The reason for this `discrepancy' is that the squared ratio of the averaged SNR in the first and second columns of Table 2
, SNR2[β0, Xineos]/SNR2[β0, Eiger], was closer to the ratio of the fourth powers of the corresponding spatial resolutions, , as expected in CT (Howells et al., 2009
; Nesterets & Gureyev, 2014
), rather than to the ratio of the third powers of the resolutions, , as expected in the NRU and in the definition of the biomedical X-ray imaging quality.
7. Conclusions
We have analysed the performance of PB-CT in a general theoretical context using the NRU as the theoretical basis and the intrinsic imaging quality characteristic, QS, as a measure of a system's performance. Recall that QS 2 is equal to the ratio of SNR2 to the spatial resolution in the appropriate power n that corresponds to the dimensionality of the images (n = 2 for projection imaging and n = 3 in the case of CT), normalized by the incident According to NRU, QS cannot be larger than unity or, more precisely, cannot be larger than the Epanechnikov constant, that is just slightly larger than 1 (de Hoog et al., 2014
; Gureyev et al., 2016
), for conventional absorption-based X-ray imaging systems. However, we have shown that in PBI setups the intrinsic imaging quality characteristic can be as high as cn(γ/NF)n/4, where γ is the ratio of the real decrement to the imaging part of the of the imaged object, NF is the minimal Fresnel number and cn is a dimensionless quantity of the order of 1. This means that in PBI the maximum number of bits of information about the imaged sample delivered by each detected photon can be close to (γ/NF)n/4. In hard X-ray PBI of biological samples, the ratio γ/NF can be significantly larger than unity, which indicates potential for a beneficial violation of the NRU.
In order to objectively assess the improvement of biomedical X-ray imaging quality in PBI, one typically needs to evaluate the key image parameters, such as contrast, SNR and spatial resolution, as a function of not only the imaging setup but also with respect to the radiation dose delivered to the sample in the process of imaging. The latter is particularly important in the case of medical X-ray imaging. For this purpose, in the present paper we introduced the biomedical X-ray imaging quality characteristic, QC, which was obtained by incorporating the contrast and radiation dose into the intrinsic imaging quality characteristic. Essentially, while QS evaluates the objective performance of the imaging system alone, QC allows one to assess the imaging quality of the system with respect to a particular sample or type of samples, by appropriately weighting the imaging quality metric in proportion to the contrast and in inverse proportion to the dose. We showed that the biomedical X-ray imaging quality in CT is directly proportional to the image contrast and is maximized with respect to the X-ray absorption when the average transmission through the sample is equal to 1/e2 ≅ 0.135. We also showed that the biomedical X-ray imaging quality in PB-CT is typically much smaller than unity because of the mathematical ill-posedness of the CT reconstruction (Natterer, 2001
). We have previously compared objective assessments of imaging setups in terms of the measured QS with the systematic subjective evaluation of the same images by groups of radiologists and medical imaging specialists (Baran et al., 2017
; Donato et al., 2025
; Pakzad et al., 2025
). In future work, it will be important to also compare the (objective) biomedical X-ray imaging quality characteristic with the assessment of images by radiologists.
We have also analysed the results of an experiment conducted at IMBL (Australian Synchrotron) in which two full fresh mastectomy samples were imaged at 4 mGy MGD in PB-CT mode with 5 m sample-to-detector distance, using monochromatic parallel 32 keV X-rays and two different detectors. One of the detectors was a flat-panel detector with a pixel size of 99 µm and a spatial resolution of approximately 150 µm, while the other detector was a photon-counting one with 75 µm pixels and an approximately single-pixel PSF. We showed that the difference in the quality of the CT-reconstructed images obtained from the scans under the same conditions using the two different detectors was consistent with the theoretical results obtained earlier in this and related papers. In particular, the gain in SNR in PB-CT in comparison with conventional CT at the same dose and spatial resolution, was significantly larger in the 3D case, compared with the SNR gain in the individual 2D projections. In the case of the detector with 75 µm pixels, the gain in SNR in PB-CT was demonstrated to be larger than a factor of 10. Note that, since the radiation dose is proportional to SNR2, a gain in SNR by a factor of 10, without loss of spatial resolution or contrast, corresponds to a potential reduction of the radiation dose by a factor of 100, without a loss of image quality. Even higher SNR gains (of the order of 102) and dose reduction factors (of the order of 104 to 105), compared with equivalent conventional absorption-based imaging, were demonstrated previously (Kitchen et al., 2017
) in a synchrotron-based small animal X-ray imaging experiment with a spatial resolution of approximately 40 µm. The present results, obtained for full intact human mastectomy samples under clinically acceptable radiation dose and spatial resolution, is an important step in our ongoing work towards breast cancer imaging of human patients at the Australian Synchrotron in the near future.
Note that the theoretical expressions for the 3D gain factor in equations (22)
–(23)
and for the biomedical X-ray imaging quality characteristic in equation (24)
are based on a number of approximations and are subject to certain validity conditions (such as, for example, the approximate of the sample and the near-Fresnel imaging conditions) which may or may not be satisfied in a particular experiment. Therefore, it is advisable to check the status of these conditions in a given experiment before attempting to apply these equations for estimation of the 3D gain coefficient or biomedical X-ray imaging quality in a particular setup. Furthermore, even when the validity conditions are satisfied, the approximate nature of these theoretical formulas means that they are likely to provide only qualitatively accurate predictions for the results of the relevant measurements.
APPENDIX A
Calculation of the gain factor in PB-CT using NRU
Applying equation (6)
and then equation (8)
to the case of the filter function O(r) = T3,inv(r), we can obtain
The first equality in equation (27)
is the result of the assumption about having a fixed incident fluence and using `full' phase retrieval with a = [δRλ/(4πβ)]1/2, as discussed in Section 5
in the paragraph following the one containing equation (20)
. The second equality in equation (27)
follows from equation (6)
. The third equality in equation (27)
is based on the fact that, due to the conservation in free-space propagation, SNR2 is the same in (flat regions of) CT-reconstructed volumes in the object and image planes: SNR2[P−1CR] = SNR2[P−1C0]. Finally, the last equality in equation (27)
is the consequence of equation (8)
, with O = T3,inv and I = P−1CR.
Using the definition of and the fact that the convolution with O = T3,inv does not change the first integral norm, we can express the right-hand side of equation (27)
as . We assume that the width of the PSF in the image plane after the CT reconstruction is much narrower than T3,inv, or, equivalently, that the Fourier transform of P−1CR is much broader and is almost constant on the support of the Fourier transform of T3,inv. Then using the Parseval theorem, we obtain
where
=
, f(A) =
and A
. Note that f(∞) = 1. From equations (27)
and (28)
, we now obtain equation (22)
.
APPENDIX B
Calculation of the gain factor in PB-CT using noise power spectrum
It is known (Goodman, 2000
) that the noise variance is equal to the integral of the noise power spectrum (NPS). In particular, =
, where WR,β(q) is the NPS of β(r) reconstructed from a PB-CT scan collected at z = R. An explicit expression for WR,β(q) can be obtained in conjunction with equation (6)
, in the case of a PB-CT scan with Ma flat-field corrected projections uniformly distributed over the 180° of sample rotation and uncorrelated with photons on average per detector pixel (Nesterets & Gureyev, 2014
),
where is the linear size of the (square) detector pixels. We assumed for simplicity that the detector has a rectangular-shape modulation transfer function χ with the width equal to
, where
is the width of the detector PSF. By integrating equation (29)
with respect to ρ and η over the range of all detectable frequencies, we obtain
where
is the mean number of photons registered in one detector pixel during the CT scan. Unfortunately, we are not aware of a method that would allow one to calculate the double-integral in equation (30)
exactly analytically. Therefore, the following results rely on approximations, similar to the results given by Nesterets & Gureyev (2014
). The behaviour of equation (30)
is determined primarily by the ratio =
. When
<< 1, the denominator in equation (30)
can be approximated by 1, and the simple integration results in
The latter result becomes exact in the case of a = 0, i.e. in the case of conventional CT with R = 0. On the other hand, when >> 1, equation (30)
can be approximated as
Taking the square root of the ratio of equations (32)
and (31)
, we obtain equation (23)
.
APPENDIX C
Biomedical X-ray imaging quality in CT
Here we investigate why the measured values of SNR[β0] in Table 2
were significantly lower than the measured values of SNR[I] in Table 1
. A related question is about the measured values of being significantly smaller than unity. It can be shown from equation (31)
, or from the results of Gureyev et al. (2016
), that
where =
. Note that the mean number of photons collected in a single detector pixel with the size
during the CT scan is
=
. Equation (33)
implies, in particular, that SNR[β0] is proportional to the fourth power of the spatial resolution, as mentioned earlier in the paper,
In the experiment considered in Section 6
, the effective spatial resolution of the Eiger detector was = 83.2 µm, Ma = 600 and the 2D was Iin,2D ≅ 0.19 µm−2. Using the known values for the linear absorption coefficient of adipose tissue at 32 keV, μ ≅ 2.62 × 10−5 µm−1 (TS-Imaging, 2025
), we get
2.62 × 10−5 µm−1 × 8.64 cm
2.26. It then follows from equation (34)
that SNR[β0] ≅ 0.69. This theoretical value is close to the average measured value SNR[β0] ≅ 0.52 for the Eiger detector in Table 2
and is also substantially lower than the values of SNR[I] in Table 1
. The last fact does not have a fundamental nature, because, at a given radiation dose, the ratio of SNRs in the CT projections and in the reconstructed slices of β depends on the number of projections in the CT scan.
Substituting equation (34)
and the expression for the MGD
into equation (17a), we obtain
where MR = is equal to π/2 times the number of spatial resolution units in the radius of the CT reconstruction cylinder. In the experiment described in Section 6
, we had Cm ≅ 0.27,
2.26,
0.5 and MR = π × 5.5 cm/(2 × 83.2 µm) ≅ 1038. Substituting these numbers into equation (35)
, we obtain
. The latter theoretical value is close to the value
measured in the experiment in Section 6
. This theoretical estimation confirms that in the present configuration should indeed be expected to be much smaller than unity. It was already noted in Section 4
that ≃ 1 requires SNR ≅ 103, while in the present configuration we had SNR[β0] ≅ 0.52. The key reason for the relatively small value of
in CT reconstruction is the ill-posedness of the CT reconstruction, which leads to the appearance of the (typically, large) numerical factor MR in the denominator of equation (35)
. At the same time, the term in the numerator of equation (35)
has the theoretical maximum equal to approximately 0.54, which is achieved at = 2 (the latter value is quite close to the value
≅ 2.26 in our experiment described in Section 6
). Note that the SNR for `projected' values of β, e.g. for β values integrated or averaged along one of the spatial coordinates, can be several orders of magnitude higher than for the 3D distribution of β. In other words, the biomedical X-ray imaging quality of the synthetic 2D mammographic images, obtained from CT-reconstructed slices, can be close to unity. This shows that the `dose fractionation theorem' (Howells et al., 2009
; McEwen et al., 1995
) does not hold for the types of CT imaging considered in the present paper, which is again explained by the mathematical ill-posedness of CT reconstruction; see also a related discussion by Gureyev et al. (2018
).
Acknowledgements
We gratefully acknowledge the two donors of mastectomy samples, without whom this work would not have been possible. The experimental part of this research was undertaken on the Imaging and Medical beamline at the Australian Synchrotron, part of ANSTO.
Conflict of interest
The authors declare no conflict of interest in regards to the present paper.
Data availability
The experimental data used in the present paper cannot be made publicly available due to the patient data privacy considerations.
Funding information
The following funding is acknowledged: National Health and Medical Research Council (grant No. APP2011204).
References
Arhatari, B. D., Stevenson, A. W., Abbey, B., Nesterets, Y. I., Maksimenko, A., Hall, C. J., Thompson, D., Mayo, S. C., Fiala, T., Quiney, H. M., Taba, S. T., Lewis, S. J., Brennan, P. C., Dimmock, M., Häusermann, D. & Gureyev, T. E. (2021). Appl. Sci. 11, 4120. Web of Science CrossRef Google Scholar
Baran, P., Pacile, S., Nesterets, Y. I., Mayo, S. C., Dullin, C., Dreossi, D., Arfelli, F., Thompson, D., Lockie, D., McCormack, M., Taba, S. T., Brun, F., Pinamonti, M., Nickson, C., Hall, C., Dimmock, M., Zanconati, F., Cholewa, M., Quiney, H., Brennan, P. C., Tromba, G. & Gureyev, T. E. (2017). Phys. Med. Biol. 62, 2315–2332. CrossRef CAS PubMed Google Scholar
Beltran, M. A., Paganin, D. M. & Pelliccia, D. (2018). J. Opt. 20, 055605. CrossRef Google Scholar
Bezak, E., Beddoe, A. H., Marcu, L.G., Ebert, M. & Price, R. (2021). Johns and Cunningham's The Physics of Radiology, 5th ed. Springfield: Charles C. Thomas. Google Scholar
Brombal, L. (2020). J. Instrum. 15, C01005. CrossRef Google Scholar
Brombal, L., Donato, S., Dreossi, D., Arfelli, F., Bonazza, D., Contillo, A., Delogu, P., Di Trapani, V., Golosio, B., Mettivier, G., Oliva, P., Rigon, L., Taibi, A. & Longo, R. (2018). Phys. Med. Biol. 63, 24NT03. CrossRef PubMed Google Scholar
Bronnikov, A. V. (1999). Opt. Commun. 171, 239–244. Web of Science CrossRef CAS Google Scholar
Dectris AG. (2025). Eiger2 detectors, https://www.dectris.com/en/detectors/x-ray-detectors/eiger2/eiger2-for-synchrotrons/eiger2-x-cdte/#Specifications (accessed 12 May 2025). Google Scholar
de Hoog, F., Schmalz, G. & Gureyev, T. E. (2014). Appl. Math. Lett. 38, 84–86. CrossRef Google Scholar
Donato, S., Caputo, S., Brombal, L., Golosio, B., Longo, R., Tromba, G., Agostino, R. G., Greco, G., Arhatari, B., Hall, C., Maksimenko, A., Hausermann, D., Lockie, D., Fox, J., Kumar, B., Lewis, S., Brennan, P. C., Quiney, H. M., Taba, S. T. & Gureyev, T. E. (2025). Med. Phys. 52, e17950. CrossRef PubMed Google Scholar
Goodman, P. (2000). Statistical Optics. New York: Wiley. Google Scholar
Gureyev, T. E., Brown, H. G., Quiney, H. M. & Allen, L. J. (2022). J. Opt. Soc. Am. A 39, C143–C155. CrossRef CAS Google Scholar
Gureyev, T. E., Kozlov, A., Nesterets, Y. I., Paganin, D. M., Martin, A. V. & Quiney, H. M. (2018). IUCrJ 5, 716–726. CrossRef CAS PubMed IUCr Journals Google Scholar
Gureyev, T. E., Kozlov, A., Paganin, D. M., Nesterets, Y. I. & Quiney, H. M. (2020). Sci. Rep. 10, 7890. Web of Science CrossRef PubMed Google Scholar
Gureyev, T. E., Nesterets, Y. I. & de Hoog, F. (2016). Opt. Express 24, 17168–17182. CrossRef PubMed Google Scholar
Gureyev, T. E., Nesterets, Y. I., de Hoog, F., Schmalz, G., Mayo, S. C., Mohammadi, S. & Tromba, G. (2014). Opt. Express 22, 9087–9094. CrossRef PubMed Google Scholar
Gureyev, T. E., Nesterets, Y. I., Kozlov, A., Paganin, D. M. & Quiney, H. M. (2017). J. Opt. Soc. Am. A 34, 2251–2260. CrossRef Google Scholar
Gureyev, T. E., Nesterets, Y. I., Stevenson, A. W., Miller, P. R., Pogany, A. & Wilkins, S. W. (2008). Opt. Express 16, 3223–3241. CrossRef PubMed Google Scholar
Gureyev, T. E., Paganin, D. M., Myers, G. R., Nesterets, Ya. I. & Wilkins, S. W. (2006). Appl. Phys. Lett. 89, 034102. Web of Science CrossRef Google Scholar
Gureyev, T. E., Paganin, D. M. & Quiney, H. M. (2024). J. Synchrotron Rad. 31, 896–909. CrossRef CAS IUCr Journals Google Scholar
Gureyev, T. E., Stevenson, A. W., Nesterets, Ya. I. & Wilkins, S. W. (2004). Opt. Commun. 240, 81–88. CrossRef CAS Google Scholar
Higgins, P. D., Attix, F. H., Hubbell, J. H., Seltzer, S. W., Berger, M. J. & Sibata, C. H. (1992). Mass Energy-Transfer and Mass Energy-Absorption Coefficients, including In-Flight Positron Annihilation for Photon Energies 1 keV to 100 MeV, NISTIR 4812. NIST, Gaithersburg, MD, USA (https://nvlpubs.nist.gov/nistpubs/Legacy/IR/nistir4812.pdf). Google Scholar
Howells, M. R., Beetz, T., Chapman, H. N., Cui, C., Holton, J. M., Jacobsen, C. J., Kirz, J., Lima, E., Marchesini, S., Miao, H., Sayre, D., Shapiro, D. A., Spence, J. C. H. & Starodub, D. (2009). J. Electron Spectrosc. Relat. Phenom. 170, 4–12. Web of Science CrossRef CAS Google Scholar
Hubbell, J. H. & Seltzer, S. W. (2004). X-ray Mass Attenuation Coefficients, NIST Standard Reference Database 126. NIST, Gaithersburg, Maryland, USA (https://dx.doi.org/10.18434/T4D01F). Google Scholar
Johns, P. C. & Yaffe, M. J. (1985). Med. Phys. 12, 289–296. CrossRef CAS PubMed Google Scholar
Kato, H. (2014). Jpn. J. Radiol. Technol. 70, 684–691. CrossRef Google Scholar
Kitchen, M. J., Buckley, G. A., Gureyev, T. E., Wallace, M. J., Andres-Thio, N., Uesugi, K., Yagi, N. & Hooper, S. B. (2017). Sci. Rep. 7, 15953. Web of Science CrossRef PubMed Google Scholar
Liu, Q., Suleiman, M. E., McEntee, M. F. & Soh, B. P. J. (2022). J. Radiol. Prot. 42, 011503. CrossRef Google Scholar
Mandel, L. & Wolf, E. (1995). Optical Coherence and Quantum Optics. Cambridge University Press. Google Scholar
McEwen, B. F., Downing, K. H. & Glaeser, R. M. (1995). Ultramicroscopy 60, 357–373. CrossRef CAS PubMed Google Scholar
Momose, A. (2020). Phys. Med. 79, 93–102. Web of Science CrossRef PubMed Google Scholar
Natterer, F. (2001). The Mathematics of Computerized Tomography. Philadelphia: SIAM. Google Scholar
Nesterets, Y. I. & Gureyev, T. E. (2014). J. Phys. D Appl. Phys. 47, 105402. Web of Science CrossRef Google Scholar
Nesterets, Y. I., Gureyev, T. E., Mayo, S. C., Stevenson, A. W., Thompson, D., Brown, J. M. C., Kitchen, M. J., Pavlov, K. M., Lockie, D., Brun, F. & Tromba, G. (2015). J. Synchrotron Rad. 22, 1509–1523. Web of Science CrossRef IUCr Journals Google Scholar
Nielsen, M. A. & Chuang, I. L. (2010). Quantum Computation and Quantum Information, 10th Ann. ed. Cambridge University Press. Google Scholar
Paganin, D., Gureyev, T. E., Mayo, S. C., Stevenson, A. W., Nesterets, Ya. I. & Wilkins, S. W. (2004). J. Microsc. 214, 315–327. CrossRef PubMed CAS Google Scholar
Paganin, D., Mayo, S. C., Gureyev, T. E., Miller, P. R. & Wilkins, S. W. (2002). J. Microsc. 206, 33–40. Web of Science CrossRef PubMed CAS Google Scholar
Pakzad, A., Turnbull, R., Mutch, S. J., Leatham, T. A., Lockie, D., Fox, J., Kumar, B., Häusermann, D., Hall, C. J., Maksimenko, A., Arhatari, B. D., Nesterets, Y. I., Entezam, A., Taba, S. T., Brennan, P. C., Gureyev, T. E. & Quiney, H. M. (2025). arXiv:2505.05812. Google Scholar
Ronneberger, O., Fischer, P. & Brox, T. (2015). U-Net: Convolutional Networks for Biomedical Image Segmentation in Medical Image Computing and Computer-Assisted Intervention (MICCAI2015), edited by N. Navab, J. Hornegger, W. M. Wells & A. F. Frangi, pp. 234–241. Springer International Publishing. Google Scholar
Sakurai, J. J. (1967). Advanced Quantum Mechanics, Section 4.6. Massachusetts: Addison-Wesley. Google Scholar
Tavakoli Taba, S., Baran, P., Lewis, S., Heard, R., Pacile, S., Nesterets, Y. I., Mayo, S. C., Dullin, C., Dreossi, D., Arfelli, F., Thompson, D., McCormack, M., Alakhras, M., Brun, F., Pinamonti, M., Nickson, C., Hall, C., Zanconati, F., Lockie, D., Quiney, H. M., Tromba, G., Gureyev, T. E. & Brennan, P. C. (2019). Acad. Radiol. 26, e79–e89. CrossRef PubMed Google Scholar
Thompson, D. A., Nesterets, Y. I., Pavlov, K. M. & Gureyev, T. E. (2019). J. Synchrotron Rad. 26, 825–838. Web of Science CrossRef IUCr Journals Google Scholar
TS-Imaging (2025). X-ray complex refraction coefficient calculator, https://ts-imaging.science.unimelb.edu.au/Services/Simple/ICUtilXdata.aspx (accessed 12 May 2025). Google Scholar
This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.
access
journal menu



