Carbonate texture identification using multi-layer perceptron neural network

2.1 Image acquisition

Image acquisition is the first step in an image processing system. The general aim of any image acquisition was to transform an optical image into an array of numerical data which can manipulate on a computer. Suitable digital cameras achieved image acquisition. A digital image is a two-dimensional array of pixels, and each pixel characterizes by its (x, y) coordinates and its value. Digital images are characterized by matrix size, pixel depth, and resolution, containing a gray tone compared with real images.

The matrix size is determined from the number of columns (m) and the number of rows (n) of the image matrix (m × n). The acquired images are in jpg images format with 8-bit, 256 linear gray tones, and the size of a matrix selected by us is 680 × 700 pixels, acquired by digital cameras with different dots per inch. The utilized magnification was about ×4–10, which is acceptable for acquiring smaller and larger grains. Since the color variability of the studied thin section isn’t distinctive, all images are converted to grayscale for easy calculations and to get a homogeneous database. The proposed methodology counts only the shape of white and gray areas, considered cement, matrix (micrite), and grains, respectively.

2.2 Image preprocessing

Image preprocessing is done to improve the image’s quality and upgrade the image highlights for further handling. In the gained images, the following image preprocessing techniques in Python are as follows:

1. grayscale conversion,

2. normalization,

3. data augmentation, and

4. image standardization.

To use this preprocessing technique, we have created a Google Colab account and Google Drive data onto Colab. The next step was importing the data set to Colab. We loaded an image file from the data set and randomly printed it (Figure 1).

Figure 1

Flowchart of the proposed methodology.

2.2.1 Grayscale conversion

Grayscale is converting images from colored to black and white. Since most pictures do not need color to be recognized, we use grayscale, which reduces the number of pixels in an image, thus, reducing the computations required in machine learning algorithms [27]. From executing the code, the image was converted as below (Figure 2).

Figure 2

(a) Printed image of mudstone from the data set; (b) conversion of the original image into grayscale. In both figures, units of X and Y axes represent pixels.

2.2.2 Normalization

This process consists of projecting image data pixels (intensity) to a predefined range (usually (0, 1) or (−1, 1)). It is commonly used on different data formats, and we want to normalize all of them to apply the same algorithms to them. Its benefits include fairness across all images and a standard learning rate (Figure 3).

Figure 3

Output image after normalizing our data, units of X and Y axes represent pixels.

2.3 Image segmentation

Image segmentation helps to get the region of interest (ROI) from the image. This is the process of dividing an image into separate areas. The part where the image is divided is called an image object. This is based on image properties such as similarity and discontinuity. The purpose of image segmentation is to simplify the image for better analysis.

This is the process of assigning a label to each pixel in the image. There are several segmentation algorithms like- Otsu’s segmentation, edge-based segmentation algorithms, region-based segmentation algorithms, clustering-based segmentation algorithms etc. We find it most appropriate to use Otsu’s Image segmentation (threshold-based segmentation) since some of them can take more time during training if the data set is huge and take huge computation time (Figure 4).

Figure 4

Input image segmentation according to Otsu’s segmentation: (a) input image; (b) grayscale converted image; (c) threshold output image. Units of X and Y axes represent pixels.

2.4 Feature extraction

Features are the primary attributes that help us identify a particular image and are the properties that make it unique. While working on the image data set, you need to extract various image features. This helps to separate images based on specific features or aspects. Image features cannot be defined exactly, but shape, size, and orientation make up image features. This step represents a sequence of operations on the image matrix to measure the features of the objects. These values are given as input layers (i.e., pixels) to the neural network. Marmo et al. and Gonzales have performed various test to select different features in carbonates and related images [9, 28]. According to the authors, the final set of determined features is listed below:

1. percentage of the most frequent gray tones;

2. number of edge pixels in the gray areas;

3. number of edge pixels between the gray and the white areas;

4. number of white areas;

5. number of pixels of white areas;

6. number of white areas composed of more than 150 pixels;

7. number of pixels within the white areas consisting of more than 150 pixels.

Figure 5

(a) Foraminiferal grainstone; (b) image histogram of a grainstone texture shown in (a). Units of X and Y axes represent pixels.

Figure 6

(a) Canny edge detection output in a mudstone image; (b) Canny edge detection output in a grainstone image. Units of X and Y axes represent pixels.

Figure 7

(a) Mudstone composed of a few skeletal grains (ostracods) in a micritic matrix; (b) image histogram of a mudstone texture shown in (a).

Figure 8

(a) and (b) Perimeter pixels in a grainstone.

Features (1)–(7) were extracted from the entire image while the remaining 16 features were extracted from each of the four selected largest white areas (4 features (8)–(11) per area):

1. number of pixels of the white area;

2. length in pixel of the white area;

3. number of small convex deficiencies;

4. number of large convex deficiencies.

They represent 23 features, and each of them is described further in this article.

1. Percentage of the most frequent gray tones

   In an image processing context, an image’s histogram refers typically to a histogram of the pixel intensity values. This histogram shows the number of pixels in an image at each different intensity value found in that image. If we use an 8-bit gray scale image, there are 256 different possible intensities, so the histogram will graphically display 256 numbers showing the distribution of pixels amongst those grayscale values. The exact output from the implementation operation is an image of the required histogram. It is composed of equally spaced bins, each representing a range of gray tones. It is necessary to sum the number of pixels of the bins chosen around the bin that corresponds to the most frequent gray tone. The value of this feature is the ratio between this sum and the sum of pixels in all bins. Since the image of a grainstone consists of many shades of gray, the histogram is distributed almost across all possible grayscale values (Figure 5). On the other side, an image of a mudstone shows a narrower and higher histogram because a very fine-grained mudstone is represented by many pixels distributed on a few shades of grayscale values (Figure 6).

2. Number of edge pixels in the gray areas

   In Image Processing, an edge is defined as continuous pixel positions where an abrupt change of intensity (gray or color) values occurs. Edges represent boundaries between objects and backgrounds. We used the Canny Edge detector algorithm to identify edges in the whole image. The detector looks for places where the intensity changes rapidly. The resulting image is a black/white image in which a black pixel corresponds to an edge pixel in the original image. We choose this edge detector because it is more likely to detect weak edges. Canny recommended an upper and lower ratio of 2:1. The low and high thresholds are 200 and 100, respectively. The value of the feature is the number of black pixels. The resulting output image of mudstone is a few black pixels because it is composed of very fine-grained particles (<0.05 mm) as the micritic matrix (Figure 7). At the same time, a grainstone has many black pixels because it represents a grain-supported texture. It is composed of numerous grains of skeletal and non-skeletal origin that generate many edges (Figure 8).

3. Number of edge pixels between the gray and the white areas

   A pixel is part of the perimeter if it is nonzero, and it is connected to at least one zero-valued pixel. The pixel connectivity used is an 8-connected neighborhood. In other words, a pixel is considered as a perimeter pixel if it satisfies both criteria, a white pixel in the starting image that shows one or more gray pixels in its 8 × 8 neighborhood. The value of the feature is the number of perimeter pixels. A grainstone has many perimeter pixels due to the high content of grains due to many white areas, while a mudstone shows few white areas and is characterized by few perimeter pixels (Figure 8).

4. Number of white areas

   Grayscale images are converted to black and white, and gray pixels in the original image are converted to black. The algorithm finds and counts white areas. The feature value is the number of white areas, and grainstone has far more than mudstone.

5. Number of pixels of white areas

   The value of this feature is the number of pixels that make up the white area used to extract the previous feature.

6. Number of white areas composed of more than 150 pixels

   The previous black-and-white image and the image obtained from the inscription are taken into account. The number of pixels belonging to each white area is estimated, and the feature value is the number of white areas exceeding 150 pixels (Figure 8). This threshold was experimentally evaluated on thin-section without significant post-depositional features (dolomite crystals, fractures, and veins filled with sparry calcite).

7. Number of pixels within the white areas composed of more than 150 pixels

   The value of this feature is the number of pixels in the white area that has already been used to calculate the previous feature. Mudstone and wackestone do not show large white areas.

8. Number of pixels in the white areas

   The results of the previous labeling algorithm are considered here. The white areas are ordered according to their number of pixels. The four most significant white areas are chosen by estimating the number of pixels for each. In this way, four features are obtained. Different types of grainstones (foraminiferal grainstones, peloidal grainstones, etc.) have been used to count the number of white areas to determine the above threshold.

9. Length in pixel of the white area

   The feature value for each of the previous four most significant white areas is the number of columns used to represent the area in the image matrix.

10. Number of small and large convex deficiencies

    The convex hull and the corresponding convex deficiencies are computed for each of the previous four white areas [29]. A labeling algorithm locates and counts the convex deficiencies. The value of the features (j) and (k) corresponds to convex deficiencies composed of less than 100 pixels and convex deficiencies consisting of more than 1,000 pixels, respectively. The values 100 and 1,000 were determined by analysis of a wide variety of white areas in images of grainstone and packstone. Due to intergranular calcite cement, a white area has a large convex hull, a low number of minor convex deficiencies, and a high number of large convex deficiencies (Figure 9). On the contrary, a white area due to intragranular calcite cement has a small convex hull, a high number of small convex deficiencies, and a low number of large convex deficiencies (Figure 9).

Figure 9

(a) Determination of convex deficiencies; (b) large convex deficiencies.

3.1 Neural network model and architecture

Our study restricts the number of output classes to four (mudstone, wackestone, packstone, and grainstone). Due to this, we found the most appropriate to use the multi-layer perceptron neural network (DLP). A multi-layer perceptron (MLP) represents a feedforward artificial neural network. It is composed of an input layer to receive the signal. This output layer decides or prediction about the input. Between those two, an optional number of hidden layers are the true computational engine of the MLP. In our case, we have used just one hidden layer. Each neuron relates to the neurons of the following layer (Figure 10). According to PCA (principal component analysis), which is used to reduce the number of parameters in the neural network, the input layer has 13 neurons; the hidden layer has nine neurons, while the output layer has four neurons. Each one of them corresponds to a Dunham textural class. Other parameter configurations show less than 91.0% of correct classification.

Figure 10

The architecture of the proposed MLP where each layer consists of 13, 9, and 4 nodes.

3.2 Activation functions for each neuron

The activation function of a neural network defines how the weighted sum of inputs is transformed into outputs from one or more nodes in a layer of the network. The activation function used in the hidden layer is usually chosen based on the nature of the neural network architecture. We used the Rectified Linear Unit (ReLU) activation function, the most common function for hidden layers for neural network models with MLP architectures [32]. The activation function for the output layer is based on the type of prediction problem that we are solving. Since our problem is a classification problem and more than two mutually exclusive classes (multiclass classification), our output layer will have one node per class, and a softmax activation should be used [32].

3.3 Loss function

Cross-entropy is the default loss function that is used for multi-class classification problems. Cross-entropy will calculate a score that summarizes the average difference between the actual and predicted probability distributions for all classes in the problem.

3.4 Training and validation of the data set

Once configured, the neural network model needs to be trained on our data set consisting of input features and the desired class as output. The back-propagation method used the training procedure: we give the input neurons the primary pattern, and therefore, the network provides its output. If this is not equal to the desired output pattern, then the procedure computes the difference (mean square error) between these two values and changes the weights to minimize them. These operations are repeated for each. The data set was divided into training, validation, and test data set. We assigned small random values except 0 to weights utilizing 500 gray-scale images of thin sections from the data set to train the neural network, 125 for each Dunham class. The image data sets represent micro photos from previous studies carried out by different authors in the shallow-water carbonates of the Upper Cretaceous, Kruja zone (Albania) [33,34]. Thin sections affected by dolomitization and recrystallization were not taken into count by us. Our model, after each epoch, learns the data and updates the weights properly. After some iterations, the model starts memorizing the data, and even though the training error decreases, the validation error increases, causing the overfitting of the model. We used early stopping as a form of regularization to avoid overfitting (Figure 11). The validation data set is represented by 154 images (38 images for each texture class), and it was tested every 20 epochs. During the training and validation tests, the training error in continuity decreases as training proceeds, while the validation set error only decreases up to a certain point (two points). It starts to increase (Figure 11). As shown in Figure 11, after some iteration, the training is stopped early to prevent model overfitting at the minimum of the validation error that corresponds to 0.63. The related training error is 0.34.

Figure 11

Training and validation error every 20 epochs.

3.5 Test of the data set

The test set is a separate data set used to test the trained model. We utilized 348 images (87 for each class) taken from the same data set to test the technique. The MLP has classified exactly 315 (90.5%) images, and 33 (9.5%) images were incorrectly classified. Regarding mudstone class type, two images were classified as wackestone and packstone due to the secondary porosity of rounded shape filled with sparry calcite (Table 1). As for the wackestone class type, three images were classified as mudstone due to the very fine-grained bioclasts that cannot be distinguished at the used magnification. One of the wackestones was classified as packstone because of irregularly distributed empty fenestrae (Table 1). If we look at the packstone class type, we can notice that one of the images was classified as wackestone. Very fine peloidal grains are not detectable at the used magnification.

On the other hand, six packstones images were classified as grainstones due to the coarse bioclast (mainly rudist fragments) occupying the entire field of view (Table 1). The last class type corresponds to grainstones, and it is the most precise classification since only three of the images were incorrectly classified as packstones. This is because of the presence of intraclasts that are mainly reworked packstones. To summarize the visualization of the algorithm’s performance, we have used a confusion matrix. It is a table where the columns represent the true classes, and the rows represent the four predicted classes (Table 1). From this table, every error that occurs with the adjacent with the correct one, except mudstone class where only one of the images is misclassified as packstone, but we can see that a grainstone is never classified as mudstone and vice-versa (Table 1 and Figure 12).

Figure 12

Misclassified image in test set, green points show correct class and red points corresponding neural network output.

To fulfill Table 1, we used a classification scheme where the highest output suggests the choice of the class. The problem is when the MLP’s four output values are very similar, and the highest value differs from 1. In such cases, the MLP classification of the image has a high probability of being uncertain. Only with the specialist’s help is it possible to classify the aforementioned images precisely. Those were chosen using a simple threshold rule. The specialist should classify the relevant image if the highest output value is less than or equal to the rejection threshold. Otherwise, the proposed MLP classification can be accepted (Figure 13).

Figure 13

The classification percentages for each value of rejection threshold in range [0.0, 1.0].

The optimal rejection threshold value is significant. If the rejection threshold value increases, more images are rejected, making fewer errors.

Each of the highest output values is compared to a rejection value in the range [0.0, 1.0] using the MLP output values for the test set. According to Figure 13, the number 0.55 appears to be the rejection threshold for classifying an appropriate percentage of images as correct, wrong, or uncertain. The threshold distinguishes 4% of the images as uncertain (14 images), 13.0% of the images as wrong (45 images), and 83% as correctly classified images (289 images). Calculating the difference between the first and second highest values is another criterion used to distinguish between uncertain classes. By additional testing, the MLP is unsure between the two classes corresponding to these two values if the difference is smaller than 0.4. Only 13 images need to be classified by direct observation under the microscope.

3.6 Geological setting

In the Triassic period, a general break-up of Adria, which was followed by the Early Jurassic disintegration of the Southern Tethyan Mega platform and the formation of the Adriatic Basin, resulted in the creation of the Apulian Carbonate Platform, Apennine Carbonate Platform, and Central Mediterranean Carbonate Platform [35]. During the Mesozoic and Paleogene, the passive margin of the Apulian Plate was bordered by broad carbonate platforms, known in Albania as the Kruja Zone or Kruja Platform. The Kruja Zone represents the most western coastal zone in the Outer Albanides. It lies in the center of Albania, and it extends northward into the Adriatic Dinaric platform (Montenegro, Croatia) and southward into the Gavrovo platform (Northern Greece) [36]. In Albania, only upper Mesozoic Paleogene deposits are exposed [37,38,39,40].

Sedimentologic and petrographic studies have been carried out on an Upper Cretaceous shallow-water carbonate succession near Burizana. This section lies northwest of the anticline structure, part of the Kruja Zone (Makareshi carbonate anticline structure) that represents the Apulian eastern passive margin [41] (Figure 14). The stratigraphic succession of Kruja is mainly composed of thick deposits of Cretaceous carbonates (about 1.5 km thick), confined in a subtidal to the supratidal environment.

Figure 14

Map facies from middle-upper Lias to upper Cretaceous in the peri Adriatic zone (after Zapaterra, 1994).

This succession is interrupted by a stratigraphic gap that begins during the late Maastrichtian and continues (depending on the zone) until the Early Eocene [42]. During this gap, bauxite formation was formed [43]. After this gap, the series continues with bioclastic limestones (Middle Eocene). Above the limestone formation, the “passing marl” horizon of Late Eocene marks the transition to the Oligocene flysch, 1500 m thick. The Upper Cretaceous carbonate deposits of the Burizana section display eight typical microfacies, ranging from subtidal to supratidal environments [44] (Figures 15 and 16).

Figure 15

Lithological column of the Burizana section (Qorri, [13]).

Figure 16

Images from different microfacies of the Burizana section: (a) mudstone with sparse small benthic foraminifera (mainly miliolids) and problematic algae Thaumatoporella; (b) wackestone with horizontally oriented ostracods; (c) packstone with large bivalve sections (rudists); (d) grainstone with benthic foraminifera (Moncharmontia sp., Quinqueloculina sp.).

3.7 Burizana section as a key study

The test set of the Burizana section (Makareshi structure, Albania) consists of a total of 250 thin sections, which include 75 mudstones (30%), 70 wackestones (28%), 46 packstones (18.4%), and 59 grainstones (23.6%). The above four Dunham textural classes contain transition textures like mudstone-wackestone, wackestone-packstone, and packstone-grainstone. There are approximately 25 thin sections that were excluded from the test set because the original texture was affected by meteoric diagenesis. Since the thin sections of the Burizana section were previously studied under the microscope, it was possible to correlate manual and automatic classification.

The MLP has correctly determined the texture with 91.3% accuracy, which means that 229 of the entire images are correctly classified. The proper confusion matrix is displayed in Table 2.

Sedimentologic and petrographic characteristics can complicate the primary depositional texture, which is not captured by our technique. Neomorphism and partial dolomitization have been identified in the micritic matrix of wackestones. These diagenetic processes assign higher values to the packstone texture, confusing classification. For example, a small number of well-defined dolomite crystals in mudstone caused by diagenetic dolomitization were identified as grains during feature extraction, increasing the percentage of grainy material artificially.

On the other side, even the grain-supported texture shows some complications. Micritization, fine, and densely packed grains in some packstones are other complex features that our neural network cannot determine. As a result, the fine-grained packstones appear to have mud-supported texture and can be considered wackestone. In addition, since the studied section consists of shallow-water carbonates, the presence of intraclasts and peloids is very common, and it can be considered mud. As a result, the contrast between the grains (peloids and micritized grains) and the matrix is not always visible. Increasing the contrast may not personalize the different grayscale values during feature extraction. The grains should have a well-defined internal structure and clear boundaries so that the grayscale values contrast with the values in the matrix. These conditions are necessary for all packstones to be accurately determined by the neural network.

If we compare the classification accuracy between two textures, the mud-supported texture has been more accurately determined (94.3% accuracy) than the grain-supported texture (85.1 precision), and only 21 images were error classified. It is important to remember that the accuracy of determining texture depends on the difference between the two highest values. If the difference is less than 0.2, generally, the technique determines a mixture of two textures.