Accuracy of the independent atom approximation in digital tomosynthesis Monte Carlo simulations

Monte Carlo (MC) codes serve as the gold standard simulation tool during design and optimisation of x-ray imaging systems. Such simulations often model Rayleigh scattering based on the Independent Atom Approximation Model (IAM). This model neglects the low range molecular interference (MI) effects of non-crystalline materials such as human tissues. Previous work has found discrepancies in the simulated images of planar x-ray images between IAM and MI models. However, insignificant differences were found for computed tomography (CT) reconstructions. In this work we present Geant4 MC simulations of a flat panel source digital tomosynthesis (DT) system for human extremities. Results show that with a 1:9 scatter to primary ratio (SPR) in the x-ray projections, the DT reconstructions are insensitive to the differences of the IAM and MI models. Therefore, MC codes that use the IAM model are sufficient for the study of DT systems. That is because DT algorithms have a larger effect on image quality than the few percent change in the noise due to a physical model and noise suppression methods make this change even less important. Dependency of this conclusion on SPR must be considered in other DT modalities where SPR might be larger.


Introduction
There has been a recently renewed clinical and commercial interest in digital tomosynthesis (DT) systems. It is driven by the lower patient dose, due to fewer projections compared to computed tomography (CT), along with increased concerns about cumulative dose from recurring medical diagnostic procedures. This interest is also complemented by the lower manufacturing and maintenance costs of DT compared to CT as well as by the portability of the former which could allow imaging at the patient's bedside for more personalised health care or deployment at the site of the accident for a true point of care imaging.
As a result, DT is under extensive research with a variety of tools, the gold standard of which being Monte Carlo (MC) simulations. Regardless of architectural differences, all MC codes employed must accurately model physics at the diagnostic energies of ∼10s keV. In most MC codes the default model for the simulation of photon coherent scattering, also known as elastic or Rayleigh scattering, is the Independent Atom Approximation Model (IAM). Its accuracy has been tested against more realistic models that include molecular interference (MI) effects caused by soft tissues (Tartari et al 2001), (Poludniowski et al 2009) and it was found that discrepancies propagate to the final image quality of 2D x-ray images, but changes in CT slices are minor. 2D x-ray images of chest, pelvis and skull had 5% discrepancies on scatter, mostly on anatomical edges but CT reconstruction algorithms were insensitive to this difference and produced slices with discrepancies less than 10 H (Poludniowski et al 2009).
That is because the IAM neglects MI effects of human tissues. The Rayleigh scattering differential cross section is calculated assuming tissues are composed of free atoms, like a gas. This leads to a maximum scattering probability at 0  scattering angle. However, human tissues and even liquids such as water have a low range order that causes MI effects and thus scattering is not peaked at 0  (4  for 60 keV photons in liquid water) (Narten and Levy 1971). Literature is rich in experimental data of human (Poletti et al 2001) and animal tissues (Peplow and Verghese 1998)  Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. demonstrating a total or multiple local peaks at scattering angles other than 0. A compilation of references to such data can be found in a recent work by Paternò (Paternò et al 2020).
In 2D x-ray imaging, these effects are direct on the images but in CT the reconstruction process minimises their effect. However, DT is an incomplete reconstruction problem due to the limited angle of projections used. As such, the conclusions about 2D x-ray imaging or CT might not be valid for DT. Therefore, a comparison of the IAM and MI models in the simulation of DT image stacks is necessary. To the authors' knowledge the literature lacks such a comparison.

Imaging geometry
There are multiple DT irradiation geometries available. Examples include systems with a stationary or mobile detector paired with a single scanning x-ray source, or systems with arrays of stationary x-ray sources paired with a stationary detector. The latter design spares the engineering challenges and image artefacts of movements, in either the system or the patient. However, regardless of the irradiation geometry, the DT reconstruction principles remain the same. Therefore, this study utilised the irradiation geometry of the orthopaedic DT system manufactured by Adaptix Ltd, shown in figure 1.
The Adaptix orthopaedic DT system includes a compact, stationary flat panel source at the top and a stationary flat panel detector on the bottom, parallel and centred to each other and at a source to image distance (SID) of 20 cm. The source is made of a silicon wafer on which 45 electron field emitters have been etched in a7 7 grid without the corners and with 1 cm pitch. The wafer acts as the cathode held at a 60 kV potential difference from the anode containing a transmissive x-ray target. The anode also acts as an x-ray filter and as a vacuum to air interface. A multi hole collimator sits parallel to and in-between the cathode and anode, with the holes at a lateral shift from the emitters. When the system is turned on, all electron beams are absorbed in the collimator, but magnetic coils can be used to steer individual electron beams through dedicated apertures in the collimator to interact with an x-ray target. Finally, downstream from the anode and in air, lies a printed circuit board (PCB) that controls the system.
The flat panel source array allows the user to irradiate the detector with each individual source sequentially and thus create x-ray projections from different angles. These 45 projections are then normalised on a pixel-by-pixel basis with a set of 45 air shots. Details of the normalisation method are given later. Finally, the normalised images are processed with a DT reconstruction algorithm (Soloviev et al 2020) to produce a stack of DT slices.
2.2. Monte Carlo simulation 2.2.1. The Monte Carlo codes The Adaptix DT system is simulated in Geant4 (Agostinelli et al 2003) using a flat panel detector model and an analytical description of the x-ray source array. The analytical source model is a fit from results of a FLUKA (Ferrari et al 2005) Monte Carlo simulation of the Adaptix DT x-ray source array (Primidis et al 2021). Both Geant4 and FLUKA use IAM by default. Geant4 was chosen for the rest of the simulation as FLUKA does not allow the necessary source code adaptations to implement MI effects, whereas a Geant4 MI extension has been published (Paternò et al 2018), (Paternò et al 2020). Therefore, it was deemed most efficient that the available analytical model is used to sample the x-ray photons and Geant4 is used to model their transport through the imaging subject and the detector, with and without MI effects. A range cut of 1 mm was used in Geant4 and the simulations run for 7 h in a high-performance computing facility.

Source array model
For the scope of this study, all 45 x-ray sources of the Adaptix DT system are considered equivalent and with their beam axis normal to the source plane. Therefore, the same analytical description is used to sample photons from each source, with only a lateral shift of the photon position required depending on the source that is simulated. The x-ray beam, as described by the analytical functions, is a rotationally symmetric gaussian cone beam, including a penumbra and with an energy spectrum between 20-60 keV. Finally, it is further described by a scalar parameter which expresses the average photon yield per mAs.
The photon yield per mAs is assumed identical for all sources. It had been simulated previously in FLUKA and Geant4 with agreeing results that were also comparable to experimental values from a prototype machine (Primidis et al 2021). Also, the mAs of each emitter is, by Adaptix's specifications, a constant value and the same for all emitters. So, the total number of photons emitted by each x-ray source is calculated by multiplying the simulated photon yield with the mAs per emitter. This results in 2.5 billion photons per x-ray image that are sampled in Geant4 using the analytical source model. Due to the source equivalence mentioned above, the same number of photons is used for all 45 images.

Phantoms
Two phantoms are used in this study. The first phantom is the ICRP110 reference adult male voxel phantom (International Commission on Radiological Protection 2009); specifically, the left hand starting above the wrist and ending after the fingertips. The 53 tissues of the voxel phantom have been grouped down to 5 based on their density as shown in table 1, a method inspired by a similar 5-material voxel phantom in the literature (Poludniowski et al 2009). The phantom's hand does not contain tissues 1 and 50 so they are not used, and all material compositions are taken from tables or assumptions in the literature (Paternò et al 2020). It must be noted that the original voxel dimensions 2.137´2.137´8 mm 3 were decreased to 2´2´8 mm , 3 shrinking the phantom from 87.617´89.754´200 mm 3 to 82´842 00 mm .
3 For the interests of this study, the reduction of the size is deemed unimportant. Also, this work is not about the image resolution but about the distribution of primary and scattered radiation. Therefore, the large voxel size is not of concern. Nevertheless, the complexity of the phantom is representative of the complexity of human tissue and that alone is satisfactory for the purpose of this study.
The second phantom is comprised of basic Geant4 shapes and resembles a severely broken extremity as shown in figure 2. The phantom is made of an elliptical water cylinder with length, minor axis and major axis equal to 10 cm, 5 cm and 7 cm, respectively. The water cylinder acts as the soft tissue inside of which there are two annular cylinders and two spheres made of bone. The outer diameter of the two cylinders is 2 cm and the inner diameter is 1 cm. Both cylinders are made of healthy bone and their insides are made of osteoporotic bone. The bone cylinders are diagonally cut at one end to mimic the fracture and the two bone spheres between them act as fragments of 6 mm and 2 mm diameter. The two bone compositions are shown in table 1 and they are those from literature (Paternò et al 2020) that best match experimental scattering data (Royle and Speller 1995) using a method suggested by Tartari (Tartari et al 2001) More details of that method are presented later.

Detector model
The flat panel detector is modelled as a plane with 100 μm pixel size and resolution of 1500´1000 and 1500´1500 pixels for the primitive and for the voxel phantom, respectively. Therefore, the detector is either15 10 cm 2 or15 15 cm .
2 Each pixel sums the number of the incident photons in it, weighted by the cosine of their angle of incidence. The detector plane is at a 5.6 mm distance from the phantoms to simulate the distance between the detector surface and the detector active layer. Scoring is done using a concrete class derived from the G4UserSteppingAction class.

Molecular interference effects in Geant4
To investigate the difference on the image caused by more realistic MI effects on diagnostic x-ray coherent scatter, the user source code is linked with a recently published extension that implements these effects in Geant4 (Paternò et al 2020). This extension is a modification to the Geant4 Penelope interaction models that originate from the PENELOPE code (Salvat et al 2019). It allows the user to either activate coherent scattering MI effects on a set of specific tissues and materials for which MI data are available or to define a tissue as a compound of water, fat, bone matrix and hydroxyapatite. The MI data of the compound are automatically calculated using the MI data of the four materials with a method described by Tartari (Tartari et al 2001). This method changes the coherent scattering differential cross sections but the extension by Paternò (Paternò et al 2020) also integrates them to further recalculate the total scattering cross section. Finally, the MI effects can be turned on and off by means of a single variable which offers a user-friendly way to replicate the same simulation with the IAM and with the MI coherent scattering models and compare the result.

Image processing
Prior to the image acquisition, 45 images without the phantoms in place are simulated with the same number of photons as would be simulated with the phantoms in place. These air shots are used to calculate the pixel-wise attenuation matrices, μ. These are defined as the natural logarithm of the ratio of the air shot, I , 0 divided by the phantom shot, I, according to the Beer-Lambert law, These 45 attenuation matrices are the input for the DT reconstruction algorithm.
μ is a real number. It is negative in pixels that collect more photons after the beam scatters in the phantom, as photons are redirected to pixels where they would not normally arrive from an air shot; < I I . 0 μ is positive in pixels that, due to scattering and absorption, receive fewer photons with the phantom than with the air shot; > I I . 0 Theroretically, μ can also be zero in pixels where the beam intensity has not changed with the inclusion of a phantom, so = I I , 0 but practically this never happens due to the stochastic nature of photon emission and interaction in air and tissues.
It is worth noting that in the clinic, air shots are preloaded on the system or taken during the various standard calibration procedures, without interfering with the patient's scan.
The DT reconstruction algorithm accepts various data formats, but this study uses uncompressed PNG images. Therefore, the 2D attenuation matrices with negative and positive float numerical data have to be converted into 16-bit grey scale PNG images. At the time of this work, the actual imaging system extracted the 2D attenuation images in 16-bit grey scale PNG format and then reconstructed them, so we decided to maintain this practice in the simulation. To do so, a user Python script scans all 45 attenuation matrices and creates a single histogram with the frequency of their pixel values. This is done separately for the IAM and the MI matrices. Then, the lower 1% and higher 99% bins of the two histograms are used to define the grayscale's white and black ends. But, to have the same grayscale for both IAM and MI, if between the IAM and MI histograms, the minimum 1% mark is a positive number, this number will be the black end of the grayscale for both models otherwise zero is used. For the white end, the 99% margin caused brightness saturation on the projections from the outer sources, so we decided to use1.2 99% margin as the value for the white end. The 99% margin used is the maximum between that of IAM and MI histograms. With this method, all 45 attenuation matrices created with both scattering models are converted to grey images using the same grey scale. It is only between the two phantoms that the grey scales are different.

Digital tomosynthesis algorithm
The reconstruction algorithm used in this study is a modification of the filter of back projections (Soloviev et al 2020). In contrast to the filtered back projection algorithm, the back projection comes first and the filtering follows. It differs from the original filter of back projections suggested by Bates and Peters (Bates and Peters 1971) in that respect that the ramp filter is applied over planes (slices) parallel to the flat panel detector, while the original one performs filtering over the planes perpendicular to the detector with parallel rays illumination. In this case the ray tracing technique, which is usually employed for computing back projections, is replaced with an equivalent image mapping procedure. Thus, attenuation images are mapped at a particular height above the detector as described in Soloviev et al (2020) and their average is fed to the ramp filter. The cone beam geometry is taken into account by the image mapping procedure.
This approach is well suited to the DT geometry. It does not require three-dimensional grid or mesh allocation for storing reconstruction results. Moreover, it performs a user interactive, slice-by-slice reconstruction where each slice position can be chosen at runtime. This allows the display of reconstructed slices on the fly as soon as they are available. The reconstruction slice order can be arbitrary, thus, allowing the most relevant slices to be reconstructed and displayed first.
Contrary to CT, DT reconstructs slices from a set of projection images taken over a limited angle range. Therefore, the inverse problem suffers from incompleteness of data. This manifests itself in the appearance of negative numbers and poor image quality. To address this issue, the reconstruction problem must be regularized (Engl et al 2000).
When the reconstruction is performed in the Fourier domain, an addition of a small number to the ramp filter |k| (truncated at the maximum allowed value of the wave number) has a regularization effect. Its value is about a few sizes of the reciprocal lattice cell. Alternatively, filtering can be done by computing a convolution of the average of mapped images with the spatial analogue of the ramp filter over the finite window. Reduction of the size of the convolution window has a regularization effect (Soloviev et al 2020).
The reconstruction algorithm also allows the user to choose a resolution for the DT slices different to that of the input 2D x-ray images. In this study, the resolution of the slices is reduced by a factor of 2 to speed up the reconstruction.

Comparison of results
First a Python script reads the raw projection matrices generated by the simulation and calculates the pixel-bypixel relative difference caused by the two scattering models on the distribution of primary and scatter radiation. This is done to compare these results with literature on the differences found in planar x-ray images (Poludniowski et al 2009). Then, another Python script calculates the pixel-by-pixel relative percentage difference between the two scattering models for each DT slice. For completion, the same comparison is done on the grey scale attenuation images.

Results
The difference between IAM and MI in the spatial density of primary and scattered photons for one of the 45 projections of the voxel phantom is shown in figure 3. The figure presents both single and multiple scatter components of both coherent and incoherent scattering. The difference between the IAM and the MI models in scatter density is within 5% and is most profound at the edges, which agrees with literature (Poludniowski et al 2009). The structure of the phantom is visible in the difference of the single coherent scatter, but it is lost in all other scatter components. Also, the highest relative differences are found in the coherent scatter, both single and multiple, and the recalculated total scattering cross section has induced a small difference in the primary radiation too. However, the difference in the total radiation remains low. This can be explained by the fact that with both models the primary (not scattered) radiation forms 90% of the detected radiation, the single coherent scatter forms less than 3% and the multiple coherent less than 0.2%. Finally, the projection in figure 3 is from a source on the top right (A2, figure 2), so the higher difference at the bottom left of the image, mostly visible in the primary radiation, is due to the high noise of the penumbra.
The same information for the primitive phantom is shown in figure 4. As above, results agree with the literature regarding the difference in the scatter density. The difference in the primary radiation is higher than before, however the difference in the total radiation remains close to 3%. In this case the primary radiation with both models makes up 90%, the single coherent scatter is less than 3% and the multiple coherent scatter forms less than 0.2%, which explains the small difference in the total radiation. Moreover, this detector is 5 cm shorter on the vertical axis so the penumbra on the bottom left is not as visible as before. However, the larger height of this phantom in the beam direction causes cropping on the bottom left of the difference in the scattered radiation.
Histograms with the raw pixel values of all attenuation matrices after using the IAM and MI scattering models are presented in figures 5 and 6 for the voxel and for the primitive phantom, respectively. Both models produce a practically identical distribution of pixel values apart from a small difference after 1.5 for the primitive phantom. In both figures, the peak at zero is for pixels not behind the phantom and the negative values are caused by scattered radiation. In figure 6, the clearly defined peaks of the primitive phantom at 1 and 1.9 are from pixels behind the water cylinder and behind the bone cylinders, respectively. On the contrary, the higher complexity of the voxel phantom produces convolved peaks.
From figure 5, the voxel phantom attenuation matrices resulting from both scattering models are converted to 16-bit grey scale PNG images using a linear scale with white and black ends at 0 and 2.43, respectively. Similarly, from figure 6, the attenuation matrices of the primitive phantom are converted with ends at 0 and 2.52.
Examples of PNG attenuation images of both phantoms and with both Rayleigh scattering models  are shown in figures 7 and 8. All images are illustrated with their previously and respectively allocated grey scales. As intuitively expected from the radiation distributions in figures 3 and 4, the IAM and the MI scattering models do not cause a visually obvious difference in the x-ray attenuation images. Calculating this difference reveals a discrepancy less than 4%.
Moreover, the SNR of a200 200 pixel region was measured just above the wrist and at the fingers in figure 7. Figures 7(a) and (b) have SNR equal to 11.30 and 11.35 above the wrist and 8.21 and 8.22 at the fingers. Figures 7(d) and ( Similarly, the SNR on the top left side of the water and in the centre of the left bone was measured in figure 8 using80 80 pixel regions. Figures 8(a) and (b) have SNR equal to 11.35 and 11.29 in the water and 19.73 and 19.65 in the bone region. Figures 8(d) and (   All SNR results are shown in figure 9. Each projection is named after the source that generated it with a letter-number combination, with letters denoting horizontal and numbers vertical position of the source relative to the source plane as illustrated in figure 2. It is clear from figure 9 that the two models do not change the SNR for any of the images. The SNR difference among the sources is caused by the change in projection angle and the manual positioning of the regions during the calculation.
Finally, slices from the reconstructed volumes of both phantoms are shown in figure 10 and figure 11. Again, there is no visually obvious difference. By calculation we find this difference lower than 4% and only dark areas show values above that, an obvious result of numerical precision effects. Also, as shown in figure 12, the SNR of various areas on slices of the two phantoms are affected to less than 2% by adding or removing MI effects. These areas are on structures that are in focus in the selected slices. Each slice brings different structures in focus therefore intra-slice comparison is meaningless.

Discussion
In these simulations we decided to use a pixel size of 100 μm, a low flux per projection and an arbitrary way (a) using IAM and (b) using MI are images generated by the A2 source, (c) relative difference of (a), (b). (d) using IAM and (e) using MI are images generated by the G6 source, (f) relative difference of (d)-(e). Dark spot in (a) and (b) is due to the alignment of the x-ray beam with a gap within the curled fingers. Images are binned with a 4´4 kernel prior to calculation of their difference. Differences appear minor and only black areas cause large deviations due to operations with very small numbers. For true colours see the online version. Figure 8. Attenuation images of the primitive phantom and their absolute percentage difference caused by using either IAM or MI. (a) using IAM and (b) using MI are images generated by the A2 source, (c) relative difference of (a)-(b). (d) using IAM and (e) using MI are images generated by the G6 source, (f) relative difference of (d)-(e). Images are binned with a4 4 kernel prior to calculation of their difference. Differences appear minor and only black areas cause large deviations due to operations with very small numbers. For true colours see the online version. of converting photon flux to pixel value for the following reasons.
There is a minimum distance of 5.6 mm between the phantoms and the active layer of the detector. With a Rayleigh scattering peak at 4°for 60 keV photons in liquid water (Narten and Levy 1971), these scattered photons would be expected at a minimum offset of 392 μm. Also, bone has multiple peaks in a wide range of momentum transfer values that includes that of liquid water (Paternò et al 2020). Moreover, the mean energy of the spectrum is well below 60 keV (Primidis et al 2021) so scattering will be at larger angles causing lateral offsets that are wider than 392 μm. Therefore, 100 μm pixels are small enough to resolve these photons when binned by 4 (effective size 400 μm).
The photon yield was simulated in FLUKA and Geant4 and it was compared to experiments in our laboratory. The simulations generate the same estimation for the photon flux which is about 4 times the experimental values of a prototype machine (Primidis et al 2021). We have identified collimator misalignments in the machine as the reason for that and we are optimising the manufacturing to minimise it. Therefore, there is enough confidence that the expected photon flux is estimated well. Nevertheless, the projections are noisy due to the inherently low flux of the system that stems from the mAs specifications and the simulated photon yield per mAs. This could add doubts regarding the validity of our results in a higher flux and thus less noisy procedure. However, in figure 3(a) the colour map is uniform close to 0% difference while in figure 4(a) the cold area of the bone cylinders in the total radiation is below 5%. This simple geometry is unrealistic therefore based on figure 3 with the complex and more realistic voxel phantom, increasing the photon flux should not change the conclusions.
As was mentioned above the reconstruction problem is incomplete. In general, the regularization affects the SNR and the image contrast (Engl et al 2000). If we consider the reconstruction in the Fourier domain, the value of regularization parameter smaller than optimal usually provides reconstructed images of higher contrast but with lower SNR. For example, a small value amplifies the Figure 9. SNR of various areas in the attenuation images with both scattering models. Images are named after the source that produces them. The effect on the SNR due to a change in the scattering model is smaller than the symbol size. Changes on an image by image basis are due to the different projection angles and the manual way of placing the area in which the SNR is calculated. Figure 10. (a), (d) two tomosynthesis reconstructed slices at different depths in the voxel phantom using the IAM model. (b), (e) the same slices using the MI model, (c) absolute percentage difference of (a)-(b), (f) absolute percentage difference of (d)-(e). Voxels that contain the phantom have a discrepancy mostly below 4% while dark voxels have discrepancies above that. A zone of zero difference surrounding the phantom comes from the reconstruction algorithm. The same applies for the rectangular frame of zero difference on the outermost regions of the slices. The reconstruction algorithm reduces the image dimensions by half to speed up the process and resulting images are binned with a2 2kernel prior to calculation of their difference. For true colours see the online version. noise, which is intrinsically present in attenuation images, whilst a larger than the optimal value produces smoother images, but with a stronger background component. In many cases only an approximation for the optimal value can be found. The same is true for the width of the window when the reconstruction is performed spatially by convolution. Nevertheless, the image quality is usually acceptable when the SNR is greater than 10, which implies that the noise level in reconstructed images can even reach 10%.
Therefore, variation of the noise level below a few per cent due to approximations made in a physical model are not important in terms of the reconstruction. In addition, computational techniques for noise suppression make small variations in the noise level caused by the choice of a physical model to be even less significant.
Although arbitrary, it is not uncommon in the literature to use cosine weighted photon fluence as the signal for simulated images and the decision to use the 99% margin in figure 5 and figure 6 multiplied by 1.2 as the white end avoids brightness saturation in the projections. The shade of these projections as shown in figures 7 and 8 are similar to what we see in the laboratory and what we have seen in the literature. As such, we think our methodology is reasonable and straightforward to understand and implement. Figure 11. (a), (d) two tomosynthesis reconstructed slices at different depths in the primitive phantom using the IAM model, (b), (e) the same slices using the MI model, (c) absolute percentage difference of (a)-(b), (f) absolute percentage difference of (d)-(e). Voxels that contain the phantom have a discrepancy mostly below 4% while dark voxels have discrepancies above that. A zone of zero difference surrounding the phantom comes from the reconstruction algorithm. The same applies for the rectangular frame of zero difference on the outermost regions of the slices. The reconstruction algorithm reduces the image dimensions by half to speed up the process and resulting images are binned with a2 2kernel prior to calculation of their difference. For true colours see the online version. Figure 12. SNR of various areas on the slices from figures 10 and 11. The first four column pairs are from the voxel phantom and the remaining are from the primitive one. The areas on the voxel phantom slices are 50-pixel wide squares, placed on two positions on the left fingers (LF1, LF2) and two positions on the right finger (RF1, RF2). The areas on the primitive phantom are circles with 40-pixel diameter placed in the water above the left bone (WA), in the centre of the left bone (LB) and in the large bone sphere (BS). The selection of all regions is based on the fact that these are in focus. SNR remains unaffected within 2% by the inclusion of MI effects. However, digital breast tomosynthesis (DBT) uses much lower energies than this work. At those energies and with vast ranges of breast thickness and composition, the SPR can be different to the 1:9 of this work. More importantly, the portion of the single coherent scatter could be higher and as shown in figure 3 and  figure 4, it holds spatial information about the phantom. So, there could be situations where there is a larger component of single coherent scatter which carries structural information that is inaccurately modelled and with this inaccuracy being propagated to the image quality even after a DT reconstruction. This could also deteriorate scatter correction methodologies based on MC scatter estimations. Therefore, our conclusions cannot be intuitively extrapolated to DBT. Thankfully, this is trivially addressed on an ad hoc basis by simply identifying combinations of energy spectrum, breast thickness and composition that have higher ratios of single coherent scatter to total radiation. If a manufacturer or medical physicist finds their equipment to be in that situation, they can test whether the IAM produces different image estimations compared to the MI model. Otherwise, the random nature of the other scatter components diffuses the information carried by the MI corrected coherent scatter making the corrections obsolete and the IAM model adequate.

Conclusion
The recent heightened academic and commercial interest in developing DT imaging systems has brought attention to MC radiation transport literature where the accuracy of commonly used MC photon scattering models such as the IAM had been challenged. Conclusions so far had been drawn for planar x-ray imaging and for CT, the former requiring more accurate MI effects to be included in the models and the latter being sufficient with the IAM. With this report we have complemented the above conclusions by proving that DT algorithms are also insensitive to physical effects like their CT counterparts. Therefore, multipurpose MC codes that use the IAM are adequate for the simulation of DT imaging systems. This conclusion is limited to DT with low SPR. However, it is trivial to estimate the SPR of different DT scenarios such as DBT. If SPR is similar to or lower than this report, then the IAM should be sufficient. Otherwise, the Geant4 extension by Paternò (Paternò et al 2020) can give estimates of model discrepancies on an ad hoc basis.