Extraction of three-dimensional fracture trace ... - GeoScienceWorld

0 downloads 111 Views 42MB Size Report
Jun 16, 2016 - east, and west were marked on the outcrop with the aid of a laser ...... and connectivity of joint system
Research Paper

GEOSPHERE GEOSPHERE; v. 12, no. 4

Extraction of three-dimensional fracture trace maps from calibrated image sequences Thomas D. Seers1,2 and David Hodgetts1 1

doi:10.1130/GES01276.1 14 figures; 2 interactive models; 1 table; 1 ­supplemental file CORRESPONDENCE:  thomas​.seers@​qatar​.tamu​ .edu CITATION: Seers, T.D., and Hodgetts, D., 2016, Extraction of three-dimensional fracture trace maps from calibrated image sequences: Geosphere, v. 12, no. 4, p. 1323–1340, doi:10.1130/GES01276.1. Received 28 September 2015 Revision received 27 February 2016 Accepted 6 May 2016 Published online 16 June 2016

OL D G

2

Department of Petroleum Engineering, Texas A&M University at Qatar, Qatar School of Earth, Atmospheric and Environmental Sciences, University of Manchester, Manchester M13 9PL, UK

ABSTRACT 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 meshand 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 pre­sent 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 varia­ bility 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.

1. INTRODUCTION OPEN ACCESS

This paper is published under the terms of the CC‑BY license.

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

petro­leum 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 out-

© 2016 The Authors

GEOSPHERE  |  Volume 12  |  Number 4

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1323

Research Paper crop surface may not be diagnostic of the presence of fractures. Many dis­conti­ nuity 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 deriva­ tive 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 recog­nized 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 indi­vidual 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

GEOSPHERE  |  Volume 12  |  Number 4

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 diag­ nostic 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 threedimen­sional 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

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1324

Research Paper mini­mal 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.

2. CALIBRATED IMAGES 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, [xw yw zw ]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):

m = K[R t]M ≡ PM, (1)

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



 r1,1 r1,2 r1,3 t x    [R t] =  r2,1 r2,2 r2,3 t y  . (2) r   3,1 r 3,2 r 3,3 t z 

C is the camera intrinsic matrix, containing the internal parameters of the camera:



 αx C =  0   0

γ αy 0

u0  v 0  , (3)  1

where g 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 ax and ay are products of Euclidean focal length f and scaling factors mx and my (ax = f · mx , ay = f · my ), which relate f to distance in pixels (Hartley and Zisserman,

GEOSPHERE  |  Volume 12  |  Number 4

Figure 1. Components of the pinhole camera model. The transform [R|t] maps the world coordinate (M) = [xw yw zw ]T to its equivalent point on the image plane m = [u v]T.

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):

C = R−1t ≡ −Rt t. (4)

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

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1325

Research Paper 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 (multiview 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.

3. THREE-DIMENSIONAL TRACE-MAP EXTRACTION AND ANALYSIS 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

GEOSPHERE  |  Volume 12  |  Number 4

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 compre­ hensively 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,

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1326

Research Paper

A

C

Trowbarrow Quarry

Phase congruency image 0

2.5 m

B

image intensity

1

2.5 m

Non-local means filtered image 0

phase congruency

1

D

Unsupervised image edges

2.5 m

2.5 m

Figure 2. Pixel-based feature delineation (Trowbarrow Quarry, Northwest UK; photogrammetry data set provided by M. James). (A) Input image. (B) The image is de-noised using a non-local means filter. (C) Phase congruency image edge detection. (D) Putative edges extracted from a phase congruency image of the quarry using hysteresis thresholding and binary skeletonization.

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 con-

GEOSPHERE  |  Volume 12  |  Number 4

nected 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

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1327

Research Paper

A

Unsupervised curves

C

Connected components

2.5 m

B

Unsupervised and manually fitted curves

2.5 m

D

E

2.5 m

Extracted 3D trace map: manually fi fitted and unsupervis unsupervised curves

4m

Figure 3. (A) Unsupervised curves are fitted using line-segment extraction and growing. (B) Trace objects mapped using unsupervised curve extraction are edited, merged, and supplemented with manually delineated trace features. (C) Connected two-dimensional (2D) trace components are identified using region labeling. Different colors indicate distinct label groups. (D) Extracted three-dimensional (3D) trace map from the Trowbarrow Quarry outcrop.

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 ema­ nat­ing from the origin of the world coordinate system, representing the opti­ cal 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

GEOSPHERE  |  Volume 12  |  Number 4

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 ˆ (where ωℝ ˆ ∈ℝ 3 and camera is achieved by locating the rotational axis ωℝ S(wx , wy, wz ) = 1) by finding the cross product between the vectors, a = [0 0 1]T, and the directional vector, b, obtained from the digitized trace (ωˆ = a × b). The matrix, Rωˆ , describing the rotation between a and b around ωˆ by the angle, h, is then obtained by solving the Rodrigues rotational formula (e.g., Murray et al., 1994):  cos η + ω 2x (1− cos η) ω x ωy (1− cos η) − ω z sinη ωy sinη + ω x ωy (1− cos η)    Rωˆ (η) = ω z sinη + ω x ωy (1− cos η) cos η + ωy2 (1− cos η) −ω x sinη + ωy ω z (1− cos η)   −ω sinη ω ω (1− cos η) ω sinη ω ω (1− cos η) cos η + ω 2z (1− cos η) . (5) x z x y z  y 

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1328

Research Paper

A

Optical Center

B

Camera Position C

R

Ray Tree

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 desira­ ble 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 poten­tial 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,

GEOSPHERE  |  Volume 12  |  Number 4

Masked Pixels

Figure 4. (A) Ray tree generated in nadir view. (B) Transformed ray tree with non-highlighted pixels masked. Individual curves are indexed.

Indexed Pixels

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

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1329

Research Paper

A

B

0.5 m

C

Fracture trace curve vertices

0.5 m

D

Nearest neighbor = end nodes

Non-overlapping

Figure 5. (A) Pixel resolution renders of overlapping images aid the identification of matched trace segments. (B) Merged trace segments displayed on a textured mesh. (C) Several potential overlap configurations complicate the merger of three-dimensional curves. (D) If curves are non-overlapping, nearest-neighboring end nodes are linked (smoothed using spline interpolation). In the overlapping case, vertex nearest-neighbor relationships are harnessed to locate overlapping portions of the curves to be merged.

Non-overlapping Nearest neighbor stabilization = point of overlap

Edge overlap

Center overlap

Nearest neighbor stabilization

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 tracemap 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

GEOSPHERE  |  Volume 12  |  Number 4

Overlapping

Point of overlap

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

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1330

Research Paper 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 (l1, l2, l3 ) of covariance matrix of each trace feature’s point set as a measure of vertex coplanarity: M = ln(l1 /l3 ), and collinearity: K = ln(l1 /l2 )/ln(l2 /l3 ) (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 ap-

Outcrop normal

Cluster centroid vMF c. parameter

N

Set 1 κ = 43 Dip = 71.1 Dip Az. = 39.3 n = 150

θ

f(θ | μ_

θ,k)

θ

B

Set 1 κ = 16 Dip = 13.4 Dip Az. = 306.4 n = 150

Outcrop best fit plane

C

N

45

29 26

8 1 5 4 17 318 2 97

Eigensolution

GEOSPHERE  |  Volume 12  |  Number 4

Trace collinearity (K)

28

20 34 30 42 37

11

vMF Distributed Fractures N

D

28

20 34 30 25 33 37 51 41 38 49 27 45

29 26

8 1 5 4 10 17 318 132 199 7

0.5

Random n = 150

Probabilistic

Seers and Hodgetts  |  Three-dimensional trace-map extraction

Trace collinearity (K)

A

plied 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 digi­ tized 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 (b) 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.

Figure 6. (A) Factors determining the proba­ bility density characteristics of a point, θ, on the 2-sphere around the arbitrary axis of rotation set by the end nodes of a threedimen­sional structural lineament (see text for details). (B) Two approximately orthonormal synthetic fracture sets drawn from the von Mises Fisher (vMF) distribution. (C) The analytical (eigensolution) to orientations of traces extracted through fractional Brownian surface intersections, the planes of which correspond to a subset of the vectors displayed in (B). (D) Equivalent prob­ abilistic solution improves the precision of the orientation estimates, particularly for traces with high collinearity values (K).

0.5

1331

Research Paper 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

GEOSPHERE  |  Volume 12  |  Number 4

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 ­ pper crust controls permeability and mechanical competence (Dershowitz u 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

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1332

Research Paper

Intersection length = P21 3D Trace map

Outcrop surface normal

Outcrop surface

θ Spherical kernel

Orientation set mean vector

Mesh vertex

Figure 7. A measure of P21 (areal intensity) is calculated at each point on the mesh surface by locating the intersection length of the trace map within a surrounding spherical kernel. A localized measure of P32 (volumetric intensity) is obtained using the conversion factor of Wang (2005), which requires that the angle θ between the outcrop normal and set vector is constrained.

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.

4. APPLICATION TO SHEAR-BAND CHARACTERIZATION 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 faultbounded, north-south–trending sedimentary basin located within Northwest England (Fig. 8A). The basin comprises a mixed sequence of predominantly

GEOSPHERE  |  Volume 12  |  Number 4

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 milli­meterscale 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, GUIbased application for 3D reconstruction using structure from motion estimation (Wu, 2011). The sparse point cloud and camera intrinsic-extrinsic parame­ ters generated by this procedure were then used as the basis for dense scene reconstruction using the patch-based, multi-view stereo tool (PMVS2) of Furu­ kawa 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 sur-

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1333

Research Paper

Arm

sth

te Dyke

Eden Shales

c ills

in nn Pe

hm atc Scr

Vale of Eden Basin

Sherwood Sst.

wai

e

ar Fa u

Fa

ult

lt

ill ult Fa

Carbonif.

tem Sys

Riv

er Ede n

A A

0

50 km

B

PENRITH

Penrith Sst.

Fault

nh ide Ma

ald

sw

ko

Kir

LACY’S CAVES

Faults

o 0

2 km

Figure 8. (A) Onshore and offshore outcrop of Permo-Triassic rocks in central and northern England and Wales, with the location of the Vale of Eden Basin indicated (redrawn from Fowles and Burley [1994]). (B) The central Vale of Eden Basin, with major fault trends and the present study area (Lacy’s Caves) indicated (adapted from Mountney and Thompson [2002]).

face 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). Supporting materials for Seers, T. D. and Hodgetts, D. Extraction of threedimensional fracture trace maps from calibrated image sequences The supporting materials consists of three 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): 1) Lacys_mesh_colour_unreg.ply - This is the unregistered coloured triangulated mesh generated from the raw point cloud data output by PMVS2 using Poisson surface reconstruction (see Section 4). This surface reconstruction can be viewed in the open source mesh viewer Meshlab (http://meshlab.sourceforge.net/) 2) Lacys_connectivity_map.ply -3D trace map displaying connectivity attribute. Note that this model has been scaled, rotated and translated into a local coordinate frame (i.e. it is true scale and orientation but not geolocated). This model can also be viewed using Meshlab. 3) Lacys_fault_intensity_all_sets.ply - P21 (areal intensity) and P32 (volumetric intensity) for all sets identified in Lacy's Caves displayed on a decimated mesh. P21 and P32 attributes are displayed in Fig. 14c&d of the manuscript. The base surface model can be also be viewed using Meshlab, though to view the P21/P32 attribute maps, users will need to download Virtual Reality Geological Studio (VRGS) - a multo modality mesh and image analysis software developed at the University of Manchester by Dr David Hodgetts (free to academia: http://www.click2go.umip.com/i/software/vrgs.html) For more information or support using these materials please contact: [email protected]

Supplemental 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.

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).

1

GEOSPHERE  |  Volume 12  |  Number 4

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 m = 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 m = 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.

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1334

Research Paper

A

B

View angle Fig. 9C

0.5 m

0.5 m

C Cataclastic shear bands

See Interactive Model 1 for a fully manipulable version of the Lacy’s Caves outcrop model.

0.5 m Figure 9. (A–C) Different views of a textured triangular irregular mesh-based surface model of the Lacy’s Caves outcrop constructed from a structure from motion estimation–multi-view stereo (SfM-MVS)–derived point cloud using the Poisson surface reconstruction algorithm. See Interactive Model 1 for a fully manipulable version of the Lacy’s Caves outcrop model.

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

GEOSPHERE  |  Volume 12  |  Number 4

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.

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1335

Research Paper

A

B Figure 10. (A) Curves (n = 802) extracted from the digital surface model shown in Figure 9. (B) The raw curve data set is merged to form a three-dimensional trace map composed of 645 trace objects. See Interactive Model 2 for a fully manipulable version of the Lacy’s Caves 3D trace map.

See Interactive Model 2 for a fully manipulable version of the Lacy’s Caves 3D trace map.

00.75 75 m

00.75 75 m N

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 goodnessof-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

GEOSPHERE  |  Volume 12  |  Number 4

Set 1: n = 508

Set 2: n = 66

Dip = 85.3 Dip Az. = 136.6

Dip = 87.6 Dip Az. = 231.9

Figure 11. Fault orientation sets identified using the K-means method and a loss function to determine the non-systematic component (displayed in blue) for cataclastic shear-band traces extracted from the Lacy’s Caves outcrop. Two major fault sets are identified from the cluster analysis: Set 1: n = 508; direction of μ = 085/137; k = 9.8; Set 2: n = 66; direction of μ = 088/232; k = 13.4. Blue measurements are classified as non-systematic by the cost minimization method.

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1336

Research Paper TABLE 1. SUMMARY STATISTICS FOR CATACLASTIC SHEAR BANDS MAPPED FROM LACY’S CAVES*

Set 1 Set 2

n

*µdip

µazimuth

Fisher’s k†

µlength (m)

σlength (m)

tmin (m)

tmax (m)

µP21 (m/m2)

µP32 (m2/m3)

508 66

85.3 87.6

136.6 231.9

9.8 13.4

0.77 0.70

0.49 0.54

0.10 0.09

2.60 3.02

41.14 4.75

35.56 3.76

*µ indicates mean values; σ indicates standard deviation; tmin and tmax are the minimum and maximum trace lengths. Fisher (1953).



A

BT

60%

20%

100 % 100

80%

60%

40%

20%

% T Type

Figure 12. (A) Three-dimensional trace map displayed with global connectivity relations indicated: vertices within each label group are displayed using the same color. Note that the trace network is almost entirely connected. (B) Styles of observable trace terminations. See Figure 11 for the color classification.

μlength = 0.70

σ length = 0.49

60

40%

80%

0.75 0 75 m

12

μlength = 0.77

e yp

60%

~90% abutting

B

XT

40%

yp e

80%

σ length = 0.54

AD Test: Lognormal AD = 0.54 p = 0.15

40

AD Test: Lognormal AD = 0.63 p = 0.10

8

n faults

100%

:

n faults

B

20%

A

4

20

0 0

1

2

3

0

Trace length (m)

0

1

2

3

Trace length (m)

Figure 13. Trace-length distributions for (A) Set 1 and (B) Set 2. Anderson-Darling (AD) test statistics for goodness-of-fit to the lognormal distribution are also displayed (α = 0.05).

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).

GEOSPHERE  |  Volume 12  |  Number 4

5. DISCUSSION AND CONCLUSIONS 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.

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1337

Research Paper

A

P21

176

Set 1 P21

1m

C

m/m²

P21

0

185

1m

m/m²

P32

D

0

171

All Sets

P21

Figure 14. Vertex attribute maps displaying values of cataclastic shear band P21 (­ areal fault intensity: m/m2) calculated across the exposure face. The map was produced using a spherical kernel ­radius of 0.25m. (A) Set 1 P 21 . (B) Set 2 P 21 . (C) Set 1 and Set 2 P21. (D) Set 1 and Set 2 P 32 derived from values of P21 using the conversion factor of Wang (2005).

P32

m/m²

0

0

1m

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, petro­leum 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.

ACKNOWLEDGMENTS 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

GEOSPHERE  |  Volume 12  |  Number 4

33

Set 2 P21

All Sets

1m

P21

B

m²/m³

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. REFERENCES CITED Anderson, T.W., and Darling, D.A., 1954, A test of goodness of fit: Journal of the American Statistical Association, v. 49, no. 268, p. 765–769, doi:​10​.1080​/01621459​.1954​.10501232​. Aydin, A., 2000, Fractures, faults, and hydrocarbon entrapment, migration and flow: Marine and Petroleum Geology, v. 17, no. 7, p. 797–814, doi:​10​.1016​/S0264​-8172​(00)00020​-9​. Baecher, G., Lanney, N., and Einstein, H., 1977, Statistical description of rock properties and sampling: Proceedings of the 18th U.S. Symposium on Rock Mechanics: American Rock Mechanics Association, no. 5C1-8. Barton, C.A., Zoback, M.D., and Moos, D., 1995, Fluid flow along potentially active faults in crystalline rock: Geology, v. 23, no. 8, p. 683–686, doi:​10​.1130​/0091​-7613​(1995)023​2​.3​.CO;2​. Buades, A., Coll, B., and Morel, J.-M., 2005, A non-local algorithm for image denoising: IEEE Computer Society Conference on Computer Vision and Pattern Recognition, v. 2, p. 60–65. Canny, J., 1986, A computational approach to edge detection: IEEE Transactions on Pattern Analy­ sis and Machine Intelligence, no. 6, p. 679–698. Carrivick, J.L., Smith, M.W., Quincey, D.J., and Carver, S.J., 2013, Developments in budget remote sensing for the geosciences: Geology Today, v. 29, no. 4, p. 138–143, doi:​10​.1111​/gto​.12015​. Casini, G., Gillespie, P., Vergés, J., Romaire, I., Fernández, N., Casciello, E., Saura, E., Mehl, C., Homke, S., and Embry, J.-C., 2011, Sub-seismic fractures in foreland fold and thrust belts: insight from the Lurestan Province, Zagros Mountains, Iran: Petroleum Geoscience, v. 17, no.  3, p.  263–282, doi:​10​.1144​/1354​-079310​-043​.

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1338

Research Paper Chandler, J., and Fryer, J., 2013, Autodesk 123D catch: How accurate is it?: Geomatics World, v. 2, no. 21, p. 28–30. Clemen, R.T., and Winkler, R.L., 1999, Combining probability distributions from experts in risk analysis: Risk Analysis, v.  19, no.  2, p.  187–203, doi:​10​.1111​/j​.1539​-6924​.1999​.tb00399​.x​. Cook, R.L., Porter, T., and Carpenter, L., 1984, Distributed ray tracing: Computer Graphics, v. 18, no. 3, p. 137–145. Debevec, P.E., Taylor, C.J., and Malik, J., 1996, Modeling and rendering architecture from photo­ graphs: A hybrid geometry-and image-based approach: Computer Graphics Proceedings, ACM SIGGRAPH 96, p. 11–20. Dershowitz, W.S., and Herda, H.H., 1992, Interpretation of fracture spacing and intensity: The 33th U.S. Symposium on Rock Mechanics (USRMS), 3–5 June, Santa Fe, New Mexico: Rotter­dam, A.A. Balkema, p. 757–766. Dershowitz, B., LaPointe, P., Eiben, T., and Wei, L., 2000, Integration of discrete feature network methods with conventional simulator approaches: Society of Petroleum Engineers, Reservoir Evaluation & Engineering, v. 3, no. 2, p. 165–170, doi:​10​.2118​/62498​-PA​. Di Ruberto, C., 2004, Recognition of shapes by attributed skeletal graphs: Pattern Recognition, v.  37, no.  1, p.  21–31, doi:​10​.1016​/j​.patcog​.2003​.07​.004​. Dyer, R., 1988, Using joint interactions to estimate paleostress ratios: Journal of Structural Geology, v.  10, no.  7, p.  685–699, doi:​10​.1016​/0191​-8141​(88)90076​-4​. Faugeras, O.D., Luong, Q.-T., and Maybank, S.J., 1992, Camera self-calibration: Theory and experiments, in Sandini, G., eds., Proceedings Computer Vision—ECCV’92 1992: Berlin, Heidel­ berg, Springer, p. 321–334. Favalli, M., Fornaciai, A., Isola, I., Tarquini, S., and Nannipieri, L., 2012, Multiview 3D reconstruction in geosciences: Computers & Geosciences, v. 44, p. 168–176, doi:​10​.1016​/j​.cageo​.2011​ .09​.012​. Fernández, O., 2005, Obtaining a best fitting plane through 3D georeferenced data: Journal of Structural Geology, v.  27, no.  5, p.  855–858, doi:​10​.1016​/j​.jsg​.2004​.12​.004​. Fischler, M.A., and Bolles, R.C., 1981, Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography: Communications of the Association for Computing Machinery, v. 24, no. 6, p. 381–395, doi:​10​.1145​/358669​.358692​. Fisher, Q., and Knipe, R., 2001, The permeability of faults within siliciclastic petroleum reservoirs of the North Sea and Norwegian Continental Shelf: Marine and Petroleum Geology, v. 18, no.  10, p.  1063–1081, doi:​10​.1016​/S0264​-8172​(01)00042​-3​. Fisher, R., 1953, Dispersion on a sphere: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 1953, Volume 217: The Royal Society, p. 295–305. Fitzgibbon, A.W., 2001, Simultaneous linear estimation of multiple view geometry and lens distortion: IEEE Computer Society Conference on Computer Vision and Pattern Recognition, v. 121, p. I-125–I-132. Florez-Niño, J.-M., Aydin, A., Mavko, G., Antonellini, M., and Ayaviri, A., 2005, Fault and fracture systems in a fold and thrust belt: An example from Bolivia: American Association of Petroleum Geologists Bulletin, v. 89, no. 4, p. 471–493, doi:​10​.1306​/11120404032​. Fowles, J., and Burley, S., 1994, Textural and permeability characteristics of faulted, high porosity sandstones: Marine and Petroleum Geology, v. 11, no. 5, p. 608–623, doi:​10​.1016​/0264​ -8172​(94)90071​-X​. Fujimoto, A., Tanaka, T., and Iwata, K., 1986, Arts: Accelerated ray-tracing system: IEEE, Computer Graphics and Applications, v. 6, no. 4, p. 16–26, doi:​10​.1109​/MCG​.1986​.276715​. Furukawa, Y., and Ponce, J., 2010, Accurate, dense, and robust multiview stereopsis: Pattern Analysis and Machine Intelligence: IEEE Transactions on Pattern Analysis and Machine Intelligence, v. 32, no. 8, p. 1362–1376. García-Sellés, D., Falivene, O., Arbués, P., Gratacos, O., Tavani, S., and Muñoz, J.A., 2011, Supervised identification and reconstruction of near-planar geological surfaces from terrestrial laser scanning: Computers & Geosciences, v. 37, no. 10, p. 1584–1594, doi:​10​.1016​/j​.cageo​ .2011​.03​.007​. Ghosh, K., and Mitra, S., 2009, Structural controls of fracture orientations, intensity, and connectivity, Teton anticline, Sawtooth Range, Montana: American Association of Petroleum Geologists Bulletin, v. 93, no. 8, p. 995–1014, doi:​10​.1306​/04020908115​. Gillespie, P., Walsh, J., Watterson, J., Bonson, C., and Manzocchi, T., 2001, Scaling relationships of joint and vein arrays from The Burren, Co. Clare, Ireland: Journal of Structural Geology, v.  23, no.  2, p.  183–201, doi:​10​.1016​/S0191​-8141​(00)00090​-0​. Gillespie, P., Monsen, E., Maerten, L., Hunt, D., Thurmond, J., and Tuck, D., 2011, Fractures in carbonates: from digital outcrops to mechanical models: Outcrops revitalized: Tools, techniques and applications: Tulsa, SEPM (Society for Sedimentary Geology), p. 137–148.

GEOSPHERE  |  Volume 12  |  Number 4

Grasmueck, M., 1996, 3-D ground-penetrating radar applied to fracture imaging in gneiss: Geophysics, v. 61, no. 4, p. 1050–1064, doi:​10​.1190​/1​.1444026​. Haralick, R., and Shapiro, L., 1992, Robot and Computer Vision (vols. 1 and 2): Reading, Massachusetts, Addison-Wesley. Hartley, R., and Zisserman, A., 2003, Multiple view geometry in computer vision: Cambridge, UK, Cambridge University Press, 657 p. Hatton, C., Main, I., and Meredith, P., 1994, Non-universal scaling of fracture length and opening displacement: Nature, v. 367, no. 6459, p. 160–162, doi:​10​.1038​/367160a0​. Hill, R., 1963, Elastic properties of reinforced solids: some theoretical principles: Journal of the Mechanics and Physics of Solids, v. 11, p. 357–372, doi:​10​.1016​/0022​-5096​(63)90036​-X​. Hitchmough, A.M., Riley, M.S., Herbert, A.W., and Tellam, J.H., 2007, Estimating the hydraulic properties of the fracture network in a sandstone aquifer: Journal of Contaminant Hydrology, v.  93, no.  1, p.  38–57, doi:​10​.1016​/j​.jconhyd​.2007​.01​.012​. Holliday, D., 1993, Mesozoic cover over northern England: Interpretation of apatite fission track data: Journal of the Geological Society of London, v. 150, no. 4, p. 657–660, doi:​10​.1144​ /gsjgs​.150​.4​.0657​. Horn, B.K., 1987, Closed-form solution of absolute orientation using unit quaternions: Journal of the Optical Society of America A, v. 4, no. 4, p. 629–642, doi:​10​.1364​/JOSAA​.4​.000629​. Hua, W., Zhang, H., Lu, Y., Bao, H., and Peng, Q., 2004, Huge texture mapping for real-time visuali­zation of large-scale terrain: Proceedings of the Association for Computing Machinery Symposium on Virtual Reality Software and Technology 2004: Association for Computing Machinery, p. 154–157. Kazhdan, M., Bolitho, M., and Hoppe, H., 2006, Poisson surface reconstruction: Proceedings of the Fourth Eurographics Symposium on Geometry Processing 2006, Volume 7, p. 61–70. Kemeny, J., and Post, R., 2003, Estimating three-dimensional rock discontinuity orientation from digital images of fracture traces: Computers & Geosciences, v. 29, no. 1, p. 65–77, doi:​10​.1016​ /S0098​-3004​(02)00106​-1​. Kemeny, J., Turner, K., and Norton, B., 2006, LIDAR for rock mass characterization: Hardware, software, accuracy and best-practices, in Tonon F., and Kottenstette, J., eds., Laser and Photogrammetric Methods for Rock Face Characterization: Alexandria, Virginia, American Rock Mechanics Association, p. 49–62. Kim, G.-B., Lee, J.-Y., and Lee, K.-K., 2004, Construction of lineament maps related to ground­ water occurrence with ArcView and Avenue™ scripts: Computers & Geosciences, v. 30, no. 9, p.  1117–1126, doi:​10​.1016​/j​.cageo​.2004​.09​.002​. Kovesi, P., 1999, Image features from phase congruency: Videre: Journal of computer vision research, v. 1, no. 3, p. 1–26. La Pointe, P., 2002, Derivation of parent fracture population statistics from trace length measurements of fractal fracture populations: International Journal of Rock Mechanics and Mining Sciences, v.  39, no.  3, p.  381–388, doi:​10​.1016​/S1365​-1609​(02)00021​-7​. Le Garzic, E., de L’Hamaide, T., Diraison, M., Géraud, Y., Sausse, J., De Urreiztieta, M., Hauville, B., and Champanhet, J.-M., 2011, Scaling and geometric properties of extensional fracture systems in the Proterozoic basement of Yemen: Tectonic interpretation and fluid flow implications: Journal of Structural Geology, v. 33, no. 4, p. 519–536, doi:​10​.1016​/j​.jsg​.2011​.01​.012​. Lemy, F., and Hadjigeorgiou, J., 2003, Discontinuity trace map construction using photographs of rock exposures: International Journal of Rock Mechanics and Mining Sciences, v. 40, no. 6, p.  903–917, doi:​10​.1016​/S1365​-1609​(03)00069​-8​. Lowe, D.G., 1999, Object recognition from local scale-invariant features: Proceedings of the IEEE International Conference on Computer Vision, v. 2, p. 1150–1157. Mauldon, M., Dunne, W., and Rohrbaugh, M., 2001, Circular scanlines and circular windows: New tools for characterizing the geometry of fracture traces: Journal of Structural Geology, v.  23, no.  2, p.  247–258, doi:​10​.1016​/S0191​-8141​(00)00094​-8​. Meagher, D., 1982, Geometric modeling using octree encoding: Computer Graphics and Image Processing, v.  19, no.  2, p.  129–147, doi:​10​.1016​/0146​-664X​(82)90104​-6​. Möller, T., and Trumbore, B., 1997, Fast, minimum storage ray-triangle intersection: Journal of Graphics Tools, v. 2, no. 1, p. 21–28, doi:​10​.1080​/10867651​.1997​.10487468​. Mountney, N.P., and Thompson, D.B., 2002, Stratigraphic evolution and preservation of aeolian dune and damp/wet interdune strata: An example from the Triassic Helsby Sandstone Formation, Cheshire Basin, UK: Sedimentology, v. 49, no. 4, p. 805–833, doi:​10​.1046​/j​.1365​-3091​ .2002​.00472​.x​. Murray, R.M., Li, Z., Sastry, S.S., and Sastry, S.S., 1994, A mathematical introduction to robotic manipulation: Boca Raton, CRC Press, 463 p.

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1339

Research Paper Neuman, S.P., 1988, Stochastic continuum representation of fractured rock permeability as an alternative to the REV and fracture network concepts, in Custodio, E., et al., eds., Groundwater Flow and Quality Modelling: Netherlands, Springer, p. 331–362, doi:​10​.1007​/978​-94​ -009​-2889​-3_19​. Nixon, M.S., and Aguado, A.S., 2012, Feature extraction and image processing for computer vision: Oxford, UK, Academic Press, 609 p. Odling, N., Gillespie, P., Bourgine, B., Castaing, C., Chiles, J., Christensen, N., Fillion, E., Genter, A., Olsen, C., and Thrane, L., 1999, Variations in fracture system geometry and their implications for fluid flow in fractured hydrocarbon reservoirs: Petroleum Geoscience, v. 5, no. 4, p.  373–384, doi:​10​.1144​/petgeo​.5​.4​.373​. Odling, N.E., 1997, Scaling and connectivity of joint systems in sandstones from western Norway: Journal of Structural Geology, v. 19, no. 10, p. 1257–1271, doi:​10​.1016​/S0191​-8141​(97)00041​-2​. Otoo, J., Maerz, N., Xiaoling, L., and Duan, Y., 2011, 3-D discontinuity orientations using combined optical imaging and LiDAR techniques: Proceedings of the 45th U.S. Rock Mechanics/­ Geomechanics Symposium 2011, American Rock Mechanics Association, https://​www​ .onepetro​.org​/conference​-paper​/ARMA​-11​-367. Pandey, G., McBride, J.R., Savarese, S., and Eustice, R., 2012, Automatic Targetless Extrinsic Cali­ bration of a 3D Lidar and Camera by Maximizing Mutual Information, in Proceedings of the Twenty-Sixth AAAI (Association for the Advancement of Artificial Intelligence): Conference on Artificial Intelligence, p. 2053–2059. Pearce, M.A., Jones, R.R., Smith, S.A., and McCaffrey, K.J., 2011, Quantification of fold curvature and fracturing using terrestrial laser scanning: American Association of Petroleum Geologists Bulletin, v. 95, no. 5, p. 771–794, doi:​10​.1306​/11051010026​. Priest, S., and Hudson, J., 1981, Estimation of discontinuity spacing and trace length using scanline surveys: Proceedings International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 1981, Volume 18, Elsevier, p. 183–197. Priest, S.D., 1993, Discontinuity analysis for rock engineering: Springer Science & Business M ­ edia, doi:​10​.1007​/978​-94​-011​-1498​-1​. Reid, T., and Harrison, J., 2000, A semi-automated methodology for discontinuity trace detection in digital images of rock mass exposures: International Journal of Rock Mechanics and Mining Sciences, v. 37, no. 7, p. 1073–1089, doi:​10​.1016​/S1365​-1609​(00)00041​-1​. Renshaw, C.E., 1997, Mechanical controls on the spatial density of opening-mode fracture networks: Geology, v. 25, no. 10, p. 923–926, doi:​10​.1130​/0091​-7613​(1997)025​2​ .3​.CO;2​. Rogers, S., Elmo, D., Webb, G., and Catalan, A., 2013, Volumetric fracture intensity measurement for improved rock mass characterisation and fragmentation assessment in block caving ­operations: Rock Mechanics and Rock Engineering, v. 48, p. 1–17. Roncella, R., and Forlani, G., 2005, Extraction of planar patches from point clouds to retrieve dip and dip direction of rock discontinuities: Proceedings of Laser Scanning, p. 162–167. Sahini, M., and Sahimi, M., 2003, Applications of Percolation Theory: London, CRC Press, 276 p. Scaramuzza, D., Harati, A., and Siegwart, R., 2007, Extrinsic self calibration of a camera and a 3D laser range finder from natural scenes: San Diego, California, IEEE/RSJ International Conference on Intelligent Robots and Systems, p. 4164–4169. Seers, T.D., and Hodgetts, D., 2013, Comparison of digital outcrop and conventional data collection approaches for the characterization of naturally fractured reservoir analogues: Geological Society of London, Special Publications, v. 374, p. 374–313. Seers, T.D., and Hodgetts, D., 2016, Probabilistic constraints upon structural lineament orientation constrained using numerical analysis: Journal of Structural Geology, v. 82, p. 37–47, doi:​ 10​.1016​/j​.jsg​.2015​.11​.004​.

GEOSPHERE  |  Volume 12  |  Number 4

Smits, B., 2002, Efficient bounding box intersection: Ray Tracing News, v. 15, no. 1. Snavely, N., Seitz, S.M., and Szeliski, R., 2006, Photo tourism: Exploring photo collections in 3D, in Proceedings ACM transactions on graphics (TOG), Volume 25, Association for Computing Machinery, p. 835–846. Sternlof, K.R., Karimi-Fard, M., Pollard, D., and Durlofsky, L., 2006, Flow and transport effects of compaction bands in sandstone at scales relevant to aquifer and reservoir management: Water Resources Research, v. 42, no. 7, W07425, doi:​10​.1029​/2005WR004664​. Stone, M., 1961, The opinion pool: Annals of Mathematical Statistics, v. 32, no. 4, p. 1339–1342, doi:​10​.1214​/aoms​/1177704873​. Sturzenegger, M., and Stead, D., 2009, Close-range terrestrial digital photogrammetry and terrestrial laser scanning for discontinuity characterization on rock cuts: Engineering Geology, v.  106, no.  3, p.  163–182, doi:​10​.1016​/j​.enggeo​.2009​.03​.004​. Terzaghi, R.D., 1965, Sources of error in joint surveys: Geotechnique, v. 15, no. 3, p. 287–304, doi:​ 10​.1680​/geot​.1965​.15​.3​.287​. Thiele, S.T., Micklethwaite, S., Bourke, P., Verrall, M., and Kovesi, P., 2015, Insights into the mechan­ics of en-échelon sigmoidal vein formation using ultra-high resolution photogrammetry and computed tomography: Journal of Structural Geology, v. 77, p. 27–44, doi:​10​.1016​ /j​.jsg​.2015​.05​.006​. Umili, G., Ferrero, A., and Einstein, H., 2013, A new method for automatic discontinuity traces sampling on rock mass 3D model: Computers & Geosciences, v. 51, p. 182–192, doi:​10​.1016​ /j​.cageo​.2012​.07​.026​. Underhill, J., Gayer, R., Woodcock, N., Donnelly, R., Jolley, E., and Stimpson, I., 1988, The Dent fault system, northern England—Reinterpreted as a major oblique-slip fault zone: Journal of the Geological Society of London, v. 145, no. 2, p. 303–316, doi:​10​.1144​/gsjgs​.145​.2​.0303​. Van Knapen, B., and Slob, S., 2006, Identification and characterisation of rock mass discontinuity sets using 3D laser scanning, in Culshaw, M.G., Reeves, H.J., Jefferson, I., and Spink, T.W., eds., Engineering Geology for Tomorrow’s Cities: Geological Society, London, Engineering Geology Special Publications, IAEG2006 Paper 438, CD-ROM insert, 10 p. Vasuki, Y., Holden, E.-J., Kovesi, P., and Micklethwaite, S., 2014, Semi-automatic mapping of geological structures using UAV-based photogrammetric data: An image analysis approach: Computers & Geosciences, v.  69, p.  22–32, doi:​10​.1016​/j​.cageo​.2014​.04​.012​. Wang, L., Kang, S.B., Szeliski, R., and Shum, H.-Y., 2001, Optimal texture map reconstruction from multiple views: IEEE Computer Society Conference on Computer Vision and Pattern Recognition, v.  341, p.  I-347–I-354, doi:​10​.1109​/CVPR​.2001​.990496​. Wang, X., 2005, Stereological interpretation of rock fracture traces on borehole walls and other cylindrical surfaces: Blacksburg, Virginia Polytechnic Institute and State University, 123 p. Woodcock, N., 1977, Specification of fabric shapes using an eigenvalue method: Geological Society of America Bulletin, v. 88, no. 9, p. 1231–1236, doi:​10​.1130​/0016​-7606​(1977)88​2​.0​.CO;2​. Wu, C., 2011, VisualSFM: A visual structure from motion system, http://​ccwu​.me​/vsfm/ (accessed 9 April 2016). Yamaguchi, F., 2012, Curves and surfaces in computer aided geometric design: Heidelberg, Springer Science & Business Media, 378 p. Zhang, L., 2014, Camera calibration, in Ikeuchi, K., eds., Computer Vision: A Reference Guide, p.  76–77, doi:​10​.1007​/978​-0​-387​-31439​-6_164​. Zhang, Z., 2000, A flexible new technique for camera calibration: IEEE Transactions on Pattern Analysis and Machine Intelligence, v. 22, no. 11, p. 1330–1334.

Seers and Hodgetts  |  Three-dimensional trace-map extraction

1340

Research Paper

Interactive Model 1: Lacy’s Caves outcrop model. Return to Figure 9.

GEOSPHERE

|

Volume 12

|

Number 4

Seers and Hodgetts

|

Three-dimensional trace-map extraction

Interactive Model 1

Research Paper

Interactive Model 2: Lacy’s Caves 3D trace map. Return to Figure 10.

GEOSPHERE

|

Volume 12

|

Number 4

Seers and Hodgetts

|

Three-dimensional trace-map extraction

Interactive Model 2