Piecewise smooth reconstruction of normal vector field on digital data

We propose a novel method to regularize a normal vector field defined on a digital surface (boundary of a set of voxels). When the digital surface is a digitization of a piecewise smooth manifold, our method localizes sharp features (edges) while regularizing the input normal vector field at the same time. It relies on the optimisation of a variant of the Ambrosio‐Tortorelli functional, originally defined for denoising and contour extraction in image processing [ AT90 ]. We reformulate this functional to digital surface processing thanks to discrete calculus operators. Experiments show that the output normal field is very robust to digitization artifacts or noise, and also fairly independent of the sampling resolution. The method allows the user to choose independently the amount of smoothing and the length of the set of discontinuities. Sharp and vanishing features are correctly delineated even on extremely damaged data. Finally, our method can be used to enhance considerably the output of state‐of‐the‐art normal field estimators like Voronoi Covariance Measure [ MOG11 ] or Randomized Hough Transform [ BM12 ].


Introduction
Processing discrete 3D data to enhance their geometric quality is an important task in many shape modeling and computer graphics applications. For instance, it influences deeply object rendering, shape matching, shape compression, salient part extraction or the stability of numerical simulation performed on object boundary. We focus here on processing digital surfaces, which are boundaries of sets of voxels (i.e. points of Z 3 ). These data naturally come from the segmentation or binarization of MR, X-ray tomographic or micro-tomographic images. Most of the time, the original shape of interest has a piecewise smooth boundary whose digitization introduces arithmetic artefacts and noise around its boundaries. A typical application is illustrated in Figure 2. Starting from a volumetric tomographic images of a snow sample microstructure ( Figure 2-(a)), a binarization process extracts the ice/air interface. Thermo-mechanical properties of the snow sample at this micro- scale level are related to geometrical properties of the ice/air interface [CLB * 00]. Reconstructing smooth ice crystals facets while preserving sharp edges at crystal facets boundaries ( Figure 2-(b)) is a crucial step in snow metamorphism modeling [FB08]. Preserving the digital nature of the input surface is important here to alleviate potential approximations when converting it to a triangular mesh.
In this paper, we propose a variational approach to recover such underlying piecewise smooth object. More precisely, we regularize a raw input normal vector field on the digital surface while preserving sharp features defined as loci of high discontinuities in the vector field. Although many works have been proposed for estimating the normal vector field on point clouds or digital shapes [CP05,MOG11,BM12,CLL14], their theoretical guarantees require smoothness on the original shape. These methods are often unreliable around sharp features of the surface. Statistical and voting strategies can then be used to control the anisotropic regularization of the estimation [BM12, ZCL * 13]. Note that getting a consistent normal vector field considerably facilitates further postprocessing, notably surface reconstruction (as stated in [BTS * 14]).
In contrast with former approaches, our method reconstructs a consistent and smooth normal field, and delineates discontinuities of this field at the same time. We use a variational formulation that incorporates the fit to an input normal data, its piecewise regularization, and a penalty on the regularity and length of discontinuities. When input data is an image, this functional is known as Ambrosio-Tortorelli (AT) functional [AT90]. We adapt it to the piecewise smooth reconstruction of a normal vector field. To avoid known difficulties related to the discretization and optimization of AT functional, we reformulate it in a discrete calculus setting. The model is thus able to locate precisely discontinuities as 1-dimensional curves on the digital surface, and to smooth at the same time the input normal field in directions without discontinuity. Experiments show that sharp and vanishing features are correctly delineated even on extremely damaged data.

Related works
Estimating normals on point clouds. There has been many works published on this topic in the last twenty years. We mention here only recent methods that either present theoretical robustness guarantees or incorporate sharp features in their formulation.
By formalizing and extending former works using the Voronoi diagram for normal estimation (e.g. [ACSTD07]), Mérigot et. al. estimate the geometry of point clouds with a covariance matrix analysis of the distance function [MOG11], the so-called Voronoi Covariance Measure (VCM). The VCM matrix is proved robust to Hausdorff perturbation. Another interest of this method is that edges/sharp features in data are recognizable in eigenvalues. Cuel et. al. have generalized the VCM to make it resilient to outliers [CLMT15]. However, when data is corrupted, the smoothing parameter affects both the normal estimation and the feature detection, which is either unstable or produces wide and fuzzy features.
Boulch and Marlet follow a statistical approach to estimate normals on point clouds [BM12], called Randomized Hough Transform (RHT). The idea is to consider many random triples of input points and to bin their normal direction into a spherical histogram.
Maximal vote gives the local normal. The method offers some (statistical) guarantees and the output normal field is not smoothed around edges. However, the normal field is not perfectly smooth along smooth parts of the shape. Zhang et. al. adopt a low-rank representation to segment locally the feature points into subspace [ZCL * 13]. This approach handles better point clouds with variable sampling around features, at the price of a very high computation cost and parameter tuning. Liu et. al. improve the computation time of the preceding method, essentially by restricting the set of points where heavy computations are done [LZC * 15]. Both methods use thresholds to discriminate between feature and common points.
Normal estimation on digital surfaces. Digital data are specific point clouds whose coordinates are integers (subsets of Z 3 ) associated to a combinatorial structure. As said in the introduction, such data often come from the binarization or segmentation of 3D images (e.g. thresholding, Bayesian classification, watershed). Formally, given a compact shape X ⊂ R 3 , its digitization at gridstep h is {z ∈ Z 3 , hz ∈ X }. Subsets of Z 3 can be seen as a collection of cubes, so their topological boundary forms a quadrangulated mesh, called digital surface, whose vertices have (half-)integer coordinates. They are thus approximations of continuous 2-manifolds with a specific isothetic noise model: samples are evenly spaced, quad normals are not informative, even perfect digitizations present an Hausdorff noise of at least h, due to arithmetic approximations. Hence local approaches to normal estimation cannot be convergent toward exact normals when increasing the resolution. Digital variants of Integral Invariants (II) [CLL14] or VCM [CLT14] have convergence guarantees for digitization of smooth shapes, but the radii must be adapted to the sampling grid step. Although robust to noise, these methods do not handle well sharp edges.
Sharp feature extraction. Many works aim at identifying specific geometric zones on shapes, like dominant points, sharp edges, salient parts, gathered under the generic name "feature". Here we are interested in sharp features defined on the original continuous shape boundary as discontinuities in the normal vector field. We present below a representative set of feature detection techniques which serve as benchmarks in Section 3.
For meshes or point clouds, many approaches are based on integral quantities computed in a local neighborhood. For instance, Pauly et. al. [PKG03] and Clarenz et. al. [CRT04] use Principal Component Analysis on nearby data points. Eigenvalues analysis, sometimes as a function of the neighborhood size, provides a feature score. Mérigot et. al. [MOG11] use also a ratio of eigenvalues of the VCM to identify large angle variations. The neighborhood size is used in these methods to regularize possible perturbations on data. Unfortunately it tends also to generate large and diffused features. Park et. al. [PLL12] have proposed a Tensor-Voting strategy on local surface patches. Edges are identified with a scale-space analysis of the tensor vote. As shown in [LCL15], this technique does not provide sufficiently robust results on digital surfaces. Mellado et. al.
[MGB * 12] have introduced a fast least square spherical fitting approach to a point cloud to create a multi-scale feature score. Although qualitatively relevant, experiments show that it fails to provide a precise localization of sharp features. Some authors propose a statistical approach to extract feature lines in point clouds and we mention the very recent work of Zhang et.  al.
[ZGW * 16]. Authors identify potential features as points which vary too much from growing faces. Consistency between feature points, lines, and faces is then used to eliminate false positives, and provides a reconstruction of feature lines and not solely a score. However, regularization of point positions and normals, feature detection, and topological reconstruction are performed as a pipeline involving several parameters. Finally, features can be extracted following a spectral analysis of the shape from eigenvalues of the surface Laplacian matrix [GBAL09,SOG09,SLMR14]. In this context, features are characterized by spectral quantities which are locally stable and distinguishable from its neighborhood. However, these features or saliency maps are not related to sharp edges, but more to a behaviour with respect to vibrations.
Very few works exist for extracting sharp features on digital surfaces. We mention the recent approach of Levallois et. al. [LCL15], which detects features in scale-space, using multigrid convergence properties of curvature estimation with digital integral invariants. It will be compared to our method in section 3.
Variational models for piecewise smooth reconstruction. The Mumford-Shah (MS) model is a classical variational formulation of image restoration and segmentation [MS89]. It allows the piecewise smooth reconstruction of a possibly perturbated function, while controlling the amount of discontinuities. Notably difficult to solve directly, this formulation has inspired many works in image analysis as well as more tractable variants. Most of them transform the original formulation so that discontinuities are the boundary of regions, and the result is a partition. Since we would like to process shapes with fading sharp features, such relaxations are not suitable. Our approach follows the Ambrosio-Tortorelli (AT) relaxation of MS [AT90], which can handle dangling or fading discontinuities. However, AT is not so easy to solve numerically, and we discuss this issue with more details in Section 2.4.
We also mention also the convex relaxation of MS functional proposed by Alberti et. al. [ABDM03], by means of a calibration (some well-chosen vector field). It leads to algorithms targeting global optimum of the MS for scalar functions [PCBC09] and for vectorial functions [SCC12]. However, global optimality is not always guaranteed and these approaches are computationally very demanding (more than 600 seconds for a 128 × 128 images with 32 greylevels with a fast GPU implementation reported in [PCBC09]). They are not usable on meshes with hundred of thousands faces.
Variational models for mesh denoising. Note that other image processing techniques have proven to be useful in the context of mesh processing. Apart from the ubiquitous spectral approaches, we can mention methods that regularizes meshes while trying to preserve features. Fleishman et. al. [FDCO03] propose to use bilateral filtering for this purpose. This approach is not very well adapted to digital surfaces, which are too perturbated at a local scale for this local processing. Along the same lines, Jones et. al. [JDD03] smooth a mesh by projecting vertices onto predicted local tangent planes. Both approaches are interesting on meshes whose perturbation is small, about 1/5 of the average edge length, but would fail on our data. More recently, several authors have considered l 0 (sparsity) or l 1 minimization (total variation) schemes to denoise a mesh while preserving its features [HS13, WYL * 14, WZCF15, ZWZD15]. Compressed sensing optimization is used to detect features in [WYL * 14]. All these recent approaches denoise triangulated meshes and use vertex positions as input. They achieve impressive results especially on piecewise flat shapes, when the topology of the initial mesh is correct. Although our formulation is specific to normal vector field reconstruction on digital surfaces and does not use vertex positions, reformulating our AT model for mesh processing allowing comparisons with these literature is a challenging future works.
Other works with AT functional. Chambolle and Dal Maso [CDM99,BC00] and Bourdin and Chambolle [BC00] have proposed to use finite elements to solve AT functional for restoring 2D images, with a refinement procedure valid only for planar domains. Shah uses AT functional to determine 2D and 3D shape skeletons [Sha05]. The skeleton is determined by the discontinuity set of the gradient of the squared distance to the boundary of the shape. Hence AT functional can recover smooth skeletons of controled length, even for noisy shape. Shah emphasizes the fact that the skeleton may be thick. In a more recent work, Rosman et. al.
[RBB * 11] have proposed a motion-based segmentation process of a 3D articulated shape into rigid parts. Assuming that motion transformation should vary smoothly into rigid parts and vary suddenly into non-rigid parts, they regularize a map between shapes describing different motions with AT functional to recover a motion estimation. The numerical scheme used in these two works cannot detect thin discontinuities and the resulting reconstructions are smooth everywhere, in opposition to our formulation. Finally, Kee and Kim [KK14] have studied several convex relaxations of AT functional, by using the factorability of some non-convex problems. Although interesting, their best relaxation presents some blurring around edges on given images, as well as some false contours.
We discuss more about it in the supplementary material.

Contributions
Our method takes into account the fact that, on a piecewise smooth shape, sharp features lie exactly where normal vectors are not continuous. The AT functional defines two functions, a function v describing the set of discontinuities, a function u that is a smoothed version of the raw input normal field except at places where v indicates a discontinuity. Thanks to discrete calculus formulation, we propose a global optimization approach that minimizes AT functional on digital surfaces. We obtain a robust regularization of the normal vector field and a precise localization of sharp features as a set of edges. As a side product, our method can enhance the output of state-of-the-art anisotropic normal estimators. We describe our model in Section 2 as well as its numerical implementation. Section 3 presents an experimental evaluation of our method as well as many comparisons with existing techniques. Section 4 concludes with a discussion on the pros and cons of the method as well as possible future directions of research.
2. Variational reconstruction of a piecewise-smooth normal vector field

Ambrosio-Tortorelli approximation of MS functional
Let Ω be a 2-dimensional compact domain (for instance the plane R 2 or a surface S). Let g : Ω → R d be the function we wish to approach. The continuous formulation of the Ambrosio-Tortorelli functional can be written as: (1) where u, v ∈ SBV (Ω) (continuous functions in Ω with bounded variations). Function u is designed to be a piecewise smooth approximation of g, while function v is meant to approach discontinuities. More precisely, the first term expresses the fit to the input function g, the second term tells that u must be smooth except at places where v is zero. The third term measures the smoothness of function v while the last term forces v to stay close to 1. The second term is critical since it decides between a smooth part of the reconstruction (v ≈ 1) and a discontinuity (v ≈ 0).
Ambrosio and Tortorelli have proved that functional AT ε Γconverges toward the Mumford-Shah functional as ε tends toward 0 [AT90]. This convergence means that discontinuities are approached by stiffening progressively function v around them [Bra98]. At first, a large ε makes the function v be a diffuse approximation of discontinuities. Then, parameter ε is decreased so that the strip where v is close to 0 is made more and more narrow. Ultimately, in the continuous case, the set {v = 0} is exactly the set of discontinuities. Using dimensional analysis, ε is a distance defining the width of diffusion of discontinuities. 1/α is the area over which input data is smoothed. The smaller α is, the stronger the smoothing is. Finally, 1/λ is a distance, proportional to the length of discontinuities we allow in the piecewise smooth reconstruction. The smaller λ is, the longer discontinuities are.
At first glance, this functional perfectly fits our purpose: if we define the domain Ω to be our input surface and if we set g as the raw normal vector field, considered as a normal direction field, then minimizing AT ε with ε → 0 provides (1) the piecewise smooth approximation u of this input field and (2) the set of discontinuities, i.e. sharp features, as {v ≈ 0}. Unfortunately, it is not easy to optimize this energy in the discrete world. We discuss this problem in details in Section 2.4. Just note for now that the optimal v has almost everywhere value 1, and has value 0 only on a set of null 2-Hausdorff measure: singularities are thin. Therefore it is difficult to capture v with a sampling, and most numerical schemes fail.

Digital formulation of AT functional
Our input data is a set of voxels V (i.e. a subset of Z 3 ) and the domain Ω will be the boundary of V , seen as a collection of closed unit cubes. Then Ω is decomposed as a 2-dimensional closed cubical complex K, composed of unit cube faces, edges and vertices. Note that we add a notion of adjacency between face elements to get a digital surface [Her92]. From this primal cubical complex, a dual complex is constructed by considering dual vertices at the center of each primal square face, dual edges orthogonal to primal edges joining adjacent primal square face centers and dual faces as umbrellas centered around primal vertices with borders composed of dual edges (see Fig. 3).
Eq. (1) is discretized by considering discrete k-forms as scalars associated with k-cells of K. This approach follows the discrete calculus approach of [GP10] and shares also several characteristics of the discrete exterior calculus scheme [DHLM05], modified to handle cubical cells instead of simplexes. It induces the discretization of linear operators between k-forms as matrices.
The input field g of raw normal vectors is defined by its three dual 0-forms g 1 , g 2 , g 3 associating to each primal face the three components of the normal vector. Similarly, estimated normal u is defined as three dual 0-forms u 1 , u 2 , u 3 associating to each primal face the components of the regularized normal vector. Finally v is a primal 0-form associating the feature selection scalar to each primal vertex. This definition allows us to make ε as small as we want (typically 1 4 of the distance between two voxels), and v become identically 1 everywhere except at some vertices lying on discontinuities where it will be able to retain value 0.
Let d k (resp. dk) be the primal (resp. dual) discrete exterior derivative operator from primal (resp. dual) k-forms to primal (resp. dual) k + 1-forms. Discrete inner product between primal (resp. dual) k-forms is denoted by ·,· k (resp. ·,· k ). The wedge product between a primal 0-form γ and a dual 1-form β is a dual 1-form whose discretization can be computed as where diag (·) denotes the matrix whose diagonal elements are the components of the given vector and the operator M transports a primal 0-form to a dual 1-form, and, for each edge, simply averages the values at the extremities of this edge without taking into account orientation. In our setup, it is defined as M := 1 2 |d 0 |. With u i as dual 0-forms and v as a primal 0-form, we formulate a digital version of AT ε as follows: A possibility would be to inject metrics into the discrete calculus, by means of inner products (as in [GP10]). This is difficult in our case since the input data V may be very noisy. We therefore choose identity metrics for inner products. Denoting A the matrix form of d 0 (the incidence matrix sending vertices toward edges) and B the matrix form of d0 (the incidence matrix sending faces toward edges), we rewrite (2) in matrix/vector form: where column vectors u i and v contain the scalars defining the discrete k-forms u i and v, given a numbering of the k-cells. All the matrix/vector computations in (3) and (4) correspond to very simple and local operations on the quad mesh, as illustrated on Fig. 3.

Optimization of digital variant of AT
Energy AT d ε should have null derivatives with respect to u and v at optimum. Consequently we alternately solve the following systems until equilibrium: The three linear systems in Eq.(5) are identical and estimate {u i } for a given v. Eq.(6) is a linear system that resolves v for fixed {u i }. All involved matrices are symmetric, definite, positive, so we use a Cholesky LL factorization to solve the systems. Both systems involve sparse matrix whose number of non-null coefficients is proportional to the number of faces or vertices of the input mesh. They are almost equal for digital surfaces with bounded genus.
Ambrosio and Tortorelli have shown the Γ-convergence of AT ε toward the MS functional. This means that a sequence of minimizers for decreasing ε tends toward the MS minimizer. Since our problem is convex for fixed u or fixed v, we compute at each iteration a unique minimum. A result of convex analysis related to block coordinate descent (see [Ber99], Prop. 2.7.1) ensures that the minimizers at each iteration must converge to a stationnary point of AT d ε . Remember that a large ε convexifies the AT energy. Our optimization process thus starts with a relatively large ε (generally 2 in our examples), and after convergence for a given ε, we decrease ε (typically by dividing it by 2) and the process is restarted from this initial guess. We provide further illustration in supplementary material with a visual feedback of this progressive optimization.
The input dual 0-forms (g 1 , g 2 , g 3 ) corresponds to the input raw estimation of the normal vector field on the digital surface. In Section 3 and supplementary material, we discuss about the influence of the normal vector estimator on the AT minimization. When the minimization stops, the dual 0-forms {u i } define the regularized normal vector field and the 0-form v is a scalar value in [0, 1] which is close to zero on features. The final feature extraction evaluated in Section 3 corresponds to the collection of 1-cells such that the values of v at their vertices are both smaller than 1 2 . This naive thresholding on v values could be improved but it already leads to very stable and precise results.

Discussion about discrete formulation
There are many alternate ways of discretizing functional (1): some of them should definitely be avoided, while others are acceptable but remain behind our formulation. We sum up in this section our experience on several alternative discretizations.
General comments. The continuous formulation of Ambrosio-Tortorelli involves two functions u and v defined over some domain. When ε tends to 0, the two last terms tends toward the perimeter of the set of discontinuities. The main problem when discretizing this functional comes from the fact that v tends toward a function that cannot be approached by a standard sampling. More precisely v tends toward a function that is almost everywhere 1, meaning except on sets with null Hausdorff-2 measure.
Discrete calculus versus other numerical schemes. A discretization with standard finite differences fails for this functional. Indeed, as noted for instance by Bourdin and Chambolle [BC00], if h is the sampling step, then both ε, h and h/ε must tend toward zero to get a correct approximation of AT ε, otherwise results are poor. In practice this means that ε should be between 5-10 times greater than h, leading to "features" more than 10 voxels wide. Therefore, several people have proposed to use a finite-element approach to solve this problem [CDM99,BC00] in the case where Ω is a planar square (e.g. an image domain). The mesh should then be adapted progressively to the set of discontinuities, by refining triangles so that edges are tangent to discontinuities. This is required to get a correct approximation of the gradient of u close to discontinuities. Such approaches are complex to implement. We cannot follow this strategy for two reasons. First it would be extremely costly on our input meshes, which have often 1 million faces. Second, our input meshes may be strongly perturbated and do not constitute an accurate tiling of a smooth domain. It would be very difficult to refine such mesh along discontinuities. In supplementary material, we provide a comparison between our AT formulation and classical finite difference or finite element methods on flat domains.
Possible variants in discrete calculus. From (1), it is natural to consider u as dual 0-forms and v as a 0-form. Different formulations could have been considered, for instance with both u and v being 0-forms. However, the crucial point in the formulation is the discretization of the term v 2 |∇u| 2 since it relates the two functions u and v. We observe that this term must be seen as a 1-form living on edges. It is consistent with the fact that the antiderivative of u lives on edges. Form v must thus either be defined on 1-cells or be transported to it. We have chosen u to be a dual 0-form, since the normal field is associated to faces of digital meshes, and v is a 0form so that it lives in between normals. Form v is transported by simple averaging onto edges as a 1-form.
Γ-convergence towards a unitary normal vector field. In the vectorial case, Focardi et. al. proved that the Ambrosio-Tortorelli functional Γ-converges toward a kind of Mumford-Shah functional, where the length term is equal to one to three times the perimeter [FI14]. Optimizing a unitary tangent vector field may be difficult [KCPS13], we prefer a simpler component-wise approach. In smooth parts, each component converges independently and since the input field g is unitary, the overall vector is close to be unitary almost everywhere.

Experiments
One of the main originality of AT functional is its ability to regularize the normal vector field and to localize features at the same time. For comparisons to state of the art methods, we evaluate them independently in the following sections. Our approach is parametrized by an input raw normal vector field and α, λ and ε values. In the following experiments, we discuss the impact of these parameters on the regularized field and on the feature extraction. Extended analysis is provided in supplementary material. For some experiments, noisy digital objects are considered. The noise model adapted to digital data consists in flipping the grid point value at p with probability defined by a power law k 1+dt(p) for some user-specified k ∈ [0, 1] (dt(p) being the distance of p to the boundary voxels of the original digital shape).

Feature extraction
In Figure 4, we compare the feature extraction ability of AT functional with methods working at a single scale like VCM [MOG11] or with methods using a scale-space (Pauly et al [PKG03], Mellado et al [MGB * 12] and Levallois et al [LCL15]). Except for the latter one which is specific to digital surfaces, all methods have been originally designed for point cloud processing. The digital surface is thus seen as a point cloud defined by surfel centers. Such point cloud induced by the digitization process can be seen as an Hausdorff sampling of the underlying continuous object, and is consistent with sampling hypotheses usually considered in point cloud processing (see [MOG11] for example).
For the VCM approach [MOG11], each point p of the surface is associated with a Voronoi Covariance Measure, which depends on two parameters R and r: the offset radius R dilates the input set while the convolution radius r defines the Voronoi cells that are integrated to smooth the measure. Both parameters allow to control the impact of noise in the input data while preserving geometrical information. To extract a feature score, authors compute a ratio of the eigenvalues of the convolved VCM at each point p: r(p) := λ1 λ0+λ1+λ2 . The point is considered as a sharp edge when the ratio is greater than a threshold parameter T . Value of r(p) greater than T is thus the feature score of p. Pauly et al. [PKG03] use the eigenvalues of the covariance matrix of the point cloud at each point on the shape surface in a given neighborhood. More precisely, they consider the values of τ i (p) := λ0 λ0+λ1+λ2 for a range of radii {R i } i = 0..n (λ 0 ≤ λ 1 ≤ λ 2 being the eigenvalues of the covariance matrix). Since these eigenvalues decrease as the curvature increases, τ i will be higher on edges than on flat parts of the surface. To enhance this classification, the weight ω(p) is defined as the number of times τ i is greater than τmax on the range of radii (ω(p) := Card{τ i ≥ τmax | 0 ≤ i < n}). This quantity is then used as a feature score.
Mellado et al [MGB * 12] consider a least squares spherical fitting approach at surface points. The scale-space parameter is the neighborhood size considered in the fitting. Then, following their notations, for each scale t, they fit an algebraic hyper-sphere and get the algebraic offset distance τ between p and the 0-isosurface, the unit normal η and the signed curvature κ of the hyper-sphere. Then, they compute a geometrical variation ν(p,t) at a point p defined from fitted sphere parameters. For a range of scale t ∈ [r min , rmax], the authors define a continuous feature score function f (p) := tanh(ν(p,t))dt that differentiates regions with no geometrical variations from those with high variations.
Finally, Levallois et al [LCL15] is based on a scale-space analysis of Integral Invariant based curvature estimators [PWY * 07, CLL14]. More precisely, these estimators are defined from the volume, or covariance matrix eigenvalues, of the intersection between a ball of radius R and the object surface. For a given range of radii [r min , rmax], an analysis is performed on the estimated values in order to classify each surfel of the digital surface into three classes, Edge part, Smooth part and Flat part (or zero mean curvature part).
In Figures 4, we have considered the same set of parameters for all shapes (with or without noise). Even if some existing approaches provide robust feature selection (VCM or [LCL15]), the global minimization of AT functional allows us to have more precise and thin delineation of sharp features. Furthermore, our features are outputed as a collection of 1-cells (quad edges), compared to a feature score that needs to be post-processed to localize a onedimensional feature (see zoom on the noisy OctaFlower shape). From the set of parameters chosen for our method in these experiments, we can observe that the smooth edge vanishing on the front part of fandisk object has not been reconstructed by AT compared to [LCL15] or [PKG03]. Changing the parameters would allow us to capture such feature. In supplementary material, we provide more experiments to show the influence of the parameters (α, λ and ε) on the results.

Regularization of normal vector field
We then evaluate the quality of the regularized normal vector field. To compare normal vector fields, we use a rendering of quads associated with their normal vector and a specular material. In Figure 5, we consider several noisy versions of the Fandisk object and the estimated normal vector field using VCM, RHT, II and our approach where the input normal vector field is II with radius 4. Both VCM and RHT provide robust and moderately smooth normal vector reconstructions with respect to the noise level. However, edges are also smoothed. Since II has not be designed to be anisotropic with respect to discontinuities, large integration radii lead to smooth reconstruction but all edges disappear (see row "II, r = 8"). If the integration is small (r = 4 for instance) edges are preserved but noise becomes visible. If we feed the AT functional with this II input estimation (input functional g), we obtain a smooth reconstruction with preserved edges. In supplementary material, we provide additional results when combining RHT or VCM as input to AT functional. We observe that AT functional always improves the input normal vector field while preserving the features.

Discussions and limitations
This paper proposes a new method for the piecewise regularization of normal vector fields onto digital surfaces. Its originality lies in its global variational formulation that addresses at the same time the regularization of normals and the delineation of sharp features as the locii of vector discontinuities. Experiments show that sharp features are correctly localized even on damaged data while keeping smooth normal vectors except across features. Our method can also be used to enhance the normal estimation of several state-ofthe-art methods (see also supplementary material for a quantitative evaluation).
Limitations. As many methods developed for this purpose (regularization in presence of sharp features), our formulation requires parameters (α, λ, ε). First note that parameter ε is related to the optimization process, and we have used the same sequence 2, 1, 1 2 , 1

Input data
Mérigot et al.
Mellado et al.

[PKG03]
Levallois et al.  based analyses. Note that we have used the same α and λ in all presented experiments. The supplementary material also confirms the relative invariance of the method to large ranges of parameters.
Another limitation is that the method requires a combinatorial surface structure (digital surface), and a lot of the robustness of the method is related to this topological information. Hence the method cannot be used as is on an arbitrary cloud of points. Last, the method performs a global analysis of possible discontinuities on the surface. This analysis is thus performed at a given scale (a balance between α and λ) and not in an adaptive way. Therefore the distinction between sharp or smooth feature is related to an angle variation and it cannot directly handle well variable noise over the input normal vector field.
Finally, our discrete AT formulation is not convex, so our alternate optimization procedure, while being convergent, may stop on a local minima. Our experiments nevertheless show that the multiscale approach (decreasing sequence of ε) limits this issue. Other convex relaxations of AT model (like in [KK14]) would be interesting to cast in a discrete calculus framework. The classical convex relaxation of MS model of [ABDM03, PCBC09, SCC12] seems for now too computationnaly costly for practical geometry processing.
Future works. First it would be interesting to introduce metrics in our discrete calculus formulation in order to remove some metrication artefacts introduced by our staircase surfaces. One possibility is to inject the projected area of quads along the current estimated normal into the hodge star operator . We intend to see if some results can be enhanced with this new metric.
It is clear also that our method is rather straightforwardly extensible to triangulated meshes. DEC [DHLM05] offers exactly the tools to transpose this formulation to such meshes. This reformulation to meshes would allow us to address specific mesh processing tasks such as mesh reconstruction with optimization of the vertex position, similarly to [HS13, WYL * 14, WZCF15, ZWZD15].
A more challenging task would be to transpose this work on point clouds. Approximated Laplacian operators exist on such data and should help in designing an AT functional over point clouds.
Our piecewise smooth regularization model does not require at all that g (and u) be a normal vector field. In fact, the AT functional is suited to the piecewise reconstruction of any function. We intend to check its ability to regularize arbitrary functions, vector fields or tensors, in problems involving unknown discontinuities. Last, this model is fully compatible with user-interaction either for selecting smooth zones (forcing v to be one in these places) or for delineating manually some sharp features (forcing v to be zero).