Characterization of breast tissues in density and effective atomic number basis via spectral X-ray computed tomography

Tissue differentiation in mammography is challenging due to the overlap of tissues with similar linear attenuation coefficients μ in the diagnostic energy range. A comparison of the attenuation of glandular tissue and close-by malignant tumors found that the average attenuation values are almost identical even in the low-energy range (Chen et al 2010, Fredenberg et al 2018). Computed tomography (CT) solves tissue overlapping issues and provides additional diagnostic information based on tissue morphology in 3D volumes. On the other hand, breast-dedicated CT uses higher X-ray energies than mammography due to the imaging of uncompressed breasts, thus reducing attenuation contrast between similar tissues. Tissue separation can be significantly improved by acquiring data at multiple energy levels using spectral CT imaging (Alvarez and Macovski 1976). New generations of whole-body spectral CT systems obtain source-generated spectral information from two (or more) scans at different energy levels (dual-energy), or detector-generated spectral response with layered and photon-counting (PC) detectors. PC detectors are particularly interesting due to the availability of multiple energy thresholds (up to 12 (Danielsson et al 2021)), flat spectral response, and low imaging noise. The first breast CT scanners with PC detectors are being introduced in clinics (Kalender et al 2017, Berger et al 2019, Schmidt et al 2022, Zellweger et al 2022) and synchrotron facilities (Longo et al 2016, 2019), where spectral separation can be obtained using multiple monochromatic x-ray beams. Despite the available technology for spectral imaging in both systems, the benefits of such an approach in breast CT imaging have not been explored yet.

Given two (or more) scans at different energies, the content of each voxel can be described as a linear combination of the attenuation of two basis materials in the process known as material decomposition (Alvarez and Macovski 1976). Decomposition to two basis materials is the most common approach because two physical effects contribute to image formation in the diagnostic energy range: the photoelectric effect and Compton scattering. The photoelectric-Compton basis itself is very convenient because of the well-defined dependencies of these effects on density and atomic number. However, such basis spans an infinite range of physical materials, given that a 'purely Compton' material would be the element in the limit Z → 0, while a 'purely photoelectric' one would correspond to Z → ∞ . Such a broad range would degenerate the separation of similar materials, such as soft tissues. To address this issue, the use of a pair of physical (or even virtual) materials that span the range of materials of interest has been shown to significantly reduce decomposition uncertainty (Champley et al 2019). A popular choice of physical materials for decomposing biological tissues are polymethyl-methacrylate (PMMA) and aluminum (Al) basis pair (Lehmann et al 1981), also used in this study. Estimated weights of the linear combination of two basis materials carry physical information but have no particular comprehensive meaning except when one of the basis materials is the same as the material of interest (e.g. iodine quantification or calcium scoring). However, using the weights and the known energy dependence of basis attenuation coefficients virtual monochromatic images (VMIs) can be extrapolated at an arbitrary energy value within the diagnostic energy range (McCollough et al 2015). VMIs have simple interpretation and clinical case studies show that VMIs are diagnostically valuable for certain tasks (Albrecht et al 2019), and they are a standard output in the first clinical whole-body PC-CT (Naeotom Alpha, Siemens) (Rajendran et al 2021). Decomposition to ρ and Z_eff is a step further toward extracting the physical information about tissue composition from spectral measurements. While material decomposition describes the attenuation of one material as a linear combination of the attenuations of other known materials, ρ/Z_eff decomposition exploits the underlying physics of x-rays interaction with materials to extract information about their physical and chemical properties. Several works have been published exploring this approach (Heismann et al 2003, Torikoshi et al 2003, Szczykutowicz et al 2011, Azevedo et al 2016, Busi et al 2019, Champley et al 2019, Vrbaski et al 2022). Density and Z_eff are intuitive units that can be easily interpreted and correspond to physical properties that can be measured by other techniques, allowing for easy comparison. Previous work showed the benefits of ρ/Z_eff for tissue differentiation (Pascart et al 2019), radiotherapy planning (Hudobivnik et al 2016), and interventional radiology (Liu et al 2023).

In the present paper, a spectral study was specifically designed to characterize breast CT images in terms of material density ρ and the effective atomic number Z_eff. While ρ/Z_eff decomposition has been implemented in some clinical CT scanners (Rajiah et al 2020), no studies specifically focused on breast tissue characterization were performed. In section 2.1, a theoretical approach to computing ρ and Z_eff values from physical material decomposition instead of the photoelectric-Compton decomposition was developed to improve the accuracy of results. Using synchrotron x-ray beams at several energies and a high-resolution PC detector described in section 2.2, CT scans of a dedicated phantom and 5 mastectomy samples were acquired. Due to the superior imaging quality and spectral separation available with synchrotron setup, the work investigates the feasibility and the potential diagnostic benefit of ρ/Z_eff decomposition in breast CT imaging. This information is of interest since PC detectors matured enough to operate in clinical conditions with the ability to obtain better spectral separation than other approaches mentioned earlier (Danielsson et al 2021). The introduced concept of effective atomic number was further investigated in section 2.3 and several approaches to the computation of Z_eff were analyzed in appendix C. Section 2.4 describes tissue preparation and section 2.5 provides a detailed explanation of the practical realization of the analysis with careful consideration of the well-known problem of the decomposition noise (Dong et al 2014, Niu et al 2014, Mechlem et al 2016). Finally, the data analysis approach is given in section 2.6. We also showed that virtual monochromatic μ values at desired energy can be computed using the known dependence of x-ray attenuation on density and atomic number.

2.1. Theoretical model

In the material decomposition approach, the attenuation coefficient μ of a given material is expressed as a linear combination of the (known) attenuation coefficients of a pair of basis materials, here labeled μ₁ and μ₂

$$\begin{eqnarray}&&\mu ={x}_{1}\,{\mu }_{1}+{x}_{2}\,{\mu }_{2},\end{eqnarray} \tag{ 1 }$$

where x₁ and x₂ are the coordinates of the material in the reference frame identified by the selected basis. In the first step of our method, we used equation (1) to determine the coefficients x₁ and x₂, which represent the (energy-independent) relative concentrations of each basis material. Because more than two spectral scans were available for each sample (Piai et al 2019), instead of using traditional matrix inversion (Zhang et al 2019), x₁ and x₂ coefficients were calculated by using a least-square fit of the form:

$$\begin{eqnarray}&&\displaystyle \sum _{i}{\left(\mu ({E}_{i})-({x}_{1}\,{\mu }_{1}({E}_{i})+{x}_{2}\,{\mu }_{2}({E}_{i}))\right)}^{2}.\end{eqnarray} \tag{ 2 }$$

The fit procedure consists of a voxel-by-voxel minimization of the sum over the energies E_i of the squared residuals. From such a point of view, including more images of different energies adds further points to the plane, which is expected to increase the accuracy of the decomposition.

The attenuation of X-rays is influenced by several physical effects, which are dependent on the energy and material properties of the attenuator. In the range of photon energies useful for medical CT imaging, the attenuation coefficient μ of given material of density ρ, the atomic number Z, and atomic mass A are approximated as the sum of the two contributions

$$\begin{eqnarray}&&\mu (\rho ,Z,A,E)=\displaystyle \frac{\rho \,Z}{A}\left[{K}_{1}\,{Z}^{n-1}\,{f}_{P}(E)+{K}_{2}\,{f}_{C}(E)\right],\end{eqnarray} \tag{ 3 }$$

where the functions f_P (E) and f_C (E) encode the energy dependencies of the photoelectric and Compton effects, respectively, and K₁, K₂ and n are constants. Using equation (3) to describe basis materials in equation (1) and assuming A = 2Z, as it is almost true for any chemical element with Z ≲ 20, the linear combination coefficients of equation (1) read

$$\begin{eqnarray}&&{x}_{1}=\displaystyle \frac{\rho \,{Z}_{1}\left({Z}^{n}{Z}_{2}-Z\,{Z}_{2}^{n}\right)}{{\rho }_{1}Z\left({Z}_{1}^{n}{Z}_{2}-{Z}_{1}{Z}_{2}^{n}\right)},\end{eqnarray} \tag{ 4a }$$

$$\begin{eqnarray}&&{x}_{2}=\displaystyle \frac{\rho \,{Z}_{2}\left(Z\,{Z}_{1}^{n}-{Z}^{n}{Z}_{1}\right)}{{\rho }_{2}Z\left({Z}_{1}^{n}{Z}_{2}-{Z}_{1}{Z}_{2}^{n}\right)}\end{eqnarray} \tag{ 4b }$$

ρ_j and Z_j being the density and the atomic number of the jth basis material. It is straightforward to notice that equations (4a ) and (4b ) depend on both the density and the atomic number of the material. In order to decouple the two dependencies, it is necessary to rotate the reference frame (by an angle $\phi =\arctan \,\left({\rho }_{2}/{\rho }_{1}\right)$) and then rescale the second coordinate dividing it by the first one. The resulting coordinates, which read

$$\begin{eqnarray}&&{x}_{\rho }=\displaystyle \frac{\rho }{\sqrt{{\rho }_{1}^{2}+{\rho }_{2}^{2}}},\end{eqnarray} \tag{ 5a }$$

$$\begin{eqnarray}&&{x}_{Z}=\displaystyle \frac{{Z}^{{\ell }}\left({\rho }_{1}^{2}+{\rho }_{2}^{2}\right)-\left({Z}_{1}^{{\ell }}{\rho }_{1}^{2}+{Z}_{2}^{{\ell }}{\rho }_{2}^{2}\right)}{({Z}_{2}^{{\ell }}-{Z}_{1}^{{\ell }}){\rho }_{1}{\rho }_{2}}\end{eqnarray} \tag{ 5b }$$

(with ℓ = n − 1) are labeled x_ρ and x_Z because of their exclusive dependencies on the variables mentioned in the subscripts. Equations (5a ) and (5b ) represent the expected relationships between x_ρ and ρ, and between x_Z and Z, for a given choice of basis materials. Such expressions are used in the calibration procedure described in section 3.1, where the (known) densities and atomic numbers of reference materials will be fitted against the corresponding (measured) values of x_ρ and x_Z using the functional forms

$$\begin{eqnarray}&&{x}_{\rho }(\rho )=\kappa \,\rho ,\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{x}_{Z}(Z)=p\,{Z}^{\lambda }+q\end{eqnarray} \tag{ 6b }$$

stemming directly from equations (5a ) and (5b ) with the (physically motivated) coefficients replaced by the effective fit parameters κ, λ, p, and q. The resulting calibrated relationships map the measured x_ρ and x_Z of any imaged material onto its actual values of density and effective atomic number.

2.2. Scan setup

2.3. Calibration phantom and effective atomic number

2.4. Breast mastectomy samples

In addition to phantom scanning, post-mastectomy breast tissue images were analyzed in a retrospective study. The analyzed surgical samples (N = 5) were fixed in formalin, sealed in a vacuum bag, and conserved at 4 ^◦C. The preliminary analysis of the same data was published (Piai et al 2019) as a feasibility study of the synchrotron breast CT approach. All the procedures adopted in this work followed Directive 2004/23/EC of the European Parliament and of the Council of 31 March 2004 on setting standards of quality and safety for the donation, procurement, testing, processing, preservation, storage, and distribution of human tissues. In the present work, we further processed the data to extract ρ and Z_eff of breast tissues. Tomographic reconstructions of selected samples are given in figure 1. They all contained adipose, fibro-glandular, and tumorous tissue, but only in sample 4 existed a region of glandular tissue clearly separated from the fibrous. It also contained calcification regions that were not evaluated in this study. All samples contained some type of malignant tissue: samples 1, 2, and 3 contain infiltrating ductal carcinoma, sample 4 contains infiltrating ductal carcinoma with a core of desmoplastic tissue, and sample 5 contains vastly differentiated infiltrating ductal carcinoma. Samples 1, 2, 4, and 5 also contained portions of the skin.

Zoom In Zoom Out Reset image size

The mean glandular dose of 5 mGy delivered per scan was computed according to a dedicated Geant4 Monte Carlo simulation (Fedon et al 2015, Mettivier et al 2016). The phantom scans were acquired at 25, 28, 32, and 35 keV, and mastectomy samples at several monochromatic beam energy levels in the range of 24–38 keV, from which three scans were selected for the spectral analysis.

2.5. Practical implementation

The theoretical model described in section 2.1 was implemented with Python (Python 3.10.0) and an interactive delineation tool was developed to select an arbitrary region of interest within a reconstructed spectral data set (Matplotlib, Python 3.10.0) containing a single tissue type. In the first step of the process (framed with a dashed line in figure 2), the voxel-to-voxel material decomposition method using the PMMA-Al basis is applied to the selected regions (figure 2(a)), resulting in PMMA and Al maps for each ROI (figure 2(b)). In the next step, the material maps are combined into 2D histograms, one for each tissue type (figure 2(c)). In literature, this approach to material visualization is often referred to as CT fingerprinting (Rajiah et al 2020). Due to the presence of noise, the obtained histograms are blurred and elongated.

Zoom In Zoom Out Reset image size

The correct decomposition coefficients x₁ and x₂ are considered to be the centers of obtained distributions. They are extracted using a 2D Gaussian fit method of the form:

$$\begin{eqnarray}&&G({x}_{1},{x}_{2})=A\,{e}^{\left(-\tfrac{{u}^{2}}{2\,{S}_{u}^{2}}+\tfrac{{v}^{2}}{2\,{S}_{v}^{2}}\right)}\end{eqnarray} \tag{ 7 }$$

A is the intensity of the peak and S_u , S_v the spreads (that is, the standard deviations of the associated distributions) along the major and minor axes, respectively. The peak coordinates ${\bar{x}}_{1}$ and ${\bar{x}}_{2}$ relative to (x₁, x₂)-plane of the histogram are contained in the quantities u and v

$$\begin{eqnarray}&&u=({x}_{1}-{\bar{x}}_{1})\cos \theta +({x}_{2}-{\bar{x}}_{2})\,\sin \theta ,\end{eqnarray} \tag{ 8a }$$

$$\begin{eqnarray}&&v=({x}_{2}-{\bar{x}}_{2})\,\cos \theta -({x}_{1}-{\bar{x}}_{1})\,\sin \theta \end{eqnarray} \tag{ 8b }$$

θ is the tilt of the major axis of the spot with respect to the horizontal direction. After the fitting, the peak coordinates ${\bar{x}}_{1}$ and ${\bar{x}}_{2}$ found according to equations (8a ) and (8b ) are used in equations (4a ) and (4b ) to compute the final output of the procedure: the single values for densities and effective atomic numbers of selected tissues (figure 2(d)). The measurement uncertainty, estimated as the standard error on the center of fitted 2D Gaussian distribution, is propagated through all mathematical transformation steps and measurement calibration. For the sake of conciseness, the complete analysis is given in appendix A.

2.6. Data analysis

In addition to the experimental data obtained using the calibration phantom, the procedure described in the theoretical model (section 2.1) was also applied to published μ ('true') values (Schoonjans et al 2011) of phantom material inserts at the energy levels used in the experiment. The true data points are used to evaluate the accuracy of the model by calculating the percentage error between the experimental and the ground truth data points, as given in equation (9)

$$\begin{eqnarray} &&\% \,{error}=\displaystyle \frac{| \mathrm{ground}\,\mathrm{truth}-\mathrm{measured}| }{\mathrm{ground}\,\mathrm{truth}}\times \,100.\end{eqnarray} \tag{ 9 }$$

The segmentation of the breast CT reconstructions was performed by an experienced radiologist, who selected ROIs containing adipose, fibro-glandular, glandular, skin, and tumorous tissue already knowing the mastectomy content from specimens sampled for the histological examination. Considering magnification, slices were reconstructed with 0.057 × 0.057 mm pixel size and slice thickness of 0.049 mm. The segmentation was done slice by slice to extract values from a 3D volume, using at least 3 pixel thick (0.17 mm) margins to avoid partial volume effects. For each sample and tissue type, density and Z_eff values and their uncertainties were estimated respectively as the mean and standard deviation evaluated over 10 consecutive CT slices. To estimate the discrimination power of ρ/Z_eff and x_PMMA and x_Al decomposition, hence diagnostic potential, a mean Silhouette score (MSS) using Euclidean distance as a metric was computed on mean values of all tissue types collected from the mastectomy samples (Rousseeuw 1987). The MSS is a tool to quantify how similar are samples within the same cluster and how separated they are from other clusters, defined as:

$$\begin{eqnarray}&&\mathrm{MSS}=\frac{1}{n}\displaystyle \sum _{k=0}^{n}\frac{{\mathrm{NC}}_{k}-{\mathrm{IC}}_{k}}{max({\mathrm{NC}}_{k},{\mathrm{IC}}_{k})},\end{eqnarray} \tag{ 10 }$$

where IC is the mean intra-cluster distance and NC is the mean nearest-cluster distance for each sample k, and n equals the number of selected tissue types × number of selected mastectomy samples. Negative values of MSS indicate cluster overlap, values around 0 signal that samples are on or close to the boundary between clusters, and positive values up to 1 indicate increased cluster separation.

Finally, extracted density and effective atomic numbers from the calibration phantom and sample 4 were fed to the mathematical relationship given in equation (3) to extrapolate linear attenuation coefficients of tissues for arbitrary (virtual) monochromatic energy levels. The uncertainty was propagated from the standard deviation of obtained ρ and Z_eff values.

3.1. Calibration phantom results

The accuracy of material decomposition in phantom materials following the Gaussian fitting in the histogram space is shown in figure 3(a). Quantities x₁ and x₂ in equations (4a ) and (4b ) are replaced with x_PMMA and x_Al because of the particular basis choice. Figure 3(b) corresponds to the decoupling of density and effective atomic number in equations (5a ) and (5b ). The fitting parameters A, S_u , S_v defined in equation (7), and angle θ for phantom materials are given and further discussed in appendix B.

Zoom In Zoom Out Reset image size

The decoupled values from figure 3(b) were calibrated to absolute density and effective atomic numbers (table 1) applying the least-squares fitting method to the functional forms given in equations (6a ) and (6b ). Both calibration curves were based on the phantom material inserts in the range of interest for soft tissue imaging and the mapping to correct density and the effective atomic number was obtained at R² = 0.998 and R² = 0.997, respectively. Calibration curves are given in figure 4 and corresponding density and Z_eff values are given in table 2. Using equation (9), % errors in calibration materials were estimated to be below 3% and 1.5%, respectively, after the density and effective atomic number calibration.

Zoom In Zoom Out Reset image size

Table 2.  Z_eff and ρ with the corresponding standard errors σ_Z and σ_ρ obtained from the calibration phantom and %err computed using equation (9).

Material	ρ ± σρ (g cm−3)	Zeff $\pm {\sigma }_{{Z}_{\mathrm{eff}}}$	%err ρ	%err Zeff
PE	0.963 ± 0.004	5.36 ± 0.02	2.4	1.3
Water	0.971 ± 0.001	7.51 ± 0.02	2.9	0.9
PA	1.149 ± 0.004	6.10 ± 0.01	0.7	1.0
PMMA	1.198 ± 0.004	6.40 ± 0.01	0.7	1.4
POM	1.411 ± 0.004	7.04 ± 0.01	1.0	0.4
PTFE	2.204 ± 0.003	8.55 ± 0.01	0.2	0.2

3.2. Breast mastectomy results

Tissues inside the mastectomies delineated by a radiologist were quantitatively analyzed in terms of their density and effective atomic number and presented in figure 5, with average values given in table 3. For comparison, decomposition to basis material coefficients x_PMMA and x_Al are also given in figure 5. Skin tissue was not included in the graph as it is anatomically well separated from other tissues of the breast.

Zoom In Zoom Out Reset image size

Table 3. Average Z_eff and ρ from 5 mastectomy samples.

Tissue type	ρ ± SDρ (g cm−3)	Zeff $\pm {\mathrm{SD}}_{{Z}_{\mathrm{eff}}}$
Adipose	0.90 ± 0.02	5.94 ± 0.09
Fibro-glandular	0.96 ± 0.02	7.03 ± 0.12
Tumor	1.07 ± 0.03	7.40 ± 0.10
Skin	1.08 ± 0.02	7.31 ± 0.06

The level of separation between adipose, fibro-glandular, and tumor tissue in ρ and Z_eff space was found to be 0.31 (on a scale of −1 to 1) using MSS. Pairwise comparison of adipose and fibro-glandular (MSS = 0.59), adipose and tumor (MSS = 0.74), and fibro-glandular and tumor (MSS = 0.17) showed that adipose tissue can be distinguished from the other two, while fibro-glandular and tumor are closer together in ρ/Z_eff space. In comparison, MSS in x_PMMA and x_Al space was −0.02. From figure 5 it can be seen that adipose tissue can be distinguished from fibro-glandular and tumorous purely based on effective atomic number values. On the other hand, fibro-glandular and tumorous tissues are overlapping both in their effective atomic numbers and densities and can be distinguished only if both quantitative values are observed together.

In addition to adipose, fibro-glandular, and malignant tissue, the mastectomy labeled as sample 4 contains well-separated pure glandular tissue and skin. Thus, an extensive investigation was conducted on this tissue. Extracted density and effective atomic number were 0.98 ± 0.01 and 5.42 ± 0.06 for adipose, 1.07 ± 0.02 and 6.56 ± 0.20 for fibro-glandular, 1.18 ± 0.02 and 6.88 ± 0.09 for glandular, 1.19 ± 0.01 and 7.04 ± 0.07 for tumorous, and 1.19 ± 0.02 and 6.87 ± 0.04 for skin tissue, respectively. In this particular mastectomy sample, mean density and Z_eff are not distinguishable between the skin and glandular tissue. Glandular tissue has higher Z_eff and density than fibro-glandular tissue, but almost the same density and slightly lower Z_eff than tumorous tissue.

Virtual μ values of phantom materials and tissues in sample 4 are given in figure 6. The obtained ρ and Z_eff values correctly map to experimentally measured μ values and provide information about tissue separation at lower energies with respect to those used for tissue scanning.

Zoom In Zoom Out Reset image size

Conventional CT scans offer a relatively low specificity when distinguishing between soft tissues of slightly different compositions. On the other hand, spectral imaging advances the ability to distinguish such tissues by probing their attenuation properties at several energy levels. Spectral information is used for the decomposition of data in basis materials, such as PMMA-Al. The decomposed material maps can be used as an intermediate step to estimating the uniquely defined physical quantities of imaged tissues. Utilizing these quantities, virtual monochromatic μ values at arbitrary energy levels can be further extrapolated. In this work, a new approach to material decomposition and ρ/Z_eff estimation was presented to characterize breast tissues.

A least-square fitting approach to two-basis material decomposition was adopted for an over-determined system when μ values can be measured at several energy levels. This was the case in our retrospective study with monochromatic beams, but the same approach can be applied to multi-threshold photon counting detectors. The procedure for decoupling density from Z_eff depends on the quality of performed material decomposition, but is at the same time independent of the method itself, and can be applied to other methods published in the literature. Experimentally computed concentrations x_PMMA and x_Al of basis materials and related basis x_ρ and x_Z in the phantom were found to be in good agreement with the ground truth values. Moreover, small associated standard errors in figure 3 show that plastic inserts can be separated in both material and x_ρ /x_Z basis. The theoretically derived functions in equations (6a ) and (6b ) were fitted to 6 data points corresponding to tissue-equivalent plastics in the phantom. The high accuracy (R² > 0.99) of the fitting procedure justifies the assumptions made in theoretical derivation. The systematic errors observed in figure 3 for some materials (e.g. PTFE and PE), which can probably be ascribed to slight differences between the composition of the phantom materials and the ones published in Schoonjans et al (2011), do not impact the process of calibration, leading to an agreement with the ground truth data at an average % error of 1.34% for ρ and 0.89% for Z_eff, as shown in table 2. These results demonstrate that soft-tissue-equivalent plastic materials of similar composition can be distinguished from spectral CT data using estimated effective atomic number and density values.

The preliminary study on the breast cancer mastectomy samples was an attempt to demonstrate the feasibility of our method to distinguish between fibro-glandular and tumorous tissues inside the breast. Based on the available samples, it was shown that starting from spectral data it is possible to separate adipose, fibro-glandular and tumorous tissues based on their physical characteristics. This might be useful in risk assessment, cancer diagnosis, and the assessment of the status of the disease. Although pure glandular and tumorous tissue in Sample 4 had almost the same density and slightly different effective atomic numbers, this could be due to the desmoplastic core present inside the tumor and no conclusions could be made based on a single piece of evidence. The importance of Z_eff to x-ray attenuation can be observed when comparing adipose and fibro-glandular tissue clusters. It can be seen that lower μ values for adipose tissue are driven by lower Z_eff, rather than significantly lower density. Distinguishing tumorous and pure glandular tissue is challenging because only slight differences exist in both density and Z_eff. In our study, we observed that it is not possible to distinguish between different tissues solely on density or effective atomic number or μ value alone while reasonable discrimination (MSS = 0.31) can be obtained considering 2D clusters in ρ/Z_eff space. A worse separation (MSS = −0.02) can be obtained using just material coefficients x_PMMA and x_Al. Therefore, using density and effective atomic number maps in the diagnostic workflow could be beneficial, potentially allowing the identification of the tissue type based on quantitative measurements. Because of improved tissue separation, interpretability, as well as the ability to measure these quantities with other experimental techniques, density and Z_eff should be preferred over a simple material decomposition approach.

VMIs are usually computed directly from decomposed basis (e.g. μ₁ and μ₂) using x₁ and x₂ in equation (1). Equivalently, physically relevant ρ/Z_eff space can be used to compute other quantitative maps established in clinical practice. Linear attenuation values in figure 6 were calculated using equation (3) at energies within and outside the energy range used in the experiment. Virtual μ values calculated using the Z_eff defined by Champley et al (2019) were in agreement (within the measurement error) with experimentally measured μ coefficients in both the phantom and sample 4 and they enabled comparison with experimentally obtained μ values in other studies. For sample 4, which contained tissues examined by Fredenberg et al (2018), μ values were compared at energies of 20, 30, and 40 keV with an average % error of 3.3, 7.3, and 3.4%, for glandular, adipose, and tumor tissue, respectively.

Our approach to tissue analysis could be directly applied to state-of-the-art synchrotron radiation breast CT setups currently developed in Trieste at Elettra Sincrotrone SYRMEP beamline (Longo et al 2019), and at ANSTO Imaging and Medical beamline in Australia (Tavakoli Taba et al 2021). With these experimental settings, it should be mentioned that multiple energy acquisitions would be required, thus resulting in an increased dose to the breast. On the other hand, due to the high contrast-to-noise ratio of phase-contrast images, the dose per scan could be reduced thus bringing overall acceptable radiation exposures. Considering clinical systems in hospitals, the advent of spectral CT paves the way to material decomposition following a single shot acquisition without a significant increment of the dose. Photon-counting breast CT systems are in clinical practice (Berger et al 2019) and an extension to spectral applications is expected.

Studies estimating both ρ and Z_eff of human tissues using synchrotron CT systems are almost nonexistent. The study by Torikoshi et al (2003) introduced a method to compute these quantities avoiding the material decomposition task. An average accuracy of 0.9% for ρ_e and 1% for Z_eff is comparable to our method, but the dose level used was not reported. The analysis was performed on low-Z_eff plastic materials using a pair of monochromatic acquisitions. Considerably more papers using conventional systems have been published, but not focusing on breast tissues or breast dedicated scanners. Szczykutowicz et al (2011) performed a methodologically similar approach to the one presented in this paper using a clinical scanner and test object without any noise remedying approach. They successfully decompose electron density and Z_eff but at the cost of a significant reduction in signal-to-noise ratio. Lalonde and Bouchard (2016) developed a model in which materials are decomposed in a compressed basis with principal component analysis, using the fact that human tissues are composed of a very limited number of elements. Then, the first principal components are virtual materials containing a certain fraction of those elements. From there Z_eff values were computed. Azevedo et al (2016) implemented their System-Independent-Rho-Z (SIRZ) method to obtain physical quantities of phantom materials independent of the shape of the X-ray spectrum. The photoelectric-Compton decomposition is performed in sinogram space and absolute ρ_e and Z values are obtained after the calibration procedure. While physical characterization was successfully described for materials of fairly distinguishable compositions, the noise behavior was also not described in this work. Champley et al (2019) released a follow-up paper focusing on the optimization and simplification of spectral modeling as a new method called SIRZ-2. Most recently, Busi et al (2019) developed a physical characterization method using spectral detectors. They claim higher robustness and increased estimation accuracy (25%) compared to SIRZ methods. Extended work from the same group was published in Jumanazarov et al (2021) to optimize the computation speed. Machine learning solutions to ρ/Z_eff extraction were also tested by Su et al (2018) using dedicated phantoms with several tissue-equivalent inserts. Good results in computing Z_eff were obtained using artificial neural networks and the random forest method with a relative error between 1% and 2% at clinically relevant doses. However, a low-dose scanning and evaluation of the model on materials that were removed from the training set led to errors of up to 6%. Nonetheless, the authors showed that the machine learning approach is robust and computationally efficient.

Considering specifically the estimation of ρ and Z_eff for breast tissues, only a few studies exist and most have been made with experimental setups not used in diagnostic radiology. Berggren et al (2018) performed clinical evaluation of breast skin Z_eff obtained from 709 screening patients using planar spectral mammography. They reported slightly higher values of 7.365 (95% confidence interval: 7.364, 7.366) comparable to our findings of 7.31 ± 0.06, with a difference of less than 1% that might be due to the fact that tissues in our experiment were formalin-fixed. Gobo et al (2020) used a combination of transmission and scattering measurements with ²⁴¹Am source and an x-ray tube, while Antoniassi et al (2011) performed scattering measurements at 90 degrees by using low-energy beams. Given diverse formulations of Z_eff across the literature, we made a comparison in relative change using Z_eff of nylon as a reference value in table 4.

Our study shows potential for quantitative breast imaging by translating spectral information into the computation of physically relevant quantities, but it also has some limitations. The accuracy of the presented method is mainly governed by the quality of material decomposition, which is highly dependent on the denoising of decomposed data. We gave up the spatial information inside the ROI to obtain quantitatively correct material decomposition during the 2D Gaussian denoising approach. Inaccurate tissue segmentation would result in an erroneous Gaussian fitting procedure as the number of peaks in 2D histogram space corresponds to the number of tissue types being evaluated at once. The appearance of a histogram containing several plastic materials was published in our previous study (Vrbaski et al 2021). Thus, either correct tissue segmentation or an algorithm capable of correctly fitting multiple Gaussian functions in case no segmentation is performed are critical aspects of the present model. For the proof of concept, we relied on high-quality synchrotron beam radiation, but to make this approach broadly used, we plan to apply the method to more accessible polychromatic sources. Finally, the conclusions drawn in this study were based on the analysis of 5 mastectomy samples fixed in formalin which could slightly bias the measured μ values (less than 0.5% (Fredenberg et al 2018)). More samples will be evaluated to further confirm these findings. Despite the mentioned limitations, given the general validity of the proposed decomposition model and the foreseeable use of spectral detectors in breast CT scanners, the present feasibility study paves the way for its application to clinical spectral breast CT data.

A model incorporating CT reconstructions of an arbitrary number of spectral energy channels was developed to compute material density and effective atomic number. The density and effective atomic number of soft-tissue-equivalent plastic materials were computed with an average accuracy in the order of 1%, and the same approach was applied to the set of 5 mastectomy samples. The quantitative analysis presented here suggests that adipose, fibro-glandular, and tumorous tissues can be distinguished, given the MSSs obtained for each tissue pair. Density and effective atomic number can also be used for physics-based extrapolation of virtual monochromatic linear attenuation coefficients outside of the experimental energy range. Breast CT is an emerging technology that is capable to provide three-dimensional imaging of the breast and its extension to spectral imaging is certainly desirable for exploring the potential of material characterization in clinical trials.

This work has been supported by the Istituto Nazionale di Fisica Nucleare (INFN, National Scientific Committee 5 for Technological and Interdisciplinary Research) and Elettra-Sincrotrone Trieste S.C.p.A. Authors would like to acknowledge the whole SYRMA-3D collaboration members. We also would like to thank Dr Cristina Marrocchio, Radiologist and Research fellow at the University of Parma, Italy supporting us in tissue identification from CT and histology data sets. Sandro Donato has been supported by the 'AIM: Attraction and International Mobility'-PON R&I 2014-2020 Calabria and 'Progetto STAR 2'—(PIR01_00008)—Italian Ministry of University and Research.

All data that support the findings of this study are included within the article (and any supplementary information files). Data will be available from 1 June 2023.

As it was shown in section 2.5, regions of uniform composition result in scattered values of elongated shape (i.e. clusters) in the (x₁, x₂)-plane histogram. Such behavior was ascribed to the unavoidable amount of noise carried by the tomographic reconstructions on which the whole decomposition method is based. Detailed analysis of the propagation of measurement uncertainties, related to the size of the cluster, was carried out. The Gaussian fitting procedure introduced in section 2.5 individuates a preferred direction at angle θ aligned to the major axis u of elliptical Gaussian. The length in this direction and the orthogonal one (minor axis) are identified as the Gaussian spreads in a reference frame (x_u , x_v ), which is parallel to the major and minor axes of the spot. The standard errors of the centroid coordinates can be approximated by the ratios between the Gaussian spreads S_u and S_v and the square root of the volume under the Gaussian surface

$$\begin{eqnarray}&&{\sigma }_{{x}_{u}}=\displaystyle \frac{{S}_{u}}{\sqrt{V}},\end{eqnarray} \tag{ 1.1a }$$

$$\begin{eqnarray}&&{\sigma }_{{x}_{v}}=\displaystyle \frac{{S}_{v}}{\sqrt{V}},\end{eqnarray} \tag{ 1.1b }$$

where the volume is given as

$$\begin{eqnarray}&&V=2\pi \,A\,\displaystyle \frac{{S}_{u}\,{S}_{v}}{{b}_{1}\,{b}_{2}}\end{eqnarray} \tag{ 1.2 }$$

that it is equal to the total number of voxels contributing to the corresponding two-dimensional Gaussian function. Quantities b₁ and b₂ are bin sizes in both directions and A is the peak value. Calculated errors are then translated into uncertainties of coordinates x_ξ and x_ζ in the reference frame rotated by the angle ϕ defined in section 2.1. Rotation from the frame identified by angle θ to the frame identified by angle ϕ is defined as

$$\begin{eqnarray*}&&\left[\begin{array}{c}{x}_{\xi }\\ {x}_{\zeta }\end{array}\right]=\left[\begin{array}{cc}\cos \alpha & \sin \alpha \\ -\sin \alpha & \cos \alpha \end{array}\right]\left[\begin{array}{c}{x}_{u}\\ {x}_{v}\end{array}\right],\end{eqnarray*}$$

where α = ϕ − θ. Propagating uncertainty through the rotation of the reference frame leads to

$$\begin{eqnarray}&&{\sigma }_{{x}_{\xi }}^{2}={\cos }^{2}\alpha \,{\sigma }_{{x}_{u}}^{2}+{\sin }^{2}\alpha \,{\sigma }_{{x}_{v}}^{2},\end{eqnarray} \tag{ 1.3a }$$

$$\begin{eqnarray}&&{\sigma }_{{x}_{\zeta }}^{2}={\sin }^{2}\alpha \,{\sigma }_{{x}_{u}}^{2}+{\cos }^{2}\alpha \,{\sigma }_{{x}_{v}}^{2}.\end{eqnarray} \tag{ 1.3b }$$

From the rotated frame, to the rescaled x_ρ , x_Z frame, the division of the second coordinate by the first one implicates

$$\begin{eqnarray}&&{\sigma }_{{x}_{Z}}^{2}=\displaystyle \frac{{x}_{\zeta }^{2}}{{x}_{\xi }^{4}}{\sigma }_{{x}_{\xi }}^{2}+\displaystyle \frac{1}{{x}_{\xi }^{2}}{\sigma }_{{x}_{\zeta }}^{2}\end{eqnarray} \tag{ 1.4 }$$

while propagating x_ρ remains trivial

$$\begin{eqnarray}&&{\sigma }_{{x}_{\rho }}^{2}={\sigma }_{{x}_{1}}^{2}.\end{eqnarray} \tag{ 1.5 }$$

Finally, shifting to the pair ρ, Z requires equations (6a ) and (6b ), leading to the uncertainties

$$\begin{eqnarray}&&{\sigma }_{\rho }={\kappa }^{-1}\,{\sigma }_{{x}_{\rho }}\end{eqnarray} \tag{ 1.6a }$$

$$\begin{eqnarray}&&{\sigma }_{Z}=\displaystyle \frac{{\left({x}_{Z}-q\right)}^{\tfrac{1-\lambda }{\lambda }}}{\lambda \,{p}^{1/\lambda }}\,{\sigma }_{{x}_{Z}}.\end{eqnarray} \tag{ 1.6b }$$

In addition to the estimated centers of the distributions in figure 3(a), the 2D Gaussian fit method outputs several other parameters relevant for the accurate estimation of basis material concentrations summarized in table B1. The angle θ at which the cluster is extended remains constant (≃ −0.12) for all materials indicating that blurring is not dependent on the material type. It is rather a result of combined contributions of image acquisition, reconstruction, and material decomposition noise. The amount of blurring in the major direction S_u and the direction orthogonal to it S_v is of the same order of magnitude for all materials. The S_u values are around 2 orders of magnitude larger than S_v values. The amplitude A of the Gaussian function remains constant for the same size of the ROI for all plastic inserts. These properties are important for the Gaussian fitting procedure because an initial guess for fitting parameters can be given, improving the robustness of the method and increasing the computational speed. A more rigorous statistical description of these features will be the subject of a forthcoming standalone communication.

Effective atomic number as a property of a material depends on its chemical composition, but it is not uniquely defined in the literature. Very often different definitions are adopted by researchers depending on the experimental setup, energy range, and type of compounds. Several definitions that have been proposed to compute this quantity can be divided into: (i) methods that compute Z_eff as a weighted average of composing elements in the compound, (ii) and methods that rely on mass attenuation coefficients of composing elements to compute Z_eff (Bonnin et al 2014). The second group of methods was developed to solve the fundamental problem: Z_eff as the quantity defined in the first category is not specifically tied to the absorption property of a material. Thus, weighted sums of elements in a compound used in the first category are exclusively valid for certain energy ranges and often to a certain set of elements contained in the mixture, while methods in the second category rely on tabulated attenuation properties of materials to inherently define Z_eff as a quantity tied to attenuation property of the material.

In this section, we compared several available methods for the materials in the calibration phantom. ZcompARE (Vrbaski 2022) is user-friendly software with a graphical interface that can be used to compare several methods most often used. The comparison of the methods for the calibration phantom materials in the energy range 20–40 keV is given in table C1.

It can be seen that methods by Spiers et al, Tsai and Cho, Glasser et al, Gowda et al, and Champley et al provide very similar Z_eff numbers for our experimental conditions. Methods by Hine et al and Puumalainen et al gave considerably lower Z_eff numbers for all materials. Calibration functions of the form defined in equations (6a ) and (6b ) were applied to all definitions for Z_eff and results are given in figure C1. The MSS defines the level of separation in the range −1 to 1 was found to be approx. 0.31 for all methods used. Thus, the distinguishment of tissues could be obtained irrespective of the proposed definition.

Zoom In Zoom Out Reset image size

The method defined by Champley et al (2019) showed slightly better agreement of virtual monochromatic μ values computed using equation (3) with experimentally measured data. Constants n, K₁, and K₂ used in equation (3) were estimated to be 4.44, 9.4, and 1.6 using the least-square fit method for elemental materials found in the human body in the Z range of 1–20 across the diagnostic energy range (20–200 keV).