Adaptive Segmentation for Discontinuity Detection on Karstified Carbonate Outcrop Images From UAV-SfM Acquisition and Detection Bias Analysis

The identification of fractures and discontinuities has great importance on the fluid flow estimation in hydrocarbon reservoirs since they influence the properties of porosity and permeability. Due to the inaccessibility and sparsity of reservoir data, the fracture characterization is generally assessed through the study of outcrop analogues using remote sensing or in situ observations by a specialist. Considering the remote sensing methods, the unmanned Aerial Vehicle (UAV) acquisition combined with Structure from Motion (SfM) photogrammetry is a low-cost way to generate products like orthorectified images, allowing manual and automated methods of fracture trace detection. Automatic approaches, commonly used to address this problem, present some known limitations and disadvantages due to the nature of the outcrops and weather conditions during UAV acquisitions. In this work, we focus on fracture detection over karstic regions that are highly fractured. For this, we evaluated a series of adaptive segmentation methods based on thresholding. The Sauvola local adaptive segmentation presented the best result when compared to a manually annotated ground truth. The segmentation results were further improved by the use of the binary denoising method Non-Local means. We also carried an evaluation of the influence of the sun position in the fracture detection, and to reduce this inherent bias we combined three UAV acquisitions done over the karstic carbonate outcrop, namely Rosário pavement in the Jandaíra formation northeast Brazil. With the proposed methodology we acquired more accurate fracture data over the study area, which follows the directional statistics of previous works carried out in the region.


I. INTRODUCTION
Carbonate reservoirs concentrate important hydrocarbon reserves located at great depths, as in the case of pre-salt in the Brazilian offshore. These reservoirs are mainly affected by diagenesis and post-deposition dissolution caused by karstification [1], which leads to highly fractured rocks contributing The associate editor coordinating the review of this manuscript and approving it for publication was Gerardo Di Martino.
to the heterogeneity of the carbonate reservoirs [2]. The karstification process cause fragility and fracturing of the reservoir rocks changing the inherent properties of porosity and permeability, increasing the importance of the characterization of fracture (and overall discontinuities) in reservoirs to estimate fluid flow [3].
The identification of discontinuities and characterization of the discrete fracture network (DFN) in oil reservoir rocks is generally done by seismic or well logs, where the first suffers from low spatial resolution, and the latter, while having a much better resolution, suffers from the sparsity of the collected data [4]. To overcome the inherent uncertainty of these kinds of data acquisition, the fracture characterization can be carried using analogue outcrops i.e. outcrops with similar characteristics to the reservoir [4], [5].
The fracture characterization in analog outcrops is usually carried out in situ. However, remote sensing techniques powered by LiDAR (Light Detection And Ranging) [6], [7] and/ or Unmanned Aerial Vehicles (UAV) image acquisitions, combined digital photogrammetry algorithms like the Structure from Motion (SfM) and Multi-view Stereo (MSV) [8]- [10], can be used to create digital representations of the outcrop such as orthorectified images and digital outcrop models (DOM) that enable geoscientists to perform fracture characterization in a virtual environment [5], [8], [11], [12]. The UAV-SFM-MVS photogrammetry processing is a low-cost method for outcrop geometry reconstruction. It identifies key point pairs in overlapping images to estimate camera positions using the bundle adjustment.
Regardless of the selected surveying approach, the fracture characterization for DFN modeling requires the determination of fracture network geometric parameters like length distribution, spatial density, fracture aperture, and orientation [13]. Additionally, fracture network topology parameters like joints and nodes are also assessed to estimate fracture network connectivity [14].
Sampling strategies either in the field or in digital models include the scanline one-dimensional sampling, and/or circular or window sampling [15], where the fracture characterization can be carried out manually, or by assisted or automated tools. Manual characterization is performed by visual interpretation and digitization of fractures over an outcrop or images and 3D models. However, manual interpretation can be time-consuming [16] and be affected by ambiguity and subjectivity biases [17], being the main motivation for assisted and automated methods [18], which detect and quantify fractures using 3D geometry and feature detection in images.
Fracture detection in images is an active field of research, taking advancements in computer vision and image processing techniques also applied to other field areas like medical applications, road identification, and civil engineering. For this, signal analysis methods like Phase Congruence [19], Wavelet [20]- [22], Shearlet [23], [24], line detection with Hough transform [25], [26], ridge and edge detection filters with Canny [19], [25]- [29] and Steger algorithm [26], [29] have also being applied in fracture detection. With most of these methods generating noisy images when detecting discontinuities in highly fractured outcrop or with irregular shadowed areas and vegetation. Deep learning methods based on Convolutional Neural Networks (CNNs) for semantic segmentation were recently developed to successfully solve some of these problems when detecting fractures in outcrop and rock images [30]- [34], however, extensive and specialized training datasets are needed to obtain the best results. Another caveat is the interpretation bias in the automated fracture detection caused by the sunlight direction when detecting shadow areas as fractures, a problem also pointed by [26] but, to our knowledge, no other work has considered this when generating the DFN. Aiming to tackle these drawbacks, we propose the use of a simple adaptive thresholding method alongside the fast denoising technique Non-Local means [35]. We also evaluated and combined three UAV acquisitions over the Rosario pavement, a karstic carbonate outcrop in the municipality of Felipe Guerra in the State of Rio Grande do Norte, Brazil. This combined acquisition and following fracture detection are then used to estimate fracture intensity, fracture length distribution, and fracture aperture. With this strategy, we aim to capture a more accurate representation of the complex fracture network of the karst environment.

II. METHODOLOGY
The methodology employed in this work and the general workflow is presented in Figure 1 and better detailed below. The workflow consists of UAV acquisitions, photogrammetry processing, data acquisition and validation of adaptive segmentation methods, and the use of a binary denoising method. Once the segmentation method is validated, we perform the fracture characterization and DFN generation using multiple UAV acquisitions to analyze the bias in the fracture detection. Finally, we combine the fracture segmentation from each UAV acquisition to properly characterize the DFN, evaluating the fitted distribution of the fractures features and computing spatial fracture indices.

A. STUDY AREA AND GEOLOGICAL SETTINGS
The sudy was carried out in the municipality of Felipe Guerra in the Rio Grande do Norte state, northeast Brazil ( Figure 2). The area of interest consists in an extensive outcrop, the Lajedo Rosário, where rocks of the Jandaíra Formation (Potiguar Basin) crop out [36], [37].
The Jandaíra Formation consists of an extensive succession of carbonatic rocks, mostly of limestones, such as grainstones, packstones, and secondarily mudstones, with variable degrees of dolomitization [38], [39]. The strata present a variable bioclastic content of bivalve mollusks and algae [38], [40]. The corresponding sediments deposited from the Turonian (± 93.9 a 89 Ma) to the Campanian (± 83.6 Ma) [41] as sea carbonate sediments within turbulent to shallow water environments [42]. The strata show a gentle dip (2-5 • ) towards NNE and surface as prolonged NE-SW plateaus in a regional level [39]. These carbonate banks occur with distinct thicknesses across the basin, from 30m in the Upanema Region and above 600m in the Apodi region [36]. The basal parts of the formation show a transition of the calcareous sediments to the sandstones of the underlying Açu Formation, which indicates a marine transgression in the basin. Thus, the Jandaíra formation relates to the transitional rift stage that ends the continental phase [43]. The outcropping Jandaíra Formation presents extensive deformation by fractures [39], [44], [45]. The literature reports main fracture families running N-S and, E-W directions, with subsidiary NE-SW and NW-SE-striking families, in both cases presenting sub-vertical dip, and, lastly, a set of sub-horizontal, bedding-parallel fractures [44], [45]. Calcite veins and open fractures with N-S, NNE-SSW, and E-W directions have been observed [44], and E-W and WNW-ESE directed stylolites [39]. These structures relate to diverse stages of subsidence and regional lateral movement related to the context of the Atlantic margin opening [39].

B. IMAGE ACQUISITION
Image acquisition was carried using UAVs DJI Mavic Pro 2 and a DJI P4 multispectral in multiple flights over the study area shown in Figure 3 to make it possible to analyze the impact of the time of the day in the acquisition and fracture detection using the proposed method in this work.
The DJI Mavic Pro 2 is a UAV with a 20M effective pixel embedded RGB sensor, 29 minutes maximum flight time and 907g (2lbs), while the DJI P4 multispectral is a UAV with 6 embedded cameras, 1 RGB sensor and 5 monochromatic sensors capturing five distinct wavelengths with 2M effective pixels in each sensor, weighing 1.47kg (3.2lbs). The DJI P4 also includes a Real Time Kinematic (RTK) Global Navigation Satellite System (GNSS).
The UAV image acquisition is carried at three different times of the day respecting the 60% to 80% overlap between image acquisitions. The SfM-MVS processing is carried using the Agisoft Metashape photogrammetry software version 1.6 considering medium settings of quality for photo alignment and cloud densification.
The main products used in this work are the georeferenced orthorectified images RGB data in the region shown in detail on Figure 3. These products are aligned and resized to the same resolution using visual key points to allow the fracture detection comparison. The sun position during each flight is registered using an online tool 1 that computes the sun position given the date-time and global position.

C. ADAPTIVE THRESHOLDING SEGMENTATION
Thresholding techniques, also known as binarization, are well-known techniques in character recognition and segmentation of Magnetic Resonance Imaging (MRI), being regionbased segmentation methods that aim to separate objects in the foreground from the background. Fracture identification can make use of these methods as fractures and discontinuities naturally appear as dark areas in outcrop images.
A binary or thresholded image (b) is given by where (x, y) represents the pixel position, I (x, y) is the grayscale intensity value and T is the threshold value that is determined by a chosen method. For adaptive thresholding, mean and standard deviation estimation, as well as contrast estimation, are computed over a moving window that is convolved over the image. In this work, our strategy is to apply local adaptive thresholding to maximize the fracture segmentation in images with a high level of detail or arbitrary shadowed areas like images of karstified carbonate outcrops with vegetation.
For this, we evaluate several adaptive thresholding techniques using those proposed by Niblack, Sauvola, and Phansalkar, which are common techniques for segmentation of biomedical images and character recognition. For comparison, the local Gaussian filter and the global thresholding filter Otsu are evaluated. These thresholding techniques are shortly described.
The Niblack local thresholding [47] T (x, y) is given by where m(x, y) is the mean value and δ(x, y) is the standard deviation of a window of pixel values over a convolution window of size w, while k is a bias value set to 0.5. Adapted from [40].
The Sauvola adaptive thresholding [48] is given by where m(x, y) and δ(x, y) are the same aforementioned. R is a standard value for the grayscale midpoint normally set to 128, and k is a bias value in the range 0.2 to 0.5. The Phansalkar adaptive thresholding [49] updates the Sauvola technique by trying to better estimate the threshold in dark regions with low contrast by adding two new constants and an exponential operation that add more impact to the mean. This local thresholding is given by where the constants p is set to 3, q is set to 10, and k is set to 0.25 following the original paper.
Otsu's method [50] uses the intensity distribution (histogram) to iteratively find the threshold value. In the first step, the algorithm tries to find two classes where the mean for each class is calculated and used to find two threshold values and three classes: background, foreground, and ''to be determined'' (TBD). The method is applied again in the TBD class and in the last step, the new TBD area is divided into two classes finding the global threshold value.
For the adaptive local thresholding techniques, the size of the moving is set to a fixed size of 91 × 91 pixels aiming to obtain the threshold for an area large enough to cover the fracture width in our images.

D. BINARY DENOISING METHODS
After segmenting, the outcrop images with the adaptive thresholding methods, some noise related to fracture areas incorrectly detected or to other artifacts are expected. In this stage, we evaluated common techniques for noise reduction in binary images like Graph cut, Markov Random Fields (MRF), and Non-Local Means (NLM).
The Non-Local means (NLM) [51] algorithm is a denoising method that uses the similarity within the image space, inferring the true value of a pixel in a region analyzing other regions in the image that have similar pixel distribution. Its equation considers the noisy image as an element v, while x is the element in the image that we want to determine the average, based on Gaussian neighborhoods similar to the neighborhood of x. In its base form of the NLM is given by where G α is a Gaussian kernel with a standard deviation alpha, C(x) is a normalizing factor, and h is a filtering parameter. In a discrete image v = {v(i)|i ∈ I } the NLM of an element i in the noisy image v is given by where the weights {w(i, j)} j depend on the neighborhood similarity of the pixels i and j.  The NLM denoising algorithm we employed is based on the work of [35] with computing optimizations made available by the OpenCV library.

E. SKELETONIZATION/FRACTURE CENTER PROFILES
To identify discontinuity traces we need to extract the centerline or medial axis of the segmented fractured areas using a process called skeletonization. The skeleton of an object can retain its topological, geometrical, or both forms of the original object, depending on the technique used which is based on distance transformation, Voronoi diagram, or thinning.
Algorithms like the developed by [52] and [53] are traditional skeletonizing algorithms based on iterative thinning, e.g. iteratively remove pixels around objects until they have a 1-pixel width. In this work, however, we employed the distance transformation called Medial Axis Transform (MAT) to obtain the fractures aperture.
The MAT algorithm uses a 3 × 3 filter pass evaluating how many pixels are connected to the center pixel, and depending on this evaluation, the pixels are marked for removal or not. In the next iteration of the algorithm, the pixels are removed and the distance transform is computed for each object center pixel to its border, where the Euclidean distance is computed.  if neighborhood center pixel plus another pixel are marked then 4: Mark the pixel as a termination 5: end if 6: if neighborhood center pixel plus another tree pixels are marked then 7: Mark the pixel as a cross node 8: end if 9: end for 10: for every node in terminations do 11: Vectorize fracture pixels using region growing 12: end for 13: for every node in cross nodes do 14: Vectorize fracture pixels using region growing 15: end for 16: Simplify vectorized fracture traces to lines The filter pass is illustrated in Figure 6, where we captured two common situations when convolving the 3 × 3 neighborhood. After marking these node regions, we use them as initial points to a region-growing algorithm that will walk along with the center fracture profile.

G. IMAGE COMPARISON METRICS AND VALIDATION DATA
For image segmentation evaluation, this work uses orthorectified outcrop images available in [54] that were obtained from UAV acquisition (height of 60 to 70m) and SfM photogrammetry processing generating products with 1cm/pixel resolution. A region of this outcrop image was manually annotated considering the most prominent fractures and discontinuities both pixel-wise segmenting its shape and interpreting its delineations with GIS tools in a primary study carried by the Vizlab 2 laboratory aimed to solve the lack of public datasets of annotated fractures for appropriate comparison between works. The pixel-wise annotated images and fracture delineations available on GitHub 3 are used in this work for direct comparison of the delineaments and validation of the segmentation methods.
For comparison, a set of evaluation metrics were used to verify the proximity of the segmented images to the image manually segmented, which was used here as the ground truth. This set of metrics consists of the Structure Similarity Index Measure (SSIM), the mean Intersection Over Union (mIoU), the F1 score or F-measure, and the Peak signal-tonoise ratio (PSNR).
where x and y are the pixel values in both images, µ x and µ y are the average of x and y respectively, σ 2 x and σ 2 y are the variance of x and y respectively, σ xy is the co-variance of x and y, C 1 and C 2 are the function stabilizers given by (K 1 L) 2 and (K 2 L) 2 respectively where K 1 and K 2 are constants and L is the dynamic range of the image.
The PSNR measures how close two images are to each other, where higher values indicate higher similarity. The PSNR is given by and M and N are the pixel size dimensions in the x and y axis, and I and I are the first and second image, respectively.
The F1 score or F-measure and the mIoU are based on the confusion matrix elements: TP for true positives; FP for false positives; TN for true negatives and; FN for false negatives.
The mean Intersection over Union (mIoU), also known as Jaccard similarity index, is given by The F1 score is given by where and

H. DISCRETE FRACTURE NETWORK AND FRACTURE STATISTICS
Once the fracture network is mapped and digitized, the properties of strike/dip orientation, fracture intensity, length distribution, aperture distribution, and topology that define the geological DFN model can be estimated.
Considering fracture characterization based on fracture traces, the most common index for outcrop fracture characterization is the one based on the P index [15], [56]. These properties are directly related to the dimensions of the observed data and the number of dimensions considered for its estimation, given by the notation P xy where x is the dimension size of the sample and y is the number of dimensions used in the analysis. Figure 7, based on the work of [57], brings the different indices (P xy ) used to determine the fracture intensity from the measurement and sample dimensions.
From the discrete fracture network, the attributes like length distribution can be sampled for the entire fracture network or subsets, depending on the orientation of the main fractures previously identified or on the multi-modal analysis of the rosette histogram and additional circular statistics [14].
The distributions of the lengths of natural fractures generally follow the Power Law [13], [15], [59] distribution, although the log-normal, negative exponential and gamma distributions can also be found [14], [15], [60], [61]. In addition, the fracture aperture is also assessed by the distribution of the sample measurements following log-normal, power-law, and exponential distributions as in the work of [60]. In this work, we consider the distributions powerlaw, exponential, and log-normal in its probability density function (PDF) form to evaluate the best fit in a log-log plot.
The PDF of the Exponential distribution is given by where λ is estimated byλ = 1/x, with x ∈ X , and X is the set of attributes measured. The PDF of the Log-normal distribution is given by where µ is the mean value of the distribution, and σ is the standard deviation. The PDF of the Power-Law distribution is given by where α is the estimated fractal dimension [21] given bŷ To evaluate the distribution of the lengths and aperture data, a common method is to plot the fracture data in the log-log plot as an inverse y-axis against the cumulative PDF of the distribution that we want to measure the similarity in the x-axis. The x and y positions in this plot are used to determine the best fit of the data applying the least-squares linear regression algorithm, and the best fit is evaluated using the coefficient of determination R 2 .

I. IDENTIFYING THE SEGMENTATION BIAS AND GENERATING THE FINAL DFN MODEL
To analyze a possible bias from the UAV acquisition and fracture detection, considering different conditions or acquisition times, we analyze each fracture segmentation product related to each product. The analysis is carried firstly by analyzing the weighted rosette diagram, with each bar being a sum of lengths in a given direction, where we expect to find slight differences in the resultant directional statistics.
To obtain a more accurate representation of the fracture network from the outcrop images, we combine multiple fracture segmentations. This is done by accumulating the denoised fracture segmentations generating a new binarized image. This binarized image will then follow the process of skeletonization/discontinuity detection to obtain the resultant DFN model with length and aperture statistics.

III. RESULTS AND DISCUSSION
Initial evaluation, as presented in methodology, was performed in the outcrop region using the annotated ground truth. Preliminary results, as in Figure 8 and Table 1, showed better segmentation results for the Sauvola local threshold method and will be used as base adaptive segmentation method.
These results were obtained using a window size of 91 × 91 pixels, as we found out that, using the standard  . Evaluation of arbitrary window sizes for the Sauvola. In (b) using a small window the fracture detection is noisy and with artifacts when the fracture areas surpasses the window size, while we detect fracture elements with geometry closer to the areas in black (a) as we increase the window sampling size.
values for the window size, the fractures aperture could surpass this window in width, segmenting fracture regions incorrectly. This is exemplified in Figure 9, where we evaluate distinct window size widths for the Sauvola technique.
The reason for this is due to the impact of the standard deviation in the threshold level, areas with low standard deviation  conduce, in general, to a low threshold level, setting the area as foreground (white). This also explains why vegetation areas were suppressed from the segmentation as they are more uniform having low standard deviation than the fractured areas.
Although presenting the best results, the Sauvola technique seems to be still noisy in visual comparison to the annotated ground truth. As we presented in the methodology the Non-Local means method was the chosen technique for this task. Figure 8f shows the thresholded image after applying the denoising techinique showing a closer image to the ground truth than the other methods evaluated. The values of mIoU, SSIM, F-score, and PSNR increased to 0.91, 0.85, 0.95, and 59.13 respectively, as shown in Table 1 for the Sauvola/NLM approach. VOLUME 10, 2022  These results show the efficiency of the proposed method when detecting fractures in outcrop images. Although not directly compared with other methods already explored in related works, they showed promising results and a research opportunity when visiting computer vision methods often overlooked by recent works.
In terms of segmentation results against related works that used images of highly fracture karstified outcrop, the work of [62] presents an approach using global thresholding which, as shown in our results using Otsu, mistakenly segments trees as fracture areas, forcing the research of [62] to work in a small area of the outcrop without vegetation. Also, the trace detection is carried using Hough transform, which requires several initial parameters, while our work requires only the window sample size.
We also carried a preliminary study in this region employing semantic segmentation methods using deep learning [34], that although it presented lower results of mIoU (0.79 against 0.91 of this work), as pointed in the works of [30]- [33], this results might have been influenced by the quality and quantity of the training dataset.
After evaluating the proposed method using the validation set presented in Section II-G, we proceed with the evaluation of the bias of the sunlight position and its influence on fracture detection. In this step, multiple flights over the Rosario pavement were executed using the DJI Mavic Pro 2 and the DJI Phantom RTK as detailed in Section II-B.
The flights with the DJI Mavic Pro 2 were carried out in the morning (10 AM) and in the afternoon (2:34 PM) lasting 35 to 40 minutes, recovering 443 and 431 images respectively. The photogrammetric processing on the Agisoft Metashape (illustrated in Figure 10) generated orthophotos with 2.27 cm/pixel and 2.58cm/pixel. The flight with the DJI Phantom, made in the afternoon (1:46 PM), recovered 605 images, with photogrammetric processing generating an orthophoto with 6.69 cm/pixel resolution.   Figure 11 shows the orthophotos and fracture characterization on a section of the outcrop considering the three UAV acquisitions, showing the differences in color tone, light, and shadows. These difference also influences the fracture characterization and directional statistics in each case, showing that more fractures were detected in the perpendicular direction to the sun direction, which is observed when analyzing the weighted rosette diagram both in the number of discontinuities detected and main directions highlighted.
This confirms the hypothesis that the position of the light source can inflict an undesired bias in the fracture detection as also pointed in the work of [26], however, this bias is only observed in a vertical outcrop. Other works like [63] point the importance of capturing the outcrop images with the UAV during sunny weather opposing the recommended cloudy conditions ideal for general purpose photogrammetry. In addition, [28] highlights the importance of the light source being in an oblique position to the outcrop face that we want to detect discontinuities using shadowed areas.
To overcome the inherent bias of choosing one of these acquisitions for the DFN modeling, the combination of these fracture segmentation results as we expect, could lead to a more accurate representation of the fracture network. Figure 12 shows the result of the detection and fracture characterization to obtain the DFN using the combined segmentation of the three UAV flights after the procedures of skeletonization by medial axis transform for aperture analysis (Figure 12a) and line detection (Figure 12b). The weighted rosette diagram shows four main directions of fracture discontinuities: N-S, E-W, NE-SW, and NW-SE.
These main directions identified align with the main fracture directions observed in the Jandaíra formation (N-S and E-W) and its striking families (NE-SW and NW-SE) [44], [45].
The log-log distribution fitting for both aperture data and fracture length data are presented in Figures 13 and 14, respectively. As expected for this kind of fracture data, the length distribution is closer to the power-law distribution where the coefficient of determination R 2 is 0.97. For aperture data, the fitting against the log-normal distribution presented a better R 2 score, 0.97.
As final results, the Figure 15a shows the estimated fracture intensity index P 21 using sample areas of 5m, and Figure 15a shows an interpolation of this fracture intensity data.

IV. CONCLUSION
Fracture segmentation in images is still a challenging task in geology applied computing, which can be observed by the different number of methods developed over the years. Due to the nature of the techniques and the great differences between scenarios like outcrop composition and geometry, light conditions, and accessibility of the chosen sensing technique, the automated fracture detection and statistic interpretation of fracture attributes can suffer from the existence of false negatives and false positives during the segmentation process.
In this work, we proposed the use of simple adaptive thresholding methods combined with a fast binary denoising method. In our evaluation, the Sauvola technique performed better than other adaptive thresholding methods, by being more accurate in the fracture segmentation also avoiding noisy areas like areas with vegetation. We also presented a vectorization method capable of extracting fracture lines using topology elements in the fracture network.
Additionally, the weather conditions can greatly influence the fracture detection inflicting an undesired bias in the fracture attribute analysis hypothesis that was confirmed by the evaluation of multiple UAV acquisitions. The combination of these acquisitions led to more accurate data for the analysis of the fracture attributes and a better framework for the fracture network characterization.
Concerning more recent approaches to detect fractures in images, like the deep learning methods using CNN, this work and related works that rely on traditional computer vision methods, have some advantages like the possibility of processing the entire image without needing to divide the images into small patches due to memory and processing limitations when training CNNs. Albeit showing better results in most cases the CNNs for fracture segmentation need datasets that must encompass most of the expected scenarios, and the image processing and computer vision methods can contribute to increasing usable datasets while providing ready to use results for fracture characterization.
Future works in this research aim to expand the DFN characterization by analyzing the fracture detection bias regarding topology while exploring stochastic modeling and reservoir fluid flow simulation. We also point that an extensive evaluation of fracture segmentation methods is yet to be seen in the literature.