Haralick texture feature analysis for characterization of specific energy and absorbed dose distributions across cellular to patient length scales

Objective. To investigate an approach for quantitative characterization of the spatial distribution of dosimetric data by introducing Haralick texture feature analysis in this context. Approach. Monte Carlo simulations are used to generate 3D arrays of dosimetric data for 2 scenarios: (1) cell-scale microdosimetry: specific energy (energy imparted per unit mass) in cell-scale targets irradiated by photon spectra (125I, 192Ir, 6 MV); (2) tumour-scale dosimetry: absorbed dose in voxels for idealized models of 125I permanent implant prostate brachytherapy, considering ‘TG186’ (realistic tissues including 0% to 5% intraprostatic calcifications; interseed attenuation) and ‘TG43’ (water model, no interseed attenuation) conditions. Five prominent Haralick features (homogeneity, contrast, correlation, local homogeneity, entropy) are computed and trends are interpreted using fundamental radiation physics. Main results. In the cell-scale scenario, the Haralick measures quantify differences in 3D specific energy distributions due to source spectra. For example, contrast and entropy are highest for 125I reflecting the large variations in specific energy in adjacent voxels (photoelectric interactions; relatively short range of electrons), while 6 MV has the highest homogeneity with smaller variations in specific energy between voxels (Compton scattering dominates; longer range of electrons). For the tumour-scale scenario, the Haralick measures quantify differences due to TG186/TG43 simulation conditions and the presence of calcifications. For example, as calcifications increase from 0% to 5%, contrast increases while correlation decreases, reflecting the large differences in absorbed dose in adjacent voxels (higher absorbed dose in voxels with calcification due to photoelectric interactions). Significance. Haralick texture analysis provides a quantitative method for the characterization of 3D dosimetric distributions across cellular to tumour length scales, with promising future applications including analyses of multiscale tissue models, patient-specific data, and comparison of treatment approaches.


Introduction
Across medical physics, energy deposited in tissues is critically important. This provides the impetus for advancements in absorbed dose determinations which, in turn, enable the evolution of radiation medicine. In radiation therapy contexts, state-of-the-art absorbed dose calculations in patient-specific models are becoming standard (Chetty et al 2007, Das et al 2008, Beaulieu et al 2012, Farr et al 2021. Monte Carlo (MC) simulations provide detailed absorbed dose distributions within realistic virtual patient models that enable new radiobiological modelling and outcomes analyses (Jin et al 2011, Miksys et al 2017, Vigneault et al 2019. Multiscale MC simulations are bridging macro-and microscopic regimes, considering both specific energy (z, defined as the energy imparted ò per unit mass m) in cellular targets as well as absorbed dose to tissue (¯= =  D d dm z, where d is the mean specific energy imparted to an infinitesimal volume of mass dm andz is the mean specific energy (ICRU 2011)), in diverse contexts, e.g., radiation therapy (Oliver and Thomson 2018a), experiments probing radiation response (Underwood and McMahon 2019), novel treatment approaches including nanoparticles in cancer radiation therapy (Zygmanski and Sajo 2016), and diagnostic imaging (Oliver and Thomson 2019).
While exquisite dosimetric distributions are being determined across length scales and are essential to advancement of the field, traditional analysis approaches rely on relatively rudimentary statistics (e.g., mean, standard deviation) and clinical metrics (e.g., D 90 , minimum absorbed dose received by hottest 90% of the target) to compactly summarize dosimetric data. These reductionist approaches result in loss of information regarding the spatial distribution of dosimetric data. This obscurement of spatial information in analyses may represent an important limitation.
We present herein an approach to address this shortfall in the spatial characterization of 3D (micro) dosimetric distributions by introducing Haralick feature texture analysis (Haralick et al 1973) in this context. Texture analysis quantitatively characterizes the spatial arrangement of intensity values in a region and has been widely applied in fields such as image processing (Soh and Tsatsoulis 1999) and computer vision (Dollar et al 2012). It has seen application in medical physics and bioinformatics, e.g., skin imaging (Ou et al 2014), radiomics (Vallières et al 2015), Raman spectroscopic mapping of radiation response (Vrbik et al 2019). Inspired by radiomics methods, which involve extraction of numerous spatial/texture features from medical images and the use of these features to predict therapeutic response, 'dosiomics' was recently proposed. Dosiomics focuses on extraction of spatial features from absorbed dose distributions for radiotherapy response prediction, with early applications focusing on toxicity rates in radiotherapy (Rossi et al 2018). These early studies have involved the calculation of many textural features and use of sophisticated mapping techniques for prediction of a subset of treatment outcomes (Gabryś et al 2018, Ebert et al 2021. Recent literature has demonstrated that considering a joint dosiomics and radiomics approach provides increased predictive capability over dose volume histogram derived metrics when attempting to correlate extracted features with treatment outcomes (Wu et al 2020). While dosiomics is gaining momentum for its predictive capabilities, the use of texture analysis for quantitative characterization and interpretation of 3D dosimetric distributions has not been considered to our knowledge.
The present work focuses on applying texture analysis to characterize 3D dosimetric distributions across cellular to tumour length scales, showing its interpretability and potential across traditional dosimetry and microdosimetry contexts. Haralick et al (1973) described readily computable textural features based on greylevel intensities of pixels in images. By considering either voxel absorbed dose or specific energy as analogous to pixel intensity, we demonstrate the application of Haralick analysis to dosimetric data generated from MC simulations and discuss trends in extracted textural features in terms of radiation physics and dosimetry.

Methods
The present work considers two scenarios to illustrate the approach and demonstrate its potential across microscopic-to-tumour length scales, as follows: 1. Cell-scale microdosimetry: 3D distributions of specific energy in μm-scale targets for photon sources; and 2. Tumour-scale dosimetry: 3D distributions of dose in mm-scale voxels for permanent implant prostate brachytherapy (PIPB).
The flowchart in figure 1 provides an overview of our methods.     scatter phantom, all composed of 0.998 g cm −3 water. The voxels have side lengths of 3 μm, comparable in size to cellular targets (Oliver and Thomson 2018a); previous work demonstrated the reliability of EGSnrc for scoring in micron-scale geometries . The transport cutoff and production threshold for the kinetic energy of electrons and photons is 1 keV, with interaction scoring used to compute specific energy in each voxel. Three source spectra are used: 125 I, 192 Ir, and a 6 MV medical linear accelerator photon spectrum. The minimum, mean, and maximum energies of each of the sources considered are presented in table 2. Phantom dimensions are motivated by considering the electron continuous slowing down approximation range for each source (Oliver and Thomson 2018a) to ensure build-up and scatter, and have side lengths of 1.0 cm ( 125 I and 192 Ir), and 6.0 cm (6 MV). The source is a parallel beam of square cross section with side length taken to be half that of the scatter phantom. For each source, sufficient histories are modelled to achieve a mean specific energy (dose),z, of approximately 40 mGy. This process is repeated with different random number seeds to obtain 15 different realizations of the specific energy distribution (z between 37.4 and 41.6 mGy) for analysis. As an example of the calculation time: simulation of an 125 I source requires 4 hours to achieve 3.5 × 10 9 histories (40 mGy) on a single Intel Xeon 3.0 GHz CPU (model: 5160). For all sources, the statistical uncertainty is less than 1.5% (k = 1) on this dose (40 mGy) which is estimated using history-by-history statistics considering a single scoring volume occupying the same space as the microdosimetric model (120 × 120 × 120 μm 3 ). The choice of 40 mGy is used to illustrate the application to cell-scale microdosimetry and is somewhat arbitrary: Haralick textural features may be explored over a wide range of mean specific energies (absorbed dose; results not presented herein) and further investigation is left to future work (section 4), including varying voxel size.

Tumour-scale dosimetry
The tumour-scale simulation geometry is shown in figure 3. The PIPB model is an idealized rectilinear prostate scoring volume with prostate volume, number of seeds, and seed spacing from Haidari et al 2019. The prostate model uses (2 mm) 3 scoring voxels and has dimensions 3.2 × 3.6 × 4 cm 3 (total volume 46.08 cm 3 ). A triangular arrangement is used to place 84 Amersham OncoSeed model 6711 125 I seeds with midpoint separation of 0.9 cm inside the prostate volume; the seed model was previously benchmarked (Safigholi et al 2020). For all PIPB calculations, absorbed dose is approximated as collision kerma which is calculated in each voxel using a tracklength estimator relying on tabulated values of mass-energy absorption coefficients distributed with egs_brachy (Chamberland et al 2016). Intraprostatic calcification modelling is motivated by Chibani and Williamson (2005) with (2 mm) 3 calcifications scattered randomly throughout the prostate volume. Intraprostatic calcifications from 0% to 5% by volume (%IC) are considered. For each level of %IC, 10 different realizations of calcification distribution are considered. To ensure sufficient scatter, 2 cm of male soft tissue surrounds the prostate volume. Elemental compositions and mass densities (table 3) of tissues follow the TG186 (Beaulieu et al 2012) recommendations and are from ICRU Report-46 (ICRU 1992) and Woodard and White (1986).
Two sets of calculation conditions are investigated. Calculations using these idealized models with tissuebased media assignments are referred to as 'TG186' in the following. Calculations based on the TG43 (Nath et al 1995) formalism (in which each voxel's elemental composition, mass density, and scatter material is assigned to that of 0.998 g cm −3 water and interseed attenuation is not modelled using egs_brachy's 'superposition' run mode) are referred to as 'TG43'. Absorbed dose distributions are normalized to achieve a nominal prescription dose D 90 = 140 Gy for simulations under TG43 conditions; the same scaling factor is used for TG186 simulations. Simulations of 10 9 histories are used for both sets of calculations and achieve average statistical uncertainties of less than 0.2% (k = 1) in the prostate volume; each calculation can be completed in under two hours on a single Intel Xeon 2.4 GHz CPU (model: E5-2683V4).

Haralick analysis
To compute the Haralick features, all dosimetric data generated must first undergo a 'quantization' process, consistent with Haralick analysis of image grey-levels (Haralick et al 1973). Effectively, this process maps each dosimetric datum (: specific energy or absorbed dose in a voxel) to an integer value ( ¢  : 1, K, N) representing quantization levels (or greyscale levels). The quantization approach implemented herein uses the Python digitize() function (Van Rossum and Drake 2009, Harris et al 2020) with the following bin delimiters, where  1 ,  2 are determined based on the scenario under investigation. For the cell-scale microdosimetry scenario,  1 is the smallest and  2 is the largest magnitude of specific energy observed over all microdosimetric datasets. For the tumour-scale dosimetric scenario,  1 is the smallest and  2 is the largest magnitude of absorbed dose (over all datasets) for a voxel assigned to prostate tissue which does not contain a seed; all voxels with absorbed dose higher than  2 are assigned to N (to limit the effect of very high values of absorbed dose that are greater than three times the prescription dose). The value of N is set to 100 for all analyses presented herein, following the approach and considerations of efficiency and sufficiency described by Soh and Tsatsoulis (1999).
The haralick() function available in the Python (Van Rossum and Drake 2009) package mahotas (Coelho 2013) is used to extract Haralick features. These calculations involve the construction of the grey level co-occurrence matrix (GLCM), P θ . Each element P θ (i, j) in this N × N matrix counts the number of times an adjacency pair, adjacent elements with values i and j, appears in the quantized distribution ¢  within a displacement of one voxel and for a particular adjacency direction, θ. There are 26 possible adjacency directions between neighbouring elements when considering the 3-dimensional rectilinear systems analyzed within this work. To increase efficiency, the routine used herein computes symmetric GLCMs; which are semi-directional . Homogeneity increases as adjacency pairs become more common and the distribution becomes more uniform. Local homogeneity is a measure of difference between adjacency pairs, and increases as these differences become smaller and the distribution generally becomes more uniform. Contrast increases as the difference between adjacency pairs increases. Correlation measures the level of linear predictability between adjacency pairs. Entropy quantifies the distribution of adjacency pair combinations, increasing as more occupied adjacency pair combinations are observed.

Results
In the following, we present Haralick textural measures extracted from specific energy distributions generated on cellular length scales considering multiple initial spectra with varying mean energies (sec. 3.1) and absorbed dose distributions generated for tumours with varying levels of prostatic calcification undergoing PIPB considering both TG43 and TG186 conditions (sec. 3.2). Figure 4 presents a subvolume of 3D specific energy distributions with mean specific energy of 40 mGy generated using 125 I, 192 Ir, and 6 MV photon spectra. These figures demonstrate the considerable differences in patterns of energy deposition within micron-scale voxels for different photon source spectra but for the same absorbed dose (mean specific energy). For example, relatively higher energy sources (e.g., 6 MV) produce more uniform distributions of specific energy, with smaller differences in specific energy deposited between adjacent voxels, compared to lower energy sources (e.g., 125 I).

Cell-scale microdosimetry
The differences in specific energy distributions generated using different initial spectra are quantified via the Haralick texture features extracted, depicted in figure 5. For example, homogeneity is largest when extracted from specific energy distributions generated using a 6 MV source, reflecting the relative smoothness and recurring patterns of specific energy deposited across adjacent voxels, observed in the specific energy distribution shown in figure 4. Contrast is highest when extracted from specific energy distributions generated using 125 I, reflecting the large variations in specific energy in adjacent voxels with many voxels receiving no specific energy as well as many near 300 mGy (upper end of the colour scale shown in figure 4); contrast decreases when extracted from specific energy distributions generated with increasing source mean energy. Correlation is largest when extracted from specific energy distributions generated using the 192 Ir source, reflecting a higher linear relationship between adjacent voxels. Local homogeneity extracted from specific energy distributions generated using the 6 MV source is higher than for either 192 Ir or 125 I, indicating that in addition to the overall smoothness observed in figure 5, specific energy distributions generated using a 6 MV source have relatively small differences in specific energy deposited across adjacent voxels. Entropy is largest when extracted from specific energy distributions generated using an 125 I source. The higher measure of entropy indicates large variations in magnitudes of specific energy deposited within voxels, producing many distinct adjacency pair values. Figure 6 presents 2D arrays of the 3D absorbed dose distributions generated using TG43 and TG186 conditions considering 0, 3, and 5 %IC. Comparison of the panels reveals differences in absorbed dose deposited due to calculation conditions (TG43 compared with TG186) as well as the effects of intraprostatic calcifications, reflected in the Haralick measures shown in figure 7.

Tumour-scale dosimetry
Dosimetric distributions generated under TG43 and TG186 conditions with 0%IC differ only slightly, with slightly higher values of absorbed dose observed in the centre of the slice for TG43 (figure 6(a)) compared with TG186 ( figure 6(b)); correspondingly, the Haralick measures are comparable, with relatively small differences observed: homogeneity and local homogeneity are higher for TG186; entropy is lower. Focusing on TG186 distributions with varying calcifications (%IC), there are relatively higher magnitudes of absorbed dose deposited in regions composed of calcification and lower doses in surrounding regions; correspondingly, the Haralick measures change markedly with increasing %IC. As %IC increases from 0 to 5%, contrast increases while correlation decreases. Homogeneity and local homogeneity initially decrease as %IC rises above zero, but then rise after 3%IC (homogeneity) or 4%IC (local homogeneity). Entropy increases as % IC changes from 0 to 3%, but then decreases once %IC reaches above 4%. Furthermore, relatively larger levels of variation in extracted textural measures are observed when considering patients with higher %IC.

Discussion
The Haralick textural measures extracted for both cell-scale and tumour-scale scenarios quantitatively characterize the corresponding 3D microdosimetric and dosimetric distributions. These textural features reflect differences in dosimetric distributions due to source spectra, calculation conditions, and the presence of calcifications within the patient model, all interpretable in terms of the underlying radiation physics, further described below. To our knowledge, the present work is the first time that texture analysis has been used for quantitative characterization and interpretation of 3D absorbed dose distributions across micro-to macrolength scales. Other work involving the application of texture analysis to dosimetric data (dosiomics) has employed an approach which emphasizes a general collection of textural features and mapping to treatment outcomes (Gabryś et al 2018, Rossi et al 2018, Wu et al 2020, Ebert et al 2021, Placidi et al 2021, without investigation of the interpretation of the textural features and connecting with underlying fundamental radiation physics and dosimetry. The Haralick texture features extracted from 3D specific energy distributions for the cell-scale scenario are sensitive to microdosimetric differences due to source spectra. For these microdosimetric distributions and associated Haralick measures, the predominant interactions and ranges of electrons set in motion therein may be used to interpret results. For the 125 I photons, photoelectric and Compton interactions dominate, and the most energetic electron set in motion by a photon with energy equal to the mean energy of the spectrum (mean Given the tortuous path of the electron in condensed media such as water in conjunction with this relatively short CSDA range, the energy of a primary photon will generally be deposited in only one or a few voxels. This results in large variations in specific energy across the 3D distribution, corresponding to high measures of contrast and entropy. On the other hand, Compton interactions dominate for the 6 MV photons (mean energy 2.0 MeV) and generate relatively high energy electrons with CSDA ranges large in comparison with the dimensions of the entire 3D scoring volume for the cell-scale microdosimetry scenario. Thus, energy is deposited across many voxels by slowing down electrons (Villegas et al 2013, Oliver andThomson 2018b). With relatively similar magnitudes of specific energies in adjacent voxels and a higher frequency of repeating adjacency pairs of (quantized) specific energy, higher measures of homogeneity and local homogeneity are observed for 6 MV compared with either 125 I or 192 Ir. Entropy is lowest, reflecting the fact that few distinct energy states are occupied. For 192 Ir photons (mean energy 350 keV), Compton interactions dominate but electrons have shorter range than for 6 MV and so deposit their energy amongst fewer scoring voxels. In comparison, distributions generated using 125 I and 6 MV sources are observed to have either highly sporadic or largely similar distributions of specific energy respectively. These differences in the distribution of specific energy are what contribute to relatively lower measurements of correlation for distributions generated using 125 I and 6 MV sources.
For the tumour-scale scenario, the Haralick measures quantify observed differences in 3D absorbed dose distributions due to simulation conditions (TG43 versus TG186) and %IC. Comparing Haralick measures for TG43 and TG186 with 0%IC, only relatively small intrafeature textural differences are observed and these correspond to the relatively small dosimetric differences observed (Miksys et al 2017). For TG186 (0%IC), interseed attenuation reduces the magnitude of absorbed dose deposited in the immediate vicinity of seeds, compared with TG43 simulations in which interseed attenuation is not modelled, resulting in higher measures of homogeneity and local homogeneity, and lower entropy for TG186.
The introduction of intraprostatic calcifications markedly changes absorbed dose distributions calculated under TG186 conditions (Miksys et al 2017), with corresponding changes in Haralick measures. Due to the relatively high effective atomic number of calcifications compared with prostate tissue (ICRU 1992), calcifications have enhanced photoelectric interactions that result in relatively high magnitudes of absorbed dose deposited within IC voxels but lower magnitudes of absorbed dose deposited within surrounding prostate tissue (due to attenuation; differences in fluence). The presence of these higher absorbed dose IC voxels markedly changes all Haralick measures. Extracted measures of contrast and correlation are observed to linearly increase and decrease respectively for any increase in %IC introduced. Increasing contrast and decreasing correlation indicates that the presence of calcifications causes large differences in absorbed dose deposited between adjacent voxels and a lower linear predictability between adjacent voxels for higher %IC respectively.
When considering extracted measures of homogeneity and entropy, the largest differences are observed when comparing textural features extracted from models which have 0% and 3% IC. Over the 0% to 3% IC range, the decreasing homogeneity reflects how the addition of calcifications causes relatively large changes of absorbed dose deposited across adjacent voxels (Miksys et al 2017), thus disrupting patterns of adjacency pairs that previously existed. Furthermore, increasing entropy over the 0% to 3% IC range indicates that the addition of calcification introduces a greater magnitude of distinct adjacency pairs observed.
As %IC increases from 0% to 4%, extracted measures of local homogeneity from dose distributions are observed to decrease; indicating that patterns of adjacency pairs are being disrupted, and the differences in dose across adjacent voxels are increasing. When considering higher levels of %IC, trends in extracted measures of homogeneity, entropy, and local homogeneity are observed to reverse. The reversal of all trends within extracted features reflects how the dosimetric changes caused by the presence of calcifications, such as adjacency pairs with relatively large differences in dose deposited, are becoming prevalent for patients with higher levels of %IC. Furthermore, when considering patients with higher levels of %IC, the disruption in distributions of absorbed dose caused by the random scattering of calcified voxels is observed to result in an increase in the relative variation of extracted measures of homogeneity, contrast, local homogeneity, and entropy.
The present work, focusing on Haralick texture analysis for quantitative characterization of dosimetric distributions and interpretation of extracted textural features in terms of radiation physics, is distinct from recent applications of texture analysis to dosimetric data in the emerging field of dosiomics. Dosiomics involves correlating measures extracted from dosimetric distributions (traditional Dose Volume Histogram (DVH) measures, textural, and other high order measures) to radiotherapy outcomes (Gabryś et al 2018, Ebert et al 2021, Placidi et al 2021, generally without presentation and interpretation of the textural measures collected from dosimetric distributions, as done in the present study. For example, Rossi et al (2018) compared prediction modelling for late toxicity rates in prostate cancer radiotherapy using a conventional DVH approach and a joint approach using DVH and dosimetric texture features (considering 42 textural measures). Wu et al (2020) investigated whether dosiomic features (considering 297 first order and textural dosiomics features) can be used to predict locoregional recurrences in patients undergoing intensity-modulated radiotherapy treatments. Adachi et al (2021) considered 6,808 dosiomic features (both first order and texture features) extracted from calculated absorbed dose distributions to predict radiation pneumonitis after lung stereotactic body radiation therapy. The present work illustrates that Haralick texture analysis can characterize dosimetric distributions across cellular-to-patient length scales, and that texture measures may be interpreted in terms of underlying physics.
Future work will build on the simplified scenarios considered herein to apply Haralick texture analysis in meaningful contexts across cellular and patient length scales. Microdosimetric studies will explore Haralick measures extracted from microdosimetric distributions that consider various target sizes, absorbed dose levels, and source spectra, as well as using dosimetric data from MC simulations involving realistic microscopic modelling of cell populations Thomson 2018a, 2018b) and novel treatment strategies, e.g., involving gold nanoparticles Thomson 2017, 2020). In tandem with these recent computational developments enabling microdosimetric modelling of cell populations in tissues, experimental techniques are being developed for high-spatial resolution dosimetry, e.g., Raman spectroscopy in conjunction with radiochromic film ( In the context of tumour-scale dosimetry, future work will apply Haralick texture analysis to characterize 3D realistic patient-specific absorbed dose distributions. Building on the simplified prostate model used herein, future work will use patient-specific dosimetric data coming from MC simulations that include realistic intraprostatic calcifications from CT images and contours. For example, Miksys et al used patient CT images and contour data to perform MC simulations, investigating different approaches to model intraprostatic calcifications; they observed varying calcification volumes within the prostate, with calcifications in some patients visible as more randomly-dispersed/distributed smaller volumes (similar to those considered herein) and in others as larger contiguous masses (Miksys et al 2017). An interesting future application could involve comparing different treatment modalities e.g., considering brachytherapy combined with external beam therapy (Hoskin et al 2013, Spratt et al 2017 and investigating how the collected textural measures may be able to inform treatment planning. While the examples in this paragraph are for prostate treatments, they are generalizable to other treatments and, indeed, to radiation medicine in general. For example, Haralick texture analysis could be used to quantitatively characterize and compare 3D absorbed dose distributions for different treatment modalities such as different external beam techniques, proton therapy, brachytherapy, and targeted radionuclide therapy. Extending those studies to then consider biological response might yield further insights, and could lead to practical considerations for clinical applications. Future applications of quantitative Haralick analysis in such contexts may necessitate investigation of robust quantization techniques and parameters which are being developed in other contexts (Brynolfsson et al 2017, Löfstedt et al 2019.

Conclusion
The present work demonstrates that Haralick texture analysis can quantitatively characterize 3D dosimetric distributions across cellular to tumour length scales. Haralick measures of homogeneity, contrast, correlation, local homogeneity, and entropy quantify patterns of specific energy characteristic of the sources of radiation ( 125 I, 192 Ir, 6MV) used in the microdosimetric scenario. On tumour scales, these same features characterize differences in 3D absorbed dose distributions due to simulation conditions or absorbed dose calculation methodology (TG43 versus TG186) as well as the effects of intraprostatic calcification. The trends in Haralick measures considered, namely homogeneity, contrast, correlation, local homogeneity, and entropy, are interpretable in terms of the underlying radiation physics. Interesting possibilities remain for future work, including assessment of 3D absorbed dose distributions for different radiotherapy treatments; applications to cellular dosimetry, including novel treatment strategies such as nanoparticles for radiation therapy; assessing patterns of double strand breaks or other biologically-motivated quantities and ensuing biological response.