Luca Calatroni
email: calatroni@i3s.unice.fr

Martin Huska
email: martin.huska@unibo.it

Serena Morigi
email: serena.morigi@unibo.it

Giuseppe Recupero
email: giuseppe.recupero3@unibo.it

Keywords: Variational surface restoration, Willmore energy, sparse non-convex optimization, Surface denoising, Surface inpainting, Context-aware mesh completion

Technologies for 3D data acquisition and 3D printing have enormously developed in the past few years, and, consequently, the demand for 3D virtual twins of the original scanned objects has increased. In this context, feature-aware denoising, hole filling and context-aware completion are three essential (but far from trivial) tasks. In this work, they are integrated within a geometric framework and realized through a unified variational model aiming at recovering triangulated surfaces from scanned, damaged and possibly incomplete noisy observations. The underlying non-convex optimization problem incorporates two regularisation terms: a discrete approximation of the Willmore energy forcing local sphericity and suited for the recovery of rounded features, and an approximation of the 0 pseudo-norm penalty favoring sparsity in the normal variation. The proposed numerical method solving the model is parameterization-free, avoids expensive implicit volume-based computations and based on the efficient use of the Alternating Direction Method of Multipliers. Experiments show how the proposed framework can provide a robust and elegant solution suited for accurate restorations even in the presence of severe random noise and large damaged areas.

A unified surface geometric framework for feature-aware denoising, hole filling and context-aware completion

Introduction

In spite of the remarkable progresses achieved in the fields of 3D scanning and 3D printing technologies, the available digitizing techniques often produce defective data samples corrupted by random noise and often subject to a local lack of data. Typically, this is mainly due to occlusions, surface reflection, scanner placement constraints, etc. In the context of digital restoration of cultural heritage art-works, for instance, the scanned object itself (e.g., the archaeological findings) may be incomplete and damaged due to the fact that some of its (missing) parts have been ruined over time due to wear and tear. In these cases, to facilitate the downstream processing of its digital content, the object shape needs to be denoised and repaired. Generally speaking, desirable properties of a surface repair toolkit shall include:

• Feature-aware denoising: remove undesirable noise or spurious information from the data, while preserving original features, including edges, creases and corners, with special care on the robustness for defective and incomplete point sets. • Smooth hole filling/Inpainting: inpainting is the process of recovering a missing or damaged region in the surface by filling it in a plausible way using available information. The result of a surface-inpainting operation depends on specific application considered. In digital cultural heritage restoration, for instance, surface inpainting is understood as the recovering of the holes in the data or the removal of the scratches/cracks possibly present in the scanned objects. In prototype manufacturing, the goal is shifted towards a waterproof virtual reconstruction, so that the related operation is rather interpreted as a smooth hole filling. Either case, all damaged areas should be filled in a seamless way that is minimally distinguishable from their surrounding regions. • Context-aware completion: when a priori knowledge on the missing/damaged parts of the scanned model is known, it is desirable that the completion of the damaged areas occurs by pasting known data -such as template patchesautomatically or semi-automatically under user guidance. This allows, for example, to repair a damaged part of an artifact by filling the region of interest with a patch taken from a valid/undamaged region of the model itself or even from other 3D geometric models.

An example of the three geometric tasks processed by the proposed geometric variational framework is illustrated in Fig. 1. The original noisy and incomplete scanned angel mesh (see Fig. 1(left)) is denoised while keeping all the holes, see Fig. 1(center, first row). Then the inpainting tool filled the holes smoothly, as shown in Fig. 1(center, second row) driven by the inpainting mask illustrated on the left. The large damaged region on the head is recovered by replacing a hair curl patch selected from a different, undamaged, mesh, see the recovered mesh in Fig. 1(center, third row). Finally, the completion of the damaged part together with hole filling is performed and illustrated in Fig. 1(right).

We propose a unified approach for these challenging geometric tasks by defining a variational problem encoding a priori knowledge of the particular problem (i.e. the mask operators) directly in the cost functional.

Here, it is assumed that a corrupted surface S embedded into R 3 , and possibly characterized by the presence of a damaged (incomplete) region S D ⊂ S, is represented by a triangulated mesh M 0 = (V 0 , T 0 ) with V 0 ∈ R n V ×3 being the set of n V vertices, and T 0 ∈ R n T ×3 being the set of n T triangle faces, sharing n E edges.

The aforementioned three geometry processing tasks are addressed by means of the following unified variational formulation

V * ∈ arg min V J (V ; M E ), (1)
J (V ; M E ) := λχ S\S D F(V ; V 0 ) + + R 1 (V ; M E ) + R 2 (V ; M c E ),
where χ S\S D : S → {0, 1} denotes the characteristic function of the subset S \ S D , while the binary mask operators M E ∈ {0, 1} n E and M c E = 1 n E -M E , characterize the specific surface geometry considered. As a result of the discretization on the triangulated mesh, the role of the characteristic function is played by a mask operator M V ∈ {0, 1} n V whose zero values identify the region S D . The set of vertices V * solution of the unconstrained optimization problem (1) defines a restored triangulated surface M * = (V * , T * ) which provides a solution of the three surface geometry tasks, depending on the particular setup considered.

In case of surface completion, the region identified by S D is replaced by a given template patch P with boundary b P . In this case, the only compatibility assumptions required is that the boundary of S D in M 0 , named b 0 , and b P have the same number of vertices. If this is not the case, a suitable subdivision process can be preliminarily applied. The template patch P can be identified on the object itself as well as on other objects, thus allowing a mesh editing process.

The proposed approach does not need any global or even local 2D parameterization, nor any sophisticated octree data structures to efficiently solve implicit volumetric computations [START_REF] Sharf | Contextbased surface completion[END_REF][START_REF] Park | Surface completion for shape and appearance[END_REF]. The data are explicitly treated as connected samples of a surface embedded in R 3 .

The functional J (V ; M E ) in ( 1) is characterized by the presence of the sum of two regularization terms: the sparsity-promoting term R 1 (V ; M E ) and the sphericity-inducing penalty R 2 (V ; M c E ). Furthermore, a fidelity term F(V ; V 0 ), weighted by the scalar parameter λ ≥ 0 is used to control the trade-off between fidelity to the observations and regularity in the solution V * of (1).

The regularizer R 1 favours solutions with piece-wise constant normal map and sharp discontinuities. Such regularizer can be designed in a way to penalize a measure of the "roughness" or bumpiness (curvature) of a mesh, or, equivalently, to promote sparsity on this measure. A natural bumpiness measure for a surface is the normal deviation measuring the normal variation between adjacent triangles. The ideal sparse-recovery term one would like to consider is the non-convex, noncontinuous 0 pseudo-norm, but its combinatorial nature makes the minimization of (1) an NP-hard problem. We then rather consider as regularizer R 1 (V ; M E ) a sparsity-promoting parametrized non-convex term, whose form provides an effective control on the sparsity of the normal deviation magnitudes being more accurate than the 1 norm, while mitigating the strong effect and the numerical difficulties of 0 pseudo-norm. Numerical experiments will show its efficiency in handling high levels of noise, producing good-shaped triangles, and faithfully recovering straight and smoothly curved edges.

As far as the R 2 regularization term is concerned, we choose it so as to encode a geometric energy, aimed to force local sphericity in correspondence of rounded regions. We considered here the Willmore energy, which has to be preferred over standard approaches based on mean curvature flow due to its scale invariance nature. Such energy is a quantitative measure of how much a given surface deviates from a round sphere and it is defined by the following curvature functional

E w (S) = 1 2 S (h 2 -k)dA, (2)
where dA is the area element, and h and k are the mean and Gaussian rigidities, respectively. The Willmore energy E w (S) is non-negative, and vanishes if and only if S is a sphere [START_REF] Bobenko | Discrete willmore flow[END_REF]. For compact and closed surfaces, and surfaces whose boundary is fixed up to first order, i.e. positions and normals are prescribed, finding the minima of ( 2) is equivalent to minimize the Willmore bending energy E h (S) = 1 2 S h 2 dA since the two functionals differ only by a constant (the Euler characteristic of the surface S), [START_REF] Willmore | Surfaces in conformal geometry[END_REF]. In this paper we present a discrete Willmore energy, which, in contrast to traditional approaches, follows an edge-based discrete formulation.

Compared to the earlier version of this work [START_REF] Huska | Sparsity-aided variational mesh restoration[END_REF] which focused only on surface denoising, this work further integrates the Willmore energy term in [START_REF] Sharf | Contextbased surface completion[END_REF] to promote fairness and extends the model usability to more tasks.

From an algorithmic point of view, we solve the (non-convex) problem (1) by means of an Alternating Direction Method of Multipliers (ADMM) scheme. This allows us to split the minimization problem into three more tractable sub-problems. Closed-form solutions for two of these problems can be found, while for the third, nonconvex, one different optimization solvers can be used. For this substep, we compare standard gradient descent, with heavy ball and Broyden-Fletcher-Goldfarb-Shanno (BFGS) schemes, endowed with suitable backtracking strategy applied to guarantee the convergence to stationary points of the sub-problem considered.

Numerical experiments will demonstrate the effectiveness of the proposed method for the solution of several exemplar mesh denoising, inpainting and completion problems.

The rest of the paper is organized as follows. In Section 2, we review the main mathematical approaches related to the three geometry tasks considered. In Section 3 the proposed geometric variational model is presented; details on its numerical optimization by means of the ADMMbased scheme are described in Section 4. In Section the three task can be realized by solving the optimization problem [START_REF] Sharf | Contextbased surface completion[END_REF]. Experimental results and comparisons are given in Section 6. We draw the conclusion in Section 7.

Related works

Standard numerical approaches solving the mesh denoising problem can be, essentially, divided into three classes. The first class inherits PDE-based techniques from analogous problems arising in image processing and addresses the task by using linear/nonlinear diffusion equations, see, e.g., [START_REF] Tasdizen | Geometric surface processing via normal maps[END_REF][START_REF] Lysaker | Noise removal using smoothed normals and surface fitting[END_REF], with particular care to preserve local curvature features [START_REF] Morigi | Nonlocal surface fairing[END_REF]. Recently, also thanks to their considerable impact in the image processing field, two further major approaches have started to be investigated: data-driven and optimization-based methods. Approaches belonging to the former class aim to learn the relationship between noisy geometry and the ground-truth geometry from a training dataset, see, e.g., [9]. Optimizationbased mesh denoising methods formulate the mesh restoration problem as a minimization problem where a denoised mesh best fitting to the input mesh while satisfying a prior knowledge of the ground-truth geometry and noise distribution is sought. These approaches grew their popularity more and more also thanks to rapid development of studies on sparsity-inducing penalties. Among them, penalty terms aimed at approximating the 0 pseudo-norm have been directly applied for denoising mesh vertices in [START_REF] He | Mesh denoising via l0 minimization[END_REF] and noisy point clouds in [START_REF] Sun | Denoising point sets via l0 minimization[END_REF]. However, the strong geometric bias favoured by the use of the 0 pseudo-norm can produce spurious overshoots and fold-backs, hence these methods may become computationally inefficient, even under small amounts of noise. Alternatively, an 1 penalty can be used. This is very frequent in recent sparse image/signal processing problems as well as in the mesh processing community, see, e.g., [START_REF] Avron | 1 -sparse reconstruction of sharp point set surfaces[END_REF], where 1 -sparsity was adopted to denoise point sets in a two-phase minimization strategy. As it is well-known, the 1 norm tends to underestimate high-amplitude values, thus struggling in the recovery under highlevel noise and presenting undesired staircase and shrinkage artifacts.

In this work, we focus on an optimizationbased approach for mesh denoising, and exploit a non-convex penalty approximating the 0 pseudonorm so as to induce sparsity without artifacts and promoting fairness.

As far as the hole filling or repairing problem is concerned, standard approaches range from fourth-order surface diffusion PDE methods [START_REF] Clarenz | A finite element method for surface restoration with smooth boundary conditions[END_REF], to volumetric approaches mainly based on signed distance functions to implicitly represent the surface [START_REF] Davis | Filling holes in complex surfaces using volumetric diffusion[END_REF][START_REF] Argudo | Àlvar Vinacua: Biharmonic fields and mesh completion[END_REF]. Other popular non-polygonal methods rely on Radial Basis Functions implicit interpolations [START_REF] Casciola | Fast surface reconstruction and hole filling using positive definite radial basis functions[END_REF], and Moving Least Squares [START_REF] Wang | Filling holes on locally smooth surfaces reconstructed from point clouds[END_REF]. A commonly-adopted fairness prior is ∆V 2 2 , proposed for smooth hole filling with the so-called least squares meshes, see [START_REF] Sorkine-Hornung | Leastsquares meshes[END_REF].

Our penalty function also favors smoothness. However, unlike least squares meshes, we adopt a nonlinear curvature measure -the Willmore energy -which leads to a more rounded shape filling.

When some a priori knowledge on the missing part is available, we can do more than simply fill the hole as we can complete/repair the hole geometry with a context-aware template patch, with a minimally distinguishable transition zone. Many efforts have been devoted to the automatic selection of the template patch in the object itself, under a best-matching assumption in a contextsensitive manner [START_REF] Sharf | Contextbased surface completion[END_REF][START_REF] Park | Surface completion for shape and appearance[END_REF], or by similarity between synthesizing geometry features [START_REF] Harary | Contextbased coherent surface completion[END_REF].

Our completion proposal is rather based on the assumption that the patch to be pasted has been pre-selected and placed at the desired position by the user.

Variational Recovery Model

Solving the variational problem (1) on surfaces requires the definition of the discrete manifold representing the underlying object of interest as well as the discrete approximation of the first-order differential operators involved.

We thus assume M := (V, T ) to be a triangulated surface (mesh) of arbitrary topology approximating a 2-manifold S embedded in R 3 , with

V = {v i } n V i=1 ∈ R n V ×3
being the set of vertices, and T ∈ N n T ×3 the set of face triangles,

T = {τ i } n T i=1 . Implicitly, we further denote by E ⊆ V × V ∈ N n E ×2 , E = {e j } n E
j=1 the set of edges. We denote the first disk, i.e. the triangle neighbors of a vertex v i , by Let N : R n V ×3 → R n T ×3 be the mapping computing the piecewise-constant normal field over the triangles of the mesh, where the m-th element is the outward unit normal at face

D(v i ) = {τ m | v i ∈ τ m }.
τ m = (v i , v j , v k ), defined as N m (V ) := (v j -v i ) × (v k -v i ) (v j -v i ) × (v k -v i ) 2 T ∈ R 3 , m = 1, . . . , n T . (3)
Notice that the normal vector's sign depends on the orientation of the face. The desire for consistently oriented normals is that adjacent faces have consistent orientation. Under this discrete setting, the scalar functions x, y, z : Ω ⊂ R 2 → R defined on S are sampled over the vertices v i = (x i , y i , z i ) ∈ V of the mesh M, and are understood as piecewise linear functions. We now introduce the discretization of the gradient operator on a 3D mesh. Since the normal field is piecewise-constant over the mesh triangles, the gradient operator vanishes to zero everywhere but the mesh edges along which it is constant. Therefore, the gradient operator discretization is represented by a sparse matrix D ∈ R n E ×n T defined by

D ij =    l i if τ j τ k = e i , k > j, -l i if τ j τ k = e i , k < j, 0 otherwise , (4)
where

l i = e i 2 , i = 1, . . . , n E is the length of i-th edge.
The matrix D can be decomposed as D = L D, with L = diag{l 1 , l 2 , . . . , l n E } being the diagonal matrix of edge lengths, whose values may be updated during the iteration scheme considered, and D ∈ R n E ×n T an edge-length independent sparse matrix.

Key ingredients of the proposed formulation (1) are the two operator masks M V and M E . The role of the mask M E is to adapt the recovery according to the surface morphology, while M V selects the region to be preserved in the inpainting and completion tasks. M E is a sharp detection mask represented by a binary vector M E ∈ {0, 1} n E which has 1s in correspondence with sharp edges. Recalling that the dihedral angle associated to the edge e i is the angle between normals to the adjacent triangle faces τ and τ s which share e i , we classify e i as a sharp edge if the dihedral angle θ s ∈ [0, 360) is greater than a given threshold th. In formulas

(M E ) i = 1 if (θ s > th) 0 otherwise . (5)
Given M E , its complementary mask is the vector M c E = 1 n E -M E . Fig. 2 shows M E for three different surface meshes, where we empirically set th = 30, which typically produces good results.

The influence of the choice of the mask M E in realizing the denoising task is shown in Fig. 3. The perturbed sharp sphere is illustrated in Fig. 3 on the left panel, and the denoised meshes on the right panel, obtained by applying the proposed method under the choice 3(b) and Fig. 3(e), respectively. The spacevariant mask M E obtained with th = 30, and illustrated in Fig. 3(c) is applied to obtain the denoised mesh in Fig. 3(d).

M E = 0 n E , M E = 1 n E , in Fig.
Stemming from the consideration by which a general scanned surface is characterized by sharp as well as rounded features, we specify the form of problem (1) to determine solutions V * which are close to the given data V 0 according to the observation model,

V * ∈ arg min V ∈R n V ×3 J (V ; M E , λ, a) J (V ; M E , λ, a) := λ 2 n V i=1 ((M V ) i (V i -(V 0 ) i )) 2 + + n E j=1 (M E ) j φ (DN (V )) j 2 ; a + +(M c E ) j (DN (V )) j 2 2 , (6)
where • 2 denotes the Frobenius norm. The functional in ( 6) involves three terms designed to meet three different and competing requirements that arise quite naturally from the intuitive concept of surface recovery: (1) fidelity to the known data;

(2) a parametric (defined in terms of the parameter a ∈ R + ) discontinuity-preserving smoothing favoring piece-wise constant normals; (3) smooth connection between parts and inside unknown regions. The functional is composed by the sum of smooth convex (quadratic) terms and a nonsmooth non-convex regularization term, thus the functional J in ( 6) is bounded from below by zero, non-smooth and can be convex or non-convex depending on the values of M E and a.

Sparsity-inducing penalty

We aim at constructing a parameterized sparsitypromoting regularizer characterized by a tunable degree of non-convexity a ∈ R + inducing sparsity on the vector of components (DN ) i 2 , i = 1, . . . , n E , which represent the normal variation between adjacent triangles sharing the i-th edge.

A substantial amount of recent works has studied the class of sparsity-promoting parametrized non-convex regularizers, given their provable theoretical properties and practical performances [START_REF] Lanza | Convex image denoising via non-convex regularization with parameter selection[END_REF][START_REF] Soubies | A unified view of exact continuous penalties for 2 -0 minimization[END_REF]. We consider here one of the most effective representative of this class, i.e. the Minimax Concave (MC) penalty φ(•; a) : [ 0, +∞) → R, introduced in [START_REF] Zhang | Nearly unbiased variable selection under minimax concave penalty[END_REF] and used previously in [START_REF] Huska | Convex non-convex segmentation of scalar fields over arbitrary triangulated surfaces[END_REF] applied to (DN ) i 2 in the context of mesh editing, defined by:

φ(t; a) = - a 2 t 2 + √ 2a t for t ∈ 0, 2/a , 1 for t ∈ 2/a, +∞ (7)
which, for any value of the parameter a, satisfies the following assumptions:

• φ(t; a) ∈ C 0 (R) ∩ C 2 (R \ {0}) • φ (t; a) ≥ 0 , φ (t; a) ≤ 0, ∀ t ∈ [0, ∞) \ { 2/a} • φ(0; a) = 0, inf t φ (t; a) = -a.
We denoted by φ (t; a) and φ (t; a) the first-order and second-order derivatives of φ with respect to the variable t, respectively.

The parameter a allows to tune the degree of non-convexity, such that φ(• ; a) mimics the asymptotically constant behaviour of the 0 pseudo-norm for a → ∞, while behaves as an 1 regularization term, for values a approaching to zero. For values of a in between, the MC penalty function in ( 7) is a sparsity-inducing penalty which preserves sharp features in normal variations better than 0 -pseudo-norm regularizer, and more accurate than 1 regularizer which tends to produce shrinkage effects.

This motivated us to use it in the construction of the regularizer R 1 (V ; M E ).

Edge-based discretization of the Willmore energy

Numerical approximations of the Willmore energy in digital geometry processing and geometric modeling are mainly based either on finite element discretization and numerical quadrature [START_REF] Clarenz | A finite element method for surface restoration with smooth boundary conditions[END_REF][START_REF] Hari | Computation of open willmore-type surfaces[END_REF], or on discrete differential geometry approaches ab initio. Discrete isometric bending models, derived from an axiomatic treatment of discrete Laplace operators [START_REF] Bergou | Discrete quadratic curvature energies[END_REF], the discrete conformal vertex-based energy well-defined for simplicial surfaces using circumcircles of their faces [START_REF] Bobenko | Discrete willmore flow[END_REF][START_REF] Bobenko | A conformal energy for simplicial surfaces[END_REF], and the integer linear programming approach [START_REF] Schoenemann | On a linear programming approach to the discrete willmore boundary value problem and generalizations[END_REF] all fall into the latter class.

Here, we consider an alternative edge-based discrete approximation of the Willmore energy [START_REF] Park | Surface completion for shape and appearance[END_REF] for open triangulated surfaces M represented by polygonal meshes. This energy is a sum over contributions from individual edges

E(M) = 1 2 n E j=1 e j 2 (DN ) j 2 , (8)
where (DN ) j measures how the surface "curves" near e j . To derive the continuum limit of (8) in the limit of vanishing triangle size, we assume that S is a 2-dimensional manifold of arbitrary topology embedded in R 3 and parameterized by (X, Ω) with Ω ⊂ R 2 , an open reference domain, and define

X : Ω → S; ξ → X(ξ),
the corresponding coordinate map (that is, the parametrization of S at a given point). We denote the local coordinates in Ω as (ξ 1 , ξ 2 ). For a given point x ∈ X(Ω) ⊂ S, the tangent space T x S at x is spanned by

r 1 := ∂X(x) ∂ξ1 , r 2 := ∂X(x) ∂ξ2
, the induced metric is given by g ij = r i • r j , its inverse is denoted by g ij , so that g ik g kj = δ j i , or in matrix notation [g ij ] = [g ij ] -1 , and its determinant is defined as

det(g) ≡ |g| = 1 2 ik jl g ij g kl = 1 2 (g ij g kl -g ik g jl ).
The second fundamental form II : T x S ×T x S → R is the symmetric bilinear form represented by the coefficients

L ij = -r i • ∂ j n, 1 ≤ i, j ≤ 2.
When the grid size of the triangulation M is sent to 0, the energy (8) approximates the Willmore energy as stated by the following Prop.1.

Proposition 1 Let S ⊂ R 3 be a 2-dimensional manifold, M an underlying flat triangulated approximation of S. Let M j be regular flat triangulated surfaces M j ⊂ R 3 with size(M j ) → 0 and M j → S for j → ∞. Then, the discrete energy (8) approximates the Willmore energy of S, i.e.

lim j→∞ E(M j ) = 1 2 S (h 2 -k)dS. (9)
Proof Let us first consider the integrand of (2) in the continuum, with h = κ 1 +κ 2 = tr(L i k ) being the mean curvature and k = 1 2 κ 1 κ 2 = det(L i k ) the Gaussian curvature where κ 1 , κ 2 represent the principal curvatures. The second fundamental form with components L ij relates with the linear map L k i with respect to the basis of TxS, according to the matrix equation:

[g ij ][L ij ] = [L i j ]
, and we denote L ij = k g ik L k j , 1 ≤ i, j ≤ 2. Following notations in [START_REF] Seung | Defects in flexible membranes with crystalline order[END_REF], we use the identity g ij g kl = g ik g jl + il mng mj g nk in the integrand of (2) as

h 2 -2k = (L i i ) 2 + il mnL m l L n i = (g ik L ik ) 2 + il mn(g mj L jl )(g nk L ki ) = g ij g kl L ik L jl . (10
) Substituting in [START_REF] He | Mesh denoising via l0 minimization[END_REF] the Weingarten equations

∂ i n = L k i r k and L k i , i = 1, 2, we have g ij g kl L ik L jl = L j k L k j g kj g jk = L j k r j L k j r k g kj = ∂ k n • ∂ j n g kj (11)
which is the gradient of the normal vector field. Therefore, replacing ( 10)-( 11) in (2), we get

S ∂ k n • ∂ j n g kj dS. ( 12
)
For sufficiently fine, non-degenerate tessellations M j approximating S, we consider a partition of the undeformed surface S into the disjoint union of diamond-shaped tiles, T , associated to each mesh edge e. Following Meyer et al. [START_REF] Meyer | Discrete differential-geometry operators for triangulated 2-manifolds[END_REF], one can use the barycenter of each triangle to define these regions or, alternatively, the circumcenters. Over such a diamond partition, the integral ( 12) is defined as the sum over all the diamond tiles, which reads

S ∂ k n • ∂ j n g kj dS = n E i=1 ( T )i |∂ i n| 2 d T . ( 13
)
If the triangles do not degenerate, we can approximate the area of the diamond related to the edge e i in M j by e i 2 , i.e. d T ≈ e i 2 , which implies that

|∂ i n| 2 d T ≈ (DN ) i 2 e i 2 .
The result of the limiting process depends on the triangulations considered. In particular, we assume the triangulations of S consist of almost equilateral triangles. For our purposes, the discrete Willmore energy will be used based on the observation that (2) is invariant under rigid motions and uniform scaling of the surface, which implies that E(S) itself is a conformal invariant of the surface S, see [START_REF] Bobenko | Discrete willmore flow[END_REF].

Remark 1 Even if the introduced discrete formulation is very simple when compared with the ones introduced in [START_REF] Bobenko | Discrete willmore flow[END_REF][START_REF] Schmidt | Universal formulae for the limiting elastic energy of membrane networks[END_REF], it practically produces good results. In order to validate the effective applicability of the proposed discrete Willmore energy, we evaluated E(M) in ( 8) on a uniformly tessellated sphere, for decreasing average edge-size h = {0.1208, 0.0308, 0.0076}. The achieved energy values E(M h1 ) = 0.1182, E(M h2 ) = 0.0075, E(M h3 ) = 0.00047, tend to zero, as theoretically expected from (2).

Numerical solution of the optimization problem

In this section, we illustrate the ADMM-based iterative algorithm used to compute the numerical solution of [START_REF] Tasdizen | Geometric surface processing via normal maps[END_REF].

In order to define the ADMM iteration on triangular mesh surfaces, we first consider a matrix variable N ∈ R n T ×3 with row components defined as in [START_REF] Bobenko | Discrete willmore flow[END_REF], and resort to the variable splitting technique by defining t ∈ R n E ×3 as t := DN , where D is defined in (4). The optimization problem ( 6) can be thus reformulated as

{V * , N * , t * } ∈ arg min V,N,t λ 2 n V i=1 ((M V ) i (V i -(V 0 ) i )) 2 + + n E j=1 (M E ) j φ ( t j 2 ; a) + (M c E ) j t j 2 2 , s.t. t = DN, N = N (V ) . (14)
We define the augmented Lagrangian functional associated to problem [START_REF] Davis | Filling holes in complex surfaces using volumetric diffusion[END_REF] as

L(V, N, t, ρ 1 , ρ 2 ; λ, β 1 , β 2 , a) := λ 2 n V i=1 ((M V ) i (V i -(V 0 ) i )) 2 + (15)
+

n E j=1 (M E ) j φ t j 2 ; a + (M c E ) j t j 2 2 - -ρ 1j , t j -(DN ) j + β 1 2 t j -(DN ) j 2 2 + + n T m=1 -ρ 2m , N m -N m (V ) + β 2 2 N m -N m (V ) 2 2 ,
where β 1 , β 2 > 0 are scalar penalty parameters, and

ρ 1 ∈ R n E ×3 , ρ 2 ∈ R n T ×3
represent the matrices of Lagrange multipliers associated to the constraints. We now consider the following saddle-point problem:

Find (V * , N * , t * , ρ * 1 , ρ * 2 ) ∈ R n V ×3 × × R n T ×3 × R n E ×3 × R n E ×3 × R n T ×3 s.t. L(V * , N * , t * , ρ 1 , ρ 2 ) ≤ ≤ L (V * , N * , t * , ρ * 1 , ρ * 2 ) ≤ L (V, N, t, ρ * 1 , ρ * 2 ), ∀(V, N, t, ρ 1 , ρ 2 ) ∈ R n V ×3 × × R n T ×3 × R n E ×3 × R n E ×3 × R n T ×3 . ( 16
)
An ADMM-based iterative scheme can now be applied to approximate the solution of the saddlepoint problem ( 15)-( 16). Initializing to zeros both the dual variables ρ

(0) 1 , ρ (0) 2 and setting N (0) m = N m (V (0)
) , m = 1, . . . , n T , the k-th iteration of the proposed alternating iterative scheme reads:

t (k+1) = arg min t∈R n E ×3 L(V (k) , N (k) , t; ρ (k) 1 , ρ (k) 2 ) , (17)
N (k+1) = arg min N ∈R n T ×3 , Nτ =1 L(V (k) , N, t (k+1) ; ρ (k) 1 , ρ (k)
2 ) ,

V (k+1) = arg min V ∈R n V ×3 L(V, N (k+1) , t (18)
2 ) ,

ρ (k+1) 1 = ρ (k) 1 -β 1 t (k+1) -DN (k+1) , (19)
ρ (k+1) 2 = ρ (k) 2 -β 2 N (k+1) -N V (k+1) . (20)
The updates of Lagrangian multipliers ρ 1 and ρ 2 have closed form. In the following we show in detail how to solve the three minimization sub-problems ( 17),( 18) and ( 19) for the primal variables t, N and V , respectively.

Sub-problem for t.

The minimization sub-problem for t in ( 17) can be explicitly rewritten as:

t (k+1) = arg min t∈R n E ×3 n E j=1 (M E ) j φ t j 2 ; a + +(M c E ) j t j 2 2 -ρ 1j , t j -(DN ) j + + β 1 2 t j -(DN ) j 2 2 , (22)
where we omitted the constant terms in [START_REF] Argudo | Àlvar Vinacua: Biharmonic fields and mesh completion[END_REF]. Due to the separability property of φ(•; a), problem ( 22) is equivalent to n E independent, threedimensional problems for each t j , j = 1, . . . , n E in the form t (k+1) j

= arg min tj ∈R 3 φ t j 2 ; a +

α 2 t j -r (k+1) j 2 2 , (23) where r
(k+1) j := 1 β 1 + 2(M c E ) j β 1 DN (k) j + ρ (k) 1 j and α = β1+2(M c E )j (M E )j
. where we conventionally set x 0 = 0.

Necessary and sufficient conditions for strong convexity of the cost functions in ( 23) are demonstrated in [START_REF] Huska | A convex-nonconvex variational method for the additive decomposition of functions on surfaces[END_REF]. In particular, problems [START_REF] Huska | Convex non-convex segmentation of scalar fields over arbitrary triangulated surfaces[END_REF] are strongly convex if and only if the following condition holds:

β 1 + 2(M c E ) j (M E ) j > a, ∀j = 1, . . . , n E =⇒ =⇒ β 1 = ε max j {(M E ) j a -2(M c E ) j } , for ε > 1. ( 24
)
We noticed that the sub-problem is always convex when t j has associated (M E ) j = 0, as it eliminates φ(•; a) from the sub-problem. Whenever [START_REF] Clarenz | A finite element method for surface restoration with smooth boundary conditions[END_REF] holds, the unique minimizers of (23) can be obtained in closed form as

t (k+1) j = min(max(ν -ζ/ r j 2 , 0), 1) r j , where ν = α α -a and ζ = √ 2a α -a .
We remark that the condition on β 1 in (15) only ensures the convexity conditions (24) of t-subproblem (23), but does not guarantee convergence of the overall ADMM scheme.

Sub-problem for N . The minimization subproblem [START_REF] Sorkine-Hornung | Leastsquares meshes[END_REF] for N can be reformulated as:

N (k+1) = arg min N ∈R n T ×3 , Nτ =1 β 1 2 t (k+1) -DN 2 2 + + ρ (k) 1 , DN -ρ (k) 2 , N + + β 2 2 N -N V (k) 2 2 .
The first optimality conditions lead to the following three linear systems, one for each spatial coordinate of

N ∈ R n T ×3 D T D + β 2 β 1 I N = ( 25
)
β 2 β 1 N V (k) + ρ (k) 2 β 1 + D T t (k+1) - 1 β 1 ρ (k) 1
.

Since β 1 , β 2 > 0, the linear system coefficient matrix is sparse, symmetric, positive definite and identical for all three coordinate vectors. The systems can thus be solved efficiently by applying, e.g., a unique Cholesky decomposition. At each iteration, the edge lengths diagonal matrix L in D = L D, defined in (4), needs to be updated as the vertices V move to their updated position. For large meshes, an iterative solver warm-started with the solution of the last ADMM iteration, is rather preferred. A normalization is finally applied as N represents a normal field.

The reconstructed normal map N * obtained by solving ( 14) via the proposed ADMM, satisfies the orientation consistency, as proved in [START_REF] Huska | Sparsity-aided variational mesh restoration[END_REF], thus reducing the foldovers issue. This property is not trivially satisfied by most of the two-stage mesh denoising algorithms (normal smoothing and vertex update). They present the normal orientation ambiguity problem in the vertex updating stage, which provokes ambiguous shifts of the vertex position due to direction inconsistency of the normal vectors [START_REF] Zheng | Bilateral normal filtering for mesh denoising[END_REF][START_REF] Sun | Fast and effective feature-preserving mesh denoising[END_REF]. In [START_REF] Liu | Triangulated surface denoising using high order regularization with dynamic weights[END_REF], this issue is solved by an orientation-aware vertex updating scheme.

Sub-problem for V . Omitting the constant terms in [START_REF] Argudo | Àlvar Vinacua: Biharmonic fields and mesh completion[END_REF], the sub-problem for V reads

V (k+1) = arg min V ∈R n V ×3 J V (V ) ( 26
)
J V (V ) := λ 2 n V i=1 ((M V ) i (V i -(V 0 ) i )) 2 + + n T m=1 ρ (k) 2m , N m (V ) + + β 2 2 N (k+1) m -N m (V ) 2 2
.

The functional J V (V ) is proper, smooth, nonconvex and bounded from below by zero. A minimum can be obtained applying the gradient descent (GD) algorithm with backtracking satisfying the Armijo condition or using the BFGS algorithm. In order to balance between the slow convergence properties of GD and the high computational costs required to compute the operators involved in the BFGS method, we also considered a heavy-ball type rule, following [START_REF] Zavriev | Heavy-ball method in nonconvex optimization problems[END_REF], and its extension with backtracking (covering also nonsmooth problems) given in [START_REF] Ochs | ipiano: Inertial proximal algorithm for non-convex optimization[END_REF]. In particular, the heavy-ball method is a multi-step extension of gradient descent, which, starting from

V (0) = V (k) ,
iterates over V as follows

V (j+1) = V (j) -α j ∇J(V (j) ) + δ j (V (j) -V (j-1) ), j = 1, 2, . . . ( 27
)
where α j > 0 is a step-size parameter and 0 ≤ δ j < 1. Note that for δ j = 0, ( 27) reduces to the gradient descent method. In [START_REF] Ochs | ipiano: Inertial proximal algorithm for non-convex optimization[END_REF], convergence of the scheme above to stationary points is proved in the context of non-convex cost functions as the one in [START_REF] Bergou | Discrete quadratic curvature energies[END_REF], with an extension also to non-smooth scenarios.

All the numerical optimization methods here considered rely on a easily computable formula for the gradient of the functional J V in [START_REF] Bergou | Discrete quadratic curvature energies[END_REF], which is derived in the following.

Proposition 2 Let sτ m := (v j -v i ) × (v k -v i ) 2 /2 be the area of the triangle τm = (v i , v j , v k ) with updated vertices in V , and Nm(V ) = ((v j -v i ) × (v k -v i ))/(2sτ m ). For all triangles m = 1, . . . , n T , ∇v i J V (V ) = λ(M V ) i (v i -v 0 i )+ + τm∈D(vi) 1 2sτ m ρ (k) 2m -β 2 N (k+1) m - -ρ (k) 2m -β 2 N (k+1) m , Nm(V ) Nm(V ) × (v k -v j ). ( 28
)
Proof The gradient of J V (V ) in ( 26) w.r.t vertex v i ∈ V , i = 1, . . . , n V is non-zero only over the triangles sharing v i which are contained in the first disk D(v i ). Therefore the sum in ( 26) is reduced to

∇v i J V (V ) = λ(M V ) i (v i -v 0 i )+ τm∈D(vi) τm=(vi,vj ,v k ) ∇v i      z q , (v j -v i ) × (v k -v i ) gi      , where z = ρ (k) 2m -β 2 N (k+1) m , q = (v j -v i )×(v k -v i ) 2
, and the third term in [START_REF] Bergou | Discrete quadratic curvature energies[END_REF] reduces to the scalar product g i since both N (k+1) m and Nm(V ) have unitary norm. In order to compute ∇v i (g i ), we resort on the following two properties, which hold for every constant vectors w, u ∈ R 3 and can be easily proved:

1. ∇ vi w, (v j -v i ) × (v k -v i ) = w × (v k -v j ) ; 2. ∇ vi w (v j -v i ) × (v k -v i ) 2 , u = -w, u (v j -v i ) × (v k -v i ) × (v k -v j ) (v j -v i ) × (v k -v i ) 3 2 . ( 29
)
To evaluate the product rule derivative, we apply property 1, with w = z/ (v j -v i ) × (v k -v i ) 2 for the left-side term constant, while property 2 is applied with w = z and u = (v j -v i ) × (v k -v i ) for a rightside term kept constant. Combining the results leads to the explicit formula for ∇v i g i :

∇v i g i (V ) = ρ (k) 2m -β 2 N (k+1) m × (v k -v j ) (v j -v i ) × (v k -v i ) 2 - ρ (k) 2m -β 2 N (k+1) m , (v j -v i ) × (v k -v i ) (v j -v i ) × (v k -v i ) × (v k -v j ) (v j -v i ) × (v k -v i ) 3 2 , ( 30
)
which reduces to [START_REF] Schoenemann | On a linear programming approach to the discrete willmore boundary value problem and generalizations[END_REF].

In Figure 4 (first and second rows) we report the graphs showing both the energy decay and the gradient norm decay for the three different algorithms used, i.e. GD (with and without backtracking), BFGS and heavy-ball with backtracking. The plots are related to the meshes twelve (first column) and block (second column) as representative of the entire set of meshes analyzed in the experimental section.

We remark that the use of Armijo-type backtracking rule is justified by the difficult expression [START_REF] Schoenemann | On a linear programming approach to the discrete willmore boundary value problem and generalizations[END_REF], which makes the accurate estimation of the Lipschitz constant L J V of ∇J V (V ) quite challenging. In the proposed strategy, a (typically) initially large step-size α 0 is then reduced depending on whether the following inequality is verified:

J V (V (j+1) ) ≤ J V (V (j) ) -c 1 α ∇J V (V (j) ) 2 2 (31)
with c 1 ∈ (0, 1) and where V (j) denotes the j-th update of V given by [START_REF] Bobenko | A conformal energy for simplicial surfaces[END_REF].

From the convergence plots we notice that upon a manual selection of a sufficiently small constant step-size α the convergence of plain GD without backtracking is as good as the one of the heavy-ball algorithm combined with backtracking. However, the former choice is problem-dependent, hence a backtracking strategy automatically adjusting the value of α to an appropriate size is preferred.

The graphs in Figure 4 (third row) show the robustness to the initialization V (0) , and show that both GD and heavy-ball with backtracking are consistent, regardless of the chosen initialization. However, the natural and most efficient choice for V (0) is a warm start given by the matrix V (k) obtained as a solution of the problem [START_REF] Bergou | Discrete quadratic curvature energies[END_REF] in the previous ADMM iteration.

The rigorous analysis of the convergence properties of our proposed three block ADMM scheme following, e.g., [START_REF] Wang | Global convergence of admm in nonconvex nonsmooth optimization[END_REF] is not easy to derive. However, we will provide some evidence of the numerical convergence in Section 6.

A practical use of the geometry repair framework

In the following we provide some details for the practical use of the geometric framework introduced above in view of its application to the three different tasks we are interested in.

Feature-aware mesh denoising. The goal of a surface denoising algorithm is to remove undesirable noise or spurious information on a 3D mesh, while preserving original features, including edges, creases and corners. The restored surface is a 3D mesh that represents as faithfully as possible a piecewise smooth surface, where edges appear as discontinuities in the normal field. To achieve this goal, the natural choice is to set in (6) M V = 1 n V and define M E as in [START_REF] Huska | Sparsity-aided variational mesh restoration[END_REF], so as to distinguish salient edges from smooth regions. In case of severe noise, Figure 4 First and second rows: plots of the energy J V (first row) and gradient norm (second row) decay for subproblem [START_REF] Bergou | Discrete quadratic curvature energies[END_REF] with GD (with and without backtracking), heavy-ball update with backtracking and BFGS . Third row: energy decay for GD with backtracking and heavy-ball with backtracking, for three different initializations V (0) . First column mesh twelve and second column mesh block corrupted by noise level γ = 0.2. the estimate of the mask M E may be affected by false edge detections. In such case, we suggest to recompute the edge mask M E along the ADMM iterations.

Smooth hole filling/inpainting. In contrast to techniques for image inpainting, which make use of the given spatial structure of the data (the regular grid of an image), surfaces lack a natural underlying spatial domain, which brings an additional degree of freedom in the setting of problem. At the same time, vertices' positions encode both function values and the domain of the function to be reconstructed. The initial mesh, M 0 = (V 0 , T 0 ) thus has to be set as the original (possibly noisy) incomplete mesh with trivially enclosed and labeled disconnected holes -region S D -marked as zeros in M V . On the other hand, the mask M E can still be defined as in [START_REF] Huska | Sparsity-aided variational mesh restoration[END_REF], by additionally forcing zero values on the edges in S D . The proposed geometric repair algorithm then performs simultaneously denoising, outside the holes, and smooth filling in the internal part of the holes, through the regularizer R 2 .

Context-aware completion. In some applications smooth filling of holes is not sufficient: this is the case in archaeology and in general cultural heritage applications where the main goal is the reconstruction of a digital twin of a cultural heritage object. Some parts of the original 3D model can be damaged or missing but can be completed by means of characteristic parts taken from the object under consideration or from others. Given the original incomplete mesh M 0 with a region of interest bounded by a curve b 0 and characterized by vertices V ⊂ V 0 and triangles T ⊂ T 0 , together with a template patch P = (V P , T P ), bounded by a curve b P , we build a repaired mesh M * by replacing ( V , T ) by (V P , T P ) and blending the two parts through the proposed variational model.

Note that, in case the region of interest on M 0 that has to be completed is a hole, then trivially ( V , T ) are empty sets.

We assume that the template patch P is properly aligned in the correct position and that both polygonals b 0 and b P are approximants of oriented, closed, simple curves in R 3 with the same same number of vertices. The correct positioning can thus be performed either automatically (by rigid body transformation algorithms) or through user interaction.

A narrow band around b 0 , named strip(b 0 ), containing at least 2-disk of triangle neighbors adjacent to b 0 , plays the role of S D . Hence M V is the characteristic function of M 0 \ strip(b 0 ), i.e. is zeros only on strip(b 0 ).

The operator mask M E has values one for each sharp edge in both M 0 \ strip(b 0 ) and P. According to the user desiderata, the blending can be performed in three different ways: edges in strip(b 0 ) all zeros in M E to force a smooth joint with the template; edges in strip(b 0 ) all ones, to keep a sharp connection; edges in strip(b 0 ) defined by the spatially adaptive M E in (5) to maintain geometric continuity G 0 /G 1 over the blended region.

The vertices V * of the completed surface M * are obtained by minimizing [START_REF] Tasdizen | Geometric surface processing via normal maps[END_REF], properly initialized with V (0) = (V 0 \ V ) ∪ V P , while maintaining the connectivity defined by T * = (T 0 \ T ) ∪ T P . The connectivity T * is automatically achieved as we imposed b 0 ≡ b P .

We refer the reader to Fig. 1 for visual representation of the three different tasks performed. Moreover, Section 6 offers additional insights.

Numerical Examples

We validate the proposed geometric framework both qualitatively and quantitatively on a variety of benchmark triangulated surfaces characterized by different sharpness and smoothness features and on some real datasets.

At the aim of a quantitive validation, meshes M 0 = (V 0 , T 0 ) have been synthetically corrupted. The noisy vertices in V 0 correspond to underlying noise-free vertices V GT by the following additive degradation model

V 0 = V GT + η d , (32)
where the product η d accounts for the noise perturbations. Namely, η ∈ R n V is assumed to be at each vertex independently and identically distributed as a zero-mean Gaussian random variable, i.e. η i ∼ Gauss(0, σ 2 ), i = 1, . . . , n V , with known variance σ 2 , and d ∈ R n V ×3 is a vector field of noise directions with elements d i ∈ R 3 , i = 1, . . . , n V , which can be either random directions or the normals to the vertices. The perturbations are thus characterized by a noise level γ ∈ R + defined by σ = γ l, with l representing the average edge length. Quantitative evaluation is done in terms of the following error metrics, which measure the discrepancy of the computed V * , N * w.r.t. the noise-free mesh V GT , N GT :

• Mean squared angular error (MSAE)

M SAE = E[∠(N GT , N * ) 2 ] , (33)
• L 2 vertex to vertex error (E V )

E V = V * -V GT F n V . (34)
For all the tests, the iterations of the ADMM algorithm are stopped as soon as either of the two following conditions is fulfilled:

k > 200 , V (k+1) -V (k)
2 / V (k) 2 < 10 -6 .

(35) Fig. 5 shows the energy decay curve versus the number of iterations for some of the meshes reported in this section. We observe that for all meshes considered the energy converges to a stationary value. This represents an empirical validation on the numerical convergence of the proposed ADMM-based minimization scheme. Having performed a comparative analysis between inner solvers in Section 4, we used the GD algorithm with backtracking for solving the subproblem for V , with warm start strategy allowing us to restrict to a few number (three in our experiments) of GD iterations while achieving good relative accuracy.

With respect to the comparisons showed with competing approaches for mesh repairing, we remark that most of them are based on hierarchical data structures and combined with various heuristic algorithms. On the contrary, the results presented in the following are directly derived from the solution of the proposed unified mathematical optimization problem and do not require any heuristic post-processing procedures.

All the meshes are rendered in flat shading model and visualized using ParaView software. Example 1: feature-aware denoising. To evaluate the performance of the proposed method for mesh denoising, we compared the results with other state-of-the-art variational methods for mesh denoising, namely the methods introduced in [START_REF] He | Mesh denoising via l0 minimization[END_REF][START_REF] Zheng | Bilateral normal filtering for mesh denoising[END_REF][START_REF] Sun | Fast and effective feature-preserving mesh denoising[END_REF][START_REF] Zhang | Guided mesh normal filtering[END_REF], which have been kindly provided by authors of [START_REF] Zhang | Guided mesh normal filtering[END_REF] at https://github.com/bldeng/GuidedDenoising, and a learning-based approach, presented in [9]. For each method, we show their best results achieved by tuning the corresponding set of parameters.

Fig. 6 shows the denoised meshes colored by their mean curvature scalar map, with fixed range, together with zoomed details on mesh edges. From ours [START_REF] Sun | Fast and effective feature-preserving mesh denoising[END_REF] [33]

[10] [START_REF] Zhang | Guided mesh normal filtering[END_REF] [9] a visual inspection we notice remarkable overlaps in the denoised meshes obtained from the other compared methods, and severe perturbations of the triangle shapes in the reconstructed meshes. To further demonstrate how robust our approach is w.r.t. to increasing noise perturbation, in Fig. 7 we reported qualitative and quantitative results for noise levels γ = {0.2, 0.2, 0.3, 0.4, 0.5, 0.6}from top to bottom. In the last row the mesh has been corrupted by arbitrary perturbations on the noise directions (d i ) in [START_REF] Huska | A convex-nonconvex variational method for the additive decomposition of functions on surfaces[END_REF]. Below each recovered surface, we report the quantitative evaluations according to the two error metrics (M SAE × 10 2 , E V × 10 6 ). Both quantitatively and qualitatively the results confirm the effectiveness of the proposed variational model in preserving sharp features while smoothly recovering rounded parts. Finally, we can comment on the efficiency of our algorithm which computational time is, on average, one order less than the 2 -0 denoising method [START_REF] Sun | Denoising point sets via l0 minimization[END_REF] which is the slowest, while it is comparable to the other compared methods.

To improve the estimation of mask M E for severe noise, we dynamically updated the edge mask M E every three ADMM iterations.

Example 2: hole filling/inpainting. We applied our geometric framework for the recovery of various meshes M 0 which exhibit holes or damaged parts. Fig. 2 illustrates the basic workflow for the inpainting task on angel mesh which takes as input the original eventually noisy mesh M 0 (Fig. 2, left) and the inpainting mask M V , which can be of arbitrary topology, in the figure the holes to be filled are marked as 0 in M V and blue colored. The recovery of angel mesh using smooth hole filling is illustrated in Fig. 2(second row).

Fig. 8 (first row), shows the challenging Igea mesh which presents a deep groove on the left side of the mouth and a shallower one on the right cheek. Our geometric framework was able to inpaint the shallower hole perfectly, while the deep one was filled in a satisfactory, even if not complete, way. This is justified by the different contribution of Willmore vs sparsity-inducing penalties. The latter acts more strongly with respect to the former, especially for high levels of noise. Hence, adding suitable weights to the two penalties could overcame to this disparity.

The data set minerva shown in Fig. 8 (second row, first column) presents a few holes caused by the scanner acquisition, in the head and under the nose. Moreover, a vertical strip has been intentionally added to the inpainting region S D in order to remove the groove provoked by the gluing of the two parts of the minerva's face. This dataset has been provided by ENEA, Bologna, Italy, and acquired by a VIVID laser scanner. The dataset presents inherent noise due to the optical acquisition system. The result of repairing the damaged geometry and filling surface holes is illustrated in Fig. 8(c).

In Fig. 8 (third row) the inpainting framework has been applied to repair a shard from neolitic pottery received by the CEPAM laboratory (CNRS France), obtained by fusion of more fragments. The inpainting region, shown in Fig. 8 (third row, second column) has been intentionally imposed to eliminate obvious fractures between joined fragments.

Example 3: context-aware completion. We finally applied context-aware completion as an editing tool for seamless object fusion. Completion results for the meshes lion, screwdriver, and igea are illustrated in Fig. 9-10. The templates P smoothly complete the original surfaces.

A critical aspect in context-aware completion is the continuity imposed in the joint region, which we denoted by strip(b 0 ). Conditions for geometric continuity between parametric surfaces are well assessed, while for meshes a rigorous treatment on this topic is still missing. In our framework, according to the user's desiderata, the template P can be joined to M 0 , both smoothly, by setting M E (strip(b 0 )) ≡ 0, in a sharp manner by setting M E (strip(b 0 )) ≡ 1, or in a blended fashion by simply using the M E mask of one of the two meshes (or even a combination of them). Therefore, imposing different continuity conditions for strip(b 0 ) means to define in a different way the mask M E in correspondence to the strip(b 0 ).

A typical example is shown in Fig. 11(left panel) where a synthetically created hole on the fandisk mesh M 0 is filled with a similar corner patch -template P cyan colored. In the right panel, we report details onto the completion area

Conclusions

We presented a novel geometric framework for denoising, inpainting and context-based completion for the recovery of damaged and incomplete scanned data. In contrast to volumetric approaches which use complex data structures and sophisticated procedures, we formulate the solution of the three tasks in terms of a single variational problem which is parameterization-free and normal consistent. The proposed approach is intended to repair damaged and incomplete meshes resulting from range scanning as well as for all modeling operations aiming at replacing damaged or missing parts of the surface. Future investigations will focus on the study of the theoretical convergence of the proposed numerical algorithm, which minimizes a functional that is spatially variant and characterized by a convex-non-convex structure: it varies spatially from being convex (due to the presence of the Willmore energy) to non-convex (MC penalty) according to the mask operator M E . Nevertheless, the algorithm demonstrates empirical convergence and very satisfying practical performance. The encouraging results can further be extended to the completion of missing part of objects and template patches with boundaries of different topology in order to validate the process on more realistic cases. Finally, a future direction will be to couple the proposed approach with image inpainting models favoring the completion of texture-like regions.

Figure 1

 1 Figure 1 Applications of the proposed surface geometry framework: incomplete and noisy surface input M 0 (left); denoised surface M * (center, first row), inpainting mask M V and inpainting result M * (center, second row), context-aware completion of the curly hair detail without inpainting (center, third row); context-aware completion result M * with preliminary inpainting (right). The S D region is represented in blue in the masks.

Figure 2

 2 Figure 2 Examples of M E mask for three different meshes: blue colors represent values 1, while red colors represent 0 values.

Figure 3

 3 Figure 3 Effect of the mask M E on the denoising task: original noisy mesh (a); setting M E = 0n E (b); using a spacevariant M E mask (c)-(d); setting M E = 1n E (e). The perturbed sharp sphere on the left panel has been corrupted according to (32) with γ = 0.15.

Figure 5

 5 Figure 5 Empirical convergence of ADMM algorithm for some reconstructed meshes.

Figure 6

 6 Figure 6 Examples of denoising: results of noisy-free input meshes (first row) corrupted by noise levels γ = {0.15, 0.3, 0.3, 0.2, 0.2}, from left to right.

Figure 8

 8 Figure 8 Examples of surface inpainting: (a) original damaged object; (b) inpainting mask M V ; (c) inpainted surface.

Figure 9

 9 Figure 9 Examples of context-aware completion: M 0 with template patch P, output M * for two different incomplete meshes.

Figure 10 Figure 11

 1011 Figure 10 Example of context-aware completion: M 0 with template patch P (left); mask M V (middle); output M * from two different camera points of view (right).

Research is supported in part by INDaM-GNCS research project 2020. LC and GR acknowledge the support received by the Laboratory I3S, the UCA JEDI Academy on Complex Systems and the project Arch-AI-story at UCA. LC acknowledges D. Binder, L. Cassard and V. L. Coli (CEPAM, CNRS) for providing the data use for the examples reported in Figure 8.


François Laurent
email: flaurent@pasteur.fr&jbmasson@pasteur.fr

Hippolyte Verdier

Maxime Duval

Alexander Serov

Christian L Vestergaard

Jean-Baptiste Masson

TRamWAy: Mapping physical properties of individual biomolecule random motion in large-scale single-particle tracking experiments

Motivation: Single-molecule localization microscopy allows studying the dynamics of biomolecules in cells and resolving the biophysical properties of the molecules and their environment underlying cellular function. With the continuously growing amount of data produced by individual experiments, the computational cost of quantifying these properties is increasingly becoming the bottleneck of single-molecule analysis. Mining these data requires an integrated and efficient analysis toolbox. Results: We introduce TRamWAy, a modular Python library that features: (i) a conservative tracking procedure for localization data, (ii) a range of sampling techniques for meshing the spatio-temporal support of the data, (iii) computationally efficient solvers for inverse models, with the option of plugging in user-defined functions and (iv) a collection of visualisation tools and a simple web-based interface.

I. INTRODUCTION

Rapid progress in single-molecule localization microscopy (SMLM [START_REF] Betzig | Imaging Intracellular Fluorescent Proteins at Nanometer Resolution[END_REF][START_REF] Manley | High-Density Mapping of Single-Molecule Trajectories with Photoactivated Localization Microscopy[END_REF]) has allowed the recording of moving biomolecules in cells with spatio-temporal precision of tens of nanometre-milliseconds. The biomolecules' random walks allow inferring local biophysical and biochemical properties of their cellular environments (Fazel et al. 2020, Sgouralis andPressé 2017). Solving this task is a challenging inverse problem, notably due to spatio-temporal heterogeneities of both the environments and biomolecule populations, as well as the inherent stochasticity of biomolecule motion, experimental noises and short durations of most trajectories. Bayesian inference (Beheiry et al. 2016[START_REF] Bryan | Inferring effective forces for langevin dynamics using gaussian processes[END_REF] provides a principled framework to address this inverse problem. Continuous increase in both the number of trajectories and particle densities has recently paved the way to global cellular mapping (Beheiry et al. 2016[START_REF] Briane | An overview of diffusion models for intracellular dynamics analysis[END_REF][START_REF] Hoze | Heterogeneity of AMPA receptor trafficking and molecular interactions revealed by superresolution analysis of live cell imaging[END_REF][START_REF] Salomon | Dense mapping of intracellular diffusion and drift from single-particle tracking data analysis[END_REF]. The massive increase in data throughput induces new computational challenges however. Common experiments can yield millions of recorded points and hundreds of thousands of trajectories [START_REF] Giannone | High-content super-resolution imaging of live cell by uPAINT[END_REF], and models can have hundreds of thousands of parameters [START_REF] Laurent | Mapping spatio-temporal dynamics of single biomolecules in living cells[END_REF]. This imposes significant constraints on algorithmic design and computing infrastructures. In this article, we introduce TRamWAy, The Automatic Random Walk Analyzer, a Python library for Bayesian spatio-temporal analysis of dynamic localization microscopy data. TRamWAy implements the full processing chain from raw SMLM data to biophysical landscapes ready to be visually or programmatically explored. This includes simplified and robust particle tracking, tessellation and segmentation of the spatio-temporal support of the data, and inference of the parameters of the local dynamics such as diffusivity, drift, force and effective interaction energy (Fig. 1(a)).

II. SOFTWARE IMPLEMENTATION AND METHODS

In this section, we introduce a selection of key methods provided in TRamWAy. The full list of features, tutorials and notebooks can be found at https://tramway-tour.readthedocs.io (i) Particle tracking: TRamWAy can take as input both trajectories or time-stamped localizations. For treating localization data, a non-tracking graph-based algorithm is featured to extract an optimal pairing between molecules in successive frames. Trajectories can be reconstructed by linking consecutive elementary displacements. The procedure, (iv) Post-processing: Due to the large variety of possible post-analyses, TRamWAy focuses on making the outcome of these analyses robust to various sources of noise and bias. TRamWAy stores all data, inferences and analyses in HDF5-based files to be manipulated in Python with the RWA-python library, or in text files. Inferred parameter landscapes can be rendered in many ways, in Python or using a minimalist web interface. The graphical output options include movies for displaying dynamic landscapes. Furthermore TRamWAy features a Bayesian evidence test for force/drift detection in heterogeneous diffusion landscapes [START_REF] Serov | Statistical Tests for Force Inference in Heterogeneous Environments[END_REF].

TRamWAy matches the current demand for automated pipelines and high flexibility to programmatically handle multiple and large datasets. Furthermore, observed molecules are often scattered into distinct spatial domains, which may be studied independently and in parallel as regions of interest (ROIs). While exploratory data analysis in single ROIs can be done with the InferenceMAP software (Beheiry et al. 2016) in 2D and in 3D with Genuage [START_REF] Blanc | Genuage: visualize and analyze multidimensional single-molecule point cloud data in virtual reality[END_REF], TRamWAy features sets of building blocks to automate parallel processing of large-scale data with large amounts of ROIs (Fig. 1(b)). The workflow is specified in a declarative way, i.e. in arbitrary order. TRamWAy can parallelize it on the local host and onto a remote host such as a high-performance computing (HPC) cluster. On remote hosts, TRamWAy automatically dispatches user scripts after minimal parsing for execution as job arrays in a container, which is also automatically fetched if missing. Input data are not transferred; only outputs are retrieved once the jobs are complete. User-defined functions can be added to the inference.

III. SUMMARY

The TRamWAy library allows spatio-temporal mapping of biomolecules physical properties in the cell. Examples of application to Virion assembly in HIV can be found in [START_REF] Laurent | Mapping spatio-temporal dynamics of single biomolecules in living cells[END_REF] and [START_REF] Floderer | Single Molecule Localisation Microscopy Reveals How HIV-1 Gag Proteins Sense Membrane Virus Assembly Sites in Living Host CD4 T Cells[END_REF]. The library provides tools to post-process and analyse the biophyscal maps in order to link them to biological results. TRamWAy is modular and can be rapidly repurposed and is associated to numerous notebooks to ease future users'

  FIG. 1. Local (a) and cloud (b) processing chains for SMLM data: Tracking: Graph-based sequential linking between consecutive frames with tracking if needed (see §2i). Segmentation: Multi-constrained spatio-temporal meshing of trajectories. Inference: Parameter inference is performed and physical parameters are associated to the mesh. Feature extraction: Postprocessing of the spatio-temporal maps to extract relevant features, e.g. depth of effective confining potential, heterogeneity of diffusion at the cell membrane, etc. (b) Computational scalability: The inference is the main computational bottleneck. Most analysis are focused in ROIs. TRamWAy features building blocks to help the user parallelize the analysis and run it automatically on a remote HPC host with the resulting files retrieved back to the local computer.

FUNDING

This study was funded by the French Agence Nationale de la Recherche, grants ANR-17-CE23-0016, PIA/ANR-16-CONV-0005 & ANR-19-P3IA-0001.

DATA AVAILABILITY

Manual and tutorials are available at https://tramway-tour.readthedocs.io. Beheiry, M. E. et al. (2016). A primer on the bayesian approach to high-density single-molecule trajectories analysis. Biophys

TRamWAy is a Pythonlibrary and can be installed with pip and conda. The source code is available at https://github.com/DecBayComp/TRamWAy.

developments.


Dr François Fillaux
email: fillaux@glvt-cnrs.fr

Dr Alain Cousson
email: cousson@llb.saclay.cea.fr

Where are protons and deuterons in KH p D 1-p CO 3 ? A neutron diffraction study

Keywords: Neutron diffraction, Isotope mixture, Nonlocality, Quantum superposition, Hydrogen bonding

The crystals of potassium hydrogen carbonate (KHCO 3 ) and the KDCO 3 analogue are isomorphous. They are composed of hydrogen or deuterium bonded centrosymmetric dimers (HCO - 3 ) 2 or (DCO - 3 ) 2 . The space group symmetry of KH p D 1-p CO 3 (p ≈ 0.75) determined with neutron diffraction is identical to those of KHCO 3 and KDCO 3 . This is at variance with a random distribution of H and D nuclei. These crystals are macroscopic quantum systems in which protons or/and deuterons merge into macroscopic states.

Widmung

It is a great pleasure to dedicate this paper to Professor Hans-Heinrich Limbach on the occasion of his sixty fifth birthday, and to wish him many more happy and productive years of activity in chemistry.

Introduction

The formalism of quantum mechanics extrapolated from the level of electrons and atoms to that of everyday life leads to conclusions totally alien to commonsense, such as Schrödinger's Cat in a superposition of "alivedead" states and nonlocal observables [START_REF] Laloë | Quantum mechanics, strange correlations; paradoxes and theorems[END_REF][START_REF] Leggett | Testing the limits of quantum mechanics: motivation, state of play, prospects[END_REF]. Such conflicts lead to a dichotomy of interpretation consisting in that, while at the microscopic level a quantum superposition indicates a lack of definiteness of outcome, at the macroscopic level a similar superposition can be interpreted as simply a measure of the probability of one outcome or the other, one of which is definitely realized for each measurement of the ensemble [START_REF] Leggett | Macroscopic quantum systems and the quantum theory of measurement[END_REF][START_REF] Leggett | Quantum mechanics versus macroscopic realism: Is the flux there when nobody looks?[END_REF][START_REF] Leggett | Realism and the physical world[END_REF]. For open systems, this dichotomy can be legitimated by decoherence [START_REF] Zurek | Decoherence, einselection and the quantum origin of the classical[END_REF] stipulating that an initial superposition state should lose its ability to exhibit quantum interferences via interaction with the environment. However, since the quantum theory does not predict any definite dividing line between quantal and classical regimes, macroscopic quantum behaviour is possible for systems decoupled from, or very weakly coupled to, the surroundings [START_REF] Caldeira | Quantum tunnelling in a dissipative system[END_REF]. In principle, there is no upper limit in size, complexity, and temperature, beyond which such systems should be doomed to classicality.

As a matter of fact, defect-free crystals are macroscopic quantum systems with discrete phonon states at any temperature below melting or decomposition. This is an unavoidable consequence of the translational invariance of the lattice. Yet, the dichotomy of interpretation arises for O-H• • • O hydrogen bonds, when the coexistence of two configurations, say O1 -H • • • O2 and O1 • • • H -O2, is conceived of as "disorder" [START_REF] Schuster | The hydrogen bond. Recent developments in theory and experiments[END_REF]. Protons are thought of as dimensionless particles, with definite positions and momenta, moving in a double-well coupled to an incoherent thermal bath. By contrast, vibrational spectra provide unquestionable evidences of phonon states showing that the translational invariance and the quantum nature of lattice dynamics are not destroyed by proton transfer. It is therefore compulsory to elaborate a purely quantum rationale avoiding any classical ingredient [START_REF] Fillaux | Proton transfer in the KHCO 3 and benzoic acid crystals: A quantum view[END_REF]. Such a theoretical framework has been elaborated for the crystal of potassium hydrogen carbonate (KHCO 3 ) made of hydrogen bonded centrosymmetric dimers (HCO - 3 ) 2 separated by K + entities [START_REF] Thomas | Hydrogen-bond studies. LXXXIV. An x-ray diffraction study of the structures of KHCO 3 and KDCO 3 at 298[END_REF][START_REF] Thomas | Hydrogen bond studies. XCII. Disorder in (HCO 3 ) 2-2 and (DCO 3 ) 2-2 dimers: A neutron diffraction study of KHCO 3 and KDCO 3[END_REF]. The coexistence of two configurations for protons was first regarded as statistical in nature [START_REF] Benz | Jump motion of deuterons along hydrogen bonds in KDCO 3 . A deuteron relaxation study[END_REF][START_REF] Eckold | Proton disorder and phase transition in KHCO 3[END_REF]. However, vibrational spectra (infrared, Raman, inelastic neutron scattering-INS) and neutron diffraction evidence the macroscopic quantum behaviour of protons and theoretical arguments suggest that this behaviour is intrinsic to the crystal state [START_REF] Fillaux | The Pauli principle and the vibrational dynamics of protons in solids: A new spin-related symmetry[END_REF][START_REF] Ikeda | Incoherent-elastic-neutron scattering study of the vibrational dynamics and spin-related symmetry of protons in the KHCO 3 crystal[END_REF][START_REF] Fillaux | Observation of the dynamical strucuture arising from spatially extended quantum entanglement and long-lived quantum coherence in the KHCO 3 crystal[END_REF][START_REF] Fillaux | Macroscopic quantum entanglement and "super-rigidity" of protons in the KHCO 3 crystal from 30 to 300 K[END_REF][START_REF] Fillaux | A neutron diffraction study of macroscopically entangled proton states in the high temperature phase of the KHCO 3 crystal at 340 K[END_REF].

In the present paper, we examine whether macroscopic states survive in a mixed crystal KH p D 1-p CO 3 , grown from a water solution composed of H 2 O, HDO, D 2 O molecules. Since there is no segregation of H and D atoms into separate domains or clusters, it might seem straightforward to suppose that the crystal is a mixture of HH, HD, DD dimers, with probabilities p 2 , 2p(1p), (1p) 2 , respectively [START_REF] Takeda | An indication of slowing down of hydrogen atom transfer in isotopically mixed hydrogen bonds of benzoic acid crystals[END_REF][START_REF] Xue | Isotope effects associated with tunneling and double proton transfezr in the hydrogen bonds of benzoic acid[END_REF]. A random distribution of such dimers should destroy the translational invariance of the crystal, as well as-needless to say-the space group symmetry, and this should be evidenced by neutron diffraction.

We report below neutron diffraction measurements showing that a mixed crystal (p ≈ 0.75) is isomorphous to the pure analogues KHCO 3 and KDCO 3 . This is at variance with any random distribution of H and D nuclei, so the concept of macroscopic states previously elaborated for pure crystals should hold for isotope mixtures, as well.

The organization of this paper is as follows. The crystal structure of KHCO 3 at 300 K is presented in Sec. 1 and the theoretical framework for H or D macroscopic states is presented in Sec. 2. In Sec. 3 we present evidences of macroscopic proton states obtained with neutron diffraction. Finally, some consequences of the P 2 1 /a space group symmetry determined with neutron diffraction for the mixed crystal are discussed in Sec. 4.

1 The crystal structure of KHCO 3 and KDCO 3 The crystal of KHCO 3 is monoclinic, space group P2 1 /a (C 5 2h ), with four equivalent entities per unit cell (Fig. 1). Centrosymmetric dimers (HCO - 3 ) 2 linked by moderately strong hydrogen bonds, with lengths R OO ≈ 2.58 Å, are well separated by K + ions. All dimers lie practically in (103) planes, hydrogen bonds are virtually parallel to each other, and all protons are crystallographically equivalent (indistinguishable). This crystal structure is unique to probing proton dynamics along directions x, y, z, parallel to the stretching (νOH), the in-plane bending (δOH), and the out-of-plane bending (γOH) coordinates, respectively.

From 14 K to 300 K, there is no structural phase transition. The increase of the unit cell dimensions and of the hydrogen bond length are marginal. Only the population of proton sites changes significantly. Below ≈ 150 K, all dimers are in a unique configuration, say L (see arrows in Fig. 1). At elevated temperatures, protons are progressively transferred along the hydrogen bonds to less favored sites (configuration R) at ≈ 0.6 Å from the main position. The center of symmetry is preserved and all proton sites remain indistinguishable. The population of the less favored site (or interconversion degree ρ ≈ 0.22) is determined by an asymmetric double-well [START_REF] Fillaux | Calculation of infrared and raman band profiles of strong hydrogen bonds. OH stretching band and proton dynamics in crystalline potassium hydrogen carbonate[END_REF][START_REF] Fillaux | Proton dynamics in the hydrogen bond. The inelastic neutron scattering spectrum of potassium hydrogen carbonate at 5 K[END_REF], but there are controversies as to whether this leads to statistical disorder [START_REF] Benz | Jump motion of deuterons along hydrogen bonds in KDCO 3 . A deuteron relaxation study[END_REF][START_REF] Eckold | Proton disorder and phase transition in KHCO 3[END_REF] or quantum delocalization [START_REF] Fillaux | Macroscopic quantum entanglement and "super-rigidity" of protons in the KHCO 3 crystal from 30 to 300 K[END_REF][START_REF] Fillaux | A neutron diffraction study of macroscopically entangled proton states in the high temperature phase of the KHCO 3 crystal at 340 K[END_REF][START_REF] Fillaux | Proton transfer across hydrogen bonds: From reaction path to Schrödinger's cat[END_REF][START_REF] Fillaux | Proton transfer in the KHCO 3 and benzoic acid crystals: A quantum view[END_REF].

The isomorphous P 2 1 /a crystal of KDCO 3 is made of a very similar arrangement of centrosymmetric dimers. The hydrogen bond length is slightly longer (R OO ≈ 2.61 Å) and the population of the less favored site is smaller (ρ ≈ 0.12) [START_REF] Thomas | Hydrogen bond studies. XCII. Disorder in (HCO 3 ) 2-2 and (DCO 3 ) 2-2 dimers: A neutron diffraction study of KHCO 3 and KDCO 3[END_REF]. This suggests a slightly greater potential asymmetry for KDCO 3 than for KHCO 3 [START_REF] Fillaux | Calculation of infrared and raman band profiles of strong hydrogen bonds. OH stretching band and proton dynamics in crystalline potassium hydrogen carbonate[END_REF].

Theory

In this section we present the bases of the theoretical framework leading to macroscopic states for protons or deuterons [START_REF] Fillaux | Macroscopic quantum entanglement and "super-rigidity" of protons in the KHCO 3 crystal from 30 to 300 K[END_REF][START_REF] Fillaux | A neutron diffraction study of macroscopically entangled proton states in the high temperature phase of the KHCO 3 crystal at 340 K[END_REF].

The adiabatic separation

Within the framework of the Born-Oppenheimer approximation, the vibrational Hamiltonian can be partitioned as

H v = H H/D + H at + C (H/D)at , (1)
where H H/D and H at represent the sublattices of light nuclei (H + or D + ) and heavy atoms, respectively, while C (H/D)at couples the subsystems. For OHO or ODO bonds, coupling terms between OH/OD and O• • • O degrees of freedom are rather large [START_REF] Schuster | The hydrogen bond. Recent developments in theory and experiments[END_REF][START_REF] Novak | Hydrogen bonding in solids. Correlation of spectroscopic and crystallographic data[END_REF], hence beyond the framework of the perturbation theory. Two approaches, either semiclassical or quantum, are envisaged. Semiclassical protons or deuterons are dimensionless particles moving across a potential hypersurface [START_REF] Benderskii | Tunneling splittings in vibrational spectra of non-rigid molecules: IX. Malonaldehyde and its isotopomers as a test case for fully coupled multidimensional tunneling dynamics[END_REF][START_REF] Benderskii | Tunneling splittings in vibrational spectra of non-rigid molecules. X. Reaction path Hamiltonian as zero-order approximation[END_REF][START_REF] Giese | Multidimensional quantum dynamics and infrared spectroscopy of hydrogen bonds[END_REF][START_REF] Smedarchina | Tunneling dynamics of double proton transfer in formic acid and benzoic acid dimers[END_REF][START_REF] Tautermann | The ground-state tunneling splitting of various carboxylic acid dimers[END_REF]. Complex trajectories involving heavy atom coordinates lead to mass renormalization, and to incoherent phonon-assisted tunnelling [START_REF] Benz | Jump motion of deuterons along hydrogen bonds in KDCO 3 . A deuteron relaxation study[END_REF][START_REF] Eckold | Proton disorder and phase transition in KHCO 3[END_REF][START_REF] Skinner | Proton transfer in benzoic acid crystals: A chemical spin-boson problem. Theoretical analysis of nuclear magnetic resonance, neutron scattering, and optical experiments[END_REF]. This approach is quite natural when the Born-Oppenheimer surface is known, but quantum effects can be severely underestimated.

Alternatively, if the classical concept of "trajectory", totally alien to quantum mechanics, is abandoned, adiabatic separation of the two subsystems, namely H H/D and H at , may lead to tractable models [START_REF] Schuster | The hydrogen bond. Recent developments in theory and experiments[END_REF][START_REF] Fillaux | Macroscopic quantum entanglement and "super-rigidity" of protons in the KHCO 3 crystal from 30 to 300 K[END_REF][START_REF] Giese | Multidimensional quantum dynamics and infrared spectroscopy of hydrogen bonds[END_REF][START_REF] Fillaux | Extended tunnelling states in the benzoic acid crystal: Infrared and raman spectra of the OH and OD stretching modes[END_REF][START_REF] Witkowski | Infrared spectra of the hydrogen-bonded carboxylic acids[END_REF][START_REF] Maréchal | Infrared spectra of H-bonded systems[END_REF]. Then, protons or deuterons in a definite eigen state should remain in the same state in the course of time, while heavy atoms oscillate slowly, in an adiabatic hyperpotential depending on the proton/deuteron state, through the coupling term. This separation is relevant for KHCO 3 and KDCO 3 because adiabatic potentials for different proton or deuteron states do not intersect each other. In fact, the separation is rigorously exact in the ground state, since protons or deuterons should remain in this state for ever, if there is no external perturbation. Then, protons are bare fermions, while deuterons are bosons, so quantum correlations should be different for the two crystals.

Macroscopic states

Consider a crystal composed of very large numbers N a , N b , N c (N = N a N b N c ) of unit cells labelled j, k, l, along crystal axes (a), (b), (c), respectively. The two dimers per unit cell are indexed as j, k, l and j ′ , k, l, respectively, with j = j ′ . For centrosymmetric dimers, there is no permanent dipolar interaction, so interdimer coupling terms and phonon dispersion are negligible [START_REF] Fillaux | Proton dynamics in the hydrogen bond. The inelastic neutron scattering spectrum of potassium hydrogen carbonate at 5 K[END_REF][START_REF] Ikeda | Inelastic neutron scattering study of the localized dynamics of protons in KHCO 3 single crystals[END_REF][START_REF] Kashida | Inelastic neutron scattering study of KHCO 3[END_REF]. The eigen states of the sublattice of protons can be therefore represented in a rather simple way with the basis sets of eigen states for isolated dimers A dimer (H1,H2) or (D1,D2) is modelled with coupled centrosymmetric collinear oscillators in three dimensions, along coordinates {α 1jkl } and {α 2jkl } (α = x, y, z). The center of symmetry is at {α 0jkl }. The mass-conserving symmetry coordinates independent of j, k, l, and their conjugated momenta,

α s = 1 √ 2 (α 1 -α 2 + 2α 0 ), P sα = 1 √ 2 (P 1α -P 2α ), α a = 1 √ 2 (α 1 + α 2 ), P aα = 1 √ 2 (P 1α + P 2α ), (2)
lead to uncoupled oscillators at frequencies ω sα and ω aα , respectively, each with m = 1 amu. The difference ( ω sαω aα ) depends on the coupling term (say λ α ). The wave functions, {Ψ a njkl (α a )}, {Ψ s n ′ jkl (α s -√ 2α 0 )}, cannot be factored into wave functions for individual particles, so there is no local information available for these entangled oscillators. Then, the wave functions of protons (fermions) or deuterons (bosons) are subject to the symmetrization postulate of quantum mechanics [START_REF] Cohen-Tannoudji | Mécanique Quantique[END_REF].

Protons

The degenerate ground state of indistinguishable fermions must be antisymmetrized with respect to site permutation. For this purpose, the wave function is written as a linear combination

Θ 0jkl± = 1 √ 2 α Ψ a 0jkl (α a ) Ψ s 0jkl (α s - √ 2α 0 ) ± Ψ s 0jkl (α s + √ 2α 0 ) (3)
and the antisymmetrized state vectors with singlet-like |S or triplet-like |T spin symmetry are:

|0jkl+ ⊗ |S = Θ 0jkl+ ⊗ 1 √ 2 [| ↑ 1 ↓ 2 -| ↓ 1 ↑ 2 ] ; |0jkl-⊗ |T = |Θ 0jkl-⊗ 1 √ 3 | ↑ 1 ↑ 2 + | ↓ 1 ↓ 2 + 1 √ 2 [| ↑ 1 ↓ 2 + | ↓ 1 ↑ 2 ] . (4)
The oscillators are now entangled in position, momentum, and spin. In contrast to magnetic systems [START_REF] Cowley | Quantum entanglement and neutron scattering experiments[END_REF], there is no level splitting, so the symmetry related entanglement is energy-free. It is also independent of λ α . Furthermore, as there is no significant exchange integral for protons separated by ≈ 2.2 Å [START_REF] Fillaux | Comment on "quantum correlations between protons in potassium bicarbonate[END_REF], protons are not itinerant particles and there is no sizeable energy band width. Consider now the whole sublattice of protons. The spatial periodicity leads to collective dynamics and nonlocal observables in three dimensions. With the vibrational wave function for the unit cell j, k, l, namely Ξ 0jklτ = Θ 0jklτ ± Θ 0j ′ klτ , where τ = "+" or "-" for singlet-like or triplet-like symmetry, respectively, trial phonon waves are written as

Ξ 0τ (k) = 1 √ N Nc l=1 N b k=1 Na j=1 Ξ 0jklτ exp(ik • L), (5)
where k is the wave vector, L = ja + kb + lc is the lattice vector, a, b, c, are the unit cell vectors. This equation could represent collective dynamics of (H,H) dimers if they were composed bosons in a crystal made of indistinguishable dimer entities (KHCO 3 ) 2 . However, x-ray and neutron diffraction show that the structure is composed of monomer entities (KHCO 3 ) related to each other through symmetry operations. The probability density of each atom is equally distributed over all equivalent sites and, conversely, the probability density at each site includes contributions from all indistinguishable nuclei of the same kind. Consequently, antisymmetrization applied to the sublattice of indistinguishable fermions leads to k • L ≡ 0 modulo 2π.

This means that there is no phonon (no elastic distortion) in the ground state. This symmetry-related "super-rigidity" [START_REF] Fillaux | Macroscopic quantum entanglement and "super-rigidity" of protons in the KHCO 3 crystal from 30 to 300 K[END_REF] is totaly independent of any proton-proton interaction. The state vectors in three dimensions can be then written as:

|H + = |Ξ 0+ (k = 0) ⊗ |S ; |H -= |Ξ 0-(k = 0) ⊗ |T . (7)
The wave functions Ξ 0τ (k = 0) represent collective oscillations of the super-rigid lattice as a whole, with respect to the center of mass of the crystal. Finally, the ground state of the sublattice is a superposition state as:

√ N |Ξ 0+ (k = 0) ⊗ |S ; √ N |Ξ 0-(k = 0) ⊗ |T . (8)
This ground state is intrinsically steady against decoherence. Irradiation by plane waves (photons or neutrons) may single out some excited states. Entanglement in position and momentum is preserved, while the spinsymmetry and super-rigidity are destroyed. However, the spin-symmetry reappears automatically after decay to the ground state, presumably on the time-scale of proton dynamics. Consequently, disentanglement reaches a steady regime such that the amount of transitory disentangled states is determined by the ratio of density-of-states for the surroundings (atmosphere, external radiations...) and for the crystal, respectively. This ratio is so small that disentangled states are too few to be observed. Nevertheless, they allow the super-rigid sublattice to be at thermal equilibrium with the surroundings, despite the lack of internal dynamics.

The main source of disentanglement is the thermal population of excited proton states. However, even at room temperature, the thermal population of the first excited state (< 1% for γOH ≈ 1000 cm -1 ) is of little impact to measurements.

Deuterons

For the sublattice of bosons in the isomorphic crystal of KDCO 3 , Eqs (3), ( 4) and ( 6) are not relevant. There is neither spin-symmetry nor superrigidity. Dynamics are represented with symmetry coordinates (2) leading to wave functions

Θ 0jkl = α Ψ a 0jkl (α a )Ψ s 0jkl (α s - √ 2α 0 ), (9)
and phonon states

|D k± = 1 √ N Nc l=1 N b k=1 Na j=1 |Ξ 0jkl± exp(ik • L), (10)
where Ξ 0jkl± = Θ 0jkl ± Θ 0j ′ kl . Needless to say, the H and D nuclei have the same number of degrees of freedom (12 per unit cell), but the symmetrization postulate shrinks the size of the allowed Hilbert space from ∼ 12 N for bosons to ∼ 12N for fermions.

Probing quantum entanglement with neutrons

Neutrons (spin 1/2) are unique to probing the spin-symmetry of macroscopic proton states [START_REF] Caldeira | Quantum tunnelling in a dissipative system[END_REF]. However, quantum entanglement is extremely fragile, as it is not stabilized by any energy. Only "noninvasive" experiments, free of measurement-induced decoherence, are appropriate [START_REF] Leggett | Quantum mechanics versus macroscopic realism: Is the flux there when nobody looks?[END_REF]. For neutron scattering, this means (i) no energy transfer (ii) no spin-flip and (iii) particular values of the momentum transfer vector Q.

(Q = k i -k f ,
where k i and k f are the initial and final wave vectors, respectively.)

Figure 2: Diffraction patterns of KHCO 3 at 30 K (left) and 300 K (right) in the (a * , c * ) reciprocal plane at k = 0, after Ref. [START_REF] Fillaux | Macroscopic quantum entanglement and "super-rigidity" of protons in the KHCO 3 crystal from 30 to 300 K[END_REF]. The arrows emphasize ridges of intensity parallel to Q z and perpendicular to the dimer plane (dash lines along Q x ).

Noninvasive diffraction events preserving the super-rigidity occur when components of the momentum transfer along proton coordinates (Q x , Q y , Q z ) match a node of the reciprocal sublattice of protons, so the momentum transfer is the same at all sites. Then, neutrons probe super-rigid states without any induced distortion. The only information conveyed by such events is the perfect periodicity of the sublattice, so the Debye-Waller factor is equal to unity at any temperature. In addition, thanks to the spin-symmetry, the scattered intensity is proportional to the total crosssection σ H ≈ 82.0 b [START_REF] Fillaux | Macroscopic quantum entanglement and "super-rigidity" of protons in the KHCO 3 crystal from 30 to 300 K[END_REF][START_REF] Lovesey | Nuclear scattering, Theory of Neutron Scattered from Condensed Matter[END_REF]. Otherwise, if the matching condition is not realized, (i) the spin-symmetry is destroyed, (ii) the intensity scattered by protons is proportional to the coherent cross-section σ Hc ≈ 1.76 b, (iii) Bragg-peaks are depressed by the usual Debye-Waller factor for nonrigid lattices. Quantum correlations can be therefore clearly observed, thanks to the dramatic enhancement factor σ H /σ Hc ≈ 45. Furthermore, the intensity scattered by heavy atoms, proportional to σ cKCO3 ≈ 27.7 b, is depressed by the Debye-Waller factor: exp -2W KCO3 (Q). Therefore, the contribution at large Q-values is much weaker than that of the entangled sublattice and the contrast of intensity is enhanced at elevated temperatures.

The dashed lines in Fig. 1 show that proton sites are aligned along x and y, but not along z. Consequently, the noninvasive condition can be realized for Q x and Q y , whereas Q z does not coincide, strictly speaking, with a nod of the reciprocal lattice of protons. The diffraction pattern is composed of rods of diffuse scattering parallel to Q z , cigar-like shaped by the Debye-Waller factor. A detailed analysis of the scattering function reveals a complex pattern of rods of intensity, in full agreement with observations [START_REF] Fillaux | Macroscopic quantum entanglement and "super-rigidity" of protons in the KHCO 3 crystal from 30 to 300 K[END_REF][START_REF] Fillaux | A neutron diffraction study of macroscopically entangled proton states in the high temperature phase of the KHCO 3 crystal at 340 K[END_REF]. For example, Fig. 2 shows slices of the diffraction patterns of KHCO 3 at 30 and 300 K, in the (a * , c * ) reciprocal plane, at Q y = 0 (k = Q y /b * = 0). The cigar-like shaped rods are observed at Q x = 0 and ±(10.00 ± 0.25) Å-1 , in accordance with the distance of ≈ 0.6 Å between double lines of proton. These rods are evidences of macroscopic quantum correlations in the sublattice of protons. As anticipated from [START_REF] Thomas | Hydrogen-bond studies. LXXXIV. An x-ray diffraction study of the structures of KHCO 3 and KDCO 3 at 298[END_REF], these rods are not observed for KDCO 3 [START_REF] Fillaux | Observation of the dynamical strucuture arising from spatially extended quantum entanglement and long-lived quantum coherence in the KHCO 3 crystal[END_REF]. At low temperatures, the rods are partially hidden by the anisotropic continuum, centered at Q = 0, due to incoherent scattering. At 300 K, this continuum is dramatically depressed by the Debye-Waller factor, while diffraction by the super-rigid sublattice is largely unaffected. Quite remarkably, quantum correlations are not destroyed by interconversion. The sublattice is therefore a superposition of macroscopic states, |L and |R , corresponding to the structures sketched in Fig. 3.

Theory and experiments show that protons (deuterons) in the crystal field are not individual particles possessing properties in their own rights. Instead of that, they merge into macroscopic states. The theory suggests that such matter fields, intrinsic to translational invariance, should occur in many crystals. However, the diffraction patterns in Fig. 2 are easy to interpret because protons are indistinguishable and hydrogen bonds are aligned (Fig. 1). For more complex structures, seeking rods of diffuse scattering in diffraction patterns could be problematic. In the next section, we propose a complementary approach to evidence the nonlocal nature of protons and deuterons.

Isotope mixtures

In the present section we examine whether macroscopic states exist in isotope mixtures KH p D 1-p CO 3 . We anticipate that any superposition of ( 7) and ( 10) should preserve the P 2 1 /a space group symmetry of the crystal, whereas any statistical distribution of isotopes, here or there, should destroy the translational invariance. Because neutrons are scattered by nuclei, the diffraction pattern should discriminate between these two possibilities. By contrast, the P 2 1 /a symmetry evidenced by x-ray would not exclude, per se, a random distribution of isotopes, for they have identical electronic structures.

A nearly cubic specimen (3 × 3 × 3 mm 3 ) was cut from a large crystal grown by cooling slowly a saturated solution in a mixture of H 2 O (≈ 75%) and D 2 O (≈ 25%). Measurements were conducted with the four-circle diffractometer 5C2 at the Orphée reactor (Laboratoire Léon-Brillouin) [40]. The crystal wrapped in aluminum was maintained at (300 ± 1) K. Data reduction was carried out with CRYSTALS [START_REF] Watkin | CRYSTALS. issue 11[END_REF]42]. A complete determination of the positional and thermal parameters is currently in progress. The P 2 1 /a symmetry is at variance with a distribution of H and D nuclei at random. It is consistent with the existence of four sublattices, H L , H R , D L , D R , analogous to those sketched in Fig. 3, such as:

HD(p, ρ H , ρ D , ) = p[(1 -ρ H )H L + ρ H H R ] + (1 -p)[(1 -ρ D )D L + ρ D D R ]. (11)
With a diffractometer based at a reactor source one measures a limited range of the reciprocal space, typically |Q| 5 Å-1 . The matching condition for noninvasive measurements is realized only at Q x = Q y = 0, so for almost all measured Bragg-peaks, the spin symmetry is destroyed, the cross-section is σ Hc , and the sublattice of protons behaves exactly the same as a non-rigid lattice of bosons [START_REF] Thomas | Hydrogen-bond studies. LXXXIV. An x-ray diffraction study of the structures of KHCO 3 and KDCO 3 at 298[END_REF]. The four sublattices are therefore consistent with a superposition of separable macroscopic states such as: |H k± L , |H k± R , |D k± L , |D k± R . There is no direct information as to whether quantum correlations for proton states hold in the mixed sublattice. This point deserves further investigations.

Conclusion

The P 2 1 /a space group symmetry of the KH p D 1-p CO 3 (p ≈ 0.75) crystal determined with neutron diffraction shows that H and D nuclei are not individual particles located at definite sites and mutually exclusive. Wave-like protons and deuterons are equally distributed over all indistinguishable sites and, conversely, each site is occupied by all wave-like H and D nuclei with total occupancies p and (1p), respectively. This nonlocal behaviour allows the sublattice H p D 1-p to match the symmetry of the electronic structure that is unaffected by isotope substitution. The superposition of macroscopic states sounds as weird as a Cat in a superposition of |alive and |dead states.

To the question "Where are protons and deuterons?", our answer is that the particle-like representation, implicitly meant by "Where", is irrelevant. Only the wave-like, or matter-field, representation is permitted by the crystal structure. The theory is based on fundamental principles of quantum mechanics and it is free of any ad hoc hypothesis or parameter. The nonlocal nature of protons and deuterons could be therefore of general significance. Neutron diffraction by hydrogen bonded crystals with partial isotope substitution could be a well suited and powerful approach to study nonlocal behaviors.

Figure 1 :

 1 Figure 1: Schematic view of the crystalline structure of KHCO 3 at 300 K. The arrows point to the sites occupied at low temperature. Dashed lines through protons are guides for the eyes. The ellipsoids represent 50% of the probability density for nuclei.

Figure 3 :

 3 Figure 3: Projections on the (a, b) plane of the L (bottom) and R (top) configurations.

Table 1 :

 1 Unit cell parameters in Å and • units determined with singlecrystal neutron diffraction at 300 K. λ = 0.8305 Å, space group P 2 1 /a. The variance for the last digit is given in parentheses. = 2n and h0l for h = 2n) confirms the P 2 1 /a space group assignment. The unit cell parameters in Table1are identical to those of KHCO 3 .

		KH 0.75 D 0.25 CO 3 KHCO 3 [17]
	a( Å) 15.180(1)	15.180(1)
	b( Å) 5.620(1)	5.620(4)
	c( Å) 3.710(1)	3.710(4)
	β	104.67(1) •	104.67(5) •
	A preliminary inspection of intensities at the absent reflections (0k0 for