Petrography is the branch of petrology focused on classifying and describing rocks from both a microscopical and a megascopical point of view [ 3 ]. With the difusion of Scanning Electron Microscopes (SEM), optical analyses for petrography reached good quality results due to the high resolution of the images. Indeed, SEM images can reveal features that are invisible with optical microscopy [ 8, 12, 15 ]. In particular, they allow distinguishing different mineral phases based on variations of gray intensity. Other information such as the size, distribution, shape, orientation and textural relationships of minerals can be retrieved.

Specifically, SEM images are highly useful for pore characterization in rocks [ 2, 16 ]. Petrographic Image Analysis (PIA) studies the presence and distribution of pores, which is important for analyzing the permeability of rocks to natural gas and oil [ 6, 9 ]. Within this process, SEM images are retrieved from ground samples called thin sections. After acquisition, image processing algorithms can extract diferent characteristics of pores, for example their size and shape. [ 4, 13 ]. Finally, domain experts use these results to inspect rock properties such as their permeability.

Since each thin section possibly contains from thousands to millions of pores, their manual categorization is applicable only to a small subset. This entails an approximated result that highly depends on the sampling performed over the initial data. An automatic categorization/clusterization method may help avoiding most of the manual efort that is required to perform the analysis. Both classification and clustering methods could be taken in consideration. However, classification requires having a suficient amount of training data with manually generated labels. This implies that domain experts should label large quantities of pores to model the characteristics of diferent categories. Since there are no public datasets with labeled pores and generating new labeled data is very time consuming, we adopted a clustering based approach. Clustering methods are able to learn the data distribution and automatically derive a set of groups, which then can be more easily and quickly mapped to geological categories by domain experts. Clustering algorithms do not need training labels and can be applied to thousands of samples. With the proposed methodology, domain experts can concentrate their eforts only on the geological interpretation of each cluster, while leaving to the algorithm the most time consuming task.

In this paper, we propose the Big Petro pipeline, a clustering methodology specifically designed for automatically categorizing pores inside thin sections. The proposed pipeline has been designed and evaluated on real datasets with the support of the domain experts of ENI, a leading company in oil and gas extraction. The contribution of our paper can be summarized as follows: • Development of a customized semi-automatic process for clustering geological pores. This methodology helps domain experts to avoid manually categorizing all pores. • Definition of a variant of the standard usage of hierarchical clustering. Specifically, we propose to obtain clusters with an Adaptive Multi-level Dendrogram Cut, being guided by both automatic metrics and domain experts knowledge. • Inspection of diferent methodologies to describe the generated clusters in terms of characterizing attributes distribution for each cluster. The description approaches allow domain experts to interpret the diferent groups and derive interesting insights from the analyzed thin sections. 2

Clustering techniques [ 11 ] are widely used in many application ifelds. To the best of our knowledge, this is the first attempt to apply clustering techniques to automatize pore analysis. During our preliminary study of geological pores categorization, we inspected the behavior of the following algorithms: DBSCAN [ 5 ], KMeans [ 7 ], and agglomerative hierarchical clustering [ 10 ].

DBSCAN [ 5 ] is a density based algorithm. Its driving idea is joining points in the same cluster when they exceed a minimum neighborhood density. The main issue with this technique is that the distribution of our data shows variable densities. Setting a low minimum density results in a single cluster containing all the pores, while when setting a high density many pores are considered as noise. Multiple runs of DBSCAN with diferent configurations could potentially allow identifying clusters with diferent densities. However, also the self-tuning techniques that use multiple runs of DBSCAN [ 1 ] did not yield good quality results on our data.

K-means [ 7 ] is a well-known centroid based clustering algorithm, which requires specifying the (fixed) desired number of clusters. The results with this technique were not satisfactory, as this method builds globular clusters with a balanced number of points, but the types of pores in our datasets present more complex shapes and are very imbalanced in number. For example the category of small pores is typically an order of magnitude more numerous than the one containing bigger pores.

Hierarchical clustering [ 10 ] techniques allow defining clusters by inspecting the cluster hierarchy with a structure called dendrogram. By cutting the dendrogram at a fixed level, k clusters are generated. The generated clusters can be characterized by diferent densities and imbalanced cardinalities. This technique, already in its standard version, allowed us to obtain the best results, according to the analysis performed by ENI domain experts. For this reason, we selected hierarchical clustering as preferred method to be customized and integrated in our Big Petro pipeline. 3

The objective of our work is to automatically discover diferent categories of geological pores in SEM images of rock samples by means of the proposed Big Petro pipeline. The Big Petro pipeline is composed of diferent blocks, depicted in Figure 2. In the following paragraphs we describe the characteristics of each block.

Data extraction. Each dataset to be analyzed consists of a single thin section acquisition, represented as a grayscale image. Figure 1a shows an example of a small portion of thin section, where the black irregular spots are the pores to be categorized. In this first step, a tool developed by ENI processes the grayscale image that automatically selects the darker regions representing pores. Afterwards, the tool analyzes each pore and extracts a set of geometrical features describing the pore shape and size. The result of this operation is a structured dataset, with a record for each pore, characterized by 32 numerical attributes. Section 5 describes more in detail the characteristics of the considered datasets, extracted from diferent thin section samples.

Data preparation. The second step of the pipeline involves the preprocessing of the datasets obtained in the data extraction phase. Each dataset is typically characterized by the presence of many small pores represented by only few pixels. The low pixel resolution of these pores does not allow computing most of the geometrical characteristics necessary for the analysis. Hence, we filter out the smaller pores by applying a domain provided threshold on their size. This operation also yielded smaller, more manageable datasets (from millions of pores to thousands), which allowed a reduction of the clustering algorithm processing time. In Section 5, Table 1 shows further details of the dataset cardinality before and after this filtering procedure.

Next, attribute values are normalized by applying the z-score normalization, after which each attribute domain range is characterized by mean equal to zero and unitary variance. This step is necessary to guide the distance metric used for clustering to treat all the attributes in the same way, independently of their value range.

Adaptive Multi-level Dendrogram Cut. Pore clustering takes place in this phase. It is composed of four steps: (i) dendrogram construction, (ii) super-clusters generation, (iii) domain experts analysis, and (iv) sub-clusters generation. The first task is addressed by running an agglomerative hierarchical clustering algorithm, which builds the dendrogram describing the hierarchy of the pore aggregations. The second step consists of generating a ifrst set of clusters by cutting the dendrogram at a specific height. These clusters correspond to macro-groups of pores, called superclusters and provide an initial coarse-grain aggregation of pores. Their number is chosen by analyzing the silhouette score [ 14 ], as described in Section 5.

Next, each super-cluster is manually analyzed by domain experts to decide if further partitioning is needed. The standard approach consisting of one single cut of the dendrogram may generate some clusters already containing homogeneous pores and some other clusters grouping heterogeneous pores. Simply cutting the dendrogram at a diferent height will split also the already high-quality macro-clusters. For this reason, only for the subset of super-clusters that are deemed by the domain experts to need further split in sub-clusters, the sub-cluster generation step is executed. Specifically, given a cluster to be further divided, we analyze its corresponding dendrogram sub-tree and cut the hierarchy at a new level of depth. The process continues with a loop involving the domain expert analysis and the sub-clusters generation. At each iteration, new levels of sub-clusters are added until the results are satisfactory. Further details on this process are provided in Section 4.

Cluster description. In this step, the obtained clusters are analyzed with diferent descriptive techniques. In Section 5, we show how, by inspecting the attribute distribution in the obtained clusters and the aggregation hierarchy in the dendrogram, interesting insights are obtained. We also consider PCA plots to show the distribution of the cluster labels along the directions where data have maximum variance. 4

The dendrogram is generated by an agglomerative hierarchical clustering algorithm. Since our datasets consist of continuous attributes, the euclidean metric is used for computing the distance between pore pairs. After inspecting the results obtained with diferent linkage techniques, we selected Ward’s method [ 17 ]. In particular the single linkage criterion did not fit our purposes, as it is sensible to noise and the points belonging to lower density clusters are not merged properly in the dendrogram. More specifically, single linkage connects two clusters when they present a pair of points which are very close together. For this reason higher density clusters are merged together in the first stages of the dendogram and points in lower density regions appear in different isolated clusters. Instead, complete linkage tends to break bigger clusters to have equal-sized groups. The algorithm behaviour with this metric is approximately similar to K-means and does not produce good quality results. Average linkage produces a trade-of result between the single and complete methods. However, by inspecting the first levels of the dendrogram together with ENI domain experts, we observed that Ward’s method yields a better separation of the main pore groups.

As mentioned before, agglomerative hierarchical clustering yields a binary tree structure, called dendrogram. The leaves of the tree represent single pores to be clustered, while nodes describe the merging points of clusters at diferent hierarchical levels. The height of each merging point is defined by the linkage distance between the two nodes to be merged.

Let Ch be the set of k clusters obtained by cutting the dendrogram at a specific height h:

Ch = {c0h , ..., cih , ..., ckh−1} The dendrogram cut at height h is exploited to generate the first level of groups, denoted as super-clusters. For example in Figure 3 the obtained super-clusters are tiny-pores, small-pores, big-pores1 and big-pores2 (these names have been conventionally assigned by domain experts after the analysis of the macro-clusters). The value of h is selected by inspecting the silhouette value of the obtained groups. In Section 5 we further discuss how this value is chosen for the diferent analyzed datasets.

When the domain experts deem the granularity of some macroclusters as insuficient, the iterative generation of sub-clusters is activated. Sub-clusters are generated by inspecting the descendants of their corresponding super-cluster. The set of sub-clusters of a cluster cih obtained by cutting the dendrogram at height h′ is defined as:

subclustersh′ (cih ) = {chj′ : chj′ ∈ descendants(cih )} where descendants(cih ) is the set of descendants of cih in the tree hierarchy and chj′ are the nodes obtained by cutting the dendogram at the new height h′. For example, in Figure 3 the set subclustersh′ (tiny- pores) is composed of sub-clusters tiny-a, tiny-b, tiny-c. The definition of subclustersh′ (cih ) can be applied recursively to obtain multiple sets of subclusters, at diferent levels of the hierarchy for diferent cluster subsets, until the results are satisfactory (see Figure 2). 5

ENI provided us three diferent datasets ( Dset1, Dset2 and Dset3), whose samples were extracted from three thin sections mostly composed by carbonates. As described in Section 3, the smallest pores in each dataset are filtered out using a threshold specified by ENI domain experts. Consider the left part of Table 1. The biggest dataset, Dset3, presents 1.3M pores, while after filtering its size is 4K rows. The smallest dataset, Dset2, changes from 300K to 5.9K pores. Note that, even if the original size of the datasets can be quite diferent, the number of pores after filtering takes similar values, ranging between 4K and 6K . 1

Super-clusters generation. After building the dendrogram on the three datasets, we selected the proper values of h to cut the hierarchy tree and obtain the super-clusters. In particular, Figure 4 plots the silhouette values against the number k of super-clusters in each dataset. For all datasets the silhouette score decreases when k grows. However, even if the silhouette is higher with lower values of k, a very low number of superclusters does not provide meaningful groups from a geological point of view. Datasets Dset1 and Dset3 show stable silhouette values in the range k = 3 to k = 4. Hence, we selected k = 4, the maximum value of k before the silhouette rapidly decreases, for these two datasets. Instead, dataset Dset2 shows stable silhouette values for k = 2, 3 and k = 4, 5. For the same motivation, we set k = 5, i.e., the highest number of clusters in the region of stable silhouette values.

Table 1 shows the size of the obtained super-clusters. All datasets are characterized by two big clusters (tiny-pores, smallpores) and two smaller ones (big-pores1, big-pores2). Moreover, Dset2 shows a cluster (big-pores3) with a single pore, which can be considered an outlier. This confirms that the choice k = 5 for dataset Dset 2 produces similar results to the ones obtained with the other two datasets.

Figures 6a, 6b, and 6c show the results of a Principal Component Anaysis (PCA) on the pores of the three datasets by visualizing the projection of the points on the first three components. The diferent colors represent points belonging to the generated superclusters. For all three datasets, the two biggest clusters (tiny-pores, small-pores) are characterized by the highest density, while the others are sparser. Furthermore, the cluster structure is similar for all datasets. The biggest clusters are characterized by lower values of component pc1, while the smaller ones by higher values of pc1. Finally, the outlier in Dset2 is clearly visible as an isolated point in Figure 6b (big-pores3).

Sub-clusters generation. According to the analysis performed by ENI domain experts, clusters big-pores1 and big-pores2 are pure and well characterized for all datasets. Instead, clusters tiny-pores and small-pores need to be further separated. Being guided by ENI domain experts we generated k ′ = k ′′ = 3 sub-clusters from tiny-pores and small-pores respectively. Figure 3 depicts the resulting dendrogram for dataset Dset1. The names of the subclusters are obtained by appending a sufix ( -a,-b,-c) to the name of their corresponding super-cluster (e.g., tiny-a, tiny-b). Note that the cut point for the sub-clusters of small-pores is deeper in the dendogram than the one for the three sub-clusters of tinypores. This result cannot be obtained with a single dendrogram cut and motivates the need for the proposed Adaptive Multi-level Dendrogram Cut. To reach the granularity of the small-pores subclusters with a single dendrogram cut we need to choose a height H corresponding to K = 22 clusters. However, this K value would break the clusters big-pores1, big-pores2, tiny-a, tiny-b, tiny-c and would produce a lower quality clustering.

Figure 5 shows the silhouette value for the final clusters obtained with Adaptive Multi-level Dendrogram Cut and those generated with a single dendrogram cut. The number of clusters for the single dendrogram cut is K = 22, 19, 13 for Dset1, Dset2 and Dset3 respectively. The silhouette of the clusters obtained with our method is clearly higher than the one based on one single cut. This entails that our method not only produces a better partition from a geological point of view, but also a better domainagnostic quality. A further discussion on the results shown in Figure 5 is provided in Section 6.

Cluster description. We provide a description of the obtained clusters by inspecting the attributes that allow characterizing them. The PCA representation in Figure 6 shows the distribution of the diferent super and sub-clusters. We inspected the attributes which are most representative of each of the three principal components, by considering the ones with highest eigenvalues. For example, the first component, pc1, is characterized by size attributes. In particular, higher values of pc1 entail higher size of the pores. In Figures 6a, 6b, 6c, the super-clusters tinypores and small-pores are characterized by lower values of pc1, while big-pores1 and big-pores2 present higher values. Figure 7a shows the distribution of one of the size attributes in Dset1 for the diferent super-clusters. It confirms that the distinction between super clusters is partially driven by the size of the pores. (a) Super-clusters of Dset1. (b) Super-clusters of Dset2. (d) Final clustering of Dset1. (e) Final clustering of Dset2. (f) Final clustering of Dset3.

Figures 6d, 6e, 6f show the final clustering, including both the sub-clusters and the super-clusters that have not been further divided. It can be noticed that sub-clusters present a finer separation that also involves components pc2 and pc3. For example in Dset1 clusters tiny-b and tiny-c are well separated by pc2. The most relevant attributes for pc2 describe how much pores are stretched, obtaining a non circular shape. Looking at Figure 7b it is possible to see that pores in tiny-c show a higher stretching than those in tiny-b that are more circular. The third principal component, pc3, is related to pore irregularity, instead. From Figure 6d pores in cluster small-a present a lower shape irregularity than those in small-c. This is also confirmed by the histograms in Figure 7c. Following this methodology, ENI domain experts are able to derive interesting insights, describing the geological characteristics of the analyzed thin section. 6

The analysis conducted in the previous sections highlights the following results. First, the silhouette score, a domain-agnostic quality index for evaluating clusters, is not always a good metric from a domain driven point of view. Consider Figure 5. The silhouette score computed on the super-clusters is higher than the one obtained with the Adaptive Multi-level Dendrogram Cut. For example, in Dset3 super-clusters reach a silhouette score of 0.73, while Adaptive Multi-level Dendrogram Cut clusters score 0.2. However, the quality of the clusters from the domain experts point of view is higher for the proposed method, as the super-clusters show an excessively coarse subdivision. In fact, super-clusters have a strong correlation with the pore size, but cannot capture more complex characteristics of the pore shape.

Second, the cluster description confirmed that sub-clusters represent a finer subdivision of the initial groups according to more complex characteristics, such as the pore stretching and shape irregularity. According to the domain experts, this grouping has a good quality from a geological point of view, because it generates groups that are distinguished by both geometrical and genetic characteristics. These characteristics are fundamental for geologists to inspect the structure of the analyzed thin section and its permeability. Hence, the role of domain knowledge is fundamental in driving the clustering process to obtain higher quality results. 7

In this paper, we proposed the Adaptive Multi-level Dendrogram Cut method for clustering geological pores in thin sections and we evaluated its results on three diferent real datasets provided by ENI. We demonstrated that the Adaptive Multi-level Dendrogram Cut allows obtaining sub-clusters that capture more complex characteristics of the pore shape and have higher silhouette values than the sub-clusters that would be obtained with the corresponding single level dendrogram cut. Currently, the Big Petro pipeline is under further evaluation by ENI domain experts, which are planning to deploy it shortly.

As future work, we plan to integrate this pipeline with domain driven self-tuning techniques that will allow fully-automatic pore clustering. Furthermore, we will work on improving the cluster description step to help domain experts in interpreting clustering results.

The research leading to these results has been supported by the SmartData@PoliTO center for Big Data and Machine Learning technologies.