Haralick texture analysis for microdosimetry: characterization of Monte Carlo generated 3D specific energy distributions

Objective. Explore the application of Haralick textural analysis to 3D distributions of specific energy (energy imparted per unit mass) scored in cell-scale targets considering varying mean specific energy (absorbed dose), target volume, and incident spectrum. Approach. Monte Carlo simulations are used to generate specific energy distributions in cell-scale water voxels ((1 μm)3–(15 μm)3) irradiated by photon sources (mean energies: 0.02–2 MeV) to varying mean specific energies (10–400 mGy). Five Haralick features (homogeneity, contrast, entropy, correlation, local homogeneity) are calculated using an implementation of Haralick analysis designed to reduce sensitivity to grey level quantization and are interpreted using fundamental radiation physics. Main results. Haralick measures quantify differences in 3D specific energy distributions observed with varying voxel volume, absorbed dose magnitude, and source spectrum. For example, specific energy distributions in small (1–3 μm) voxels with low magnitudes of absorbed dose (10 mGy) have relatively high measures of homogeneity and local homogeneity and relatively low measures of contrast and entropy (all relative to measures for larger voxels), reflecting the many voxels with zero specific energy in an otherwise sporadic distribution. With increasing target size, energy is shared across more target voxels, and trends in Haralick measures, such as decreasing homogeneity and increasing contrast and entropy, reflect characteristics of each 3D specific energy distribution. Specific energy distributions for sources of differing mean energy are characterized by Haralick measures, e.g. contrast generally decreases with increasing source energy, correlation and homogeneity are often (not always) higher for higher energy sources. Significance. Haralick texture analysis successfully quantifies spatial trends in 3D specific energy distributions characteristic of radiation source, target size, and absorbed dose magnitude, thus offering new avenues to quantify microdosimetric data beyond first order histogram features. Promising future directions include investigations of multiscale tissue models, targeted radiation therapy techniques, and biological response to radiation.


Introduction
Microdosimetry, the systematic study and quantification of the spatial and temporal distribution of absorbed energy in irradiated matter (Rossi and Zaider 1996), is used for numerous experimental and computational investigations designed to advance the field of Medical Physics.For example, recent studies have demonstrated that Monte Carlo (MC) simulations are capable of providing 3D specific energy (z, defined as the energy imparted ò per unit mass m) distributions using realistic microscopic modelling of cell populations (Villegas et al 2013, Oliver and Thomson 2018a, 2018b).In parallel, experimental techniques are being developed for highspatial resolution dosimetry, e.g.Raman spectroscopy in conjunction with radiochromic film (Mirza et al 2016, Mcnairn et al 2021), or microdetector arrays (Bachiller-Perea et al 2022).In general, considering microdosimetric quantities provides a complementary, microscopic perspective compared to conventional and ubiquitous dosimetric quantities such as absorbed dose ( = = ¯ D d dm z, where  d is the mean specific energy imparted to an infinitesimal volume of mass dm and z is the mean specific energy (ICRU 2011)).The significance of microdosimetric investigations is demonstrated by its potential for supporting development of novel treatment modalities (e.g.nanotherapy involving gold nanoparticles (Zygmanski andSajo 2016, Sakata et al 2018)), and new understandings of radiobiological effects of radiation exposure (Pater et al 2016, Villegas et al 2016, DeCunha et al 2021, Bertolet et al 2021).
Despite the importance of microdosimetry, analyses of microdosimetric distributions often rely on firstorder statistics (e.g.mean and standard deviation) with little-to-no consideration of the spatial distribution of data (Villegas et al 2013, Oliver andThomson 2018b).However, recent work introduced Haralick texture analysis (Haralick et al 1973)for the quantitative characterization of dosimetric distributions across micro-to macroscopic length scales (Mansour and Thomson 2023), focusing on examples to illustrate proof-of-concept.Initially proposed as a method to quantify the relation of neighbouring pixels within images, Haralick analysis has seen broad application for the classification of multidimensional data in fields such as image processing (Soh andTsatsoulis 1999, Clausi 2001) or computer vision (Dollar et al 2012) and with recent applications in medical physics, e.g.radiomics (Vallières et al 2015), dosiomics (Rossi et al 2018).
Haralick analysis uses a grey level co-occurrence matrix (GLCM) to quantify the spatial arrangement of intensity values within a distribution, often pixels within an image.To generate a GLCM, observed intensity values are typically mapped to a finite number of bins, often referred to as 'grey levels' with a maximum number of grey levels (N), in a process that is referred to as 'quantization'.Many projects have been looking for optimal quantization approaches in various applications (Soh and Tsatsoulis 1999, Clausi 2002, Gomez et al 2012) and no standardized method exists.Furthermore, it has been demonstrated that Haralick measures depend on the quantization method (Brynolfsson et al 2017) and researchers have investigated the use of different quantization approaches for use within Haralick texture analyses.For example, quantization methods and sensitivity to the number of grey levels have been investigated in different contexts e.g.image classification (Soh andTsatsoulis 1999, Clausi 2002) or radiomics (Vallières et al 2017).To address the sensitivity to quantization, Löfstedt et al (2019) developed a modified set of Haralick features with the intention of determining textural measures that are asymptotically invariant to N, while preserving interpretations of many texture descriptors.While Löfstedt et al (2019) demonstrated that the modified Haralick measures were asymptotically invariant to quantization for two datasets (comprising images from T1-weighted MRI volumes of the brain and histology images of colorectal cancer glands), this is not always the case (Garpebring et al 2018).
The present work explores the application of Haralick texture analysis to characterize 3D microdosimetric distributions of specific energy in cell-scale targets.We build on previous work (Mansour and Thomson 2023) which illustrated the potential for Haralick features to characterize microdosimetric distributions considering a single target size (cubic 3× 3 × 3 μm 3 volumes) and one magnitude of absorbed dose (40 mGy) for three source spectra ( 125 I, 192 Ir, and a 6 MV medical linac).In particular, we investigate how Haralick features can quantify changes in 3D microdosimetric distributions considering varying cell-scale target sizes, 6 photon spectra, and varying magnitudes of absorbed dose.Given that previous investigations have demonstrated considerable variation in specific energies in populations of cell-scale targets with varying target volume and dose (Villegas et al 2013, Oliver andThomson 2018b), we adopt the approach of Löfstedt et al (2019) to compute Haralick measures and assess sensitivity of computed features to the number of grey levels used in quantization.

Methods
The flowchart in figure 1 provides an overview of our methods to generate 3D specific energy distributions and then calculate modified Haralick texture measures.Similar to the work by Mansour and Thomson (2023), textural measures are calculated from quantized MC generated specific energy distributions.However, the work presented herein considers a wider microdosimetric parameter space with varying scoring voxel volumes (1-15 μm side lengths), mean specific energies (10-400 mGy) and incident photon source spectra (mean energy: 20 keV-2 MeV).The resultant expanded parameter space necessitates the use of modified Haralick analysis to reduce sensitivity to quantization.Section 2.1 describes the MC simulations used to generate the dosimetric data; the implementation of Haralick textural analysis is described in section 2.2.

Monte Carlo simulations of 3D microdosimetric specific energy distributions
All simulations are performed using the EGSnrc (Kawrakow et al 2021b) application egs_brachy (Chamberland et al 2016), pulled from the EGSnrc_CLRP githubegs_brachy branch with most recent commit 134c731.The XCOM photon cross section database (Berger et al 2010) and the NRC bremsstrahlung cross section database (Kawrakow et al 2021a) are used.Rayleigh scattering, bound Compton scattering, and electron impact ionization are turned on.The high-resolution random number generator option is used.All other transport parameters are EGSnrc defaults.The transport cutoff and production threshold for the kinetic energy of electrons and photons is 1 keV.Table 1 provides an overview of MC simulations performed (see figure 1 of the supplementary materials for a schematic diagram).A water phantom is modelled, consisting of a 40 × 40 × 40 array of cubic voxels as the scoring volume embedded in a larger scatter phantom, all with uniform mass density 0.998 g cm −3 .Cell-scale scoring voxels are considered, with side lengths of 1, 2, 3, 5, 7, 9, 11, and 15 μm.Previous work has demonstrated the reliability of EGSnrc for scoring in micron-scale geometries (Oliver andThomson 2018b, Martinov andThomson 2020).
Six photon source spectra with keV to MeV source energies are considered: 103 Pd, 125 I, 120 kVp, 192 Ir, 60 Co and a 6 MV medical linac photon spectrum, with minimum, mean, and maximum energies summarized in table 2. The 103 Pd, 125 I, and 192 Ir spectra are distributed with egs_brachy (Chamberland et al 2016) while the 60 Co and 6 MV spectra are distributed with EGSnrc (Kawrakow et al 2021b).The 120 kVp spectrum corresponds to an x-ray source with a tungsten target and 2 mm of aluminum filtration and is obtained from the Siemens x-ray spectrum simulation tool (https://bps.healthcare.siemens-healthineers.com/booneweb/index.html) (Boone   , 192 Ir and 120 kVp), 3.0 cm ( 60 Co) 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.
Interaction scoring is used to compute specific energy in each voxel.The mean specific energy over all scoring voxels, z, corresponds to the absorbed dose in the volume made up of all the scoring voxels.Absorbed doses between 10 and 400 mGy are achieved by varying the number of histories simulated.To estimate statistical uncertainty on absorbed dose, we replace the (40) 3 array of target voxels by a single scoring volume spanning their collective volume, computing the uncertainty using history-by-history statistics (Walters et al 2002).Statistical uncertainties vary with scoring voxel volume, source spectrum, and absorbed dose magnitude, and for k = 1 range from 10% ( 103 Pd, (1 μm) 3 scoring voxel volume, D = 10 mGy) to <0.3% (all spectra, (15 μm) 3 scoring voxel volume, D 150 mGy).
For each spectrum and absorbed dose magnitude considered, multiple realizations of the 3D specific energy distribution are determined by repeating MC simulations with different random number seeds; these realizations are used to assess corresponding variation in Haralick measures.The number of different realizations of specific energy distributions considered varies as a function of z and is motivated by the observed variation in statistical uncertainty of absorbed dose.At z=10 mGy, 50 realizations are considered.At 400 mGy, at least 4 realizations are considered.As an example of the calculation time: an 125 I simulation requires 4 h to achieve 3.5 × 10 9 histories (D = 40 mGy, (3 μm) 3 scoring voxel volume) on a single Intel Xeon 3.0 GHz CPU (model: 5160).

Haralick analysis
As part of the process to compute Haralick features, the specific energy distributions undergo a 'quantization' process (figure 1).This involves mapping (or binning) each value z from an MC-determined 3D specific energy distribution (denoted  in figure 1) to an integer value to yield a 'quantized' 3D array (denoted ¢  ).The quantization approach implemented herein uses the digitize() function  3.
As noted in section 1, we adopt the approach developed by Löfstedt et al (2019) to compute Haralick textural measures.Their approach was developed with the goal of obtaining modified Haralick measures that are asymptotically grey level invariant, which we investigate in section 3.4.The haralick() function available in the Python (Van Rossum and Drake 2009) package mahotas (Coelho 2013) is used in conjunction with the source code provided by Löfstedt et al (2019), available on tomlof github repository with most recent commit f83acbf.These calculations involve the construction of the GLCM, P θ .Each element in P θ (i, j) 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, θ.Thus, P θ (i, j) is an N × N matrix, where N represents the largest value in the quantized array ¢  (N £ N b ).The computation of a modified GLCM ( q p ), is outlined in Löfstedt et al (2019), and involves the explicit consideration of the number of quantization levels, N, for each respective GLCM.The mahotas package is used to calculate symmetric GLCMs for 13 of the 26 possible adjacency directions required to compute directionally-invariant features (Mansour andThomson 2023, Löfstedt et al 2019).The source code provided by Löfstedt et al (2019) is used to calculate modified GLCMs and textural measures.Due to a high level of correlation between textural measures (Vrbik et al 2019), we focus on five textural features: homogeneity ( H), contrast ( CON), correlation ( COR), local homogeneity ( LH), and entropy (  E), where the tilde notation follows the convention of Löfstedt et al (2019):

Results
In the following, we present Haralick features extracted from 3D specific energy distributions considering scoring voxels with side lengths between 1 and 15 μm (section 3.1), mean specific energies (absorbed dose) between 10 and 400 mGy (section 3.2), and 103 Pd, 125 I, 120 kVp, 192 Ir, 60 Co, and 6 MV incident photon spectra (section 3.3).The sensitivity of textural measures extracted from specific energy distributions to quantization (examining sensitivity to N B and corresponding Δz) is investigated in section 3.4.The supplementary materials for this work include additional specific energy distributions and figures of extracted Haralick measures to complement those presented herein.

Sensitivity to voxel volume
Figure 2 presents specific energy in the central plane of voxels generated using 1, 3, and 15 μm side length voxels, an 125 I source, and z equal to 10, 50 and 100 mGy.These figures illustrate variability in 3D specific energy distributions over the populations of voxels of different volumes.In particular, specific energy distributions generated for larger voxels (e.g. 15 μm side length) tend to be relatively uniform, with smaller differences in the magnitude of specific energy deposited across adjacent voxels, in contrast with smaller voxels (e.g. 1 μm side length) which have relatively large variations in specific energy deposited across adjacent voxels and many voxels receiving no specific energy.These trends are accentuated as the magnitude of z decreases.Figure 3 depicts Haralick textural measures extracted from specific energy distributions with 1, 2, 3, 5, 7, 9, 11, and 15 μm side length voxels considering 103 Pd, 125 I, 192 Ir, and 6 MV incident source spectra.These textural measures quantify the differences in the spatial distribution of specific energy induced by varying scoring voxel volumes.Different trends in extracted textural measures as a function of voxel volume are observed when considering distributions with different mean specific energies, and so Haralick features for 10, 50, and 100 mGy are presented as examples to illustrate trends.
Relatively large measures of homogeneity and local homogeneity in conjunction with relatively small measures of contrast and low measures of entropy (large and negative) are observed when extracted from distributions generated using 1 μm side length voxels for all magnitudes of z and incident source spectra considered.These trends in textural measures reflect the many voxels with zero specific energy.As voxel volumes increase (e.g. 3 and 5 μm side lengths), measures of homogeneity and local homogeneity decrease while measures of contrast and entropy increase, due to specific energy being deposited in more voxels but with relatively large differences in the magnitude of specific energy deposited.As larger voxel volumes are considered, trends in local homogeneity and contrast as a function of voxel volume are observed to reverse (contrast decreasing; local homogeneity increasing).These trends reflect a reduction in relative difference in the magnitude of specific energy deposited across adjacent volumes and a general increase of recurring patterns of specific energy deposited across adjacent voxels.
Correlation is observed to decrease as a function of voxel volume when considering relatively low energy sources (e.g. 103Pd, 125 I) reflecting a reduction in the linear relationship of specific energy deposited across adjacent voxels.For the higher energy sources ( 192 Ir, 6 MV), correlation varies little over the range of voxel sizes considered.

Sensitivity to mean specific energy
Figure 4 presents specific energy distributions for varying mean specific energy z (absorbed dose) and voxel volumes with a fixed incident source spectrum.Specific energy distributions for z equal to 10, 50 and 250 mGy with 1, 5, and 15 μm side length voxels and a 120 kVp source are presented.Specific energy distributions with large magnitudes of z in conjunction with relatively large voxel volumes (e.g.250 mGy, 15 μm side length) have relatively small differences in the magnitude of specific energy deposited across adjacent voxels.Specific energy distributions with lower magnitudes of z are observed to have relatively large variations in specific energy deposited across adjacent voxels; these effects are pronounced when considering smaller voxel volumes.For example, specific energy distributions with small magnitudes of z in conjunction with relatively small voxels (e.g. 10 mGy, 1 μm side length) have many voxels with zero specific energy with few non-zero values interspersed somewhat sporadically.Figure 5 presents textural measures extracted from specific energy distributions for 1, 5, and 15 μm side length voxels and 103 Pd, 125 I, 192 Ir, and 6 MV incident sources as a function of z.Measures of homogeneity and local homogeneity are largest while measures of contrast and entropy are smallest when extracted from specific energy distributions with relatively low z in conjunction with small voxels (e.g. 1 μm side length-figure 5 left hand panels), reflecting the many voxels with zero specific energy (figure 4(a)).As z increases and more voxels have non-zero specific energy (figures 4(b,c)), homogeneity and local homogeneity decrease while entropy and contrast increase.These trends differ from those for larger voxels (e.g. 15 μm side length-figure 5 right hand panels): homogeneity and local homogeneity increase with increasing z, while contrast and entropy decrease, reflecting decreasing variations in specific energy across voxels (figures 4(g,h,i)).The trends for more intermediate volume voxels (e.g. 5 μm side length) are similar to those for smaller voxel sizes for lower z and larger voxels for higher z (figure 5 centre column panels).
Correlation is observed to increase as a function of z when considering relatively low energy incident spectra (e.g. 103Pd and 125 I) in conjunction with large voxel volumes (e.g. 15 μm) reflecting an increase in linear predictability.For smaller voxel volumes or higher energy sources ( 192 Ir, 6 MV), correlation varies little over the range of z considered.

Sensitivity to incident spectrum
Figure 6 presents a subvolume of 3D specific energy distributions generated using 103 Pd, 120 kVp and 60 Co photon spectra, with voxel volumes of (3 μm) 3 , (9 μm) 3 , and (15 μm) 3 and z equal to 30 and 100 mGy.These figures illustrate variations in specific energy distributions for different sources.For lower magnitudes of absorbed dose and smaller voxels (e.g. 30 mGy, 3 μm side length voxels: figures 6(a,b,c)), the sources with the lowest energy have the largest number of voxels with zero specific energy and also large differences in the magnitude of specific energy in adjacent voxels (figure 6(a)).As mean source energy increases, specific energy is deposited across more voxels and fewer receive no specific energy (figures 6(b,c)).Overall, simulations with higher energy sources result in specific energy deposited more evenly across many voxels compared to lower energy sources.Figure 7 presents textural measures extracted from specific energy distributions for 3, 5, and 9 μm side length voxels and z of 30, 50, and 100 mGy as a function of the incident source mean energy.Extracted textural measures quantify differences in specific energy distributions due to different source spectra.For example, relatively high values of homogeneity and contrast in conjunction with relatively low measures of entropy are observed when evaluating specific energy distributions generated using a low energy incident spectrum (e.g. 103Pd) in addition to low magnitudes of z and small voxels (e.g. 30 mGy and 3 μm).These extracted textural measures reflect large variations in specific energy deposited across adjacent voxels, with many volumes receiving no specific energy as well as many near 150 mGy (upper end of the colour scale shown in figure 6).As the mean photon energy of the incident spectrum increases ( 125 I, 120 kVp, and 192 Ir) with these relatively low magnitudes of z and small voxels, we observe a reduction in measures of homogeneity and increasing entropy; this reflects specific energy being distributed across more voxels with varying magnitudes and fewer voxels receiving no specific energy (figure 6(b)).As even higher energy incident spectra are considered (figure 6(c)), or if specific energy distributions have relatively high magnitudes of z and larger voxels (figures 6(d,e,f,g,h,i)), decreasing measures of contrast and entropy in conjunction with increasing measures of homogeneity and local homogeneity are observed.These trends reflect increasingly uniform specific energy distributions, with relatively small differences in the magnitude of specific energy deposited, and recurring patterns of specific energy deposited across adjacent voxels.
For the voxel volumes and magnitudes of z investigated herein, measures of correlation are generally largest for 192 Ir source specific energy distributions.

Sensitivity to quantization
While previous sections considered quantization of specific energy distributions using N b = 2048 and Δz = 2.6 mGy for the Haralick analysis, the present section investigates sensitivity to these parameters (table 3).For each simulated specific energy distribution, distinct quantized distributions are evaluated with Δz between 1.3 and 11 mGy and Haralick texture measures are determined.Figure 8 presents Haralick textural measures extracted from 3D specific energy distributions generated using 103 Pd, 125 I, 192 Ir, and 6 MV sources, with 1, 5, and 15 μm side length voxels and mean specific energy of 50 mGy as a function of Δz.All extracted textural measures are observed to be sensitive to values of N b , Δz used for quantization, however sensitivity depends on the Haralick feature considered as well as the magnitude of mean specific energy and voxel size.For smaller voxel volumes (e.g. 1 μm side length), measures of homogeneity decrease as Δz increases (figure 8(a)) while entropy increases (figure 8(m)); however, contrast, correlation, and local homogeneity show less variation with Δz (figures 8(d,(g,j)).For larger voxels (centre and right-hand columns in figure 8), the variations in Haralick measures with changing Δz are less pronounced until the magnitude of Δz becomes large relative to variations in specific energies deposited across all voxels in the 3D scoring volume.Figure 9 presents the textural measure homogeneity extracted from 3D specific energy distributions generated using 103 Pd, 125 I, 192 Ir, and 6 MV sources, with 1, 5, and 15 μm side length voxels and Δz of 1.3, 2.6, 5.2, and 42 mGy as a function of of z.These figures illustrate that the overall trends in homogeneity versus z are consistent for the smaller values of Δz, whereas oscillatory-type behaviour with increasing z is observed with larger Δz.

Discussion
Recent microdosimetric studies considering specific energy distributions (Villegas et al 2013, Oliver and Thomson 2018a, 2018b) have focused on histograms and first-order statistics.However, histogram features ignore patterns in 3D specific energy distributions which are characteristic of radiation source, target size, and absorbed dose magnitude, making these blunt tools for analysis of exquisite microdosimetric data.Fortunately, Haralick textural analysis offers new avenues to quantify the 3D spatial distribution of specific energy.Our results show that the Haralick textural features of homogeneity, contrast, correlation, local homogeneity, and entropy successfully quantify spatial trends in the broad parameter space considered with keV to MeV photon source energies, cell-scale targets of varying size ((1μm) 3 -(15 μm) 3 ), and absorbed doses spanning an order of magnitude ( z ranging from 10 to 400 mGy).Trends in the Haralick measures for MC-determined 3D specific energy distributions are interpretable in terms of the underlying radiation physics, e.g.predominant interactions, electron ranges.
For all source spectra considered, the stochastic nature of radiation transport and energy deposition is apparent when evaluating specific energy distributions with relatively small voxels (e.g.(1 μm) 3 , (3 μm) 3 ) and low magnitudes of mean specific energy ( z, absorbed dose) for which many voxels have zero specific energy.Correspondingly, there are relatively higher measures of homogeneity and local homogeneity, and relatively low measures of contrast and entropy, compared with distributions for larger voxels and/or magnitudes of absorbed dose (e.g. 100 mGy and (15 μm) 3 ).With increasing target size, trends in Haralick measures-such as decreasing correlation ( 125 I, 103 Pd), homogeneity, and local homogeneity, increasing contrast and entropy (e.g.figure 3)reflect characteristics of each 3D specific energy distribution such as more voxels having non-zero specific energy.On the other hand, for larger voxels and increasing z (absorbed dose), we see homogeneity and local homogeneity increasing as specific energy is more evenly-shared amongst targets, as well as contrast and entropy decreasing; correlation varies with source type and voxel size (e.g.figure 5).Trends in Haralick measures characterize differences in specific energy distributions with source energy.For higher energy sources, Compton interactions predominate and set in motion energetic electrons with ranges large compared to the volume of targets (table 2), resulting in more sharing of energy amongst targets in the scoring volume.For lower energy sources, photons predominantly undergo photoelectric interactions, setting in motion electrons with R CSDA small in comparison with the overall 3D scoring volume, resulting in the energy of an incident photon being deposited in one or a few voxels, depending on voxel dimension.Specific energy distributions corresponding to the sources of different mean energies are characterized by the Haralick measures which vary accordingly.While detailed observations are provided in section 3.3 (e.g.figure 7), overall we see contrast decreasing with source energy; correlation and (local) homogeneity often (not always) higher for higher energy sources, and entropy often becoming larger and more negative with increasing source energy.We adopted the modified Haralick analysis approach developed by Löfstedt et al (2019) which they developed to reduce sensitivity to quantization method with the goal of 'asymptotic' grey level invariance (seen in some but not all applications (Garpebring et al 2018)).Mansour and Thomson (2023) considered a linear quantization scheme for their specific energy analysis using original Haralick measures, however, they considered a limited parameter space with only a single magnitude of absorbed dose (40 mGy) and one (3 μm) voxel size.The present work has a much broader parameter space with varying 1-15 μm target sizes, 6 sources with keV to MeV mean photon energies, and a correspondingly large range of magnitudes of absorbed dose (10-400 mGy).In preliminary investigations (results not presented), we considered the original Haralick measures in conjunction with various quantization methods, linear and non-linear schemes (Soh and Figure 7. Haralick features extracted from 3D specific energy distributions as a function of the mean energy of photons for 103 Pd, 125 I, 120 kVp, 192 Ir, 60 Co, and 6 MV sources for z equal to 30 mGy (a,d,g,j,m), 50 mGy (b,e,h,k,n), and 100 mGy (c,f,i,l,o) and voxels with side lengths of 3, 5, and 9 μm.The Haralick measures presented are mean values calculated from multiple realizations considering the same incident source, voxel volume, and z; the error bars displayed show the Type A standard uncertainty of these extracted features.
Tsatsoulis 1999, Di and Gao 2017), but we found that results were highly sensitive to quantization.We thus adopted the the modified Haralick texture analysis approach of Löfstedt et al (2019), using uniformly-spaced quantization levels with global minimum-maximum limits to enable comparisons of texture measures for different microdosimetric distributions.
While Löfstedt et al (2019) found that their modified Haralick features were asymptotically invariant to grey level for the two sets of images they considered, other work demonstrates that the features are not always invariant to quantization (Garpebring et al 2018).We observe some sensitivity to quantization (section 3.4).Haralick measures of homogeneity and entropy when extracted from specific energy distributions with small scoring volumes (e.g.(1 μm) 3 or (5 μm) 3 ) in conjunction with relatively low z and low energy sources (e.g. 103Pd or 125 I) are observed to be the most sensitive to quantization.Despite some sensitivity in extracted textural measures to quantization, trends in extracted measures as a function of z, voxel volume, and incident spectrum are observed to be subject to relatively small variations once Δz is sufficiently small (figure 9).Future work should continue to investigate quantization scheme in conjunction with Haralick analysis approach, with considerations specific to each application and context.
Interesting possibilities remain for applications of Haralick analysis to microdosimetric data.Building on the studies presented herein, future work will characterize spatial distributions of energy deposited within more sophisticated models and relate extracted textural measures to biologically relevant measures.Recent work investigating specific energy deposited within realistic cellular models which consider nuclei and cytoplasm compartments composed of representative media have demonstrated the sensitivity of resultant specific energy distributions to changes in source and geometric model of targets, and that specific energy distributions calculated with absorbed dose ∼mGy and ∼μm sized targets do not resemble normal distributions (Oliver and Thomson 2018a).Additionally, recent work has demonstrated that the spatial energy deposition clustering of ionizing radiation on nanometer length scales can be used as a surrogate for radiobiological parameters (Villegas et al 2016).Haralick texture analysis may be explored for spatial characterization of not only specific energy distributions, but also other microdosimetric quantities and biological response measures.

Conclusion
The present work demonstrates that Haralick texture analysis can quantitatively characterize 3D distributions of specific energy in populations of cell-scale targets.Measures of homogeneity, contrast, correlation, local homogeneity, and entropy quantify patterns of specific energy characteristic of the sources of radiation ( 103 Pd, 125 I, 120 kVp, 192 Ir, 60 Co, and 6 MV), voxel volumes (1, 2, 3, 5, 7, 9, 11, and 15 μm side length voxels), and mean specific energy (absorbed dose, between 10-400 mGy).Interesting possibilities remain for future work including application to cellular dosimetry and considering novel treatment modalities such as nanoparticles for radiation therapy.
Figure 9. Homogeneity extracted from 3D specific energy distributions with 1 μm (a,b,c,d,e), 5 μm (e,f,g,h), and (i), and 15 μm (i,j,k), l) side length voxels with Δz of 1.3 mGy (a,e,i), 2.6 mGy (b,f,j), 5.2 mGy (c,g,k) and 42 mGy (d,h,l) as a function of z .The Haralick measures presented are mean values calculated from multiple realizations considering the same voxel volume, incident source, and z; the error bars (often too small to see) represent the Type A standard uncertainty of these extracted features.

Figure 1 .
Figure1.Overview of the framework to generate microdosimetric distributions, process, and extract directionally invariant modified textural features.
(from Python (Van Rossum and Drake 2009) package numpy (Harris et al 2020)) with bin delimiters = i = 0, 1, 2, ..., N b .The values of z min and z max are set to the smallest and largest magnitudes of specific energy observed over all the microdosimetric datasets considered: integer N b defines the maximum number of quantization levels (or bins).Table 3 presents the different values of Δz, the quantization level 'width', alongside corresponding values of N b .Sections 3.1, 3.2, and 3.3 present Haralick features determined with N b = 2048 whereas section 3.4 assesses sensitivity to the value of N b , considering all values in table et al(2019)  note that the interpretation of the modified texture features does not generally change when compared to the original measures (mainly because computation of the measures involves intermediate rescalings/renormalizations).Entropy is an exception as the modified entropy measures can be negative(Löfstedt et al 2019).A concise interpretation of the original Haralick measures is provided elsewhere(Kumar et al 2014, Vrbik et al 2019, Mansour and Thomson 2023).

Figure 3 .
Figure 3. Haralick textural features extracted from 3D specific energy distributions as a function of voxel volume with the mean specific energy of 10 mGy (a,d,g,j,m), 50 mGy (b,e,h,k,n), and 100 mGy (c,f,i,l,o) for 103 Pd, 125 I, 192 Ir, and 6 MV sources.The Haralick measures presented are mean values calculated from multiple realizations considering the same voxel volume, incident source, and the error bars displayed show the Type A standard uncertainty of these extracted features.

Figure 5 .
Figure 5. Haralick textural features extracted from 3D specific energy distributions as a function of mean specific energy (z ) with 1 μm (a,d,g,j,m), 5 μm (b,e,h,k,n), and 15 μm (c,f,i,l,o) side length voxels for 103 Pd, 125 I, 192 Ir, and 6 MV sources.The Haralick measures presented are mean values calculated from multiple realizations considering the same z , incident source, and voxel volume; the error bars displayed show the Type A standard uncertainty of these extracted features.

Figure 8 .
Figure 8. Haralick textural features extracted from 3D specific energy distributions with nominal z equal to 50 mGy, 1 μm (a,d,g,j,m), 5 μm (b,e,h,k,n), and 15 μm (c), (f,i,l,o) side length voxels as a function of quantization level width (Δz).The Haralick measures presented are mean values calculated from multiple realizations considering the same voxel volume, incident source, and z; the error bars displayed show the Type A standard uncertainty of these extracted features.

Table 1 .
Overview of MC simulations: source, phantom, scoring voxels, and quantity scored.
 z m

Table 2 .
(Oliver and Thomson 2018b)with minimum, mean, and maximum photon energies; CSDA range is indicated for an electron of kinetic energy equal to the source mean photon energy (ICRU 2016).CSDA , for each source(Oliver and Thomson 2018b)to ensure build-up and scatter, and have side lengths of 1.0 cm ( 103 Pd, 125 I

Table 3 .
The maximum number of quantization levels (N b ), and corresponding quantization level width (Δz), considered for the investigation of textural measures to quantization (section 3.4).