Modern computer vision and machine learning techniques, when applied in Fractogra- phy bare the potential to automate much of the failure analysis process and remove human induced ambiguity or bias. Given the complex interaction between intrinsic (e.g. microstructure) and extrinsic (e.g. environment, loading history) factors leading to failure, deep learning methods, which exhibit very high eciency in establishing complex interconnections between the input data, may end up revealing new correla- tions and information that is encoded onto the complex geometries of fracture surfaces and remained hidden from us so far. In this work, we examine the potential use of an unsupervised learning pipeline to classify fracture surfaces of ve tungsten heavy alloys following their chemical content (i.e. Tungsten percentage). Encouraged by the success of the algorithms, we move on and analyze the features on the fracture surfaces which are governing the decision process of the algorithms. The fractographic interpretation of these features shows that the extent of plasticity on the fracture surface serves as a measure for the classi cation process. The examined pipeline can be used to identify failures originating from erroneous manufacturing processes, leading to locally varying Tungsten concentrations and ultimately premature failure. Corresponding author. Tel. +30-6932913594; E-mail address: tsopanidisstelios@gmail.com Preprint submitted to Elsevier May 8, 2020 Keywords: Quantitative Fractography, Convolutional Neural Networks, t-SNE, Unsupervised Learning, Data Clustering, Machine Learning Interpretability 1. Introduction Materials science has heavily relied on visual inspection since its early days. The analysis of microstructural data, obtained through various microscopy techniques has become a corner stone in materials science and engineering, allowing to establish process- microstructure-properties correlations. Similarly, fractographic examination of frac- tured specimens and components is extensively used to understand the mechanisms driving the fracture process and their relation with the operating conditions and un- derlying microstructure. It is no surprise then that the material science community has readily adopted emerging technologies in the eld of computer vision, and more speci cally recent advances in Machine Learning (ML). Mainly, these approaches im- plement supervised learning methods, where a neural network, after being trained on a large dataset of images (usually obtained by di erent microscopy methods), performs predictions on a test dataset. The focus of the majority of published works in this eld, is on the identi cation and classi cation of the microstructure of di erent material sys- tems, aiming to establish microstructure-property linkages or serve as a tool for quality 2 assurance (DeCost et al., 2017). Kondo et al. (2017) used seven ceramic micrographs, synthesized with varying sintering temperature and time, to train a Convolutional Neu- ral Network (CNN) aimed at identifying the corresponding ionic conductivities. In Az- imi et al. (2018), a semantic segmentation CNN is trained in extracting and classifying microstructural constituents of low carbon steel in Scanning Electron and Light Optical Microscopy images, while DeCost et al. (2019) implement a similar pixelwise CNN to segment four principal micro-constituents in a heat-treated Utrahigh Carbon Steel mi- crostructure dataset. Chowdhury et al. (2016) implemented di erent machine learning methods (visual bag of words, texture and shape statistics, and pretrained CNNs) to extract features from micrographs that depict dendritic morphologies and consequently classify them with support vector machine, voting, nearest neighbors and random forest models. Recently, a signi cant contribution to the eld of crystallography was made by Kaufmann et al. (2019), where a collection of Electron BackScattered Di raction (EBSD) patterns with characteristic crystal structures was used to train a machine learning algorithm in identifying the origin Bravais lattice or point group of a sample from raw EBSD data. While the papers cited so far are only a small portion of the pub- lished literature where machine learning techniques are employed for microstructural characterization, much less progress was done in the eld of fractography. The major- ity of published works where computer vision techniques are employed for automated fractography are based on classical image processing techniques, e.g. edge detection, texture analysis etc. To the authors best knowledge, Bastidas-Rodriguez et al. (2020); Konovalenko et al. (2018); Kitahara and Holm (2018) and Tsopanidis et al. (2020) , are the only published works where ML techniques have been successfully applied for frac- tographic analysis which goes beyond feature enhancement (e.g. automatic extraction of dimple size). Although supervised machine learning approaches are proven very ecient, the requirement of an extensive and time-costly training process constitutes a drawback and the reliance of the annotation process on human input can potentially introduce bias or errors. Furthermore, these algorithms usually fail in performing predictions 3 on microscopy images of material systems with morphological and micro-structural characteristics that deviate substantially from the training dataset. Thus, recent pub- lications propose unsupervised learning approaches, which combined with pre-trained CNNs ( transfer learning ) intent to overcome this limitation and enable inference on microscopy images from broader range of material systems; or even achieve mate- rial independent performance. The methods introduced in Lubbers et al. (2017), in which the Gram matrices, constructed from the activation maps of the hidden layers of a pre-trained CNN, are used for microstructural characterization, appear to be very promising. Finally, Kitahara and Holm (2018) propose an unsupervised learning algo- rithm that implements a pre-trained CNN to extract features from two di erent image datasets of metal alloys and combines it with dimensionality reduction algorithms to group them into clusters according to their visual similarities. Most of the aforemen- tioned publications present machine learning methods that exhibit very high accuracy, but only few exceptions (Kondo et al., 2017; Kaufmann et al., 2019) go one step further and investigate the internal operations that enable this high ecacy or interpret the performance of the algorithm. Kondo et al. (2017) is using a variation of Class Activa- tion Maps (CAM) (Zhou et al., 2016) and Grad-CAM (Selvaraju et al., 2017) methods to visualize the locations that are activated by the intermediate layers of the neural network and provide a physical interpretation of the results. Similarly, Kaufmann et al. (2019) compute the Grad-CAM visualizations and their interpretation shows that the neural network bases its predictions on the same features a crystallographer would use to identify the crystal structure. Despite the great promise and proven successes of machine learning techniques listed above, some degree of skepticism is often met when discussing their contribution and future perspectives with scientists from the fracture mechanics community. Many voices are heard, calling for a more rigorous evaluation of the proposed methods and better interpretation and understanding of the internal functions that enable the high ecacy of these algorithms (Holm, 2019). The second, perhaps more alluring question, relate to the extent to which machine learning can add to our current understanding of fracture 4 processes and to what degree it is over-hyped or truly a new tool in our arsenal that can lead us to new, exciting revelations. Thus, two main questions regarding machine learning methods need to be answered: {What are they good for? In computer science and information theory, machine learning methods have shown very promising results and their eciency is proven, however in fracture mechanics only few very recent publications have started exploring their potential and the results reported, although very encouraging, require more extensive research as to their advantages and limitations. An im- portant question in this regard is whether machine learning is only useful to make previously tedious tasks more readily achievable (e.g. Tsopanidis et al. (2020) ) or it can allow us to go beyond what is conceived feasible today. Consider for example the on going quest to correlate fracture surface roughness with fracture toughness (Barak et al., 2019; Mandelbrot et al., 1984; Srivastava et al., 2014), can new correlations be found? What other information could be hidden in the complex geometry of fracture surfaces we have not even considered looking for so far? {How do they work? A better understanding of the functionality of the Machine Learning methods implemented in Fracture Mechanics will enable more ecient evaluation of their performance and allow the research community to explore new applications and realize their full potential. To some extent machine learning methods have been treated as a "black box" and this approach does not allow the researchers to properly evaluate their applicability and creates a reasonable incertitude. The inferability of machine learning techniques, given that new cor- relations can be found, are the basis for nding causalities and expanding our knowledge as to the mechanistic origin of those presumed correlation. The work presented here is aimed at scratching the surface of the two questions posed above through a simple example. We will demonstrate that while the topography of 5 a fracture surface is traditionally considered a poor indicator to its chemical content, a data pipeline can be constructed to assist in clustering fracture surfaces of tungsten heavy alloys (WHA), following their tungsten content. More over, we identify the features on which the used machine learning algorithm is basing it predictions. Our research builds upon the work of Kitahara and Holm (2018) and presents two algorithms that enable clustering and classi cation of an SEM fracture images dataset, according to the composition of the alloy material samples that are used for acquir- ing the fracture surfaces. The rst part of this data clustering pipeline constitutes a CNN, which is using pre-trained set of weights to extract meaningful features from the input fracture images, and the second part is a sequence of dimensionality reduction algorithms that cluster the images according to the similarity of the extracted features. Subsequently, this pipeline is extended with the addition of a minimally supervised algorithm, which enables the classi cation of the clustered data points. Finally, the visualization of the activation maps of the last convolution layer of the neural network enables the analysis of the internal operations that enable the pipeline to cluster the fracture images and the interpretation of the neural network's functionality. 2. Material and Method 2.1. Dataset Five di erent samples of tungsten heavy alloys (WHA), with di erent composition of tungsten (90wt%, 92wt%, 95wt%, 97wt% and 99wt%) were produced. The WHA at tensile specimens at varying tungsten density were loaded at a strain rate of 1000 s1, using a split Hopkinson tension bar apparatus. Subsequently, the fracture surfaces were scanned using a Tescan MIRA3 FEG-SEM and a dataset of 10 SEM images with size of 40964096 pixels was created. This dataset contains 2 SEM images for each tungsten composition category (each image taken from a di erent specimen). Finally, the images were cropped to the size of 448 448 pixels and the nal dataset is composed by 810 SEM images; 162 images for each tungsten composition. 6 90% 92% 95% 97% 99% Figure 1: SEM fracture images from WHA samples with di erent mass percentage of tungsten. The images in each row belong to samples with the same tungsten mass percentage. Fracture in WHA is mainly attributed to W-W boundary brittle fracture, W-matrix de-cohesion, ductile matrix rupture and W grain to grain cleavage. As a result, one fracture surface contains many fracture mechanisms within it, thus making it com- plicated to infer. The fracture surfaces are very similar and even for an experienced 7 fractographer will be very dicult to realize any di erences. In Fig. 1, four randomly selected images of size 448 448 pixels from each sample group are presented. With exception of the samples with 99% mass percentage of tungsten, the rest of the SEM images exhibit very similar morphology. 2.2. Clustering data pipeline Following the Kitahara and Holm (2018) approach, a clustering algorithm is devel- oped. The main objective of this algorithm is to group the dataset of WHA fracture images into di erent clusters according to the tungsten mass composition of the sample they belong to. Figure 2: Schematic owchart of the clustering data pipeline. The clustering algorithm is a 4-part data pipeline, where the dimensionality of the input data is reduced after each stage. A schematic owchart of the algorithm is presented in Fig. 2 and the function of each stage is de ned as follows: VGG16 : the VGG16 neural network architecture (Simonyan and Zisserman, 2014), where the nal Fully-Connected layers are replaced by a Global Average Pooling (GAP) layer, is used to extract features from the input SEM fracture images. The network, using a set of pre-trained weights on the ImageNet dataset (Russakovsky et al., 2015), performs predictions on the SEM images. Each frac- ture image of the dataset is imported into the network and, after passing through all the consecutive layers, the last convolution layer outputs a sequence of 512 two-dimensional feature maps for each image. The GAP layer at the end of the 8 network computes the mean value for each feature map and exports a feature vector in R512for each input image. This modi ed network architecture enables the extraction of 512 features from each input image. Thus, a set of weights, although being trained on an extensive dataset of images that have no relevance to fracture images, when loaded to the modi ed VGG16 architecture allows the network to reduce the dimensionality of the input images from 224 224 to 512 dimensions and extract meaningful information from the input data; the fracture images were resized to 224 224, before inserted into the VGG. PCA : a Principal Component Analysis (Karl Pearson F.R.S., 1901) on the feature vectors extracted by VGG allows the computation of the 50 principal components that maximize the total variance of the data. This linear-dimensionality reduction algorithm expresses each feature vector, which corresponds to a certain fracture image, with respect to 50 principal axes that maximize the total variance of the data. The algorithm enables the estimation of the "importance" of each principal component, by computing the eigenvalues of the covariance matrix of the set of feature vectors and furthermore the computation of the weight that each feature extracted by VGG has on each principal axis, by computing the eigenvectors of the covariance matrix. t-SNE : this algorithm (van der Maaten and Hinton, 2008) performs a non-linear dimensionality reduction on the set of feature vectors in R50that is exported by the PCA. It computes similarities between the data points in the 50-dimensional space and projects them onto the 2D space, according to these similarities. The nal output of this stage is a 2D plot. By the end of this third part, the pipeline achieves to project the initial fracture images dataset onto data points on a 2D plot, in positions that enable clustering according to similarities of the fracture images. k-Means : the k-Means algorithm (Lloyd, 1982) groups the data points, exported by t-SNE into 5 clusters, according to their Euclidean distances in the 2D space. 9 As a result, k-Means assigns a label to each data point according to the cluster that it belongs. This 4-part data pipeline is an unsupervised machine learning algorithm that groups any dataset into a prede ned number of clusters, without requiring any training process; even the rst part of the pipeline that is a neural network does not require additional training, since it implements the transfer learning technique. The computer code for this algorithm is developed in Python with the implementation of Keras (Chollet et al., 2015) and scikit-learn (Pedregosa et al., 2011) libraries; Keras is used for building the VGG16 model and subsequently the inference operation of the rst part of the pipeline, while the scikit-learn functions are implemented for the other parts of the pipeline. The nal output is a 2D scatter plot where the position of each data point is de ned according to similarities on the high-dimensional space of the input dataset. The performance of this algorithm on the dataset of the WHA fracture image is evaluated in the next section. 2.3. Classi cation data pipeline The previous clustering pipeline can be extended into a classi cation algorithm, with the addition of a k-Nearest Neighbors (KNN) algorithm (Cover and Hart, 2006) at the end of the pipeline. The k-Nearest Neighbors is referred to as a minimal supervision algorithm because the training process only involves saving the training data points and the corresponding labels. During the inference operation, the algorithm calculates theEuclidean distances between each test data point and all the training data points, which enables the computation of knearest neighbors in the training dataset for each test data point. Finally, KNN assigns to each test data point the label that has the most representatives among the kneighbors. The de nition of a computational framework that combines the training and in- ference operations of the algorithm enables the classi cation of the fracture images according to the 5 tungsten composition labels. A schematic owchart of this compu- tational framework is presented in Fig. 3 and the outline of this computation process 10 Figure 3: Schematic owchart of the computation framework of the classi cation algorithm. follows the next steps: 1. Initially, 90% of the WHA fracture images dataset is set as training images, while the rest 10% is set as test images. 2. The entire dataset is imported into the VGG -PCA -t-SNE pipeline and the result is a 2D scatter plot. This plot is separated into a plot that contains only the training data points and another one with the test data points. 3. The training data points, projected onto the speci c locations on the 2D plot by the VGG -PCA -t-SNE pipeline, are imported into another pipeline con- structed by k-Means and KNN algorithms. 4. A mesh grid with dimensions large enough to accommodate every training and test data point is created and KNN enables the classi cation of every grid point into one of the 5 tungsten composition labels, using the training data points and their annotations. The result of this step is a colormap, where each area is assigned to a di erent tungsten composition label. 5. The nal step involves plotting the test data points, with the positions predicted by the VGG -PCA -t-SNE pipeline, onto the colormap. The label of the area that each test data point is placed de nes the classi cation of the test point. This computation process enables the classi cation of a fracture images test dataset according to pre-de ned labels, when a training dataset with the corresponding anno- tations are provided. 11 3. Results and discussion 3.1. Clustering of fracture images The clustering data pipeline, presented in the previous section, is implemented on the dataset of WHA fracture images and the results are presented in this section. Di er- ent parameters for the intermediate algorithms of the pipeline have been investigated. More speci cally, the PCA algorithm was tested with 50 or 100 principal components and t-SNE was implemented with Barnes-Hut approximation for the gradient compu- tation (Van Der Maaten, 2013), and it was tested for perplexity values in the range of (5, 50), learning rate in range (100, 1000) and the number of iterations was varied between 1000 to 10000 iterations. Finally, k-Means algorithm was tested with ran- dom ork-means++ centroid initialization; kmeans++ is an initialization scheme for the cluster centroids, which has been implemented in scikit-learn according to Arthur and Vassilvitskii (2006). The combination of the parameter values that resulted into the best performance for this clustering pipeline is presented in Table 1. Algorithm Parameters PCA Number of principal components = 50 t-SNE Perplexity = 40, Learning rate = 200, Iterations = 3000 k-Means Centroid initialization = k-means++ Table 1: Parameters of the clustering pipeline that enable the highest eciency clustering of the WHA fracture images dataset. The dataset of 810 SEM images, obtained by fracture experiments on the WHA samples with di erent tungsten composition, is imported to the clustering data pipeline and the resulted 2D scatter plot is shown in Fig. 4. The initial qualitative assessment of the result shows that the algorithm succeeds in creating three distinct clusters (green, blue and red data points in the plot), while k- Means manages to separate the main cluster in the middle by assigning di erent labels 12 Figure 4: Clustering results on the WHA fracture images dataset. to the data points. Subsequently, in order to quantify the performance of the algorithm the data points are plotted with their ground truth labels (Fig. 5). The ground truth plot shows that the algorithm fails to cluster correctly the data points that belong to the rst SEM fracture image of the sample with 90wt% tungsten composition; for each tungsten composition sample two SEM images were obtained and then cropped before added to the dataset, as it was explained in the beginning of the previous section. The erroneous predictions are marked (see crossed data points in Fig. 6) and the accuracy of the algorithm is computed to 81.2%. Additionally, the confusion matrix presented in Fig. 7 veri es that the main source of clustering errors is the data points that originate from the rst SEM fracture image of the sample with 90wt% tungsten composition. In numerous tests on the same dataset, the accuracy of the algorithm in clustering the WHA fracture images was never exhibited a value below 73%. Subsequently, the dataset images that belong to the rst SEM image of the sample with tungsten composition of 90wt% are removed and the accuracy of the algorithm in 13 Figure 5: Ground truth of the clustering results on the WHA fracture images dataset. Figure 6: Clustering results with errors in predictions marked by crossed data points. 14 Figure 7: Confusion matrix of the clustering results shows that the rst SEM image from the sample with 90wt% tungsten composition is the main source of clustering errors. Figure 8: Clustering results, after removing the data points of the rst SEM image of the sample with 90wt% tungsten composition, with errors in predictions marked by crossed data points. clustering the remaining dataset images is computed. In Fig. 8 the clustering result is presented, with the errors being crossed as before. The improvement is signi cant and the confusion matrix (Fig. 9) shows a uniform distribution of errors. Finally, the 15 accuracy is raised to 91%, after removing the data from this particular SEM image that the algorithm failed to cluster correctly. Figure 9: Confusion matrix of the clustering results, after removing the data points of the rst SEM image of the sample with 90wt% tungsten composition. The visual di erences between the two SEM images of the sample with 90wt% tungsten composition are obvious, as it can be seen in Fig. 10, and this explains the failure of the algorithm in grouping the data points that originated from these two images in the same cluster. This observation initiated the inquiry of the internal functions that enable the clustering operation of the algorithm. The question that arises is: {Which features of the fracture images "attract the attention" of the neural net- work and from these extracted features which weigh more in the nal result of the dimensionality reduction algorithms? The answer to this question is the subject of the last part of this section, which focuses on the visualization of the activation maps of the neural network and the inter- pretation of the functionality of the clustering algorithm that enables its ecacy. 16 (a) A102: rst SEM image of the 90wt% tungsten compo- sition sample. (b) A103: second SEM image of the 90wt% tungsten com- position sample. Figure 10: The two SEM images obtained after scanning the sample with the 90wt% tungsten com- position. Note that both images are taking under the same magni cation and eld of view 3.2. Fracture images classi cation The WHA fracture images dataset, without including the dataset images from the A102 SEM image (in Fig. 10a), is imported into the classi cation algorithm. 675 images of this dataset are assigned to training and the rest 54 are considered as the test images. Following the computational framework for the classi cation of the dataset, the algorithm is plotting the data point that correspond to the test images onto a colormap that de nes the labels for the di erent tungsten compositions. For the k- Nearest Neighbors algorithm in the classi cation data pipeline, the number of nearest neighbors is set to k= 15 and uniform weight distribution is used. The result of this classi cation process is presented in Fig. 11a. The position that each test data point on the colormap, which is computed by the classi cation data pipeline, de nes the label that is assigned to the data point. Knowing the ground truth 17 (a) The result of the classi cation for the WHA dataset. (b) The ground truth labels of the classi cation result. Figure 11: Classi cation results and ground truth annotations, with classi cation errors inside the red circles. annotations for the test data points allows the evaluation of the classi cation errors (data points inside the red circles in Fig. 11b) and enables the computation of the classi cation accuracy. The accuracy of the classi cation predictions is computed to 93%, while in numerous tests performed on the same dataset the lowest accuracy was 85%. 3.3. Comparison with Haralick texture descriptors method The evaluation of the clustering and classi cation results shows that the two al- gorithms presented here exhibit very high accuracy when performing on a dataset of fracture images with very similar structural and morphological characteristics. In order to have a reference and better understanding of the algorithm eciency, the WHA dataset is imported into an algorithm that extracts di erent texture descrip- tors and consequently enable image classi cation, following the method proposed by Haralick et al. (1973). The central component of this statistical methodology, which is commonly used for images classi cation tasks where texture is of importance (e.g. fractographs), is the Gray Level Cooccurrence Matrix (GLCM). The Haralick method 18 de nes 14 texture descriptors, which after being trained on a training dataset, enable classi cation on the fracture images of a test dataset. The performance of the Haralick texture descriptors algorithm on the WHA fracture images dataset is very poor, with classi cation accuracy of 21%. Practically, this evalu- ation shows that the algorithm, which is based on a well-established method commonly used for image classi cation, fails to classify the fracture images of the WHA dataset and reveals the complexity of the task at hand. Additionally, it shows the great po- tential of the Clustering and Classi cation algorithms that are presented in this paper and makes more evident the importance of their high eciency performance. 3.4. Features visualization and Interpretability The objective of this section is to understand the internal functions that enable the clustering, and as an extension the classi cation, pipeline in performing with the reported high ecacy on the task at hand. More speci cally, the study focuses on the visualization and consequently the interpretation of the features that the algorithm is identifying in the input fracture images and through this information aims to gain a better understanding on how these extracted features enable accurate predictions. Considering the architecture of the VGG16 network, the output of the last convo- lution layer is a set of 512 activation maps, with dimension of 7 7, for each input fracture image. This set is denoted as: fAc kjk= 1;2;:::M andc= 1;2;:::Ng, where N= 810 is the total number of images in the dataset and M= 512 the number of activation maps. When this set of activation (or else feature) maps is inserted into the Global Average Pooling (GAP) layer, it exports a vector ( ~ ac2R512) with the mean values of the extracted activation maps from each image. This set of feature vectors, f~ acjc= 1;2;:::Ng, is imported into PCA. Initially, the PCA algorithm creates a data matrix X, with dimensions NM, by placing the feature vectors ~ acin the rows of the matrix, and then computes the covariance matrix S, with dimensions MM, as follows: 19 S=XTX Finally, the computation of the eigenvalues ( k) and the eigenvectors ( ~ uk) of the covariance matrix are computed with the Singular Value Decomposition (SVD) method. The eigenvalues express the total variance of the linear combinations of the feature vectors (~ ac) with the computed eigenvectors ( ~ uk), which are called Principal Compo- nents and their mathematical de nition is: ~Pk=X~ uk This enables the representation of the 512-dimensional dataset of the feature vec- tors with respect to the p= 50 principal components with the highest eigenvalues. Furthermore, the eigenvalues express the importance of each principal component on the reduced dimensionality representation and the eigenvectors that correspond to each principal components express the importance of each of the 512 extracted features by the VGG network on the speci c component. Thus, selecting the rst principal component it becomes possible to weigh its in uence in the reduced dimensionality dataset, through the computation of the corresponding eigenvalue, and even assign a weight to each feature extracted by VGG on this speci c component, by computing the 512 coecients of the corresponding eigenvector. The coecients of the corresponding eigenvector are the weights of the VGG extracted features. The computation of the eigenvalues from the PCA algorithm of the clustering pipeline on the WHA dataset shows that the rst principal component expresses 41.6% of the information in the reduced dimensionality representation of the dataset. Note that the second and the third principal components express 16.5% and 8.7% of the information, respectively, while the in uence of the rest of the principal components is even less signi cant. Following this analysis becomes apparent that for each input fracture image cof the WHA dataset, the weighted, by the coecients of the eigenvector of the rst principal 20 component, sum of all the activation maps exported by the nal convolution layer will show the exact features of the input image that are utilized by the clustering algorithm in the process of grouping this image with the rest of the images in the WHA fracture images dataset. This conclusion is mathematically de ned as follows: Ac heatmap =512X k=1jwkjAc k where the setfAc kjk= 1;2;:::512 andc=constantgis the set of the activation maps for the speci c input image and jwkjare the absolute values of the coecients of the eigenvector ~ uthat corresponds to the rst principal component of the PCA: ~ u=2 6666666666664w1 w2    w5123 7777777777775 The resulted activation heatmap Ac heatmap has dimensions of 7 7 and when resized to the original dimensions of the input image (448 448 pixels) enables the localization of the features that the clustering algorithm is identifying in the speci c image cof the WHA dataset. Implementing this computational framework in a python algorithm enables the vi- sualization of the set of activation heatmaps: fAc heatmapjc= 1;2;:::Ng, computed for each fracture image of the WHA dataset. Fig. 12 presents some representative activation heatmaps from each group of tungsten composition, while the entire set of activation heatmaps, which counts 810 images, together with the corresponding WHA fracture images dataset are published in Materials Data Facility (MDF) with DOI: https://doi.org/10.18126/aph0-olbz . 21 Figure 12: Activation heatmaps for a representative sub-set of the WHA fracture dataset. 22 A careful examination of the activation heatmaps concludes that the algorithm is focusing on the ductile fracture of the alloy material, since it activates the areas between the tungsten grains that accommodate deep and clear dimples. This ductile behavior is attributed to the fracture of nickel-Ferrous (NiFe) matrix and it can be hypothesized that the identi cation of the ductility of the fracture surface enables the algorithm to cluster the images according to the composition of nickel in the alloy material. Apparently, the algorithm does not focus on the tungsten's fracture and this very interesting observation can be justi ed by the fact that the fracture of tungsten can have di erent manifestations (intergranular or transgranular micro-fracture modes) in the alloy's fracture surface, which makes the identi cation more complex. This hypothesis provides a very reasonable interpretation of the internal operation of the clustering algorithm that enables the accurate clustering and elucidates the algorithm's ecacy in identifying the di erent fracture images. Although this analysis interprets the operation of the algorithm and the internal processes that enable the high eciency clustering, it does not provide an explanation on the fact that the algorithm fails to group the two SEM fracture images from the sample with 90wt% tungsten composition. In order to address this issue, the activation heatmaps from the fracture images, cropped from the two original SEM images obtained by scanning this particular sample (see Fig. 10), are used to reconstruct the activation heatmaps that correspond to these two initial SEM images. In Fig. 13 the activation heatmaps are shown. It is evident from these images that the algorithm has activated fewer areas in the activation heatmap of the image A102 compared to the activated areas in the heatmap of the A103 image. Following the previous analysis, the algorithm considers that the surface area of image A102 presents less ductility than the surface in the image A103. It is clear that a WHA sample produced with a speci c tungsten composition - in this case 90wt% - does not provide a fracture surface with uniform tungsten composition and it can be perfectly reasonable that the di erent areas that have been scanned by the SEM contain di erent ratios of tungsten and nickel, hence accommodating di erent amount of areas with ductile 23 fracture. Thus, the activation heatmaps in Fig. 13 can provide an interpretation of the failure of the algorithm in grouping the data from these two SEM images. (a) Activation heatmap corresponding to A102. (b) Activation heatmap corresponding to A103. Figure 13: Activation heatmaps that correspond to the two SEM images obtained after scanning the sample with the 90wt% tungsten composition. Finally, the confusion matrix in Fig. 7 shows that the data points that belong to the images cropped from the A102 SEM image are mainly clustered together with the data points obtained by the sample with 97wt% tungsten composition and secondarily with the data points from the 92wt% tungsten composition sample. Fig. 14 presents the activation heatmaps of the SEM images whose data points are grouped together by the clustering algorithm. The similarities between the activation heatmaps of the A102 SEM image and the maps from the other two samples are signi cant and certainly more pronounced than the similarities between the A102 and A103 SEM images, which are obtained from the same sample. Undoubtedly, a more systematic study of this clustering algorithm with di erent fracture datasets is required to prove or disprove the assumptions that are made in this interpretation of the algorithm's internal operation analysis. This paper encourages further investigation and the authors are convinced that a profound understanding 24 (a) A102 (b) 92wt% (c) 97wt% Figure 14: Activation heatmaps that correspond to samples with (a) 90wt% (b) 92wt% (c) 97wt% tungsten composition. of the functionality of the recently introduced Machine Learning algorithms can open up new possibilities with novel ecient characterization tools that optimize current methods in Fracture Mechanics. 25 Summary This paper evaluates the performance of an Unsupervised Machine Learning data pipeline in the clustering and classi cation of the SEM fracture images obtained from WHA samples with di erent tungsten compositions. The algorithm achieves to cluster the fracture images of this dataset according to their tungsten composition with accu- racy that reaches to 81.2%. Additionally, the data points that are erroneously clustered belong mainly to a certain SEM image and when they are removed from the dataset the accuracy of the algorithm raises to 91%. Subsequently, by adding a minimally supervised Machine Learning algorithm at the end of the clustering pipeline and de ning a new computational framework, the classi cation of the fracture images according the tungsten composition of the sample that they belong is enabled. The classi cation accuracy is computed to 93%, while the Haralick texture descriptors fail to classify this fracture images dataset, since they achieve accuracy of 21%. Finally, in an attempt to interpret the functionality of the clustering algorithm, this research work introduces a method that enables the visualization of the features in the input fracture images that are activated by the algorithm and enable the identi cation and clustering of the images according to the tungsten composition. More speci cally, a set of activation heatmaps - one for each input fracture image - is computed to enable the visualization of the exact locations in each input fracture image that are activated by the algorithm and consequently are used to enable the accurate clustering of the dataset. The interpretation of these activation heatmaps shows that the algorithm is using the ductile micro-fracture modes (dimples), which are the product of the ductile fracture of the Nickel-Ferrous matrix, in order to group the WHA dataset according to the tungsten composition. This estimation of ductility of the fracture surfaces is correlated to both the percentage of the NiFe presence in the alloy as well as to the relative amount of matrix failure exhibited by the specimen and this creates a criterion for the clustering process. One can speculate, that since di erent loading scenarios 26 (eg. loading rate, triaxiality etc.) will result in di erent crack paths and will a ect the relative amount of ductile vs brittle fracture indicators in the surface, a similar pipeline could be used to classify specimens following the conditions at which they failed. This hypothesis will be tested in a future work. Data availability The dataset used for this study and the corresponding Activation Maps are pub- lished in Materials Data Facility (MDF) with DOI: https://doi.org/10.18126/aph0-olbz . The source code is available at https://github.com/SteliosTsop/WHA_Clustering_ Classification_Visualization . Acknowledgments The nancial support provided by the Pazy foundation young researchers award Grant # 1176 and the European Union's Horizon2020 Programme (Excellent Science, Marie-Sklodowska-Curie Actions) under REA grant agreement 675602 (Project OUT- COME), is gratefully acknowledged. 27