The world’s largest fossil oyster reef, formed by the giant oyster Crassostrea gryphoides and located in Stetten (north of Vienna, Austria), is studied in this article. Digital documentation of the unique geological site is provided by terrestrial laser scanning (TLS) at the millimeter scale. Obtaining meaningful results is not merely a matter of data acquisition with a suitable device; it requires proper planning, data management, and postprocessing. Terrestrial laser scanning technology has a high potential for providing precise 3D mapping that serves as the basis for automatic object detection in different scenarios; however, it faces challenges in the presence of large amounts of data and the irregular geometry of an oyster reef. We provide a detailed description of the techniques and strategy used for data collection and processing. The use of laser scanning provided the ability to measure surface points of 46,840 (estimated) shells. They are up to 60-cm-long oyster specimens, and their surfaces are modeled with a high accuracy of 1 mm. In addition, we propose an automatic analysis method for identifying and enumerating convex parts of shells. Object surfaces were detected with a completeness of 69% and a correctness of over 75% by means of a fully automated workflow. Accuracy of 98% was achieved in detecting the number of objects. In addition to laser scanning measurements, more than 300 photographs were captured, and an orthophoto mosaic was generated with a ground sampling distance (GSD) of 0.5 mm. This high-resolution 3D information and the photographic texture serve as the basis for ongoing and future geological and paleontological analyses. Moreover, they provide unprecedented documentation for conservation issues at a unique natural heritage site.

The fossil oyster reef in Stetten, Austria, is the world’s largest excavated fossil oyster reef, formed by large seashells. It is part of geo-edutainment park, “Fossilienwelt Weinviertel,” and is a protected site. About 50,000 shells up to 60 cm long cover a 459 m2 area (Harzhauser et al., 2015). The reef consists primarily of shells of Crassostrea gryphoides (Schlotheim, 1813), which is an ostreid bivalve with largely calcitic shells. Each individual developed two convex, strongly elongated valves of different sizes (Fig. 1). A ligament in the hinge area, together with a single adductor muscle located posteroventrally in the interior shell cavity, kept the valves articulated during the bivalve’s lifespan. After death, the valves usually became disarticulated.

The densely packed shell bed represents an event layer that was formed ∼16.5 million years ago in a tropical estuary (Sovis and Schmid, 1998, 2002; Latal et al., 2006; Harzhauser et al., 2009). This unique fossil accumulation was formed at the onset of the last global climate maximum, known as the Middle Miocene Climate Optimum (Zachos et al., 2001). The shell bed is preliminarily interpreted to have formed during a single hydrodynamic event, such as a tsunami or a major storm (Harzhauser et al., 2009, 2015). A detailed analysis is required to elucidate the formation of the structure. For instance, the number of shells, the ratio between left and right valves, their orientation, and the degree of fragmentation provide bases for further interpretation. Values for those parameters are described in Harzhauser et al. (2015) on the basis of a subset of the oyster reef. However, such data cannot be acquired on site for a larger area; it is impossible to step on the shell bed without destroying the fragile fossils. Moreover, the field measurements might be biased by subjectivity, increasing the difficulty of cross-checking the data afterwards. Finally, traditional onsite measurements consume a great deal of time. Non-destructive georeferenced data acquisition represents a major breakthrough for the analysis of such geological structures. The extent of the reef and the good condition of the fossils qualify the oyster reef for a detailed documentation for current and future scientific investigations. Thus, geospatial and georeferenced information needs to be acquired for modeling its objects. Therefore, one aspect of the study is the technology transfer between the disciplines of photogrammetry and paleontology. That requires testing already existing methods from photogrammetry and laser scanning and evaluating their suitability for this particular application.

The very irregular shell surface not only poses a problem in interpretation (manually as well as automatically) but also in aligning the orientation of the individual acquisitions relative to each other. In both methods (laser scanning point clouds and photogrammetric images), manually identified corresponding points needed to be identified. Neither an iterative closest point (ICP) algorithm nor an automatic bundle block adjustment, based on automatically extracted correspondences alone, provided acceptable results. Thus, surface reconstruction based on images was not further explored, due to the additional consideration that the main interest lies in the analysis of the oyster field rather than in a determination of the best performing acquisition method for a unique experiment. To minimize risk, the most reliable 3D acquisition method applicable was selected (see Acquisition Strategy section).

Documentation, archiving, and analysis of natural and man-made objects and scenes are ongoing endeavors. If the focus is on geometrical aspects, 3D point clouds (Otepka et al., 2013) are often the primary choice for the description of objects, along with raster digital elevation models (DEMs) and triangulated models (meshes). For smaller scenes with complex geometry, terrestrial laser scanning (TLS) is a well-established method to describe objects with centimeter precision for extended areas of up to 100 m2. However, it is challenging to achieve mm accuracy over a large site with irregular shapes. That level of accuracy requires not only one scan but, over a larger scene, coverage by several scans. The surfaces of the shells are very irregular, feature high curvature and roughness, and include overhangs. Therefore, studying fossil reefs from a geological point of view by analysis of their surface models requires millimeter resolution because the objects range in size from a few millimeters to a few decimeters. Respective approaches that provide very high resolution and precise 3D models are needed. In this paper, the technique comprises geometric, as well as photographic, documentation with a resolution of 1 mm and 0.5 mm, respectively. Furthermore, features such as measures of concavity and convexity can be extracted automatically, allowing information extraction and geological interpretation.

In summary, the contributions of this article are: (1) a report documenting the uniqueness of this geological site; (2) a strategy for capturing irregular surfaces by high-resolution point clouds; (3) an analysis method to detect convex parts of shells; and (4) a method to count and detect the number of shells automatically.

We considered the state of the art in surface modeling related to: geometrical shapes such as cylindrical, conical, or spiral surfaces; 3D surfaces in geology obtained by surveying technologies (photogrammetry and laser scanning); surface representations and different feature extraction techniques, primarily in geology; and studies focused on modeling small objects or structures similar to our fossilized objects, such as grains, river bed rocks, rough terrain, or relief with convex features. We also introduced TLS studies related to current living oyster and coral reefs. With different surface representations, accuracy is dependent on scale, occlusions due to acquisition geometry, point density, and ruggedness of the original surface. Surface modeling has been investigated for decades (Hoppe et al., 1992; Edelsbrunner and Mücke, 1994), taking sets of unordered 3D points as input and providing a simplified surface as output. Automatic recognition of objects with simple shapes, such as planes, was described by Schnabel et al. (2007) and Pottmann et al. (2004). Recognition of a wide variety of surfaces (spiral surfaces, cones, etc.) from 3D data was studied by Pottmann et al. (2002) and Hofer et al. (2005). In addition to those studies, the source of the geometric 3D point information was considered. For instance, detection of cylindrical, conical, or spherical surfaces from range images was described by Marshall et al. (2001) and Han et al. (2004). Data sources affect properties such as noise characteristics and occlusions. Shape analysis restricted to 2D in a photographic image was studied as well, for purposes such as extracting lines and curves based on the Hough transformation (Duda and Hart, 1972). In Heijmans (1994), Boykov and Jolly (2001), and Chan and Vese (2001), the analysis of scientific images to approximate edges of objects was studied. Determined pixels of “object” or “background” provide hard constraints for segmentation. Both sources, geometry and radiometry, can be used in a combined 3D analysis.

As described by Pesci and Teza (2008) and Burton et al. (2011), typical surfaces acquired in a geological survey are irregular. To identify features, or, more generally, for pattern recognition, intensity data (including spectral imaging), as in Kurz et al. (2013), can also be used as input to evaluate the reflectance and absorption properties of the material. Geological features in 3D, such as surfaces or geological bodies, are traditionally represented by polygons, polyhedrons, or voxels (Jessell, 2001). However, 3D triangulated surfaces, surface normals, and shaded relief models allow for a more precise interpretation (Samson and Mallet, 1997; Caumon et al., 2009). Detailed 3D information provides a virtual copy of the studied object and is more intuitive to interpret than traditional data from geological fields (Buckley et al., 2010). That information has demonstrated the use of laser scanning and photorealistic modeling for mapping and modeling the configuration of geological surfaces in exposed rock outcrops, as well as the ability to capture the geometry of subtle and small-scale features. Zakšek et al. (2011) present relief characteristics with hypsometric colors (graduated bands of color), elevation contours, and details by certain standard spatial analyses (e.g., slope, aspect, curvature, or local relief model), but the best visual impact is obtained with relief shading—a representation of the relief in a natural and intuitive manner. First approaches of modeling complex 3D geological surfaces were proposed by Mallet (1989) in the Paradigm GOCAD research project and by Bellian et al. (2005). Buckley et al. (2008) and Sankey et al. (2011) presented an overview of interpreting 3D data from TLS in a geological context. Measurements of surface roughness have examined elevation differences of the study site to document various surface conditions at submeter scales. Geometric documentation of natural surfaces for geological analyses with high resolution and precision, in the decimeter to centimeter range, became possible with advances in surveying technology (photogrammetry) and laser scanning (Buckley et al., 2008; Tarolli et al., 2009; Pfeifer et al., 2011). With close-range sensing technologies, with limits below 100 m, another order of magnitude in resolution and precision was reached (Hoffmeister et al., 2012; Milenković et al., 2015). Terrestrial laser scanning in complex and rough-surface terrain with millimeter-scale resolution was recently documented by Nield et al. (2013) and Arav et al. (2014). However, these studies considered only a few scan positions and a limited number of 3D points. Although surface conditions are quite different between oyster shells and grain or rock modeling applications (Feng and Röshoff, 2004; Heritage and Milan, 2009; Hodge et al., 2009b; Wang et al., 2013), there are still certain methodical similarities between our approach and theirs to evaluate shape, orientation, distribution, and size.

Process analysis of surface in geology using high-resolution TLS was described in a multi-scale approach by Brasington et al. (2012) on a river bed surface, which was comparable to the work of Baewert et al. (2014). The goal of Baewert’s research was to examine the influence of the number of scan positions and grid cell size on roughness calculations during postprocessing. In addition to data acquisition and surface representation from full-waveform TLS data, Di Salvo and Brutto (2014) extracted individual objects. The extraction of a number of geometric primitives, such as rock blocks, served to calculate the volumes of individual elements and to measure distance between them. In order to improve automatization of object extraction, Brodu and Lague (2012), Otepka et al. (2013), Blomley et al. (2014), and Weinmann et al. (2015) suggested deriving features for each point and classifying the point cloud with respect to those features. They described extracting optimized 3D or 2D features from the derived optimal neighborhood size to optimally describe the local structure for each 3D point. Eigen features from a covariance matrix of a point set with the sample mean are commonly used geometric features that can describe the local geometric characteristics of a point cloud and indicate whether the local geometry is linear, planar, or spherical (Lin et al., 2014). In general, estimating geometric features from point clouds has many applications in computer graphics (Merigot et al., 2011); however, this technique is a novelty in geology and, more specifically, in paleontology.

Recently, close-range TLS was applied for the first time for monitoring intertidal and subtidal oyster reefs in estuarine science (Rodriguez et al., 2014) due to their economic importance. Acoustic techniques are state of the art for these purposes (Allen et al., 2005; Grizzle et al., 2005). Similarly, coral reef surfaces were investigated by the use of remote-sensing approaches with a resolution below 1 m by Goodman et al. (2013) and Hamylton et al. (2014). All these studies tend to focus on the overall structure of the reef and do not aim to detect single shells. Consequently, millimeter resolution has not, to our knowledge, been achieved to date.

Nothegger (2011) proposed to use phase-shift scanners for the application of high-resolution and highly accurate scanning of large sites, and among them was the oyster reef. Such a laser scanner measures a cloud of 3D points. An analysis of the point cloud to improve precision (e.g., filtering) was described by Dorninger and Nothegger (2009). Therefore, we present an improved 3D data set of the oyster reef with ∼100 times more points compared with a previous campaign (Haring et al., 2009); our data set resulted in approximately one billion points, representing the entire site of ∼459 m2. This was achieved by the ability to acquire many more scanning positions with a better configuration. Image data were acquired as well. The photos were used to generate high-resolution orthophotos, which were primarily used for the interactive verification of the automatically generated results and also for the generation of textured 3D models for visualization purposes (Harzhauser et al., 2015, 2016; Djuricic et al., 2016).

Study Site

The excavated shell bed (27 m × 17 m) is protected by an indoor hall as part of the geo-edutainment park, “Fossilienwelt Weinviertel.” The hall comprises the oyster reef, surrounded by two natural walls of sand, two construction walls, a mobile bridge, a lower visitors’ platform with reflectors below, an upper platform with stairs and projectors, and reflectors on the ceiling. In August 2012, the oyster reef was declared a protected site under the nature conservation laws of Lower Austria (§12 NÖ Naturschutzgesetz) due to its exceptional scientific importance.

Control Points

The control points were distributed regularly (Fig. 2A), and they served the purpose of registering 3D laser scans, as well as being a link to an external coordinate system, as discussed in the Coordinate Systems section. Control points were placed on the pillars around the oyster reef (15 points in Fig. 2B), to provide the reference point network. In addition, 18 control points (from point number 23 to 40) were placed directly on the oyster reef from a mobile bridge spanning across the oyster reef (Figs. 2A and 2C), and seven more points were placed on the visitors’ platform around the reef (points 16, 17, 18, 19, 20, 21 and 22). Control points on the reef were distributed in a grid of 6 m × 3 m (Fig. 2A) as paper marks. All 40 points were measured by means of a total station (Leica 1200 series) from two different positions from the platform and from one position at the opposite corner; their precision is better than 2 mm.

Points 1, 2, 4, 6, and 10 (see Fig. 2A and Supplemental Fig. S11 of the reference network were permanently marked by holes embossed in the steel pillars of the construction walls (see Supplemental Table S1 [footnote 1]). The data collection was carried out in January 2014 and might be repeated in the future with new technology such as hyperspectral laser scanning (Hakala et al., 2012; Kurz et al., 2013; Puttonen et al., 2015) or imaging.

Data Acquisition and Platform

There is no possibility of positioning an instrument or a tripod on the oyster reef because of conservational concerns. To overcome this limitation, a long mobile bridge was constructed, which can be moved in a south/north direction (Figs. 3 and 4). Most data were acquired from it, but the reef borders were also visible from the visitors’ platform. To provide an unobstructed view, a long, stiff (and heavy) metal structure was built to carry the scanner (Fig. 3A) and the camera on the top (Fig. 3B), under a 3.1 m × 2.1 m tent.


We used a FARO Focus3D laser scanner that was controlled remotely from the stable visitors’ platform. This is a high-speed 3D laser scanner for detailed measurement and documentation that collects 3D points with a rate of up to 1 MHz and a range of up to 100 m with a 3.8-mm-diameter laser beam (FARO, 2016). Experimental analysis has shown that, over short distances, the single-point measurement precision of FARO TLS data is actually ∼2 mm up to 10 m range (FARO, 2016).

The surface of the exposed oyster reef was surveyed by TLS, and texture capturing was done by photographic imaging. The precision of the individual measurement is defined by the scanner specification. To increase the precision of the result, we applied point-based filtering (i.e., averaging). Hence, we are able to achieve 1 mm precision. Occlusions were minimized by the constellation of the scanner during scanning (high overlap and various incidence angles). For photogrammetric data acquisition, a uniform illumination was necessary to enable an image acquisition that was almost shadow-free.

A Canon 60D camera with a Canon EF 20 mm f/2.8 lens was used to capture images (each 5184 × 3456 pixels). On average, the ground sampling distance (GSD) was 0.6 mm, and the footprint of the image was ∼3.1 m × 2.1 m on the reef. Images were taken with 80% overlap over the longer image edge and 50% overlap over the shorter side. The camera was mounted approximately orthogonal to the oyster reef plane. Due to the low ambient light, artificial lighting (close to studio conditions) was necessary for image acquisition. Particular emphasis was taken to ensure homogeneous, diffuse lighting conditions to minimize shadows. Therefore, the structure for image acquisition was surrounded by reflectors and a tent (Fig. 3B).

Acquisition Strategy

Scanning the fossil shell bed required multiple viewpoints to ensure minimization of scan shadows and to provide homogeneous point coverage. The scanning height of 1.5 m above the reef also led to the capturing of overhangs. Thus, scans were acquired in a regular grid of 2 m × 2 m, with the distance between the laser scanning points ranging from 0.7 mm (center) to 1 mm (outer edges) per scan. In total, 83 scan positions were sufficient to collect data covering the entire reef. From all those positions, approximately one billion points were acquired at the site, corresponding to ∼150 points per square centimeter. The high number of nearly orthogonal scanning positions enabled a homogenous point density on the surface of the oyster reef. Regarding investigation of the quality of measurements obtained with laser scanners (Pfeifer et al., 2007; Kaasalainen et al., 2011; Soudarissanane et al., 2011; Milenković et al., 2015), it is recognized that not all points measured by the laser scanner can serve for modeling of the irregular surface. Hence, the most reliably scanned area of each scan is taken into consideration by restricting the distances to less than 15 m.

The scans were oriented in the acquisition coordinate system (ACS) (introduced in more detail below in the Coordinate Systems section) with the control points on the reef, the walls, and the visitors’ platform by means of the scanner software FARO Scene. In each scan, 15 or more control points (i.e., planar photogrammetric targets; Figs. 2B and 2C) were identified, covering the entire field of view. A global accuracy of better than 3 mm was achieved, but local discrepancies between overlapping scans are smaller, in the region of 1–2 mm.

Bundle Block Adjustment

Photogrammetric images were oriented in the ACS by means of a bundle block adjustment (BBA) with control points (photogrammetric targets and distinct points in the TLS data) on the reef. We used Agisoft PhotoScan (PhotoScan software, for automatic feature extraction, matching, and relative orientation of all images. Measurements of photogrammetric targets and distinct points were added manually. Subsequently, a datum transformation was applied with the control points. All observations were exported (image coordinates of feature points and control points and control-point coordinates in ACS) and used as input for a BBA in software developed in-house, thus allowing better control over camera calibration and control-point consideration.

The distinct points in the TLS data were necessary to counteract distortions in areas without targets (i.e., the lower part of the reef, which was less accessible). The distinct points were measured directly in the TLS data and were necessary for automatic detection of wrong tie-point correspondences. To avoid tensions between different coordinate systems (Luhmann et al., 2014) and realizations of coordinates, all existing stochastic quantities (accessible from the total station campaign and the registration of the TLS data) of the control points were introduced into the BBA. The camera calibration was performed by self-calibration (also in the BBA). The mean error of tie-point coordinates close to the reef center was ∼1.5 mm (in X and Y directions) and 2.4 mm (in Z direction); and close to the edge of the reef, it was 2 mm (X and Y) and 4 mm (Z).

Coordinate Systems

During its formation in Early Miocene time, the living oyster reef as well as the shell bed were horizontal. During the Middle Miocene, the entire basin fill, including the shell bed, was tilted tectonically by ∼25° to the west (Zuschin et al., 2014). We thus established the following coordinate systems to reconstruct the original situation during shell bed formation for subsequent geological and paleontological analyses: (1) acquisition coordinate system (ACS); (2) local analysis coordinate system (LACS); and (3) Universal Transverse Mercator (UTM) system.

The ACS is a local coordinate system used during the measurements with a total station and a TLS (Fig. 2A). Polar measurements (angles and distances) were used to determine the 3D positions of points on the object relative to a well-defined point where the total station was placed, which was the origin of the acquisition coordinate system (point “a” in Fig. 2A and Supplemental Fig. S1 [see footnote 1]). A second well-defined point on the visitors’ platform was selected for purposes of orientation (point “b” in Fig. 2A). The horizontal direction from the first to the second point was set to zero.

The LACS (Fig. 4) is a coordinate system defined by a Euclidean transformation from the Acquisition Coordinate System. It represents the “horizontal geological” coordinate system. The axes follow the dip and strike directions, and the third axis, Z, is orthogonal to the adjusting plane through the field. Axis X extends downward in the field, and Y is horizontal. Thus, the Y axis is basically parallel to the length axis of the hall. The origin of the local analysis coordinate system was placed behind the corner of the oyster field, which is diagonally opposite the entrance of the hall.

Transformation from the ACS to the LACS and from the ACS to the UTM

The oyster reef is situated on a plane with small deviations relative to the total extent of the reef. An orthogonal regression plane was determined from a subset of the TLS point cloud. This subset was a sub-sampled point cloud that ensured the global approach (working on the entire reef at once) was successful by avoiding numerical problems (covariance matrix of one billion points). The plane was determined by a principal component analysis (PCA) of the point cloud. From this, the rotation matrix to transform the points from the acquisition coordinate system into the local analysis coordinate system was determined. The entire area of interest was rotated by 20.26° to define a locally horizontal geological coordinate system, the LACS.

The transformation from the ACS into the LACS is described by a transformation comprising rotation, translation, and scale identical to 1.0. The local analysis coordinate system has the same scale as the range measurements of the total station (Supplemental Table S2 [see footnote 1]).

Finally, a transformation to the UTM System was also performed to ensure compatibility with future campaigns and to allow for comparison with other geological data from the region. The transformation parameters are given in Supplemental Table S2 (see footnote 1).

The registration method used control points (only the reference points on the pillars obtained from the total station measurements) to transform multiple scans into the local analysis coordinate system. These control points are transformed to UTM33N coordinates on the basis of four east, north, up (ENU) control points observed by GPS measurements outside the oyster hall (Supplemental Fig. S1 [see footnote 1]). When the parameters obtained from the spatial similarity transformation are applied, the points fit with a root mean square (rms) coordinate difference of 5 mm. Height is obtained from GPS measurements based on ellipsoidal height above WGS84 per the GRS80 ellipsoid (Hofmann-Wellenhof et al., 2011).

Data Management

The large data volume, 50.3 GB, required piece-wise processing of the data. Therefore, the data were organized into 81 rectangular tiles, with 2–17 million points per tile (110 MB to 1 GB). In total, 29 tiles are in a range from 110 MB to 500 MB; 21 tiles are in a range from 500 to 900 MB; and 32 tiles are in a range from 900 MB to 1GB. Not all of the 81 tiles are covered entirely with shells; some also contain plain sand areas. The area covered by shells and visible from the scan positions is 367 m2.

The tiles were defined with an extension of 2.1 m (east-west) by 3.1 m (north-south). An overlap of 5 cm (i.e., 10 cm in total) was chosen to avoid border effects during processing, thus resulting in 81 tiles in a grid structure, with a starting point xo = 30 m, yo = 0 m, and the grid offset 3 m in Y and 2 m in X.

In this section, we describe digital surface model and orthophoto generation, as well as manual reference data collection. We also present the method for automatically detecting shell surfaces and numbering them (see workflow diagram in Fig. 5).

Point-Cloud Manipulation

To improve the point cloud, outliers were removed, and noise was filtered (Nothegger and Dorninger, 2009). The improved point cloud was input into the meshing process, where a polyhedral mesh was created (see below). The consistency of the point cloud was aided by exploiting overlapping scans—up to five, depending upon the location within the reef, thus minimizing errors caused by autocorrelated noise in the distance sensor of the FARO scanner. These errors are not reduced by filtering a single scan only. Due to the complex surface structure, it often occurs that the emitted laser beam hits two or more distinct surfaces simultaneously. It is impossible for scanners utilizing the phase-shift measurement technique to discriminate between those signals solely on the basis of the integral signal received by the scanner (Nothegger and Dorninger, 2007). Considering the scanning position for each individual scan, locally estimated surface normals can be analyzed to eliminate those points being characterized by normals facing almost orthogonal to the direction of the scanning beam. Additionally, single outlying points were eliminated as outliers after multiple scans were merged, thus identifying apparent (not actual) surfaces that were not facing the scanner viewpoint but lay inside the beam direction. Registration of the individual scans must be sufficiently accurate to avoid additional errors. That goal was realized by using at least 15 control points per scan, with a precision of better than 2 mm (see above).

The 3D model based on the TLS data was generated using the Poisson surface reconstruction method (Kazhdan et al., 2006). The result is a triangulated surface with an approximately homogeneous edge length of 0.7 mm on average. The method also reduced the number of points to ∼40% of the original TLS points on the reef. The Poisson surface reconstruction algorithm is especially suitable for merging point clouds or meshes because of its implicit averaging, resulting in a surface that appears smooth even in the presence of noise in the individual point clouds. The reason for this is that, internally, the surface is reconstructed as a twice-differentiable scalar field, and the triangulation—the output of the algorithm—is an approximation of that smooth scalar field. Due to this averaging of up to five individual point clouds, the localized upper limit of the standard deviation between the triangulated mesh and the generated surface mesh is up to 2 mm, given that the surface is sufficiently flat. In areas with high local curvatures, the area of the laser beam footprint (3 mm) limits the achievable accuracy.

By using the nodes of the triangulated model, a homogeneous, high-resolution 3D point cloud is generated for further analysis. A high-resolution and dense 3D point cloud of a small part of the oyster reef is shown in Figure 6.

DSM Generation

Surface models may be derived in the form of a regular grid. The digital surface model (DSM) describes the top surface of the terrain (see Animation 1,2 and Fig. 7). The regular DSMs of the oyster field, rather than the unstructured point clouds, are appropriate in this stage of the work to obtain the relief of the reef surface—estimating concavity, convexity, or flatness of the locally selected neighborhood region. Detailed high-resolution (millimeter) surface estimation allows analysis of the surface model for features at the centimeter to decimeter scale (Hodge et al., 2009a; Bertin et al., 2016). However, concerning our application, work with point clouds will be used in future research.

Because the oyster reef surface includes overhangs, not all details can be maintained when converting from the point cloud to a scalar-valued function (2.5D representation) parameterized over the reference plane. In areas of overlap, the uppermost surfaces were chosen to be represented in the height model, as they were always scanned, whereas the lower surfaces were not necessarily scanned due to occlusions.

The points below overhangs were excluded by using point-cloud processing on a cell basis. In cells of 1 mm × 1 mm, the point with the maximum height value was included in the interpolation process. The interpolation method was moving least squares (Lancaster and Salkauskas, 1981), with a tilted plane used as the surface model and interpolating the eight nearest points at the grid post. Grid spacing was 1 mm. The maximum radius from the grid post to the points was 6 mm, but, with the exception of border areas, the radius was rarely reached by the eight nearest points.

Computation was performed using scientific software for orientation and processing of airborne laser scanning data (OPALS; Pfeifer et al., 2014). The number of points for the tiles to be interpolated was 14 million points on average, and this was reduced to 6 million grid posts, excluding the overlapping area.

The digital surface model and corresponding hillshade data are publicly available at PANGAEA (see footnote 1).


Texture is assigned to the DSM for more realistic visualization of the oyster reef. With the knowledge of exterior and interior orientation of the recorded images that have been computed by the bundle block adjustment, it was possible to project the RGB information onto the DSM surface, resulting in an orthophoto (OP). In total, 300 individual images with a nominal ground resolution of ∼0.6 mm per pixel were used in the OP generation process. The high-resolution OP has a pixel size of 0.5 mm. A radiometric correction (contrast and brightness adjustment) was applied on the images to minimize the remaining influences of artificial lighting, while enhancing contrast by histogram stretching. Thus, the image colors do not appear natural, but they are optimized for visual interpretation in the context of scientific paleontological analysis (Fig. 8). The 3D distance between any single DSM point and the camera projection center (depth) was in the range of 2.4–3.6 m, due to reef topography.

The orthophoto data are publicly available at PANGAEA (see footnote 1).

Reference Data

For automated analysis of the oysters in the reef, the detection of individual objects is necessary. For the purpose of evaluation of our method, the automatic detections need to be compared against corresponding reference data. Reference data for object delineation was obtained manually for 10,284 objects. The majority of objects represent Crassostrea gryphoides, which is also the largest species found on the oyster reef. It is up to 60 cm long, according to data analysis by Harzhauser et al. (2015). The brownish and massive calcitic oyster shells are represented exclusively by disarticulated single valves. Expert knowledge is required to detect and outline the individual shell correctly and to determine its parameters—species; level of overlap; shell side (left, right, unknown); orientation (convex side up/down, not known); fragmentation (complete, low, moderate, high); and tile number. Digitization of the validation data was performed manually in ArcGIS, producing vector data stored as shape files. The fastest digitization was achieved with the shaded relief of the DSM as the bottom layer, a semi-transparent orthophoto superimposed onto it, and the previously digitized outlines on the topmost layer. One row and one column of tiles were digitized manually (Fig. 9); these 13 tiles form a cross in the center of the oyster reef and allow computing gradients based only on reference data.

Automatization Method for Individual Shell Detection

A major goal of this study is developing a method for the automatized detecting and counting of individual shells. Mature Crassostrea valves have a surface area ranging between 2000 mm2 and 60,000 mm2 in this geological exposure, while juvenile or fragmented specimens can be smaller. Specimens with the convex side facing upward have a rough surface and are distributed over the entire reef. A single Crassostrea valve typically has a thickness of ∼25 mm, often elevating the top surface of the shell entirely from the reef. To detect those specimens, a measure of convexity was targeted; specifically, openness was computed for each grid point (Figs. 10 and 11).

The aim of the openness feature is to provide a measure of visibility of the neighborhood from the viewing point. It includes all surrounding points that are in line of sight with that viewing point and excludes points that are blocked by local terrain or any obstacles. It was introduced by Yokoyama et al. (2002) as a cone-fitting approach to measure the opening angle of a cone centered at a grid point and constrained by the neighboring elevations within a specified radial distance (Fig. 11). It is basically a convexity/concavity measure calculated from local terrain profiles in principal compass directions (Mandlburger et al., 2009). In most cases, that measure works well in distinguishing convex surfaces from the rest of the terrain. There was a further investigation to find which size of neighborhood is optimal for analysis of the variation in the rougher surfaces. Thus, the feature provides an angle of the fitted cone, which is derived for each point, with a preset kernel (k) radius (e.g., neighborhood size is: 2k + 1; k = 1 mm, 3 mm, 5 mm, . . . , 27 mm). Openness is expressed in two models—positive for concave features and negative for convex features.

The following quantities are introduced to compute openness (see also Fig. 11): A (Xa, Ya, Za)—reef view point; B (Xb, Yb, Zb)—point restricting the view in the neighborhood of A due to topography; C (Xc, Yc, Zc)—any point in the neighborhood of A; D—horizontal distance between A and B; φ—azimuth angle (0°, 45°, 90°, . . . , 315°); r = 2k + 1—kernel radius (i.e., neighborhood size) in pixels; θ—elevation angle between reef view point (A) and each grid point (B) located along four sections (axis parallel and diagonal; see Fig. 11) with specified azimuth angle φ and kernel radius k; α—zenith angle at a reef DSM grid point calculated along one of eight azimuths φ and specified kernel radius k; —mean of angles along r; n—number of section points dependent on k; point B is the point that has the minimum α for each section.
The applied workflow is based on a raster of negative openness, i.e., a raster f(x,y). The threshold t is set to zero, and values are converted into binary values (blue and white in Fig. 12D), where value 1 is assigned to foreground and 0 refers to background, i.e., g(x,y).

Single pixels and small areas are excluded for further analysis based on a minimum area threshold (Fig. 12E). Mathematical morphology is applied to obtain smoother outlines and to fill in small gaps (few pixels). A majority filter replaces the neighboring pixel values, depending on the size of the kernel (Fig. 12F); for example, m = 1 (3 × 3 mm), m = 2 (5 × 5 mm), up to m = 6 (13 × 13 mm). Further processing includes connected-components analysis on the basis of the 8-neighborhood array. Pixels in the neighborhood are connected in regions. Each component obtains a unique ID. Figure 10A shows automatically detected convex surfaces overlapping the mesh model of the oyster reef. In the first step, the area of connected components is incompletely filled due to the diversity of concavity and convexity. To achieve higher completeness detection, a flood-fill operation is performed (Fig. 12G). This operation changes connected background pixels (0s) to foreground pixels (1s) and stops when it reaches object boundaries. Correspondingly, the method is applied to the whole data set to assess the applicability of the method to larger areas, and the results are analyzed.

Figure 13A shows the processed results of the tile adjacent to the north of the center tile (see Fig. 9) as an example of the applied method: true positive (TP) pixels are colored in blue; false negative (FN) pixels are colored in dark yellow; false positive (FP) are presented with red color; and true negative (TN) features are in gray. Results are evaluated against the manual reference data (Fig. 13B), where the results are assessed for more than 10,000 specimens by comparison with the manually determined shell polygons within the central area of the reef (Fig. 9). The surface matching validation method uses completeness (reddish background in Supplemental Table S33), correctness (bluish background in Supplemental Table S3 [see footnote 3]), and quality assessment (greenish background in Supplemental Table S3 [see footnote 3]) as proposed by Heipke et al. (1997); see also Supplemental File. For the given scenario, the detection accuracy, based on validated surface matching and visual inspection of obtained results, reveals the feasibility of the proposed methodology.

The result shown in Figure 13 is an optimized result based on test runs using different kernel sizes: k = 1 (3 × 3 mm), k = 3 (7 × 7 mm), k = 5 (11 × 11 mm), and up to k = 27 (55 × 55 mm). The test runs showed that completeness and quality decreased from k = 1 to k = 3 and start increasing from k = 3 up to k = 25. Correctness is maximal at k = 3, but quality is minimal. The quality results show independence from majority filters, indicated by only minor changes (Fig. 14 and Supplemental Table S3 [see footnote 3]). Therefore, parameters k = 25 and m = 3 (majority filter size 3) are chosen as the best parameters for the automatic surface detection of shells. In addition to the surface matching, the total number of detected objects was examined.

Thus, the number of detected objects was validated by means of the following parameters: (1) the best parameters estimated for surface matching, k = 25 and m = 3, and another pair of parameters were determined to be optimum for the number of detected objects; (2) k = 7 and m = 6. Surfaces of the shells were detected effectively for 69% of the data based on surface matching using k = 25 and m = 3 (Fig. 15, boldface in Table 1 and Supplemental Table S3 [see footnote 3]) and for 50.25% of the data based on the detected number of objects. Overall, an accuracy of over 98% was achieved in detecting the number of objects by the use of another pair of parameters, k = 7 and m = 6 (Supplemental Table S43); but then surface matching completeness decreased slightly to 67%. According to the reference data (13 tiles), the number of shells (excluding highly fragmented shells smaller than 17 cm2) ranges from 322 to 466 shells per tile with a mean of 387. Within those tiles, the number of automatically detected shells (based on k = 25 and m = 3) ranges from 170 to 210 shells per tile with a mean of 195. On the basis of k = 7 and m = 6, shells range from 310 to 410 with a mean of 381.

The automatically detected shells and the reference data were also used to estimate the total number of shells on the reef. A total of 10,284 shells served as reference data (corresponds to 13 tiles or 78 m2). In the reference data, the number of individual shells larger than 17 cm2 was 6148. If overlapping shells are counted as one shell (as they typically appear in the automatic procedure), the number was reduced to 5035, the factor 1/1.23. The number of automatically detected objects over the region of the reference data was 4958 when the optimal numbers for object detection (k = 7, m = 6; boldface result in Table 1) were used. The inner reef area is completely covered with shells and comprises 49 tiles, corresponding to 294 m2. The outer tiles are not entirely covered by shells, since shells are only found at an area below 6 m2, i.e., the area of tile. The number of automatically detected shells larger than 17 cm2 in the inner reef area is 18,264. Multiplied by a factor of 1.23 to account for overlapping shells, this leads to a number of 22,432 shells. In comparison, multiplying the number of manually delineated shells in the 13 reference tiles (6148) by 49/13 would correspond to 23,173 shells. Thus, the automatically detected number of shells is 97% of what is computed from the reference data and the area ratio, which compares well with the 98% (see above; Table 1 boldface result). Extending this estimation from the 49 tiles (294 m2) to the total area of the reef that is covered with shells, 367 m2, the estimated number of shells, based on the automatic detection and the factor 1.23 to account for overlaps, is 28,002—that is, the estimated number of shells larger than 17 cm2 for the entire reef area. When the ratio of all shells to large shells in the reference area, 10,284/6148, is applied, the estimated total number of shells is 46,840. This number includes complete but also fragmented specimens.

Principally, openness was chosen as the most reliable geometrical feature for the analysis. On the basis of our observations, the feature openness describes very well the oyster reef surface roughness and supports visualization of various similar reliefs, as described in Yokoyama et al. (2002), Zakšek et al. (2011), and Doneus (2013). In general, surface roughness plays a key role in determining characteristics of the oyster reef terrain, where shells of different sizes are widely distributed. Convex-up surfaces provide information about the respective exterior parts of the shells, which are more dominant on the reef than interior parts (Harzhauser et al., 2015). The term “convex-up,” indicating the upward position of the convex part of a curved bivalve shell, is standard in analyses of shell beds (Martin, 1999; Harzhauser et al., 2015, 2016). Consequently, convex-up surfaces appear in processed data as negative openness (Fig. 11), which was used in our method as a main feature to differentiate between objects.

After applying the negative openness feature, the method worked insufficiently for convex-down oriented shells, because the interior parts of those shells may have concave geometries. Convex-down shells were partly detected via a border band, while information from the interior part was missing. Therefore, mathematical morphological operations such as flood-fill and majority filter improved the completeness of the results and helped to better structure connected components. Thereby, the outlines of the oyster shells were smoothed, noise was minimized, and small remaining gaps were filled exactly inside the concave shells and not in the other concave parts of the surrounding terrain (Figs. 12F and 12G). That procedure was effective, since concave (i.e., convex-down) shells have noticeably narrow convex borders that form mostly closed components (approximately elliptical rings). Those narrow borders were also possible due to the change of curvature, which was still clearly noticeable on the margins for unburied and flat shells.

Kernel-size optimization has a large influence on the results. Loss of information might occur if the kernel size is enlarging (Fig. 15) because increasing the neighborhood size corresponds to the scale of surface roughness. The value of the kernel size should not be too large (e.g., not larger than the width of the shell), if the aim is to capture details on the reef, such as low or moderate fragmentation of specimens. Herein, we developed a computational tool for selecting the kernel size to include enough information to correspond as closely as possible to the reference. If the kernel size ranges between k = 15 and k = 21, noise and numbers of small polygons decrease. In this analysis, the focus was on complete oyster shells and slightly and moderately fragmented specimens. Therefore, an optimal kernel had to be selected to minimize noise in the surroundings. Moreover, polygons smaller than 17 cm2 (or highly fragmented specimens) were eliminated from automatic detection. The same criterion was applied to the reference data (Fig. 9), where minimal polygons are excluded, because complete Crassostrea valves have a minimum surface area of 17.98 cm2.

With regard to validation, results from automatic object detection were compared with ground truth data under the same conditions (i.e., two overlapping shells in the manual data set are presented as one object due to rasterization and likewise in automatic detection). High-resolution 3D point-cloud data show high potential for surface categorization and modeling. Our final results document the success in both surface matching and object detection. Thus, irregular shapes may be reliably recognized by means of convex surface information, although the object shapes are challenging for automation. Still, separation of overlapping parts is a work in progress due to the virtual interruption of the borders of the underlying shells and also due to many irregular protruding parts. That difficulty may be handled in future work by investigating and including more features, such as radiometrically corrected intensity data.

Our study aimed at high-resolution digital documentation of the world’s largest known fossil oyster reef, comprising thousands of shells of Crassostrea gryphoides. The entire site has been surveyed by terrestrial laser scanning and digital imagery at resolutions of 1 mm and 0.5 mm, respectively. This paper describes data acquisition and processing to maximize scanning efficiency of the very complex surface. The unfavorable data acquisition conditions made the task even more challenging in terms of low visibility, tilted geometry of the site, and restricted access to the surface. Nevertheless, our study has shown that, with a sophisticated configuration of TLS, image acquisition, and processing, we were able to deal with the given conditions and to significantly improve the current scientific understanding of the oyster reef. In total, 83 scan positions, 300 images, and 40 ground control points (GCPs) were acquired and measured and were transformed from the acquisition coordinate system into a local analysis coordinate system, as well as to UTM. The resulting filtered point cloud consists of one billion points with an average resolution of 1 mm. It represents a reliable and objective documentation of the entire excavation site, providing data for current and future paleontological and geological analyses. To make the captured surface accessible to the scientific community, various visualization techniques have been applied including 3D triangulated surfaces, shaded digital surface models, and color-coded relief models. Furthermore, an experiment on automatic analysis of high-resolution laser scanning data was performed, in which geomorphometric derivatives were calculated, such as openness (convexity/concavity). We developed a method to identify and enumerate convex parts of the shells automatically. The method is based on geometric feature extraction, morphological operations, and connected component analysis. This method could also be applied to studies related to geological reefs, volcanoes, or rocks where high-resolution surface data are required. The performed quality assessment resulted in 69% successful detections based on shell surface matching (k = 25 and m = 3) and over 98% of detected objects with the applied neighborhood size 15 × 15 mm (k = 7) and majority filter m = 6. The manually derived reference data were used to cross-check and validate these values.

Our data set is the largest of its kind in geology and allows for a broad range of scientific analyses, including taphonomy, structural geology, tectonics, and paleoecology (e.g., Harzhauser et al., 2015). Aside from these purely scientific aspects, the fossil oyster reef is a protected natural heritage site. It is the main attraction of a geopark, and accidental damage or vandalism by visitors, as well as destruction by natural hazards, are constant threats to this unique site. The digital models are a preeminent and unprecedented documentation of the shell bed and may serve as a basis for restoration and maintenance. Accordingly, this study represents a best-practice example of documentation in earth sciences.

The Smart-Geology project is financed by the Austrian Science Fund (FWF) project 25883-N29. The authors are grateful to the management and employees of the “Fossilienwelt Weinviertel” (Stetten, Austria) for providing access and organizational and practical help during the data acquisition. We would like to thank Andreas Roncat for his kind support during data acquisition. We want to thank Sharif Hasan for participation in reference database development and to thank Francisco Arroyo Pinilla, whose activities have been carried out within the European Union Program Leonardo da Vinci. We also want to thank native English speakers Michael Hornáček and Glen Foss for valuable comments and proofreading of the manuscript. Balázs Székely contributed partly as an Alexander von Humboldt Research Fellow. We thank the Associate Editor Richard Jones and two reviewers, Simon Buckley and anonymous, for their excellent efforts to improve this paper.

1Supplemental File: Figure S1, Tables S1 and S2, digital surface model and corresponding hillshade data, and orthophoto data can be accessed via PANGAEA at
2Animation 1. Visualization of entire oyster reef using shaded digital surface model with moving virtual sun every 27 degrees (azimuth) and with fixed zenith angle of 45°. The animation appears in the PDF and is also available at, or go to the full-text article on to view Animation 1.
3Supplemental File (including Supplemental Tables S3 and S4). The file is available at, or the full-text article on to view the Supplemental File.
Gold Open Access: This paper is published under the terms of the CC-BY license.