Identification of Salt Deposits on Seismic Images Using Deep Learning Method for Semantic Segmentation

: Several areas of Earth that are rich in oil and natural gas also have huge deposits of salt below the surface. Because of this connection, knowing precise locations of large salt deposits is extremely important to companies involved in oil and gas exploration. To locate salt bodies, professional seismic imaging is needed. These images are analyzed by human experts which leads to very subjective and highly variable renderings. To motivate automation and increase the accuracy of this process, TGS ‐ NOPEC Geophysical Company (TGS) has sponsored a Kaggle competition that was held in the second half of 2018. The competition was very popular, gathering 3221 individuals and teams. Data for the competition included a training set of 4000 seismic image patches and corresponding segmentation masks. The test set contained 18,000 seismic image patches used for evaluation (all images are 101 x 101 pixels). Depth information of the sample location was also provided for every seismic image patch. The method presented in this paper is based on the author’s participation and it relies on training a deep convolutional neural network (CNN) for semantic segmentation. The architecture of the proposed network is inspired by the U ‐ Net model in combination with ResNet and DenseNet architectures. To better comprehend the properties of the proposed architecture, a series of experiments were conducted applying standardized approaches using the same training framework. The results showed that the proposed architecture is comparable and, in most cases, better than these segmentation models.


Introduction
Seismic imaging enables the visualization of underground structures and it is used in the discovery of hydrocarbon fuel reserves. It is based on the emitting of sound waves that reflect on underground structures and are detected on a surface using receiver devices called geophones (see Figure 1). The reflected sound signal is recorded for further processing that leads to a threedimensional (3D) representation of underground rocks structure [1]. Seismic images show the boundaries of different types of rocks. Theoretically, the strength of the rejected signal is directly proportional to the difference in the physical properties of the rocks at the point of contact. This practically means that seismic images contain information about boundaries between rocky deposits, while they expose very little about the rocks themselves [2]. Seismic images are used in the exploration of hydrocarbon fuel reserves by helping detect potential reservoir rocks and that is why identification of salt deposits plays an important role. The development of salt domes (see Figure 2) can deform surrounding rocks forming traps that hold oil and natural gas. Salt sediments have characteristics that make them both simple and difficult for identification. Salt density is usually about 2.14 g/cm3, which is less than most nearby rocks. Salt seismic velocity of around 4.5 km/s is usually bigger than the surrounding rocks. This difference causes a sharp reflection at the salt sediment boundary. Salt sediment is usually an amorphous form without some special internal structure, which means that there is typically not much reflection within the sediment itself unless it comes from some other rock that is trapped within [2].  The goal of the TGS-NOPEC Geophysical Company (TGS) sponsored Kaggle competition [2] was to develop an algorithm that can classify each pixel in a 101 x 101 seismic image patch as salt or not salt. The data that was provided for training consists of 4000 seismic image patches along with appropriate binary masks that depict salt regions. For testing and scoring purposes, an additional 18,000 seismic image patches were provided. In addition, for each seismic image patch, depth information (in feet) is provided. Illustration of the dataset is depicted in Figure 3, where several pairs of seismic image patches and corresponding masks, along with depth information, are shown.
The way it is stated puts the problem into the semantic segmentation category. Semantic segmentation is one of the key problems in computer vision and represents a natural next step beyond image classification and object localization. Deep convolutional neural networks (CNN) [3] have been successfully applied to all these three problems. CNNs combine three architectural ideas that make them, to a certain degree, invariant to translation and distortion: local receptive fields, weights sharing, and spatial subsampling [4]. Deep CNNs are capable of building a hierarchy of features that makes them suitable for classification tasks, so as for the localization and semantic segmentation of objects in images. Although CNNs were known since the 1990s, deep CNNs are coming into the focus of interest in 2012 when a network named AlexNet [5] won the ImageNet Large Scale Visual Recognition Challenge [6] (hereinafter referred to as ImageNet). AlexNet achieved a top-five error of 15.3% which was 10.8% better than the second-placed solution. This result was made possible by training using graphic processing units (GPU) which is considered a turning point for the progress of deep learning. In the following years, there has been a significant improvement in deep CNN architectures that lead to much better classification results in ImageNet challenge. In 2013, the winner was ZFNet [7] with a top-five error of 14.8%. The main contribution of the authors of this network is not that much in the classification result but in the development of visualization technique by mapping learned filters to image patches. The following year brought a huge advance and two significant solutions. VGGNet [8] got a top-five error of 7.3% while promoting simplicity and depth. The winner in 2014 was GoogLeNet [9] with a top-five of 6.7%. It is interesting to mention that GoogLeNet is not just much more accurate, but also has 12 times fewer parameters than AlexNet. In the next year, the error dropped to 3.6% which is almost twice as good result compared to the previous year. It was thanks to ResNet [10] architecture that is inspired by VGGNet simplicity with the introduction of residuals that enabled efficient training of deep CNNs. The 2015 winner was ResNet variation with 152 layers. CNN architectures for classification are especially important because they can be relatively easily extended and used for semantic segmentation. There are many different architectures [11][12][13][14][15][16][17][18] that arise in the past couple of years to tackle this problem. In this paper, a novel architecture is proposed that came as a result of the author's participation in the TGS competition [2]. The architecture represents the evolution of the U-Net [14] model integrating some ideas from ResNet [10] and DenseNet [19] architectures. Since the introduction of this architecture provided significant improvement of the competition score compared to vanilla U-Net, a series of post-competition experiments were conducted employing standardized segmentation models with the existing training framework. The results showed that the proposed architecture is comparable and, in most cases, better than tested segmentation models.
The paper is organized as follows: Section 2 presents an overview of related work concerning semantic segmentation and salt deposits identification. In Section 3, the proposed method is described presenting the used CNN architecture and some implementation details. Obtained results are presented and discussed in Section 4. Finally, Section 5 presents the conclusions.

Related Work
The problem of seismic image analyses and salt identification attracts many researchers. Typically for all computer vision problems, the traditional approach to seismic image analyses was based on hand-crafting different feature extractors and later processing of appropriate responses. In one of the first papers in this area, Pitas and Kotropoulos [20] suggested a method based on texture analyses for semantic segmentation of seismic images. Methods based on texture attributes remain feasible in recent years [21,22]. Harper and Clapp [23] suggested ways in which different seismic image attributes can be calculated and used to identify salt deposits. A method based on the calculation of new seismic image attributes is also used by Shafiq et al. [24]. Analysis of 3D seismic images was the research subject in [25][26][27]. Amin and Deriche [25] propose a method based on using a 3D multi-directional edge detector. Wu [26] relies on probability calculation to identify salt sediment boundaries. Di et al. [27] propose multi-attribute clustering using the k-means algorithm to identify salt sediment boundaries. The use of machine learning techniques to identify four characteristic seismic structures based on different seismic attributes extracted from the image was suggested by Wrona et al. [28].
Deep CNN architectures for classification, that arise in recent years, can be relatively easily extended and used for semantic segmentation. The fully convolutional network (FCN) [11] represents one of the first solutions following that path. An FCN is constructed by removing the top classifier layers of CNN, adding a convolutional layer to reduce the number of filters to the number of output classes, and upsampling the obtained feature map to match the input size of the network. Since the spatial component of the last CNN feature map is several times smaller than input size, to achieve better results it is possible to form a hypercolumn [12] using feature maps obtained from several stages of the CNN. The number of filters in each feature map is equalized with additional convolutional layers and, after upsampling to a certain size, an addition operation is used to construct a hypercolumn that is used for pixel-level classification.
More complex CNN architectures for semantic segmentation are DeconvNet [13] and U-Net [14]. Both networks share a similar architectural idea: the network consists of encoder and decoder parts. In the case of DeconvNet, they are called convolution and deconvolution networks respectively. The encoder part is basically a classification CNN with a removed top, while the decoder has a mirrored structure with blocks that consist of upsampling or transpose convolution layers followed by convolutional layers. U-Net architecture, in addition, introduces skip connections that enable copying and concatenating earlier encoder feature maps to corresponding upsampled decoder feature maps, so decoder convolution layers process them both. A very similar approach, with slight differences, is promoted by SegNet [15] architecture.
Deep learning models for semantic segmentation become increasingly popular in the automatic processing of remote sensing images for building footprint extraction [29][30][31], road extraction [32], land use detection [33], etc. Motivated by good results in many areas, CNNs are becoming researchers' default choice for the segmentation of seismic images and identification of salt deposits. A huge number of papers [34][35][36][37][38][39][40][41] in 2018 and 2019 supports the claim. Dramsch and Lüthje [34] evaluated several classification deep CNNs with transfer learning to identify nine different seismic textures from 65 x 65 pixel patches. Di et al. [35] addressed the same problem using a deconvolutional neural network with three convolutional and three deconvolutional layers to speedup seismic interpretation. The same authors [36] addressed a problem of salt body delineation using a classification approach. They relied on the CNN with two convolutional and two fully connected layers, that was fed with 32 x 32 image patches. Waldeland et al. [37] proposed a method for the identification of salt bodies from 3D seismic data by classifying 65 x 65 x 65 data cubes using a custom CNN. Zeng et al. [38] proposed a method for salt body identification using the U-Net segmentation network in combination with the ResNet classification network to delineate salt bodies with high precision. Shi et al. [39] approached the problem of salt body extraction as 3D image segmentation and proposed an efficient approach based on deep CNN with an encoder-decoder architecture. They fed a CNN with 128 x 128 patches and predicted a per-pixel salt mask. Wu et al. [40] extended this approach to 3D and directly processed 128 x 128 x 128 data cubes to predict salt bodies. Finally, a paper by Babakhin et al. [41] should be emphasized since it describes the first-place solution on the TGS's challenge [2]. Their key 'ingredient to success' was a semi-supervised method which utilizes unlabeled (test) data for multi-round self-training. In the first round of training, they used only available labeled data, after that, in following rounds they trained the model on the available labeled data and predicted confident pseudo-labels for the unlabeled data. The architecture they employed used ResNet-34 [10] and ResNeXt-50 [42] encoders with scSE modules [43], feature pyramid attention (FPA) module [44] after the last encoder block, and hypercolumns [12] for predicting the segmentation mask.
The approach proposed in this paper uses segmentation CNN, unlike classification approaches presented in [34,36]. The prediction is performed on 2D images instead of approaches [37,40] that process 3D data cubes. It is based on a novel network architecture, that is presented in the next section, which distinguishes it from [35,38,39,41].

Method Description
This section presents the proposed network architecture that came out as a result of the author's participation in TGS sponsored Kaggle competition. Technical details regarding the implementation of the network and the training procedure that was applied are also included in this section.

Network Architecture
The architecture of the network that produced the final, i.e., the best solution provided by the author's participation in TGS's competition, is shown in Figure 4. The presented architecture originates from U-Net [14] and it includes certain modifications made in a search to improve the score.
As can be seen from Figure 4, the network mainly consists of three types of blocks (C, D, and U). C-block is the most common and the most complex block in the network. It consists of 5 convolutional layers with a kernel size of 3 x 3. The block is generated using two parameters: input number of filters (f), and an output number of filters (p). The first 4 layers have the same number of filters (f) and they are 'closely coupled' using addition before ReLU activation. The final, fifth, layer has p filters and its main use is to adjust the number of output filters to the desired value. The specific coupling of the layers' output maps was inspired by ResNet [10] and DenseNet [19] architectures, although the particular 'wiring' is the original contribution. In addition to previously mentioned convolutional layers, C-block also includes batch normalization and ReLU activation layers. The number of filters in convolutional layers of C-blocks is defined using parameter n that represents the number of filters in the first convolutional layers. The experiments were conducted with the following values for n: 16, 24, and 32.
D-block is the second block type that is used, and it can be found in the encoder part of the network always following C-block. The responsibility of D-block is the downsampling of the spatial dimension of the feature map by a factor of 2 using the MaxPool layer. It also includes the Dropout layer with a 20% dropout rate that is used to increase the robustness of the trained model and prevent overtraining. Dropout is only applied during the training phase and it represents annulment (setting values to 0) of a certain percent (dropout rate) of the input feature map. This forces the network to take more input features into account when building higher-level features.
U-block is found in the decoder part of the network and it corresponds in a way to D-block in the encoder part. U-block is used for the upsampling spatial dimension of the input feature map by a factor of 2 which is achieved by duplicating each column and row. The second responsibility of Ublock is a concatenation of a previously upsampled feature map with an output feature map from the appropriate C-block of the encoder (shown by the dashed line in Figure 4).
The final network output is obtained applying a 1 x 1 convolution, with 1 filter, and sigmoid activation to the output of the last C-block. The 1 x 1 convolution is used to reduce the number of filters to the desired network output, while sigmoid activation limits the output values per pixel to the interval (0, 1). The corresponding values are rounded to values 0 or 1 to produce the final output mask.

Implementation and Training
For implementing the proposed architecture, the Python programming language and the Keras [45] library were chosen. Keras is a high-level library that specifies a simplified interface for implementing deep neural networks and in this case, it relies on TensorFlow [46] as a backend engine. The particular network implementation uses Keras functional API, and each of the proposed architectural blocks is implemented in a separate Python function. The source code for the project can be found at https://github.com/a-milosavljevic/tgs-salt-identification.
The training was carried out on 80% of the initial training set, while the remaining 20% was used for validation. More precisely, the data were divided into 5 folds, so 5 networks were trained using the different fold for the validation. The final prediction was obtained using an ensemble of 5 networks were corresponding outputs were summed, divided by the number of the networks in the ensemble, and then rounded to produce a binary mask. To keep validation sets representative compared to the training set, before splitting data into 5 folds, samples were sorted by the number of salt pixels in the mask. After that, images were split by taking 5 samples at a time and assigning them to each of the 5 allocated folds.
The data augmentation technique was applied in both training and test phase (results preparation phase). In the training phase, data augmentation was applied on two levels. The first level represents a doubling of the number of training samples by horizontal flipping all the images and corresponding masks. The second level of data augmentation used random image transforms that included translation, and intensity scaling and shifting. The translation transform relied on the fact that input images are 101 x 101, while the network input is 128 x 128. In 25% of the cases, the image was centered, i.e., it had a translation of 13 pixels in both directions. In the remaining 75% of the cases, the image was translated by random values between 0 and 27 in both directions. The appropriate translation was also applied to the corresponding mask. Besides translation, the unfilled part of the network input (and expected output) was filled by image and mask mirroring. Transformations related to image values include intensity scaling by multiplying image values with a random value between 0.8-1.2, such as intensity shifting by adding random value from the range ±0.2.
Test time augmentation (TTA) works as follows. For each image in the test set network output is evaluated for original and horizontally flipped variation. In the second case, the output is also flipped and averaged with the original output. When not only one, but the ensemble of networks, is used, averaging is done for all network outputs before thresholding operation that produces the final segmentation mask.
For the network training, Adam [47] optimizer (module optimizers) and binary_crossentropy loss function (module losses) were used. As the metric to select the optimal model binary_accuracy function from the metrics module was applied. This metric shows the percentage of correctly classified pixels. Also, the intersection over union (IoU) metric has been tracked and is calculated by dividing the number of pixels in the intersection of ground truth and predicted masks, and the number of pixels in the union of two masks. The IoU metric is typically used for segmentation problems, and the competition metric was based on it, but slightly better results were achieved using binary accuracy for selecting the optimal model.
To control training processes, Keras' callback objects were applied. Appropriate classes are part of the callbacks module and the proposed solution used: ReduceLROnPlateau, EarlyStopping, ModelCheckpoint, and CSVLogger. ReduceLROnPlateau is used to reduce learning rate using a certain factor when in a certain number of epochs certain parameter that is being tracked is not improved. In the proposed solution, binary accuracy on the validation set was tracked, and if there were no improvements after 10 epochs of training, the learning rate was reduced by 0.1. The initial value of the learning rate was 10-3. EarlyStopping callback, as the name implies, were used to end training using a similar mechanism to the previously described one. In the proposed solution, early stopping was activated after 30 epochs when binary accuracy on the validation set does not improve. CSVLogger callback was used to record training and validation losses and accuracies during the training process. Finally, ModelCheckpoint is used to record the best model measured using binary accuracy on the validation set. Each time the accuracy improves a model gets saved, so when the training ends, the model that had the best performance on the validation set is obtained.

Results and Discussion
All Kaggle hosted competitions have clearly defined duration, rules, datasets, and a metric for submission evaluation. In the case of TGS's competition training dataset consists of 4000 seismic image patches and corresponding segmentation masks that are 101 x 101 pixels in size. The test set that is used for model evaluation consists of 18,000 seismic image patches. For all seismic image patches, depth information (in feet) of the sample location is provided. The competition results consist of masks that are estimated for every test image patch. The binary segmentation masks are encoded using the run-length encoding (RLE) and submitted using the standard comma-separated values (CSV) file.
The competition lasted from July 19 to October 20, 2018. Each day it was possible to submit a maximum of five solutions that are scored using a competition-specified metric. The preliminary (public) score is given based on 34% of the test data. The best result of a contestant determines its placement on the public leaderboard. Final results are announced after the competition ends and they are based on the remaining 66% of the test data. The final score is calculated for only two solutions selected by a contestant and the better value is used to place the contestant in a final (private) leaderboard. One week before the competition ends it is possible to merge teams which led to the merging of their individual results.
The main advantage of a competition organized in such a manner is the ability to try out and get an immediate assessment for many ideas in a short time, as well as to measure the quality of oneʹs own solution compared to the other participants with a high level of objectivity.

Competition Results
The results achieved with the proposed network were ranked in the top 14%, i.e., as 446th of 3221 contestants. Before presenting the results, it is necessary to introduce the metric used by the competition to evaluate submissions. The metric was based on averaging IoU values on 10 thresholds ranging from 0.5 to 0.95 with a step of 0.05. For example, on threshold 0.5, the predicted mask is counted as a hit if the IoU value is greater than 0.5.
If T represents ground truth mask, while Y represents predicted mask, on each threshold value t precision P(t) is calculated using the following rule: (1) Finally, evaluation metric M of the predicted mask is calculated by averaging precisions calculated at all 10 thresholds: An overview of the results for different configurations regarding the number of networks in the ensemble and the number of filters in the first convolutional layer (parameter n) are presented in Table 1. Along with private scores that were the basis for the final ranking, corresponding public scores, that were visible during the competition, are also shown. For comparison purposes, the winning solution [41] score is also shown. The total number of submitted solutions for the winning team was 316, while the author had 42 total submissions. In the presented results two trends can be spotted. The first trend is score increase with the number of convolutional filters (parameter n), while the second one suggests that adding more networks to an ensemble also pays off. The best author's submission used 25 networks which may sound surprising. It was a last-minute attempt to increase the score by combining all previously trained models. As a comparison, the winning solution came from the ensemble of 40 networks.
The things that are not visible from the presented results are countless attempts that did not prove successful. These attempts include other architectures, different loss functions and metrics for picking the best model, some postprocessing efforts involving morphology operations, other methods for data augmentation, etc.
As an illustration of training processes, Figure 5 depicts two charts showing the accuracy and IoU metrics on the training and validation sets during one round of training.

Post-Competition Analyses
Judging the model architecture just by comparing the competition results is not the best approach since many other things influence the final score. For that reason, after the competition ended, a few experiments were conducted to see how different segmentation models, using existing training framework, would rank compared to the proposed model. Luckily, thanks to Pavel Yakubovskiy's Github library called Segmentation Models [48], that was fairly easy. Segmentation models are a Keras-based Python library that implements four popular segmentation architectures and dozens of ImageNet pretrained backbones that can easily be combined into the segmentation model of choice. The supported architectures include previously mentioned U-Net [14], such as the feature pyramid network (FPN) [16], LinkNet [17], and pyramid scene parsing network (PSPNet) [18].
When it comes to a backbone selection, the segmentation models library supports a range of models pretrained on the 2012 ILSVRC ImageNet dataset. For the experiments, ResNet-34 [10] backbone is used. Since ResNet network decreases the size of feature maps four times in first two layers, to preserve the precision of predicted output map, input images were upscaled two times by adding UpSampling2D layer at the input and adjust output size by applying AvgPool2D on the output (see Figure 6 for the illustration). A similar method of doubling the size of input images was applied in the winning solution, as described in [41].
A comparison of different segmentation models is presented in Table 2. For each of the models, a total number of convolutional layers and filters-public and private scores are shown. For scoring purposes, ensembles of five networks were used where each of the networks was trained using the different fold as the validation set. The training process, including data augmentation during both training and results preparation phase, was identical to the one used in the competition and previously described in detail in Section 3.2. The first two entries correspond to the results obtained using the network described in this paper. The first entry represents the result obtained during the competition, while the second one was produced by a slightly modified network that excludes dropout layers in D-blocks. Finally, the last four entries correspond to the models implemented using Yakubovskiy's library. FPN, LinkNet, and PSPNet were constructed using default parameters, while the U-Net model was constructed with 512, 256, 128, 64, and 32 filters in each of the decoder blocks which is double the default. This change was introduced to give U-Net similar 'power', in terms of the number of convolutional filters, to the one the original model had. Figure 6. An illustration of the experiment framework based on Yakubovskiy's segmentation models [48] and input image upscaling. Please note that in the case of PSPNet, due to network limitations, an input size of 120 x 120 (upscaled to 240 x 240) was used. Looking at the scores presented in Table 2 it can be noticed that the values do not differ much. The best public score is achieved by the FPN model while the winner in the private score section was the proposed model without dropout layers. Interestingly, looking at public scores, the unchanged original model had a slightly better score. The results also showed that original U-Net architecture proves to be less powerful than the proposed architecture, even though it used a pretrained ResNet-34 encoder as the backbone. The FPN model showed comparable-and even the best result in the public section-employing slightly more convolutional filters that both the proposed and U-Net models. Finally, LinkNet and PSPNet performed slightly worse than the proposed model, but it is worth mentioning that PSPNet uses about two times fewer convolutional filters than the other models.
To further explore how these segmentation models behave, the outputs were visualized and compared. Figure 7 depicts six sample images along with ground truth salt deposits mask and output masks for each of the models. The examples were handpicked to illustrate some easy and some hard cases. The samples were chosen from the validation sets, so the outputs are produced by the networks that have not seen these inputs during the training.
The first three samples can be treated as the easy ones, i.e., all models' outputs show very good matching with the ground truth image. The fourth sample is slightly more complex, so differentiation between the models can be observed, however mostly depicting the 'right shape'. The last two samples are representatives of hard cases where the true mismatch occurs. The ground truth images for both samples show no presence of salt deposits, but the output images have plenty of it. Finding a way to reduce the occurrence of such cases seems like a good research direction. For the comparison with the original model training charts depict in Figure 5, corresponding accuracy and IoU metrics charts are presented in Figures 8-12. Figure 8 depicts charts for the original model without dropout, while Figures 9-12 depicts charts for Yakubovskiy's segmentation models U-Net, FPN, LinkNet, and PSPNet, respectively. Again, the chart data are taken from the training of the fifth model (last partition) and the red line in (a) depicts the best accuracy epoch when the model is saved. Since the same training framework is used, all these charts show similar behavior. The oscillations of the validation metrics are stronger initially, and then stabilizes when the learning rate decreases. It appears that the oscillations are smallest for the original model ( Figure 5) and that is probably due to the dropout that has been applied. If we compare the IoU metric for the training set, it can be noticed that, after the first epoch of training, our models get ~50%, while Yakubovskiy's models go to ~60%. This behavior can be explained by the fact that Yakubovskiy's models use a pretrained ResNet-34 encoder that gives them an initial boost.

Conclusions
In this paper, an approach for the identification of underground salt deposits on seismic images using deep CNN for segmentation is presented. Seismic imagining and salt identification play a significant role in oil and gas discovery, while automatization of the process of their analyses is highly important for the companies involved in hydrocarbon fuel exploration.
This paper came as a result of the author's participation in the TGS sponsored Kaggle competition [2]. The solution that has a score of 0.85241 on the private leaderboard was based on the original architecture of segmentation CNN that was described in the paper. The proposed architecture was derived from the U-Net model that was responsible for the authorʹs initial score of 0.76996. Since the improvement was significant, to better comprehend the properties of the proposed architecture, a series of post-competition experiments were conducted to compare obtained results with the results achieved by standardized approaches using the same training framework. The results show that the proposed approach is comparable and, in most cases, better than these segmentation models. The five models' ensemble of the proposed architecture without dropout layers had the best private score of 0.85219, while the FPN-based ensemble recorded the best public score of 0.83623. It is worth mentioning that segmentation models used in the comparison employed transfer learning, i.e., they used ResNet-34 encoders pretrained on the ImageNet 2012 dataset.
For future work, several directions can be examined. The first one is to apply transfer learning to the proposed architecture by replacing the current encoder with some pretrained network and only keep the decoder part. The second direction would be to try to improve salt identification metrics by training additional classification models that would predict if a patch has salt or not, so the segmentation model would be used conditionally. Finally, the third direction would be to test the architecture on some other dataset and see if it holds up.
Author Contributions: Conceptualization, methodology, validation, formal analysis, supervision, writingreview and editing, data curation, Aleksandar Milosavljević. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.