An automated fracture trace detection technique using the complex shearlet transform

Representing fractures explicitly using a discrete fracture network (DFN) approach is often necessary to model the complex physics that govern thermo-hydro-mechanical–chemical processes (THMC) in porous media. DFNs find applications in modelling geothermal heat recovery, hydrocarbon exploitation, and groundwater flow. It is advantageous to construct DFNs from the photogrammetry of fractured outcrop analogues as the DFNs would capture realistic, fracture network properties. Recent advances in drone photogrammetry have greatly simplified the process of acquiring outcrop images, and there is a remarkable increase in the volume of image data that can be routinely generated. However, manually digitizing fracture traces is time-consuming and inevitably subject to interpreter bias. Additionally, variations in interpretation style can result in different fracture network geometries, which, may then influence modelling results depending on the use case of the fracture study. In this paper, an automated fracture trace detection technique is introduced. The method consists of ridge detection using the complex shearlet transform coupled with post-processing algorithms that threshold, skeletonize, and vectorize fracture traces. The technique is applied to the task of automatic trace extraction at varying scales of rock discontinuities, ranging from 10⁰ to 10² m. We present automatic trace extraction results from three different fractured outcrop settings. The results indicate that the automated approach enables the extraction of fracture patterns at a volume beyond what is manually feasible. Comparative analysis of automatically extracted results with manual interpretations demonstrates that the method can eliminate the subjectivity that is typically associated with manual interpretation. The proposed method augments the process of characterizing rock fractures from outcrops.

1 Introduction

Naturally fractured reservoir (NFR) modelling requires an explicit definition of fracture network geometry to accurately capture the effects of fractures on the overall reservoir behaviour. The National Research Council (1996) suggested the idea of using geologically realistic outcrop fracture patterns to guide subsurface fracture modelling. In recent work, the use of deterministic discrete fracture networks (DFNs) based on trace digitization from the photogrammetry of outcrop analogues was investigated by Bisdom et al. (2017) and Aljuboori et al. (2015) for reservoir fluid flow simulation and well testing. Outcrop-derived DFNs encapsulate 2-D fracture network properties at a scale that cannot be characterized using either standard surface approaches (scanlines and satellite imagery) or subsurface techniques (seismic imaging/borehole imagery/core sampling). Ukar et al. (2019) suggested a comprehensive set of protocols to select fractured outcrops that are representative of the subsurface. Stochastic and geomechanical DFNs are alternatives to outcrop-derived DFNs for fractured reservoir modelling. Stochastically generated DFNs have the disadvantage that they cannot replicate the spatial organization of fracture network patterns observed in nature (Thovert et al., 2017). Geomechanically derived DFNs are based on the physics of fracture propagation (e.g. Olson et al., 2009; Thomas et al., 2018) and can reproduce realistic fracture patterns provided the complex paleostress field and paleo rock properties are known; however, they are also computationally intensive and hence have limited applicability. A carefully chosen fractured outcrop that is relatively free of noise (fractures resulting from exhumation and weathering and not hidden too much by vegetation) may be used to interpret realistic fracture networks which are geometrical inputs used in simulating various subsurface thermo-hydro-mechanical–chemical (THMC) processes.

Recent advances in unmanned aerial vehicles (UAVs) and stereo-photogrammetry has dramatically simplified the acquisition of georeferenced datasets of fractured outcrop images (e.g. Bemis et al., 2014; Harwin and Lucieer, 2012; Turner et al., 2012). Photogrammetry using the structure from motion (SfM) principle is a relatively inexpensive and rapid technique by which 3-D outcrop models are built by identifying, extracting, and positioning common points in georeferenced outcrop images (Donovan and Lebaron, 2009). Images are captured using a camera-equipped UAV that is capable of following pre-programmed flight missions where flight path, altitude, velocity, and overlap are specified. The images undergo further processing steps that include generating sparse point clouds of common points, aligning the images, generating dense point clouds (3-D representation of outcrop geometry), and generating meshed surfaces (Bisdom et al., 2017). Interpreting fractures on the image orthomosaics with conventional Geographic Information System (GIS) software completes the outcrop-based DFN workflow.

Manually interpreting fractures is time-consuming and forms a bottleneck in an outcrop-based DFN workflow. A manual interpretation has a fair degree of associated subjectivity, and interpreter bias may take the form of specific scales of features being inadvertently omitted or deliberately ignored (Bond et al., 2007; Scheiber et al., 2015). Manual interpretation also suffers from a lack of repeatability owing to the level of expertise of the interpreter, and the interpretation criteria followed (Hillier et al., 2015; Sander et al., 1997). Reproducibility may not be guaranteed even with the same interpreter in multiple trials (Mabee et al., 1994). According to Bond et al. (2015), quantifying the magnitude and impact of subjective uncertainty is difficult. Long et al. (2018) conducted a study on the variability of fracture interpretation in which geologists with varying levels of expertise interpreted a single image. They found considerable variation in fracture topology, orientation, intensity, and length distributions in the interpretations. Andrews et al. (2019) made a detailed quantification of subjective bias in scan-line-based fracture data collection, the associated effects on derived fracture statistics, and suggested protocols for managing the variations. Peacock et al. (2019) delved into the multiple reasons for bias and the resulting implications for modelling. Given the amount of data generated in short UAV flight missions, man-hours spent in interpretation, and the need to de-bias interpretation as much as possible, automatic feature detection techniques may be considered. Automated approaches can speed up the process, improve accuracy, and exploit the acquired data to the fullest possible extent.

In this paper, we apply an automated method to extract digitized fracture traces from images of fractured rocks. The method utilizes the complex shearlet transform measure to extract fracture ridge realizations from images. Post-processing image analysis algorithms are coupled with the ridge extraction process to vectorize fracture traces in an automated manner. The complex shearlet transform was introduced by Reisenhofer (2014) and King et al. (2015) and previously applied to problems such as detecting coastlines from synthetic aperture radar (SAR) images (King et al., 2015) and propagating flame fronts from planar laser-induced fluorescence (PLIF) images (Reisenhofer et al., 2016). We present automatic fracture extraction results from drone images of two carbonate outcrops (Parmelan, France and Brejões, Brazil) and station-scale images of igneous dyke swarms.

2 Background

2.1 Review of automated and semi-automated fracture detection approaches

Rapid digitization of geological features from photogrammetry is challenging owing to issues like spatially varying image resolution, inadequate exposure, the presence of shadows due to the effects of topography on illumination conditions, and chromatic variations in essential features. False positives are non-geological features (such as trees, shrubbery, and human-made structures) that are detected using semi-automated or automated approaches (Vasuki et al., 2014). The removal of false positives is time-consuming. On the other hand, essential features that are not detected at all (referred to as false negatives) by an algorithm further complicate the task of automated feature extraction. Automated methods, in general, detect more features than what is present in the image (Abdullah et al., 2013). In this section, we review some approaches for automatic fracture detection based on the class of algorithm used.

Automated fracture detection utilizing higher-dimensional data such as point clouds, digital elevation models (DEMs) and digital terrain models (DTMs) have an advantage in that depth variations are captured and can be used to extract features. Thiele et al. (2017a) presented an approach based on a least cost function algorithm applicable to ortho-photographs of jointed fracture sets and 3-D point cloud data. Masoud and Koike (2017) introduced a software package to detect lineaments from composite grids derived from gravity, magnetic, DEMs, and satellite imagery. Bonetto et al. (2015, 2017) presented semi-automatic approaches that extract lineaments from DTMs utilizing the curvature of geological features. Hashim et al. (2013) presented an edge detection and line linking method using enhanced thematic mapping (ETM).

Colorimetry of an image can be used to detect features. By partitioning features in the image, e.g. matrix rock, as lighter shades of grey and fractures as darker shades of grey, fracture pixels may be extracted separately from matrix rock using pixel values. Vasuki et al. (2017) developed an interactive colour-based image segmentation tool using superpixels (Ren and Malik, 2003), which are groupings of pixels that are perceptually similar.

Edge detection techniques identify points in images where sharp changes in image intensity occur. Some of commonly used edge detection techniques in image processing are Canny, Sobel, Prewitt, Robert, Kuwahara, and Laplacian of Gaussian filters. Alternatively, edges may be detected using methods that are invariant to contrast and illumination in images. Phase symmetry and phase congruency algorithms (Kovesi, 1999, 2000) fall under this category. Phase symmetry is an edge detection technique that is invariant to local signal strength. The method works identifies the axis of a feature by isolating pixels symmetric along profiles that are sampled from all orientations except parallel to the feature. The axes of symmetry are regions where frequency components either approach a maximum or minimum. The phase congruency method is another edge detection method that detects features by identifying points where Fourier components are maximally in phase. This approach is also invariant to the magnitude of the signal. The property of invariance enables the identification of structures within the image even in the presence of noise. Vasuki et al. (2014) utilized an edge detection algorithm using the phase congruency principle coupled with a multi-stage linking algorithm for the detection of fault maps.

The Hough transform (Duda and Hart, 1972) is another technique that has been used to detect lineaments in images. The Hough transform identifies pixels in binary images that are likely to represent rock fractures using a voting procedure. Each pixel in a binary image is represented as a sinusoidal curve in a 2-D parametric space (or a Hough space). The voting procedure accumulates a vote for each curve in the parametric space corresponding to each non-zero pixel in the binary image. The curves with the highest votes are selected as probable fractures since they correspond to the largest number of non-zero pixels. Results by Callatay (2016) using the Hough transform for fracture detection report the following limitations. Firstly, the detection is limited to a given fracture orientation set owing to the definition of the Hough transform parameter space. Secondly, the issues of false positive detection and discontinuities persisted. The method is also limited by the fact that it needs a binarized image to start.

The development of wavelet theory in the field of harmonic analysis has led to applications in edge detection (Daubechies, 1992; Heil et al., 2006). Mallat and Hwang (1992) proposed wavelet-based approaches for edge detection. Wavelet-based methods differ from gradient-based edge detection methods that search for local maxima of the absolute value of the gradient. Felsberg and Sommer (2001) introduced monogenic wavelets for the purpose. Tu et al. (2005) considered the use of the magnitude response of complex wavelet transforms. Wavelets, owing to their isotropic properties, cannot extract curve-like features due to the lack of directional information (Labate et al., 2005). A number of wavelet-based approaches that have been proposed to overcome this lack of directional information such as curvelets (Candès and Donoho, 2005), ridgelets (Candès and Guo, 2002), contourlets (Do and Vetterli, 2005), bandlets (Le Pennec and Mallat, 2005), wedgelets (Donoho, 1999), shearlets (Guo et al., 2005), and band-limited shearlets (Yi et al., 2009).

2.2 The complex shearlet transform

In images of fractured outcrops, the presence of discontinuous gaps due to rupture within the rock mass, which occur naturally and which maybe enlarged through weathering processes, is commonly used as a defining criterion by interpreters to digitally trace and classify fractures within the rock mass. Fractures may also be partially or completely sealed by the presence of infilling material that maybe mineralogically different from the adjacent rock material. In such a case, the contrast in colour and texture of the infill material provides an interpretative criterion for classification of these material regions as fractures. The presence of such prominent discontinuities within otherwise smooth regions of rock images can be exploited by the complex shearlet transform to precisely identify positions in the form of edges and ridges.

The basis of the complex shearlet transform applied to fracture extraction from images emanates from wavelet theory. Wavelets are rapidly decaying wavelike oscillations possessing a finite duration. Wavelet transforms are routinely used in digital signal processing applications, which are often time–domain signals. They can also be applied to image data which can be regarded as 2-D functions. Wavelet transforms are not able to detect the directionality of structural features in image data since they may only be dilated or translated. Shearlets that were introduced by Labate et al. (2005) as a new class of multidimensional representation systems overcame a major shortcoming of wavelets by enabling dilation, shear transformation, and translation operations. The isotropic dilation of wavelets was replaced with anisotropic dilation and shearing in the case of shearlets. These modifications have resulted in shearlets possessing a number of properties that make them better suited to handling sparse, geometric features in 2-D image data compared to traditional wavelets (Kutyniok and Labate, 2012).

The complex shearlet transform is a complex-valued generalization of the shearlet transform that was developed by Labate et al. (2005) to handle geometric structures in 2-D data. Reisenhofer (2014) and King et al. (2015) proposed the idea of creating complex shearlets by modifying the shearlet construction so that real parts of the generating function are even-symmetric and imaginary parts of the generating function are odd-symmetric. They used the Hilbert transform to convert even-symmetric functions into odd-symmetric ones and vice versa. The complex shearlet measure for ridge and edge detection implemented in Reisenhofer (2014), Reisenhofer et al. (2016), and King et al. (2015) merged the ideas of phase congruency (Kovesi, 1999) and complex shearlets.

The complex shearlet measure first introduced by Reisenhofer (2014) and improved by King et al. (2015) was used for applications like coastline detection (King et al., 2015), flame front detection (Reisenhofer et al., 2016), and feature extraction from terrestrial lidar inside tunnels (Bolkas et al., 2018). Karbalaali et al. (2018) used the complex shearlet transform for channel edge detection from synthetic and real seismic slices. Reisenhofer et al. (2016) presented a comprehensive comparison of complex shearlet-based feature detection compared with conventional edge detectors such as Canny (Canny, 1986), Sobel (Sobel and Feldman, 1973), phase congruency (Kovesi, 1999), and another shearlet-based edge detector (Yi et al., 2009). Bolkas et al. (2018) also made specific comparisons between the performance of Canny, Sobel, and Prewitt (Prewitt, 1970) edge detection methods versus space–frequency transform methods such as wavelets, contourlets, and shearlets. A detailed overview of the complex shearlet transform is provided in Appendix A for the interested reader.

3 Methods

3.1 The automatic detection process

The automated fracture trace detection method that we present has five main steps (see Fig. 1). The first step of the method uses the Complex Shearlet-Based Ridge and Edge Measure (CoShREM), a MATLAB implementation by Reisenhofer et al. (2016) that utilizes functions from Shearlab3D developed by Kutyniok et al. (2016) and Yet Another Wavelet Toolbox (Jacques et al., 2011). The first step, namely the ridge detection, is dependent on a number of input parameters tabulated in Tables 1 and 2. Equation (A27) gives the expression for the ridge measure.

Figure 1 The automated fracture trace detection workflow.

Table 1 Shearlet system parameters.

Table 2 Detection parameters.

An optimal set of deterministic parameter values which can extract features on all scales is not known a priori. Therefore, we vary the input parameters corresponding to the construction of the shearlet system and the ridge detection parameters within user-defined ranges to compute multiple ridge realizations. A ridge ensemble map is obtained by superposing the ridge images and normalizing. A simple sigmoid function is applied to the normalized ridge ensemble to non-linearly scale and thereby isolate higher image intensities. A user-defined threshold is then applied to the intensity values of this non-linearly scaled, normalized ridge ensemble image to extract a highly probable, binarized ridge network. The threshold is set by a visual comparison of the input image with the extracted ridges. The range for each parameter in Tables 1 and 2 is ascertained by first testing the effect of variation in each parameter with respect to a chosen base case image. This approach to automated detection captures features of multiple scales and highlights regions of uncertain feature extraction within the image.

The second step is the segmentation of the detected ridges using Otsu thresholding (Otsu, 1979). This operation removes small, disconnected, and isolated ridge pixel clusters. The third step is a skeletonization procedure where clusters of pixels representing the segmented ridges are thinned into single pixel representations. For intersecting fractures, the skeletonization procedure preserves the topology of the fracture network by recognizing and splitting the frame at the branch point. This step ensures that in subsequent DFN representation, there is no further effort expended on manually connecting the detected segments.

The fourth step involves piecewise linear polyline fitting to the skeletonized clusters. By default, our code attempts to fit polylines rather than lines to the pixel clusters. Polyline fitting retains geologically realistic veering and curvature of fractures in the vectorized result. We use functions from the Geom2D toolbox (Legland, 2019) for the skeletonizing and polyline fitting. The fifth step is a line simplification procedure applied to the piecewise linear polyline clusters. A large number of polyline points would increase the size of vectorized files; hence, we use the Douglas–Peucker line simplification algorithm (Douglas and Peucker, 1973) implemented by Ahmadzadeh (2017). The algorithm simplifies a piecewise linear polyline into one which has fewer segments. The number of polyline points assigned to each skeletonized cluster is set constant in the code, but this may be modified to be a linear function of the cluster size measured in pixels. If the image is georeferenced or the image scale is known, the code georeferences the simplified polylines and writes to a vectorized shapefile format.

The DFN in the vectorized shapefile format may now be used for any application that requires explicit fracture network geometry. An example of a fractured Posidonia shale micro CT (computed tomography) image slice from Dwarkasing (2016) (see Fig. 2) illustrates the effects of each of the steps involved.

Figure 2 Illustration of the steps involved in the automatic fracture extraction using a 40×34 mm fractured shale core image. (a) CT scan core image from Dwarkasing (2016). (b) Normalized ridge ensemble. (c) Segmentation applied to the ridge ensemble. (d) Skeletonization applied to the segmented ridge. (e) Vectorized polylines fitted to the skeletonized clusters. (f) Effect of line simplification applied to a single vectorized segment.

3.2 Sensitivity analysis of parameters on extraction results

Since the detection results may vary owing to different parameter combinations, we conducted a sensitivity analysis to investigate the ridge extraction output with variation in parameter input. An example of a fractured image sample representing Mesoproterozoic sandstone from the Tomkinson Province, Northern Territory, Australia (Fig. 3a), is chosen to study the effect of shearlet parameter variation. The image dimensions are 1313×1311 pixels, and it has four prominent fractures with two of them forming an intersection. A subtler fracture is present towards the top left and a thick fracture located at the bottom left of the image. A base case set of parameters for constructing a shearlet system and for ridge identification is set up in the table next to Fig. 3a. We vary all parameters one by one with respect to this base case. Ridge extraction using the base case shearlet system shows that the major intersecting fracture system is identified; however, the largest fracture is detected only partially and that, too, only at the peripheries (see Fig. 3b). The subtle fracture is detected but disconnected. A large amount of noise is also present.

Figure 3 Effects of variation in ridge parameters on extracted ridges and the corresponding vectorizations using a fractured siliciclastic example. A constant greyscale threshold is applied to the ridge map and all other parameters with respect to post-processing are kept constant. (a) Fractured rock image (5 cm×5 cm) and tabulated base case parameters. (b) Ridge maps and vectorized traces for base case. (c) Effect of a higher Gaussian effect support compared to wavelet support. (d) Effect of a large difference in wavelet effective support with respect to Gaussian support. (e–o) Lower bounds of parameters with respect to base case, corresponding ridge maps and traces. (f–p) Upper bounds of parameters with respect to base case, corresponding ridge maps, and traces.

The complex shearlet system is constructed by the tensorial product of a Mexican hat wavelet and a Gaussian filter. The first two parameters waveletEffSupp and gaussianEffSup refer to the pixel widths over which the wavelet amplitudes sharply change from zero. The even- and odd-symmetric elements of the constructed shearlet system using the base case parameters for the siliciclastic example are depicted in Figs. 4i and ii. We chose to maintain a ratio of 2 between waveletEffSupp and gaussianEffSup. The effect of increasing the effective support on the complex shearlet system is shown in Fig. 4xvii–xix. Figure 4xx and xxi indicate the effects of large ratios between the wavelet effective support and Gaussian support. The second parameter is the scalesperOctave, which determines the number of intermediate scales per octave. An octave is the interval between two frequency peaks. For example, we may consider a wavelet that is scaled by a factor of 2. Physically, this means a stretching of the wavelet, thereby decreasing the frequency. The base-2 logarithmic ratio of the reduced frequency with respect to the original frequency is the number of octaves by which the frequency has reduced. We set the number of octaves as a constant value of 3.5. This implies that there are seven scales for the complex shearlet system as can be seen in Fig. 4iii–ix. The shearLevel parameter indicates the discrete number of orientations that the complex shearlet system can assume. The selected value of 3 indicates that there are 2³+2 (10) orientations possible for the complex shearlet system (see Fig. 4x–xvi) and 2(2³+2) (or 20 shearlets). For large images and large number of shearlets, computational effort is quite expensive. The alpha parameter is the degree of anisotropy induced by scaling with a null value of alpha maximizing the degree of anisotropy. We vary alpha, shearLevel, and the scalesperOctave, but the effects on the constructed complex shearlet system are minimal as can be seen from Fig. 4xxii–xxx.

Figure 4 Effects of parameter variation on the constructed complex shearlet system for the fractured siliciclastic example. (i) Even-symmetric elements of the complex shearlet system constructed using the base case parameters in Fig. 3. Full system is 1313×1318 pixels. (ii) Odd-symmetric elements of the complex shearlet system using the base case parameters in Fig. 3, (iii)–(ix) depiction of seven scales, (x)–(xvi) depiction of seven orientations (out of possible 10) for the odd-symmetric elements of the complex shearlet system, (xvii)–(xix) effect of wavelet effective support and Gaussian effective support on the even-symmetric elements of the complex shearlet system, (xx) effect of Gaussian effective support double that of wavelet effective support, (xxi) effect of wavelet effective support very large than Gaussian effective support, (xxii)–(xxiv) effect of scales per octave on the even-symmetric elements, (xxv)–(xxvii) effect of anisotropy parameter on the even-symmetric elements, and (xxviii)–(xxx) effect of shear levels on the even-symmetric elements.

The effects of variation in the parameters on ridge extraction are depicted in Fig. 3c–p. Decreasing the value of the support by half identifies finer features, but then the largest fracture is completely missed (Fig. 3e). When the support is doubled, the emphasis on larger features is more pronounced (Fig. 3f). The effects of increasing and decreasing scalesperOctave are depicted in Fig. 3g and h with a higher value resulting in a finer ridge map. The effect of an increase and decrease in the number of shear levels on the final ridge map is quite minimal as can be seen from Fig. 3i and j. The effect of anisotropy parameter alpha is depicted in Fig. 3k and l with minimal anisotropy resulting in a finer ridge map. The minContrast parameter is a greyscale threshold (values from 0 to 255) applied to Eq. (A27) to extract ridges. A larger value suppresses noisy features as can be seen from the comparison between Fig. 3m and n. The offset parameter is a scaling offset between odd-symmetric and even-symmetric shearlets quantified in octaves. Reducing the value of this parameter results in a coarser ridge map with enhanced connectivity (Fig. 3o) compared to the larger value, which results in a finer map (Fig. 3p).

From an interpreter's point of view, three different scales of fracturing need to be identified and false features also need to be suppressed. From the sensitivity analysis, the parameters that are most important to generate high-probability ridge maps are the wavelet supports (required to capture multiple scales of fracture), greyscale contrast (suppressing noise and thereby false features), and even–odd offset (which suppresses ridge detachments). This example illustrates the necessity of computing a ridge ensemble instead of searching for an ideal parameter combination.

3.3 Shearlet parameter selection

To decide upon the shearlet parameter space to generate multiple ridge realizations, we chose one sample image (see Fig. 5a). Base case parameters are chosen based on recommendations underlined in Reisenhofer et al. (2016) for shearlet construction and ridge detection, and these are tabulated in Table 3. The use of these results in the overlay depicted in Fig. 5b. As can be observed from visual inspection of the overlay of the detected ridges over the original image, the automatic method can extract a large number of fractures. However, there are still some false positives (features detected on the trees and inside the large karstic cavities) and false negatives (undetected small-scale fractures).

Table 3 Shearlet system and detection parameters used to extract ridges for the base case.

Figure 5 Effect of multiple ridge realizations on a sample image from Parmelan anticline, France (Prabhakaran et al., 2019b). (a) Base case image used for testing the effect of multiple ridge realizations. (b) Ridge map obtained using the base case shearlet parameters in Table 3. (c) Overlay of ridges obtained using base case shearlet parameters over the test image. (d) Normalized ridge intensity ensemble map obtained after 1050 ridge realizations using the parameters in Table 4. (e) Threshold ridge intensity map that enhances features. (f) Overlay of ridges using the threshold ridge intensities over the test image.

To select the parameter ranges, we vary parameters with respect to the base case ridge image, thereby generating multiple ridge images. We use the structural similarity measure (Wang et al., 2004) to quantify the difference between the base case ridge image and other ridge images. Structural similarity (SSIM) is a measure commonly used in image quality assessment that returns one value as a measure of similarity between two images, where one image is the reference image. The SSIM is calculated for each ridge realization image corresponding to each parameter with respect to the base case ridge image. The SSIM for variation in scaling offset, anisotropy scaling α, Mexican hat wavelet support, Gaussian filter support scales, minimum contrast, scales per octave, and number of shear levels are depicted in Fig. A3 according to the range of parameters in Table 4. From the analysis of the effects of parameters, we decided to vary the shearlet construction parameters so that we have 70 complex shearlet systems (see Table A1 for the parameters used to construct the 70 complex shearlet systems).

Table 4 Ensemble for parameter variation.

The total number of stochastic runs for the ridge detection is the number of combinations of shearlet systems and ridge specification parameters. Using such an approach, a probability map of detected features may be obtained, based on which cut-off thresholds can be defined to remove false positives. The result of such a stochastic run with 1050 realizations is depicted in Fig. 5. From this result, the utility of the method is evident; the features that are obscured by shadows and the shrubbery have a low-strength signal which can then be filtered away, thus reducing the number of false positives. Another advantage is that both large-scale and fine features are captured, which may not be possible using a single set of shearlet parameters.

4 Results

4.1 Trace extraction results from Parmelan, France

4.1.1 Geological setting of the Parmelan plateau

We tested the automated fracture extraction method on an example from a carbonate outcrop from the Parmelan plateau in the Bornes Massif, France. The Bornes Massif is a northern subalpine chain in the western French Alps. The method was applied to a photogrammetric orthomosaic derived from a 3-D outcrop model. The outcrop model was built from source photos acquired using a DJI Phantom 4 UAV. The image resolution is 18.6 mm per pixel. Processing of the drone images and generating the orthomosaic was done using AgiSoft PhotoScan Professional (Version 1.2.6) (2016) software. The Parmelan anticline in France (see Fig. 6) is situated in the frontal part of the Bornes Massif and consists of Upper Jurassic to Cretaceous rocks of the European passive margin (Huggenberger and Wildi, 1991; Gidon, 1996, 1998; Berio et al., 2018).

Figure 6 Location of the Parmelan plateau in France within the Bornes Massif depicting drone flight base points for six drone missions. This map was generated using satellite imagery obtained from ESRI World Imagery (https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer, last access: 12 February 2019) and modified using ArcGIS 10.3, ArcMap 10.3 software by ESRI (http://www.esri.com/, last access: 12 February 2019). Service layer credits: ESRI, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User Community.

This NE–SW trending anticline consists of a wide, flat crestal plateau bounded by steeply dipping limbs. Carbonates form the roof of a kilometre-scale box fold formed during the Alpine orogeny (Bellahsen et al., 2014). On the crestal plateau, a 1.7 km by 2.3 km large pavement of flat-lying shallow-water carbonates is exceptionally well exposed. The Parmelan outcrop is a good example of fracture patterns formed in a fold-and-thrust setting. We applied the automatic fracture detection technique on an orthomosaic that has been stitched together from drone photogrammetry over six different drone missions over the Parmelan. The combined extent of the six orthomosaics is depicted in Fig. 7a, and the areal extent of each orthomosaic is depicted in Fig. 7b.

Figure 7 Drone photogrammetry coverage area from the Parmelan (a) Region within the Parmelan plateau highlighting the areal extent of the drone photogrammetric orthomosaics which are projected over the base map. Manually traced large-scale faults are depicted in red. This map was generated using satellite imagery obtained from ESRI World Imagery (https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer, last access: 12 February 2019) and modified using ArcGIS 10.3 and ArcMap 10.3 software by ESRI (http://www.esri.com/, last access: 12 February 2019). Service layer credits: ESRI, DigitalGlobe, GeoEye, Earthstar Geographics, CNES/Airbus DS, USDA, USGS, AeroGRID, IGN, and the GIS User Community. (b) Spatial extent of the drone coverage of each of the six UAV flight missions in different colours.

4.1.2 Automatic extraction results on the Parmelan orthomosaic

Considering memory requirements and for faster computation, the image domain was divided into georeferenced sub-tiles using the Grid Splitter plugin in QGIS software. Visual filtering was carried out to remove tiles that did not have exposed rock, had a large degree of shrubbery, and which were at the orthomosaic edges where image resolution is poor. A total of 1000 tiles were chosen for the automated interpretation process. The areal extent of the orthomosaic covered 0.589 km², and this region is depicted in Fig. 7. The region covered by the tiles amounts to 0.379 km², and this is shown as an overlay of the selected tiles in Fig. 8a. Structural measurements were recorded at four small-scale stations (around 2–5 m² per station) depicted in Fig. 8c–f.

Figure 8 Trace extraction results from the Parmelan. (a) Selected tiles from the orthomosaic (Prabhakaran et al., 2019b) spatial extent are highlighted. Points where structural measurements were collected (stations) are depicted using red squares. The region where there is a comparison between manual and automatic interpretations is highlighted by the green square. (b) Spatial variation in the fracture intensity depicted as a P₂₁ plot computed using the box counting method. (d) Rose and stereoplot of Station 1. (e) Rose and stereoplot of Station 2. (e) Rose and stereoplot of Station 3. (f) Rose and stereoplot of Station 4. The Parmelan dataset (Prabhakaran et al., 2019b) is available under a CC-BY-NC-SA license.

An ensemble of 1050 ridges was computed using a set of shearlet parameters. A threshold for the ridge intensity was chosen to filter out the false positives. The threshold was determined by a visual examination of the overlay of detected ridges over the original images. The subsequent post-processing steps yielded features in each tile. These were georeferenced and stitched back into a single vectorized file representation. Around 3 million features were extracted from the Parmelan orthomosaic. The P₂₁ fracture intensity was computed using the box-counting method by dividing the tile into a 25×25 (pixels) regular grid. The P₂₁ fracture intensity plot highlights the spatial variation in fracturing over the Parmelan plateau (see Fig. 8b). The vectorized fracture shape files along with the Parmelan basemap are presented as a public dataset (see Prabhakaran et al., 2019a).

4.1.3 Comparison with manual interpretation and structural observations

To compare results of the automated approach to a manual interpretation, we chose a subregion within the Parmelan orthomosaic. The selected subregion depicted in Fig. 9a consists of a 24 m×24 m tile of the Parmelan orthomosaic. The image indicates fractures that seem to be isolated, without a well-connected topology, and which are predominantly aligned along a NW–SE direction. The fracturing intensity is variable across the tile. The contrast between fractures and the host rock fabric is intensified by the karstification of the fractures, which can be attributed to weathering and dissolution. Figure 9b depicts an overlay of the automatically interpreted fractures overlain over the original tile. A total of 2910 features was extracted in this tile. This example highlights some of the technical challenges associated with automated fracture trace detection. Shrubbery is present in the image, which obscures certain relevant features. The north-western corner of the image is blurred since it forms the extent of the orthomosaic.

Figure 9 Comparison between automatic and manual interpretation on a tile from the Parmelan. (a) Tile from the Parmelan orthomosaic (Prabhakaran et al., 2019b) depicting intense fracturing with an organization along the NW–SE corridors. (b) Overlay of fractures traced using the automatic detection method. (c) Overlay of fractures manually traced for the tile at a zoom of 1:2000. (d) P₂₁ Fracture intensity for automatic extracted fractures. (e) P₂₁ Fracture intensity for manually extracted fractures. (f) Length-weighted rose plot and cumulative trace length frequency distribution for the automated result. (g) Length-weighted rose plot and cumulative trace length frequency distribution for the manual result. The Parmelan dataset (Prabhakaran et al., 2019b) is available under a CC-BY-NC-SA license.

The image also depicts open cavities or blobs, which could be the result of localized weathering. The effect of the cavities on the feature extraction is that only an edge is detected. Overall the fracture extraction efficiency is quite dependent on the resolution and quality of images. In the case of the Parmelan data acquisition, the UAV was flown at an altitude of 50–70 m above the pavement; therefore, features such as closed veins and slightly open fractures are below the resolution of the drone camera. A higher image resolution is necessary to extract such features. In our specific case study, good lighting and exposure during the UAV flight mission prevented shadows from obscuring the imagery. Figure 9c depicts a manually performed interpretation at a zoom level of 1:2000 on the raster image with a total of 341 features. P₂₁ fracture intensity comparisons of both automatic and manual traces are shown in Fig. 9d and e. The difference between the automatic and manual interpretation highlights the inclination of the interpreter to neglect small-scale features. Based on geological experience and prior knowledge of the field area, there is a tendency to interpret and link together disconnected features from the original raster image. The closest small-scale station to the sub-tile depicted in Fig. 9a is station 2. There is agreement between the rose plots of station 2 (see Fig. 8c) and the interpretations (Fig. 9e and f). The observed fractures in both cases are predominantly sub-vertical.

4.1.4 Application to mineralized fractures

We now showcase an example of a close-range image containing mineralized veins that are invisible to photogrammetry at altitudes of 40–70 m. The resolution of this image is 0.18 mm per pixel and was taken using a handheld DSLR camera. In this high-resolution image, the fracture infill has similar colour to the host rock as can be seen in Fig. 10a. A manual interpretation of the veins (at a zoom of 1:750) is depicted in Fig. 10b. Using a well-tuned set of parameters with reduced wavelet effective supports, it is possible to extract the much thinner and subtle features, as depicted in Fig. 10c. It can be observed from a comparison between Fig. 10b and c that a large number of false features are also highlighted alongside the features of interest. The main contributors to the extraction of these non-fracture features are the natural rugosity of the rock face, the presence of pebbles, pockmarks, and erosion features. The arrangement of these artefacts displays a very different pattern: small lines with random direction compared to the fractures which are consistently oriented and quite continuous. The veins are also of different thicknesses, with a few veins anastomosing and some branching in a horsetail manner. Some of the thicker veins also exhibit a microstructure within the mineral infill. Further tuning of parameters in order to capture all the veins while also suppressing false features is quite challenging, and hence we do not explore this in further detail. Despite the noise, the automated method is not limited to capturing only open fractures but can also extract closed fractures.

Figure 10 Extension of the automated method to extract mineralized fractures. (a) Image from Parmelan depicting mineralized fractures. (b) Manual interpretation of mineralized fractures. (c) Ridge ensemble.

4.2 Trace extraction results from Brejões, Brazil

4.2.1 Geological setting of the Brejões pavement

The second case study for the automated extraction method is a carbonate outcrop from the Irecê Basin, central Bahia, Brazil (see Fig. 11a, b). The Irecê Basin is located within the northern region of the São Francisco Craton. The Brejões pavement study area is within the Irecê Basin and consists of Neoproterozoic platform carbonates of the Salitre Formation (750–650 Ma). The Neoproterozoic cover was affected by the Brasiliano Orogeny (750–540 Ma) in two separate folding events resulting in fold belts around edges of the São Francisco Craton (Ennes-Silva et al., 2016). The Brejões pavement UAV imagery that we used for our analysis was acquired by Boersma et al. (2019). Structural measurements from Boersma et al. (2019) are shown in Fig. 11c. The orthomosaic covers an area of 0.81 km² and consists of fractured, black oolitic limestones that correspond to Unit A1 of the Salitre stratigraphy (Guimarães et al., 2011). The resolution of the Brejões orthomosaic is 20.3 mm per pixel.

Figure 11 Trace extraction results from the Brejões outcrop. (a) Bahia state in NE Brazil. (b) Location of the Brejões outcrop in the state of Bahia. (c) Structural measurements from Brejões outcrop, adapted from Boersma et al. (2019). (d) Selected tiles from the Brejões orthomosaic (Prabhakaran et al., 2019a) for the automated extraction. (e) Spatial variation in the fracture intensity depicted as a P₂₁ plot computed using the box counting method. Maps depicting administrative boundaries of Brazil and Bahia state was modified from free vector spatial data downloadable at DIVA-GIS (https://www.diva-gis.org/Data, last access: 14 April 2019). The Brejões dataset (Prabhakaran et al., 2019a) is available under a CC-BY-NC-SA license.

4.2.2 Automatic extraction results on the Brejões orthomosaic

The Brejões orthomosaic is split into 222 tiles for the analysis, and this region is shown in Fig. 11d. The Brejões example has a different fracturing style than the Parmelan and consists of an intricate pattern of multi-scale conjugate fractures. The shearlet combinations utilized in the case of the Parmelan were insufficient to capture this variation in scales. Specifically, in the Brejões case, the large-scale features were not captured. A visual inspection of the ridges was necessary to identify the shearlet combinations that amplified the large-scale features. The contribution of these ridges was increased (factor of 8) in the ridge ensemble to highlight these large deformation features. Figure 11e depicts the P₂₁ fracturing intensity computed using the box-counting method by dividing each tile into a 25×25 (pixel) regular grid. The vectorized fracture shape files along with the Brejões basemap are presented as a public dataset (see Prabhakaran et al., 2019b).

4.2.3 Comparison with manual interpretation and structural observations

The automatically extracted features from the Brejões image data were compared with manual interpretations performed by and obtained from Boersma et al. (2019) at seven stations. The automatic interpretations were trimmed to the peripheries of the manual interpretations for a fair comparison between both vectorizations. The location of these stations alongside the automatic versus manual interpretations is shown in Fig. 12. A comparison of the rose plots and cumulative length distributions of the manual and automatic interpretations is depicted in Fig. 13. A few observations can be made from the comparison. Firstly, similar to the Parmelan case, the interpreter picks a lower number of features. Secondly, there is a tendency to extend fractures across image regions where there is no real evidence of rock failure. Thirdly, there is an inconsistency in specifying the connecting topologies between the interpreted traces.

Figure 12 Comparison between manual (left) and automatic (right) interpretation on seven stations within the Brejões outcrop (Prabhakaran et al., 2019a). The manual interpretations were obtained from Boersma et al. (2019). The Brejões dataset (Prabhakaran et al., 2019a) is available under a CC-BY-NC-SA license.

In some stations (see Mid no. 2, Mid no. 3, and North in Fig. 12), the automated interpretation suffers from a large number of false positives. A close examination indicates that the presence of shadows and eroded, undulating topography of the rocks is the main reasons for these false positives. In the Brejões case, the drone was flown at around 10:00, and hence the exposure of the outcrop face was not optimal. The inclined illumination enhances shadows on the rugged topography, which are then seen as false positives in the automatic interpretation. False positives due to shrubbery are minimal in the station regions considered.

Figure 13 Comparison of trace length-weighted rose plots and cumulative trace length distributions for automatic and manual trace interpretations from Brejões outcrop stations. The Brejões dataset (Prabhakaran et al., 2019a) is available under a CC-BY-NC-SA license.

4.3 Benchmarking with data from Thiele et al. (2017)

We further tested the automated trace detection on a recently published case study from Thiele et al. (2017a). The images selected are orthophotographs of two 10 m×10 m areas from Bingie Bingie Point, New South Wales, Australia (see Figs. 14a and 15a). The exposed rocks are Cretaceous to Paleogene dykes, intruding diorites, and tonalities cross-cut by joint sets (Thiele et al., 2017a). The images are complex as they contain both open and closed fractures of different scales, distributed between multiple lithological layers. The images also contain water, shadows, and debris, which makes it even more challenging. We chose this dataset to benchmark the quality of our results with those presented using the semi-automatic cost-function-based trace mapping approach of Thiele et al. (2017a).

Figure 14 Comparison of benchmark image 1. (a) Bingie Bingie area 1 from Thiele et al. (2017a). This image is available for download at https://doi.org/10.4225/03/5981b31091af9 (Thiele et al., 2017b) under a CC-BY-4.0 license. (b) Normalized ridge map using complex shearlet automatic extraction. (c) Threshold applied to the normalized ridges. (d) Binarized ridges map. (e) Vectorized automatic traces. (f) Length-weighted rose plot of the automatic interpreted traces. (g) Cumulative frequency plot of the automatic interpreted traces. (h) Assisted cleaned-up trace map for area 1. (i) Length-weighted rose plot of the assisted traces. (j) Cumulative frequency plot of the assisted traces. (k) Assisted trace interpretation modified from Thiele et al. (2017a) for comparison.

The variation in fracture scales implied that similar to Brejões, a different set of shearlet combinations was needed. We generated 2700 ridge realizations which were used to construct a normalized ridge ensemble map for both images (see Figs. 14b and 15b). A simple, non-linear sigmoid function was applied to the normalized ridge intensity to enhance ridge strength (see Figs. 14c and 15c) and a threshold was chosen based on visual comparison with the source image to yield highly probable, binarized ridge images (see Figs. 14d and 15.d). The subsequent workflow steps, as described in Sect. 3.1, were followed to obtain vectorized traces (see Figs. 14e and 15e). The vectorized traces were used to render assisted interpretations depicted in Figs. 14f and 15f which are comparable in quality to the assisted interpretation of Thiele et al. (2017a).

Figure 15 Comparison of benchmark image 2 (a) Bingie Bingie area 2 from Thiele et al. (2017a). This image is available for download at https://doi.org/10.4225/03/5981b31091af9 (Thiele et al., 2017b) under a CC-BY-4.0 license. (b) Normalized ridge map using complex shearlet automatic extraction. (c) Threshold applied to the normalized ridges. (d) Binarized ridges map. (e) Vectorized automatic traces. (f) Length-weighted rose plot of the automatic interpreted traces. (g) Cumulative frequency plot of the automatic interpreted traces. (h) Assisted cleaned-up trace map for area 2. (i) Length-weighted rose plot of the assisted traces. (j) Cumulative frequency plot of the assisted traces. (k) Assisted trace interpretation modified from Thiele et al. (2017a) for comparison.

In the published results of Thiele et al. (2017a), assisted interpretations of both areas are achieved in 37 and 34 min. We can report better performances of 27 and 32 min for the same areas. The time does not include computing of the ridge realizations. Once the high-probability trace map was generated, the subsequent steps of the automated detection workflow took around 3 min. The remainder of the time was used to perfect the assisted interpretation. The post-processing tasks performed in this second step were the removal of false positives owing to shadows, water, and debris and the joining of segments which were disjointed due to poor resolution within the image. Although we have performed a benchmarking exercise with the data from Thiele et al. (2017a) and also compared our results with manual interpretation, it would be useful to compare them with more manual interpretations to further validate the accuracy of the technique. Such comparison, however, can be done only on networks which are either limited in their spatial extent or in the number of features interpreted. For large orthomosaics, a benchmarking exercise can be challenging as few manually rendered datasets are comparable in network size.

5 Discussion

The extraction of fracture traces from photogrammetric data is a necessary processing step to construct DFN representations. DFNs created using fracture patterns that are directly extracted from rock images are advantageous as they honour the spatial architecture of fracture networks. Automated extraction methods reduce the human component in data processing, and we have achieved this using the complex shearlet transform ridge detection method accompanied by post-processing steps. The complex shearlet method can detect both edges as well as ridges in fractured rock images. We find that the ridge measure works very well for the extraction of fractures, and we use the ridge measure in all our case studies. Although the method performs very well and can extract much more traces than is possible manually while reducing interpreter bias, there are some issues that need to be mentioned. In this section, we discuss the validity and limitations of the technique, areas where there is scope for further development, and also describe some potential extended applications of the method.

5.1 Validity and limitations

• Detection of mineralized features. The method works well when the features of interest are barren and prominent. When fractures are closed and filled, then they are generally harder to detect and require high-resolution images (<1 mm per pixel), which can be recorded only at very close ranges at very low UAV flight altitudes. Recent outcrop studies (Ukar et al., 2019) indicate that many of the barren features in outcrop are absent within the same subsurface lithological unit while maintaining a good correspondence between mineralized features in both outcrop and subsurface. When mineral fill has a marked colour contrast with respect to the host rock (as in vein data published recently by Meng et al., 2019), then superpixel segmentation algorithms can be successful (Vasuki et al., 2017). In the case of poor contrast, the complex shearlet transform would require a great deal of manual tuning of detection parameters to extract reliable results. At such close ranges, as is needed for veins extraction, it is also likely that many more noisy features unrelated to fracturing would arise.

  Since the mineral fill of fractures can provide a clearer picture of the evolution, timing, and stress history of fractures, identifying them on an outcrop scale is important. This is doubly significant, when the goal is to directly extrapolate fracture patterns from a particular outcropping to the same subsurface target. In such a case, close range UAV-mounted hyperspectral data acquisition would be better suited than conventional imaging and image processing methods. With hyperspectral imaging, data are collected in near-continuous spectral bands. The spectral response of minerals constituting the rock is owed to atomic molecular-level processes triggered on interaction with a light source (active or passive) and this may be utilized to identify mineral composition. Since the mineral fill of veins is likely to have a different spectral response from the mineralogy of the host rock, this variation may be used to isolate the pixels that correspond to veins.

  A recent review on close range hyperspectral imaging for mineral identification identifies various previous studies performed for specific minerals (Krupnik and Khan, 2019). It would be interesting to observe, identify, and distinguish between mineralized sequences based on the differences in spectral response of the fracture infill material. Since hyperspectral data are much more voluminous and with significantly more complex image processing than conventional photogrammetry, such an analysis could be confined to selected regions within the outcrop. In conjunction with conventional UAV photogrammetry that covers a larger spatial area, laboratory-based geochemical studies, and outcrop observations (scanline sampling, abutting relations, etc.), a more detailed fracture characterization may be conducted.

• Detection of large cavities and false features. Both the Parmelan and Brejões pavements exhibit karstification, with the Parmelan containing many more collapsed karstic regions. The presence of such low-aspect ratio discontinuities is quite rare in siliciclastic and volcanic outcrops but can prove problematic in the application of the method in carbonate outcrops where karstification is severe. Both the ridge and edge measures would fail in identifying such blobs or would, at best, extract the periphery of the cavity. In recent work by Reisenhofer and King (2019), blob detection measures have been developed within the shearlet framework and could potentially solve this issue.

  Another issue is the effect of undulating topography and shrubbery in generating false positives. False positives generally appear when there is shrubbery, shadows, very rugged terrain, and non-fracture bedding planes. In the case of the Parmelan, the use of multiple ridges was successful in suppressing the false positives owing to shrubbery. However, in Brejões, false positives due to underbrush were more difficult to suppress because they shared the same scale as that of the fractures. In Brejões, shrubbery was also present within some of the wider fractures causing false negatives. In such cases, manual interference is necessary to either mask the regions of shrubbery before the automatic extraction or to remove (or connect) the vectorized traces after the automated extraction. Additionally, carbonate outcrops are prone to widespread erosion owing to exposure to meteoric water from precipitation cycles and air corrosion. Geomorphological features owing to these erosive processes may also play a role in the generation of false positives.

• Parameter selection. A significant difference in fracture scales within an image of interest can prove problematic for the method. In such a case, a vast number of ridge detection runs and an associated increase in computational time is needed to construct a ridge ensemble that takes into account all scales of discontinuities and yields a satisfactory result. When such variation is localized, the image could be segmented into regions that correspond to varying fracture intensities and processed separately. This may be difficult to assess a priori and in such cases would require trial runs. In the Brejões outcrop example and the close range Parmelan vein example, this difference in fracture scales was ubiquitous throughout the exposure and more pronounced than the Parmelan outcrop. Using a visual comparison with the original image, the effect of ridges resulting from certain shearlet parameter combinations was enhanced, so that the ridge ensemble is improved. In Brejões, it was the large-scale features that needed to be strengthened, while in the case of the Parmelan vein example, the smaller features needed sharpening. Since parameter selection is still done manually, a more comprehensive way of arriving at the optimal shearlet combination is desirable. An algorithm that automatically optimizes for shearlet parameters corresponding to each individual scale of fracture is worthy of attention.

• Artificial fragmentation of traces. Manual fracture interpretation from images often involves the step of classifying fracture traces into separate sets based on ground truth observations or with respect to fracture strike. The automated method described here in its current form can only extract traces and cannot distinguish between or classify traces as belonging to separate sets. When fractures intersect each other, the issue of artificial fragmentation of seemingly continuous traces arises. If an image consists of two orthogonally intersecting fractures, the automated method would result in four traces intersecting at a single branch point, even though a manual interpretation would only identify two fracture traces belonging to two different geometric sets. This type of fragmentation would result in different length distributions as observed in (insert figures that compare automatic and manual interpretations). This kind of fragmentation is not an issue if such an outcrop DFN is used for geometric input for a flow or geomechanics simulation. This is because the process of meshing models with explicitly specified DFN geometry would, in any case, require the specification of all intersection points (or forced fragmentation of long intersecting fractures). Therefore, the practitioner must exercise caution when using cumulative length distributions derived from outcrop DFNs that are automatically extracted.

  A single fracture could also be fragmented without being cut by other intersecting fractures. This may happen in the case of false negatives (due to shadows falling over part of a fracture, debris, or shrubbery within an open fracture and when fracture opening is very thin at some regions along the fracture length) that cause fragmentation of fractures with gaps in between them. This kind of fragmentation affects the topology of the network in addition to depressing the height cumulative length distribution. It maybe noted that the specification of fracture endpoints manually is also fraught with bias (Peacock et al., 2019). A solution would be to use a range of linking thresholds to connect traces and study the effects of threshold values on network topology and length distribution.

5.2 Recommendations for future work

• Link between extractable P₂₁, drone flying altitude, and camera resolution. From the P₂₁ analysis on the Parmelan and the Brejões automatically extracted fractures, the maximum value P₂₁ was around 8 m⁻¹. The same drone model was used in both cases (DJI Phantom 4), and the flying altitude was also similar (between 40 and 70 m). Although such a conjecture needs further verification, there could be a relation between the resolution of imagery and maximum extractable fracture intensity. Often flight altitudes are chosen by drone pilots depending upon considerations such as local topography, weather conditions, and the presence of impediments (such as trees, electricity poles, and telecommunication towers). A detailed analysis of the relation between flying altitude (and consequently image resolution) and extracted fracture intensity could provide drone pilots with insights and guidelines for UAV-based outcrop analysis. The ideal flying resolution to identify features of interest may be ascertained by carrying out a series of acquisitions at a location where ground truth is known.

• Generating data for fractured reservoir modelling workflows. Fractured reservoir characterization workflows in the oil and gas industry have traditionally used stochastic techniques that attempt to extrapolate averaged fracture statistics (either from borehole imagery, core data, or outcrop analysis) to reservoir volumes. The use of multiple point statistics (MPS) for fracture network generation was highlighted by Bruna et al. (2019) as an alternative approach to DFN modelling. MPS uses training images of realistic fracture networks to learn patterns and then generate non-stationary fractured reservoir models. Corrected for false positives and noise, the automated method can produce accurate, geologically realistic, and unbiased training images that can feed into the MPS workflow. Since our method can extract large-scale fracture networks (millions of features from sub-square kilometre regions), it is also well suited to providing training data for deep-learning workflows. Recently, the use of generative adversarial networks (GANs) for geological modelling at the reservoir scale was proposed by Dupont et al. (2018), Zhang et al. (2019) as an alternative to conventional geostatistics, MPS, and object-based modelling. GANs form a subset of deep-learning architectures that is used for generative modelling (Goodfellow et al., 2014). GANs that are trained on realistic data can then generate geologically realistic, non-stationary models.

• Deep-learning methods for trace extraction. Deep-learning methods have revolutionized computer vision applications. Various neural architectures have documented high degrees of accuracy in machine vision tasks such as the overall image classification, identification and classification of objects within an image images, localization of objects, extraction of regions of interest (semantic segmentation), and extraction of regions corresponding to individual objects (instance segmentation). The problem of fracture trace extraction falls within the problem category of region extraction of individual objects and hence may be attempted using techniques such as mask region-based convolutional neural networks (R-CNNs) (He et al., 2017). Deep-learning methods, in general, require large amounts of labelled data to train. In the case of a mask R-CNN, the library of training images must contain marked regions (or overlays) indicating pixels of interest that correspond to individual fractures. The automated method described in this paper can be used to rapidly generate a large number of overlay images that can be used as training data for mask R-CNN architectures.

6 Conclusions

This paper presents a method to automatically detect and digitize fracture traces from images of rock fractures using the complex shearlet transform. The technique replaces the task of manually interpreting fractures, which is time-consuming, prone to interpreter bias, and suffers from a lack of repeatability. The case studies that are presented highlight the utility of the complex shearlet-based measure for automatically detecting fracture traces from 2-D images. The automatic trace detection method combines the complex shearlet ridge measure with a series of post-processing steps that include image segmentation, skeletonization, polyline fitting, and polyline simplification. We tested the method at different scales of rock displacement, at outcrop scale (∼10² m), and at station scale (<10 m), using two orthomosaics reconstructed from drone photogrammetry and two rock pavement images. We have considered carbonate and igneous rock lithologies in the case studies. Using the method, we have extracted millions of 2-D features from outcrop-scale drone orthophotos. The processing time of the technique depends upon the intensity of fracturing and the complexity of the fracture networks contained within the image. The automatic trace extraction results are quantitatively compared with manually interpreted fractures on selected subsamples of the image domain using fracture trace density metrics. The automated technique is capable of extracting a much larger number of features, with a marked reduction in bias. The method outlined in this paper greatly simplifies the process of generating deterministic, outcrop-based DFNs. The automatically extracted fracture patterns can be used by structural geologists to link deformation features to tectonic history and by geo-modellers in subsurface NFR modelling.

Appendix A: Overview of the complex shearlet transform

A1 The continuous shearlet system

A shearlet-generating function consists of an anisotropic scaling matrix and a shear matrix. Let the shearlet-generating function be

$$\begin{array}{}\text{(A1)}& \mathit{\psi }\in L^{\mathrm{2}}\left({\mathbb{R}}^{\mathrm{2}}\right).\end{array}$$

The admissibility criterion for the shearlet-generating function is

$$\begin{array}{}\text{(A2)}& \underset{R^{\mathrm{2}}}{\int }{{\displaystyle \frac{\left|\widehat{\mathit{\psi }}\left({\mathit{\xi }}_{\mathrm{1}}{\mathit{\xi }}_{\mathrm{2}}\right)\right|}{{{\mathit{\xi }}_{\mathrm{1}}}^{\mathrm{2}}}}}^{\mathrm{2}}\mathrm{d}{\mathit{\xi }}_{\mathrm{2}}\mathrm{d}{\mathit{\xi }}_{\mathrm{1}}<\mathrm{\infty },\end{array}$$

where $\widehat{\mathit{\psi }}$ is the 2-D Fourier transform of ψ.

A shearlet satisfying Eq. (A2) is an admissible shearlet or a continuous shearlet (Kutyniok and Labate, 2012). The admissibility condition implies that a reconstruction formula exists for the associated continuous shearlet transform. In order to achieve an optimally sparse approximation of an image that possesses anisotropic singularities, the analysing elements must consist of waveforms that range over several scales, orientations, and locations with the ability to become very elongated. To this end, a combination of a scaling operator to generate elements at different scales, an orthogonal operator to change orientations, and a translation operator to displace elements over the 2-D plane is used. The scaling matrix A_a is defined as (Labate et al., 2005)

$${\mathbf{A}}_a=\left(\begin{array}{cc}a& \mathrm{0}\\ \mathrm{0}& a^{\mathit{\alpha }}\end{array}\right),\mathit{\alpha }\in \left[\mathrm{0},\mathrm{1}\right].$$

The value of α controls the degree of anisotropy. (For more information on the anisotropy scaling molecules or α molecules, see Grohs et al., 2016.) The scaling matrix is parabolic when $\mathit{\alpha }=\frac{\mathrm{1}}{\mathrm{2}}$.

An orthogonal transformation is applied to change the orientations of waveforms. Rotation operators are not preferred as they destroy the structure of the integer lattice Z² whenever the rotation angle is different from $\mathrm{0},\pm \frac{\mathit{\pi }}{\mathrm{2}},\pm \mathit{\pi },\pm \frac{\mathrm{3}\mathit{\pi }}{\mathrm{2}}$. Changes in the structure of integer lattice are problematic when transitioning from continuum to a digital setting. Hence, a shearing transformation is used where the anisotropic shearing transformation matrix S_s is defined as

$${\mathbf{S}}_s=\left(\begin{array}{cc}\mathrm{1}& s\\ \mathrm{0}& \mathrm{1}\end{array}\right)\text{, where the parameters }a\in {\mathbb{R}}^{+},s\in \mathbb{R}.$$

The shearing matrix S_s preserves the structure of the integer grid for any s ∈ℕ. The shearing matrix parameterizes orientations using the variable s associated with slopes rather than angles and leaves the integer lattice invariant, provided s is an integer. The difference between isotropic and anisotropic dilation with shearing is depicted in Fig. A1a and b.

Figure A1 (a) Isotropic elements capturing a discontinuity curve (b). Sheared, anisotropic elements capturing a discontinuity curve (modified from Kutyniok and Labate, 2012).

A shearlet system is defined as (Kutyniok and Labate, 2012)

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}\mathrm{SH}\left(\mathit{\psi }\right)=& \left\{{\mathit{\psi }}_{a,s,t}=a^{{}^{-\mathrm{3}}/{}_{\mathrm{4}}}\mathit{\psi }\left({{\mathbf{A}}_a}^{-\mathrm{1}}{{\mathbf{S}}_s}^{-\mathrm{1}}\left(\cdot -t\right)\right)\right.\\ & \left.a\in {\mathbb{R}}^{+},s\in \mathbb{R},t\in {\mathbb{R}}^{\mathrm{2}}\right\},\end{array}\end{array}$$

where $(\cdot -t)$ denotes the translation by a point t.

The corresponding shearlet transform for mapping a function f ∈ L²(ℝ²) into coefficients, SH${}_{\mathit{\psi }}f\left(a,s,t\right)$ specified by scaling a, shearing s, and translation t is given by

$$\begin{array}{}\text{(A4)}& f\to {\mathrm{SH}}_{\mathit{\psi }}f\left(a,s,t\right)=f,{\mathit{\psi }}_{a,s,t}.\end{array}$$

A2 Cone-adapted continuous shearlet systems

Equation (A4) renders horizontal shearlets elongated at very fine scales, which is problematic in digital implementations. Because the shearing operator can range over a non-bounded interval, directions are not treated uniformly. To overcome this drawback of shearing, the cone-adapted shearlet system was introduced in which the frequency plane is split into a horizontal and vertical cone that restricts the shear parameter to bounded intervals (see Fig. A2a). Dividing the frequency plane in such a manner ensures uniform treatment of directions (Guo et al., 2005; Kutyniok and Labate, 2012). A cone-adapted shearlet system can be tiled by further division of the frequency domain. Such a tiling configuration (see Fig. A2b) ensures that all directions are treated “almost equally” (Kutyniok and Labate, 2012). There is still small, but controllable bias in the coordinate axes directions. The cone-adapted shearlet systems can therefore be expressed as the union of a horizontal cone, a vertical cone, and a low-frequency centre component. The frequency plane is thus split into four horizontal and vertical cones with a low-frequency square region in the centre. The low-frequency region is given by the relation (Kutyniok and Labate, 2012)

$$\begin{array}{}\text{(A5)}& \mathcal{R}=\left\{\left({\mathit{\xi }}_{\mathrm{1}},{\mathit{\xi }}_{\mathrm{1}}\right):\left|{\mathit{\xi }}_{\mathrm{1}}\right|,\left|{\mathit{\xi }}_{\mathrm{2}}\right|\le \mathrm{1}\right\}.\end{array}$$

Figure A2 The cone-adapted continuous shearlet system. (a) Bias in directions is handled by dividing the frequency plane into four cones 𝒞₁, 𝒞₂, 𝒞₃, 𝒞₄, and a square low-frequency box region in the centre  ℛ. (b) Trapezoidal-shaped wedge tiling of the frequency-induced domain induced by the shearlet transform (modified after Kutyniok and Labate, 2012).

Figure A3 Variation in the structural similarity index of base case ridges with shearlet and detection parameters. (a) SSIM vs. anisotropy exponent. (b) SSIM vs. minContrast. (c) SSIM vs. wavelet effective support. (d) SSIM vs. scaling offset. (e) SSIM vs. scales or octave. (f) SSIM vs. shear levels.

Inside each cone, the shearing variable s is only allowed to vary over a finite range. This produces elements with uniformly distributed orientations. The union of the generating functions for the horizontal cones ψ ∈ L²(ℝ²), vertical cones $\stackrel{\mathrm{̃}}{\mathit{\psi }}\in L^{\mathrm{2}}\left({\mathbb{R}}^{\mathrm{2}}\right)$, and for the square low-frequency region φ ∈ L²(ℝ²) is expressed as (Kutyniok and Labate, 2012)

$$\begin{array}{}\text{(A6)}& \mathrm{SH}(\mathit{\phi },\mathit{\psi },\stackrel{\mathrm{̃}}{\mathit{\psi }})=\mathrm{\Phi }\left(\mathit{\phi }\right)\cup \mathrm{\Psi }\left(\mathit{\psi }\right)\cup \stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\stackrel{\mathrm{̃}}{\mathit{\psi }}\right),\end{array}$$

where

$$\begin{array}{}\text{(A7)}& {\displaystyle }\mathrm{\Phi }\left(\mathit{\phi }\right)=\left\{\mathit{\phi }{}_t=\mathit{\phi }\left(\cdot -t\right):t\in {\mathbb{R}}^{\mathrm{2}}\right\},\text{(A8)}& {\displaystyle }\begin{array}{rl}& \mathrm{\Psi }\left(\mathit{\psi }\right)=\left\{{\mathit{\psi }}_{a,s,t}=a^{-\frac{\mathrm{3}}{\mathrm{4}}}\mathit{\psi }\left({{\mathbf{A}}_a}^{-\mathrm{1}}{{\mathbf{S}}_s}^{-\mathrm{1}}\left(\cdot -t\right)\right)\right.\\ & \left.:a\in \left(\mathrm{0},\mathrm{1}\right],\left|s\right|\le \mathrm{1}+a^{\frac{\mathrm{1}}{\mathrm{2}}},t\in {\mathbb{R}}^{\mathrm{2}}\right\},\end{array}\text{(A9)}& {\displaystyle }\begin{array}{rl}& \stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\mathit{\psi }\right)=\mathit{\{}{\stackrel{\mathrm{̃}}{\mathit{\psi }}}_{a,s,t}=a^{-\frac{\mathrm{3}}{\mathrm{4}}}\stackrel{\mathrm{̃}}{\mathit{\psi }}({{\stackrel{\mathrm{̃}}{\mathbf{A}}}_a}^{-\mathrm{1}}{{\mathbf{S}}_s}^{-\mathrm{1}}(\cdot -t))\\ & :a\in \left(\mathrm{0},\mathrm{1}\right],\left|s\right|\le \mathrm{1}+a^{\frac{\mathrm{1}}{\mathrm{2}}},t\in {\mathbb{R}}^{\mathrm{2}}\mathit{\}}.\end{array}\end{array}$$

The scaling matrix for vertical cones, ${\stackrel{\mathrm{̃}}{\mathbf{A}}}_a$ is expressed as

$$\begin{array}{}\text{(A10)}& {\stackrel{\mathrm{̃}}{\mathbf{A}}}_a=\left(\begin{array}{cc}{a}^{\mathit{\alpha }}& \mathrm{0}\\ \mathrm{0}& a\end{array}\right).\end{array}$$

The cone-adapted continuous shearlet transform is expressed as the mapping

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}& f\to {\mathrm{SH}}_{\mathit{\phi },\mathit{\psi },\stackrel{\mathrm{̃}}{\mathit{\psi }}}f\left(t^{\prime },(a,s,t\right),\left(\stackrel{\mathrm{̃}}{a},\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t}\right))=\\ & \left(f,{\mathit{\phi }}_{t^{\prime }},f,{\mathit{\psi }}_{a,s,t},f,{\stackrel{\mathrm{̃}}{\mathit{\psi }}}_{\stackrel{\mathrm{̃}}{a},\stackrel{\mathrm{̃}}{s},\stackrel{\mathrm{̃}}{t}}\right).\end{array}\end{array}$$

A3 The discrete cone-adapted shearlet system

A discrete version of the cone-adapted shearlet system may be defined with the scaling parameter j, shearing parameter k, and translation parameter m for a sampling factor of $c=\left(c_{\mathrm{1}},c_{\mathrm{2}}\right)\in {\left({\mathbb{R}}_{+}\right)}^{\mathrm{2}}$. Similar to Eq. (A6) this is a union of the generating functions for the vertical, horizontal, and low-frequency central region.

$$\begin{array}{}\text{(A12)}& {\displaystyle }\mathrm{SH}(\mathit{\phi },\mathit{\psi },\stackrel{\mathrm{̃}}{\mathit{\psi }};c)=\mathrm{\Phi }(\mathit{\phi };c_{\mathrm{1}})\cup \mathrm{\Psi }(\mathit{\psi };c)\cup \stackrel{\mathrm{̃}}{\mathrm{\Psi }}(\stackrel{\mathrm{̃}}{\mathit{\psi }};c)\text{(A13)}& {\displaystyle }\mathrm{\Phi }\left(\mathit{\phi };c_{\mathrm{1}}\right)=\left\{\mathit{\phi }{}_m=\mathit{\phi }\left(\cdot -c_{\mathrm{1}}m\right):m\mathit{ϵ}{\mathbb{Z}}^{\mathrm{2}}\right\},\text{(A14)}& {\displaystyle }\begin{array}{rl}& \mathrm{\Psi }\left(\mathit{\psi };c\right)=\left\{{\mathit{\psi }}_{j,k,m}={\mathrm{2}}^{\frac{\mathrm{3}}{\mathrm{4}}j}\mathit{\psi }\left(S_k{A}_{{\mathrm{2}}^j}\cdot -M_{\mathrm{c}}m\right):\right.\\ & \left.j\ge \mathrm{0},\left|k\right|\le ⌈{\mathrm{2}}^{\frac{j}{\mathrm{2}}}⌉,m\in {\mathbb{Z}}^{\mathrm{2}}\right\},\end{array}\text{(A15)}& {\displaystyle }\begin{array}{rl}& \stackrel{\mathrm{̃}}{\mathrm{\Psi }}\left(\stackrel{\mathrm{̃}}{\mathit{\psi }};c\right)=\left\{{\stackrel{\mathrm{̃}}{\mathit{\psi }}}_{j,k,m}={\mathrm{2}}^{\frac{\mathrm{3}}{\mathrm{4}}j}\stackrel{\mathrm{̃}}{\mathit{\psi }}\left({S_k}^T{\stackrel{\mathrm{̃}}{A}}_{{\mathrm{2}}^j}\cdot -{\stackrel{\mathrm{̃}}{\mathbf{M}}}_{\mathrm{c}}m\right):\right.\\ & \left.j\ge \mathrm{0},\left|k\right|\le ⌈{\mathrm{2}}^{\frac{j}{\mathrm{2}}}⌉,m\in {\mathbb{Z}}^{\mathrm{2}}\right\},\end{array}\end{array}$$

with

$${\mathbf{M}}_{\mathrm{c}}=\left[\begin{array}{cc}{c}_{\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& c_{\mathrm{2}}\end{array}\right];{\stackrel{\mathrm{̃}}{\mathbf{M}}}_{\mathrm{c}}=\left[\begin{array}{cc}{c}_{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& c_{\mathrm{1}}\end{array}\right],$$

(M_c and ${\stackrel{\mathrm{̃}}{\mathbf{M}}}_{\mathrm{c}}$ are sampling matrices for horizontal and vertical cones)

$${\mathbf{A}}_{{\mathrm{2}}^j}=\left[\begin{array}{cc}{\mathrm{2}}^j& \mathrm{0}\\ \mathrm{0}& {\mathrm{2}}^{j/\mathrm{2}}\end{array}\right];{\stackrel{\mathrm{̃}}{\mathbf{A}}}_{{\mathrm{2}}^j}=\left[\begin{array}{cc}{\mathrm{2}}^{j/\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& {\mathrm{2}}^j\end{array}\right],$$

(${\mathbf{A}}_{{\mathrm{2}}^j}$ and ${\stackrel{\mathrm{̃}}{\mathbf{A}}}_{{\mathrm{2}}^j}$ are dyadic scaling matrices for horizontal and vertical cones)

$${\mathbf{S}}_{\mathbf{k}}=\left[\begin{array}{cc}\mathrm{1}& k\\ \mathrm{0}& \mathrm{1}\end{array}\right]\text{ (shearing matrix)}.$$

The discrete cone-adapted shearlet transform associated with ϕ, ψ, and $\stackrel{\mathrm{̃}}{\mathit{\psi }}$ is given by the mapping

$$\begin{array}{}\text{(A16)}& \begin{array}{rl}& f\to {\mathrm{SH}}_{\mathit{\phi },\mathit{\psi },\stackrel{\mathrm{̃}}{\mathit{\psi }}}f\left(m^{\prime },(j,k,m\right),(\stackrel{\mathrm{̃}}{j},\stackrel{\mathrm{̃}}{k},\stackrel{\mathrm{̃}}{m})=\\ & \left(f,{\mathit{\phi }}_{m^{\prime }},f,{\mathit{\psi }}_{j,k,m},f,{\stackrel{\mathrm{̃}}{\mathit{\psi }}}_{\stackrel{\mathrm{̃}}{j},\stackrel{\mathrm{̃}}{k},\stackrel{\mathrm{̃}}{m}}\right).\end{array}\end{array}$$

A4 The complex discrete cone-adapted shearlet system

Taking the complex valued wavelet of a real valued even-symmetric wavelet generator ψ^even ∈ L²(ℝ²) and using the Hilbert transform operator (ℋ), a complex valued shearlet generator is obtained (from Reisenhofer, 2014; King et al., 2015):

$$\begin{array}{}\text{(A17)}& {\mathit{\psi }}^{\mathrm{c}}={\mathit{\psi }}^{\mathrm{even}}+i{\mathit{\psi }}^{\mathrm{odd}}.\end{array}$$

The complex valued function can be written in terms of a Hilbert transform pair of an even-symmetric real valued shearlet and an odd-symmetric real valued shearlet (from Reisenhofer, 2014; King et al., 2015):

$$\begin{array}{}\text{(A18)}& {\mathit{\psi }}^{\mathrm{c}}={\mathit{\psi }}^{\mathrm{even}}+i\mathcal{H}{\mathit{\psi }}^{\mathrm{even}}.\end{array}$$

The Hilbert transform operator is written as

$$\begin{array}{}\text{(A19)}& \mathcal{H}\left(f\right)\left(t\right)=\underset{a\to \mathrm{\infty }}{lim}\underset{-a}{\stackrel{a}{\int }}{\displaystyle \frac{f\left(\mathit{\tau }\right)}{t-\mathit{\tau }}}\mathrm{d}\mathit{\tau }.\end{array}$$

The discrete cone-adapted complex shearlet system is given as (Reisenhofer, 2014; King et al., 2015)

$$\begin{array}{}\text{(A20)}& \mathrm{SH}(\mathit{\phi },\mathit{\psi },\stackrel{\mathrm{̃}}{\mathit{\psi }};c)=\mathrm{\Phi }(\mathit{\phi };c_{\mathrm{1}})\cup \mathrm{\Psi }(\mathit{\psi };c)\cup \stackrel{\mathrm{̃}}{\mathrm{\Psi }}(\stackrel{\mathrm{̃}}{\mathit{\psi }};c)\end{array}$$

and

$$\begin{array}{}\text{(A21)}& {\mathrm{SH}}^{\mathrm{c}}(\mathit{\phi },\mathit{\psi },\stackrel{\mathrm{̃}}{\mathit{\psi }};c)=\mathrm{\Phi }(\mathit{\phi };c_{\mathrm{1}})\cup {\mathrm{\Psi }}^{\mathrm{c}}(\mathit{\psi };c)\cup \stackrel{\mathrm{̃}}{{\mathrm{\Psi }}^{\mathrm{c}}}(\stackrel{\mathrm{̃}}{\mathit{\psi }};c),\end{array}$$

where

$$\begin{array}{}\text{(A22)}& {\displaystyle }\mathrm{\Phi }\left(\mathit{\phi };c_{\mathrm{1}}\right)=\left\{\mathit{\phi }{}_t=\mathit{\phi }\left(\cdot -c_{\mathrm{1}}m\right):m\in {\mathbb{Z}}^{\mathrm{2}}\right\},\text{(A23)}& {\displaystyle }\begin{array}{rl}& \mathrm{\Psi }{}^{\mathrm{c}}\left(\mathit{\psi };c\right)=\mathit{\{}{{\mathit{\psi }}^{\mathrm{c}}}_{j,k,m}={{\mathit{\psi }}^{\mathrm{c}}}_{j,k,m}+i({\mathcal{H}}_{\left(\mathrm{1},\mathrm{0}\right)}\text{T }\mathit{\psi }{)}_{j,k,m}:\\ & j\ge \mathrm{0},\left|k\right|\le \lceil {\mathrm{2}}^{\frac{j}{\mathrm{2}}}\rceil ,m\in {\mathbb{Z}}^{\mathrm{2}}\mathit{\}},\end{array}\text{(A24)}& {\displaystyle }\begin{array}{rl}& \stackrel{\mathrm{̃}}{\mathrm{\Psi }}{}^{\mathrm{c}}\left(\stackrel{\mathrm{̃}}{\mathit{\psi }};c\right)=\mathit{\{}{{\stackrel{\mathrm{̃}}{\mathit{\psi }}}^{\mathrm{c}}}_{j,k,m}={{\stackrel{\mathrm{̃}}{\mathit{\psi }}}^{\mathrm{c}}}_{j,k,m}+i({\mathcal{H}}_{\left(\mathrm{0},\mathrm{1}\right)}\text{T }\stackrel{\mathrm{̃}}{\mathit{\psi }}{)}_{j,k,m}:\\ & j\ge \mathrm{0},\left|k\right|\le \lceil {\mathrm{2}}^{\frac{j}{\mathrm{2}}}\rceil ,m\in {\mathbb{Z}}^{\mathrm{2}}\mathit{\}}.\end{array}\end{array}$$

Correspondingly the discrete complex cone-adapted shearlet transform is given by the mapping,

$$\begin{array}{}\text{(A25)}& \begin{array}{rl}& f\to {\mathrm{SH}}_{\mathit{\phi },\mathit{\psi },\stackrel{\mathrm{̃}}{\mathit{\psi }}}^{\mathrm{c}}f\left(m^{\prime },(j,k,m\right),(\stackrel{\mathrm{̃}}{j},\stackrel{\mathrm{̃}}{k},\stackrel{\mathrm{̃}}{m}))=\\ & \left(f,{\mathit{\phi }}_{m^{\prime }},f,{{\mathit{\psi }}^{\mathrm{c}}}_{j,k,m},f,{{\stackrel{\mathrm{̃}}{\mathit{\psi }}}^{\mathrm{c}}}_{\stackrel{\mathrm{̃}}{j},\stackrel{\mathrm{̃}}{k},\stackrel{\mathrm{̃}}{m}}\right).\end{array}\end{array}$$

A5 Edge and ridge detection using the complex shearlet transform

The behaviour of the coefficients of the even-symmetric and odd-symmetric shearlets can be used to detect edges and ridges. An edge measure for an image f∈L²(R²), a location x∈R², and a shear parameter s are given as (Reisenhofer, 2014; King et al., 2015)

$$\begin{array}{}\text{(A26)}& \begin{array}{rl}& {\mathbf{E}}_{\mathit{\psi }}\left(\mathbf{f},\mathbf{x},\mathbf{s}\right)=\\ & {\displaystyle \frac{\left|{\sum }_{a\mathit{ϵ}A}\mathrm{Im}\left(f,{{\mathit{\psi }}^{\mathrm{c}}}_{a,s,x}\right)\right|-{\sum }_{a\mathit{ϵ}A}\left|\mathrm{Re}\left(f,{{\mathit{\psi }}^{\mathrm{c}}}_{a,s,x}\right)\right|}{\left|A\right|{max}_{a\mathit{ϵ}A}\left|\mathrm{Im}(f,{{\mathit{\psi }}^{\mathrm{c}}}_{a,s,x})\right|+\mathit{\epsilon }}},\end{array}\end{array}$$

where $A\subset {\mathbb{R}}^{+}$ is a set of scaling parameters, ψ is a real valued symmetric shearlet, and ϵ prevents division by zero. The complex shearlet-based edge measure can give approximations of the tangential directions of an edge. A line measure or ridge measure is obtained by interchanging the role of the even-symmetric and odd-symmetric shearlets (Reisenhofer, 2014; King et al., 2015):

$$\begin{array}{}\text{(A27)}& \begin{array}{rl}& {\mathbf{L}}_{\mathit{\psi }}\left(\mathbf{f},\mathbf{x},\mathbf{s}\right)=\\ & {\displaystyle \frac{\left|{\sum }_{a\mathit{ϵ}A}\mathrm{Re}\left(f,{{\mathit{\psi }}^{\mathrm{c}}}_{a,s,x}\right)\right|-{\sum }_{a\mathit{ϵ}A}\left|\mathrm{Im}\left(f,{{\mathit{\psi }}^{\mathrm{c}}}_{a,s,x}\right)\right|}{\left|A\right|{max}_{a\mathit{ϵ}A}\left|\mathrm{Re}(f,{{\mathit{\psi }}^{\mathrm{c}}}_{a,s,x})\right|+\mathit{\epsilon }}}.\end{array}\end{array}$$

Both the edge and ridge measures given above are inspired from the phase congruency measure of Kovesi (2000). The edge and ridge measures are almost contrast invariant.

Table A1 Shearlets.