In-plane image quality and NPWE detectability index in digital breast tomosynthesis

A rigorous 2D analysis of signal and noise transfer applied to reconstructed planes in digital breast tomosynthesis (DBT) is necessary for system characterization and optimization. This work proposes a method for assessing technical image quality and system detective quantum efficiency (DQEsys) for reconstructed planes in DBT. Measurements of 2D in-plane modulation transfer function (MTF) and noise power spectrum (NPS) were made on five DBT systems using different acquisition parameters, reconstruction algorithms and plane spacing. This work develops the noise equivalent quanta (NEQ), DQEsys and detectability index (d’) calculated using a non-prewhitening model observer with eye filter (NPWE) for reconstructed DBT planes. The images required for this implementation were acquired using a homogeneous test object of thickness 40 mm poly(methyl) methacrylate plus 0.5 mm Al; 2D MTF was calculated from an Al disc of thickness 0.2 mm and diameter 50 mm positioned within the phantom. The radiant contrast of the MTF disc and the air kerma at the system input were used as normalization factors. The NPWE detectability index was then compared to the in-plane contrast-detail (c-d) threshold measured using the CDMAM phantom. The MTF and NPS measured on the different systems showed a strong anisotropy, consistent with the cascaded models developed in the literature for DBT. Detectability indices calculated from the measured MTF and NPS successfully predicted changes in c-d detectability for details between 0.1 mm and 2.0 mm, for DBT plane spacings between 0.5 mm and 10 mm, and for air kerma values at the system input between 157 µGy and 1170 μGy. The linear Pearson correlation between the detectability index and threshold gold thickness of the CDMAM phantom was −0.996. The method implements a parametric means of assessing the technical image quality of reconstructed DBT planes, providing valuable information for optimization of DBT systems.


Introduction
The development of digital breast tomosynthesis (DBT) has opened the way to various implementations of plane-based quasi-volumetric imaging in mammography. As with all medical imaging equipment, quality control (QC) and routine performance testing is required at various stages of the equipment lifecycle (European Commission 2018). Guidance on QC testing suitable for DBT system is available from the European Reference Organisation for Quality Assured Breast Screening and Diagnostic Services (EUREF) in the form of version 1.03 of the DBT protocol (van Engen et al 2018). The protocol covers a range of technical tests that give some insight equipment factors that can be tracked over time and that influence system imaging performance. The works of Rodriguez-Ruiz et al (2016) and Marshall and Bosmans (2012) describe the application of a number of these tests, including the modulation transfer function (MTF), the artefact spread function (ASF) and simple noise analysis. In a more detailed approach-that goes beyond QC- Mackenzie et al (2017) describe Fourier based metrics measured on DBT projection images, useful for characterization of these devices, however the reconstructed planes are not assessed. An alternative is to directly measure imaging performance using phantoms with some degree of realism, containing some relevant imaging tasks (Kiarashi et al 2015, Cockmartin et al 2017, Ikejimba et al 2017, Glick and Ikejimba 2018. These methods provide direct insight into the visibility of critical structures such as microcalcifications and mass-like lesions, however the technical aspects of system performance are inevitably tied up or combined in a single quality figure or metric of performance.
A complementary approach uses technical image quality metrics to assess global system performance that includes scatter and reconstruction. Metrics such as presampling MTF, noise power spectrum (NPS), noise equivalent quanta (NEQ) and detective quantum efficiency (DQE) describe the performance of linear imaging systems (ICRU Report 54 1996, Cunningham 2000 and are used in mathematical model observers for task-based image quality assessment. In recent work, we have used cascaded linear system theory to include the influence of scattered radiation on signal and image noise on signal-to-noise ratio (SNR) transfer through an imaging system, including the scatter reduction device, in the form of the system DQE (DQE sys ) (Monnin et al 2017). Intended as an extension of the DQE concept to the whole imaging chain, DQE sys is evaluated in the presence of typical scattering phantoms, with the system antiscatter device and the detector considered in a cascaded model. This work however, was limited to projection images. With regard to DBT systems employing image reconstruction, cascaded system models in combination with the filtered back projection (FBP) reconstruction algorithm have been developed in the frequency (Fourier) domain Siewerdsen 2008, Zhao andZhao 2008). These models show good agreement with image characteristics for different test objects and quantitative measurements of MTF and NPS (Zhao et al 2009, Zhao et al-2017, Marshall and Bosmans 2012, Calliste et al 2017. Detectability of targets of different shapes and sizes estimated from Fourier-based imaging metrics have been shown to correlate well with observer data for a range of imaging systems and conditions in DBT (Gang et al 2010, Reiser and Nishikawa 2010, Richard and Samei 2010. These results provide an initial validation of Fourier-based imaging metrics for modeling and predicting imaging quality and detection performance for DBT systems. In this work, we develop a method for measuring the imaging performance metrics in the frequency (Fourier) domain for DBT, based on the approximation of a stationary and spatially invariant imaging system considered linear in a small-signal approximation. NEQ, system DQE and detectability index for the non-prewhitening model observer with eye filter (NPWE) were calculated for five DBT systems, from measurements of 2D MTF and NPS. The DBT systems involved in this study differ in terms of the scan angle, number of projections, tube motion during image acquisition, reconstructed plane spacing, and reconstruction algorithm. The NPWE detectability index calculated for these systems was compared to the contrast-detail (c-d) thresholds obtained with the CDMAM phantom.

Theory
We consider an imaging system composed of a grid, a detector and a reconstruction process. The mean photon fluence and scatter fraction (SF) are respectively given byq in and SF in at the system input (i.e. at the breast support table), and byq and SF out at the detector input (or image) plane. For all the theoretical developments below,q =q in and SF out = SF in for systems not using a grid. For the present study, we have assigned the coordinates in the reconstructed volume as z for the direction perpendicular to the reconstructed planes, x for the tube travel direction and y for the perpendicular front-back direction. Fourier-based metrics are based on the approximation of a stationary and spatially invariant imaging system. Although logarithmic transforms, iterative reconstructions and image processing algorithms applied to 2D projections are not linear processes, DBT systems will be considered linear in a small-signal approximation for a small range of exposure variations in the images Siewerdsen 2008, Zhao andZhao 2008). Such systems can be described as products of transfer functions for signal and noise in the Fourier space.

Impulse response function (IRF)
We consider a linear and spatially invariant imaging system with pixel values d that depend linearly on the input photon fluence q, with a constant offset a and a gain ∂d/∂q. Signal scattering is described by the convolution product of the input signal with the point spread function of the whole imaging system that includes detector blur, geometrical blur and scatter (PSF in equation (1)).
For a linear system, the gain ∂d/∂q is constant (invariant to exposure). Equation (1) can be expressed in the frequency domain (noted in capital letters) with the MTF, the modulus of the Fourier transform of the PSF (equation (2)).
The IRF is defined as the response of the imaging system to a Kronecker delta function. The system MTF (MTF sys ) is the IRF measured with scatter and normalized to 1.0 at its maximal value (Monnin et al 2017).
where IRF d and IRF q denote the IRF with amplitudes expressed in image pixel values d and photon fluence values q, respectively. The quantity ∆q is the amplitude of the impulse signal, i.e. the number of photons attenuated in the object used to produce the impulse signal (a 0.2 mm aluminum disc in our study). Since the maximal value of MTF sys is equal to 1.0, the gain of the imaging system is the maximal IRF d value (IRF d,max ) normalized by the impulse signal (equation (4)).
The following steps assume that the image used for IRF calculation and the region of image used for NPS calculation have the same system gain.

NEQ
Except for the very low frequency range (below 0.1 mm −1 ), the NEQ at the system input (NEQ in ) is given by equation (5) (Monnin et al 2017).
In the image plane, the NPS can be expressed in pixel values d (NPS d ) or in q values (NPS q ). The link between the two NPS is the gain squared, as shown in equation (6).
For imaging systems with small signal linearity, the NEQ in the image (noted NEQ) is expressed by equation (7a) (ICRU Report 54 1996). In the presence of scatter, the detector MTF is replaced by the system MTF measured with scatter, andq is the total (primary + scattered) photon fluence at the detector (Monnin et al 2017).
Equations (3) and (7a) give a generalized expression for the NEQ with IRF d and NPS d , without any signal linearization or normalization (equation (7b)).
It can be seen that the NEQ self-normalizes, thus measurement of the zero-frequency signal for MTF normalization is not needed, and equation 7(b) works for any system with a linear response or linear in a small-signal approximation. This requires that the object used for the IRF calculation and the region of image used for the NPS calculation have the same system gain. For pre-processed images in projection imaging ('For Processing' images), the gain is held constant and MTF and NPS can be calculated from different images. For DBT images, the gain applied can vary from image to image, depending on the image content, and therefore the IRF d and NPS d have to be calculated from the same image. For images acquired with scattered radiation, MTF sys , IRF d and NEQ will vary as a function of (1 − SF out ) 2 (Monnin et al 2017).

System DQE
The system DQE (DQE sys ) measures the detective efficiency of the whole system in the presence of scattered radiation (Monnin et al 2017). It is the ratio between the NEQ in the image (equation (7b)) and the NEQ at the system input (equation (5)).
An equivalent expression for DQE sys involving the MTF sys and the total grid transmission (T t ) is given in equation (8b).
Equation (9) defines the system NPS (NPS sys , expressed in mm −2 ), that would be equal toq · MTF 2 sys for a perfect imaging system.

NPWE model observer
We consider an object with a Fourier spectrum S obj f x , f y , f z that creates a peak radiant contrast ∆q at the detector. The peak radiant contrast is defined as the x-ray attenuation in the object point of maximal thickness. The detectability index (d ′ ) of the NPWE model observer for this object in a DBT plane is expressed by equation (10a), where VTF is the visual transfer function of the human eye , Prakash et al 2011. (3) and (10a) give equation (10b), a formulation for the detectability index for DBT reconstructed planes that uses the in-plane IRF d and NPS d .
The VTF used in this study is an approximation of the Barten contrast sensitivity function of the human eye (Barten 1990) with a general form given in equation (11) (Burgess 1994, Segui andZhao 2006).
The values of c and n depend on the detection task and observer viewing conditions. In this study, a VTF consistent with that used by Segui and Zhao (2006) was used, with n = 2.2, c = 0.02 mm 2 , and a viewing distance of 400 mm without magnification. The maximum VTF response occurs at 1.06 lp mm −1 (figure 5).
The ratio ∆q/∆q represents the proportion between the peak radiant contrast of the object to detect and the flat radiant contrast of the disc used for IRF d calculation. In this study, this ratio will be equal to 1.0 because d' will be calculated for aluminium discs of various diameters (diameters equal to the CDMAM discs) with a fixed thickness equal to that of the object used to produce the impulse signal (a 0.2 mm aluminum disc in our study). This ratio can be calculated from the properties of the known object (thickness and attenuation coefficient) or measured using a dosimeter or alternatively estimated from the 'For Processing' projections using the pixel values re-expressed in air kerma or mAs values (i.e. from linearized images). For a thin flat object (for example a disc) with a thickness T much thinner than the plane spacing ∆z (T >> ∆z), the Fourier spectrum of the object can be decomposed into an in-plane xy-shape function and a thin z-aperture given by sinc (πTf z ) ∼ = 1, as shown in equation (12).
The 3D IRF d integrated along the z-frequency axis gives the 2D in-plane IRF d , measured in the plane that contains the whole object and provides the maximum IRF d (Zhao and Zhao 2008).
The d' index given in equation (15) can be expressed as a product between a function of the disc diameter f (D) and a radiant contrast-to-noise ratio (∆q/σ q ) that depends on the air kerma (equations (16) and (17)). Where: (17) is a frequency-dependent factor that involves the signal (IRF d and S obj ), noise (NPS) and VTF, without any dependence on air kerma for quantum limited systems. The gain ∂d ∂q in the denominator is introduced via σ d = ∂d ∂q σ q .

C-d analysis
Considering a thin flat object of thickness T and linear attenuation coefficient µ in a homogeneous material, the radiant contrast (absolute difference in photon fluence) ∆q between the object and the homogeneous background at the detector input is given in equation (18).
Where γ is a power parameter that has to be adjusted to the data and is expected to be close to 1.0 for thin objects. For a quantum noise limited system, the standard deviation of pixel values increases as a power function ofq, σ q ∼ = k q ·q β , where k q is a gain factor and β ≈ 0.5. Under these assumptions, the detectability index given in equation (16) increases as a power function ofq (equations (19a) or (19b)).
The c-d curves of the CDMAM phantom follow the points of constant detectability (for a constant observer performance). The threshold gold thicknesses (T Au ) for the CDMAM gold discs correspond to a fixed threshold detectability index (d ′ = d ′ T ) given in equations (20a) or (20b).
Consider that a CDMAM gold disc is replaced by an aluminium disc of same diameter but of thickness T Al . Imaging the aluminium disc in the same conditions (photon fluenceq) as the CDMAM phantom will give a detectability index d ′ given by equation (19b), closely related to the threshold thickness of the CDMAM phantom given by equation (20b). The link between T Au and d ′ is given in equation (21).
From equation (21), we expect to find a log-log linear link between the threshold gold thickness T Au of a given disc of the CDMAM phantom and the detectability index d ′ calculated with the NPWE model for an aluminium disc of same diameter and thickness T Al .

Imaging setup
Six flat panel mammography systems were involved in this study: a Fuji Amulet Innovality, a GEHC Pristina, two Hologic units (Selenia Dimensions and 3Dimensions), and two Siemens units (Inspiration and Revelation). The technical characteristics of the DBT systems are given in table 1. A homogeneous phantom made of 20 mm poly(methyl) methacrylate (PMMA) +0.5 mm Al (>99% purity) +20 mm PMMA with lateral dimensions of 180 mm × 240 mm was used for this study. The homogeneous 0.5 mm Al plate has the composition of the base plate in the CDMAM phantom (CDMAM manual, Artinis). The tube voltage, anode/filter (A/F) and grid configuration obtained when this phantom was imaged in DBT mode under automatic exposure control (AEC) were used (table 2). Only the GEHC Pristina system used an antiscatter grid. Four additional tube current-time product (mAs) values were chosen bracketing the AEC mAs value, giving five different exposure levels and a factor of five range between maximum and minimum mAs. In order to generate the c-d data required, the 0.5 mm Al plate was replaced by the CDMAM in the PMMA stack, and 8 images of the CDMAM were acquired for each mAs, with the phantom slightly moved between each DBT acquisition. The CDMAM images were scored using the CDCOM module and c-d curves were generated using the standard processing method described by Young et al (2006). The cases where the in-focus plane was not obvious, planes around the plane considered to have the best focus were scored and the plane with the lowest contrast thresholds taken to be the in-focus plane. Uncertainty on the threshold gold thickness was estimated using a bootstrap method.
For IRF d calculations, a 0.2 mm thick aluminium disc of diameter 50 mm was positioned on the 0.5 mm Al plate of the phantom, at the reference point, 6 cm from chest wall side and laterally centred (European Commission 2006). This phantom, named 'NPWE phantom' , is equivalent to the CDMAM phantom for x-ray attenuation and scatter. A typical in-plane image of the NPWE phantom is shown in figure 1. Four DBT scans of the NPWE phantom were acquired at each mAs level, with the phantom slightly moved between the four acquisitions to produce four different samplings of the disc. DBT planes of the CDMAM and NPWE phantoms were generated using the standard clinical reconstruction (table 1). The standard (ST) and high resolution (HR) modes for the Fuji used tube rotation angles of 15 • and 40 • , respectively. The ST images were reconstructed with the FBP algorithm and with the new Iterative Super-resolution Reconstruction (ISR), with pixel sizes 0.15 and 0.10 mm, respectively. The HR images could be obtained only with FBP and a pixel size 0.10 mm. Siemens systems images were reconstructed with the standard FBP (STD) method for the Inspiration and the new Enhanced Multiple Parameter Iterative Reconstruction (EMPIRE) algorithm for both Inspiration and Revelation. A 1 mm reconstruction plane spacing was available for all systems, while two additional spacings of 0.5 mm and 10 mm were available for the Pristina.

Dosimetry
The reference plane for all dosimetric quantities in this study was the detector input plane. Air kerma at the system entrance (grid + detector) was measured with a calibrated ionization chamber positioned at the reference point on the breast support table, 1 cm behind the homogeneous phantom (40 mm PMMA + 0.5 mm Al). For this measurement, the phantom was positioned in the compression plate (figure 2). The measured air kerma (cumulated air kerma for all the projections of a DBT acquisition) includes primary and scattered radiation at the system entrance plane; an inverse squared correction for distance gave the input air kerma at the detector plane, denoted K in in table 2. The product of K in with the photon fluence per air kerma unit (Boone 1998), denoted φ in table 2, gave the input photon fluenceq in . The detector air kerma (DAK) and photon fluence (q) are given, respectively, by the product of K in andq in with the grid total transmission (T t ). No correction was made for the transmission of carbon fibre cover of the breast support table.
The radiant contrast of the 0.2 mm Al disc used in the NPWE phantom (∆q in ) was measured using the homogeneous PMMA phantom positioned as in figure 2. The radiant contrast ∆q in is the difference in photon fluenceq in measured with and without an additional 0.2 mm Al plate fixed at the tube output. We assume the relative radiant contrast ∆q/q at the detector is equal to ∆q in /q in , neglecting the effect of the thin 0.2 mm Al disc on the SF at the system input (equation (23)).
The term∆q/q measured for the 0.2 mm Al disc depends on the disc thickness (T Al ) and on the linear attenuation coefficient of aluminium (µ Al ). For a given disc thickness, ∆q/q will vary only with the x-ray beam energy, and can be measured independently of the distance to the source, tube load or SF produced in the phantom.
Alternative methods can be envisaged for measuring the radiant contrast ∆q/q, for example using 'For Processing' projections with pixel values re-expressed in air kerma or mAs values (linearized images). 2D images can be used for systems that use the same A/F in 2D and DBT, while DBT 'For Processing' projections have to be used for systems where A/F is different in DBT.

SFs and grid transmission
The SFs were measured for the homogenous phantom using the beam stop method described in Monnin et al (2017). Lead discs of radii ranging between 1.5 mm and 6 mm and positioned at the reference point on the homogeneous phantom were imaged in the same conditions as the NPWE and CDMAM phantoms. The pixel values measured in the discs on the 0 • projection images (DICOM 'For Processing') were linearized to DAK values for the calculation. The SF measured with grid out gave the SF at the system input (SF in ) used for DQE sys calculations. For the Pristina, the SF was also measured with grid in, giving the SF at the grid output (SF out ). The ratio of DAK measured using a 5 × 5 mm ROI at the reference point in the linearized projection images acquired with grid in and grid out gave the total grid transmission T t . The primary grid transmission T p was calculated from T t , SF in and SF out (equation (24)).
The grid DQE (DQE grid ) is the ratio between T 2 p and T t (Monnin et al 2017).
The x-ray path length increases as a function of the inverse of the cosine of the projection angle, leading to a 10% increase for the maximum projection angle of 25 • (Siemens systems). SF thus increases with the projection angle. However, SF difference for a DBT acquisition compared to the 0 • projection will remain below a few percent for the small angles considered in DBT. Variations in SF due to x-ray path length variations as a function of the projection angle were therefore not considered in this study.

In-plane IRF d and MTF sys
The in-plane 2D IRF d was measured from the in-plane image of the 0.2 mm Al disc in the NPWE phantom with a radial version of the angled edge method described in Samei et al (1998). Calculation steps are described in Monnin et al (2016). A total of 180 radial edge spread functions (ESF) were calculated over 360 • from a square ROI 100 × 100 mm centred on the disc, using a 4 • angular aperture with a 2 • angular pitch. To reduce the noise in the measurements, for each mAs, the ESF was calculated for the four in-plane images of the NPWE phantom, with all the pixel values plotted as a function of their distance from the centre of the disc. The IRF d measured from the aluminium disc at 20 mm above the table with the SF produced by the phantom gives an estimate of sharpness for the gold discs of the CDMAM phantom.
The 0.2 mm aluminium thickness was chosen to produce a sufficiently small signal (radiant contrast below 15%) so that the system gain can be considered constant within the range of exposure variations produced by the disc in the image, and the log-normalization used before tomographic reconstruction can be considered linear in a small-signal approximation. The use of a small thickness is also necessary to avoid errors in IRF d measurement due to shadows produced at the edges of the disc at the maximal acquisition angles. As given in equation (3), IRF d is proportional to the system gain ∂d/∂q and to the radiant contrast of the disc ∆q that depends on the disc thickness and attenuation coefficient, and onq (equation (23)). The system MTF (MTF sys ) given in equation (3) is the IRF d normalized to 1.0 by its peak value. Unlike IRF d , MTF sys does not vary withq, system gain or disc properties, and was used to compare the signal frequency transfer between the different systems. The 1D MTF sys is a radial average of the corresponding 2D MTF sys .

NPS
Three dimensional NPS d data were measured from two homogeneous volumes of interest (VOIs) of 60 × 60 × 30 mm placed on both sides of the aluminium disc of the NPWE phantom, at 60 mm from the disc centre and 60 mm from the chest wall side, at in-focus plane (figure 1). The 3D NPS d for each mAs are averages of eight non-overlapping VOIs acquired at each mAs setting. No detrending correction was applied to subtract large inhomogeneities from the VOIs before noise analysis. Details of the NPS implementation are given in Monnin et al (2014). The in-plane 2D NPS d were obtained by integrating the 3D NPS along the z-axis (Zhao and Zhao 2008), giving the 2D synthesized NPS as defined by Siewerdsen et al (2002). The 1D NPS d curves are radial averages of the 2D NPS d , excluding the 0 • and 90 • axial values. The system NPS (NPS sys ) was calculated from the calculated NPS d , system gain, SF in and T t using equation (9). NPS sys was used to compare the noise frequency content of the images for the different systems, independently of the system gain, input SF and grid transmission.

NEQ, DQE sys and detectability index
The NEQ and DQE sys were calculated from equations (7b) and (8a), respectively. The detectability index d' of the NPWE model was calculated from equation (15), where ∆q = ∆q for the 0.2 mm aluminium thickness. The Fourier spectrum S obj for a disc of diameter D and thickness T can be expressed with equation (26a), where J 1 is the Bessel function of the first kind.
Because CDMAM discs are much thinner than the reconstructed plane spacing, the sinc function is approximately 1.0 below the Nyquist frequency, i.e. within the z-frequency range [0; 0.5/T z ], and therefore equation (26b) was used for d' calculation in equation (15).

X-ray beams and grid characteristics
The AEC settings obtained for the homogenous phantom for the different mammography units, the photon fluence per unit air kerma and the relative radiant contrast measured for the 0.2 mm Al thickness are given in table 2. The radiant contrast for the 0.2 mm Al disc varied between 10.7% and 13.2% for the different systems and was considered small enough such that the IRF was acquired on a linear region of the response for the different systems. The air kerma at the system input (K in ) obtained for the homogenous phantom under AEC (bold in table 2), varied between 282 µGy (Siemens Revelation) and 627 µGy (Selenia Dimensions). The SFs measured at the system (grid + detector) input (SF in ) varied between 0.413 and 0.450 for the different x-ray beams (table 2 and figure 3). Differences in x-ray spectra and in equipment can explain these differences of a few percent in SF in (range <4%) between the systems. Only the GEHC Pristina system uses a grid for DBT. The SF at the grid output (SF out ) was reduced to 0.112 for the Pristina. The corresponding primary (T p ), scatter (T s ) and total (T t ) grid transmissions were 0.746, 0.115 and 0.462, respectively, leading to a grid DQE equal to 1.20. Hence, the grid increased the NEQ and DQE sys by a factor 1.20 for the Pristina, for the phantom thickness considered.

In-plane MTF sys
For all the systems, the in-plane MTF sys results are strongly anisotropic, and act as band-pass filters along the f x -axis ( figure 4). The MTF sys measured on the images acquired at different mAs were identical. The different exposure levels were obtained from variations in tube current without changing the scan time and focus motion speed. The MTF sys for the different systems all have the same general bandpass shape, peaking between 0.7 mm −1 and 1.0 mm −1 along the x-frequency axis. These results are in good agreement with the cascaded models developed for DBT (Zhou et al 2007, Zhao andZhao 2008) and for cone beam computed tomography (CBCT) (Tward and Siewerdsen 2008). They are also consistent with in-plane MTF results previously measured in DBT (Zhao B et al 2009, Zhao C et al2017, Marshall and Bosmans 2012. MTF sys for the reconstructed planes incorporates signal blurring of the projections (detector resolution, focal spot blur and x-ray scatter) and additional contributions from the focus motion during acquisition and plane reconstruction (filtering and resampling). Flying focus systems suffer some additional blurring along the x-axis from the focus travel during acquisition of the projections while all the systems have pixel correlations introduced by the reconstruction filters and interpolations. The low-frequency MTF sys along the f x -axis is proportional to f x and sensitive to the acquisition angle (Zhao and Zhao 2008), while the high frequency response remains mainly determined by the presampling projection MTF and apodization filters. The scatter contribution can be practically considered as a constant factor equal to (1-SF out ) for spatial frequencies higher than 0.1 mm −1 (Monnin et al 2017).
The normalization of MTF sys by the peak value removes the influence of the system gain, photon fluence and edge contrast that are all multiplicative factors in IRF d (equation (3)). The radiant contrast for the 0.2 mm Al disc of ∼10% assumes the first order signal linearity remains dominant over nonlinear effects. For example, the logarithmic transform applied to the projections before reconstruction can be treated as a gain stage for small signal variations (Zhao and Zhao 2008). The thickness of the 0.2 mm Al disc represents a small fraction of the effective plane spacing (PSF in the z-direction) (Marshall and Bosmans 2012) and therefore the Al disc thickness can be considered infinitesimal and does not reduce high-frequencies in the z-direction. Consequently, high-frequency components of in-plane MTFs sys are not expected to be reduced by the object used for the calculation.
The mean radial MTF sys curves shown in figure 5 give a summary comparison of mean signal transfer between the systems and associated reconstruction parameters. The Selenia 3Dimensions and the Pristina with the 0.5 mm plane spacing have the highest MTF sys between 1.0 mm −1 and 4.5 mm −1 . For the Pristina, the thinner 0.5 mm planes give a MTF sys higher by 0.15 than the 1.0 and 10 mm planes in the frequency range 1.0-4.5 mm −1 . The Amulet ST and HR images reconstructed with FBP gave the poorest MTF sys . The new ISR algorithm improved dramatically MTF sys in the mid-range frequencies. The STD and EMPIRE algorithms give similar MTF sys on the Siemens Inspiration and Revelation.

In-plane NPS sys
In-plane NPS sys q curves, calculated using equation (9) and the input air kerma (K in ) given in table 2, are shown in figure 6. Nyquist frequencies are given in table 1. The use of the term NPS sys q removes the influence of the system gain, input SF, grid transmission and photon fluence at the detector from NPS d to give a normalized comparison (without units) of the noise for the different systems and reconstruction modes. Contributing factors to NPS sys q include the NPS of the projections (quantum, electronic and fixed pattern NPS), to which are added the pixel correlations introduced by 3D sampling, noise aliasing and filtering applied in the reconstruction processes. For all the systems, the in-plane NPS sys data show a strong anisotropic pattern in the Fourier plane, which can be expected as this parameter is mainly determined by the in-plane MTF sys squared. The NPS sys patterns for the different systems agree with NPS results previously measured in DBT (Zhao et al 2009), and are in line with the cascaded models developed for DBT (Zhao and Zhao 2008) and for CBCT Siewerdsen 2008, Baek andPelc 2011). Due to ramp filtering used for image reconstruction, the low-frequency NPS sys along the x-axis is proportional to the acquisition angle and to x-frequency squared (Zhao and Zhao 2008), and NPS sys along the y-axis is zero. Low-frequency trends in the reconstructed images are visible in figure 1 where pixel values near the centre of the PMMA phantom are higher than those near the periphery. These low-frequency trends increase the NPS below 0.2 mm −1 , as shown in figures 6 and 7, but the human VTF mitigates their influence on object detectability.
The mean radial NPS sys q in figure 7 peak between 0.5 and 1.4 mm −1 and show an averaged noise performance for the different systems. The Fuji ST and HR systems have a maximum radial NPS sys q between 3 and 3.5, higher than the other systems. The Selenia 3Dimensions and Pristina (for the 0.5 and 1 mm planes) show the lowest NPS sys q, with peak values around 2. The iterative EMPIRE algorithm gives similar NPS sys q for both the Siemens Inspiration and Revelation, with maximal value 2.2 lower by 25% than the STD reconstruction.

NEQ and DQE sys
The in-plane DQE sys data in figure 8, calculated using equation (8a), have lost the anisotropy introduced by the directional filters applied in the reconstruction processes, such as the ramp or apodization filters. DQE sys includes the loss in detection efficiency due to the detector, detector cover, grid efficiency (if used) and motion blurring (due to focus travel). Only factors that have an explicitly different effect on the IRF compared to the NPS can reduce DQE sys . An example is focus motion during the acquisition of the projections-this only reduces the signal in the travel direction (x-axis) and not the NPS. The deterministic effects of the ramp and apodization filters, present in the IRF and NPS, cancel out. The high DQE sys values along the f y -axis (higher than 1.0 for most of the systems) are due to division by values close to zero along the NPS f y -axis. DQE sys along the f y -axis have a high uncertainty and are therefore not reliable. For the Pristina, the grid increased DQE sys by 1.20 for the phantom used in this study. Because of the grid, DQE sys for the Pristina is expected to increase with SF in , and hence with breast thickness. The other DBT systems do not use a grid, and their DQE sys does not depend on SF in because IRF d is proportional to (1-SF in ) (Monnin et al 2017).
The mean radial DQE sys data in figure 9 give a global comparison of the average system efficiency. DQE sys for the different systems peaks between spatial frequencies of 0.6 mm −1 and 1.6 mm −1 . The highest and lowest peak DQE sys were for the Selenia 3Dimensions at 70% and for the Amulet with the HR mode at 25%. With a peak DQE sys value at 55%, the 0.5 mm plane spacing reconstruction for the Pristina gives DQE sys higher by 30% compared to the 10 mm plane spacing. The 1 mm plane spacing has similar DQE sys to the 0.5 mm planes below 1 mm −1 , and similar DQE sys to the 10 mm planes above 2 mm −1 . Hence, the 0.5 mm plane gives the highest DQE sys for all the frequencies. For the Amulet, ST mode is more efficient than HR mode over the whole frequency range, with peak DQE sys at 30% and 25% for the ST and HR modes, respectively. The new ISR algorithm increases the DQE sys peak from 30% to 42% for the ST mode, which is achieved at around 1.0 mm −1 for the two reconstructions. The STD and EMPIRE reconstructions available on the Siemens systems give similar maximal DQE sys around 50%. The maximal DQE sys for the Siemens is obtained at low frequency, around 0.5 mm −1 , followed by a quick falloff at higher frequency for the STD reconstruction, whereas the iterative EMPIRE algorithm maintains high DQE sys over a larger frequency bandwidth, dramatically increasing DQE sys between 0.5 and 2.5 mm −1 compared to STD. The EMPIRE algorithm gives similar DQE sys curves for the Inspiration and Revelation systems. The comparison between the FBP and new iterative reconstructions shows the importance of the reconstruction algorithm for overall system efficiency, and is a crucial step between the final (output) quality of the reconstructed DBT planes and the input photon fluence at the system entrance.
The mean radial NEQ results calculated with equation (7b) give a global comparison of the average image quality that depends on both the system efficiency and DAK chosen by the AEC for the different systems (figure 10). The 3Dimensions, Selenia Dimensions and Amulet ST with the ISR gave the three highest NEQ of the studied systems. The high DQE sys combined with a high exposure delivered by the AEC (table 2) explain these results. Despite the low DQE sys , the Amulet in HR mode gave the fourth highest NEQ peak because of a relatively high DAK compared to the other systems. The 10 mm planes for the Pristina and the STD reconstruction for the Inspiration provided the lowest NEQ, respectively below and above 1.2 mm −1 . Despite a similar DQE sys , the NEQ for the Siemens Revelation was 30% lower than for the Inspiration because of a DAK reduced by 30%.

NPWE detectability index and c-d analysis
Equation (15) gave the detectability indices (d') of the NPWE observer model for 0.2 mm aluminium discs with different diameters between 0.1 and 2 mm corresponding to the CDMAM gold discs ( figure 11(a)). Using the AEC settings, the 3Dimensions and Selenia Dimensions gave the best index for all the disc sizes, followed by the Amulet ST with ISR and the 0.5 and 1.0 mm planes for the Pristina. The Amulet ST with the FBP reconstruction gave the lowest detectability index for all the discs. The 10 mm planes for the Pristina gave a lower detectability than the 0.5 and 1.0 mm planes whatever the disc diameter. The EMPIRE outperforms the STD reconstruction of the Siemens Inspiration for all detail sizes. Compared to the Inspiration EMPIRE, the Siemens Revelation gave lower detectability because of lower DAK chosen by the AEC. These performance figures depend on the input dose used by the AEC settings and on the capabilities of the system. Figure 11(b) shows the detectability indices normalized by the square root of the entrance air kerma (d ′ √ K in ) used as a figure of merit (FOM) for the visibility of the discs of different sizes. The rank of FOM follows the DQE sys hierarchy.
The threshold gold thickness (T Au ) for CDMAM discs of different diameters corresponds to the smallest disc thickness visible in 62.5% of cases by a human observer (Young et al 2006). T Au values obtained with the CDCOM software for the in-plane images of the CDMAM phantom were plotted against d' indices calculated with the NPWE model. A negative log-log linear correlation was obtained, in agreement with equation (21) for all the systems and imaging conditions ( figure 12). This linear relationship does not depend  figure 12). The coefficient γ in equations (18)-(21) is expected to be close to 1.0, (∆q ∼ =qµT for a linear system) and the slope should therefore be close to −1.0. This difference in slope of about 15% is expected to come essentially from the radiant contrast given by the 0.2 mm Al disc on the images. The power law in equation (18) that gives the coefficient γ is only a linear approximation of contrast valuable for linear systems and thin objects. Nonlinearities in image reconstruction and processing can produce images in which the relative contrast on the image is dose dependent, and the small signal analysis used in this study accounts for a system behavior which is linear only around an operating point. The characteristics of the linear regression performed on data acquired over a dose range of five can therefore substantially deviate from the theory. As shown in equation (21), the linear link depends on the ratio of linear attenuation coefficients of the discs in the CDMAM and NPWE phantoms (gold and aluminium). Differences in x-ray beam energy can therefore shift the points and give an offset between the points for the different systems. In our study, the small differences in tube voltage (between 30 kV and 34 kV) and in A/F only caused small variations in the ratio µ Al /µ Au (variations within ±5% between 20 keV and 40 keV, below 1% in our study).
The linear correlation showed a good accuracy for all the disc diameters between 0.1 mm and 2 mm. Residuals of the linear regression did not show any bias related to disc diameter, input dose, reconstruction method or system. However, CDMAM results showed more variability for the largest discs, potentially because the threshold thickness can reach the minimum value available in the phantom for the largest disc in good imaging conditions, and the precision of the linear log-log model decreased substantially as a function of the disc diameter. Some error bars of d' values of the largest discs do not intersect the fitted line. The predictive power of the linear model is therefore substantially less good for discs of diameter 2.0 mm compared to discs of diameters 0.1 mm or 0.25 mm. Nevertheless, the experimental variations in T Au behaved with a good accuracy as predicted by the NPWE observer model, suggesting that the NPWE model is suitable for image quality characterization for low-contrast objects in reconstructed planes in DBT. The result obtained previously for 2D FFDM mammography (Monnin et al 2011) can be considered extended to reconstructed planes in DBT. The signal and noise transfers through the global system including detector, grid and reconstruction algorithm, considered as locally linear, can be calculated as functions of free parameters like DAK, beam energy, acquisition angle, breast thickness, plane spacing or reconstruction parameters.

Limitations of the study
The NPWE model gave coherent and robust results compared to the CDMAM phantom for a range of discs between 0.1 mm and 2 mm in diameter, for acquisitions angles between 15 • (Selenia Dimensions) and 50 • (Siemens Inspiration and Revelation), for plane thicknesses between 0.5 mm and 10 mm (GE Pristina), and for different reconstruction algorithms (FBP and iterative). The uncertainty in NEQ can be estimated along the uncertainties of the contributing MTF and NPS. The accuracy in MTF depends on the intrinsic measurement accuracy, estimated to be around 2.5% of the measured value (Buhr et al 2003), and on a frequency-dependent uncertainty up to 0.01 due to noise in MTF measurement. The resulting total uncertainty in MTF depends on the frequency and is at maximum 3.5%. The uncertainty in NPS is expected around 5%, as stated in IEC 62220-1. The resulting uncertainty in NEQ is at maximum 12%. The uncertainty in K in can be estimated around 3%, mainly due to uncertainty in dosemeter calibration. The resulting uncertainty in DQE is at maximum 15%. The uncertainty in d' is expected to be of the same order of magnitude as for the NEQ.
The image quality parameters (MTF sys , NEQ, DQE sys and NPWE detectability index) and CDMAM results were all assessed for objects thinner than the reconstructed plane spacing. Under this condition, the case can be treated as a 2D problem, giving in-plane imaging parameters. Hence, the disc used for the IRF measurement has to be thinner than the plane thickness. Plane spacing larger than the disc thickness can potentially make the calculated IRF d phase-sensitive, i.e. the IRF d amplitude will depend on the actual  z-position of the disc with respect to the reconstructed planes. This effect can be a source of variation in IRF d amplitude and a source of error in NEQ, DQE sys and d'. The thresholds thicknesses calculated for the thin CDMAM gold discs also have some phase dependence that does not cancel out with those of IRF d . The rate of variation of signal amplitude across the reconstructed planes is related to the z-PSF. A wide z-PSF for a DBT system will result in slowly varying signal between the reconstructed planes in the z-direction, and thus mitigate the phase-dependence of the object signal. Hence, we expected this source of error to be negligible for the systems involved in our study; this is supported by the close agreement found between CDMAM results and calculated d' results. For objects thicker than the plane spacing, the analysis must be extended to 3D. The resolution along the z-axis (z-MTF) has to be taken into account, leading to a 3D NEQ, DQE sys and NPWE model. In this study, the z-resolution was not measured, and the results in image quality and system efficiency for 3D objects could potentially differ from those obtained in this study, which only holds for thin (∼2D) objects. The z-resolution depends on both the acquisition angle and the reconstruction algorithm,  and the three-dimensional imaging performance is influenced by these parameters, but is not evaluated in this study. An extension of the model for extended objects should thus be considered. Thick objects can also give a high signal not compatible with the assumption of small signal system linearity used in this study. The challenge to this assumption of linearity comes from nonlinear iterative reconstruction techniques that can produce images in which resolution is contrast and noise (dose) dependent. As a consequence, sharpness and noise levels can vary locally in the reconstructed images, depending on the contrast of the structures present in the object being imaged. This can lead to a local spatial and signal level dependence for resolution and noise. The assessment of imaging quality using transfer functions is a linear approximation to a nonlinear system and must be ideally used under task-based conditions. The hypothesis of small signal system linearity used in this study accounts for the system behavior which is linear around an operating point (dose and contrast). This study was made under small contrast condition provided by the 0.2 mm Al disc. The good agreement between detectability indexes (d') and CDMAM threshold thicknesses (contrasts) for the different dose levels indicates that the proposed method yields a good approximation of the realized imaging performance for iterative reconstructions. This applies to the threshold details in CDMAM. For thick objects with high contrast, such an analysis with transfer functions in the frequency domain could potentially give erroneous results, and a statistical analysis performed in the spatial domain will be more appropriate (Zhao and Zhao 2008).
This study does not evaluate mean glandular dose (MGD) for these systems. MGD measurement would have required additional acquisitions using homogeneous PMMA blocks. This would also have introduced further complication in the study as the systems have different A/F combinations which in turn have different dose efficiencies, depending on spectral shape. These dose efficiencies will vary as a function of breast simulating material (PMMA in this case) and in fact on the breast dosimetry model used. By concentrating on the SNR transfer efficiency from the exit of the imaged object through to the reconstructed image, we are able to examine the influence of grid, detector and reconstruction without the added complication of spectral dose efficiency. Comparison of MGD between systems is certainly interesting, but for these reasons, out of the topic of this study.
Detectability performance using both the NPWE and CDMAM methods was established in homogeneous phantoms without anatomical noise. The ability of DBT to reduce anatomical noise is a major advantage compared to 2D mammography. A reduction in anatomical noise will tend to increase detection performance, especially for large objects like masses (Elangovan et al 2018, Hadjipanteli et al 2019. The ability of the system to reduce anatomical noise in the image plane was not taken into account in this study, as only the signal and detector noise transfer through the imaging systems were quantified. The presence of anatomical noise will reduce the detectability indices, and the magnitude of this reduction will depend on the ability of the system to reduce or restrict propagation of out of plane anatomical texture to the image plane containing the targets of interest. The comparison of detectability indices or threshold thicknesses calculated in homogeneous background is therefore not suitable for comparing projections and reconstructed DBT images. It is quite possible that the introduction of an anatomical noise term in the detectability index will modify the hierarchy of system performance found here. Systems with high z-resolution, large acquisition angle and a high number of projections are expected to produce reconstructed planes with a smaller amount of anatomical noise and this is expected to result in improved object detectability. It is also possible that the ranking of systems may depend on the amplitude of the anatomical noise in the imaged volume. The generalized NEQ (GNEQ) and DQE (GDQE) proposed by Barrett et al (1995) and Richard et al (2005) include in-plane anatomical noise and could be used for that purpose. This will be considered as an extension of the present model.

Conclusion
This study has introduced and applied a new methodology for assessing the image quality and imaging performance of reconstructed planes in DBT. Based on cascaded linear system theory, signal and noise transfer through the imaging system were described in the Fourier plane. The in-plane image quality (NEQ) and system efficiency (DQE sys ) measured with a thin (low-contrast) disc were used to compare the performance of different systems or reconstruction processes. When inserted into the NPWE observer model, these metrics correctly predicted the detectability of thin (∼2D) objects in reconstructed DBT planes. The strong correlation between NPWE indices and CDMAM threshold gold thickness demonstrated the robustness of the method for a wide range of dose, object size and reconstructed plane spacing. The measurement technique used in this work is proposed for the characterization of DBT imaging chains and the detectability assessment of thin objects with a low contrast.