Image Anal Stereol 2013;32:107-115 Original Research Paper 3D RECONSTRUCTION AND ANALYSIS OF THE FRAGMENTED GRAINS IN A COMPOSITE MATERIAL Luc Gillibertb and Dominique Jeulin Centre de Morphologie Mathématique, Mathématiques et Sytemes, Mines ParisTech, 35, rue Saint Honoré, 77305 Fontainebleau, France e-mail: luc.gillibert@cmm.ensmp.fr, dominique.jeulin@cmm.ensmp.fr (Received June 6, 2012; revised May 20, 2013; accepted June 9, 2013) ABSTRACT X-ray microtomography from solid propellant allows studying the microstructure of fragmented grains in damaged samples. A new reconstruction algorithm of fragmented grains for 3D images is introduced. Based on a watershed transform of a morphological closing of the input image, the algorithm can be used with different sets of markers. Two of them are compared. After the grain reconstruction, a multiscale segmentation algorithm is used to extract each fragment of the damaged grains. This allows an original quantitative study of the fragmentation of each grain in 3D. Experimental results on X-ray microtomographic images of a solid propellant fragmented under compression are presented and validated. Keywords: clustering, damage, mathematical morphology, reconstruction, segmentation, solid propellant. INTRODUCTION Propellants are composite materials made of linear elastic brittle grains embedded in a visco-elastic elastomer matrix. Damage can occur in the brittle grains of solid propellants under the action of a mechanical shock, inducing a possible unsafe behavior of the propellant. In order to characterize the progression of damage in such composite materials, it is interesting to analyze their microstructure at different steps of the fragmentation, as initiated in Gillibert and Jeulin (2011a). in the present paper, specimens of this material were fragmented under compression generated by the impact of a mass, and examined by means of microtomographic images obtained on a highresolution micro-CT system. From these 3D images, the goal is to estimate some statistics on each grain, that should be relevant to the progression of the 3D damage: the specific surface area of its cracks, the volume fraction of its cracks, the number of fragments and the size distribution of its fragments. For this purpose, original grains have to be reconstructed from the image of the fragmented material. Then, each fragment must be extracted, and must be associated to its initial grain. In this paper, we first introduce the type of materials and of 3D images that are studied. Then the original algorithm for the reconstruction of the particles from the observed fragmented image, based on two types of segmentation (one based on the h-minima, and one using the K-means clustering algorithm) is presented. Then a multiscale segmentation based on the stochastic watershed gives a 3D images of individual fragments. Finally, 3D image analysis measurements provide a statistical analysis of the local damage in the material, which gives a new approach to the local 3D study of damage in materials. The steps of our approach are illustrated by Fig. 1. MATERIALS AND 3D IMAGES The studied images are obtained by X-ray microtomography with a Skyscan 1172 highresolution micro-CT system at the CEA Gramat, a public laboratory affiliated with the Atomic Energy and Alternative Energies Commission. The material samples are a solid propellant in three states: the initial material, and fragmented materials with two steps of degradation generated by mechanical impacts. For all the studied images, a voxel is 3.6 jam. The original diameter of the grains is 400 jam, but there are many small fragments in the damaged specimens. The following images ares studied: - MAT1 is a 1014 x 1155 x 250 voxels image (3650.4 x 4158 x 900 jam3). A slice of this image is illustrated in Fig. 2. A mechanical impact is obtained from a 2 kg mass falling 15 cm. - MAT2 is a 1035 x 1008 x 428 voxels image (3726 x 3628.8 x 1540.8 jam3). It is a reference material without any mechanical impact (Fig. 3). However it contains some rare cracks. - MAT3 is a 1116 x 1104 x 424 voxels image (4017.6 x 3974.4 x 1526.4 jam3). A mechanical impact is obtained from a 2 kg mass falling 30 cm, with a 9.5 cm rebound (Fig. 4). Measuring the evolution of local damage in such materials is a challenge, requiring to extract the crack network in grains for their individual study. As seen in Fig. 2 and Fig. 4, the crack network is complex and many grains are highly fragmented, so that it is not easy to recover the initial grains from the image. The purpose of the algorithms developed in this paper is to give a reliable and automatic method to provide local estimates of the damage. 3D Images of a composite material • Initial state • 2 states of fragmentation Reconstruction of grains from fragmented images • Cracks closing • Watershed from h minima • Watershed from K-means Extraction of fragments by a multiscale stochastic watershed Statistical analysis of the grains fragmentation Damage characterization Fig. 1. Flowchart of the study of fragmented media. RECONSTRUCTION ALGORITHM OF THE FRAGMENTED GRAINS A classical approach is used for removing the cracks: morphological closing and a volumic opening removing connected components of the holes with a low volume (Matheron, 1967; Serra, 1982). Then, a watershed transform on the closed image is used. Introduced in 1979 by S. Beucher, the watershed is computed from a gradient image, here the inverse of the distance map of the closed image, and from a set of markers (Beucher and Lantuejoul, 1979). Two possible sets of markers for this watershed are explored. The first approach is topological and uses the h-minima filter (Soille, 2003). The minima of the inverse distance function with a depth is larger than h are used as a sets of markers. The use of the h-minima with the watershed on the distance map is very classical, but if the grains are very fragmented and if the fragments are scattered, the algorithm fails to reconstruct correctly some grains. The second approach is based on a method of cluster analysis, the K-means clustering, which aims to partition a set of observations into K clusters (Lloyd, 1982). Here, the observations are random voxels inside the mask of the grains. The kernels of the clusters, more precisely the center of gravity of the clusters, are used as sets of markers for the watershed. The number of classes for the K-means clustering algorithm is automatically computed from the initial image with a covariance-based approach. The algorithm is described in Faessel and Jeulin (2010): the authors use the covariance for estimating the average radius of the grains, and then estimate the number of grains in the image. CLOSING OF THE CRACKS The studied solid propellant has two phases: grains and matrix. Therefore, the first step of the reconstruction is to compute a binary mask for the grains. The threshold is estimated via the maximization of the interclass variance (Otsu, 1979). After the binarization, a morphological closing of the binary mask is used. The structuring element is a rhombicuboctahedron of radius 3, providing a good approximation of a sphere of small size on the digitized image. It offers a good compromise between performance and exactness. The size of the structuring element used is the same for the three images but depends on a subjective choice that is checked by visual inspection. Results with a rhombicuboctahedron of radius of radius 2 are also acceptable, but will lead to a few additional errors in the final reconstruction (the over-segmentation of a few grains). The remaining holes inside closed grains are removed with a volumic opening: all the connected components of the matrix with a small volume are removed. The threshold on the volume is the same for the three images but also depends on a subjective choice. This threshold is fixed to 113 voxels (corresponding to a sphere of radius 3). Results with a threshold corresponding to a sphere of radius 4 or 5 are also acceptable, but will lead to a few additional errors in the final reconstruction, generating the fusion of a few grains. The results of this first part of the segmentation are illustrated in Fig. 2 and Fig. 5. (A) Fig. 2. (A) 3D X-ray microtomographic image of a fragmented granular material (slice). (B) Binarization by maximization of the interclass variance (slice). Fig. 3. (A) Original image for MAT2, the reference material without any mechanical impact (slice). (B) Watershed segmentation using markers computed from the h-minima (slice). (A) (B) Fig. 5. (A) Morphological closing of the binary mask with a small rhombicuboctahedron (slice). (B) Removal of the holes with a volumic filter (slice). WATERSHED COMPUTED FROM THE MINIMA A first segmentation to reconstruct the grains is provided by markers on selected minima of the distance map. From the morphologically closed image, constructed in the section Closing of the cracks, a distance map is computed. Then, the inverse of this distance map is segmented with a watershed transform. The use of the local minima of the distance map will lead to an over-segmentation. Therefore, a prior h-minima transform is used on the distance map. The h-minima transform suppresses all minima whose depth is less than h, reducing the number of local minima. The difficulty is in the choice of h. Using a trial and error approach and a visual inspection, it is possible to achieve a good reconstruction (Fig. 6A). . li-^ÍÉÚ».-*»«. kt ** * (A) Fig. 4. (A) Original image for MAT3. An impact is obtained from a 2 kg mass falling 30 cm (slice). (B) Watershed segmentation using markers computed from the h-minima (slice). (A) (B) Fig. 6. (A) Watershed segmentation of the closed image after a h-minima filter (slice). (B) Watershed segmentation using markers computed from the K-means (slice). MARKERS COMPUTED FROM THE K -MEANS For comparison, an alternative segmentation for the reconstruction of grains is obtained from markers generated by the K-means clustering algorithm. A realization of Poisson point process is generated inside the binary mask of the image, with a low intensity. This gives a set of voxels used for a K-means clustering (Lloyd, 1982). Given a set of voxels {p1,p2,...,pn}, K-means clustering aims to partition the n voxels into K < n clusters, {S1,S2,...,SK} so as to minimize the following sum of squares in each cluster: Yf=1 LpjeSi |pj — Mi||2, where ^ is the center of mass of voxels in the cluster Si (assuming all the voxels have the same mass). There exists several algorithms for building the set of clusters minimizing this sum. Here, the MacQueen algorithm, as implemented in the software R, is used (MacQueen, 1967). The centers of mass ^ are finally used as markers for the watershed transform on the distance map of the morphologically closed image constructed in the section Closing of the cracks. The numbers K of classes used is directly estimated from the closed image. Considering for simplification the closed image as a realization of a Boolean model of spheres with a single radius R, it is deduced from the range of the covariance C(h) equal to 2R (Faessel and Jeulin, 2010), obtained by the distance h for which C(h) = Vv2. Then, the number of grains is deduced from this radius and from the Boolean model assumption (Jeulin, 1991): the volume fraction of the overlapping grains Vv is expressed as a function of the average number of grains per unit volume 6, by: Vv = 1 — exp(—64nR3) .6 and consequently the total number of grains (and of markers) is deduced from Vv and R3. The same covariance-based approach, applied on the original thresholded image, before the morphological closing, gives an estimation of the size of the fragments. From this estimation, the intensity of the Poisson point process used for the clustering algorithm is deduced. The intensity is chosen such that each fragments receives at least several points. The process achieves a satisfactory reconstruction of the grains, as illustrated in Fig. 6B. MULTISCALE IMAGE SEGMENTATION in this section, the goal is to isolate each fragment of each grain, for further morphological analysis. For the separation of the fragments, a multiscale stochastic watershed algorithm is used. The stochastic watershed segmentation was first introduced in Angulo and Jeulin (2007). The approach is based on using a large number of realizations of random markers to build a probability density function (pdf) of contours, starting from a standard watershed algorithm producing oversegmentation. The stochastic watershed was proved to be efficient for unsupervised segmentation (Noyel et al., 2007; Faessel and Jeulin, 2010). For multiscale images, the full granulometry of the image is used Gillibert and Jeulin (2011b). Using morphological openings, this granulometry can be automatically computed from the image and is used as a constraint during iterations of segmentation steps. STOCHASTIC WATERSHED The aim of the stochastic watershed Angulo and Jeulin (2007) is to estimate for each point of the contours of a standard watershed a probability (called here probability density function of contours) of detection from random markers. The first method introduced for computing the stochastic watershed is based on a large number of realizations of random markers to estimate a probability density function of contours, or of surface boundaries in 3D. The random markers are generated with a uniform distribution of their coordiantes. For the present composite material, a constant background marker is added to each set of random markers. This constant background marker is extracted by thresholding the image via the maximization of the interclass variance. For each set of markers, a constrained watershed transform is computed. Then, the Parzen window method (typically here a convolution of the probability image by a Gaussian kernel) is used to estimate the probability density function of contours from this finite set of random realizations (Parzen, 1962). A good estimation of the probability of contours generated by the stochastic watershed requires 100 to 200 realizations (Angulo and Jeulin, 2007). However, using A-flat zones to smoothen the local probability of contours, a stochastic watershed segmentation can be achieved with only 50 realizations (Faessel and Jeulin, 2010). Computing a large number of watershed transforms from simulations provides good results but is a slow process, mainly in 3D. A more efficient solution for computing stochastic watersheds is to use a graph-based approach. Probability of boundaries is directly computed with a good approximation without the use of any realization (Jeulin, 2008). As an example, the computation of the stochastic watershed with 50 realizations of watershed transforms takes 163 minutes and 20 seconds for the MAT2 sample on 3.00 GHz Pentium 4. A similar result can be achieved in 7 minutes an 8 seconds, on the same computer, using the graph-based approach. In Stawiaski and Meyer (2010) and Gillibert and Jeulin (2011b), the direct computation of the probability of the boundaries is obtained using a region adjacency graph deduced from the watershed, each vertex of the graph figuring a basin of attraction of the watershed, and each edge connecting two neigbouring basins. This graph-based approach leads to a multiscale stochastic watershed algorithm that is used now. MULTISCALE STOCHASTIC WATERSHED The main drawback of the stochastic watershed is that it is not well suited for the segmentation of objects with a wide range of scale. A variant was introduced by the authors to operate on a wide granulometric spectrum Gillibert and Jeulin (2011b). The multiscale image segmentation process is based on a simple idea: estimate the full granulometry of the image, using morphological openings, then use multiple stochastic watersheds with different numbers of markers for each size, and finally combine them to get a segmentation that is correct for each size of grains (no oversegmentation for big grains, no undersegmentation for small grains). Many hierarchical segmentation algorithms were studied, such as the waterfalls (Beucher, 1994) or the P algorithm (Beucher and Marcotegui, 2009). Here the approach is based on the merging of the watersheds basins using a minimum spanning tree (Eppstein, 2000). In the merging process, a constraint is introduced: the granulometry of the image. The algorithm is described in Fig. 7. Fig. 7. The main steps of the multiscale image segmentation process. The first step of the approach is to estimate the granulometry of the image, using morphological openings. From the granulometry, a small number of classes are chosen (3 classes in this paper). A good approach for this choice is to maximize the interclass variance. The total volume of the grains in each class x will be denoted v(x). The number of grains in each class x is deduced from v(x) and is denoted n(x). It is used to generate the corresponding number of markers in the calculation of the probability of the boundaries between grains of the segmentation. Then, the standard watershed transform is computed from the local minima of the gradient. From this watershed, the adjacency graph is constructed and a minimum spanning tree is extracted. The first class is chosen, starting from the largest grains. The stochastic watershed is computed with a number of markers equals to n(1). Based on this stochastic watershed, a first hierarchy on boundaries is computed with the merging algorithm. For each step i of the hierarchy, the granulometry of the corresponding segmentation is computed (vi(1)). In the full hierarchy, there is a size step which minimizes the difference |v;(1) — v(1)|. This step is used for the segmentation of the grains in the first class. All the segmented grains are removed from the image and added to the background mask. The minimum spanning tree is updated and the next class is chosen. The same process is applied to all the classes. When no more class is left, all the segmentations are combined together. Results are illustrated in Fig. 8. A few fragments are missing, and a few grains are oversegmented, but this errors have a small impact on the results of the measurements. (A) (B) Fig. 8. (A) Original image before multiscale segmentation (slice). (B) Final multiscale stochastic watershed segmentation (slice). IMAGE ANALYSIS AND MEASUREMENTS VOLUME FRACTION AND SPECIFIC SURFACE AREA OF DAMAGED GRAINS From the reconstructed data, the volume fraction of each grain is estimated before and after the morphological closing. From this two measurement, the volume fraction of the cracks of each grain is estimated. For each image, grains reconstructed with the h-minima markers and grains reconstructed with the K-means markers are studied for comparison. The results for the volume fraction of the cracks, presented as normalized histograms (namely the proportion of grains in % with a given property), are shown in Figs. 9, 10 and 11. The agreement between the K-means reconstruction results and the h-minima reconstruction results is excellent, showing the robustness of our segmentations for the purpose of damage measurement. With the number of intercepts (transitions from background to foreground) in 13 directions generated by a voxel and its first and second neighbours on the cubic grid, it is possible to estimate the surface area of each grain i from the closed reconstructed data (denoted Sc(i)). With the same process, the surface area of each grain before the morphological closing is estimated (denoted Sf (i)). From this two measurement, the surface area S(i) of the cracks of each grain is estimated: S(i) = Sf (i) - Sc(i) . Due to some imperfections on the original grains, the surface area estimated with this process correspond to the surface area of the cracks, the small irregularities at the surface of the grains and the porosity. The specific surface area Sspec measures the total surface area per unit of volume: Sspec(i) Si V (i) The results for the specific surface area, presented as a normalized histogram, are shown in Fig. 12, Figs 13 and 14. Once again, the agreement between the K-means reconstruction results and the h-minima reconstruction results is excellent. The damage on the grains is quantified from both specific surface area and volume fraction. On the sample MAT2, the reference material without any mechanical impact, the volume fraction of the cracks is low (mostly between 0 and 0.2) and the volume fraction is almost always 0. On the sample MAT3, impacted with a 2 kg mass falling 30 cm, with a 9.5 cm rebound, there is only 6% of the grains with a zero specific surface area of the cracks, and the volume fraction of the cracks is mostly between 0.2 and 0.5. The sample MAT1, impacted with a 2 kg mass falling 15 cm, is less damaged than the sample MAT3. This is visible on both specific surface area and volume fraction. On the sample MAT1 there is only 8% of the grains with a zero specific surface area of the cracks, and the volume fraction of the cracks is mostly between 0. 1 and 0. 4. Volume fraction of the cracks (histogram) MAT1 (HM) - MAT1 (KM) 1 '' \ ft A ¡1 \ ■i \ \ \ y \j—>T\' / v v\ /KaX Volume fraction of the cracks Fig. 9. Histogram of the volume fraction of the cracks for MAT1. HM are the grains reconstructed with the h-minima markers and KM are the grains reconstructed with the K-means markers. 4 0 Volume fraction of the cracks (histogram) Fig. 10. Histogram of the volume fraction of the cracks for MAT2. HM are the grains reconstructed with the h-minima markers and KM are the grains reconstructed with the K-means markers. Specific surface area of the cracks (histogram) Fig. 13. Histogram of the specific surface area of the cracks for MAT2 (given in ^m-1). HM are the grains reconstructed with the h-minima markers and KM are the grains reconstructed with the K-means markers. Volume fraction of the cracks (histogram) Fig. 11. Histogram of the volume fraction of the cracks for MAT3. HM are the grains reconstructed with the h-minima markers and KM are the grains reconstructed with the K-means markers. Specific surface area of the cracks (histogram) Fig. 14. Histogram of the specific surface area of the cracks for MAT3 (given in ^m-1). HM are the grains reconstructed with the h-minima markers and KM are the grains reconstructed with the K-means markers. Number of fragments (histogram) Specific surface area of the cracks (histogram) i MAT1 (HM) — MAT1 (KM) .... : 1 V>ôo< 0.02 0.03 0.05 0.06 Specific surface area of the cracks Fig. 12. Histogram of the specific surface area of the cracks for MAT1 (given in ^m-1). HM are the grains reconstructed with the h-minima markers and KM are the grains reconstructed with the K-means markers. MAT2 (HM) MAT2 (KM) MAT3 (HM) MAT3 (KM) 30 40 Number of fragments Fig. 15. Histogram of the number of fragments for MAT2 and MAT3. HM are the grains reconstructed with the h-minima markers and KM are the grains reconstructed with the K-means markers. 4 NUMBER OF FRAGMENTS With the multiscale stochastic watershed algorithm introduced in the the previous section, the separation of the fragments can be achieved. Using the labels computed with the reconstruction algorithm, each segmented fragment is associated to its original grain. Therefore, it is possible to know the number of fragments in each grain. This number includes the grain itself and is therefore larger or equal to 1. From this number, a normalized histogram is estimated, as illustrated in Fig. 15 for samples MAT2 and MAT3. For the fragmented grains (MAT3), a three scale stochastic watershed is used. For the reference material without any mechanical impact, a simple stochastic watershed, without any additional steps, gives a correct segmentation. As for the specific surface area and volume fraction, the damage of the grains can be quantified from these normalized histograms, the initial material showing almost no fragmentation, as compared to the shocked material, MAT3. The results obtained from the two types of segmentation are pretty close, showing again the robustness of the used segmentation techniques. We can observe from the normalized histogram that many grains are very fragmented (> 30 fragments per grain, as seen in Fig. 15). CONCLUSION The proposed algorithm gives a satisfactory reconstruction of the fragmented grains with both markers sets. Visual inspection reveals that the K-means markers give better results when the grains are highly fragmented and if the fragments are scattered, but the boundaries of the grains are less accurate. The h-minima markers give correct boundaries between grains but fail to reconstruct a few grains. in both cases, the closing of the grains depends on two parameters: the radius of the rhombicuboctahedron used and the volume of the holes. The h-minima markers require an additional parameter: h. The number of clusters used for the computation of the K-means is directly deduced from the image. Therefore, the K-means depend on less parameters and is less user-dependent. The two algorithms require a similar time for the reconstruction. The reconstruction of the sample MAT2, using the K-means approach, requires 18 minutes and 12 seconds on a 3.00 GHz Pentium 4. The reconstruction of the same sample requires 19 minutes and 15 seconds using the h-minima approach. Both approaches are very useful in order to generate data on the damage of materials from a fully automated 3D image analysis. The reproducibility of results on fragmentation obtained from two segmentation methods is very good. The analysis proves that the methods provide 3D damage measurements consistent with the mechanical impacts applied on the materials. They will provide useful information for the fully automatic characterization of damage in various conditions, helping to improve the reliability of solid propellants. A similar approach can be followed for the quantitative analysis of the progression of damage in materials, starting from 3D images. ACKNOWLEDGMENTS This work was supported by a grant from DGA (contract 2009 34 0006). The authors are grateful to Alain Fanget (CEA Gramat) for his advice during this study. REFERENCES Angulo J, Jeulin D (2007). Stochastic watershed segmentation. In: Banon GJF, Barrera J, Braga-Neto U, Hirata NST, eds. Proc 8th Int Symp Math Morpho (ISMM) 1:265-76. Beucher S (1994). Watershed, hierarchical segmentation and waterfall algorithm. In: Serra J, Soille P, eds. Mathematical morphology and its applications to image processing. Computational Imaging and Vision 2:6976. Beucher S, Lantuejoul C (1979). Use of watersheds in contour detection. In: Proc Int Worksh Image Process Real-Time Edge Motion Detect. Sept 17-21. Rennes, France. Beucher S, Marcotegui B (2009). P algorithm, a dramatic enhancement of the waterfall transformation. Tech Rep. CMM/Mines Paristech. Eppstein D (2000). Spanning trees and spanners. In: Sack JR, Urrutia J, eds. Handbook of Computational Geometry. Chap. 9, pp. 425-61. Elsevier. Faessel M, Jeulin D (2010). Segmentation of 3D microtomographic images of granular materials with the stochastic watershed. J Microsc 239:17-31. Gillibert L, Jeulin D (2011a). 3D reconstruction of fragmented granular materials. In: Proc 13th Int Cong Stereol (ICS-13). Beijing, China. Gillibert L, Jeulin D (2011b). Stochastic multiscale segmentation constrained by image content. In: Soille P, Pesaresi M, Ouzounis GK, eds. Proc 10th Int Conf Math Morpho Appl Image Signal Process (ISMM'11). Lect Not Comput Sci 6671:132-42. Jeulin D (1991). Modèles morphologiques de structures aléatoires et de changement d'echelle. PhD Thesis. University of Caen, France. Jeulin D (2008). Remarques sur la segmentation probabiliste. Tech. Rep. N-10/08/MM, CMM/Mines Paristech. Lloyd SP (1982). Least squares quantization in PCM. IEEE Trans Inf Theory 28:129-37. MacQueen JB (1967). Some methods for classification and analysis of multivariate observations. In: Cam LML, Neyman J, eds. Proc 5th Berkeley Symp Math Stat Probab, vol. 1. Berkeley: University of California Press. Matheron G (1967). Elements pour une theorie des milieux poreux. Paris: Masson. Noyel G, Angulo J, Jeulin D (2007). Random germs and stochastic watershed for unsupervised multispectral image segmentation. In: Apolloni B, Howlett RJ, Jain L, eds. Proc Knowledge-Based Intell Inf Eng Syst (KES). Lect Not Comput Sci 4694:17-24. Otsu N (1979). A threshold selection method from gray-level histograms. IEEE T Syst Man Cyb 9:62-6. Parzen E (1962). On estimation of a probability density function and mode. Ann Math Stat 33:1065-76. Serra J (1982). Image analysis and mathematical morphology. London: Academic Press. Soille P (2003). Morphological image analysis: principles and applications. Secaucus, NJ, USA: Springer. 2nd ed. Stawiaski J, Meyer F (2010). Stochastic watershed on graphs and hierarchical segmentation. In: Proc Eur Conf Math Indust. Wuppertal, Germany.