The routine application of digital survey technologies such as terrestrial lidar and photogrammetry to the characterization of fault and fractures in outcrop over the past decade has resulted in major advances in terms of the efficiency of discontinuity data acquisition. However, the reliance upon mesh- and point-cloud–based analysis approaches means that data sets obtained from these sources commonly offer heavily abstracted views of the measured fracture network due to the limited resolution of the input model. Here, we present an alternative approach that combines conventional two-dimensional (2D) image analysis with ray-tracing techniques to extract three-dimensional (3D) fracture trace maps from photogrammetrically calibrated image sequences. These 3D trace objects may be interrogated to obtain fracture network properties (trace length, intensity, and connectivity), with probabilistic methods used to estimate fracture orientation for high collinearity traces.

Our approach possesses a number of advantages over existing digital surface reconstruction-based methods, with the use of a 2D pixel-based approach allowing established image-processing routines (e.g., edge detection/connected components analysis) to be applied to the characterization of fracture and fault properties. Moreover, the innately high resolution of the input images results in practically lossless 3D fracture trace representation, limiting truncation effects. As a result, the method is capable of resolving local variability in higher-order fracture properties such as fracture intensity, which are difficult to derive using existing approaches. We demonstrate the approach on pervasively faulted Permian age exposures of the Vale of Eden Basin, UK.


Structural heterogeneities (e.g., faults and fractures) within the upper crust form in response to a range of stress states (e.g., isotropic, normal, shear, or a combination thereof), which may emanate from elevated fluid pressure gradients, as well as lithostatic, tectonic, or thermal loading. The preeminence of geological discontinuities in governing both geomechanical behavior and crustal circulation of geofluids means that understanding the spatial distribution of fault and fracture systems is critical to resolving numerous problems within the fields of structural geology, rock mechanics, hydrogeology, and petroleum geoscience (e.g., Aydin, 2000; Dyer, 1988; Hitchmough et al., 2007; Odling et al., 1999; Priest, 1993).

Though borehole and geophysical data sets have been employed extensively in the study of naturally fractured rocks (e.g., Barton et al., 1995; Grasmueck, 1996), outcrop data have been used more commonly as a medium for understanding the geometry and distribution of fracture and fault networks (e.g., Florez-Nio et al., 2005; Ghosh and Mitra, 2009; Gillespie et al., 2001; Odling et al., 1999). The reliance upon outcrop data sets in these investigations is due to their ability to provide information of fracture intensity, connectivity, and size, which are not directly obtainable from subsurface data sets.

Conventionally, the characterization of fractures in outcrop has been accomplished using direct measurement at the rock face using structural surveys (scanlines and/or cell mapping; e.g., Mauldon et al., 2001; Priest and Hudson, 1981) or interpretation of remotely sensed data (i.e., lineament and/or trace maps; e.g., Kim et al., 2004). However, both approaches suffer manifold disadvantages, with direct surveys being biased, inefficient, and prone to censoring (Kemeny et al., 2006), and with photogeologic interpretation of fracture networks generally being incapable of resolving discontinuity orientation (and thus size distributions and intensity segmented by identified fracture set).

Over the past decade, digital survey technologies such as terrestrial lidar and digital stereophotogrammetry have offered an attractive alternative to conventional fracture surveys, providing the outcrop coverage akin to that from photographic surveys, while producing statistically dense data sets of fracture orientation (e.g., Kemeny et al., 2006; Sturzenegger and Stead, 2009). As a result of the large data volumes such methods are capable of creating, considerable efforts have been expended by workers to develop automated and semi-automated methods of digital fracture data acquisition and analysis. In this regard, the extraction of discontinuity planes directly from point-cloud or triangular irregular network surface reconstructions has received the greatest attention, with a number of procedures being proposed (e.g., García-Sellés et al., 2011; Kemeny et al., 2006; Roncella and Forlani, 2005; Van Knapen and Slob, 2006).

While such approaches enable fracture orientation set distributions to be rapidly determined, the reliance upon point-cloud or surface reconstructions as a medium for fracture identification has major disadvantages. Given that discontinuities may be manifest either partially or completely as traces without any discernible surface topography, the coplanarity and curvature of the outcrop surface may not be diagnostic of the presence of fractures. Many discontinuity types such as veins (e.g., Gillespie et al., 2001) and deformation bands (e.g., Sternlof et al., 2006), as well as joints exposed on excavated rock cuts (i.e., tunnels and quarry walls) typically have minimal surficial expression that precludes the application of mesh- and point-cloud–based methods to their analysis. As a result, point-cloud and mesh-based outcrop surface models tend to offer an abstracted view of a target fracture network, with discontinuities becoming partially censored or omitted (equivalent to truncation and curtailment effects within conventional fracture surveys: Priest, 1993; Seers and Hodgetts, 2013). These censoring effects may impact significantly on the determined network properties, and especially upon those that are spatially dependent, such as fracturing intensity and connectivity (and thus contingent rock mass properties such as permeability and strength). For example, Seers and Hodgetts (2013) observed that reductions in fracture density related to these resolution effects resulted in broad-scale alterations in the permeability tensors of derivative discrete fracture network models.

In order to investigate such parameters, it is generally desirable to obtain a discrete and unabridged representation of all observable fractures within a surveyed domain. Given fracturing intensity and connectivity are widely recognized as being first-order controls over fracture permeability and rock mass strength (Dershowitz et al., 2000; Rogers et al., 2013), this inability to satisfactorily account for these higher-order fracture properties therefore represents a major weakness in current digital fracture analysis approaches.

1.1. Fracture Trace Maps

The most widely employed method of producing spatial representations of fracture networks in outcrop is via the creation of trace maps, whereby individual fractures are digitized from photo-imagery (e.g., Barton et al., 1995; Le Garzic et al., 2011; Odling, 1997; Renshaw, 1997). The construction of fracture trace maps from optical images provides a convenient means of representing the trace pattern and termination relationships of a given network, with the delineation of traces on an image plane being a more straightforward operation than polyline interpretation within a 3D object space. The digitization of fractures for trace-map production may be conducted manually using standard CAD-based tools (e.g., Odling, 1997) or through the application of unsupervised edge detection and segmentation (e.g., Lemy and Hadjigeorgiou, 2003; Reid and Harrison, 2000).

In addition to the relative simplicity of feature extraction, the information-rich qualities of images produced by modern digital cameras offer further advantages relating to the minimization of spatial aliasing-based truncation of fractures from digital surveys. High-megapixel imagery provides greater pixel coverage of targeted exposures when compared to the typical vertex spacing of lidar or digital stereo photogrammetry–derived point clouds captured from equivalent ranges, meaning that smaller features can be detected. Moreover, because optical images comprise multichannel displays of the visible light spectrum (typically using the RGB additive model), color, and tonal variations of the exposed fracture network may be used to guide feature identification of traces, which are expressed with minimal surface topography.

Nevertheless, the use of 2D imagery for trace-map production does have noteworthy limitations. For example, the singular viewing position and restricted field of view (FOV) of the camera may result in the partial occlusion of the surveyed fracture network (Umili et al., 2013). Furthermore, the representation of traces delineated from the image plane is purely bi-dimensional, limiting the ability to estimate true fracture trace length (i.e., chord distance). This also means that only 2D measures of fracture orientation can be obtained directly from the mapped traces (Otoo et al., 2011); although 3D orientations may be estimated if a posteriori knowledge of orientation distributions exists (i.e., from manual measurements; Kemeny and Post, 2003).

In a work pertinent to the present study, Umili et al. (2013) demonstrate a method that utilizes break lines identified using vertex principal curvature as a proxy for trace occurrence, enabling the unsupervised generation of 3D trace maps. While the improved representation of fracture traces their method offers is welcome, the reliance upon surface topography as the key diagnostic for fracture occurrence limits the usefulness of their approach. Indeed, the authors concede that their method is only applicable in situations where fractures may be identified as an edge of an exposed surface (i.e., natural exposures of rock joints). In a further work closely related to the present study, Vasuki et al. (2014) present a method for the semi-automated extraction of 3D structural discontinuities from orthorectified images and digital elevation models (DEMs), derived from unmanned aerial vehicle (UAV)–based photogrammetry. By using optical images as the primary medium for trace-map delineation, Vasuki et al. overcome some of the limitations inherent to the mesh-based approach presented by Umili et al. (2013). The reliance upon rasterized DEMs as the basis for 3D feature extraction however means that their approach is not applicable to arbitrary surfaces, effectively limiting their method to those generated using quasi-nadir view imagery (i.e., imagery captured using airborne sensory platforms).

In this work, we present an approach that retains many of the benefits associated with 3D trace-map production using photogeologic interpretation of optical images (i.e., Vasuki et al. 2014), while generating a continuous three-dimensional representation of the surveyed fracture network (i.e., Umili et al. 2013). Our method harnesses calibrated images, common to both photogrammetric and lidar digital surveys, and an optimized ray-tracing engine to extract low-aliased three-dimensional trace maps from digital surface models, with a 3D lineament merger algorithm used to marry contiguous traces from discrete images. Drawing upon the work of Seers and Hodgetts (2016), we employ a probabilistic solver to draw orientation estimates from the mapped traces. By locating the trace’s plane of best fit using a probabilistic rather than an analytical solution, we are able to obtain orientation estimates for features with vertex geometries that are not conducive to least-squares plane fitting (i.e., approximately collinear traces). By employing this framework, measures of connectivity and intensity can be derived for discrete orientation sets with minimal truncation of the surveyed fracture network. Moreover, the availability of dense, point-based representations of individual fracture traces mapped to the parent mesh allows variability in these properties to be calculated across the surface of the input digital surface model through a spherical kernel-based schema. We demonstrate the potential of our approach on pervasively faulted exposures of Permian sandstones of the Vale of Eden Basin, UK.


In contrast to previous studies, which relied upon direct mesh-based digital fracture analysis approaches, we utilize calibrated images as the proprietary medium for fracture identification. The camera projection matrix associated with detected 2D traces on the image plane is used to map extracted pixel features from digital surface models to equivalent world coordinates. The composition of the projection matrix and the procedures available for obtaining a geometric camera calibration are summarized below.

Expressed using homogenous coordinate notation, a 2D point m in pixel coordinates on the image plane, [u v]T, is related to its equivalent point M in world coordinates in the 3D object space, [xwywzw]T, by a coincidental ray passing through the camera’s optical center O, assuming a pinhole camera model (Fig. 1). Thus, the relationship between m and M can be described by (Zhang, 2014): 
where R is a 3 × 3 rotation matrix and t is a 3 × 1 translation vector describing the attitude and location of the world coordinate system axes and origin in relation to the camera coordinate system. Concatenation of R and t gives the extrinsic parameters of the camera 
C is the camera intrinsic matrix, containing the internal parameters of the camera: 
where γ is the skew coefficient between the camera coordinate system’s x and y axis (equal to 0, if pixels are quadratic) and u0, v0 are pixel coordinates of the principal point (ideally corresponding to the image center). Parameters αx and αy are products of Euclidean focal length f and scaling factors mx and myx = f · mx, αy = f · my), which relate f to distance in pixels (Hartley and Zisserman, 2003). Other significant nonlinear intrinsic parameters, such as radial distortion coefficients, cannot be encompassed by this linear camera model but may be estimated by current camera calibration procedures (e.g., Fitzgibbon, 2001). The product of C and [R|t] gives the camera projection matrix P = K[R|t] (referred to as the geometric camera calibration), with camera pose C, expressed in world coordinates being related back to the extrinsic matrix via (Hartley and Zisserman, 2003): 

In practice, the intrinsic parameters of a camera may be ascertained either through photogrammetric calibration or self-calibration (Zhang, 2000). Photogrammetric calibration relies upon the displacement of a calibration apparatus (3D reference objects, planes, and lines) whose geometry is characterized to high degree of precision, across a sequence of narrow baseline images (Zhang, 2014). Conversely, self-calibration (e.g., Faugeras et al., 1992) uses correspondences between sequences of uncalibrated images captured across static, unstructured scenes to obtain constraints upon the camera’s internal parameters, being more convenient than photogrammetric calibration at the expense of precision.

Derivation of the camera extrinsic parameters varies with the type of digital survey technique being employed. Close-range stereophotogrammetry operates upon the principles of epipolar geometry, utilizing the triangulation of intersecting back-projected rays from multiple overlapping images to estimate 3D scene structure and camera pose. Current multi-view dense stereo processing pipelines developed within the field of computer vision have largely automated this procedure. A popular approach is to use sparse image feature correspondences (e.g., using scale invariant features; Lowe, 1999) to initialize the scene structure and camera geometry (structure from motion estimation [SfM]; e.g., Snavely et al., 2006). Using the estimated camera extrinsics, sparse matches can then be iteratively expanded to nearby locations to construct a dense model, with visibility constraints used to suppress false matches (multi-view stereo [MVS]; e.g., Furukawa and Ponce, 2010).

In the case of terrestrial lidar surveys, camera extrinsic parameters may be associated to images captured from a scanner-mounted camera by applying a mount calibration, describing the known rigid transformation between the laser beam detector and the optical center of the attached camera. If this transformation is unknown (i.e., in the case of images captured off scanner), the camera extrinsic parameters may be determined by manually identifying mutual keypoints between the lidar point cloud and mapped imagery; although this process may be time consuming for large numbers of images. Unsupervised pose estimation of uncalibrated images using laser scans of equivalent scenes is, however, an ongoing area of research, largely motivated by developments in robotic navigation (Pandey et al., 2012; Scaramuzza et al., 2007).

As previously discussed, the relative resolution and richness of information contained within optical images, coupled with the highly evolved nature of 2D image-processing techniques makes digital photography attractive for mapping natural fracture networks. Exploiting the relationship between the camera and world coordinate system described by a camera projection matrix would appear to provide a suitable means of extending two- dimensional image features into the 3D domain. Indeed, the use of calibrated images to create projective textures is a well-established practice in the computer graphics and vision communities (e.g., Debevec et al., 1996), and the method has previously been used by a number of workers (Casini et al., 2011; Gillespie et al., 2011) to enhance the visualization of fracture systems on geological surface reconstructions. The method we present here is broadly analogous to texture mapping, utilizing pinhole camera projection to map detected pixel features to their location in the world coordinate system, constrained by ray-mesh intersections. In contrast to the monolithic mesh textures produced by conventional UV mapping, we are able to obtain discrete 3D lineament objects while avoiding image decimation and performance bottlenecks associated with the creation and manipulation of large projective textures (Hua et al., 2004; Wang et al., 2001). In the following sections, we will describe our implementation of 3D trace-map extraction using calibrated image sequences and present an application of the technique to the characterization of sandstone-hosted shear bands.


Developed using the MATLAB® language, our 3D trace-map extraction software comprises two graphical user interface (GUI)–based modules: (1) a 2D image editor used for the extraction of 3D fracture traces from calibrated images; and (2) a 3D viewer that facilitates the merger of traces from neighboring images, as well as the analysis of key network properties such as fracture orientation, size, connectivity, and intensity. The procedure for extracting and analyzing 3D trace maps from calibrated images is described below.

3.1. Pixel-Based Trace Delineation

To ensure that feature extraction is both efficient and comprehensive, we employ a combination of unsupervised edge detection and manual curve fitting to delineate fracture traces from input images. Following Vasuki et al. (2014), unsupervised trace extraction is based upon the phase congruency feature detector of Kovesi (1999). In contrast to commonly used edge detectors that exploit pixel intensity gradients (e.g., Canny, 1986), phase congruency is based upon the observation that in-phase Fourier components in the frequency domain of an image represent the signals of edges (Kovesi, 1999). Because the location of coincident sine waves is insensitive to shifts in illumination and contrast, phase congruency tends to outperform detectors that rely upon local gradient magnitude in natural scenes, such as rock exposures, where shifts in luminance are common (Nixon and Aguado, 2012).

The occurrence of noise is generally detrimental to the performance of edge detectors. We therefore apply the edge-preserving, non-local means (NLM) filter (Figs. 2A and 2B) of Buades et al. (2005) prior to edge extraction and then use hysteresis thresholding and skeletonization to extract a binary image from the edge map output from the phase congruency feature detector (Figs. 2C and 2D). Discontinuous fracture trace segments are initially extracted by splitting foreground pixels in the binary image at junction points (e.g., Di Ruberto, 2004), with linkage of disparate segments achieved using a novel line-growing algorithm (Fig. 3A). This line-growing operator uses a buffered search path oriented in the direction of the mean vector of a given segment, extending from its end nodes to catch potential matched edges. False positives are rejected on the basis of a user-defined angular threshold, with a weighting scheme based upon angular deviation and distance from the parent line used to determine optimum matches.

The presence of low-contrast edges in addition to the requisite filtering of images for edge detection means that unsupervised trace extraction is likely to result in a partial representation of the network. For example, only 31% of the observed trace network in Figure 2A could be mapped using the automated trace extraction procedure presented above. As a consequence, it may be necessary to manually fit curves to the image to ensure fracture traces are comprehensively delineated (Fig. 3B). In many cases, the strongest image edges are not produced by the target fracture trace network, but by extraneous features such as the outcrop edge, vegetation, shadows, or weathering patterns and/or staining. If such features are prevalent, automated trace extraction will produce a significant proportion of false positives, which must be manually identified and removed. In addition to environmental factors, the parameterization of the NLM filter, the phase congruency edge detector, hysteresis thresholding, and the line-growing operator will also impact significantly upon the results of automated fracture trace extraction. Indeed, optimum parameters used at each stage may vary widely between different images of the same outcrop and will depend upon the physical and optical characteristics of the scene, as well as the curvature and contrast of the targeted trace features. The lack of an all-encompassing parameterization framework for automated trace extraction, coupled with the considerable amount of effort required for vetting and editing its output, may severely reduce the utility of unsupervised fracture delineation. Dependent upon the level of contrast between fractures and their surrounding host rock, and the number of artifacts and unwanted edges present, manual curve fitting may prove to be the most efficient and accurate approach to fracture trace delineation for many images.

The final stage in the 2D trace-map delineation process is determining connectivity relationships between the mapped traces (Fig. 3C). We identify connected trace components by binarizing each trace map and applying the two pass region labeling algorithm of Haralick and Shapiro (1992). While ascertaining pixel connectivity between line segments in a binary image is a relatively trivial operation, the equivalent procedure in three dimensions is non-robust due to the extreme low probability of line segment intersections in a 3D object space. Later, we will describe the creation of a global 3D connectivity map from the connected region labels of individual images.

3.2. Optical Ray Tracing

In the present work, pixel-based representations of a fracture trace pattern are extrapolated into three dimensions using optical ray tracing. Such an approach was previously suggested as a possible means to map 2D features from the image to real-world coordinates by Vasuki et al. (2014) but was never implemented within their study. Not to be confused with the equivalently named technique for calculating propagating particle or wave paths, ray tracing involves the back-projection of rays from located cameras into synthetic scenes, enabling complex light-object interactions to be modeled (e.g., shadows, reflections, and refracted light; Cook et al., 1984). Though used principally in computer graphics for illumination rendering, we utilize ray tracing to determine correspondences between mapped pixel features on calibrated images and triangular irregular network-based surface models of the viewed rock exposure.

We initialize the tracer from a single ray, oriented in nadir view and emanating from the origin of the world coordinate system, representing the optical axis of the camera. Iterative replication and rotation of this initial unit vector yields a tree of rays (Fig. 4A), with each constituent ray corresponding to a single pixel of the input image (related by ax and ay). A mask is then applied, occluding rays that are not associated to indexed pixels that represent mapped fracture traces. Rotation of the ray paths to the attitude of a calibrated camera is achieved by locating the rotational axis forumla (where forumla and Σ(ωx, ωy, ωz) = 1) by finding the cross product between the vectors, a = [0 0 1]T, and the directional vector, b, obtained from the digitized trace (forumla). The matrix, forumla, describing the rotation between a and b around forumla by the angle, η, is then obtained by solving the Rodrigues rotational formula (e.g., Murray et al., 1994): 

A secondary rotation between the direction cosines perpendicular to the transformed vector T(a) and b around their shared axis is then required to remove the remaining degree of freedom between the two systems (Fig. 4B).

Because ray tracing is a computationally expensive operation, it is desirable to limit the number of potential mesh triangles included in the search for ray-triangle intersections (Fujimoto et al., 1986). To reduce the number of potential mesh triangle candidates, we use an octree-based spatial partitioning scheme to group mesh triangles into multiple bounding boxes aligned to the local coordinate axis (Meagher, 1982), with the ray-box intersection algorithm of Smits (2002) used to constrain which of these spatial bins is intersected by an incoming ray. Finally, the point of intersection between a given ray and a mesh object is determined using the fast ray-triangle intersection algorithm of Möller and Trumbore (1997), resulting in the extraction of vertex information from the input pixel-based 2D fracture trace map (Fig. 3D). Both the ray-box and ray-triangle intersection solvers are hardware accelerated for speed. Benchmarked on a HP Z420 workstation (Intel Xeon E5-1650 quad core with Nvidia K4000, 64GB RAM and running Windows 7 Enterprise Service Pack 1), the feature extractor computes the 3D locations of 150 trace objects composed of ∼400k pixels mapped on a 8.4 megapixel image from a ∼3.2 million face tetrahedral mesh (∼340k faces in the FOV) in ∼65 seconds. Computation time grows proportionately with increasing input image size and mesh density within the image FOV. For example, locating the 3D coordinates of an equivalent trace data set (150 trace objects and/or ∼800k pixels) mapped on a 16.8 megapixel image from ∼670k faces using the hardware configuration described above takes ∼147 seconds.

3.3. Nearest-Neighbor–Based 3D Lineament Merger

Due to the limited field of view of the camera and the common occurrence of line-of-sight obstructions during image capture, multiple calibrated images will typically be required to map an exposed fracture network. As a consequence, ray-tracing–based extraction of trace features is likely to result in a number of discrete 3D trace-map patches, with fractures spanning the FOV of congruent and/or neighboring images occurring as disconnected curves. We provide the means to interactively select and merge these discontinuous traces in the GUI of our 3D viewer. In cases where fracture geometries are relatively simplistic, correspondences between disconnected traces may be readily identified using the 3D curve data set alone. However, when the fracture trace pattern is complex (i.e., tightly clustered and/or anastomosing discontinuities), relationships between traces mapped from congruent images may be difficult to decipher. In such cases, we utilize pixel-resolution 3D renderings of the paired images and extracted trace objects to establish correspondences (Figs. 5A and 5B).

The merger of vertex-based representations of lineaments in a 3D space is complicated by the existence of a number of potential alignment geometries between the two curves (Fig. 5C). We utilize nearest-neighbor relationships between the vertices of the two selected curves to be merged, initially to determine whether a condition of overlap exists, and if so, where the point of overlap occurs (Fig. 5D). The simple case of non-overlapping curves is flagged when the nearest-neighboring vertices between the two lineaments are found on their respective end nodes. Should overlap be detected, the non-overlapping segments of the curve are determined by finding the point along the opposing curve where the located nearest-neighboring vertex becomes stable, with the paired vertices separated by the smallest distance providing the point of overlap.

In the case of non-overlapping lineaments, the merger is achieved by trimming the opposing end segments, with the missing center section populated by interpolating points along a parametric spline curve fitted to the trimmed portions of the input fracture traces (e.g., Yamaguchi, 2012). For overlapping lineaments, the centroid between the overlapping segments forms the center portion of the merged trace. In similitude to the non-overlapping case, trimming of opposing end segments and subsequent spline interpolation is used to create a smooth curve during the merger process. It should be noted that the merger process produces an averaged geometry for overlapping traces, effectively resulting in information loss. As a general recommendation, the length of overlap between paired traces should be minimized, which may be achieved prior to 3D feature extraction by mapping features from images that exhibit nominal overlap. In addition, basic tools are provided in the 3D viewer for selecting and trimming traces. These tools enable users to minimize overlap between paired trace features subsequent to the extraction of 3D trace-map patches. A further point to note is that the quality of the merger will depend upon the degree of spatial correspondence between the merged traces, with major deviations between the 3D curves resulting in artifact geometries after merging. The degree of spatial correspondence between paired traces is largely dependent upon the density and error characteristics of the input mesh, the accuracy of the feature delineation from the image plane, and the quality of the image calibration.

In addition to aiding the restitution of the exposed fracture trace array, this lineament merger procedure provides the means to update the connectivity indices generated during image region labeling, enabling global fracture connectivity to be quantified. The final stage of lineament merger is to update the vertex connectivity indices associated with the newly merged curves to create new global connectivity indices. After all disconnected curves have been merged, all remaining indices assigned during image region labeling are similarly updated to form a global trace connectivity map of the fracture array.

3.4. Fracture-Trace Orientation

As alluded to previously, ascertaining discontinuity orientation distributions and assigning individual measurements to orientation sets is a prerequisite stage in the analysis of many of the properties of natural fracture systems. As with any 3D point set, the orientation of a vertex-based representation of a fracture trace may, in theory, be determined by finding the best- fitting plane through its constituent points. The unconditional application of plane-fitting routines to estimate discontinuity orientation from fracture traces should however be cautioned against. Based upon the vertex properties of faults and fracture traces mapped from decimeter to regional scale, Seers and Hodgetts (2016) demonstrate that the majority of digitized traces tend to have vertices that make them unsuitable candidates for least-squares plane-fitting routines. Using the ratios between the eigenvalues (λ1, λ2, λ3) of covariance matrix of each trace feature’s point set as a measure of vertex coplanarity: M = ln(λ13), and collinearity: K = ln(λ12)/ln(λ23) (Fernández, 2005; Woodcock, 1977), they find that most structural lineaments have vertices with collinearity values that fall outside the threshold for reliable orientation estimates (K < 0.8) established by Fernández (2005). Furthermore, they identify degenerate cases whereby the least-squares plane fitted to a fracture trace will indicate the outcrop rather that discontinuity orientation, if the asperity heights of the associated discontinuity locally exceed those of the outcrop along their intersecting curve.

In previous work, iterative voting-scheme–based model estimation using random samples consensus (RANSAC; Fischler and Bolles, 1981) has been proposed as a suitable method for the estimation of discontinuity orientation from 3D fracture traces (i.e., Vasuki et al., 2014; Thiele et al., 2015). It should be noted however that even robust parametric estimators such as RANSAC or the Hough transform, in addition to analytical solutions, such as planar regression or singular value decomposition, are likely to lead to erroneous results if applied to the degenerate cases described above. The application of any fitting routine should therefore be made in conjunction with both an assessment of the vertices geometrical properties and an intersection test between the trace and the outcrop planes of best fit, if measurements are to remain credible.

Following the work of Seers and Hodgetts (2016), we rely upon probabilistic orientation estimates to determine structural attitude from 3D fracture trace maps. Informed by the results of Monte Carlo experiments, their approach relies upon the assumption that the fracture pole vector associated with a digitized trace lies upon the great circle formed perpendicular to the arbitrary axis of rotation set by the trace’s end nodes. The probability density characteristics around this rotational axis are assumed to be governed by (Seers and Hodgetts, 2016; Fig. 6A): (1) The underlying probability distribution of estimated pole vectors with respect to a given interval of M/K—assumed to be von Mises (circular normal) distributed; (2) the angle of incidence (β) between the estimated best-fit plane of a trace and that of the outcrop locally along their intersection curve (e.g., Terzaghi, 1965); and (3) the proximity of the prospective great circle to orientation cluster centroids identified from the wider discontinuity network (assumed to be von Mises-Fisher distributed).

Under this framework, discrete probability density functions (PDFs) are constructed around the great circle formed perpendicular to the trace end nodes representative of the probability characteristics described above. A linearly weighted aggregation scheme (linear opinion pool; Clemen and Winkler, 1999; Stone, 1961) is then employed to combine individual PDFs into a combined probability model from which the orientation estimate is stochastically drawn (the reader is referred to Seers and Hodgetts [2016] for a more thorough conceptual and mathematical underpinning of the approach). Extensive numerical benchmarking of the solver using synthetic lineaments extracted from intersecting fractional Brownian surfaces has shown that this probabilistic approach offers significant gains in orientation estimate precision over equivalent analytical solutions (Figs. 6B–6D; Seers and Hodgetts, 2016).

The great circle upon which the fracture plane normal is assumed to lie is initialized at its closest point to the eigensolution of the trace’s best-fit plane pole vector. Based upon the Monte Carlo experimental results of Seers and Hodgetts (2016), pole vectors along this circular arc are modeled as von Mises distributed. The concentration parameter k of the von Mises distribution, which acts as a descriptor of pole vector dispersion upon the unit circle, is selected from experimentally derived values constrained to different intervals of K and M, based upon the corresponding collinearity and coplanarity characteristics of the input trace. K-means clustering is then used to identify discontinuity clusters from fracture plane facets (sensu: Kemeny et al., 2006) and higher orientation quality traces based upon a user-defined threshold of M and K, with points along the great circle falling within the 99% spherical aperture of each cluster assumed to be von Mises-Fisher oriented. The final class of PDFs describes the probability intersection of a fracture at a given angle of incidence to the outcrop, which is relatable to the normalized Terzaghi weighting function (Terzaghi, 1965) by inverse proportionality (Seers and Hodgetts, 2016). Having constructed each discrete PDF, a hierarchical weighting system is applied prior to the aggregation of the probability density functions. The primary weighting factor is related to the collinearity characteristics of the input trace, with the clustering of discontinuity network and the angle of incidence between the orientation estimate and the outcrop making a secondary and tertiary contribution. Subsequent to the linearly weighted aggregation of each PDF, the probabilistic orientation estimate is obtained by taking the principal direction of 1000 uniformly random draws from the combined probability model.

Although the probabilistic framework described above offers marked improvements of the equivalent analytical solution in terms of orientation precision, it is still recommended that the user assess the quality of the output from the solver. This operation is conducted by selecting trace objects within the 3D viewer. Upon selection, a glyph (disk) is fitted through the trace objects centroid displaying its estimated orientation. The user may then update the orientation estimate manually by selecting a point on the arbitrary axis of rotation set by the traces end nodes, displayed around the trace feature’s centroid.

3.5. Fracture Trace-Length and Termination Styles

A major advantage of constructing a 3D digital representation of the exposed trace network is that accurate measurements of trace length can be easily derived. Within the 3D viewer, both the surficial length (the length of the discontinuity-outcrop curve of intersection) and the chord length (distance between the trace end nodes) can be computed per orientation set. Having established trace-length measurements for a given fracture subpopulation, Anderson-Darling tests (Anderson and Darling, 1954) are used to determine the goodness of fit between the empirical trace-length distributions and parametric distributional forms commonly used to describe fracture size in outcrop (lognormal, exponential, and power law; e.g., Baecher et al., 1977; Hatton et al., 1994; La Pointe, 2002; Priest and Hudson, 1981).

In addition to the determination of trace-length characteristics, termination styles of the trace network can also be assessed, with the end nodes of each trace object being assigned to one of the following classes: (1) blind termination (B)—the fracture terminates within the rock mass; (2) abutting (T) termination—the fracture terminates against another mechanical discontinuity; (3) crossing (X) termination (Anderson and Darling, 1954)—the fracture intersects another mechanical discontinuity; and (4) censored termination (C)—the fracture termination cannot be observed.

Termination styles are initially assigned using unsupervised estimation, based upon both the connectivity indices of individual traces and the proximity of their end nodes to the nearest neighboring vertices of the wider network. The user may then update the termination classification manually by selecting individual traces within the GUI. Ratios between termination classes 1–3 may then be used to quantify fracture connectivity (Sahini and Sahimi, 2003), with class 4 forming the basis for the non-parametric correction of truncated fracture trace length using Kaplan Meier estimation (e.g., Odling, 1997).

3.6. Fracture Intensity

It is widely recognized that the degree of fracturing present within the upper crust controls permeability and mechanical competence (Dershowitz et al., 2000; Rogers et al., 2013). To meet the task of quantifying fracture abundance in outcrop, we have developed a measurement framework that utilizes 3D trace maps to gain per-orientation set estimates of areal fracturing intensity (P21; Dershowitz and Herda, 1992) across the exposure face. A roving spherical kernel is used to determine fracture length per unit of outcrop surface area around each vertex of the parent digital surface model (Fig. 7), enabling the intrinsic variability in fracturing intensity to be represented as a vertex attribute map. This method is equivalent to the approach presented by Pearce et al. (2011) for calculating areal density measures for manually digitized traces extracted as polylines from laser scan data sets. Based upon the assumption that orientation clusters are von Mises-Fisher distributed (Fisher, 1953), measured values of P21 are then transformed into more geologically useful volumetric intensity measures (P32; Dershowitz and Herda, 1992) by adapting the conversion factors of Wang (2005).

It should be noted that the selection of the kernel radius will have a significant impact upon the generated vertex map, with fracture intensity values tending toward a statistically homogeneous state as the kernel volume approaches the representative elementary volume (REV) of the fractured rock mass property (see Hill, 1963). Although it is desirable to match the kernel volume to the representative element of the fractured rock mass, the existence of the latter is not guaranteed (Neuman, 1988). Due to these conceptual difficulties, we suggest that in cases where the attribute maps are to be used for property modeling for fracture network model generation (i.e., through geostatistical simulation), the kernel volume should be matched to the targeted grid volume. As a test case, within this study we set the kernel radius in proportion to the outcrop’s edge length.


To demonstrate the potential of our approach, we have extracted a three-dimensional trace map from pervasively faulted sandstone exposures from the Vale of Eden Basin, UK; these exposures have been used to assess the geometrical properties of the fault network (fault orientation, trace length, connectivity, termination styles, and intensity).

4.1. Study Area

Inferred to have formed in response to extensional faulting along the Dent and Pennine fault systems during a period of regional-scale Early Permian rifting (Holliday, 1993; Underhill et al., 1988), the Vale of Eden Basin is a fault-bounded, north-south–trending sedimentary basin located within Northwest England (Fig. 8A). The basin comprises a mixed sequence of predominantly continental red-bed Permo-Triassic deposits that rest unconformably upon the Carboniferous basement (Fig. 8B). The present study area (Lacy’s Caves) is located within the central portion of the Vale of Eden Basin, exposing Lower Permian aeolian deposits of the Penrith Sandstone Formation. The outcrop is formed as a sequence of man-made excavations exposing a dense population of cataclastic shear bands (micro faults), which occur as thin tabular zones of compacted and crushed grains within the high-porosity sandstone host rock (Fowles and Burley, 1994). The shear bands developed within the cave system are analogous to subseismic- scale faults identified within several North Sea petroleum reservoirs (i.e., the Rotliegendes of the Southern North Sea; Fisher and Knipe, 2001). The characterization of the geometry, interconnectivity, and abundance properties of the exposed fault network therefore has considerable practical application.

4.2. Digital Outcrop Data Acquisition and Processing

We employ a combination of structure from motion estimation (SfM) and multi-view stereo reconstruction (MVS) to generate the camera calibrations and digital surface model required for 3D trace-map production (e.g., Favalli et al., 2012). A photographic survey of a single cave was conducted using a Sony Nex 5N mirrorless micro-dSLR camera equipped with a 16-mm pancake-style prime lens, producing 213 eight-megapixel narrow baseline images, captured at ∼1.5 m distance from the outcrop. At this range, the camera setup produces a pixel size on the outcrop surface of ∼600 µm, enabling millimeter-scale shear bands to be readily identified. For the purpose of model registration, four orthogonal ground-control points corresponding to north, south, east, and west were marked on the outcrop with the aid of a laser level, with the distances between each key point and the center of rotation of the survey apparatus recorded. These control points were subsequently imaged within the photographic survey.

A sparse point cloud and camera pose estimates were obtained from the Lacy’s Caves photo-survey data set using VisualSFM, an open-source, GUI-based application for 3D reconstruction using structure from motion estimation (Wu, 2011). The sparse point cloud and camera intrinsic-extrinsic parameters generated by this procedure were then used as the basis for dense scene reconstruction using the patch-based, multi-view stereo tool (PMVS2) of Furukawa and Ponce (2010). Dense reconstruction produced a 6.5 million vertex RGB point cloud, covering an outcrop area of ∼25 m2, with an approximate point density of 260k points per m2. Finally, the Poisson surface reconstruction algorithm of Kazhdan et al. (2006) was used to construct the digital surface model needed for trace-map extraction from the dense point-cloud output from PMVS2 (Fig. 9 and Interactive Model 1). An octree depth of 11 was selected for the Poisson surface reconstruction, offering balance between the reduction of the surface model data volume and the preservation of the vertex density characteristics of the input point cloud. The resultant triangular irregular mesh consisted of 1.6 million vertices, 3.2 million faces, with an average nodal density of ∼64k points per m2 (a down-sampled RGB colored mesh surface reconstruction of Lacy’s Caves is made available as part of the Supplementary Files1). Finally, the Cartesian coordinates of the ground control points mapped within Lacy’s Caves were obtained using the ray-tracing techniques presented in Section 3.2. and were used to calculate the similarity transformation between the arbitrary and reference coordinate frames using the method of Horn (1987).

4.3. Fault Network Extraction

A subset of 25 partially overlapping and/or adjacent calibrated images from the photogrammetric survey was used as the basis for 3D fault trace extraction, with both manual curve fitting and unsupervised edge detection used to map the discontinuities from the image plane (see Section 3.1.). Ray-tracing–based feature extraction produced 802 three-dimensional curves, which were subsequently edited and merged within the 3D viewer, resulting in 645 trace objects, providing the basis for fault network property characterization (Fig. 10 and Interactive Model 2).

4.4. Shear-Band Orientation

Guided by the experimental results of Seers and Hodgetts (2016), thresholds of K < 2 and M > 6 were imposed upon traces selected for the determination of fault-orientation clustering, equating to a sample 95% confidence level of ∼10°–15° around the axially constrained great circle (Seers and Hodgetts, 2016). Proto-clusters identified from this higher-orientation quality data set using the K-means method were then used to condition the probabilistic solver described in Section 3.4., enabling the stochastic estimation of discontinuity orientation for the remainder of the fault trace network. Approximately 10% of orientations were manually updated during the vetting procedure described in Section 3.

Using K-means clustering, two dominant fault-orientation sets are identifiable from the aggregated analytically and probabilistically–derived orientation data (Fig. 11 and Table 1): a dominant NW-SE–oriented set (Set 1: n = 508; direction of μ = 085/137; k = 9.8) approximately aligned with the Kirkoswald fault, which is located proximal to the study area (see Fig. 8B), and a subsidiary NE-SW–oriented set (Set 2: n = 66; direction of μ = 088/232; k = 13.4) coincidental to the other principal trend of faulting within the central Vale of Eden Basin. Comparison between the fault-orientation estimates obtained from our 3D trace maps and manually derived measurements obtained from Lacy’s Caves by Fowles and Burley (1994) reveals excellent agreement, with both the conjugate pairs of shear bands (Set 1 and Set 2) being readily identifiable from both data sets (Fowles and Burley, 1994; Fig. 3). It is notable that Fowles and Burley (1994) detected a further lower-angled, eastward-dipping fault orientation set within the caves; this orientation set could not be discerned within our data set. This discrepancy may potentially have arisen due to the fact that our survey only covered a single cavern within the cave system.

4.5. Shear-Band Connectivity and Termination Styles

Summation of the number of vertices contained within each label group of the trace connectivity map (Fig. 12A) reveals that the connectivity of the fault network exposed within the study site is remarkably high, with over 98% of all trace vertices falling within a single group (a 3D discontinuity trace connectivity map from the Lacy’s Caves outcrop is made available as part of the Supplementary Files). The high degree of connectivity of the fault network is reflected further in the distribution of observable terminations, which predominantly form abutting relationships with other traces (Fig. 12B). It should also be noted that the trace network is subject to a high degree of censoring, with 29.3% and 33.0% of traces from Set 1 and Set 2 being partially or fully curtailed due to the finite window of observation offered by the outcrop.

4.6. Shear-Band Trace Length

Representing a more intrinsic measurement of fracture size in comparison to the length of the outcrop-discontinuity curve of intersection, we use the chord distance between the end nodes of each trace to serve as a measure of fault trace length in this study. The mean trace lengths of individual cataclastic shear-band segments are 0.77 m for Set 1 and 0.70 m for Set 2 (Table 1). However, these results should be treated with caution due to the large proportion of censored trace lengths reported for this data set (see the previous section). Testing per–orientation-set trace-length distributions for goodness-of-fit to parametric models commonly used to describe fracture size in outcrop (lognormal, exponential, power law, and Weibull), we find that trace-length distributions are lognormal at the 5% significance level (Set 1: p = 0.16; Set 2: p = 0.10; Fig. 13) using Anderson-Darling goodness-of-fit criteria (Anderson and Darling, 1954). Rather than reflecting the size distribution of the fault array within the rock volume, it is widely recognized that lognormal behavior in fracture trace lengths measured at the exposure face is a product of the low probability of intersection between smaller fractures and the outcrop surface (Priest and Hudson, 1981).

4.7. Shear-Band Intensity

Vertex attribute maps of fracture intensity were computed using a spherical radius of 0.25 m, which was selected on the basis of being ∼10% of the cavern wall’s edge length. The areal fault intensity vertex attribute maps generated for each orientation set indicate that fault abundance varies over several orders of magnitude within the Lacy’s Caves study area (Figs. 14A and 14B). Based upon the mean values of P21 calculated over the outcrop surface, we find that NW-SE–oriented shear bands (Set 1) make the dominant contribution to faulting intensity within the study area, with the NE-SW–oriented set (Set 2) making a subsidiary contribution. There is nearly an order of magnitude difference in areal faulting intensity between the two identified orientation sets (Set 1 mean P21: 41.14 m/m2; Set 2 mean P21: 4.75 m/m2; Table 1), reflecting the close proximity of study area to the inferred parent fault structure associated with Set 1 (i.e., the Kirkoswald fault).

Aggregated values of P21 and volumetric faulting intensity (P32) calculated for Set 1 and Set 2 reflect the remarkably high degree of fault intensity present within the study area (Set 1 and Set 2 mean P21: 45.89 m/m2; Set 1 and Set 2 mean P32: 39.32 m2/m3). Furthermore, aggregated P21 and P32 surface maps reflect magnitudinal shifts in fault abundance at meter scale (P21/P32 vertex attribute maps from the Lacy’s Caves outcrops are made available as part of the Supplementary Files), with the highest values located around dense zones of shear-band clusters (see Figs. 14C and 14D). It should be noted that due to the complex geometry of the outcrop, values of P32 within the corners of the model will prove to be unreliable. This is because the directional cosine representing the orientation of the outcrop needed to calculate the Wang conversion factor will form the average between the adjoining walls and ceiling of the outcrop model contained within the spherical kernel (see Wang, 2005).


With the advent of open-source and cloud-based photogrammetric image-processing pipelines (Chandler and Fryer, 2013; Wu, 2011) and low-cost, portable sensory platforms (e.g., multi-rotor unmanned aerial vehicles [UAVs]; Carrivick et al., 2013), digital outcrop modeling is rapidly becoming a routine tool within the geosciences. Arguably, the development of approaches to characterize structure from digital outcrops has failed to keep pace with these advances in range data acquisition, with the focus still being emphatically toward mesh- segmentation–based strategies. In this study, we have demonstrated that synergy between conventional 2D image-processing techniques and mesh-based feature extraction bridges the gap between traditional structural surveys and current strategies for digital discontinuity analysis. By generating 3D trace maps from calibrated image sequences, we are able to access fracture and/or fault network properties that traditionally require the use of inefficient manual surveys (connectivity, trace length, termination styles, and fracturing intensity), while retaining the expediency of data acquisition and broad coverage of the sampling domain common to existing digital-based methods.

The derivation of static fracture properties from rock outcrops is pertinent to a range of problems within the geosciences, with these data sets having considerable utility within the fields of structural geology, rock mechanics, petroleum geoscience, and hydrogeology. It may be surmised, however, that obtaining a digitized representation of the exposed trace network, which places constraints upon the (apparent) size, orientation, and location characteristics of individual fractures, leads the way for further applications, providing the fundamental conditioning data required to create semi-deterministic realizations of the discontinuity array. We suggest that the development of stochastic fracture network models constrained by outcrop data may have significant application to geomechanical and fluid flow problems (e.g., Gillespie et al., 2011; Seers and Hodgetts, 2013), providing a tangible link between naturally fractured outcrop analogues and their subsurface counterparts. The development of such a modeling framework is the ongoing focus of the authors.

The work is funded by the UK Natural Environment Research Council (NERC). The authors also acknowledge Total E&P UK, ExxonMobil, and the American Association of Petroleum Geologists for their financial support of the project. M. James is thanked for providing the Trowbarrow Quarry photogrammetry data set used within this study. W. Head is thanked for his assistance in collecting the Lacy’s Caves photogrammetry data set. Finally, the authors would like to thank reviewers Max Wilkinson and Steven Micklewaite whose suggestions provided significant improvements to the manuscript.

1Supplemental Files. Stanford polygon format (.ply) files demonstrating different stages of the trace map feature extraction and analysis procedure detailed within the manuscript (Lacy’s Caves, Vale of Eden, UK). Please visit http://dx.doi.org/10.1130/GES01276.S1 or the full-text article on www.gsapubs.org to view the Supplemental Files.
Gold Open Access: This paper is published under the terms of the CC-BY license.