Geomorphometric analysis of cave ceiling channels mapped with 3 D terrestrial laser scanning

The change of hydrological conditions during the evolution of caves in carbonate rocks often results in a complex subterranean geomorphology which comprises specific landforms such as ceiling channels, anastomosing half tubes, or speleothems organised vertically in different levels. Studying such complex environments traditionally requires tedious mapping, however, this is being replaced with terrestrial laser scanning technology. Laser scanning overcomes the problem of reaching high ceilings providing new options to map underground landscapes with unprecedented level of detail and 10 accuracy. The acquired point cloud can be handled conveniently with dedicated software, but applying traditional geomorphometry to analyse the cave surface is limited. This is because geomorphometry has been focused on parameterisation and analysis of surficial terrain. The theoretical and methodological concept has been based on twodimensional scalar fields which is sufficient for most cases of the surficial terrain. The terrain surface is modelled with a bivariate function of altitude (elevation) and represented by a raster digital elevation model. However, the cave is a three15 dimensional entity therefore a different approach is required for geomorphometric analysis. In this paper, we demonstrate the benefits of high resolution cave mapping and 3-D modelling to better understand the palaeohydrography of the Domica cave in Slovakia. This methodological approach adopted traditional geomorphometric methods in a unique manner and also new methods used in 3D computer graphics which can be applied to study other 3-D geomorphological forms. 20


Introduction
The caves are specific geomorphological forms typically developed in limestone rocks under appropriate hydrologic conditions.Caves contain various landforms such as speleothems or half tubes mostly organized vertically.However, more specific cave landform can form, such as ceiling channels.These channels incised in the cave roof evolve under a specific hydrological regime, which can be inferred from the channel morphology and associated features.Morphology of these speleoforms can indicate how the whole cave system developed and what the water stream parameters were at that time (Ford and Williams, 2007).Ceiling channels and other forms in the upper parts of the cave corridors are often difficult to reach and study in detail, especially, if the ceiling is several metres high.Direct cave surveying techniques such as mining compass, inclinometer, and theodolites have been traditionally used to map mainly the bottom part of caves (Mattes, 2015), but this approach is not suitable for mapping inaccessible surfaces.Therefore, groundbased remote sensing is the technological solution for such a case.Terrestrial laser scanning (TLS) is especially suitable since it uses its own source of electromagnetic energy, which makes it capable of mapping the surface with a high resolution in the darkness of the underground world (Buchroithner, 2015).Simpler methods exploiting laser distance measurers were used for capturing ceiling morphology (Lundberg and McFarlane, 2012), but the level of detail recorded does not reach the capabilities of TLS.
Laser scanning produces massive cloud of point data containing millions of three-dimensional (3-D) measurements of x, y, and z coordinates.The data can be visualized allowing for basic parameterization of the cave morphometry such as measuring distances, cutting cross sections, and colour-Figure 1.Schematic evolution of the passages in phreatic and vadose hydrological regime under conditions of high sediment flux and following abandonment (modified after Farrant and Smart, 2011).Further description of the stages is Sect.2.1.
ing.However, more information can be inferred if the point representation is converted into a continuous digital surface model.Afterwards, sophisticated methods for analysing the Earth surface morphology can be applied, which is the main domain of geomorphometry.Traditionally, a bivariate function has been sufficient to analyse the ground surface (Mitas and Mitasova, 1999).Digital elevation models (DEM) are usual representations of the Earth surface for which the methodological concept of digital terrain analysis is well developed and a plethora of geomorphometric tools exist in various kinds of software (Hengl and Reuter, 2009).In geosciences, 2-D arrays are the most popular representation of DEMs among the users.The main reason is due to the ease of manipulation and design of algorithms.However, the approach reaches its limitations for geomorphometry of 3-D volumetric surfaces such as caves (Roncat et al., 2011;Gallay et al., 2015b).
The aim of this paper is to present a novel methodology of deriving morphology information for the identification of ceiling channels using a high-resolution digital cave surface model and geomorphometric analysis.The high level of detail captured by TLS coupled with new 3-D modelling tools enable us to study specific cave landforms with unprecedented accuracy.This paper presents new aspects of studying cave genesis facilitating the traditional 2-D geomorphometric tools and innovative 3-D methods of computer graphics.The methodology is demonstrated with the data acquired by TLS in the Domica cave, Slovakia, in 2014 (Gallay et al., 2015a).This paper also provides new evidence of the hydrological regime acting during the cave formation.

Evolution of ceiling channels
Ceiling channels originate as a result of erosional action and corrosion of underground water flow, but they can evolve un-der various hydrological regimes (Ford and Williams, 2007;Pasini, 2009;Farrant and Smart, 2011).Most often there is initially a phreatic tube, which has formed along a guiding fracture in the rock (Fig. 1a).The phreatic tube grows in all directions (Fig. 1b).In case the seepage drops down, the phreatic regime changes to the vadose regime.The water stream is further incising the rock downwards.If the base level drops quickly, a keyhole passage is formed by the subsequent downward incision of the stream (Fig. 1c).The initial tube is abandoned and remaining as a half tube carved in the cave ceiling.Typically, the initial guiding fracture can be observed in the ceiling.A multi-phase evolution passage is formed if phases of base level drop and base stabilization alter (Fig. 1d).
Another scenario occurs, when the phreatic regime persists and the initial tube is filling up with sediments accumulated by the water flow but the post-depositional shrinking or erosion make enough space to sustain water flow between passage ceiling and the sediment fill (Fig. 1e).The sediment infilling protects the bottom and lower parts of the cave passage from further dissolution that is instead pushed upwards against gravity (De Waele et al., 2009).If the roof recession is compensated by sediment accumulation the upward incision continues (Fig. 1f), creating a series of typical morphologies such as ceiling half tubes (channels), wall anastomoses and half tubes, and pendants separating them (Fig. 2a-c).During such a speleogenetic process, the ceiling channel is often meandering without any relation to initial tectonic fractures (Fig. 2a).Further development of the ceiling half tube creates a paragenetic canyon (Lauritzen and Lauritsen, 1995).Presence of the specific passage morphology and scallop morphometry can be successfully used to distinguish between paragenetic and vadose regimes (Lauritzen, 1982;Lauritzen and Lauritsen, 1995).Change of the phreatic conditions to the vadose regime (e.g. by base-level rise) often causes evacuation of the sediments and therefore the paragenetic channels remain in the ceiling (Fig. 1g).Also the remnants of the sediments can be found near the ceiling vault as the evidence of the former accumulation fill.
First explanation of such formation of cave passages is attributed to Renault (1968), who termed the process as paragenesis.However, 3 decades earlier, Roth (1937) identified ceiling channels and explained their morphogenesis in the same way in his study on the evolution of the Domica cave in the former Czechoslovakia.Roth (1937) was published in Czech with a French summary, which limited its wider readership (Bella and Bosák, 2015).Pasini (2009) provided a detailed reasoning of the paragenetic formation of ceiling channels and suggest the more appropriate term "antigravitative erosion".Paragenesis is well-known and widespread in gypsum and limestone caves.Favourable environmental conditions for antigravitative erosion typically occur in the contact karst settings where there is sufficient sediment input into caves.Most sediment enters cave systems via stream sinks.Large, low gradient river caves are often subject to paragenetic development (Farrant and Smart, 2011).The sediment influx may be semi-continuous, with paragenesis occurring throughout the cave's development.The Domica-Baradla cave system is an example when allogenic drainage led to significant alluviation and paragenetic development (Bosák et al., 2004).The system evolved in the edge of a limestone block, which is in contact with gravel-clay deltalike sediments.Paragenesis is often thought to be of local interest, but increasing observations seem to indicate that many caves could be partially, if not almost completely, formed by this process (De Waele et al., 2009).Mapping the morphology of the cave ceiling thus provides essential evidence of the cave system formation.However, the morphological indicators are varying in size and abundance and they are often difficult to reach for humans.Therefore, recently emerging high-resolution mapping technologies, such as laser scanning, have potential to overcome the obstacles.

Benefits of terrestrial laser scanning and geomorphometry in cave mapping
TLS is an active remote sensing technique, which enables mapping of the cave surface with unprecedented levels of detail and accuracy.This technology has revolutionized cave mapping in the last decade (Canevese et al., 2009;Rüther et al., 2009;Buchroithner et al., 2011;Zlot and Bosse, 2014;Gallay et al., 2015a;Idrees and Pradhan, 2016).The scanner emits laser beams and records the reflected electromagnetic energy with a high frequency and accuracy.This generates a point cloud with x, y, z coordinates providing a very detailed 3-D record of the cave morphology (of the order of millimetres), which can be used for visualization, profiling, and measurements (Buchroithner, 2015).The points can be used for generation of 3-D cave surface models, which brings the possibility of detailed geometric analysis of the cave morphology.
Geomorphometry has traditionally been used in geometric analysis of terrain (land surface) (Evans, 1972;Krcho, 1973;Pike et al., 2009;Minár et al., 2013).Many morphometric parameters indicate likely occurrence of geomorphological processes, such as water erosion, landslides, or flooding (Hengl and MacMillan, 2009).However, specific forms such as caves or overhangs represent 3-D objects which require a more complex approach including new, 3-D modelling methods and tools.Traditional approximation of terrain surface with a bivariate function z = f (x, y) is limited for 3-D surfaces as there is more than one value of altitude z for a given location (x, y).Nevertheless, the bivariate approach can also be applied for modelling parts of the cave system, for example, cave ceiling or cave bottom.McFarlane et al. (2015) used this idea in identification of swiflets' nests in the Bomatong cave by extracting the upper part of the cave from a TLS data set and modelling the ceiling as a raster DEM.Another solution is to transform the orientation of the modelled TLS point cloud so that it satisfies the condition of unambiguity of the bivariate function.Starek et al. (2013) used the approach to model the lateral river bank erosion.Mahmud et al. (2015Mahmud et al. ( , 2016) ) worked with the distance from the scanner instead of the altitude to identify infiltration properties in limestone by the means of 2-D morphological analysis of the Golghota cave ceiling, Western Australia.However, it should be emphasized that generating the 2-D raster-DEM-like representations of the cave surface is possible only from a high-resolution 3-D point cloud representation of the cave acquired via laser scanning.Data surveyed by traditional cave mapping methods are not sufficiently detailed to be used in such a 2-D geomorphometric analysis.
The bivariate approach can be extended to a trivariate analysis of 3-D phenomena, where the scalar parameter w is modelled as w = f (x, y, z) (Hofierka and Zlocha, 1993).This concept is applicable for phenomena distributed in a 3-D space (e.g.air temperature, air pressure, water infiltration rate) and it is implemented in the geographic information system GRASS (Petras et al., 2015).However, the concept cannot be applied to 3-D landforms as the Earth surface is a 2-D boundary interface existing in a 3-D spatial domain.One of the solutions applicable for digital 3-D Earth surface modelling is to use polyhedral meshes (meshes).Polyhedral 3-D meshes are vector-based representations commonly used in computer graphics.Such a mesh is generated by connecting input 3-D points into a polyhedral network.Most often triangular facets are used which is similar to modelling the terrain with a triangulated irregular networks (TINs).
Currently, 3-D cave surfaces from laser scanning data need to be generated using other kinds of software, which have been recently developed in relation to the rapidly emerging 3-D data acquisition methods, for example, Blender, Meshlab, Polyworks, CloudCompare, 3-DReshaper, Geomagic Studio, VMTK.Primarily, the main application of the software is in modelling surfaces of objects other than the Earth's surface mainly for cultural heritage science, archaeology, architectural design, medicine, forensics, and civil engineering.This kind of software emerged in times of multiple-core processing units with a 64 bit architecture.Therefore, the algorithmic design is made to exploit parallel processing for fast computation and rendering of massive 3-D data sets even on a laptop.However, the tools can be applied to 3-D point clouds representing geographic landscapes.As the software is not specifically oriented to the analysis of cave morphology, specialized tools have to be used.For more details on methodologies for this issue, the reader is advised to consult works by Roncat et al. (2011), Jaillet et al. (2013), Cosso et al. (2014), Hoffmeister et al. (2015), Silvestre et al. (2015), or Gallay et al. (2015a, b), which exploited commercial software and open-source software for 3-D cave surface modelling.In this manner, the new software tools fill the gap which has emerged when the traditional 2-D geomorphometric approach adopted in GIS hit its limitations in handling 3-D surfaces.Various tools have been implemented for this task within the GIS platform, e.g.3-D analyst in ArcGIS (ESRI, 2015) or r3.modules in GRASS GIS (Neteler and Mitasova, 2008;Neteler et al., 2012), but still further research is needed to implement new 3-D analytical tools and to enable faster data rendering and processing.

Study site
The study area comprises the Domica cave situated at the south-western edge part of the Slovak Karst near the border of Slovakia and Hungary (Fig. 3; 48 • 28 40.4N, 20 • 28 12.9 E).The total length of the cave is almost 5400 m, however, its further continuation into the Hungarian Aggtelek Karst (the Baradla Cave) generates a multi-level single genetic cave system of 26 065 m (Gulden, 2016), which belongs to one of the longest in the Carpathians.The cave was inhabited by the Neolithic people but after a natural blockage of the entrance the cave had not been visited by a human until Ján Majko re-discovered it in 1926.Since that time, the cave had been widely studied from various research fields (Nováková, 2009;Papáč et al., 2014;Svitavská-Svobodová et al., 2015;Mihailović et al., 2015).The Domica-Baradla cave system is a listed UNESCO Nat- ural Heritage Site since 1997.Panorama photography views of selected parts of Domica are available as a virtual tour on http://www.panoramika.sk/firma/jaskyna-domica/10056/.Laser scanning of the cave provided non-disturbing means for highly detailed 3-D parameterization of the cave corridors enabling to study its morphology and acting as a catalyst for further research and educational initiatives.Domica is a typical example of a fluviokarst cave, which was formed by several underground water streams in limestones.The cave was formed by corrosive-erosive processes caused by superficial fluvial water and temporary streams, which sank underground at the contact of the Middle Triassic white limestones and the Pontian fluvial-lacustrine gravelsand-clay sediments (Bella, 2001).The tectonic framework was generated by south to southeast oriented pressures, which caused stress, and the compression was realized in the range of directions of north-south to northwest-southeast (Gaál and Vlček, 2011).Bella et al. (2014) performed cosmogenic nuclide dating of quartz pebbles cemented in the upper parts of Domica (340 m a.s.l.).Their results suggest that the upper level (338-340 m a.s.l.) began to form after uplifting the region above sea level earlier than the Middle Pliocene (3.47 ± 0.78 million years) when the current hydrographic network was being established.The middle level is located 10-12 m below the uppermost level.The lowest evo-lution level of Domica was found at one place by drilling at 318 m a.s.l.(Droppa, 1972).
The cave is a result of the underground flow of Styx and Domický potok shown by the oval shapes of the corridors and the quantity of allochtone pebbles found in the system.Both water streams are strongly influenced by the amount of atmospheric rainfall, which was observed according to the hydrological monitoring of water discharge of Styx and Domický potok inside the cave and simultaneously monitoring precipitation outside the cave during 1999-2001 (Peško, 2003).Average annual air temperature in the cave varies between 10-11 • C and air humidity is around 95-98 %.The average annual outdoor air temperature is between 8 and 8.3 • C and the mean annual precipitation varies around 635 mm based on monitoring during 1955-2004 at the closest meteorological station in Plešivec, 10 km north of Domica (Haviarová and Gruber, 2005).The water discharge of Styx and Domický potok is normally very low (0.4-0.5 L per second), however, it quickly increases during rapid rainfall or snow thaw when the soil is saturated or frozen.In the recent past, several rainfall events caused major flooding in the cave.This flooding was also due to inappropriate agricultural practise (Bella, 2001;Gaalova et al., 2014).
There were several phases in the cave's evolution when the streams were pushed towards the ceiling by depositing their sedimentary load in the cave what resulted in upward erosional incision of the stream, thus the formation of ceiling channels (Bella, 2000).The channels were first described by Roth (1937), who published an extraordinary study on the morphogenesis of the cave reconstructing its palaeohydrography (Bella and Bosák, 2015).His work was further extended by other followers (e.g.Kunský, 1950;Droppa, 1972), who mapped the cave and measured physically accessible parts of the cave.The last complex surveying was done by the Geological Survey, national enterprise, Spišská Nová Ves in 1975(Novoveský, 1975).However, the maps portray only the cave bottom not the ceiling.In places where the height of cave corridors restrains closer study of the ceiling, the geometric properties of channels were just estimated.In particular, the high ceiling could be only partially observed for obstacles or properly measured.The uniqueness of the Domica cave and the interest of wider research community in this area were the main driving factors to acquire a highly detailed 3-D representation of this underground landscape.

Terrestrial laser scanning
The data used in the presented research were acquired with a terrestrial laser scanner in combination with differential satellite positioning using real time kinematic method (RTK-GNSS) surveying within a 5-day mission in March 2014 in the Domica cave, Slovakia.Although the cave is drained by the permanent subsurface water flow Styx, the water level was very low and the cave was generally dry.This provided favourable conditions for scanning with an infrared laser.The survey is thoroughly described and compared with other similar surveys in Gallay et al. (2015a); hence we describe only the main facts.FARO Focus 3-D scanner was used to scan around 1600 m of the cave from 327 individual scanning positions within 40 h in total.The majority of individual scans were acquired at a scanning resolution given by point spacing of 7.6 or 6.1 mm at the range 10 m away from the scanner, respectively.In some instances the scanner was set to 3 mm spacing of measured points at 10 m range especially in large caverns with high ceilings.The final point cloud contained over 11.9 billion of measured points representing the entire show cave and some parts inaccessible by public.The total area of the scanned cave footprint was about 15 000 m 2 .
Semi-automatic registration of the individual scans was carried out using reference spheres placed in each scene and this method achieved an overall registration error of 2.24 mm (root-mean-square error, RMSE).Three RTK-GPS points surveyed outside the cave visitor centre and one point from older cave survey (Novoveský, 1975) in the rear north part of the cave were used to transform the final registered point cloud from its local coordinate system to the Slovak national cartographic system (S-JTSK) with the vertical ref-erence to the Baltic Datum after adjustment using the local geoid model.The transformation achieved the total accuracy of 21 mm (RMSE).The registration of scans and georeferencing was conducted in the SCENE8 software by FARO.The SCENE data project is used to manipulate the point data in the full spatial extent and full resolution.For other analytical tasks and 3-D modelling, the points are exported into the PXT format which is an ASCII based interchange format for point cloud data used by the Leica Cyclone software.The format stores normal vectors for each point oriented with respect to the scanning position from which the point had been located.The point cloud provides a very high detail of the order of few millimetres, which enables viewing even small geomorphological features such as soda-straw speleothems.

Generating a 3-D cave surface model
It is not convenient and efficient to analyse the entire georeferenced point cloud at the highest level of detail (full resolution).Either parts of it are selected and the selected points are further used for meshing or other analysis in full resolution or the original cloud is decimated to model the cave surface at lower resolutions.Despite some loss of information, the decimation (i.e.subsampling) homogenizes spatial distribution of the points into a semi-regular spacing thus providing control of the level of detail (spatial scale).First, for a more convenient handling, the original points were exported from the SCENE project into the PTX format by sampling every 7th point in horizontal and vertical direction with respect to their origin (i.e.scanner position).This procedure reduced the number of the exported points to almost 2 %.However, it was still sufficient to keep the level of detail required for the analysis presented in this paper (Fig. 4a).Afterwards, the cloud was subsampled in the open-source CloudCompare software (Girardeau-Montaut, 2006, 2015).Two data sets were derived by subsampling with the criterion of minimal distance between two points.This means that the points are not closer to each other than the specified value.The points were reduced having a minimal spacing of 1 and 5 cm, respectively (Fig. 4b).Before further use, the clouds were manually cleaned in Meshlab (Cignoni et al., 2008;Cignoni and Ranzuglia, 2014) from noise and outliers which originated mainly due to the presence of water.The new clouds were further used in the analysis of the cave surface morphometry based on 2-D raster representation of digital elevation models and 3-D meshes (Fig. 4c).
The acquired TLS point cloud provided a discontinuous (point) representation of the cave surface.For this reason it is necessary to generate a surface model.Computational representation of 3-D surfaces is a widely studied problem in computer graphics and disciplines where 3-D surface is the concern of research (e.g.biology, medicine, physics, geology).Surfaces are usually represented as a polygonal mesh, which comprises a collection of vertices, edges, and facets defining the shape of a polyhedral object (Tobler and Maierhofer, 2006).For 3-D modelling of the Domica cave, we adopted a triangular 3-D-mesh consisting of triangular facets.Each facet is defined by a set of three vertices and its orientation (the angle of azimuth and slope) is defined by a vector normal to the facet.Usually, it is meaningful to reduce the input point cloud before these steps are performed to test the methodological approach.
The Meshlab software (Cignoni and Ranzuglia, 2014) was used to generate a 3-D triangular mesh from the input point clouds.Meshlab is free and open-source software for mesh processing and editing capable of working with numerous 3-D file formats.Generation of a 3-D surface model requires that points are assigned with normal vectors.In our case, the points had normals assigned during scanning with respect to the scanner position.Otherwise, it would be necessary to calculate the normals for which several algorithms exist in the software.The surface model (mesh) can be reconstructed using several algorithms.We tested the Poisson surface reconstruction approach by Khazdan et al. (2006).Silvestre et al. (2015) used this procedure to model the surface of a cave chamber and identified stalactites based on local minima of the 3-D surface.The Poisson reconstruction of the 3-D surface is based on the observation that the normal field of the boundary of a solid can be interpreted as the gradient of the solid's indicator function.Therefore, given a set of oriented points sampling the boundary of a solid, a 3-D-mesh can be obtained by transforming the oriented point samples into a continuous vector field in 3-D.This is performed finding a scalar function, whose gradients best match the vector field, and extracting the appropriate isosurface.A thorough definition of the Poisson surface reconstruction can be found in Khazdan et al. (2006).The approach is also implemented in CloudCompare.It is important to mention that the vertices of the reconstructed triangular 3-D meshes do not coincide with the points of the survey points as it is in the case of the Delaunay triangulation.With this algorithm, octree depth is the key input parameter controlling the level of surface detail.It is the maximum depth of the octree hierarchy that is used for surface reconstruction by fitting an indicator function modelling the 3-D surface.The method adapts the octree to the sampling density, therefore the specified reconstruc-tion depth is only an upper bound.The number of vertices (also facets) comprised in the resulting 3-D mesh increases with the increasing value of the octree depth.After reconstructing the 3-D mesh, further processing and surface analysis was performed, such as removing non-manifold edges, duplicated vertices, decimation of the mesh, or the mesh parameterization.
The 3-D model for the full spatial extent of the scanned cave passages was generated using the Poisson surface reconstruction algorithm with the octree depth setting of 12.The data sets used originated from subsampling the original point cloud at 5 cm minimal spacing of points.The resulting model contained 4.77 million of vertices and it was stored in the binary Stanford polygon file format (PLY) having 349 MB (Fig. 4c).Selected parts of the cave were modelled in higher detail from the point cloud decimated to 1 cm minimal spacing of points and with the octree depth of 14.

Identification of ceiling channels using 2-D geomorphometry
Caves are subsurface landforms with a complex morphology.Despite the cave being a true 3-D landform, it can be partially modelled and analysed with traditional 2-D geomorphometry.The geomorphometric tools usually applied for 2-D gridded elevation data can be used to identify and measure particular features associated with the speleohydrological processes.For example, the 3-D point cloud can be filtered in order to extract the points of the highest elevation representing the cave ceiling.In the next step, DEM of the cave ceiling can be generated from the points.Similarly, points of the lowest elevation can be extracted to generate DEM of the cave bottom part.
To parameterize the ceiling and the bottom of the Domica cave, we generated continuous grid-based DEMs of the ceiling's highest altitude, the bottom's lowest altitude, passage height, and the depth of the cave ceiling below the terrain.Also, the ceiling channels were extracted as polylines by the means of traditional raster-based algorithms and map algebra (2-D geomorphometry).As the first step, the decimated and spatially homogenized point clouds were exported from CloudCompare in the LAZ format.LAZ is a compressed version of the LAS which is a public file format for the interchange of 3-D point cloud data data between data users defined by (ASPRS, 2014) format.In this way the storage demand decreases (Table 1) and the point cloud can be processed very efficiently in LAStools (Isenburg, 2014).
We used the open-source unlicensed lasgrid utility to generate the raster DEM of cave ceiling by extracting the highest altitude value of a point within a squared grid cell of 25 cm size.Similarly, the lowest altitude value of a point was extracted to generate the raster DEM of cave bottom.
The decimated point clouds of 5 cm spacing were used for this purpose.The DEMs represent the orthogonal projection of the cave morphology on a horizontal plane.The DEMs were imported into a GRASS GIS geodatabase (Neteler and Mitasova, 2008) where they were further handled by tools dedicated for raster data analysis.Relative height differences between the ceiling and the above-surface were calculated using a digital elevation model of terrain (DTM) generated from an airborne laser scanning (ALS) data set.The data originated during a mission flown in August 2014 by Photomap s.r.o., Košice over a wider area of the Domica cave.
The accuracy of measurement in open areas is reported at 0.1 m of 1σ by the data supplier.The density of laser returns classified as ground points varied between 0.5-8 points per m 2 which allowed for production of high-resolution gridded DTM of 1 m cell size.The workflow of processing the ALS data is more specifically described in Hofierka et al. (2016).
Tables 2 and 3 summarize the generated data and key parameters of the commands used.
The ceiling channels can be perceived as ridges of the cave closed volume.The algorithm used for calculating the ceiling line is based on the modified approach of Hardin et al. (2012), who extracted sand dune ridgelines.Their method is based on least cost path algorithm and traces the DEM cells of the highest elevation value between two points.We used a normalized ceiling height (DEM_CH_NORM) before calculating the cost surface model as the altitude values in our case are higher numbers (hundreds of metres) as opposed to the coastal dunes.By this means, we avoided calculation of very small cost values approaching zero which reach the limitation of the software at a certain level of decimal precision.As a result, the approach can be universally used in other kinds of landscape.The constant of 300 improves delineation of the line; however, the value can be changed ar-bitrarily.The higher it is the smoother line will be extracted.The output of this operation is called CEILING_LINE in Table 3.If channels stretch continuously across the cave roof algorithms for water flow routing could be also used.However, the course of ceiling channels in Domica is interrupted in many places by younger forms such as speleothems, tectonic ruptures, chimneys, etc.We tested the option of using r.watershed module in GRASS GIS after multiplying the DEM_CH with a constant of −1 which swapped the morphology upside-down.It provided means for comparison with the least cost path approach.However, the flow routing paths were sensitive to surface obstacles, therefore, we further worked with the least cost path algorithm.The ceiling line was delineated between several places and its geometry was parameterized by calculating length, sinuosity, and elevation summary statistics.The ceiling channels were also delineated as areal landforms by the means of geomorphometric classification using the approach of Jasiewicz and Stepinsky (2013), which is implemented in GRASS GIS by the module r.geomorphon.The method was successfully used to identify planation surfaces related to Domica-Baradla cave system in the Aggtelek Karst in Hungary by Veselský et al. (2015).The classification of the DEM_CH surface delineated the channels as ridges and some parts as summits.Also other convex features were assigned the same category although they are clearly not ceiling channels.To separate the channels form other ridge-like forms, additional processing of DEM_CH_CLASS was performed by extracting the ridges and summits with a Boolean logic and removing small clusters by a focal mode filter (r.neighbors).The remaining areas were clumped into individual objects (r.clump) of unique category value for which elevation statistics were calculated (r.stat.zonal,r.stat.quantile).Afterwards, a map algebra command was used to separate the ceiling channels from other ridge and summit forms based on similar median elevation of clumped objects.The remaining set of 77 ceiling channels objects was reduced to 26 after comparing their location with the 3-D cave surface model and observation in the cave.By this means areal segments of ceiling channels with similar morphometric properties were identified.

Analysing the cave surface using 3-D geomorphometry
The generated cave surface 3-D model enabled a new range of analysis parameterizing the cave morphology more accurately and from different perspective than with the DEM approach.The analysis was performed in Meshlab and Cloud-Compare where the suite of tools for analysing 3-D geometry is very diverse and customable.Volume and surface area were calculated for meshes generated with differing octree depths with the compute geometric measures filter in Meshlab.From the visualization point of view, just the sole viewing of the generated 3-D cave surface model is a useful means of inspecting the cave morphology in searching for ceiling channels and associated forms.Changing the light source angle reveals particular morphological features.However, more powerful measures than interactive shadowing can be calculated to improve perception of the 3-D shape and to help identifying the features.We explored various 3-D surface quality measures such as altitude, ambient occlusion and curvature parameters to colourizing the 3-D mesh based on the calculated parameter values.First, the mesh was colourized based on the z coordinate of its vertices (altitude) with the Per Vertex Quality Function filter in Meshlab.In this form, the mesh was used for interactive Web visualization described in Sect.4.6.Further, ambient occlusion was calculated with the Ambient Occlusion Per Vertex filter using default parameters.Ambient occlusion is the parameter, which improves visibility of micromorphology.It is a measure of illumination from a set of different light directions used to crudely approximate global illumination without simulating multiple reflections of light.The parts that are more occluded appear darker than the parts that are more exposed.Ambient occlusion was successfully applied not only in 3-D computer graphics for visualizing molecular surface (Tarini et al., 2006), objects of cultural heritage (Callieri et al., 2011), but also DEMs derived from airborne laser scanning data (Hinks et al., 2015).Curvature of a surface is a significant parameter commonly calculated in DEM analysis.In the presented study, we explored the multi-scale 3-D cave surface curvature parameterization with the filter for colorizing curvature of an algebraic point set surface (APSS) by Guennebaud and Gross (2007) filter in Meshlab.Gallay et al. (2015b) tested calculation of mean surface curvature at various levels of scale for a small section of the Domica cave.Their approach uses a local moving least-squares (MLS) fitting of algebraic spheres with intuitive parameters for curvature control of the fitted spheres.Computation of the curvatures requires points, mesh vertices, or facets equipped with oriented normal.The spatial scale of the fitted sphere can be controlled by setting the MLS filter scale, which is a value relative to the local point spacing of the vertices.The higher is the filter scale value the larger is the neighbourhood of the vertex for which the sphere is fitted and the curvature is calculated.The curvature value represents the inverted radius of the sphere tangent to the points in the defined neighbourhood of a vertex.We calculated the principal curvatures K1 and K2, mean curvature, and Gauss curvature.K1 curvature parameterizes the minimal curvature of the surface while K2 curvature defines the maximal curvature of the surface (Willmore, 2012).Gauss curvature is the product of K1 and K2 curvatures.The average of K1 and K2 values for a given vertex is the measure of the overall surface convexity/concavity, i.e. mean curvature parameter.As a result, generally convex features, such as stalactites, stalagmites, or pendants, have positive mean curvature values while concave features, such as sinks, cavities, half tubes, or sinks have negative values.Lai et al. (2014) used the curva-ture of a 3-D mesh as a measure of the surface roughness of rock faces.
Roughness of the 3-D cave mesh was parameterized with the roughness tool implemented in CloudCompare.For each point, the "roughness" value is equal to the distance between this point (mesh vertex) and the best-fitting plane computed on its nearest neighbours.The size of spatial scale is controlled by the kernel radius parameter, which is the radius of a sphere centred on each vertex.We tested several radius settings and finally used the radius of 0.5 and 5 m, which accurately depicted the anastomosing half tubes and the ceiling channels, respectively.

Results and discussion
The outcomes of the applied methodology demonstrate the capabilities of current state-of-the art in 3-D mapping technology and 3-D modelling software for research on 3-D hydromorphology of underground fluvial systems.We focus on explaining the benefits of the geomorphometric approach in studying ceiling of fluviokarst caves rather than providing interpretations of the generated data.The results of the presented research will be further structured in three sections according to the applied methodology in Sect. 4. Section 5.1 presents the outputs generated by 2-D geomorphometric analysis of DEM data sets.The results generated using the 3-D point cloud and 3-D cave surface model are presented in Sect.5.2.Collaboration of a geomorphologist/speleologist and GIS specialist is the key factor for efficient interpretation of the generated 2-D/3-D data sets, because in this case both disciplines require wider knowledge and experience to han-dle the data.Therefore, we generated a tool for interactive visualization and simple analysis of the generated data via the Web interface, which is presented in Sect.5.3.
The 3-D point cloud acquired by TLS is the key data source, which enabled exact quantification of the cave ceiling and from which we derived 2-D and 3-D models representing the cave surface.To familiarize the reader with the overall layout of the analysed Domica cave system, Fig. 5 provides perspective view of the 3-D cave surface model with its orthogonal projections.The presence of the ceiling channels can be clearly observed in the detail view as the elongated and horizontally winding half tubes in the cave roof.Also, we emphasize the fact that the top parts of the channels formed in horizontal levels, which, in some places, are disrupted by younger vertical or subvertical forms, such as chimneys, related to tectonics.

DEM-based analysis
The ceiling of Domica was modelled by a bivariate approach, which allowed for using diverse geomorphometric methods developed for 2-D digital terrain analysis with DEMs.The generated DEM data sets and DEM-derived data represented the cave ceiling morphology orthogonally projected on a horizontal plane (Figs. 6 and 7).Table 4 presents summary statistics of the gridded elevation models.The statistics are valid for the scanned part of the cave.
The footprint area of the cave based on DEMs is 1934 m 2 .The altitudes of the scanned cave system range between 322.51 m (minimum of DEM_CL) to 368.58 m (maximum of DEM_CH).The ceiling altitude ranges within 46 m but its height is influenced by the presence of speleothems or chimneys.Therefore, quantile values are more informative in assessing the ceiling altitude.Majority of the ceiling evolved within 328.92 to 339.93 m.However, two peaks of the statistical distribution can be observed in the histogram of the ceiling highest elevation (Fig. 6a).The lower maximum relates to the youngest evolution level of the Virgin passage developed at about 329 m a.s.l.(segment A in Fig. 8).The higher maximum is associated with the ceiling of the older evolution stages (segments B-F in Fig. 8) at about 338.75-340 m.The highest 10 % of the ceiling above 340 m is related to the concave features that originated as a result of younger ceiling collapse along tectonic fractures, which is well seen in side projections of the 3-D model in Fig. 5.
The extraction of the highest and lowest points of the original point cloud enabled calculation of a digital model of difference (DOD_C) as the relative ceiling height above the cave floor (Table 4).The morphology of the cave bottom (DEM_CL) is influenced by human interference when building infrastructure for the show cave, but the statistics provide useful data for comparison (Fig. 6b).The DOD_C model is more revealing as it quantifies the openness of the cave interior and the amount of the eroded material (Fig. 6c).The range of the heights for a certain location is 32.42 m and at 90 % of locations the cave ceiling less than 10 metres high.
The scanned cave passages have developed relatively deep below the terrain surface (Fig. 6d).Despite that, there are some chimneys almost 30 m tall, 90 % of the ceiling is more than 34 metres and 75 % more than 59 m below the terrain, respectively (DOD_CHT).No marked signs of the ceiling collapses on the above surface can be observed on the lidarbased DTM (Fig. 2).There are dolines north of the known cave passages but no connection has been proven yet.

Least cost path analysis
The cave 3-D model (Fig. 5) and the ceiling DEM (Fig. 6a) clearly reveal winding forms in the highest parts of the ceiling not noting the chimneys or stalactites.These are the ceiling channels, which are most markedly developed in the Dome of Mysteries, the Majko's Dome, and the Gothic Dome.The channels thus form ridges in the DEM_CH for which a ridgeline, i.e. ceiling line, was automatically extracted.Considering the ceiling morphology, we divided the ceiling into six contiguous segments marked A-F in Fig. 7 for which elevation, length, and sinuosity statistics were calculated (Table 5).Figure 8 shows the ceiling line and its longitudinal profile between the "dry passage" (Suchá chodba in Slovak) and Gothic Dome (segments B-D in Fig. 7).It can be seen that the profile is not smooth or flat as it would be expected in case of a generic river profile.The reason is that, besides the downwards flow of water, the formation of ceiling could be influenced by hydraulic pressure in the phreatic regime and paragenesis.After the water level decreased, and the conditions changed to vadose or epiphreatic, speleothems developed or the ceiling collapsed resulting in formation of chimneys.The ceiling line thus follows the DEM cells of the highest altitude.Nevertheless, the ceiling line very well indicates the presence of the ceiling channels.The longitudinal profile shows most of the parts preserved as horizontal levels.Interestingly in these flat parts, no obvious inclination of the profile line with respect to the flow direction can be identified.The same stands for other ceiling segments.The interruptions of the profile act as statistical outliers; therefore, quantiles are most appropriate to infer the trends in the elevation of the derived ceiling lines.
The segments can be grouped into three categories based on the median and interquartile range (IQR).The segment A (a.k.a.Virgin passage) has the lowest median elevation along the ceiling line of almost 330 m.It represents the youngest known ceiling level.Further, the segments B-D form another group having similar medians between 338.77 and 339.02 m and a relatively narrow IQR of less than a metre.These facts support they common speleogenesis related to older route of the Styx river.The third group comprises the segments E and F, which were formed by a different water flow  of the Domický potok river.Their median elevation is similar (339.78 and 339.95 m) and also it is the highest among the three groups.In comparison with the segment A, the segments B-F are mutually more similar and represent the same evolution phase.These fact were speculated in previous studies (Roth, 1937;Kunský, 1950;Droppa, 1972;Bella, 2000) but our approach provided exact measurement and evidence.Moreover, the level of meandering was quantified by the sinuosity index SI (Table 5).According to the conventional classes of the index, the ceiling lines are not straight, most of them are twisty which, however, cannot by attributed only to lateral meandering, but tectonics also affects the ceiling line.The most markedly meandering ceiling line is the D seg-ment, which has a nice marked meandering half tube channel.The meander wavelength gradually shrinks from 30 m in the Dome of Mysteries to 15 m in the Gothic Dome.

Ceiling surface classification
The visual and statistical analysis of the DEM representing a ceiling surface and the extracted ceiling lines revealed that the ceiling channels are areal speleoforms analogous to ridges.The ceiling line followed the highest parts of the cave ceiling comprising also forms which obscured the analysis of the actual channels.The current channels are remnants of more extensive ceiling channels which continuity was disrupted in subsequent evolution phases by tectonics, col- Table 5. Summary of elevation and length statistics of the ceiling line segments as shown in Fig. 8.   5.
lapses, and formation of speleothems.Their surface is gen-erally smooth and horizontal.To assess more specifically the differences within each of the ceiling line segments, we extracted the actual ceiling channels by a semi-automatic approach (Fig. 9).This analysis delineated those remaining channels, which are well identifiable by geomorphometric methods in GIS.The resulting set of the 26 channels was assigned with zonal statistics (Table 6).Most of the channels have much lower variation of elevation (SD) than in case of the ceiling lines (Table 5), indicating the smoothness and horizontal form of the channel surface.The median elevation of the channels most appropriately associates with the visual inspection of the 3-D cave surface models and observations in the cave; therefore, we classified the 26 channels into five distinctive categories accordingly (Fig. 10a).Their vertical displacement indicates vertical tectonic movement also suggested by formation of massive chimneys disrupting the hor-  izontal continuity of the channels (Fig. 10b).For example, the sequence of the channels 14, 19, 21 suggests common formation by the same palaeostream in the same time but the median of the channels 15 and 20 is 40-50 cm lower than the median of channel 22. Contiguity of the channels is disrupted by a tall chimney, which is horizontally perpendicular to the ceiling line and spans up vertically about 30 m. Similar cases exist in other parts of the cave (e.g.segments 21, 17, 16).They have been discussed since 1930s; however, no means for exact quantification were available until the emergence of high-resolution laser scanning and 3-D-modelling.

Cave ceiling analysis based on 3-D geomorphometry
The presented 2-D analysis provided diverse means of parameterization of the ceiling channels which can be sufficient for providing exact information on the cave formation.However, the analysis itself would be much more complicated without inspecting the 3-D cave surface model and comparing it with the results of the 2-D analysis.More importantly, the capability of parameterizing the 3-D surface provides new aspects for geomorphometry.We refer to such an approach of the Earth surface quantitative analysis by the   means of 3-D surface model as 3-D geomorphometry.The following sections demonstrate the results of applying 3-D surface modelling methods in geomorphometric analysis for identifying speleoforms associated with antigravitative erosion of the ceiling (Sect.5.2.1) and their parameterization at multiple levels of scale (Sect.5.2.2).

Enhanced 3-D model visualization
New possibilities for perceiving the cave morphology arise when viewing the 3-D virtual geomorphological object in an interactive fashion.Visual inspection of the features of interest from various angles and distances is very useful.In fact, it was the first tool that drew our attention to the ceiling channels.The possibility to look at the outside surface is an important benefit of the virtual 3-D model visualization as such a perspective is impossible to achieve in the real world.The ceiling channels are well depicted by the 3-D model when viewing from outside (Fig. 5).Moreover, interactive viewing of the 3-D model helped to identify morphological features associated with paragenetic formation of the ceiling channels, thus supporting existence of the process.After creation of the 3-D model, it is usually rendered in a single colour and shaded according to the orientation of surface normals with respect to the source of light.It is analogous to the 2-D shaded relief model calculated from a DEM.Instead of using the single colour, the rendering can be enhanced by defining the colours according to a certain attribute or quality.Figure 5 shows the model coloured according to the value of the z coordinate, thus elevation, and provides instant perception of the vertical position of the surface mor-phology.Ambient occlusion is another parameter which improved rendering of small-scale features and helped to reveal anastomosing half tubes based on the outside view of the model.The darker the tone the more occluded is the mesh vertex by the surrounding mesh facets.Figure 11 shows the part of the Gothic Dome as an example.Features indicating paragenesis such as anastomosing half tubes, pendants, and the distinctive meandering ceiling channel (half tube) are clearly visible from the outside view (Fig. 11a).It can seem obvious, but before seeing the model, these features remained hidden for a naked eye or obscured as they are situated high up, close to the ceiling.It is dangerous or too complicated to access a suitable position for appropriate observation.The detail view of the smaller area of the mesh (Fig. 11b) shows the cave walls from inside where pendants and the half tubes can be seen from the real-world perspective.The appearance is different, inverted.
The enhanced visualization also reveals progressive vertical and horizontal shifting of the meandering channels.Assessment of the direction of meandering provides another evidence for paragenetic or downward incision in vadose regime.Lauritzen and Lauritsen (1995) estimated the direction of the mean meandering vector and scallop morphology.In our case, no scallops were found and the morphology of the walls did not provide as many suitable places to measure the vector direction.Therefore, we used a simpler method of slicing the 3-D model perpendicularly to the z axis.The outer lines of the slices represents isolines of altitude, which were projected on the XY plane as a contour map (Fig. 12).In this way, the Dome of Mysteries (Fig. 12a) and the Gothic Dome (Fig. 12b) are portrayed.This cartographic method provides indices of the upward progressive incision of the ceiling channels in this parts of the cave.The horizontal shift of the winding contours altitude of the shift similarly as in the case of the meandering river.The contemporary position of the ceiling channel originated by lateral shift of the meanders as the arrows indicate.A river flowing in the same direction would migrate its meanders similarly.The altitude of the contours increases in the direction of the of the meander shift.Therefore, the meanders could develop only by antigravitative erosion thus incising upwards.The gradual downward widening of the passages, and absent initial guiding fractures in the ceiling, also indicate paragenetic development (Sect.2.1).

Multi-scale analysis of the 3-D cave surface
Tools for deriving metrics of a 3-D surface have been developed for the 3-D surface parameterization.We applied volumetric calculations and also derived parameters, which are analogous to the 2-D geomorphometry such as curvature and roughness.
The morphology of a natural 3-D surface is fractal to a large amount and the morphological features exist at multiple levels of scale.Therefore, the surface metrics are scale dependent.This aspect is well demonstrated by the analysis of the volume and surface area of the 3-D mesh in relation to the 3-D mesh resolution (Figs. 13 and 14).Results of the volumetric calculations need to be interpreted with care because the mesh volume and surface area change with a mesh resolution either by changing the octree depth of the Poisson surface reconstruction method or the number of input points.
The number of vertices and facets comprised in the resulting 3-D mesh increases exponentially with the increasing value of the octree depth.The time required for computation and the data size of the mesh in the PLY format are similarly related to the octree depth.On the other hand, the enclosed volume of the 3-D model decreases with increasing the mesh resolution.The area of the mesh surface rapidly increases to a certain level of detail when it tends to grow very slowly depending on the morphology represented by the mesh.We consider the volumetric parameters of the mesh with the highest resolution to be the most significant.Given its spatial extent the scanned cave surface area is about 54 300 m 2 and the mesh volume comprises approximately 54 050 m 3 .Multiple levels of scale can be explored by parameterizing morphometric properties for various local neighbourhood sizes for a mesh vertex.Such approach enabled one to control the level of scale and emphasized specific features apparent at a certain spatial resolution as sharply defined forms (Fig. 15a) or as smooth surfaces (Fig. 15a).Anastomosing half tubes, tectonic ruptures of the ceiling are clearly defined.The specific forms can be extracted based on a certain range of the parameter (Fig. 15c and d).
Surface roughness is closely related to the surface curvature.Figure 16 demonstrates the difference at two different scales.The areas of with the presence of small features such as the anastomosing half tubes and pendants appear rougher than the surface of the ceiling meanders, which is smooth at the small scale (Fig. 16a).The situation is inverted when assessing the surface at a larger scale (Fig. 16b).The roughness   values are higher for the ceiling meanders than for the areas with smaller features.The larger speleoforms clearly stand out.

Web-based tool for interactive 3-D visualization and analysis
Several tools have been recently developed enabling interactive 3-D visualization via a Web interface, thus providing appealing means for research presentation, dissemination, and further analysis (Evans et al., 2014).It is now possible to integrate 3-D content on the Web directly into the browser without plug-ins or additional components.For example, Silvestre et al. (2015) presented an approach in which X3-D, WebGL, and X3-DOM were used to enable online 3-D visualization and navigation of the interior of the Algar do Penico cave, Portugal, in several different Web browsers.Potenziani et al. (2015) introduced their 3-D Heritage Online Presenter (3-DHOP), which is an open-source software package for the creation of interactive Web presentations of high-resolution 3-D models.We used and modified the templates provided on the 3-DHOP website (http://3-Dhop.net/) to develop a tool for interactive visualization of the Domica cave 3-D model online via internet (Fig. 17).The tool is available at http://spatial3-D.science.upjs.sk/3-Dhop/SPATIAL3-D/index_domica_10x10_od13.html.The interface enables zooming, rotation, panning, changing the source light direction, and measuring Euclidean 3-D distances between two points.The model was generated in Meshlab from a reduced number of input points (3.13 million) with the octree depth of 13.Further use in 3-DHOP required conversion of the model into the compressed NEXUS format (http://www.vcg.isti.cnr.it/nexus/), which is based on a multiresolution data structure (Cignoni et al., 2005).The format allows the client to efficiently perform view-dependent visualization for faster streaming and smooth rendering in the browser.The size of the model was reduced from 148 MB in the PLY format to 20 MB after conversion.

Conclusions
The TLS technology facilitated exact and dense measurement of cave surface elevation.In combination with digital 2-D and 3-D modelling methods of surface reconstruction, it was possible to derive geomorphometric properties of the Domica cave in the West Carpathians, Slovakia.In particular, we have focused on identification and scrutiny of the cave ceiling in which the channels (half tubes) are formed.Previous studies of the Domica cave suggested a multi-level evolution of the cave with altering phases of water flow incision and sediment accumulation.As a result, there were stages of paragenetic (antigravitational) evolution during which the ceiling channels were formed.However, the previous studies were based on field observation, sediment sampling, and the conducted surveying mapped the bottom part of the cave and other physically accessible parts.The approach enabled a highly detailed and very accurate mathematical representation of the cave morphology for deriving new information on speleogenesis.Such knowledge has not been accessible prior to the advances in surveying technology and spatial data analysis.
The scientific contribution of this paper is in (i) defining a methodological approach for applying tools of 2-D geomorphometry on 3-D cave data, (ii) applying scale-dependent 3-D geomorphometric analysis of a cave, (iii) exact parameterization of the cave morphology and particularly the features associated with paragenetic formation of ceiling channels, (iv) developing a Web-based tool for 3-D interactive visualization and analysis of the Domica cave.
The ceiling channels are half tubes and the ceiling in general can be regarded as a bivariate scalar field of elevations at a certain level of scale if no overhangs are considered.Therefore, the 2-D gridded DEM representation can be used to approximate the surface of the cave ceiling.Similarly, the cave bottom can be modelled.This approach enabled (i) much simpler handling of the 3-D point cloud in the form of a DEM, (ii) exploiting the plethora of tools developed for geomorphometric analysis of DEMs in GIS software, thus (iii) deriving morphometric parameters of the cave surface.We successfully demonstrated calculation of basic morphometric parameters, scale-dependent landform classification, and least cost path approach to derive ceiling lines.This approach generated new findings on the geometry of the ceiling half tubes, their vertical position and spatial extent.
Caves are complex 3-D objects, therefore, the 2-D approach has limitations in capturing their morphology on the walls or in overhangs.We refer to the approach of the Earth surface quantitative analysis by the means of a digital 3-D surface model as 3-D geomorphometry.For this task, we applied existing tools generally developed for computer 3-D graphics to generate a 3-D cave surface model and to extract 3-D information from the model.This innovative ap-proach provided means for (i) dynamic 3-D interactive viewing of the cave system in high-resolution from arbitrary positions, (ii) revealing speleoforms associated with paragenesis such as anastomosing half tubes or pendants, (iii) estimating progress and direction of meandering of the ceiling channels.
From the data preparation point of view, we can conclude that it is important to scan at the highest possible resolution (point density) to be able to select the level of detail in the subsequent stages of the research.The original, highresolution point cloud has to be reduced and spatially homogenized for efficient computer processing but the original resolution serves as benchmark for assessing the quality of the generated model.
The methodology framework in our research can be applied in other caves with a similar configuration, but it can be also further extended.The recent trends in acquisition and processing of 3-D geographic data indicate increasing demand for development of tools capable of not only visualizing but also analysing morphology of 3-D objects of natural and urban landscapes.The future work in cave geomorphometry could exploit the multiscale dimensionality classification of 3-D point clouds (CANUPO) by Brodu and Lague (2012).From the regional perspective, the presented findings and data sets can be further used to link the cave morphology with surficial hydrography to estimate water discharge or other hydrological parameters during the cave evolution.The interpretation of the results still remains an open research area, but there are tools enabling to address questions with 2-D and 3-D geomorphometry.We hope this paper will also stimulate research on new methodologies and tools for analysis and modelling of massive 3-D point clouds in GIS.

Figure 2 .
Figure 2. Examples of ceiling channels (ch), pendants (p), and anastomosing half tubes (aht) in the Domica cave as speleoforms associated with the phreatic regime.Presence of the speleoforms in the ceiling and walls of the Gothic Dome indicate paragenetic (antigravitational) development (a, b, c), while it is more difficult to infer from the "dry passage" (Suchá chodba in Slovak) (d) from where much of the accumulated sediment was not evacuated and the passage was reshaped during multiple phases of vadose and phreatic conditions.The highest parts of the channels can be situated as high as 16 m (a, b) or 2 m above the floor on sediment accumulation (d).

Figure 3 .
Figure 3. Location of the scanned part of the Domica cave with orthogonal projection of the Domica-Baradla cave system according to Droppa (1972).

Figure 4 .
Figure 4. Point cloud from the Gothic Dome orthogonally projected on a horizontal plane representing 2 % of the full resolution as exported from the SCENE project (a), after subsampling the points at 5 cm minimal spacing (b), and as a continuous 3-D cave surface model (c) generated from (b) with the octree depth of 12.
raster values and their units DEM_CH DEM of the cave ceiling The highest altitude of points within a cell of 0.25 by 0.25 m, in metres.DEM_CL DEM of the cave bottom The lowest altitude of points within a cell of 0.25 by 0.25 m, in metres.DOD_C Model of the passage height Difference between DEM_CH and DEM_CL, cell size 0.25 m, in metres.DTM Digital terrain model Elevations of the surficial ground surface measured with airborne laser scanning, in metres.DOD_CHT Model of the depth of the cave Difference between DEM_CH and DTM, cell size ceiling below the terrain surface 0.25 m, in metres.DEM_CH_CONTOURS_05 Extraction of contour lines of the Vector contour lines of 0.5 m vertical interval ceiling extracted from DEM_CH.DEM_CH_NORM Normalized altitude of ceiling Cost surface raster based on normalized altitude values of DEM_CH, unitless.DEM_CH_NEG, Extraction of stream network from Modelling of water flow paths on ceiling model DEM_CH_NEG_STREAMS negative DEM_CH multiplied by −1.COST Cost surface model for extraction of a Cost surface raster calculated as exponential of ceiling line DEM_CH_NORM, unitless.CUM_COST Cumulative cost between two points Cumulative cost surface raster based on the COST_SURFACE raster, unitless.CEILING_LINE Extraction of ceiling line between A raster line calculated as the least cost path two points following cells of the highest altitude of ceiling, based on a cost surface calculated from DEM_CH, binary values of 0

Figure 5 .
Figure 5. 3-D model of the Domica cave oriented towards the north (a) and its detailed view towards north-west (b) showing a meandering ceiling channel.The model is coloured by the value of altitude above mean sea level which distribution is shown by the histogram.The 3-D geometry is orthogonally projected on the XY, YZ, and XZ planes.The size of the major and minor cells in the background grid is 20 and 5 m along all three axes, respectively.Note the flat, horizontal ceiling levels projected on the vertical XZ and YZ planes.

Figure 6 .
Figure 6.Raster based models with histograms representing the highest cave ceiling altitude (a), lowest altitude (b), vertical difference between the cave bottom and ceiling (c), and the ceiling depth below the terrain surface (d).The histogram y axis represents the number of raster cells in hundreds.

Figure 7 .
Figure 7. Segments of ceiling lines with assigned letters as referenced in the text and as reported in Table5.

Figure 8 .
Figure 8. Longitudinal profile along the automatically delineated ceiling line of overlaid on top of the raster DEM_CH.

Figure 9 .
Figure 9. Classification of the ceiling surface (DEM_CH) of the Majko's Dome and the Gothic Dome (Fig.3) in GRASS GIS to extract the ceiling channels as landforms.Ceiling channels are classified with the r.geomorphons module as ridges and some parts as summits (a).The two categories are extracted as binary rasters with r.mapcalc and small clusters are removed with r.neighbors (b).The areas are clumped into individual spatial units for which elevation summary statistics were calculated.Median elevation (c) most appropriately separated the ceiling channels from other ridges and summits (d).

Figure 10 .
Figure 10.Ceiling channels resulting from the semi-automatic classification of the Domica ceiling, which were grouped according to their median elevation (a) reported in Table 6 and interpretation of their vertical displacement (b).

Figure 11 .
Figure 11.Ambient occlusion of the 3-D cave surface model in the Gothic Dome rendered with grey tones enhances perception of the cave surface morphology (a).Detail view of the smaller area of the 3-D mesh (b)shows the cave walls from inside.The size of the major and minor cells in the background grid is 1 and 0.5 m, respectively.The displayed mesh was generated from the original TLS data set subsampled at 1 cm point spacing demonstrating the high resolution of cave mapping.

Figure 12 .Figure 13 .
Figure 12.Elevation contours of the 3-D cave surface projected on a horizontal plane indicate upward migration of the ceiling meanders in the Dome of Mysteries (a) and in the Gothic Dome (b).Contour interval is 0.5 m.

Figure 14 .
Figure 14.The level of detail for different 3-D meshes generated from the same input data with various octree depth parameter values of the Poisson surface reconstruction method.Example of the Gothic Dome using the input point cloud subsampled from the original at 5 cm minimal spacing.

Figure 15 .
Figure 15.Mean curvature parameterization of the 3-D cave surface model in the Gothic Dome (Fig. 7a) at the scale of 5 (a), scale of 20 (b), selecting the surface of mean curvature between −8 to −3 m −1 from (c), selected anastomosing half tube (d).The size of the major and minor cells in the background grid is 5 and 1 m, respectively.Orthogonal silhouette of the mesh is projected on the XY, XZ, YZ planes.

Figure 16 .
Figure 16.Roughness of the 3-D cave surface model in the Gothic Dome derived at 0.5 m scale (a), and 5 m scale (b), and selection of speleoforms having values between 1 and 2 (c), respectively.

Figure 17 .
Figure 17.Print-screen snapshot of the Web interface allowing for interactive viewing and basic measurements of the 3-D surface model of the Domica cave.

Table 1 .
Summary of point cloud reduction.

Table 2 .
Summary of raster data derived from 3-D point cloud in LAStools and GRASS GIS.

Table 3 .
Computer commands and their key parameters for generating the data listed in Table2.

Table 4 .
Summary statistics of the DEM data sets parameterizing morphology the scanned part of the Domica cave.

Table 6 .
Ceiling channels identified in the DEM_CH classification procedure and their summary statistics of elevation and channel halfwidth.