3 D CENTRAL LINE EXTRACTION OF FOSSIL OYSTER SHELLS

Photogrammetry provides a powerful tool to digitally document protected, inaccessible, and rare fossils. This saves manpower in relation to current documentation practice and makes the fragile specimens more available for paleontological analysis and public education. In this study, high resolution orthophoto (0.5 mm) and digital surface models (1 mm) are used to define fossil boundaries that are then used as an input to automatically extract fossil length information via central lines. In general, central lines are widely used in geosciences as they ease observation, monitoring and evaluation of object dimensions. Here, the 3D central lines are used in a novel paleontological context to study fossilized oyster shells with photogrammetric and LiDAR-obtained 3D point cloud data. 3D central lines of 1121 Crassostrea gryphoides oysters of various shapes and sizes were computed in the study. Central line calculation included: i) Delaunay triangulation between the fossil shell boundary points and formation of the Voronoi diagram; ii) extraction of Voronoi vertices and construction of a connected graph tree from them; iii) reduction of the graph to the longest possible central line via Dijkstra’s algorithm; iv) extension of longest central line to the shell boundary and smoothing by an adjustment of cubic spline curve; and v) integration of the central line into the corresponding 3D point cloud. The resulting longest path estimate for the 3D central line is a size parameter that can be applied in oyster shell age determination both in paleontological and biological applications. Our investigation evaluates ability and performance of the central line method to measure shell sizes accurately by comparing automatically extracted central lines with manually collected reference data used in paleontological analysis. Our results show that the automatically obtained central line length overestimated the manually collected reference by 1.5% in the test set, which is deemed sufficient for the selected paleontological application, namely shell age determination.


INTRODUCTION
In the past decade significant progress has been made in 3D digitizing sensor technology.To cope with the generally high amount of data, processing methods have been automated to a high degree.Therefore, 3D digitizing of shape (geometry) and appearance (images) of a huge paleontological excavation site was employed to derive digital information in form of a point cloud and orthophotos of the world's largest fossil oyster reef (Djuricic et al., 2016).The knowledge about the fossil oyster shells and digital documentation are of geological and paleontological importance and gained scientific interest because they support the evaluation of hypothesis concerning the generation of the oyster reef.The distribution and pattern of orientation and length of the individual shells can help in identifying possible causes for the current appearance of the oyster reef, e.g. if it was caused by a storm or a tsunami event (Harzhauser et al., 2015).
Often, information on, e.g., orientation and length of specimen is only available from traditional field measurements, which are not only subjective and time-consuming, but also limited in the number of measurements taken, which leads to limited spatial sampling.Furthermore, in protected sites access limitations apply, thus rendering tactile measurements impossible.Finally, scientific communication in a broader community requires original data to be available for all participants, leading to more transparent results.
The digital documentation of the oyster reef brings attention in palaeontology and sedimentology; through the restoration of the palaeosurface, age structure and size class estimations, and quantification of carbonate production per shell (Harzhauser et al., 2016).
The automatic extraction of shells from photogrammetric data, e.g., digital surface models, is studied in (Djuricic et al., 2016).There, the visualization and detection of convex shell surfaces within the complex surrounding on the oyster reef is studied, as well as the automatic quantification of the shell number.
Open questions in palaeontology concern the carbonate production of the oyster reef.In Harzhauser (2016) a relation between shell length and volume was generated.For that study the length was determined by centre line extraction.As the entire oyster reef is documented in the form of a digital surface model, the determination of the length of each specimen, and therefore the estimation of the carbon production for the entire reef becomes feasible.What is missing is a method for the reliable and accurate computation of the centre line of each oyster of the reef.
The Crassostrea oysters often provide habitat for other species, which can support an estimation of the geological age of the oyster reef.For example, the species Ostrea lives often on the exterior part of a Crassostrea oyster.Automatic derivation of encrustation of Crassostrea by Ostrea would thus further support paleontological analysis providing a methodology for documenting paleontological specimens without going to the locality.Consequently, we formulate the following aims of the research: • to derive three-dimensional (3D) shell central lines automatically by adoption existing and by developing new methods using available digital information; • to estimate shell 3D length; • to evaluate the accuracy of two-dimensional (2D) and three-dimensional (3D) shell length against manually collected reference; • to apply ratio between shell area and length with the aim of identifying encrustation or strongly curved shells.

Related Work
Generalization of object central line is studied using images, mesh surfaces, or point clouds.The state of the art presents works in a wide variety of applications of 2D and 3D approaches in computer graphics, mobile mapping and GIS, medical image analysis, etc.However, 3D central lines on surfaces in paleontological context, especially for fossil oysters using photogrammetric technology, are not yet reported in the literature.Our algorithm is similar to the algorithms that use Delaunay (1934) triangulation andVoronoi (1908) diagrams as intermediate steps to generate the vertices of the polygonal representation.
The terms central line, skeleton and medial axis are used often in literature to describe related things (Greenspan et al., 2001;Palagyi et al., 2001).The terms central line and medial axis appear most often concerning thinning of the object on the surface such as road (Doucette et al., 2004).The term skeleton is more associated to volumetric data such as 3D trees, sculptures, body (Tagliasacchi et al., 2009;Livny et al., 2010;Bucksch and Fleck, 2011;Kasap and Magnenat-Thalmann, 2011).The review paper from Cornea et al. (2005a) categorizes many of the existing skeleton algorithms for 3D shapes into classes based upon implementation and discusses how these classes achieve the various properties.The popular approaches use: i) distance field (Gagvani and Silver, 2001), where they derive three-dimensional sampled volumes and the distance field is the minimum distance from any boundary voxel; ii) thinning method as a linear process also in the number of object voxels (Palágyi and Kuba, 1998); iii) geometric method based on medial axis transform (Amenta et al., 2001) using a finite union of balls and their relationship with surfaces where the set of centrals of the balls form the medial axis of observed surface; iv) potential field such as vector field (Cornea, 2005b), i.e. a 3D array where each voxel contains a vector value (magnitude and direction) and its topological characteristics are used such as critical points and critical curves, to extract the curve-skeleton.
Moreover, pixel based methods were studied by Lam et al. (1992) when thinning methodologies based on pixel deletion or erosion criteria were needed to preserve the linear connectivity of patterns.Niblack et al. (1992) described also an algorithm for generating connected skeletons of objects in a binary image.The most dominant applications nowadays for centre line extraction based on images are road mapping and medical purposes.Reconstruction of road central lines from mobile mapping image sequences which can be used to update a Geographic Information System (GIS) database is introduced in (Mayer et al., 1998;Tao et al., 1998;Doucette et al., 2004;Miao et al., 2013).Approach for vascular medical imaging (Krissian et al., 2000) uses gradient information at a given distance of the vessel central lines tested on real X-ray images of brain vessels.An overview on medical centreline extraction is given in (Schaap et al., 2009).Furthermore, a geometric-based framework for extracting curve-skeletons is applied directly on the mesh domain and does not require voxelization (Au et al., 2008).The method extracts a 1D skeletal shape by performing geometry contraction using constrained Laplacian smoothing and it is demonstrated through various bone-skeletons examples.Raumonen et al. (2013) analyzed laser scanning data with standard fitting procedures to obtain the orientation and size of tree structures.The method presents the reconstruction of a trunk-branch skeleton based on 3D point cloud information.Dey and Zhao (2004) show a scale and density independent algorithm that approximates the medial axis from the Voronoi diagram of a set of sample points.Gorte and Pfeifer (2004) addressed in their work 3D modelling and reconstruction of a tree in terms of stem and branches, where an algorithm has been designed in 3D voxel space.The CAMPINO method is introduced by Bucksch and Lindenbergh (2008), see also (Bucksch et al., 2009), where their skeleton is represented as a graph, which can be embedded into the point clouds based on cycle elimination in a graph as derived from an octree based space division procedure.As application, Cai and Rasdorf (2008) predicted 3D lengths using LIDAR point cloud and planimetric road central line data.Tagliasacchi et al. (2009) based their approach on the notion of a generalized rotational symmetry axis for curve skeleton extraction.With proper joint handling, this leads to complete, characteristic curve skeletons even with significant missing data.Cao et al. (2010) present an algorithm based on Laplacian contraction.The method is robust to noise and sample distribution, and can handle a moderate amount of missing data.Also Huang et al. (2013) demonstrated L1-medial skeletons extraction from unorganized, unoriented, and incomplete 3D raw point clouds.

Sensors
Data was captured with the FARO Focus 3D laser scanner mounted on the mobile bridge.It collects at 1MHz 3D points in the range of 0.6 m to 100 m with a 3 mm diameter laser beam and a ranging error std.dev. is 2 mm (FARO, 2015).A Canon 60D camera with a Canon EF 20 mm f2.8 lens was used to capture images (each 5184 x 3456 pixels).In average, the ground sample distance was 0.6 mm, and the footprint of the image was about 3.1 m x 2.1 m on the reef side.Images were taken with 80% overlap over the longer image edge (comparable with the overlap in "flying direction") 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 was necessary for image acquisition (Fig. 1).

High-Resolution Point Cloud Dataset
From multiple scanning stand points around 1 billion points were acquired at the site, corresponding to about 150 points per square centimetre (Fig. 2).To improve the point cloud, outliers were removed, the amount of data was reduced and noise was filtered using the Poisson surface reconstruction method (Kazhdan et al., 2006).The Poisson surface reconstruction algorithm is especially suitable for merging point clouds (up to 5 in our case) because of its implicit averaging.

Method for Oyster Shell Central Line Extraction
Linear approximation of real objects can be simplified representation of complex object forms such as a large fossil oyster shells.Therefore, this paper introduces automatic central line determination for 1121 complete (not fragmented) Crassostrea shells.We used complete and closed boundaries as input for our method (Harzhauser et al., 2015).Oyster shell boundaries in 2D were extracted manually from digital surface models and from orthophotos.(Blum, 1967).The MAT is also a connection map between the Voronoi vertices.To determine the longest path in the connected tree (Fig. 3), the edge points (red) in it were selected, i.e. points with only a single neighbour.The longest path and its length were determined by applying Dijkstra's algorithm (Dijkstra, 1959) between all edge points in the connected tree.Furthermore, we looked for the potential intersections between the longest path end points and the border outlines by checking conditions: if beginning/end ratio threshold for the central line (of total length) and opening angle threshold were satisfied, the line is extended to boundary intersection.The junction points need to be far away from the oyster beginning/end point less than 30% (or more than 70%) of the total length of the central line; the opening angle of the junction point and two neighbourhoods surrounding points (previous and further) has to be less than 155 degrees.If both conditions are satisfied the path is extended to the outline by intersection.Otherwise, the end point was replaced with a junction point on the smallest opening angle and the line was then extended to boundary.Moreover, cubic Bézier curve is fitted to the points of the longest path.The curve is used to smooth the path model and to remove wiggling effect (Fig. 3).In our implementation it is assumed that the central line represents the largest dimension between the oyster beginning (part with hinge) and end (part with muscle), and that it is always in the middle between the shell borders independently of convex up or convex down shell positions (Fig. 4, Fig. 9).
The next operation includes embedding of the 2D central line into the corresponding 3D point cloud.That requires a relevant area of the original 3D point cloud for 3D central line extraction.Therefore, clipping of original point cloud is performed by excluding points that are outside of the shell boundary.Subsequently, assigning the third component Z to the 2D central line points is achieved by running an N-D nearest point search for every point using Quickhull algorithm (Barber et al., 1996).Finally, the closest 3D point is adopted as the most appropriate for 3D Euclidian distance calculation.The 3D central line implies tilt of the shells and curves along the top shell surface.The curves are caused by local slope change on the shell surface or by change of boundary shape which the central line follows.

RESULTS
Three-dimensional shell central lines are derived in objective and automated manner using 3D point cloud data.The length and orientation estimation of particular oyster shells are determined by using corresponding 3D central line.An example of 2D approximation and 3D central line results are shown in Fig. 4.
The described method to approximate a 2D central line produces a graph with many path options within a shell due to diverse shell shapes.The longest path information is stored in a database for 10.284 specimens, because this information was used to support length estimation of Crassostrea and other species (Harzhauser et al., 2015).However, in this paper we will limit the evaluation to the length of the 1121 complete Crassostrea shells.Fragments are not considered herein, as they are not informative for length/age calculations.In 3D, the situation becomes more difficult because of the wavy, rough and complex shell surface.It is common that Crassostrea provides habitat for other species, hence it happens that for instance specimens of Ostrea were attached to Crassostrea shells (when alive), or Venerupis shells lie inside the interior part of Crassostrea (after death) (Fig. 6).Therefore, 3D information is added to 2D line and 3D length is calculated (Fig. 5).
Two-dimensional and three-dimensional lengths of 1121 shells are compared in order to visualize possible differences between 2D and 3D lengths.Differences (residuals) support automatic discovering of strongly curved shells and potential encrustation of Crassostrea by Ostrea, Venerupis or any other type of shells or even fragments.The length differences (Fig. 5, Table 1) can be interpreted as consequence of many height surface variations.It appears that the differences are mainly under 10-20% of total shell length which most likely represents a group of oysters with no overlap (i.e.encrustation) or not strongly curved shells.Those shells are lying on the uppermost surface of the reef bed or are surrounded with sand.Interpretation of residuals over 20% is interesting as these shells may potentially overlap with another shell, barnacles, and fragments such as shown in 3D visualized cases in Fig. 6, or the strongly curved case in Fig. 8.In order to be on the safe side, a constraint is applied: if a residual is significantly large and there is no overlap and no strong curve along the surface, the distance from points of 3D central line to the shell plane is checked.If there is an unexpected or sudden jump in point height, then those points are removed, and the 3D length is recalculated.Such cases may appear due to many natural overhangs because of the shell local positions; or for the reason that irregularly shaped outlines were nearly cutting shell surrounding beside the shell edge sharply.
In order to explore multiple overlaps with more detail, we visualized the relationship between absolute residuals (3D-2D length) and existing number of overlap levels (from 0 to 3) for each shell (Fig. 6, Fig. 7, Table 1).
Among the large number of diverse cases, there are some particular shell characteristics such as by being long and thin, or short and wide, or long and wide, etc.Consequently, the ratio between area of the projected shell surface and the 3D-2D residuals was analyzed.The aim was to find a correlation between large residuals (> 0.05 m) and small areas as a confirmation of height surface variations.We focused mainly on large residuals, which allowed detecting the presence of encrustation and cases of strongly curved shells (Fig. 7).One of the most noticeable distinct differences between 2D and 3D length is the large residual of about 16 cm (Fig. 7) where the level of overlap is equal 0. Since there is no encrustation to cause surface variation, such case is interpreted as special type of shell with strong curved interior part of the surface.Projected 2D border line of this case and central line are much smaller than 3D approximation.Example of automatically discovered curved shell is visualized in Fig. 8. Based on the presented method and results, two hypotheses are validated for evaluation of expectation of potentially strongly curved or (multiple) overlapped shells in dataset: 1) Curved specimens are detectable if their residual is larger than 20% of the length and level of overlap is 0.
2) The ratio between the level number of specimen and residuals has an impact on encrustation identification.

Evaluation of automatic central line extraction
A method for evaluation of an estimated central line is implemented based on defined criteria: i) The length differences; ii) Centering of the line -exactly or approximately to within a specified error; iii) Visual appearance.
Oyster shell central line extraction may be used for various biological or paleontological applications and interpretations, and thus different evaluation measures may apply.We describe two measures.With these measures we differentiate between extraction capability and extraction accuracy.The first level of evaluation has to satisfy criteria of extraction capability and then extraction accuracy may be relevant for final evaluation.Therefore, we give priority to the length and compare the automatically extracted line with corresponding reference length.Based on comparison between automatically calculated length and 500 references (Table 2), an overestimated result of 101.5% is obtained.Accordingly, the total sum of extracted lengths is equivalent; they are well matching with each other.

Statistics
Further, in order to check the accuracy of centering, an accumulative, buffer wise approach is applied.We took 6 buffers of 0.5 cm around reference central line and checked their intersection with the extracted line.This approach enables an accurate comparison by taking a total sum of the closest distances from extracted lines to the references (blue, Fig. 9).gives information that all segments are within 0.03 m.Proportion of the length inside the closest buffers to reference line is high, about 58% first and 28% second buffer.In addition, central line extraction of 252 convex down and 248 convex up shells was separately evaluated.Obtained results show that accuracy assessment decreases in case of convex down extractions due to complexity of interior parts of the shell structure.

DISCUSSION
Delaunay triangulation was calculated between the filtered outline points to limit the cost of DT processing.The number depends on the method used for manual definition of boundaries; for instance, if picking is done manually with pointby-point selection; oysters have 25-75 points each.With free hand drawing (a tool in ArcGIS for detailed and dense outlines), an oyster outline can have from 500 to over 2000 points in it depending on its perimeter (about 1000 points on average).In general, oyster outlines are for most part smooth and welldefined.Thus, having the hundreds of points on the outline makes processing heavy without improving the central line search results.Therefore, outline point numbers were limited up to 100 and they were filtered by taking points with close to even spacing.Extracted 2D approximation was checked on the end points for two conditions: the junction distance form end point and the angle.Then, the extension was followed by a Bézier fit in order to improve visualization and minimize wiggle effect.
The presented approach is fully mathematical, dependent on the 2D geometry of shell outlines, and it is not affected from surface topology, i.e. pronounced shell concavity or convexity.
Since the oyster shell is composed of two valves, the convex up and the convex down valve can be recognized in an example of dataset (Fig. 9a).The valves are equilateral, joined proximally by a ligament and distally by the adductor muscle.After the bivalve death, the soft parts disintegrate and two valves disarticulate.In the present site no articulated valves were found (Harzhauser et al., 2015).The right and the left valve can be easily distinguished by the position of the adductor scar, always located near the posterior shell margin, at the third quarter of the former animal soft body.As a consequence of different mineralogy and diagenetic alternation the adductor scar is usually colored dark brown, deep purple or red-brown in the present specimens.The shell central line is, in its biological definition, an imaginary curve of the commissure plane, connecting oyster hinge with the distal shell margin, striking parallel to the maximum concavity zone, through the soft parts of the bivalve.It is roughly equidistant to the lateral shell margin, which may be distorted by the laterally projecting growth lamellae developing occasionally on the exterior shell surface (Fig. 9a).The shell curvature is a random phenomenon and a pure consequence of adaptation to surface conditions where shell is growing.In fossil shells, the central line may be confidently defined by checking the interior valve sides.However, in external view (Fig. 9a, b), the central line might be only approximated by averaging a width between the shell margins, being therefore not completely congruent with the biological central line.Post-mortem mechanical damage of the shell does not change the central line itself, but may complicate its reconstruction if the shell is fragmented or strongly abraded.The reference shell central line, as applied in the present study is depicted blue in Fig. 9a, b.Nevertheless, the algorithm is robust enough as it overestimates the size of both oyster shells (convex up and convex down) by 1.5% compared to 500 references.More than 86% of automatically extracted central lines are displaced less than 1 cm from reference lines.Moreover, the algorithm result may point out multiple overlaps and strongly curved shells in order to identify them automatically.Thus, this information is sufficient to support paleontological experts during interpretation, age structure calculation and to document the presence of encrustation on the oyster shells.

CONCLUSION
The relationship between residuals of 2D and 3D central line lengths and number of shell overlap levels were examined.Significant deviations of 2D/3D ratio work as a good quality estimator.Some key findings are emphasized herein: i) high resolution 3D point cloud data provides confident basis for 3D central line shell derivation, ii) the proposed 3D approach is efficient in estimating 3D length and iii) certain central line properties have direct relationship with the encrustation estimates.Findings are proven by presented analysis.Normalized residuals under 20% of normalized lengths and over 20% differences indicate possibility of multiple overlaps or strong curvatures.
We demonstrated that 3D central lines can be extracted automatically from various shapes of fossil oyster shells.Beside automatically determination of shell size, the purpose of central line extraction should support matching of 3D objects, by finding similar or identical object in a database.An extension to our work will be to obtain a set of points from a point cloud along the central line and describe relations between them in order to derive histogram description of data.Observing neighbourhood of each point will enable extraction of more information from surrounding points by combination of several geometrical features (planarity, sphericity, etc.).This will better explain how increasing neighbourhood size corresponds to amount of shell surface roughness, which is significant for analysis of bioerosion and abrasion on the world's largest fossil oyster reef.

Figure 1 .
Figure 1.The data collection setups; (a) TLS scanning from mobile bridge; (b) Camera under a special lightning tent.

Figure 2 .
Figure 2. 3D Point cloud of a few shells (convex down) originally on the oyster reef surface; 70 cm by 50 cm.The high number of almost orthogonal scanning positions enabled a homogenous point density on the surface of the oyster reef.Overhangs are captured due to the dense grid of scan positions.The scans were oriented in the Acquisition Coordinate System (ACS;Djuricic et al., 2016).A global accuracy of better than 3 mm was achieved, but local discrepancies between overlapping scans are smaller, in the domain of 1 mm to 2 mm.
The central line extraction algorithm includes following workflow of seven steps: 1: Constrained Delaunay triangulation between the fossil shell boundary points and formation of Voronoi diagram; 2: Extraction of Voronoi vertices within the oyster shell boundary and construction of a connected graph from them; 3: Reduction of the graph to the longest path via Dijkstra's algorithm; 4: Central line extension to the shell boundary.End nodes are inspected if this fulfils continuity conditions (see below: as a condition for accepting an elongation or not is given further below); 4a: If a node fills the conditions, then the central line is extended from the node to the oyster shell boundary by fitting a line.4b: If a node does not fulfil the conditions, then travel the central line to the nearest junction point and check for continuity.After a junction point is found, extend the central line from that node to the oyster shell boundary with line fitting.5: The extended longest path is smoothed by an adjusting cubic spline curve; 6: 3D clipping of the point cloud using oyster boundaries and assigning corresponding attributes; 7: Forming of 3D central line by imputing the height value of the nearest clipped point to every central line node; The mathematical definition of the 2D shell central line (steps 1-3) is the longest line which connects neighbouring Voronoi vertices inside the oyster shell boundary.The vertices were obtained from the constrained Delaunay triangulation (DT) of the boundary points (outline segments as constraining edges).To find the medial axis transform for the oyster outline, the Voronoi diagram (VD) was formed.Here, the Voronoi vertices within the oyster boundary are shown with cyan colour.By taking the edges between the Voronoi vertices within the boundary, we get a medial axis transform (MAT) for the oyster

Figure 4 .
Figure 4. Top row: Shaded DSM (2 m by 3 m) overlapped with extracted 2D lines (red) of complete Crassostrea gryphoides oysters (blue); Bottom row: visualization of fitted 3D central line on shell pair mesh model (convex down left shell and convex up right shell).

Figure 7 :
Figure 7: Left: histogram of absolute residuals between 2D and 3D lengths; Right: number of overlaps for each shell among 1121 Crassostrea oysters.

Figure 9 :
Figure 9: (Top two rows) Buffer zones growing around reference central line (blue) and visualized extracted line (red); (Bottom two rows) Quality assessment -overall, convex down, convex up accuracy of accumulative and buffer wise extracted line distribution.Based on quality assessment, 86% of extracted line lengths are within first two buffers.Beside extracted segments under 0.01 m distance from reference central line, the overall accuracy

Table 2 :
Length differences of 500 shells: reference vs. automatic extraction