Nathalie Debese

Jean-José Jacq

Thierry Garlan

Extraction of sandy bedforms features through geodesic morphometry

Keywords: Bedform, MBES, TIN, Anisotropic filter, Features extraction

State-of-art echosounders reveal fine-scale details of mobile sandy bedforms, which are commonly found on continental shelfs. At present, their dynamics are still far from being completely understood. These bedforms are a serious threat to navigation security, anthropic structures and activities, placing emphasis on research breakthroughs. Bedform geometries and their dynamics are closely linked; therefore, one approach is to develop semi-automatic tools aiming at extracting their structural features from bathymetric datasets. Current approaches mimic manual processes or rely on morphological simplification of bedforms. The 1D and 2D approaches cannot address the wide ranges of both types and complexities of bedforms. In contrast, this work attempts to follow a 3D global semi-automatic approach based on a bathymetric TIN. The currently extracted primitives are the salient ridge and valley lines of the sand structures, i.e., waves and mega-ripples. The main difficulty is eliminating the ripples that are found to heavily overprint any observations. To this end, an anisotropic filter that is able to discard these structures while still enhancing the wave ridges is proposed. The second part of the work addresses the semi-automatic interactive extraction and 3D augmented display of the main lines structures. The proposed protocol also allows geoscientists to interactively insert topological constraints.

Introduction

General context

Currently, multibeam echosounder systems (MBES) are routinely used to map fine-scale seafloor features of regional areas. By collecting dense datasets with high-resolution and accurate soundings, these acoustic sensors confirmed the ubiquity of subaqueous dunes in continental shelves and coastal areas. Since such bedforms were first described by Lüders in 1930[START_REF] Lüders | Sediments of the North Sea[END_REF], the later generation of MBES reveals a wide variety of complex sandy bedforms. Coupled with high accuracy positioning and inertial motion sensors, repeat MBES surveys are useful to measure bedform migration [START_REF] Garlan | GIS and mapping of moving marine sand dunes[END_REF].

Mobile bedforms represent vital issues for hydrographic offices as they potentially reduce water depths. Understanding sand dune formation processes is required to maintain navigation channel conditions [START_REF] Whitmeyer | Sand waves that impede navigation of coastal inlet navigation channels[END_REF]. After discovery, their migrating nature requires hydrographers to resurvey such zones to avoid navigational hazards. Because MBES surveys are costly and time consuming, the criteria for area prioritization have to be defined from the types and sizes of their dynamics [START_REF] Dorst | Improving a bathymetric resurvey policy with observed sea floor dynamics[END_REF][START_REF] Garlan | Marine dune dynamics[END_REF]. Knowledge of the behaviors of sand waves is also critical for the engineers in charge of offshore installations. Routing procedures have to be optimized as mobile bedforms may damage cables or pipelines [START_REF] Morelissen | Mathematical modelling of sand wave migration and the interaction with pipelines[END_REF]. With the substantial expansion of marine renewable energy platforms, it is crucial to predict sand wave migration as scouring may weaken established marine systems [START_REF] Németh | Offshore sand wave dynamics, engineering problems and future solutions[END_REF][START_REF] Kenyon | Sand banks, sand transport and offshore wind farms (UK Government[END_REF][START_REF] Barrie | Seabed characterization for the development of marine renewable energy on the Pacific margin of Canada[END_REF]. Moreover, in European coastal waters, mines from WWII buried in sandy environments still make the installation of platforms difficult [START_REF] Garcia | Mine burial by local scour and sand waves[END_REF]. However, sandy bedforms are of considerable interest from an economical point of view. The exploitation of marine aggregates, especially in Belgium, is increasing [START_REF] Eynde | Modelling the effects of sand extraction, on sediment transport due to tides, on the Kwinte bank[END_REF]. Nevertheless, intensive dredging activities have sedimentary and morphological impacts on marine environments due to altering the forms and volumes of the sand dunes [START_REF] De Mol | Which future for the sand extraction in the Belgian part of the North Sea? Belgium Pier -Blankenberge[END_REF].

The above examples clearly illustrate the importance of sandy bedforms for various end-users. Dunes are dynamic selforganized systems, modeled by complex interactions between flow, sediment transport and bedform [START_REF] Powell | Wave-current interaction over bedforms: observations and model predictions[END_REF]. A wide range of theoretical and numerical models have been developed to explain the generation, evolution and destruction of bedforms in both laboratory and field contexts. Up to now, processes controlling bedform dynamics still remain far from perfect [START_REF] Franzetti | Giant dune morphologies and dynamics in a deep continental shelf environment: Example of the banc du four (Western Brittany, France)[END_REF]. Model limitations are due to the complex behavior of sand dune patterns that form from nonlinear processes operating in open systems [START_REF] Werner | Modeling landforms as self-organized, hierarchical dynamical systems[END_REF]. The reduction and simplification of the physical processes involved in bedform evolution lead to morphodynamic models that are far from reality. Furthermore, the spatial variability of field observations and the scarcity of empirical analysis add further concerns to this scientific challenge [START_REF] Garlan | GIS and mapping of moving marine sand dunes[END_REF]. Because both sizes and shapes of bedforms influence and reflect hydrodynamic and sedimentological processes, another relevant approach to modelling sand dune dynamics aims to characterize their morphologies [START_REF] Lefebvre | Quantification of roughness for seabed characterisation[END_REF]. Our algorithm belongs to the latter class, as it aims to objectively describe and analyze bedforms acquired using MBES. Its originality lies in the use of a mesh data structure to represent the bedforms as surfaces, their analysis through geodesic morphometry and quantitative description using 3D parametric curves.

Related works

Sand dunes are generally described as a train of features with regular shapes, sizes and spacings. Deeper observation reveals that sand bedforms display much more complex behaviors [START_REF] Werner | Modeling landforms as self-organized, hierarchical dynamical systems[END_REF]. With the exception of Seif dunes, MBES surveys revealed subaqueous features mirroring aeolian dune shapes on a smaller scale [START_REF] Garlan | GIS and mapping of moving marine sand dunes[END_REF]. Dunes observed on the seabed can be isolated structures or organized in dune fields. Their sizes, based on spacing, cover a wide range from ripples, megaripples, and sandwaves to dunes [START_REF] Ashley | Classification of large-scale subaqueous bedforms: A new look at an old problem-SEPM Bedforms and Bedding Structures[END_REF][START_REF] Perillo | Hacia una nomenclatura y clasificación de formas de fondo transversales generadas por flujos de agua[END_REF]. The lengths, orientations and geometries of their crest lines show wide variations from linear ridges to crescent shapes [START_REF] Dyer | The origin, classification and modelling of sand banks and ridges[END_REF]. Because they are shaped by complex flow dynamics that operate on different spatial and temporal scales, sand dunes often show multiple superimposed bedform scales [START_REF] Barnard | Small-scale sediment transport patterns and bedform morphodynamics: new insights from highresolution multibeam bathymetry[END_REF][START_REF] Gutierrez | Discrimination of bed form scales using robust spline filters and wavelet transforms: Methods and application to synthetic signals and bed forms of the Río Paraná, Argentina[END_REF].

Currently, there is no consensus on how to model the morphologies and dynamics of sand bedforms [START_REF] Werner | Modeling landforms as self-organized, hierarchical dynamical systems[END_REF]. Nevertheless, numerous studies of both flume and field experiments have linked some morphological features to environmental conditions. Crest lines are generally perpendicular to the main current directions [START_REF] Garlan | GIS and mapping of moving marine sand dunes[END_REF]. Bedform height and wavelength generally increase with flow velocity [START_REF] Gutierrez | Discrimination of bed form scales using robust spline filters and wavelet transforms: Methods and application to synthetic signals and bed forms of the Río Paraná, Argentina[END_REF]. Dunes shapes, such as their degree of asymmetry or the sinuosity of their crest lines reflect hydrodynamic conditions: Bedforms with a high degree of asymmetry or those with sinuous crest lines are characterized by high migration rates [START_REF] Garlan | GIS and mapping of moving marine sand dunes[END_REF].

From these results, it makes sense to propose a geometrical reading of these bedforms to describe the local hydrodynamic and sediment dynamic processes. Up to now, geometrical features of bedforms have mainly consisted of local minima and maxima extracted from manual one-dimensional cross-sections that are drawn perpendicularly to crest lines. Manual analysis of bed elevation profiles is a tedious and time consuming task that relies on subjective measurements. This method is clearly not suitable for processing large volume of MBES datasets while taking full advantage of the detailed bathymetric information they contain.

Traditional analysis of bedforms can be improved through automatic selection of bathymetric profiles and automatic localization of their crests and troughs. When bedform profiles are not directly recorded [START_REF] Cataño-Lopera | Characterization of bedform morphology generated under combined flows and currents using wavelet analysis[END_REF][START_REF] Huang | Extraction of morphometric bedform characteristics from profiling sonar datasets recorded in shallow coastal waters of China[END_REF], they can be extracted as rows/columns of Digital Terrain Models (DTMs) provided that they are parallel to the direction of the tidal principal axis [START_REF] Knaapen | Sandwave migration predictor based on shape information[END_REF]. DTMs that do not fulfill this condition are rotated to place dune crest lines perpendicularly to raster scanning. [START_REF] Duffy | Patterns of morphometric parameters in a large bedform field: Development and application of a tool for automated bedform morphometry[END_REF] deduces the suitable rotation angle in a set of values as the one that minimizes the mean length of the dunes, while [START_REF] Cazenave | Quantitative bedform analysis using decimetre resolution swath bathymetry[END_REF] extract the bedform orientations based on 2D spectral analysis. Sand dunes are usually composed of amalgamated forms; analysis of their shapes requires the separation of their bedforms at different scales before quantifying the variabilities of their geometry. Other studies propose a bedform tracking tool to automatically find crest and trough locations from profiles smoothed using a weighted moving-average technique [START_REF] Knaapen | Quantifying bedform migration using multi-beam sonar[END_REF][START_REF] Van Der Mark | A new and widely applicable tool for determining the geometric properties of bedforms[END_REF][START_REF] Duffy | Patterns of morphometric parameters in a large bedform field: Development and application of a tool for automated bedform morphometry[END_REF]. Considering that all of the scales of sand dunes are of interest, other approaches use 1D-spectral analysis to retrieve bedform profiles of different wavelengths in the bathymetric signal and apply a detailed analysis of each Fourier-filtered bedform [START_REF] Kheiashy | Geometric and statistical characteristics of bed forms in the lower Mississippi River[END_REF][START_REF] Winter | Spectral analysis of compound dunes[END_REF][START_REF] Singh | Experimental evidence for statistical scaling and intermittent in sediment transport rate[END_REF][START_REF] Van Dijk | Quantified and applied sea-bed dynamics of the Netherlands Continental Shelf and the Wadden Sea[END_REF]. Bathymetric profiles often present significant departures from Fourier analysis assumptions. [START_REF] Cataño-Lopera | Characterization of bedform morphology generated under combined flows and currents using wavelet analysis[END_REF] propose the use of a wavelet technique to capture bedform components. Following this approach, [START_REF] Gutierrez | Discrimination of bed form scales using robust spline filters and wavelet transforms: Methods and application to synthetic signals and bed forms of the Río Paraná, Argentina[END_REF] describe bedform profiles within a hierarchical structure by combining robust spline filters to continuous wavelet transforms. As in previous morphological approaches, the geometric descriptors they propose are based on the crest and trough locations in bedform profiles. Different approaches exist to automatically detect these feature points. Crests and troughs can be selected from extremal depths obtained through bedform height and length thresholding [START_REF] Parsons | Morphology and flow fields of three-dimensional dunes, Rio Parana´, Argentina: Results from simultaneous multibeam echo sounding and acoustic Doppler current profiling[END_REF] or with the help of zero up crossings and down crossings (Van der [START_REF] Van Der Mark | A new and widely applicable tool for determining the geometric properties of bedforms[END_REF]. Different approaches exist to define and name the geometric characteristics of bedforms [START_REF] Jerolmack | A unified description of ripples and dunes in rivers, in: Taylor and Francis Group[END_REF][START_REF] Van Der Mark | On modeling the variability of bedform dimensions[END_REF][START_REF] Knaapen | Local variation in the shape of superimposed bed forms as a function of local bathymetry[END_REF]. Among the most common variables, the bedform height and length, crest and trough elevations and lee face slope are statistically analyzed to determine simple generic relationships between bedform geometries and find correlations with flow and sediment features [START_REF] Van Der Mark | Quantification of variability in bedform geometry[END_REF].

One-dimensional approaches can incorporate spatial aspects by taking into account different profile locations. [START_REF] Cataño-Lopera | Characterization of bedform morphology generated under combined flows and currents using wavelet analysis[END_REF] propose a one-dimensional wavelet transform to separate dominant wavelengths within bedform profiles. The juxtaposition of one-dimensional results enables the representation of regions with characteristic wavelength bands and analysis of their spatial-interaction behaviors.

Nevertheless, 1D approaches still remain inadequate to fully describe 3D planform patterns of bedforms. Studying the spatial distribution of co-existing bedforms of multiple scales requires a spatial approach. [START_REF] Cazenave | A method for semi-automated objective quantification of linear bedforms from multi-scale digital elevation models[END_REF] and Lisimenka and Rudowski (2013) apply a 2D Fourier transform to quantify the variability of linear bedform geometry displaying their orientation and wavelength using a vector field. [START_REF] Lefebvre | Quantification of roughness for seabed characterisation[END_REF] also apply a two dimensional spectral analysis to identify and map bedform regions according to their roughness. In addition to the Fourier decomposition technique, [START_REF] Van Dijk | Separating bathymetric data representing multiscale rhythmic bed forms: A geostatistical and spectral method compared[END_REF] propose a geostatistical filter to separate bedforms components. The crest and trough lines are then extracted as a set of points of extremal curvatures from the filtered bathymetric surface. To overcome the limitations of the 2D-Fourier transform, [START_REF] Garcia | Characterization of bed morphodynamics using Multibeam Echo Sounding (MBES) and Wavelet Transform (WT) analysis[END_REF] extend the one dimensional wavelet transform developed by Cataño-Lopera et al. (2009) to the spatial case. The morphological map leads to identification of superimposed regions at different bedform scales.

Methodological choices and paper contents

All of these studies contributed to a significant increase in our knowledge of such complex environments. Nonetheless, much remains to be done. Recent MBES surveys acquired near San Francisco and in the Irish Sea that contrast with previous observations suggest improvements to both the geometrical attributes and their associated extraction tools [START_REF] Hanes | The genesis of an inter-field marine sandwave and the associated anti-asymmetry migration of neighboring crests[END_REF][START_REF] Van Landeghem | Reversed sediment wave migration in the Irish Sea, NW Europe: A reappraisal of the validity of geometry-based predictive modelling and assumptions[END_REF]. The algorithm described below automatically extracts crest and trough lines of sandy bedforms modeled as a Triangular Irregular Network (TIN) built from MBES soundings. This algorithm requires the availability of a novel anisotropic filter. This filter, which is applied first, enhances the dominant features of bedforms while smoothing sand ripples. The resulting line features are described as 3D parametric curves and their associated domains are then subjected to deeper statistical analysis after optional expertknowledge integration. This paper is organized as follows: Section 2 describes the anisotropic filter applied to bathymetric mesh surfaces to enhance their dominant features while smoothing out their ripples. Section 3 is devoted to the automatic extraction of structural lines from filtered bedforms as 3D parametric curves. Section 4 describes the optional semi-automatic incorporation of expert knowledge into the structural lines before addressing the global statistics of their orientations. This section also covers some usual topics: accuracy, robustness, computational cost, and comparison with other line detectors. The main lines of a prospective application aiming at automatic extraction of the wave geodesic domains are then described before Section 5 that presents conclusions and perspectives.

2 De-noising and enhancement of bedform shapes through anisotropic smoothing

Data content and differential scale space

Obtaining a prior and insightful representation of the incoming data with respect to the underlying data processing framework is a key point to enable end-users to set parameters in accordance with the distribution of features. In this context, while implicitly addressing the analysis of a discrete surface, this subsection provides a short reminder of the main 2Drepresentations of the curvature pairs. For an exhaustive analysis of curvatures estimation techniques based on TIN, readers are referred to [START_REF] Gatzke | Estimating curvature on triangular meshes[END_REF]. There are two main non-hybrid families: continuous and discrete techniques. The latter have a very low-cost but provide crude estimations, whereas the former are cumbersome but provide optimal estimations. This work makes use of both types during different steps: Section 2 emphasizes the discrete approach, whereas Section 3 uses the classic continuous approach.

Let us recall relationships that exist between the common curvature pairs before emphasizing a more convenient choice: the differential scale space of Koenderink. Fig. 1 depicts the relationships between the common curvature pairs. • The principal curvatures (κ 1 ,κ 2 ) , with κ 1 ≥ κ 2 , and their orthogonal principal directions (e 1 ,e 2 ) , are the fundamental characteristics of the local shape configuration. These are the native outputs of the continuous methods. • The mean and gauss curvatures (κ H ,κ G ) can be seen as the reciprocal form of (κ 1 ,κ 2 ) where κ H = (κ 1 + κ 2 ) 2 and κ G = κ 1 κ 2 . Thus, the inverse relationship is

(κ 1 ,κ 2 ) = (κ H ± (κ H 2 -κ G ) 1
2 ) . The pair (κ H ,κ G ) is a native output of the discrete methods. • A more intuitive shape partitioning criteria can be built from (κ D ,κ H ) where κ D is the curvature mean difference:

κ D = (κ 1 -κ 2 ) 2 .
It divides the local shape domain into two half-spaces H + | H -with respect to the sign of κ H . • A major intuitive representation arises when the polar coordinates embedded in (κ D ,κ H ) are made explicit through the curvature pair (κ C ,κ P ) where

κ C = (κ D 2 + κ H 2 ) 1 2
is the curvedness, also defined as

κ C = (2κ H 2 -κ G ) 1 2
, and κ P = arctan κ H κ D is the angular phase defined on [-π/2,+π/2] . The normalized phase 2 π κ P , also known as the shape index [START_REF] Koenderink | Solid shape[END_REF], is thus defined on [-1,+1] .

The implicit (κ D ,κ H ) and explicit (κ C ,κ P ) polar coordinates are both well suited to get an intuitive overview of the morphological bedform variability while performing the anisotropic filter discussed below. However, in this data processing context, it is more convenient to make direct reference to the curvature radius: λ = κ -1 . Moreover, the datasets under study lead to the use of a logarithmic scale for the curvedness to take advantage of the Gaussian-like distribution of its logarithm. Thus, bedform shapes are described using 2D-histograms of the warped coordinates (λ D * ,λ H * ) and (λ C * ,κ P ), where

λ {D,H } * = λ C * λ C λ {D,H } , λ C * = ln(max(λ 0 ,λ C ) λ 0 )
, and λ 0 is defined as the minimal curvature radius. As stated by [START_REF] Koenderink | Solid shape[END_REF], who was the first to put forward the differential scale space ( 2 π ln κ C , 2 π κ P ) , it is worth emphasizing that the logarithm of the curvedness provides a more natural measurement scale.

In this work, the bathymetric surfaces are assumed to be provided with an extrinsic orientation. Hereafter, inner and outer

H + H - H ++ H +- H -+ H -- κ D κ H κ C κ P + π 2 - π 2 κ 12 = H + : κ 1 H - : κ 2 ⎧ ⎨ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ κ 21 = H + : κ 2 H - : κ 1 ⎧ ⎨ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ Fig. 1.
Relationship between common representations of a curvatures pair through 2D diagrams: principals κ 1 and κ 2 , mean sum κ H , mean difference κ D , gauss κ G , curvedness κ C , and angular phase κ P . The graphic design emphasizes the polar viewpoint.

regions bounded by a surface are assumed to be linked to sediments and water materials, respectively; the normal vectors are made outward-oriented. Thus, any estimated mean curvature will exhibit positive values on bedform ridges and negative values in bedform troughs or valleys. This paper focuses on ridges and valleys, i.e., structures where κ P ≈ ±π/4 . Fig. 2 gives a global overview of the test area used in this study. This bathymetric survey was carried out by the Hydrographic and Oceanographic Service of the French Navy (SHOM) in May 2005 using a SIMRAD/EM3002 MBES. Some examples of 2D-histograms shown in Fig. 3 are built from this area. These 2D-Histograms provide access to a qualitative description of these still unprocessed bedform shapes. A close look at Fig. 3a shows that some ridges are very sharp and that bedform stoss sides and ridges are superimposed with sharp-ripples waves. The existence of ridges and valleys -the former being sharper as expected -clearly comes out from 2D-histogram where the phases ±π/4 are over-represented. This accumulation effect is mainly the signature of ripples because they are statistically much more present than sand dunes. Moreover, the amplitudes of the ripples on several dunes are noticeable with respect to the dune elevation. Thus, the first step of the processing framework will need to eliminate these ripples without performing an excessive smoothing of the dune crests. To this end, the next subsections propose a filter that is not only able to smooth-out these ripples and the remnant noise, but also able to enhance the main characteristics of the bedform structures.

Surface isotropic smoothing

Let S(u,v,t) denote a 2-manifold surface evolving with re- spect to time t and embedded in the 3D-Euclidean space R 3 through a parameterization x(u,v,t) over the bounded do- main Ω . It has been shown that performing homogenous smoothing of a surface of some given amount is equivalent to letting it undergo an isotropic diffusion process during a given processing time [START_REF] Koenderink | The structure of images[END_REF]. In this framework, the initial state S 0 = S(u,v,t 0 ) is allowed to evolve under the mean curvature normal flow through the PDE ∂ t S = -κ H n , where κ H (u,v,t) is the mean curvature and n(u,v,t) is the unit normal vector. The Dirichlet conditions are currently supplemented by holding the domain borders ∂Ω frozen over time.

A popular and very general way of managing the spatial component of this PDE is to redefine S as the zero-set of an embedding function in the Euclidean plane R 3 and then to operate in this plane through well-posed numerical schemes involving relationships such as 2κ H = divn . However, when- ever the evolving surface is not allowed to undergo topological changes, which is the case for current application context, operating directly on the surface points will provide a much more efficient approach. In this way, [START_REF] Desbrun | Implicit fairing of irregular meshes using diffusion and curvature flow[END_REF] propose to directly compute the normal curvature vector as -2κ H n = lim φ A →0 ∇A A where A is a local infinitesimal ). Both plots are drawn with λ 0 set to 1 m. These 2D histograms are built on a regular grid involving 256×256 accumulation bins. A point footprint is computed as an isotropic Gaussian kernel with a 1.5-bins standard deviation and is defined on a 4-bins radius support. area of S around x which has a vanishing diameter of φ A . This vector can be seen as the result of the Laplace-Beltrami operator H(S) applied to a point x on the surface S -an operator that is also known as a generalization of the Laplacian from flat spaces to manifolds. Hereafter, the same symbol H refers to both the operator and its resulting vector. The anisotropic diffusion PDE becomes ∂ t S = 1 2 H(S) . However, whenever the application context requires a signed estimation of κ H , the normal vector will still be needed to perform the reverse factorization of H through κ H = -1 2 sgn(nH) ||H || .

H +- H + H - H ++ H -- H -+ + π 2 -π 2 λ C * κ P + π 4 -π 4 0 H + H - H ++ H +- H -+ H -- κ P + π 2 -π 2 λ C * λ

Anisotropic smoothing and shape enhancement

The methodology proposed below can be seen as a very low cost attempt to transpose the anisotropic approaches available in the flat space framework to the surface geodesic space. These approaches mainly operate through the variational minimization of an ad hoc functional and then give rise to a PDE managing a gradient descent. Intuitively, once transposed in the current application context, such a minimization process is expected to end up with an entire domain Ω spanned by a set of sharp ridges and valleys interconnected by a minimal surface. However, because the current processing objective is confined to providing a filter that is able to enhance dune shapes while filtering out both the ripples and the noise components, the proposed PDE drops any penalty management with respect to S 0 and turns to a mere iterative anisotropic quadratic regulari- zation. Thus, its main end-user parameter is the integration time.

The basic idea is to drive the local diffusion process ∂ t S = 1 2 H(S) through a scalar control function g(u,v,t) defined on [-1,+1] . For the sake of clarity, in the first step, the following discussion is limited to the quadrant H + . However, as the problem is symmetric with respect to axis κ D , an a pos- teriori transposition of the result to cope with the ridges mirrored on the quadrant H -(i.e., the valleys enhancement) is straightforward. As a basic rule, the locations where the curvedness is small enough still undergo linear diffusion, i.e., g is then set to unity. Otherwise, whenever a local configuration bears some resemblance to a ridge profile (i.e., showing both a noticeable curvedness and a vanishing Gaussian curvature) the speed of the diffusion process decreases, converging to zero for an ideal ridge. However, this adaptive iterative process must be asymmetric to drop shapes in stable configurations. If the vanishing Gaussian curvature is positive, the evolution process is damped by g while remaining diffusive. Conversely, if the vanishing Gaussian curvature is negative, the function g re- turns a negative coefficient letting inverse diffusion occurs. The remaining challenge is to define relevant thresholdsmainly, through an asymmetric transfer function with respect to the space on both sides of the κ 1 axis and a curvature magni- tude threshold along the same axis -to automatically balance between these dynamic modes.

Such an innovative and very low-cost non-linear PDE is proposed by Meyer et al. (2003) and validated by addressing the reverse engineering of noisy man-made shapes. However, unlike the application goals put forward by these authors, when addressing the current bathymetric datasets, some sharp bumps or holes may become frozen. To eliminate these artifacts or unwanted singularities, a less crude variant of their control function is described below. Let κ floor > 0 denote a curvature threshold. The control function of the new non-linear diffusion PDE is then built by cascading, in the given order, four logical cases: (iv) else g H + = 1 . Furthermore, to counteract the emergence of excessive curvatures, (v) a post-processing step progressively restores the linear diffusion while κ 1 in- creases toward a given threshold κ ceil > κ floor , and it then sets g H + = 1 if κ 1 > κ ceil . This is currently done through the update g

(i) if κ 1 ≤ κ floor then g H + = 1 , (ii) else if κ 2 > κ floor then g H + = 2κ floor (κ 1 + κ floor ), (iii) else if | κ 2 | ≤ κ H then g H + = κ 2 κ H ,
H + ← s 2 + (1 -s 2 )g H + where s ∈ [0,1] is linearly defined as s = (min(κ 1 ,κ ceil )-κ floor ) (κ ceil -κ floor )
. This new threshold should be setup with respect to the magnitude λ 0 -1 , i.e., with respect to a natural curvature bound linked to the spatial resolution of the MBES data. An overview of the action domains of g H + is depicted in Fig. 4.

If the mirror form g H -is required, it may be simply obtained by applying the swap (κ 1 ,κ 2 ) → (-κ 2 ,-κ 1 ) to these rules. However, because the non-linear PDE will still make an explicit reference to the normal diffusion vector H , g H + remains valid for the two quadrants and its counterpart g H -is not needed. Therefore, hereafter, g will implicitly denote g H + and its curvature coefficients will work simply through κ H ≡ 1 2 ||H || . In applications operating on flat spaces, it is well-known that diffusion flows evolving through an anisotropic PDE may be ill-conditioned if the operators involved in the control function are directly applied to the current state of the evolving shape. The correct numerical scheme requires applying these operators on an isotropically smoothed state of the current state of They are defined with respect to the thresholds κ floor and κ ceil and the angular phase defined through the ratio κ 1 κ 2 = -3 . The full diffusion still occurs in do- mains (i), (iv), and (v). g H + depends only on κ 1 in domain (ii), whereas it depends on both κ 1 and κ 2 in domain (iii). Sub-domain (iii+) slightly smooths the bump-biased ridges, whereas sub-domain (iii-) performs a slight inverse-smoothing on saddle-biased ridges. Data samples currently mapped on the edge between domains (iii+) and (iii-) (i.e., perfect ridges) remain strictly frozen if they are not too close to κ ceil .

κ D κ H κ 1 κ 2 κ f lo o r κ c e il + κ f lo o r - κ f lo o r (i) (ii) (iii+) (iii-) (iv) (v)
the shape [START_REF] Catte | Image selective smoothing and edge detection by nonlinear diffusion[END_REF]. In the same way, the control function discussed above has to make reference to a slightly smoothed state of S t , denoted by ! S t , obtained through the integration of the diffusion PDE over a small time interval ! t . The anisotropic PDE is ∂ t S t = 1 2 ! g t H(S t ) , where ! g t is now implicitly operated through the standard curvatures pair

(! κ H , ! κ G ) issued from ! S t , with ! S t = S t + 1 2 H(S ′ t ) ′ t =t t+ ! t ∫ d ′ t .

Discrete anisotropic diffusion

Let the TIN surface S denote a discrete counterpart of the 2- manifold S topologically built as a homogeneous simplicial 2-complex S = (V,E,T ) , where V , E , and T , denote the set of vertices, edges, and triangles, respectively, embedded in R 3 through a direct parameterization defined by the point cloud dataset {p i , i = 1,!,|V |} . Hereafter, p x and p y are the georeferenced coordinates of the MBES soundings while p z is their depths. The 1-skeleton G = (V,E) is said to be the graph of S . The k-ring of a vertex v is the set N G k (v) of the sur- rounding vertices linked to v through at least one path involv- ing at most k edges of G .

In this discrete context, the challenge is to swap from a native pointwise definition of an operator to another one based on a spatial average. Many works -mainly in the Computer Graphics domain -have been devoted to devising a discrete counterpart H for the operator H -see also Gatzke and Grimm (2006). However, [START_REF] Wardetzky | Discrete Laplace operators: No free lunch[END_REF] have shown that, from a theoretical viewpoint, all of the known discrete strategies focusing on the estimation of the normal curvature are deemed to be endowed with some intrinsic drawbacks. Meanwhile, the cotangents-based operator proposed by [START_REF] Pinkall | Computing discrete minimal surfaces and their conjugates[END_REF] and Meyer et al. (2003), currently a de facto standard that is summarized in Fig. 5, is expected to be the most versatile with respect to the current goal. The next points give the main grounds of this choice.

• Unlike many schemes that operate only on the 1-skeleton G , this operator implicitly operates through the 2-complex S , providing intrinsic estimations of H , i.e., its estimation vanishes for planar 1-rings for any relative locations of the vertices involved. Thus, a discrete diffusion flow based on this operator should not exhibit noticeable tangent drifts, which is a mandatory requirement for performing enhancements of ridges and valleys through non-linear diffusion because triangles may have to become badly shaped at these features.

• The computational cost of this scheme is very low. In addition, its vertices neighborhood is bounded to the 1-ring. This feature makes it suitable for efficiently processing an iterative integration of a PDE on a large TIN such as those built from MBES data acquired in shallow waters. • Some computational components involved in H provide, still at low cost, an additional discrete estimation of the curvature κ G through the application of the Gauss-Bonnet the- orem (see Fig. 5).

The discretization with respect to t is done through the for- ward Euler numerical scheme S t+δt = S t + 1 2 δt ! g t H(S t ) per- formed over the integration time t = mδt through m itera- tions -each of the ancillary smoothing steps being performed over ! t = ! mδ ! t . This scheme easily enables addressing non- linear PDEs. However, the limitation of such an explicit scheme is well known: At every time, δt must be kept suffi- ciently small so that 1 2 δt ! g t ||H t || does not exceed the square of the spatial granularity, i.e., a magnitude range ∝ λ 0 2 . Moreover, this theoretical step bound only holds for linear diffusion on regular grids. Thus, despite the internal regularization process, it is expected to have to proceed through a more conservative δt setup to cope with non-linear diffusion cases. It may follow that a larger number of iterations can be required to match the same space scale effect. However, because it addresses a filtering task, the current process only expects shape evolutions at the small space scale. Thus, the required integration time remains tractable.

Surprisingly, while processing datasets where λ 0 ≈ 1 m and starting from well-conditioned TIN, it is even found that working through δt = 0.5 , with δ ! t = 0.1 and ! m = 2 , still provides good results (see Fig. 6), thus saving a large amount of processing time. With this current setting of (κ floor ,κ ceil ) , the filter turns out to be almost isotropic with respect to the quadrant H -. Almost all of the ripples disappear while sharp (i.e., well localized) dunes ridges still emerge from a smooth background. Nevertheless, some ripple stamps break the continuity of the dune ridges because these intrinsic characteristics of some dune ridges can also be frozen by the filter. Thus, the selected setup of the main parameters of the filter (κ floor ,κ ceil ) is expected to provide an optimal trade-off between, on the one hand, the acceptable side-effects of the filter and, on the other hand, the robustness expected from the automatic features extraction process described below.

p i - p i p i + n = θ i + normal(T i + ) i ∑ || θ i + normal(T i + ) i ∑ || ⌢ A i ± = 0 if α i ± > π/2, 1 4 area(T i ± ) if (α i ± + θ i ± ) < π/2, 1 8 ||p i -p|| 2 cot α i ± else. ⎧ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ∇A i ± = 1 2 (p i -p)cot α i ± p ∇A = (∇A i -+ ∇A i + i ∑ ) ⌢ A = ( ⌢ A i -+ ⌢ A i + i ∑ ) κ G = (2π - θ i + i ∑ ) ⌢ A H = ∇A ⌢ A Mid point (peripheral edges)
Orthogonal (incident edges)

Contour of the influence area of 5. Computation of the normal curvature vector H through application of the discrete Meyer-Desbrun operator on the 1-ring of a vertex located at p . This graphic makes use of a circum differential indexation of the 1-ring to enhance the fragment-wise contributions in the analytical expressions. The symbol A denotes any influence area of p whose perimeter Γ always passes through the mid points of the incident edges. The actual path, denoted ⌢ Γ , follows the Voronoï influence area of p only within non-obtuse triangles. Conversely, to remain consistent with respect to the partitioning of the neighboring 1-rings, obtuse triangles require the contextual definition of their path segments. After being defined through the application of the Gauss-Bonnet theorem, the Gaussian curvature κ G at p is also provided. When the polarity of the mean curvature returned by H is required, a common way to compute the unit normal n at p is also given.

⌢ A ⌢ Γ Mid point (incident edges) α i - α i + θ i -θ i + T i - T i + ⌢ A i - ⌢ A i + p Fig.

Automatic extraction of salient linear structures as 3D parametric curves

Tracking of ridge and valley lines

The typical approaches for extracting ridge and valley lines from a TIN follow the same strategy. They track the zero-crossing of the derivative of the most significant curvature coefficient along its principal direction. These approaches are mainly of two types: global and local. A well-posed implementation of a global approach using implicit surface fitting is described by [START_REF] Ohtake | Ridge-Valley lines on meshes via implicit surface fitting[END_REF]. The local approaches, which mainly use explicit polynomial fitting, have a lower cost compared to global ones. Two popular local approaches are described by Cazals and Pouget (2005a) and [START_REF] Yoshizawa | Fast and robust detection of crest lines on meshes[END_REF], respectively. Because of the use of higher order derivatives, these approaches are very sensitive to noise. Therefore, they require isotropic smoothing of the surface prior to or during curvature estimations. In addition, principal curvature orientations become undefined at the umbilics; therefore, special care must be taken at these locations. The differential geometry classifies the ridges into several types. The salient ridges are called crest ridges. The current study focuses only on crest lines. Thus, for the sake of simplicity, subsequent references to the terms ridge and valley will implicitly refer to a crest ridge and crest valley, respectively. These standard methods are expected to provide new insights when the filter discussed in the previous section is included in their preprocessing pipeline. However, the output of these methods continues to be a mere bag of lines. Based on its filter breakthrough, the second part of this study now attempts to propose an efficient semi-global approach that can robustly provide geoscientists with a more structured output. This significantly different strategy will natively detect the main structural lines as part of an arborescent network. The strategy operates without the noisy principal directions and without performing any explicit derivation on the principal curvatures. The strategy consists of three steps. First, a scalar ridgeness indicator is built from an ad hoc analytical recombination of the curvature coefficients. Second, skeletoniza-tion fronts simultaneously start from the contours of regions where the indicator is greater than a predefined threshold. Third, they evolve simultaneously in an entropic manner with local speeds inversely proportional to the gradient modulus of the indicator level and with the constraint of not breaking the topology of the uncovered parts of their starting regions. The process stops when the fronts are frozen, and the remaining uncovered parts form the skeleton forest.

Overview of the extraction procedure

The availability of an efficient filter that is able to enhance ridge and valley structures, out the ripples, opens the door to their explicit extractions through a robust and automatic process. Its basic building blocks are high level transforms of the Mathematical Morphology -namely, reconstruction, skeletonization, tree pruning and hierarchization, influence zones, and constrained watersheds. Their theoretical foundations and their usual implementations are detailed by [START_REF] Soille | Morphological image analysis: Principles and applications[END_REF]. The entire process described below operates through its extension on the TINs. Because the internal topology of a TIN is not implicit, this choice requires additional resources to manage their structures. However, unlike the working case of the discrete regular Euclidean flat space (e.g., DTMs), where the spatial relationships are commonly described through lattice structures, the Jordan theorem still holds on TIN discrete 2-manifold. This behavior makes it easy to define vertex-based robust predicates enabling the performance of conditional erosion while preventing connected components from being split. From a methodological point of view, such a generic framework allows for addressing the segmentation of functions defined on 2-manifolds and can thus be called Geodesic Morphometry.

From an algorithmic viewpoint, the transformations listed above involve, either explicitly or implicitly, to manage the propagation of entropic fronts. Because it is based on powerful and efficient data processing structures -mainly Priority Queues (PQ) -such a dynamic model can address the pseudo-interactive analysis of very large datasets.

The aim of this work is to track both ridge and valley lines of the main sand waves. From a mathematical point of view the distinction between ridges and valleys of a surface is merely conventional . Because λ 0 = 1 , the anisotropic filter is thus locked on the curvatures pair ( κ floor , κ ceil ) = (0.2, 1.0). It is iterated with its generic internal setting -i.e., δt = 0.5, δ ! t = 0.1, and ! m = 2. The filtering scale m is set to 100 iterations.

depending on the current surface orientation. Thus, the generic process devised to track ridge lines still holds while tracking valley lines. Unless otherwise noted, the following discussion will now focus on ridges. These lines are extracted as a discrete 1-manifold network where explicit bifurcations are thus forbidden. The secondary objective is providing end-users with an automatic estimation of the volume of any sand wave; thus, the geodesic influence domain of each ridge line is made available. Although it is still a medium-term goal, some interesting results are already available and discussed in Section 4.6.

The entire procedure is summarized as a signed scalar indicator is made available to measure the saliency of the lines (Section 3.3). Then, the line extraction process is performed in four steps (Section 3.4). (i.) A mask including ridge lines is first extracted by thresholding the saliency indicator. To prevent the extraction of closed paths, any mask hole disconnected from the support borders is filled. Furthermore, to ensure compliance with the entropicwise processing context, every local minima of the indicator over the mask domain are simultaneously flooded through a geodesic reconstruction step. (ii.) This main step carries out the skeletonization of the indicator function over the mask domain, thus giving rise to a forest defined as a subgraph of the 1-skeleton G of the surface S . (iii.) Then, an ancillary step builds a structured and hierarchical description of each tree while operating through a geodesic hypermetric. As a result, each tree is broken into a primary branch (i.e., a trunk maximizing its geodesic path length) and a set of secondary branches, most of which are then pruned. (iv.) In this last step, the remaining branches are individually postprocessed. This task involves resampling, smoothing, and reparametrizing each line path. This enables access to some relevant characteristics and insertion of their augmented iconographic instances, built as generalized cylinders, into the interactive 3D-navigation scene.

Saliency indicator

This indicator provides a scalar signed value to measure how well a filtered location is linked to a sharp ridge, and a sharp valley, through large positive and large negative magnitude values, respectively. It is expected that (i) applying this indicator to the enhanced shape provided by the non-linear filter, and then (ii) taking it as the major input of the segmentation algorithm summarized above will enable the extraction of reliable high level characteristics. Let (κ C , κP ) denote the curvedness and angular phase pair computed from the pair (κ H , κG ) esti- mated from a shape Ŝ enhanced with respect to the ridges and the valleys. The saliency indicator is defined as K = κC sgn κP | sin 2κ P | n , where the strength of the smooth angular windows is currently defined with n = 4 . Hereafter, it is simply referred to as the ridgeness indicator, because it is obvious that its reverse encoding ⌣ K ≡ -K becomes the valley saliency indicator.

The pair (κ H , κG ) required to build this indicator is already available as a by-product of the filtering process. However, because the goal is now to perform a static shape analysis, it is relevant to perform this measurement more homogenously using a less crude estimator than the discrete Meyer-Desbrun operator. Thus, a classic continuous approach, operating through the fitting of an osculatory paraboloid on k-rings, is retained to obtain (κ 1 , κ2 ) first. Optionally, a regularization effect can now be managed while operating estimations with k > 1 . A mapping of this indicator on the enhanced surface discussed above is depicted in Fig. 7.

Extraction and representation of ridge lines 3.4.1 Preprocessing steps

Let M 0 denote a mask function defined on the vertices set V with M 0 = 1 if K ≥ h else M 0 = 0 , where h is a ridgeness threshold. To allow the end-user to easily express some relevant threshold h with respect to the actual distribution of the ridgeness function, this threshold can, in turn, be indirectly specified as h = K med + h * K mad through a rank-based scale h * , where K med = med i∈Ω +{K (i)} , and K mad = 3 2 med i∈Ω +{|K (i)-K med |} , are the median, and the median absolute deviation of K , respectively -both being computed over Ω + = {v i ∈V;K(i) ≥ 0} . Each Connected Component (CC) of M 0 receives a unique label identifier. Any CC of the inverse mask M 0 not connected to the domain bor- der ∂Ω is erased. This hole filling step guarantees that the linear structures extracted below will not contain circuits. Then, any CC of M 0 whose geodesic area is less than a CCmin is removed -the computation rule of a vertex area is described in Fig. 5. This cleaning step, which removes only sparse and quasi-punctual components, should not be too aggressive. Let M denote the mask issued from these preprocessing steps. The last preprocessing step performs a geodesic reconstruction of K max M over K , with K max = max i∈Ω {K(i)} . This is an iterative global transform applied until convergence. An iteration performs an erosion of K max M on the 1-ring neigh- borhood; if the eroded state of a vertex becomes lower than its reference value in K , its state is increased to this reference value. This transform results in a new reconstructed ridgeness map K R = rec K (K max M ) where any local minimum of K found in locations covered by M are now flooded. The parts of K linked to the background of M are left unchanged in K R . An example of this pair (K R ,M ) is depicted in Fig. 8a,b while focusing on one CC.

Skeleton extraction

The skeletonization of M with respect to its informative con- tent K R is done iteratively while performing successive condi- tional erosions on M until convergence. In the initial state, any vertex of M is labeled UNPROCESSED if it belongs to the mask foreground or labeled PROCESSED otherwise. The initialization step locates all of the foreground contour vertices of M and then pushes them into a waiting queue PQ with a priority equal to the negative of their ridgeness in K R . The fronts then start to progress iteratively while removing the current output vertex from PQ, whose label becomes PROCESSED. If this vertex location can be eroded with respect to its 1-ring without breaking the topology of the remaining part of its CC in M , then its mask state is zeroed. Independently of the results of this test, every neighbor of this vertex that is still UNPROCESSED and not already in PQ, is put in PQ with a priority equal to the negative of its ridgeness in K R . The process ends when PQ becomes empty. Let M SK denote the final state of M . The remaining foreground vertices of M SK form the skeleton.

However, to proscribe the occurrence of skeleton paths involving border edges, the vertices of the binary mask located on ∂Ω require a special care: At these places, the erosion step may be enabled to break the current CC topology giving rise to additional connected components within the skeleton. Thus, a final cleaning step occurs to discard from M SK spurious con- nected components with geodesic areas lower than a SKmin . The skeleton result of the tutorial CC is depicted in Fig. 8c.

Skeleton structuration into branches and pruning

Each tree of the skeleton issued from the previous step may involve many branches. Most of them account for insignificant sub-structures or for ripple residuals. Thus, a structural and hierarchical decomposition of each tree into a primary branch (i.e., its trunk) and a set of secondary branches is required to post-filter these trees while pruning some of their branches. The iterative process performing this task operates through a geodesic hypermetric based on (K R ,M SK ) . It can be described through three main steps: (i.) The initialization step marks each leaf of the skeleton with a unique label; the labels already involved in the CC identifications in M SK are reserved. Each of these vertices is registered as the tail of a starting path. At startup, these vertices are also registered as the current head of this incoming path whose initial geodesic length is thus zeroed. Each of these paths is natively labeled ONGOING and increments the ONGOING paths count of their respective CC by one unit. The vertices linked to these head locations are put in a priority queue PQ. At every step, the negative of the cumulated length of the tail-head path segment will define the priority of its current head vertex in PQ.

(ii.) The branch recovering process iterates while PQ is not empty. The current output is removed from PQ. Let i and l i denote the identifier of this vertex in S and the current path length to reach its own tail, respectively. Its processing context depends on the ongoing path count n(CC i ) of its CC. If this count is lower than two, this loop reiterates immediately. Otherwise, this unloaded vertex i is always the current head of an ONGOING path. The vertices of the undiscovered parts of the skeleton (i.e., those unprocessed that are not already registered as part of one path) found within the 1-ring neighborhood of i are then analyzed. The objective of this scan is to determine a topologically valid successor j so that the current path can progress through an additional step of length l ij . This prospective step length is quantified through a geodesic hypermetric

l ij = 1 2 (K R (i) + K R (j)) ||p j -p i || .
The topolog- ical test applied to validate j states that removing the ver- tex j must not break the topology of the unprocessed part of the skeleton. The actual processing of j still depends on n(CC i ) .

• The count n(CC i ) is greater than two. If a valid successor j has been found for i , then j becomes the new head of the path and is pushed into PQ with priority -(l ij + l i ) . Otherwise (i) the path is labeled TERMINATED and the ONGOING paths count n(CC i ) is decremented by one unit, and (ii) i is registered as the final bifurcating head of the path and receive the label of its native CC in M SK . • The count n(CC i ) equals two. Because the preprocessing step (Section 3.4.1) ensures that there are no circuits in the mask, these two remaining ONGOING paths are the two opposite half paths of the trunk of CC i . Thus, the current step always contributes to the completion of the trunk and the two remaining head points should be made to progress until they collide -this event also defines the geodesic midpoint of the primary branch. In the current context, i gets always a valid successor j . If j is not the vertex head of its half-path counterpart (i.e., j is not in PQ), the vertex j is pushed in PQ with the priority -(l ij + l i ) and the loop reiterates immediately. Otherwise, this event marks the completion of a trunk path. The two halfbranches are concatenated to form the primary branch and the ONGOING path count n(CC i ) is zeroed. This primary branch has no head and the current labels of its two tails are erased with the same number: the label of its native CC in M SK . However, a pseudo-head is still required for technical purposes, this one is made to point to the tail vertex that has the greater in S . (iii.) In this final pruning step, only secondary branches are subject to clipping. Those whose length is shorter than l BRmin are removed. Unlike the results of the two cleaning steps already described, this one should give rise to noticeable changes because many small secondary branches are indeed spurious structures. Its optimal setting encompasses expert-level trade-off (see also the discussion in Section 4.1). Hereafter, M BR denotes the cleaned and branch- structured state of M SK , where paths are now labeled with the label of their respective tail vertex. The branch-wise decomposition of the skeleton is depicted in Fig. 8c-e. The complete set of branches is displayed in Fig. 8d, where each branch is mapped with a unique color expressing the label of its tail. Only two branches -the trunk and one secondary, which are overlaid by their respective smoothed paths -remain after the pruning example shown in Fig. 8e.

Branch re-parameterization and augmented 3D representation

Each ridge line is now registered in a database as an ordered list of vertices of S describing a polygonal geodesic path along with its tail-labeled description in the mask M BR . To allow the user to perform a perceptual and quantitative analysis on the basis of a less crude set of paths, each path is then smoothed and reparameterized with respect to its curvilinear abscissa defined on [0, 1]. Only the two ends are kept frozen at their native vertex locations. A curvilinear differential analysis maps an intrinsic direct orthonormal Frenet trihedron (i.e., the tangent, normal and binormal (t F ,n F ,b F ) unit vectors-triplet along with its (curvature, torsion) scalar pair) on every locations of the path. This homogenous description enables the construction of a generalized cylinder around each path.

A last cleaning step is operated at this GUI level. To this end, a median value of K R is computed separately along each regis- tered path, and the median of this set of medians is then re-turned. The paths with median ridgeness that are lower than this value are discarded from 3D representation if their tail-head cord length is less than l THmin . This heuristic filter is mainly ex- pected to address the removal of quasi-circuits around bumps. The two paths depicted in Fig. 8e have successfully undergone this infographic-level cleaning step.

Any parametrized line can also be endowed with extrinsic attributes to provide additional visual insights. The latter can be some characteristics inherited from those found under the original discrete geodesic path (e.g., the surface normal n or the saliency indicator level K R ) as long as they have been subject to the same post-processing steps as their native vertex locations. The color texture of the cylinder encodes the amplitude of an extrinsic channel whereas its radius section is simultaneously modulated by another extrinsic channel. In the scene depicted in Fig. 8f, the underlying ridge saliency of the surface is now visually transferred over the cylinder, both through its color texture and its section radius.

Fig. 9 displays the network of ridges extracted from the Lézardrieux bank while applying this segmentation procedure to the ridgeness map in Fig. 7. In addition to the ridgeness threshold h , which is natively expressed through its rank- based relative value h * , four cleaning parameters have to be entered: the geodesic areas a CCmin and a SKmin (m 2 ) and the geodesic lengths l BRmin and l THmin (m). Because it may be convenient to express these four values as an information- , where h 1 * is the relative upper bound defining h 1 through the same rank-based rule that defines h from h * . The current setup for the ridge and valley icons is h 1 * = 10.0 and 2.5, respectively. From a perceptual standpoint, the valley lines appear with noodle-like profiles whereas in some places ridge lines appear with sausage-like profiles. The latter are not profiles with artifacts: some ridges are natively heavily notched by ripples and these configurations may have been frozen by the anisotropic filter. quanta count, the user is allowed to provide them through a number of points. These relative instances are respectively denoted {a CCmin * , a SKmin * , l BRmin * , l THmin * } . These are then internally scaled into physical units through a = 3 4 a * λ moy 2 and l = (l * -1)λ moy , where λ moy is the mean edges length of S , i.e., ~1.8 m with respect to the filtered state displayed in Fig. 6. The processing setup of the ridge lines extraction in Fig. 9 is h * = 1.0 with the cleaning parameters set {10, 10, 10, 15} giv- en in number of points. As stated previously, as long as the reassignation K ← ⌣ K is done before re-processing, the valley lines extraction is formally the same problem. The parameter set of their extraction is h * = 0.01 with the cleaning parameter set {30, 30, 40, 15} given in number of points. As their saliency levels are very poor with respect to those of the ridges, the threshold is smaller and the cleaning step becomes slightly more aggressive, especially with respect to the clipping of the secondary branches.

Discussion

High-level semi-interactive edition of topological constraints

The network in Fig. 9 still contains bifurcation errors with respect to geoscientist expertise. To scrutinize these locations, the current configurations of the bifurcations are noted in Fig. 10. These are not true segmentation errors with respect to the generic segmentation approach devised above. To successfully address these configurations, the knowledge issued from a physical model of the dune formations and interactions still has to be included in future work. However, the estimation of the path lengths through a geodesic hypermetric was already devised to counteract much of these potential geometric biases while weighting the path length by the ridgeness. However, this cannot overcome large discrepancies between a short and very salient branch and a very long branch with poor saliency. Thus, in the current situation, an expert should be allowed to proceed to an interactive topological surgery of the network while operating intuitively at the symbolic level, i.e., not requiring some time consuming interventions at the data level. This optional facility is currently addressed through a simple trick. First, all of the tail vertices involved in the structure of the trees that contain at least one secondary branch are made clickable through the insertion of a spherical pin initially shown in its neutral state, e.g., white. While navigating through the 3D scene, the user can then successively click on the pairs of tail vertices that should become trunk-linked. A color change of the pin (e.g., to red) acknowledges each selection. The segmentation process is then reiterated from the global step described in Section 3.4.3. In the initial step, any selected tail vertex now receives the same initial large length, whereas the other tails still receive an initial zero length as before. Therefore, these virtually extended paths will not begin to progress until all of the unextended paths are fully discovered. The surged network rebuild from the red pins selection, as shown Fig. 10, is depicted in the next figure.

Global statistics vs curves intrinsic orientations

Because a smoothed and re-parameterized 3D curve is mapped on each path registered in the database, some high-level ridgevalley statistics can be built at convenience. Moreover, the setup of an ad hoc iconographic instantiation can simultaneously help obtain intuitive insights while addressing the analysis of their global characteristics. An application addressing the mean orientations of the lines in spherical coordinates is depicted in Fig. 11. To this end, first, the intrinsic curvilinear frames (t F ,n F ,b F ) are simultaneously rotated around t F so that the transformed normal vectors n F belong to the plane defined by t F and the normal n of the surface. Let (t F ,n R ,b R ) denote these new orthonormal curvilinear frames. Second, the mean unit vector b R of each line is computed and encoded in spheri- cal coordinates (θ,ϕ) where θ and ϕ denote the azimuthal and polar angles, respectively. For display purposes, ridge and valley azimuths are mirrored in the West and East hemi-circles, respectively, to limit icon overlapping. The section of the iconographic paths is now made cuboidal and aligned with (n R ,b R ) . The angular dispersion of the points confirms the perceptual distribution: The azimuths of the mega-ripples (upper right part) differ from the orientation of the sinusoidal waves (diagonal part) with the mean azimuth of the barchans (lower left part) showing a greater dispersion. Because valley lines are mostly found on lee sides of sand waves, the distribution of their polar angles corroborates the main azimuth NW→SE of the current transport dynamics of this scene. . Procedure aiming at symbolic redefinition of the bifurcations that could be determined to be physically irrelevant. Considering the network in Fig. 9, seven bifurcations are identified in (h) and zoomed in (a-g), respectively. The tails of the paths involved in the bifurcations are made clickable using white pins. The user simply clicks the pin pairs he wants to become trunklinked, which turn red, and the automatic extraction process is run again using these external constraints. Some of these bifurcation hierarchies (e.g., (f,g)) are supposedly natively acceptable, whereas (a-e) need to be redefined.

Robustness and consistency

The components of the tool pair (i.e., a filter and a detector) discussed in this study are low-level measurement tools. The usage of these tools is relevant, provided that the minimum curvature radius of the structure primitives under study remains much greater than the resolution of the raw data. If this condition was not fulfilled, an alternative methodology that performs convex optimizations through a regularization model would have to be considered, thus resulting in a very different problem involving interpolation and reconstruction; then this problem would be subject to critical robustness considerations. The data resolution can be computed as the median distance between MBES soundings of the post-processed data set. This value is approximately 1 or 2 m for the currently available datasets.

Thus in this optimal setting, the robustness of the detection results is mainly linked both to the overall robustness of the complex data acquisition pipeline of an MBES survey and to its post-processing validation steps. The structures to be extracted are intrinsic geometric primitives whose inference involves the measurement of high order derivatives; therefore, for a given TIN, the accuracy of the computed results would be better than that of any manual extraction even that by a trained operator. Therefore, in the study of complex dynamic systems such as sand bedforms, the post-processing validation step of the accuracy lacks reliable prior ground truths.

Conversely, between two consecutive surveys, such bedforms can undergo drastic transformations; therefore, the consistency of the extraction of primitives should be a qualitative criterion that is more relevant than the robustness. With regard to the processing pipeline, this consistency criterion is twofold; it involves (i) the sensitivity with respect to the parameters setting and (ii) the flexibility with respect to the interpretation of the ridge network provided by the detector. The second point, already discussed in Section 3.1 and Section 4.1, relies on the ability to manage explicit bifurcations within the ridge network, thus enabling interactive insertions of high-level knowledge from the expert.

Recently, [START_REF] Debese | Ripple filtering and ridge enhancement applied to morphodynamical tracking of sand dunes[END_REF] investigate the marine inspection of another georeferenced monitoring area through the methodology discussed in Section 2 and Section 3. The challenge in their study is mainly related to the consistency with respect to the parameters setting. Their georeferenced monitoring area is a tidal bank that is significantly different from that of the banner bank described in Fig. 3. A set of 30 unevenly time-spaced surveys, ranging from one to three per year dispatched over 15 consecutive years, were conducted using two different MBES systems with similar spatial resolutions (approximately 1 m). The processing parameters of the first survey were set according to the default setup proposed in this study, and their corresponding internal instantiations were maintained static to process the other surveys. Meaningful and promising results are available within each survey because different acquisition contexts, with respect to the states of the natural system and the measurement systems, can be addressed without changing these parameters.

Computational cost

The processing tasks discussed in this study involve a sequence of four computational kernels: (i) performing anisotropic filtering, (ii) computing the ridgeness indicator map, (iii) applying the ridge detector to this map, and (iv) performing an augmented visualization of the smoothed instance of a subset of the ridge network obtained from (iii). The technical insights related to the scientific visualization of the ridge and valley network through generalized cylinders (i.e., Computer Graphics) are not within the scope of this paper. However, the last kernel (iv) uses 1D discrete primitives smoothing and drawing; therefore, its additional cost is marginal with respect to the costs of the first three kernels. Some hardware tests have been performed using a MacPro Cylinder with a Quad-Core Xeon-E5 CPU running at 3.7 GHz and having 32 GB of RAM. The overall computation cost remains mostly linear with respect to the number of vertices used by the TIN. To process a TIN that includes 1.3 million points with the default parameter settings, this computer takes 156 seconds to apply the filter, 22 seconds to compute the ridgeness map and less than 5 seconds to detect and display the ridge network. These times are obtained for the single-thread mode of operation. The most demanding computation kernel is the anisotropic filter because the relative cost of the other kernels remains marginal. Thus, a TIN issued from an MBES survey containing 2 million points could be processed with the default parameter settings in less than 300 seconds on the currently available tabletop computers. Further, this time can be easily divided by the number of CPU cores.

Comparisons with other ridge detectors

The alternative ridge extraction strategy discussed in Section 3 provides an arborescent network in which the key structural nodes are bifurcations. Thus, the latter are fundamental elements, and their occurrences are welcome, even non-ternary ones. In contrast, in typical local extraction approaches, perfect bifurcations appear as mathematical singularities (e.g., umbilics) that cannot be solved through a limited local consideration basis without some heuristics. Below, we present some comparisons with two well-known local ridge detectors: that of Cazals and Pouget (2005a), hereafter labeled "A", and that of [START_REF] Yoshizawa | Fast and robust detection of crest lines on meshes[END_REF], hereafter "B". To address path reconstructions, the detector A includes topological heuristics, whereas the detector B recovers small gaping jumps by using another heuristic rule. These detectors output a list of polygonal lines, each provided with some path integrals: PathStrength A , PathSharpness A , PathRidgeness B , PathSphericalness B , and PathCyclideness B . They can be used to post-filter the set of lines through ad-hoc thresholds. Hereafter, the detector discussed in Section 3 is labeled "C" and is used with its default setting modified with h * = 0.5 ; this detector provides its own path integral PathRidgeness C based on the saliency indicator defined in Section 3.3. [START_REF] Van Dijk | Separating bathymetric data representing multiscale rhythmic bed forms: A geostatistical and spectral method compared[END_REF] applied the detector A to automatically extract ridge lines from sand waves. However, this operation was performed through both a smoothed DEM context and a homogeneous waves orientation hypothesis; ripples and megaripples are supposed to have been be optimally discarded through an oriented 1D lowpass filter.

To perform a balanced perceptual analysis, the three corresponding post-cleaning steps are tuned to provide the network with similar complexities. This operation currently results in keeping approximately 10%, 10%, and 60% from the A outputs, B outputs, and C outputs, respectively, using the corresponding threshold settings, PathSharpness A ≥ 1E05, PathRidgeness B ≥ 1.0 and PathRidgeness C ≥ 1.0. The curvature estimations are 3-ring-based for the detectors B and C, whereas the detector A requires a 4-ring extent to provide its full capabilities through its more sophisticated estimator (Cazals and Pouget, 2005b). The resulting networks are shown in Fig. 12a. We make the following observations:

• On acute ridge locations, the detector A fails with noticeable path deviations, e.g., Fig. 12b. This result is mainly due to ill-shaped triangles and very acute ridges that are no longer compliant with the internal working hypotheses. • Conversely, the path locations provided by the detector B are very accurate. Thus, this detector can take full advantage of the filtered TIN. However, owing to its local nature, several resilient broken paths involving gaps that are larger than one edge continue to be found at locations that are heavily notched by resilient ripples, e.g., Fig. 12c. • The paths issued from C are no longer broken and are well localized. However, these paths may exhibit undersegmentation at vanishing path extremities (Fig. 12d). This behavior is not an intrinsic drawback because it depends on the input threshold setting. Now, the paths must follow the native edges of the TIN; therefore, they are less smooth than the paths provided by A or B. Nevertheless, if this path constraint resulted in an apparent limitation, it would mean that the required sampling resolution was not matched by the MBES data. Without internal changes, the detector C can be made to function in a fully global mode using a null input threshold h = 0 while continuing to accept topological constraints de- fined using the protocol described in Section 4.1. Fig. 12e shows a 3D-enhanced display of the corresponding global network result in which the branch clipping process has been disabled for didactic purposes, i.e., l BRmin * = 2 . The presence or absence of the paths labeled γ 1!4 in Fig. 12a does not indicate failure because the three detectors can provide all of these paths through their native raw outputs. The final status of these paths depends on the detection strategy during the setup of the postfiltering step, i.e., enabling or disabling false positives. The global results in Fig. 12e could act as a ground truth model versus Fig. 12a; therefore, a close inspection of Fig. 12e could result in selecting a trade-off disabling the obtention of these paths. For the 1.3 million points of the benchmark context of Section 4.4 (where the detector C requires 27 seconds), the detector A requires 357 seconds whereas the detector B requires 90 seconds; these times include the curvature estimation step.

Thus, in spite of the drawbacks and additional time consumed, performing ridge extractions by these well-known approaches may provide useful measurements while handling geomorphological structures similar to those described in this paper. However, the previous statement is valid provided that (i) the TIN is cleaned and enhanced beforehand using the anisotropic filter described in Section 2, and (ii) an ad-hoc strategy is devised to clean the over-segmented outputs found in the bag of ridges. The new detector C was just devised to address this latter point (ii) natively providing an arborescent framework. The network structured output provided by this semi-global approach should increase the flexibility and thus enable accurate interactive morphometry on complex natural structures including lineaments, the commonly found structures that typically result in both oversegmentation and topological problems. As shown in Fig. 12e, the new detection strategy could even be made global while continuing to handle topological constraints using a multi-level generalization of the protocol described in Section 4.1.

Toward automatic extraction of geodesic domains of sand bedforms

The first approach, depicted in Fig. 13a,d, transcribes the branches stored in M BR to a new binary mask M Z , where each branch still gives rise to a unique CC. Then, starting from the contours found in M Z , it simultaneously propagates geodesic fronts -one per CC -over the surface at uniform geodesic speed. If it did not collide with another before, any front stops to advance after it has travelled the geodesic distance d max . To get an optimal accuracy, the geodesic distance is computed with respect to S , not only with respect to G , through the continuous Eikonal approach described by [START_REF] Kimmel | Fast marching methods on triangulated domains[END_REF].

The second approach transposes a well-known global segmentation technique, the markers-based constrained watershed, in the geodesic morphometry context. Let ⌣ K 0 denote a copy of ⌣ K . All of the sites of ⌣ K 0 covered by M Z are set to a lower bound of ⌣ K , e.g., -K max -1 . Let ⌣ K 1 denote a copy of ⌣ K 0 . Every location of ⌣ K 1 covered by the inverse mask M Z is set to ⌣ K max . A first step performs the geodesic reconstruction of ⌣ K 1 over ⌣ K 0 . Let ⌣ K R = rec ⌣ K 0 ( ⌣ K 1 ) denote this reconstructed map where all local minima covered by M Z are now flooded. In the second step, a watershed recovers the geodesic influence zone of each branches with respect to the geodesic function ⌣ K R . However, when apply- ing a global transform to an open support cut from a larger domain, a common difficulty may arise: the choice of the boundary conditions can have drastic influences on the result relevance. To this end, the parts of the support borders ∂Ω not covered by the partition issued from the first approach are also added to the maker set M Z . The result is displayed in Fig. 13b where only the geodes- ic influence regions linked to the branches are displayed. As expected, the salient parts of the valley lines match the borders of the influence zones, especially on the lee sides of the waves. However, in other parts, some fronts have still progressed too much without colliding with a competing one.

The result displayed in Fig. 13c is the logical intersection of the two partition types. These results are promising but they are still more qualitative than quantitative. However, they will provide a relevant initial guess for a more complex algorithm involving regularization rules -theses rules also expressing theoretical or empirical knowledge about sandy bedforms.

Conclusions and perspectives

This work demonstrates that an efficient computer-based platform can be developed to semi-interactively assist geoscientists in the routine extraction of some global and intrinsic morphological characteristics of sandy bedforms from recent MBES surveys. Moreover, this could be done without smoothing their main characteristics (i.e., the ridge lines) whose dynamics could then be tracked accurately through multiple surveys, and the algorithms could work on raw bathymetric datasets while dealing with their pervasive ripples. Because these structures are often mobile with respect to both their conformations and locations, the reduction of the inter-observer variability and the preservation of the native accuracy of the data are two fundamental key points.

Recent engineering failures with respect to the premature breakdown of the basements of offshore wind farms in the North Sea, have shown that understanding these natural structures is a critical issue. After insertion into the global information pipeline, the benefits of this tool are twofold: (i) it could help marine engineers in charge of sand resource management to estimate sand volume evolution over time in a fixed area; and (ii) marine scientists could obtain reliable shape characteristics and dynamics measurements that could help them insert relevant global constraints into their complex models to calculate advanced results with tractable computational costs.

Future work will improve both the robustness and efficiency of the current approach by incorporating a more versatile anisotropic filter, a geodesic segmentation tool to quantitatively track the influence domain of each sand wave, and the semiinteractive insertion of more general topological constraints on the ridge networks. The development of a companion tool to extract the global trend of the underlying static bedforms is also required to obtain accurate measurements of sediment volume.

Notations

(κ 1 ,κ 2 ) principal curvatures, κ 1 ≥ κ 2 .

(e 1 ,e 2 ) orthogonal principal directions of (κ 1 ,κ 2 ) . (κ H ,κ G ) mean and Gauss curvatures, (κ 1 +κ 2 ) 2 , κ 1 κ 2 . (κ D ,κ H ) difference and mean curvatures, κ D =(κ 1 -κ 2 ) 2 . (H + , H - ) quadrant labels on κ P domain following sgn κ H . x(u,v,t) space and time parameterization of a surface in R 3 . S(u,v,t) continuous 2-manifold surface. n(u,v,t) unit normal vector of S . (Ω, ∂Ω) domain and Border of S . (S 0 , S t ) initial and Current states of S(u,v,t M BR arborescent setup of M SK , no branch length < l BRmin . l THmin minimum cord length of smoothed paths. .* denote a CCmin a SKmin l BRmin l THmin in number of points. (t F ,n F ,b F ) tangent, normal, binormal along smoothed path.

Fig. 2 .

 2 Fig. 2. Global overview of the location of the Lézardrieux bank. This banner bank is located in north Brittany (France). In the middle part, the 630×630 m square region visually enhances the test area used in this study.

Fig. 3 .

 3 Fig. 3. Data characteristics of the test area. (a): Shaded representation, with an exaggeration factor of 5, of the bathymetric TIN built from a set of 180000 MBES soundings acquired in the geographical area enhanced in Fig. 2. The actual bedform, whose depth ranges from 12 m to 27 m, is composed of medium, large and megaripples. All data are projected in UTM-30N. (b): Qualitative patch shape distributions through the local characteristic pair ( λ H * , λ D * ). (c): The same distribution, now shown through the pair ( κ P , λ C *). Both plots are drawn with λ 0 set to 1 m. These 2D histograms are built on a regular grid involving 256×256 accumulation bins. A point footprint is computed as an isotropic Gaussian kernel with a 1.5-bins standard deviation and is defined on a 4-bins radius support.

Fig. 4 .

 4 Fig. 4. Repartition of the six action domains of the tuning function g H + .They are defined with respect to the thresholds κ floor and κ ceil and the angular phase defined through the ratio κ 1 κ 2 = -3 . The full diffusion still occurs in do- mains (i), (iv), and (v). g H + depends only on κ 1 in domain (ii), whereas it depends on both κ 1 and κ 2 in domain (iii). Sub-domain (iii+) slightly smooths the bump-biased ridges, whereas sub-domain (iii-) performs a slight inverse-smoothing on saddle-biased ridges. Data samples currently mapped on the edge between domains (iii+) and (iii-) (i.e., perfect ridges) remain strictly frozen if they are not too close to κ ceil .

Fig. 6 .

 6 Fig. 6. Anisotropic filtering. (a): Application to the bedforms displayed in Fig. 3. (b) and (c): 2D shape distributions of the filtered TIN. A close look at Fig. 3b,c suggests to set λ C* ! 1.6 and 0., for the floor and ceil bounds of the filter, respectively -at these bounds λ D * =|λ H * |= λ C * /sqrt(2) . Because λ 0 = 1 , the anisotropic filter is thus locked on the curvatures pair ( κ floor , κ ceil ) = (0.2, 1.0). It is iterated with its generic internal setting -i.e., δt = 0.5, δ ! t = 0.1, and ! m = 2. The filtering scale m is set to 100 iterations.

Fig. 7 .

 7 Fig. 7. Mapping of the ridgeness indicator K on a filtered instance of the Lézardrieux patch. The curvatures are estimated using the ring extent order k = 3 . The angular window strength around the phases ±π/4 is defined with n = 4 . The pseudo-colors mapping focuses on the saliency level clipped with respect to the symmetric interval [-0.03, +0.03].

Fig. 8 .

 8 Fig. 8. Processing sequence showing the main steps of the extraction and visualization of the ridge lines of a sandy bedform. (a): Reconstructed ridgeness functionmap K R . (b): Cleaned and labeled CC -mask M . (c): Skeleton of M with respect to K R -mask M SK . (d):Branch-labeled decomposition of this skeleton through a geodesic hypermetric -mask M BR with all native branches preserved (nine secondaries). (e): Same mask overprinted with the basic graphic instances of the parame- terized and smoothed branches issued from the pruned skeleton. (f): Same branches displayed through their augmented 3D graphic models.

Fig. 9 .

 9 Fig. 9. Extraction of the ridge and valley network of the filtered surface, driven by the ridgeness map in Fig. 7. Concerning the iconographic working scale, the modulation of the radius of the cylinder sections operates through the range window [0.3λ moy , 2.2λ moy ] . This radius window is linearly mapped from the ridgeness interval [h * , h 1 * ], where h 1 * is the relative upper bound defining h 1 through the same rank-based rule that defines h from h * . The current setup for the ridge and valley icons is h 1 * = 10.0 and 2.5, respectively. From a perceptual standpoint, the valley lines appear with noodle-like profiles whereas in some places ridge lines appear with sausage-like profiles. The latter are not profiles with artifacts: some ridges are natively heavily notched by ripples and these configurations may have been frozen by the anisotropic filter.

  Fig.10. Procedure aiming at symbolic redefinition of the bifurcations that could be determined to be physically irrelevant. Considering the network in Fig.9, seven bifurcations are identified in (h) and zoomed in (a-g), respectively. The tails of the paths involved in the bifurcations are made clickable using white pins. The user simply clicks the pin pairs he wants to become trunklinked, which turn red, and the automatic extraction process is run again using these external constraints. Some of these bifurcation hierarchies (e.g., (f,g)) are supposedly natively acceptable, whereas (a-e) need to be redefined.

Fig. 11 .

 11 Fig. 11. Statistical analysis, in spherical coordinates, of the orientations of the ridge-valley network. The view (b), also expressing the result of the process described in Fig. 10, gives an iconographic instantiation of the network aiming at visual expression of the underlying binormals twisted with respect to the surface. The charts (a,c) depict the distribution of the mean azimuthal and polar, respectively, angle of the twisted-binormals.

Fig. 12 .

 12 Fig. 12. Comparison of three detectors aiming at extraction of structural lines from the filtered result displayed in Fig. 6a. Overprints of the three resulting networks are given in (a) with A: Cazals' detector, B: Yoshizawa's detector, and C: semi-global detector described in Section 3. The 3D zooms (b-d) point to one example of their respective drawbacks: path deviations by A in (b), path breaks by B in (c), and path under-segmentations by C in (d). The gauge network (e) displays the raw output of the detector C when it fully operates in global mode and does not perform branches clipping.

Fig. 13 .

 13 Fig. 13. Toward automatic extraction of geodesic domains of a sandy bedform. For visual purpose, these views still include the valley network, but this information is not used by the algorithm. (a): Geodesic regions recovered through uniform geodesic speed with a maximum geodesic distance d max = 45 m. (b): Geodesic regions recovered through the constrained watershed applied to the inverse ridgeness ⌣ K . (c): Intersection between these two approaches. (d): Perspective view of (a), showing that the global trend of this area is very complex and that its depth fluctuations have magnitudes similar to those of the main sand waves.

  (κ C ,κ P ) curvedness and phase, κ H κ D . λ some curvature radius -i.e., κ -1 . λ 0 actual minimal curvature radius of a surface. λ C * log-curvedness radius ln(max(λ 0 ,λ C ) λ 0 ) .

  H * ) difference and mean curvatures radius coordinates. (λ C * ,κ P ) curvedness radius and phase coordinates.

  ) . A local infinitesimal area of S around x . H(S) normal curvature vector at x . ∂ t S time partial derivative of S . g t orthogonal speed g(u,v,t) defined from (κ 1 ,κ 2 ) on S t . (κ floor ,κ ceil ) filter curvature bounds λ 0 -1 ≥ κ ceil > κ floor >0 . set of |V | vertices v . E set of |E | edges. T set of |T | triangles. p i location p x p y p z in R 3 of the i-th vertex of V . G k (v) set of neighboring vertices of v , up to its k-ring. H (S) discrete normal curvature vector at p . (δt, δ ! t ) main and secondary time steps of the filter. (m, ! m) main and secondary iteration bounds of the filter. Ŝ anisotroprically filtered state of S . denote values defined from Ŝ -e.g., κC , κP . K ridgeness indicator map κC sgn κP | sin 2κ P | n . (Ω + , Ω -) partition of Ω following sgn K . K med ,K mad MED and MAD of K over Ω + . h ridgeness threshold on K . h * indirect setting of h through K med + h * K mad . M 0 threshold vertex mask K ≥ h -inverse M 0 . M M 0 without holes and without CC area < a CCmin . K R geodesic reconstruction of K max M over K . M SK geodesic skeleton of M over K R , no CC area < a SKmin . l ij path hypermetric 1 2 (K R (i) + K R (j)) ||p j -p i || .

	! S t ! t ! S	isotropically smoothed state of S t . smoothing strength involved by ! S t . denote values defined from ! S t -e.g., ! κ 1 , ! κ H , ! g t . discrete embedding of S .

S topological structure of S -i.e., a TIN (V,E,T ) . V G graph (V,E) of S -i.e., its 1-Skeleton.

N ˆ