A Model for a Linear a-Se Detector in Simulated X-Ray Breast Imaging With Monte Carlo Software

In-silico clinical trials with digital patient models and simulated devices are an alternative to expensive and long clinical trials on patient population for testing X-ray breast imaging apparatuses. In this work, we simulated a linear-response a-Se detector as an X-ray absorber, neglecting some physical processes, such as electro-hole tracking and thermal noise. In order to tune characteristics of the simulated images toward those of the clinical scanners, the detector response curve, modulation transfer function (MTF), and normalized noise power spectrum (NNPS) were measured on a clinical mammographic unit. The same tests were replicated in-silico via a custom-made Monte Carlo code in order to define a suitable model to modify simulated images and to have realistic pixel values, noise, and spatial resolution. The proposed approach resulted to restore the slope and the magnitude of the NNPS in simulated images toward curves evaluated on a clinical scanner. Similarly, the proposed strategy for tuning noise and spatial resolution in simulated images led to a contrast-to-noise ratio (CNR) evaluated on a custom-made phantom which differed from those in measured images less than 16% in absolute value.

A Model for a Linear a-Se Detector in Simulated X-Ray Breast Imaging With Monte Carlo Software A. Sarno , R. M. Tucciariello, M. E. Fantacci , A. C. Traino, C. Valero, and M. Stasi Abstract-In-silico clinical trials with digital patient models and simulated devices are an alternative to expensive and long clinical trials on patient population for testing X-ray breast imaging apparatuses.In this work, we simulated a linear-response a-Se detector as an X-ray absorber, neglecting some physical processes, such as electro-hole tracking and thermal noise.In order to tune characteristics of the simulated images toward those of the clinical scanners, the detector response curve, modulation transfer function (MTF), and normalized noise power spectrum (NNPS) were measured on a clinical mammographic unit.The same tests were replicated in-silico via a custom-made Monte Carlo code in order to define a suitable model to modify simulated images and to have realistic pixel values, noise, and spatial resolution.The proposed approach resulted to restore the slope and the magnitude of the NNPS in simulated images toward curves evaluated on a clinical scanner.Similarly, the proposed strategy for tuning noise and spatial resolution in simulated images led to a contrast-to-noise ratio (CNR) evaluated on a custom-made phantom which differed from those in measured images less than 16% in absolute value.

I. INTRODUCTION
I N-SILICO imaging clinical trials transpose clinical tests on patient population to virtual clinical trials (VCTs), with simulated apparatuses and digital models of the patients [1].They have the potential of reducing the costs and the time needed by horizon technologies to enter the clinical practice, also avoiding issues related to volunteers' exposure to ionizing radiations.Of particular interest is their use in the evaluations of apparatuses for the breast cancer screening and diagnosis, where 3-D imaging is strongly proposed as alternative to conventional mammography whose performance are limited by the tissue superimposition in the 2-D images [2], [3], [4], [5], [6], [7], [8].
At the today state of the art, there is large availability of digital models of the breast for VCT.These present different degrees of realisms and are derived either from clinical patient data or mathematical models [9], [10]; similarly, also databases of simulated lesions have been developed [11].As regard to the simulation of the imaging apparatuses, several authors have proposed analytical calculations based on raytracing algorithms coupled with post-processing algorithms for the inclusion of Poissonian noise and scatter Compton contribution [3], [4], [5], [6].On the other hand, these two are implicitly simulated in Monte Carlo-based applications [7], [12], [13].A fundamental block of the simulation chain is represented by the model of the detector [2].Hence, it does largely contribute to the sharpness of the produced images, as well as to the noise level which affect the lesion visibility.It introduces noise to the image, mainly due to the electronic and to the digitalization processes, and reduces the overall system spatial resolution causing blurring over edges and details [14].The detector block is usually left out of the simulation platform for the needed computation complexity and its characteristics are modeled via appropriate postprocessing.For example, in the VICTRE platform [7], [8], which relies on Monte Carlo software for images computations, the detector is not simulated, but the energy released in the direct conversion flat-panel detector is converted in charge per pixel.A Swank factor of 0.99 was assumed and at each pixel an additive noise was included summing a value of charges taken from a Poissonian distribution with average values of 5200 [7].A more complex model for the electronic noise was adopted in a recent paper [8].Main simplifications in the used detector model are related to the absence of image blurring introduced by electron-hole transport within the a-Se wafer, absence of tracking of secondary electrons, as well as absence of charge drifting.These simplifications may lead to a bias between the spatial resolution of the simulated and that of the real images acquired with clinical scanners [8].Moreover, the theoretical model adopted for the noise is not benchmarked versus that in real images acquired on clinical scanners.
In this work, we propose a method for tuning noise and spatial resolution in X-ray breast images simulated via a Monte Carlo software based on experimental measurements.The employed software was developed within the advanced Geant4-based platform for VCT in X-ray breast imaging (AGATA) project [12], [13], [15] which aimed at realizing a platform for VCTs meant for the evaluation of 3-D breast imaging techniques [e.g., digital breast tomosynthesis (DBT) and dedicated breast CT (BCT)] also in comparison to 2-D digital mammography (DM).The adopted computational breast phantoms were derived from 3-D clinical images of the breast with a relatively high spatial resolution [15].The simulation software is based on the Geant4 simulation toolkit [12], [13], for a fine reproduction of the imaging physics.In the AGATA platform, the simulated detector is modeled as an absorber layer of known material and thickness; neither the optical photons (in the case of indirect-conversion detector [16]) nor the electronic readout was simulated.In this context, this work compares the imaging performance of a-Se direct conversion detector with linear response adopted on a DM/DBT clinical scanner to that of the simulated detector in the AGATA platform and proposes a post-processing approach for increasing the realisms of the produced images.

A. Simulation Platform
The simulation platform [12], [13] is based on the GEANT4 simulation toolkit and the physics list Option4.Coherent and incoherent scatter as well as photoelectric interactions were simulated.The detector was modeled as a 0.200-mm thick wafer made of Se whose density was 4.26 g/cm3; the fluorescence in the a-Se material was simulated (Se k-edge = 12.7 keV, and k α fluorescence energy = 11.2 keV, k β fluorescence energy = 12.5 keV [14], [17]).Electrons cut-off was set to 5 keV.The spectra were computed as suggested in [18].In the spectra computation, the added filtration thickness was modulated to obtain half-value layer (HVL) equal to the measured one at any of the employed tube voltage.For computing the simulated images, dose deposits were scored in the detector wafer and added up in an image matrix.The location within the matrix was selected on the basis of the location of the interaction in the detector (detector dose map) with a resolution equal to the pixel pitch of the modeled detector.The 0.300-mm focal spot was replicated in silico as a circular source.

B. Detector Model
In order to convert values in the simulated images expressed as pixel dose (in mGy), the detector response curve was evaluated.This related the pixel values to the incident air kerma evaluated at the breast support.The air kerma was calculated as suggested in [19] by scoring the photons impinging on a defined 2 cm × 2 cm surface with the center at 5 cm from the chest wall side.The backscattered photons (i.e., the scattered photons from the breast support) were not considered in the calculation.In the case of the response curve of the detector mounted on the clinical scanner-showing a linear response in the used exposure ranges, this can be modeled as a linear curve [20] as follows: Here, PV M is the pixel value evaluated in a selected image region of interest (ROI) and K M is the measured incident air kerma at the breast support.A and B are suitable fitting parameters.Similarly, for the simulated detector, a linear curve to model the detector response is introduced as follows: (2) In this case, the C computed fit parameter relates the simulated pixel value (PV S,D in mGy) and the simulated incident air kerma at the breast support K S .Simulated and measured pixel values and incident air kerma were evaluated in analogous position in the simulated and real geometry in the clinical scanner.Fitting curves were computed via Origin pro 9.0 of OriginLab Corporation.The relation between the two linear models of the detector response curve permits to recalculate the pixel values in the simulated images to the pixel values range adopted by the scanner.The simulated pixel value I S,D (x, y) is then manipulated for each of arbitrary location (x, y) in order to present values in the range of those of the modeled clinical scanner following: Thereafter, in order to tune the spatial resolution in the simulated image I S (x, y) toward that of the measured ones, on the supposition of the linear-space-invariant response of the detector, I S (x, y, ) is convolved (R-tuned) by the r(x, y) filter r(x, y) is evaluated as the inverse Fourier transform of the R(v, f ) function with (v, f ) coordinates in the Fourier domain and it relates the measured modulation transfer function MTF (MTF M ) and the simulated MTF (MTF S ) as follows: Based on the assumption that MTF M and MTF S present Gaussian shape [20], also the R function may be supposed to present Gaussian shape with standard deviation (σ R ) derived from those of the measured (σ M ) and simulated (σ S ) curves as follows: A further assumption is made on the circular symmetry of MTF M , MTF s , and R.This allows to estimate the three standard deviations from 1-D measured and simulated MTF curves and project the resulting 1-D R function in (v, f ) plane by means of circular extrusion.This filtering process leads to the I R (x, y) R-tuned image.
Finally, differences between noise in R-tuned simulated images and measured images (σ AWGN ) will permit to define the amount of additive Gaussian white noise (AWGN) necessary to reach realistic levels.With this purpose the pixel standard deviation was evaluated in I R (σ PV,R ) and in the measured image I M (σ PV,M ) at the varying of the exposure level.The variance of the AWGN (i.e., the noise power) was evaluated as follows: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

C. Detector Characterization and Image Quality Assessment
The measurements were performed using a Hologic Selenia Dimension DBT scanner.It features a 24 cm × 29 cm a-Se direct conversion flat-panel detector with a wafer thickness of 0.200 mm and a pixel pitch of 70 μm.The source to detector distance was 70.0 cm and the air gap (i.e., the distance between the superior surface of the support paddle and the detector) was 2.5 cm.The spectra were generated with W/Rh anode/filter combination for voltages of 26, 29, and 31 kV.The HVL and the air kerma measurements were obtained via a calibrate Piranha RTI multimeter [21], [22] (uncertainty = 5%).
Tests for detector characterization were replicated both on the Hologic DM/DBT physical unit and on the in-silico AGATA platform.PV S,D and PV M were evaluated as the average pixel values in ROIs of 80×80 pixels region in the projected images at 5 cm from the chest-wall side in the simulated and measured images, respectively; the same ROIs location and dimensions were also adopted for σ PV,R and σ PV,M evaluations.The incident air kerma for detector response curve computations were evaluated above the support paddle at 5 cm from the chest-wall side.
The detector MTF was evaluated via the presampled MTF.This was assessed with a 12.5-μm-diameter tungsten wire placed on the support paddle and tilted 2 • with respect to the chest wall-to-nipple direction.
In order to evaluate the impact of the manipulation of the simulated images, also the noise power spectrum (NPS) was assessed as proposed in [23] for BCT.In this case, a PMMA block was imaged twice consecutively, and the difference image projections (I DD ) were used for the computation.100 ROIs comprising 256 × 256 pixels were selected from I DD (I DD,ROI ) on the PMMA block footprint and 2-D FFT was used to compute the 2-D NPS as suggested in [23] for the 3-D case With x the pixel pitch.The normalized NPS (NNPS) was computed by dividing the NPS 2D by the square of the pixel value in a large ROI across the PMMA block projection and the 1-D NPS was computed as the radial average of the NPS 2D [23].
A customized contrast phantom was used for image quality evaluations.It was constituted by a PMMA block as the one used for the NPS evaluation, including two aluminum sheets with thicknesses of 0.025 and 1.00 mm lying on the superior surface of the PMMA block.This phantom was used for the comparison of contrast-to-noise ratio (CNR) between Al inclusions and PMMA background evaluated as follows: where (PV Al , σ 2 Al ) and (PV PMMA, σ 2 PMMA ) are the (average, variance) pixel values evaluated in a 285 × 285 pixels ROIs over the 1.00-mm Al region and PMMA region, respectively.Similarly, the pixel standard deviation in the ROIs were

III. RESULTS
Table I reports the fitting coefficients of the detector response curves as indicated in (1) and ( 2).The last column of the table also reports the R2 fitting parameters for the measured detector response curve indicating high affinity of the linear detector response.
Fig. 1 compares measured MTF to that obtained for the simulated apparatus at 29 kV, as well as the computed R function.In this case, σ S resulted 8.43 mm −1 and σ M resulted 5.95 mm -1 , with the related standard deviation of the R Gaussian filter σ R of 8.39 mm -1 .The standard deviation of Gaussian fit curves of measured and simulated MFTs (σ M and σ S ) are reported in Table II along with the estimated σ R coefficients at 26, 29, and 31 kV.Computed σ AWGN for the three investigated tube voltages are reported in Fig. 2. They are shown as function of the expected pixel value, evaluated as the average value of the same ROI used for the variance estimates in (7).These curves show a linear trend, hence linear fit presented adj-R square parameters >0.99.Slope and intercepts for computing σ AWGN as function of the expected pixel value are reported in Table III.AWGN is added to the R-tuned images, pixelby-pixel, using the pixel values as expected ones for the computation of the σ AWGN and the corresponding Gaussian noise.
Fig. 3 compares measured and simulated NNPS at 26 kV Fig. 3(a), at 29 kV Fig. 3(b), and at 31 kV Fig. 3(c).The R-filters introduced pixel signal correlation which produced a slope in the NNPS curves similar to that of the measured curves.However, the detector model simplification, which does not include the electronic chain, led to a noise underestimation and AWGN whose standard deviation is derived from ( 7) is added.Hence, for each pixel of the R-tuned image, a random number taken from a Gaussian probability density function was added.The standard deviation of such a Gaussian probability density function is calculated from the curves in Table III, based on the value of the considered pixel.The AWGN permitted to rise the NNPS curves up to the measured ones, both at 26 and 29 kV.In these two cases, the used PMMA phantom was 18 and 33 mm thick, respectively.In the case of 31 kV, the NNPS curve derived from the simulated and R-tuned images of PMMA block presented slightly higher values of the measured one after the inclusion of the AWGN [Fig.3(c)].It is worth to consider that in this case, 43 mm thick PMMA phantom was used for the NNPS evaluation; it was thinner for evaluations at 26 and 29 kV.This may have determined a more conspicuous beam hardening and a related detector response which slightly differ from that foreseen by the model.Table IV compares CNR between the 1-mm Al foil and PMMA block evaluated for simulated and acquired images on clinical scanner.The thickness of the PMMA block was 28 mm for 26-kV tube voltage and 55 mm for 29 and 31 kV tube voltages.The ratio between simulated (following the R-tuning and the AWGN inclusion) and measured CNR at 26 kV resulted 0.93 indicating a difference of 7% in the qualitative image index.In the cases of 29 and 31 kV such a ratio resulted 0.84 and 0.88, respectively.Similarly, we compared the image noise, computed as the standard deviation in defined ROIs, over a region included in the projections of the two Al foils as well as over a region behind the PMMA.
V reports the between the pixel standard deviation for the simulated image phantoms the R-tuning and the inclusion of the AWGN) and the standard deviation in phantom images acquired on the clinical scanner.The lowest values for such a ratio can be observed for the ROI comprised in the PMMA background, ranges between 1.03 at 26 and 0.83 at 31 kV.These two values are 1.05 and 0.88 for an ROI over the 0.025-mm Al foil and increase up to 1.15 and 0.88 for pixel comprised over the projection of the 1-mm Al foil.

IV. DISCUSSION AND CONCLUSION
We measured the detector response curve, the MTF, and the NNPS of an a-Se detector with a linear response adopted in a clinical mammographic unit.This detector was replicated in-silico as a 0.200-mm a-Se absorber with a pixel pitch of 70 μm.The measurements were then replicated in-silico to evaluate the differences both in the noise level and spatial resolution, with several physical phenomena neglected in the simulated platform, such as electrons and electron-hole pairs tracking and thermal noise.The differences between simulated and measured MTF curves and the assumption of the linear space invariant characteristics of the detector, allowed to define a Gaussian filter for spatial resolution tuning in the simulated images.Such a filter was modeled taking into account the differences between the simulated and measured MTF width on the hypothesis that these may be modeled as a Gaussian curve.The defined filter showed to restore the pixel signal correlation present in the measured NNPS (measured and simulated R-tuned NNPS curves presented similar slope).However, the simplification of the detector model, which does include neither the creation and tracking of the electronhole pairs nor the detector electronic chain led to simulated images with a noise level lower to that of real images acquired on physical scanner.Hence, the precise amount of additive white Gaussian noise was estimated comparing pixel standard deviation in flat-field images acquired on the clinical scanner and those simulated and tuned via the Gaussian filter.The square root of the differences of variance of pixel values in the measured image projections and that in simulated projections after the R-tuning presented a linear dependence on the expected pixel value.This permits to calculate and add to each pixel the missing amount of AWGN on the basis of its value.The limitation on this approach is related to the difficulty in knowing the expected pixel values in simulated images, and the residual AWGN is calculated from the actual pixel value.The inclusion of a suitable level of additive Gaussian noise after the image filtering showed to modify the NPS in the simulated images toward the measured one.Little differences may be ascribed to neglected noise contribution due by detector imperfections (e.g., local variations in the a-Se thickness, impurities, etc.) as well as to the beam hardening introduced by the phantom in the NNPS evaluations.CNR in a PMMA phantom with Aluminium inclusion was evaluated on images acquired via the clinical unit as well as via the insilico replication.In the second cases images were manipulated in terms of spatial resolution and noise as described by the model.CNR between 1.00-mm thick Al inclusion and PMMA background resulted 7% lower in simulated images at 26 kV; such a difference increased to 16% and 12% at 29 and 31 kV, respectively.

Fig. 3 .
Fig. 3. Measured and simulated NNPS at (a) 26 kV, (b) 29 kV, and (c) 31 kV.Dashed lines show NNPS after R-filter and green continuous lines show the same curve after the addition of white Gaussian noise.The reported dose values represent the air kerma evaluated as K S and K M .

Fig. 4
reports a test on a contrast test object comprising 28 mm of PMMA on which two Al foils are lying.These two present a thickness of 1.00 and 0.025 mm.The test is performed at 26 kV and for 1.7 mGy of air kerma estimate at Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 4 .
Fig. 4. (a) Post-processed simulated image of a contrast object phantoms (PMMA thickness = 28 mm) and (b) corresponding image acquired via a clinical scanner at 26 kV.(c) Image profiles over the yellow line indicated in (a).Air kerma at the phantom surface = 2.10 mGy.the breast support in the absence of the phantom.The ROIs in Fig.4(a) and (b) show the simulated image after R-tuning and AWGN inclusion and measured image, respectively.Fig.4(c) reports profile along the yellow line represented in Fig.4(a), which crosses the Al-PMMA edge.It can be noted that both, noise background fluctuations and sharpness of the edge are comparable.TableIVcompares CNR between the 1-mm Al foil and PMMA block evaluated for simulated and acquired images on clinical scanner.The thickness of the PMMA block was 28 mm for 26-kV tube voltage and 55 mm for 29 and 31 kV tube voltages.The ratio between simulated (following the R-tuning and the AWGN inclusion) and measured CNR at 26 kV resulted 0.93 indicating a difference of 7% in the qualitative image index.In the cases of 29 and 31 kV such a ratio resulted 0.84 and 0.88, respectively.Similarly, we compared the image noise, computed as the standard deviation in defined ROIs, over a region included in the projections of the two Al foils as well as over a region behind the PMMA.V reports the between the pixel standard deviation for the simulated image phantoms the R-tuning and the inclusion of the AWGN) and the standard deviation in phantom images acquired on the clinical scanner.The lowest values for such a ratio can be observed for the ROI comprised in the PMMA background, ranges between 1.03 at 26 and 0.83 at 31 kV.These two values are 1.05 and 0.88 for an ROI over the 0.025-mm Al foil and increase up to 1.15 and 0.88 for pixel comprised over the projection of the 1-mm Al foil.
c 2023 The Authors.This work is licensed under a Creative Commons Attribution 4.0 License.
For more information, see https://creativecommons.org/licenses/by/4.0/Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I LINEAR
FIT COEFFICIENTS OF THE MEASURED AND SIMULATED DETECTOR RESPONSE CURVE.THE LAST COLUMN REPORTS THE R 2 FITTING PARAMETER FOR THE MEASURED CURVES Fig. 1.Measured and simulated MTF at 29 kV and the calculated R filter.Continuous lines show the Gaussian fit curves.

TABLE II ESTIMATED
σ S , σ M , AND σ R COEFFICIENTS

TABLE IV SIMULATED
AND MEASURED CNR BETWEEN 1.00-MM AL FOIL AND PMMA BACKGROUND AND RATIO

TABLE V RATIO
BETWEEN PIXEL STANDARD DEVIATION IN SIMULATED IMAGES AND IN IMAGES ACQUIRED ON CLINICAL SCANNER