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

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$^\circ$ 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$^\circ$ (4$^\circ$ 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) demonstrating a total or multiple local peaks at scattering angles other than 0$^\circ .$ 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.

2.1. 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.

Zoom In Zoom Out Reset image size

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 a $7\times 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.

2.2.2. 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.

2.2.3. 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 $\times$ 2.137 $\times$ 8 ${{\rm{mm}}}^{3}$ were decreased to 2 $\times$ 2 $\times$ 8 ${{\rm{mm}}}^{3},$ shrinking the phantom from 87.617 $\times$ 89.754 $\times$ 200 ${{\rm{mm}}}^{3}$ to 82 $\times$ 84 $\times$ 200 ${{\rm{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.

Table 1. Materials used for the two phantoms and their percent decomposition in a compound of four materials (Tartari et al 2001) as matched or compared by Paternò (Paternò et al 2020) to experimental data. The 53 tissues of the ICRP110 reference adult male phantom are grouped into 5 based on their density. The left hand of the voxel phantom does not include tissues 1 and 50 so they are not used. Tissue 53 is air.

	Fat	Water	Bone matrix	Hydroxyapatite	References
Geant4 phantom
Phantom components
Soft tissue		100
Bone and bone fragments	36	15	13	36	Normal bone (Royle and Speller 1995)
Bone marrow	55	25	5	15	Osteoporotic bone (Royle and Speller 1995)
Left hand of ICRP110 reference adult male voxel phantom
Tissue number
22–25,49	83	17			Adipose tissue (Poletti et al 2001)
17,18,20,21,27–48,51,52	15	73	12		Muscle (King et al 2011)
3–16,19,26	80	20			Yellow marrow assumption matching femoral head (White, 1978)
2	22	22	19	37	Bone (King et al 2011)

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.

Zoom In Zoom Out Reset image size

2.2.4. Detector model

The flat panel detector is modelled as a plane with 100 μm pixel size and resolution of 1500 $\times$ 1000 and 1500 $\times$ 1500 pixels for the primitive and for the voxel phantom, respectively. Therefore, the detector is either $15\times 10\,{{\rm{cm}}}^{2}$ or $15\times 15\,{{\rm{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.

2.2.5. 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.

2.3. Image processing

2.4. Digital tomosynthesis algorithm

2.5. Comparison of 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Moreover, the SNR of a $200\times 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 (e) have respectively 10.82 and 10.83 above the wrist and 6.50 and 6.50 at the fingers. Among the 45 projections, using the IAM(MI) model, the average SNR of $200\times 200$ pixel regions just above the wrist is 12.76 (12.78) with standard deviation of 1.24 (1.24). The regions at the fingers have an average SNR of 8.53 (8.55) with standard deviation of 0.75 (0.75) among all 45 projections. Each region was manually shifted per image to follow the change of projection angle but regions in each image are identical between the models.

Similarly, the SNR on the top left side of the water and in the centre of the left bone was measured in figure 8 using $80\times 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 (e) have respectively 12.62 and 12.53 in the water and 21.25 and 21.18 in the bone region. Among the 45 projections, using the IAM(MI) model, the SNR of $80\times 80$ pixel regions on the top left side of the water have an average value of 12.74 (12.67) with standard deviation of 0.57 (0.55) and $80\times 80$ pixel regions in the centre of the left bone have an average SNR of 20.87 (20.80) with a standard deviation of 0.67 (0.65). These regions are shifted per projection like describe previously.

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In these simulations we decided to use a pixel size of 100 μm, a low flux per projection and an arbitrary way 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 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.

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.

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.

This work was undertaken on Barkla, part of the High Performance Computing facilities at the University of Liverpool, UK. Thomas Primidis is funded by the Accelerators for Security, Healthcare and Environment program launched by the United Kingdom Research and Innovation Science and Technology Facilities Council with reference number ST/R002142/1. Finally, the authors would like to thank Dr Joseph Wolfenden for the constructive comments during the writing of this report.

The data that support the findings of this study are available upon reasonable request from the authors.