{'fraction_non_alphanumeric': 0.03324665676077266, 'fraction_numerical': 0.011515601783060922, 'mean_word_length': 4.520892078954114, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Arbitrary order principal directions and how to compute them

Curvature principal directions on geometric surfaces are a ubiquitous concept of Geometry Processing techniques. However they only account for order 2 differential quantities, oblivious of higher order differential behaviors. In this paper, we extend the concept of principal directions to higher orders for surfaces in R^3 by considering symmetric differential tensors. We further show how they can be explicitly approximated on point set surfaces and that they convey valuable geometric information, that can help the analysis of 3D surfaces.

1. Introduction. The computation of informative tangent vector quantities on surfaces is a widely studied topic. The most standard vector quantities one can consider are the principal directions that can be estimated via differential geometry tools. However these do not necessarily serve all Computer Graphics purposes: in case of umbilical surface parts, or at monkey saddles, this vector field becomes locally irrelevant. Further, it only accounts for an edge like structure of the shape which is limited. Our approach considers per point vector quantity estimation, adopting a differential analysis point of view. We strive to go beyond the differential order 2 usually considered when analysing surfaces and extend the definition of principal directions to higher orders. We show experimentally why this definition, beyond its theoretical interest is appealing for surface analysis.
To summarize, our contributions are as follows: • The mathematical definition of arbitrary order principal directions and their link with symmetric tensor eigenvalues. • A theoretical analysis of their properties.
• A practical way of computing the directions on a sampled surface.
2. Related Work. In this section we review recent works on tangent vector quantities that can be set or estimated on a meshed or sampled surface be they guided by differential properties or designed by a global optimization process with user-prescribed constraints.
Differential quantities estimation. Estimating differential quantities on surfaces has been at the core of Geometry Processing Research. However, surface analysis restricts very often to order 1 and 2 differential properties and has seldom tackled higher order properties. Among order 2 quantities, the most famous one may be the Laplace-Beltrami operator, whose design has gathered a lot of works both from a theoretical analysis (e.g. [19,30]) and practical analysis through the Manifold Harmonics Basis (e.g. [28]). Related to the Laplace-Beltrami operator are the principal curvatures and curvature directions (or equivalently the curvature tensor) estimation, either on a point set surface [13,12] or on meshes [6], with applications to curvature lines tracking among others.
It is however possible to get access to high order derivatives of the local surface using local regression, in the context of Moving Least Squares [17]. Among those methods, Osculating Jets [5] express the surface locally as a truncated Taylor expansion wrt to a local planar parametrization. The coefficients can be estimated through a linear system solve and give then a direct access to high order cross derivatives. Interestingly, this approach proved that the error on order k differential quantities estimation in a neighborhood of radius r using a Taylor expansion of order K was in O(r K−k ). In other words, and quite counter-intuitively, to increase the accuracy of an order 2 estimation, one should still consider a large Taylor expansion order. Several other basis have been proposed following this trend, including the Wavejets basis [2] which is less sensitive to the choice of the local reference frame in the parametrization plane. When the surface is described as a point set, the regression relies on Iteratively Reweighted Least Squares in the presence of noise and/or outliers. Going further than order 2 Rusinkiewicz [26] introduces a way to compute curvatures and curvature derivatives but does not go beyond this order. All these methods essentially perform per point estimation and do not take into account any global regularization constraints. For example, on planar or spherical surfaces curvature directions are erratic in the absence of smoothness constraints, which is required by many computer graphics applications. Hence researchers have turned to the design of consistent vector fields more suited to some application purposes.
Tangent Vector Field Design. The problem of tangent vector field design is to compute a smooth tangent vector field with user prescribed constraints at given points of the surface, while optimizing some regularity criterion. We only review some of the seminal papers of this field and we refer the reader to [29] for an extensive review. Many vector field design methods focus on smooth N -symmetric vector fields, also known as rotationally-symmetric direction fields (N-RoSy). N-symmetry directions are defined as the sets of directions invariant by 2π/N rotations. Ray et al. [25] generalize the notion of singularity and singularity index to these fields, and provide a way of controlling singularities on surface meshes. Lai et al. [16] focus on casting the vector field design as a Riemannian metric design problem. Further smoothness constraints [7,14] and global symmetry enforcing constraints [21] were proposed. N-RoSy were also extended to non necessarily orthogonal nor rotationally symmetric vector fields [8] with appropriate differential operators [9] and application to Chebyshev nets computation [27].
The generalization of the Laplace-Beltrami operator to tangential vector fields and the subspace raised by its eigenvectors up to a given order may be used for regularizing vector field design [3]. Following this development of subspace methods for tangent field design [4], Nasikun et al. [20] consider tangent vector field design and processing via locally-supported tangential fields leading to fast approximation and design algorithms. Application of Tangent Vector Fields. Applications of tangential vector fields range from texture mapping and texture synthesis on surfaces (e.g. [31,15]) to fluid simulation (e.g. [1]). Shape reconstruction and quad-meshing have been tackled by combining a position field and a N-RoSy [11] yielding an extremely fast interactive algorithm.
In this paper we are also interested in computing per point sets of directions that are neither necessarily orthogonal nor rotationally symmetric but these directions stem from arbitrary order differential properties of the surface. Hence smoothness is obtained by continuity of the underlying mathematical surface.
3. Arbitrary order symmetric tensor. Our work makes extensive use of symmetric tensors and the theory developed by Qi for their spectral analysis [22,23,24]. Notice that in Qi's work, a distinction is made between a tensor and a supermatrix, that is a tensor's representation in a given basis. For clarity sake, here we rather use the tensor term for both the object and its representation in a basis.
From now on, we will always consider m = 2 since we are interested in tensors of differential properties related to surfaces of dimension 2. In this case, a symmetric tensor of order k can be seen as a k − 1-dimensional array of vectors of length 2. By convention, we define any vector to be a symmetric tensor of order 1. Given a vector v = (x, y) T , we note v k the symmetric tensor of order k generated by multiplying v k times using the Cartesian product. In particular, we set v 0 = 1 by convention.
Multiplying a symmetric tensor by a vector v produces a symmetric tensor of order lowered by 1. Let T be a symmetric tensor of order k, it is composed of two symmetric tensors of order k − 1 T x and T y and can be written T = (T x , T y ). Then T v is the symmetric tensor of order k − 1 such that : The product T v reduces the order of T by contracting on an arbitrary index. Since T is symmetric, any index used for the contraction yields the same list of numbers in T v. While eigendecomposition of matrices is a well understood theory with important applications in Geometry Processing among others, the generalization to arbitrary order tensors is highly nontrivial. Qi introduced a new definition extending eigenvalues and eigenvectors to symmetric tensors [22,23,24], as follows: [22] Given T a symmetric tensor of order k, if there exists λ ∈ C and a vector v ∈ R 2 such that: Then λ is called an E-eigenvalue of T and v is called an E-eigenvector of T . The set of λ satisfying (3.2) are the roots of a polynomial called the E-characteristic polynomial.

Arbitrary order principal directions.
Differential tensor. Tensors can be used to write Taylor expansions. As an example, one can write the two first terms of a bivariate Taylor expansion. Given v = (x, y) T an arbitrary vector, n the normal at (0, 0) and H the Hessian of a function defined on R 2 with values in R and twice continuously differentiable at (0, 0): Note that H is symmetric, and so is n since its order is 1. This expression can be generalized to higher orders Taylor expansions using tensors.
The Taylor expansion of order K of a function f from R 2 to R is: The differential tensor is now defined as: For a function f defined on R 2 with values in R, k times differentiable, the k th order differential tensor T k at (0, 0) is a k th order 2-dimensional tensor whose coefficients are: If f is differentiable then the order in which the differentiation is done is irrelevant and thus T k is symmetric. Using Definition 4.1 and v = (x, y) T , we have: and: Hence using Equation 4.2, we get the Taylor formula for a K times differentiable function: The following Lemma shows that the gradient of each of the terms involved in the Taylor expansion can be obtained by contracting the corresponding tensor k − 1 times, i.e. one time less than in the expansion. This will then allow us to search for extrema at different orders.
Proof. For k = 1, v 1 = v = (x, y) T and: By induction, the property is true for all k ≥ 1.
Theorem 4.3. Given a K times continuously differentiable function f and 1 < k ≤ K, T k is the real symmetric k th order differential tensor of f and the set of vectors v = (r cos θ, r sin θ) T such that ∂ ∂θ T k v k = 0 and v = 1 are real E-eigenvectors of T k , i.e.: Proof. First one can notice using equation 4.5, with v = (r cos θ, r sin θ) T that: Since r ∂ ∂r = x ∂ ∂x + y ∂ ∂y we get: Differentiating w.r.t angle θ to look for extrema: which yields the following equations: Since we look for real unitary vectors, we add the constraint that r = v = 1. Moreover, setting λ = T k v k , we get T k v k−1 = λv and v is a real E-eigenvalue of T k . The reverse holds using the same equations.
Definition 4.4. Given a K times continuously differentiable function f defined on R 2 with values in R, the principal directions of order k (1 < k ≤ K) of f at (0, 0) are defined as the real E-eigenvectors of its k th order differential tensor T k .
One should notice that Qi et al. defined several types of eigenvalues and eigenvectors [22,23,24]. In particular, if an E-eigenvector v is real and if its corresponding E-eigenvalue λ is also real, then v is a Z-eigenvector and λ is a Z-eigenvalue. Then our E-eigenvalues are also Z-eigenvalues in this setting. Figure 2 illustrates the principal directions of order 3 and 8 for some illustrative functions at (0, 0).
The above form is not very amenable to find the zero-crossings of the derivative of T k v k with respect to θ. Instead we propose to use its expression in the Wavejets basis (consisting of all functions (r, θ) → r k e inθ for k ∈ N, −k ≤ n ≤ k) [2], hence using polar coordinates: with φ k,n the Wavejets decomposition coefficients. Among other properties, φ k,n = φ * k,−n and φ k,n = 0 is k and n do not share the same parity (see [2] for more details).
By identification of the coefficients in front of the powers of r with the coefficients of the Taylor expansion using tensors, we have:   n(Im(φ k,n ) cos(nθ) + Re(φ k,n )sin(nθ)) = 0 A more convenient form to find roots, for example using Newton's method.
So far, we defined the principal directions as the eigenvectors of a tensor which we linked to the extrema of a function g k (θ). The eigenvalues of the tensor can be also linked to this function. Calling θ e an angle corresponding to an extremum of g k , the corresponding eigenvalue λ e can be recovered as: λ e k! This follows directly from the last equality in 4.13.
Definition 4.6. Among all principal directions, we call Maximum principal directions (resp. minimum principal directions) the directions that correspond to local maxima (resp. local minima) of g k (θ) = k n=−k φ k,n e inθ = T k v k k!r k with v = (r cos θ, r sin θ).

Properties of the principal directions.
Constraints on the principal directions. The functions g k (Definition 4.6) can be rewritten as g k (θ) = φ k,0 + 2 k n=1 (Re(φ k,n ) cos(nθ) − Im(φ k,n )sin(nθ)). From the periodicity of cosine and sine functions, we deduce that: • If k is even then g k (θ) = g k (θ + π), hence if θ 0 corresponds to a maximum principal direction, θ 0 + π is also a maximum principal direction. • If k is odd, g k (θ) = −g k (θ + π), hence if θ 0 corresponds to a maximum principal direction then θ 0 + π corresponds to a minimum principal direction. Number of principal directions. There are at most 2k principal directions of order k, since finding the zeros of ∂g k ∂θ amounts to finding the 1-norm roots of a real polynomial of order 2k (obtained by multiplying Equation 4.16 by e ikθ ). Since two maximum principal directions should be separated by one minimum principal direction and conversely, the number of maximum principal directions and the number of minimum principal directions should be equal and their maximum number is thus k. Following the parity constraints on the location of maxima and minima above, for an even order, the number of maximum principal directions is necessarily even. For similar reasons, for an odd order, the number of maximum principal directions is necessarily odd.
Order 2 principal directions. The principal directions of order 2 correspond to the classical principal curvature direction, however the maximum (resp. minimum) principal directions might not correspond to the maximum (resp. minimum) curvature directions.
Regularity and link with N-RoSy. Order k principal directions can turn into a k-RoSy, if and only if φ k,n = 0 for all n = ±k. The principal direction distribution can however not be arbitrary: this can be seen by considering order 3 principal directions, fixing their 3 angles and trying to solve for the coefficients yielding a 0 derivative for these 3 angles. This amounts to considering 4 unknowns corresponding to the real and imaginary parts of the coefficients φ 3,1 and φ 3,3 (since φ 3,n = φ * 3,−n ), linked by 6 equations given by θ i and θ i +π with i = 1 · · · 3. A rank analysis yields that the system is sometimes invertible and hence yields only the trivial solution of all coefficients set to 0. Sometimes the system has rank 2 or 3 depending on the chosen angles and hence yields a nontrivial solution. Hence not all kind of irregularities can be represented by the principal directions of the tensor. On Figure 2 we show a monkey saddle and an octopus saddle, whose respective order 3 and 8 principal directions correspond to 3 and 8-RoSy when considering only maximum (alternatively only minimum) principal directions.
6. Directions Computation per point. We now propose to compute these directions on point set surfaces by using a local parametrization around each point. In most Geometry Processing applications, surfaces are known only through a set of points, sometimes connected into a mesh, and the surface in between should be inferred by regression to estimate principal directions.
In practice, to compute principal directions, we perform a local Wavejets surface regression with a high enough order (K = 10 in our experiments). More precisely, let p be a point of the point set, let (p i ) i=1···N be its neighbors within some userdefined radius r. We compute a local tangent plane by robust PCA and deduce a local parametrization plane on which we choose an arbitrary tangent direction which serves as the origin direction for computing the local polar coordinates (r i , θ i , z i ) for each p i . Then we solve for the Wavejets coefficients φ k,n by minimizing: and σ = r 3 . This Gaussian weight avoids sharp boundary effects and makes the Wavejets estimation smoother in the ambient space.
Depending on the type of data, we can use the 1 -norm (p = 1) when there is noise and outliers, or the 2 norm (p = 2) if the data has low noise. As shown by Levin [17,18], the coefficients obtained by Moving Least Squares are continuously differentiable if a 2 norm is used. The use of the 1 norm does not provide such a guarantee. Hence depending on the required smoothness one should use a different norm in the estimation procedure of Equation 6.1.

Experiments.
Synthetic data. To illustrate the behavior of arbitrary order principal directions, we compute order 3 principal directions on a synthetic surface controlled by its Wavejets coefficients (Figure 3). The number of maximum directions is either 1 or 3 (hence either 1 or 3 minimum principal directions). On Figure 4 we show order 2 and 3 principal directions on a smooth synthetic surface evolving from a ridge to a smooth T-junction. One can see that order 3 takes slowly over order 2, with a preferred direction. Figure 5 illustrates the behavior of orders 2 to 7 principal directions on a sharp feature created by 5 intersecting planes. No order alone captures all intersection directions, but orders 5 and 7 contain the correct directions. Interestingly, order 7 degenerates to create only 5 maximum principal directions.
On Figure 6, we show the behavior of order 2 and 3 principal directions along the edges and corners of a synthetic cube. The length of the directions reflects the amplitude of the extremum. Order 3 accounting for some antisymmetry vanishes for edge points (which are symmetric) and order 2 vanishes on the corners (which are antisymmetric).
Real world models. Figures 1 and 7 show some of the principal directions computed by our method on the Armadillo model sampled with 5M points. The principal directions orders are chosen manually according to local geometric features. As expected order 2 accounts well for valleys, order 4 for valley crossings, and order 3 for some antisymmetry and monkey-saddle like features.      While it is obvious that sometimes order 2 is enough (side of the nose), order 3 is meaningful between the eyebrows and around the lips. Figures 9 and 10 show some principal directions computed on some more geometric shapes. Here order 3 becomes especially meaningful near corners, at the intersection between 3 smooth surfaces.  Robustness. To show the robustness of our principal directions estimation we add some Gaussian noise to the data. Figure 11 shows examples of principal directions estimation of order 2 and 3 on a cube with various noise levels. Importantly enough, this robustness does not stem from the principal direction decomposition itself but from the coefficients estimation. Once the coefficients are estimated the principal directions are obtained through function maximum and minimum computations, which can only introduce numerical errors. 8. Limitations. Our definition of principal directions is an extension of the principal curvature directions to higher orders, and as such are continuous on generic surfaces. For surfaces that are umbilical at every order such as a sphere or a plane, the directions will simply vanish (since no extremum will be found). Further smoothness constraints could be set locally to enforce some consistency, however this goes beyond the scope of this paper. Our approach also shares a limitation common to many local estimation methods: a radius should be chosen so that the analysis is local enough but also such that the neighborhood encloses enough points. The radius has indeed a direct impact on what is captured by the principal directions as illustrated on Figure  12. Importantly enough, our method does not perform better for curvature principal directions estimation than Osculating Jets [5] or APSS [10]. Our contribution lies rather on the extension to arbitrary orders than on the estimation accuracy itself. Finally, while it is appealing to consider that order-3 principal directions capture three ridges meeting at a single point, some precautions ought to be taken: if the ridges meet at a perfect T-junction, the maximum (or minimum) principal directions will not capture the 3 ridges directions because two maximum order-3 principal directions cannot be opposite. Figure 13 illustrates this behavior (see also Figure 5). 9. Conclusion. In this paper we introduced an extension of principal curvature directions to arbitrary differential orders and showed the link with the eigenvectors and eigenvalues of differential tensors. We showed that these new intrinsic direction fields are relevant on several shapes and can be computed efficiently even with sharp features. As a future work, some global smoothness constraints could be added, to enforce some surrogate direction computations on surfaces where the directions vanish. Many more applications of this new type of principal directions remain to be explored.