Particle shape characterizations for energetic materials by computational geometry and stereology method

Particle shapes significantly affect viscosity and flow behavior of energetic materials, and therefore affect their packability and processability. This study presents a computational geometry framework for automatically quantifying two-dimensional (2D) and three-dimensional (3D) particle shapes of energetic materials. A specimen by mixing three typical energetic materials including HMX (octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine), RDX (1,3,5-Trinitroperhydro-1,3,5-triazine) and AP (Ammonium Perchlorate) particles is used in this study. This specimen is scanned by high-resolution X-ray computed tomography (X-ray CT), yielding a volumetric image. An improved watershed analysis algorithm is used to process the volumetric image to identify individual 3D particles. The stereology sampling method is used to obtain 2D projections of 3D particles. Computational geometry techniques are developed by this study to analyze 2D particle projections and 3D particle geometries to compute seven commonly used shape descriptors, including convexity, circularity, intercept sphericity, area sphericity, diameter sphericity, circle ratio sphericity, and surface area sphericity. Results show that those different shape descriptors of energetic materials can be divided into three groups based on their numerical ranges. This study also evaluates the effectiveness and accuracy of 2D shape descriptors for quantifying the true 3D shapes. The inconsistent characterization results between 2D and 3D shape descriptors suggest that researchers should be cautious when using 2D images to characterize 3D particle shapes of energetic materials. The computational geometry framework and particle shape analysis results presented in this study can be potentially useful in numerical modeling, experimental analysis, and theoretical investigation for energetic materials. The 3D particle geometries of typical energetic materials are scanned by X-ray CT and identified by improved watershed analysis. The 2D and 3D computational geometry techniques are developed for determining seven commonly used particle shape descriptors. The stereology strategy is used to obtain 2D projections of 3D particles and evaluate the effectiveness of 2D shape descriptors to infer 3D particle shapes of energetic materials. The 3D particle geometries of typical energetic materials are scanned by X-ray CT and identified by improved watershed analysis. The 2D and 3D computational geometry techniques are developed for determining seven commonly used particle shape descriptors. The stereology strategy is used to obtain 2D projections of 3D particles and evaluate the effectiveness of 2D shape descriptors to infer 3D particle shapes of energetic materials.

1. The 3D particle geometries of typical energetic materials are scanned by X-ray CT and identified by improved watershed analysis. 2. The 2D and 3D computational geometry techniques are developed for determining seven commonly used particle shape descriptors.
3. The stereology strategy is used to obtain 2D projections of 3D particles and evaluate the effectiveness of 2D shape descriptors to infer 3D particle shapes of energetic materials.

Introduction
Energetic materials are mixtures of a large volume of hard explosive particles (90-95% by weight) and a small percentage of soft binder [1]. The key macro properties of energetic materials including sensitivity to shock, the possibility of initiation, and the subsequent chemical reaction are known to be controlled by processes occurring at their microstructure levels [2]. Under shock/impact loadings, the ignition of chemical reactions in explosive particles is controlled by localized regions of elevated temperature, i.e., hot spots [3]. A hot spot produces a small burning region that propagates in high explosives that is responsible for the buildup to detonation [4]. High explosives have complex, stochastic microstructures, which significantly affect formations of hot spots [5]. The explosive particle size and shape are the important parameters controlling the formation, ignition, and growth of hot spots in the shocked high explosives and then the macroscopic behavior of the high explosives in the shockto-detonation transition process. Therefore, the influence of explosive particle size and shape on the shock-to-detonation transition process have received increasing attention [6]. Many experimental studies have shown that shock sensitivity of heterogeneous high explosives increased when particle size was decreased [7]. Ershov et al. [8,9] investigated influences of particle sizes for HMX (octahydro-1,3,5,7tetranitro-1,3,5,7-tetrazocine), RDX (1,3,5-Trinitroperhydro-1,3,5-triazine) and PETN (pentaerythritol tetranitrate). They found that smaller particles significantly reduced reaction zone and the maximum mass velocity. In explosive formulations, shape of the particles also plays a key role in realizing better packability and processability. Kim [10] reported that spherical particles of explosives improved the mix-fluidity of formulations and had a great impact on the performance and insensitivity towards a sudden shock than non-spherical crystals. Sarangapani [11] found that spherical particles of 3-Nitro-1,2,4-triazol-5-one (NTO) reduced the viscosity by about 75% for the same explosive content to binder ratio in comparison with non-spherical NTO particles. Porosity data also supported the same conclusion that improvement in explosive loading by 10% by the usage of spherical NTO particles over non-spherical NTO particles for the given explosive to binder ratio. Particle shape and size characterizations of explosives affected their viscosity and flow behavior [12,13]. It has been reported that an increase in mix-fluidity of composite explosives resulted in improved solid explosive powder loading in liquid binder matrix [14]. Therefore, characterizations of particle shapes are important research topics to the explosive industry.
The particle shape and size characterizations significantly affect the engineering behavior of particulate materials, such as sands, powders, and pharmaceuticals. Numerous studies have shown that compared to rounded and spherical particles, angular and elongated particles exhibited larger values of strength and porosity, since angular and elongated particles did not rotate or slide as easily as rounded and spherical particles [15][16][17].
Various particle shape descriptors were proposed based on two-dimensional (2D) particle projections and threedimensional (3D) particle models [18]. Many studies focused on developing 2D and 3D image processing algorithms to automate particle shape characterizations of granular materials. These algorithms included Fourier analysis [19], angularity indices [20] and a fractal technique [21]. However, studies about particle shape characterizations of energetic materials were rare. Therefore, this study aims to present a computational geometry technique for computing 2D and 3D particle shape descriptors for typical energetic materials and evaluating the effectiveness of using 2D shape descriptors for evaluating 3D particle shapes following the stereology method.
The paper is organized in the following structure. First, seven commonly used 2D and 3D particle shape descriptors are introduced including convexity, circularity, intercept sphericity, area sphericity, diameter sphericity, circle ratio sphericity, and surface area sphericity. Next, three energetic materials including HMX (octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine), RDX (1,3,5-Trinitroperhydro-1,3,5triazine) and AP (Ammonium perchlorate) are used to prepare the specimen, which is scanned by X-ray computed tomography (X-ray CT) to obtain 3D particle geometries. The stereology sampling strategy is used to obtain 2D particle projects from 3D particle geometries. In the next section, based on the scanned particles, computational geometry techniques are developed to compute 2D and 3D shape descriptors. Finally, the 2D and 3D shape particle descriptors are compared to evaluate the effectiveness and accuracy of using 2D particle shape descriptors to infer 3D particle shapes of energetic materials. The computational geometry framework and particle shape analysis results presented in this study can be potentially useful in numerical modeling, experimental analysis, and theoretical investigation for energetic materials and other particulate materials.

Shape descriptors
There are seven commonly used particle shape descriptors for particulate materials, as shown in Table 1, which are initially proposed for characterizing 3D particle shapes. However, these 3D shape definitions were proposed around 1950's, when 3D imaging techniques were not available yet. Therefore, researchers must simplify 3D shape definitions into 2D shape definitions as shown in Table 1, because 2D particle projections can be easily obtained. The 2D and 3D shape descriptors developed for general particulate materials are used in this research to analyze particles of energetic materials.

Specimen preparation and digitization
Three energetic materials including HMX, RDX, and AP are used in this study because they the most commonly used explosive ingredients in China. The HMX, RDX, and AP used in this study are manufactured by Chinese Xi'an Modern Chemistry Research Institute. Both HMX and RDX are white crystalline powder and nitramines [22] as shown in Fig. 1a and b. The density dependence of the detonation velocity is given in the literature by the same formula [23]. RDX is slightly soluble in water (0.07 g per 100 g of water at temperature of 20°C) while HMX is insoluble (0.005 g per 100 g of water at temperature of 20°C) [24]. AP (Fig. 1c) is a colorless or white solid, soluble in water. It is a powerful oxidizer and often used as rocket propellants. The particle size distributions of three materials are determined by a laser diffraction particle size analyzer (Mastersizer 2000) as shown in Fig. 1d.
HMX, RDX, and AP particles have similar shapes as shown in Fig. 1. It is not particularly meaningful to investigate their shapes separately. Therefore, we mix them using an equal weight ratio to create a mixed sample. The mixed sample is used to fill a small plastic cylinder with a diameter of 5 mm and a height of 8 mm is to create a small specimen. This small specimen is scanned by a high-resolution X-ray CT scanner (ultra-high precision X-RAY_X-9200S manufactured by SXRAY) with a resolution of 0.98 μm, producing a grayscale volumetric image. The mean particle sizes of HMX, RDX, and AP is around 100 μm, which include about 100 voxels providing sufficient accuracy for scanning [25]. After scanning, the X-ray CT data is verified through mass fraction of the components to ensure the accuracy of 3D particle reconstruction [26].
In this volumetric image, the voxel values are proportional to densities of materials. Therefore, the explosive particles with relative high density generates brighter regions, while the air with low density generates darker regions. Then, a computer program is developed by this study to analyze the volumetric image and identify individual particles of energetic materials.
In the first step, the grayscale volumetric image is converted to a binary volumetric image by setting a threshold according to the Otsu's algorithm [27]. Then, the voxel values higher than the specified threshold is set to ones (white regions), which are corresponding to solid particles. The voxel values lower than the specified threshold is set to zeros (black regions), which are corresponding to air.
The solid particles (voxel values are ones) are used in further analysis for particle shape characterizations. However, those particles are contacting with each other. The watershed analysis is usually used to segment contacting particles [28]. The basic idea of watershed analysis is to identify constrained regions in the volumetric image, which are typically contacts between particles. Then the watershed analysis segments constrain regions to identify individual particles. However, some explosive particles also have constrained regions in their bodies, which are termed as "peanut-shaped particles" as shown in Fig. 2a. The watershed analysis cannot discriminate the constrained regions from peanut shaped particles and contacting particles. As a result, the watershed analysis mistakenly segments these peanut-shaped particles into two pieces or more pieces as shown in Fig. 2b. This paper introduces an improved watershed analysis technique proposed by Sun et al. [29] to address this issue.
The basic idea of improved watershed analysis is to analyze all the constrained regions in the volumetric image scanned by X-ray CT. Then, the constrained regions belonging to the peanut-shaped particles and belonging to the contacts between particles are discriminated by using a morphological analysis technique [29]. The improved watershed analysis only segments the constrained regions at contacts between particles to address oversegmentation issue such as Fig. 2c. The effectiveness and accuracy of the improved watershed analysis was validated by a large number of particles from six particulate specimens with different particle shape and size characterizations by Sun et al. [29]. The improved watershed analysis is used to analyze and segment the volumetric image of energetic material collected by this study to identify individual particles. A total of 940 particles are identified as shown in Fig. 3, of which three typical particles are shown in Fig. 4.
The stereology sampling method is used to obtain 2D projections of 3D particles [30][31][32]. The stereology is a branch of morphometry, which is used to infer 3D properties from 2D sections of 3D structures. A change in section direction with respect to the 3D structures generally results in a widely different set of projections, especially if the 3D structure is anisotropic and inhomogeneous in space [32]. The HMX, RDX, and AP particles all have highly irregular shapes, yielding completely different projections from different viewing angles. Following the stereology sampling strategy proposed by Cruz-Orive and Weibel [31] and Altunkaynak et al. [32], this paper designs a sampling strategy to obtain 2D sections. For each particle, a computer program is developed to automatically rotate the particle to obtain three representative projections: the maximum projection, the intermediate projection, and the minimum projection as shown in Fig. 4.

Computational geometry techniques for determining shape descriptors
The parameters including principal dimensions (d 1 , d 2 , and d 3 ), surface area (A s ), volume (V), maximum inscribing sphere, minimum circumscribed sphere, and 3D convex hull are required for computing 3D shape descriptors in Table 1. This paper develops computational geometry A is the area of the 2D particle projection A c is the area of the minimum convex boundary circumscribing the 2D particle projection Mora and Kwan [43] 3D V is the volume of the 3D particle V c is the volume of the minimum convex hull circumscribing the 3D particle A is the area of the 2D particle projection P is the perimeter of the 2D particle Altuhafi et al. [15] 3D V is the volume of the 3D particle A s is surface area of the 3D particle Intercept sphericity 2D S I = A is the area of the 2D particle projection A cir is the area of the smallest circumscribed circle of the 2D particle projection Riley [45] 3D S A = V V cir V is the volume of the 3D particle V cir is the smallest circumscribed sphere of the 3D particle Diameter sphericity 2D S D = D e D cir D e is the diameter of circle having the same area as the 2D particle projection D cir is the diameter of the minimum circumscribed circle of the 2D particle Wadell [46] 3D S D = D e D cir D e is the diameter of sphere having the same surface area as the 3D particle D cir is the diameter of the minimum circumscribed sphere of the 3D particle Circle ratio sphericity 2D S C = D ins D cir D ins is the diameter of the largest inscribing circle of the 2D particle projection D cir is the diameter of the minimum circumscribed circle of the 2D particle projection Santamarina and Cho [47] 3D S C = D ins D cir D ins is the diameter of the largest inscribing sphere of the 3D particle D cir is the diameter of the minimum circumscribed circle of the 3D particle Surface area sphericity 2D S P = P e P P e is the perimeter of the circle having the same area as the 2D particle projection P is the perimeter of the 2D particle projection Kuo and Freeman [48] 3D S P = A e A s A e is the surface area of the sphere having the same volume as the 3D particle A s is surface area of the 3D particle It is evident that only particle surfaces are useful for determining the parameters required by 3D shape descriptors in Table 1. However, after improved watershed analysis, each particle is a solid model. As such, the solid particle model is converted to "STL" particle model. The STL particle model only contains the surface geometries of particles, facilitating the shape analysis in this research. The STL is the abbreviation of Standard Tessellation Language, which is a computer graphic format widely used on a CAD (computed aided design) software, 3D printing and rapid prototyping [33].
The (Fast Fourier Transform) FFT-based 3D boundary reconstruction method proposed by Mineo et al. [34] is used to extract boundary points of solid particles and reconstruct the 3D particle surface to form the STL particle model in Fig. 5a. FFT-based 3D boundary reconstruction method has many advantages such as eliminating problem of defining a specific polynomial function order for optimum polynomial curve fitting and robustness for outof-plane noise, enabling reliable detection surface points and accurate surface reconstruction of irregular particles of energetic materials.
As shown in Fig. 5a, the STL particle models are essentially formed by small triangles and vertices. Therefore, knowing the spatial coordinates of three vertices, the surface area of the small triangle formed by three vertices can be computed, such as the triangle V 1 V 2 V 3 in Fig. 5b. Therefore, surface area of a 3D particle (S) can be determined by adding the areas of all the small triangles.
To determine volume of the 3D particle, a randomly point (O) inside particle surface can be chosen. Then, by connecting the vertices of the 3D particle and the point O, each small triangle forms a tetrahedron, such as the tetrahedron OV 1 V 2 V 3 in Fig. 5b. Knowing the coordinates of V 1 , V 2 , V 3 , and O, the volume of the small tetrahedron can be computed. Therefore, the volume of the 3D particle (V) can be determined by adding the volumes of all the small tetrahedrons. It is evident that the V is not depending on the location of point O as long as the point O is inside the particle surface. This paper suggests that O should be chosen as the centroid of the particle by computing the average of all the vertices or the center of mass of the particle to avoid possible skewed tetrahedrons.
A principal component analysis (PCA) [35,36] technique is used in this study to determine length (d 1 ), width (d 2 ), and thickness (d 3 ) of 3D particle. The PCA was firstly introduced by Pearson [37] as statistical method for analyzing 2D and 3D data points. The 3D particle vertices can be considered as 3D data points in space. Then, the PCA analyze the 3D particle vertices and find the first direction with  Fig. 6. Then, the length of the first principal component, such as the length of d 1 in Fig. 6, can be determined by projecting all the vertices on the d 1 direction. The d 1 is essentially the length of the 3D particle in Fig. 6.
In the next step, the PCA find the second direction with the second largest variance of the 3D particle vertices with constrained condition that the directions of first and second largest variances must be perpendicular to each other, such as the d 2 direction in Fig. 6. Then, the length of the second principal component, such as the length of d 2 in Fig. 6, can be determined by project all the vertices on the d 2 direction. The d 2 is essentially the width of the 3D particle in Fig. 6. Finally, the PCA find the direction have the third largest variance of the vertices with constrained condition that the directions of first, second, and third largest variances must be perpendicular to each other, such as d 3 direction in Fig. 6. Similarly, the length of d 3 can also be determined, which is essentially the thickness of the 3D particle in Fig. 6.
The 3D Euclidean distance map can be used to determine the maximum inscribing sphere of the particle as shown in Fig. 7a and b. The maximum Euclidean distance identifies the radius of the maximum inscribing sphere [38]. The coordinates of the maximum Euclidean distance values identify the center of the maximum inscribing sphere. The maximum inscribing sphere fit the original particle exactly as shown in Fig. 7c.
An iterative algorithm is developed to determine the minimum circumscribed sphere of a particle [39,40]. A 2D particle is used to illustrate the concept in Fig. 8. In the first step, only the particle surface vertices are used in Fig. 8a. Then, the minimum number of vertices bounding all the vertices is found, such as vertices 1 to 8, as shown in Fig. 8b. In the next step, two points with longest distance, such as vertex 1 and vertex 5 in Fig. 8a, are used as diameter to fit a trial circle. If this trial circle contains all the vertices, this will be the minimum circumscribed circle. If not, the points locate furthest outside of the trial circle is included in the fitting, such as the vertex 6 in Fig. 8c. Then, three vertices 1, 5, and 6 are used to fit the second trial circle. This procedure is repeated until no points locate outside the trial circle, yielding the minimum circumscribed circle as shown in Fig. 8d. The above computational process can be also applied to the 3D particle and the minimum circumscribed sphere is found as shown in Fig. 8e.
The convex hull of a particle is defined as the smallest convex shape that contains all the vertices of the particle. The 2D particle in Fig. 8a is also used to illustrate the concept. The smallest convex polygon containing all the vertices are shown in Fig. 8b, which is essentially the 2D convex hull of this 2D particle. The 3D convex hull of a 3D particle can also be identified using the same concept as shown in Fig. 9. As shown in Table 1, computations of the 2D shape descriptors require parameters including principal dimensions (d 1 , and d 2 ), area (A), perimeter (P), minimum circumscribed circle, maximum inscribed circle, and 2D convex hull. The 2D particle projections in Fig. 4 are analyzed to extract their perimeters, which is essentially a point cloud in 2D space. Then, the above computational geometry algorithms developed for 3D point cloud can also be used to analyze 2D point cloud to determine parameters required for computing 2D shape descriptors.

Relationship between 2 and 3D shape descriptors
The computational geometry techniques developed in this study are used to analyze the 3D particles and their projections to determine 3D shape descriptors and 2D shape descriptors, respectively. The results are shown in Figs. 10, 11, and 12. The 3D area sphericity (S A ), circle ratio sphericity (S c ), convexity (C x ), circularity (C), intercept sphericity (S I ), diameter sphericity (S D ), and surface area sphericity (S P ) are in ranges of 0.1 ~ 1.0, 0.2 ~ 1.0, 0.7 ~ 1.0, 0.5 ~ 1.0, 0.5 ~ 1.0, 0.7 ~ 1.0, 0.5 ~ 1.0 respectively. Therefore, these different shape definitions can be divided into three categories. The S A and S c definitions belong to the first category providing the widest numerical range for shape characterizations. The C, S I , and S P definitions belong to the second category providing the intermediate numerical range of shape characterizations. The remaining C x and S D definitions belong to the third category providing the smallest numerical range of shape characterizations.
It is evident that 2D and 3D shape descriptors give different shape evaluation results for the same specimen. The 2D S A , S c , C x , and C values underestimate the corresponding 3D image-based characterizations as shown in Fig. 10. The 2D S I values overestimate 3D S I values as shown in Fig. 11. The 2D S D and S P values vary within large ranges, causing uncertainties for analyzing 3D particle shapes, as shown in Fig. 12.

Discussion
The 2D particle shape descriptors have been widely used for describing 3D particle shapes of energetic materials and other particulate materials, although it is widely accepted that 2D particle shape descriptors may not be reliable and accurate to describe 3D particles. However, it is unknown that the accuracy levels (or error levels) of different 2D shape descriptors. This paper extensively investigates accuracy levels of seven commonly used 2D shape descriptors. The results showed the different 2D shape descriptors provide different estimation performance which could either overestimate or underestimate 3D particle shape descriptors for energetic materials. Therefore, it should cautious to use 2D images to characterize 3D particle shapes of energetic materials based on the comparisons. Today, many 3D imaging techniques are available to obtain digital 3D particle geometries, such as stereophotography, X-ray CT, neutron imaging, 3D laser scanning, white light scanning, and interferometry [41,42]. Therefore, it is necessary to use 3D geometries to accurately characterize particle shapes of energetic materials and other particulate materials. This study is limited in the following aspects and the future directions for research and recommendations should be considered. First, only three commonly used energetic materials having similar particle shapes are evaluated in this study. It is possible that other energetic materials with different particle shapes have different relationships between 2 and 3D particle shape descriptors. Therefore, the conclusions and results of this study should not be simply extended to other energetic materials. In future research, more energetic materials should be studied for extending the results of this study.
Second, the particle shapes affect the formation, ignition, and growth of hot spots in the shocked high explosives and then the macroscopic behavior of the high explosives in the shock-to-detonation transition process, but such mechanism is beyond the scope of this paper. Therefore, future research is required to investigate the fundamental mechanisms of the influence of particle shapes on macroscopic behavior of the energetic materials using experimental and numerical approaches. The computational geometry based particle shape characterization method developed by this study provide important foundation for these future studies. Third, this paper chooses three projections of each particle including the maximum projection, the intermediate projection, and the minimum projection to determine the 2D particle shape descriptors. However, the 2D particle shape descriptors computed from these three projections differ from the corresponding 3D particle shape descriptors. It is possible that such discrepancies can be minimized by choosing 2D projections from a specific orientation. Therefore, new algorithms and stereology sampling strategies should be investigated in future research to testify such hypothesis.

Conclusions
A specimen by mixing three typical energetic materials including HMX, RDX, and AP are scanned by high-resolution X-ray CT. A modified watershed analysis algorithm    is introduced by this study to process volumetric image to identify individual 3D particles of energetic materials. Then computational geometry techniques are developed to analyze 3D particles to determine seven commonly used particle shape descriptors, including area sphericity (S A ), circle ratio sphericity (S c ), convexity (C x ), circularity (C), intercept sphericity (S I ), diameter sphericity (S D ), and surface area sphericity (S P ).
Seven shape descriptors of energetic materials can be divided into three categories based on their numerical ranges. The S A and S c definitions belonging to the first category provide the widest numerical range for shape characterizations. The C, S I , and S P definitions belonging to the second category provide the intermediate numerical range of shape characterizations. The remaining C x and S D definitions belonging to the third category provide the smallest numerical range of shape characterizations.  The 2D and 3D shape descriptors are not consistent even for the same specimen of energetic material.The 2D S A , S c , C x , and C values underestimate the corresponding 3D S A , S c , C x , and C values. The 2D S I values overestimate 3D S I values. The 2D S D and S P values vary within large ranges, causing uncertainties for analyzing 3D particle shapes. The result suggest researchers should cautious to use 2D images to characterize 3D particle shapes of energetic materials.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.