SP-MORPH processing landmarking Spectral Processing

Spectral mesh processing

This task is focused on the implementation of state-of-the-art algorithms for Spectral Mesh Processing (SMP) for their application to the analysis of craniofacial geometry. SMP is a novel branch of computer geometry and vision that focuses on the analysis of surfaces embedded in 3D based on the decomposition of the geometry in the frequency domain [Zhang 2010]. Analogous to the Fourier decomposition employed in signal processing, the seminal work by Taubin [Taubin 1995] suggested that the spatial frequencies of a 3D surface (described as a triangulated mesh) could be computed from the eigendecomposition of the Laplacian operator.

The wide interest attracted by SMP soon revealed an important shortcoming. In its original formulation, Taubin shows that the Fourier decomposition of a 1D signal arises from the eigendecomposition of the Laplacian (also in 1D) and extrapolates this result to 3D using the so called combinatorial graph Laplacian. This combinatorial operator perfectly emulates its 1D counterpart, but it is based only on the connectivity of the mesh, without taking into account its geometric information. While this is acceptable in 1D, where the signals are sampled at equidistant intervals (i.e. uniform sampling), it does not suffice for surfaces embedded in 3D. In the latter, the geometry is usually represented by triangulations (simplicial complexes) which are almost invariably non-uniform regarding their size and their angles.

 Figure 1: Some eigenfunctions of the Laplace-Beltrami Operator for the mean-shape of the Basel Face Model [Paisan 2009]. Values are color coded on the surface and rescaled to highlight the zero-crossings (dark blue curves).The order of the eigenfunctions increases from left to right.

More rigorously, the Laplacian can be considered a special case of the more general Laplace-Beltrami Operator (LBO), which is defined on a manifold invariant to its parameterization, taking into account only its Riemannian metric [Reuter 2006, Bronstein 2011]. The most widely accepted formulation for the LBO on 3D meshes is based on the cotangent weights [Meyer 2002], derived from the continuous definitions of the surface curvature assuming that the geometry is reasonably approximated by the triangulation at hand.

Further research has addressed the problem from different perspectives, mainly by employing Finite Elements or Discrete Exterior Calculus (DEC). The results are still based on the cotangent weights but it has been shown that a correction factor is needed for the operator to fulfil the mathematical properties of the (continuous) LBO [Levi 2006, Vallet 2008]. Another way to intuitively understand LBO eigenfunctions is to think of them as the generalization of spherical harmonics [Li 2006, Khairy 2008] to other non-spherical geometries. Fig. 1 shows examples of some eigenfunctions of the LBO in a facial surface.

The spectrum of the LBO is an isometric invariant, and it has been shown to be a powerful descriptor as a signature for 3D shape matching and classification [Reuter 2006, Ruggeri 2010]. Among several properties, its nodal sets (the zero-sets of its eigenfunctions) are known to divide the surface in (at most) as many regions as the order of the eigenfunction, defining curves that intersect at constant angles [Levy 2006]. Thus, they are strongly linked to the geometry of the surface and its intrinsic structure [Qiu 2008]. It has been observed that level sets of LBO eigenfunctions follow geometric features [Levy 2006], highlight protrusions [Reuter 2010] and reveal (global) symmetry, which is considered very important in dysmorphology [Hammond 2008, Qiu 2009, Atmosukarto 2010].

Discrete operators using cotangent weights

The Given a surface S (2-mainold) embedded in R3, we can locally approximate it at any point by its tangent plane, orthogonal to the normal vector n. Given the mean curvature H and an infinitesimal area element A with diameter diam(A), the following holds:

where Ñ is the gradient with respect to the 3D coordinates (at the point being considered). Meyer et al. define their LBO as L = 2 H n.

The above equations regard to the continuous setting. To extend them to the discrete setting, each vertex is considered in a neighbourhood AM, which is restricted to be within its 1-ring neighbourhood, where spatial averaging is performed. The exact calculation of AM depends on whether the triangles in the 1-ring are obtuse or not but, conceptually, it is an approximation of the Voronoi region of the considered vertex appropriately corrected to deal with bad conditioned cases.

The LBO of the surface with respect to the conformal space parameters u and v, conveniently chosen as the current surface discretization (i.e. using directly the surface triangles) becomes simply the flat space Laplacian Du,v. Then:

From this point, spatial averaging is derived from finite element techniques, using a linear finite element on each triangle (i.e. to interpolate between the vertices). Using Gauss’s theorem with appropriate algebraic manipulation [Meyer 2002] leads to:

where N1(xi) is the set of 1-ring neighbours of xi, and the angles αij and βij are always those opposite to the edge formed by xi and xj, as in the example of Fig. 2. Applying the same strategy, discrete operators for the normal vector, mean curvature and Gaussian curvature are also derived (and have been implemented).

 Fig. 2- Example of the 1-ring neighborhood of vertex xi and the angles on which to compute the cotangents when considering neighbor xj. Note that both angles opposite to the edge, αij and βij, are considered in the same way so their labels might be swapped to βij and αij without any consequence.

Discrete operators using the heat kernel

In contrast to the operators described in the previous section, we also implemented an alternative approach to compute the LBO presented by Belkin et al. [Belkin 2008]. It is also based on an approximation of integrals on a discretized surface but the authors combine it with the idea of approximating the heat flow on a mesh. The most attractive feature of the algorithm is that it guarantees point-wise convergence of the derived approximation to the correct LBO.

Let M={X, T} be a triangulated mesh, obtained by sampling a surface S, composed by two sets: the vertices X = {x} and the triangles T = {t}. Given a function f : X ® R, Belkin’s algorithm approximates the LBO of the function at vertex xi as:

where h is a parameter that controls the size of the diffusion kernel, A(t) is the area of a given triangle and V(t) is the set of its vertices. Note that the neighborhood considered for vertex xi depends on h and is not limited to the 1-ring neighborhood. Thus, it might seem that the connectivity of the mesh does not play a role in the above equation, but it does in two different ways:

1.      The distances between vertices, used in the exponential, should in reality be computed as Geodesics. The use of Euclidean distances is only an approximation intended for small enough neighborhoods.

2.      For each triangle we end up with a common coefficient (or weight) that is split evenly among its vertices. Thus, the resulting weight for a given vertex depends on its valence (the number of triangles in its 1-ring neighborhood).

Apart from providing a guaranteed convergence, this approach is more flexible than the one by Meyer et al. because it allows controlling the size of the neighborhood that is considered with the diffusion parameter. However, when assessing the correctness of the operator for some synthetic cases, we found that convergence is achieved only for extremely fine tessellations and neighborhoods that are considerably larger than the 1-ring of the considered vertex.

Apart from the derivations from cotagent weights (including a few variants with the same core algorithm) and diffusion kernel, there are also interesting alternatives based on finite element approximations of higher order that bypass the computation of the LBO and address directly the generation of its spectrum [Reuter 2009], which is indeed our ultimate goal.

Comparison of discrete Laplacian operators

Different algorithms to compute the Laplacian operator are typically compared by applying them to some simple function in domains such as planes or spheres. The reason for restricting the evaluation to such simple cases is the possibility to compute the correct Laplacian analytically. Actually, an interesting aspect of the evaluation is that we cannot directly evaluate the correctness of the operator; instead, we need to apply to operator to a function and evaluate the result.

For the tests reported below we used a non-linear function f defined on the sphere, f : S2 ® R = x2, as done by Belkin et al. [Belkin 2008]. Because the function is defined only in the unit sphere, the analytic solution for the Laplacian of f is:

where θ and ϕ are, respectively, the elevation azimuth angles in spherical coordinates. Recall that we are under the simplification of a unit sphere, hence that radius is 1 and we omit it everywhere. Thus, x = sin θ cos ϕ and we get:

Fig. 3 shows the analytic result, color-coded into the surface of the sphere.

 Fig. 3- Analytic result for the Laplacian of f = x2.

To evaluate the quality of the operators, we used 3 different discretizations of the sphere: the first one was generated by a uniform sampling in the azimuth-elevation plane, which implies a higher sampling density near the poles, and the other two were generated with near-uniform sampling in the spherical domain. However, the latter two were differed in the approximation used to generate the uniform sampling: therefore, they had some differences in the underlying triangulation. In all three cases, the spheres had approximately 8000 points. Visually, both operators look very similar to the analytic result displayed in Fig. 3 on all three discretizations. However, if we plot the absolute difference between the discrete operators and the analytic result we can clearly appreciate some differences, as shown in Fig. 4.

 Figure 4: Absolute errors for Δf using discrete operators based on cotangent eights (top row) and heat diffusion (bottom row). For each row, the left-most plot corresponds to the non-uniform sampling case; the remaining two are different approximations with uniform sampling.

The main conclusion we can extract is that the operator based on cotangent weights [Meyer 2002] is more accurate and, especially, more consistent. The latter is clear when comparing the two spheres with near-uniform sampling, where the operator based on heat diffusion creates error patterns that are quite different between these two cases. It is also interesting to note that, in the non-uniform case, the cotangent-based operator has bigger errors in the poles, where the sampling density increases and the triangles become more elongated. On the other hand, the method based on heat diffusion does not have problems in the poles, because the higher sampling density helps to achieve a better approximation to the integral.

Quantitative comparisons on spheres of variable density

Our results reported above confirm the statements from the related literature, which indicate that the operator based on cotangent weights produces more accurate and consistent results than the heat diffusion operator but only for idealized cases as it could be more affected by noise and uneven triangulations.

Indeed, both approaches have some problems related to the triangulation. The cotangent operator has problems for triangles with large difference in edge sizes (i.e. with big aspect ratios) and the heat diffusion operator is sensitive to the valence of vertices (the number of adjacent triangles). In the case of the heat diffusion operator, however, this problem is partially tackled by the fact that a relatively large neighborhood is considered for each point; integration over the neighborhood tends to compensate the small deviations introduced individually by each point but that also increases the size of the neighborhood that we need to take. Larger neighborhoods are undesirable for several reasons: they take longer to process, require move memory and reduce the locality of the operator, which might therefore not be able to capture enough detail.

An interesting aspect of the diffusion operator is that its dependency with the valence of the vertices could be easily removed, or at least attenuated. The reason for this dependency is that the contribution of each vertex is weighted by a factor proportional to the area of the incident triangles: actually, each vertex of a triangle gets exactly one third of the area-weight of that triangle. This is of course suboptimal and a very reasonable approximation has been provided by Meyer et al. [Meyer 2002] to construct a Voronoi tessellation of the mesh. Hence, apart from the two original methods (cotangent and heat diffusion), we tested a third one (diffusion + Voronoi) which is the (rather trivial) extension of the heat diffusion method with the Voronoi areas from the cotangent-based method.

Fig. 5 shows the normalized L2 error on uniformly sampled spheres of variable number of points and with different neighborhood size. This error measure corresponds to the following expression:

where f is the analytically computed Laplacian of function f (on the given surface) and Δf is the discretized approximation of the same Laplacian. In all experiments of this report I used the function f(x, y, z) = x2 as done in one of the experiments from Belkin et al. [Belkin 2008]. We can see that:

·        The cotangent operator is always more accurate than the other two methods.

·        The error of the heat diffusion operator depends strongly on the neighborhood size in a valley-like shape. Thus, there is an optimal neighborhood size, in terms of error, which is dependant of the sampling density of the surface.

·        The use of Voronoi areas improves the performance of the heat diffusion operator, reducing both the error and the size of the neighborhood for which the lowest error is achieved.

·        Increasing the sampling density produces better results in all methods. Let us emphasize that we can increase the resolution in an exact manner because we know the equation of the underlying surface. This cannot be reproduced on real data.

 Figure 5: Normalized errors of Δ f on spheres with uniform sampling for different discrete approximations: cotangent cotangent (dash-dot lines), diffusion (dashed lines) and diffusion + Voronoi (solid lines). Notice that the cotangent operator does not depend on the heat diffusion parameter (horizontal axis), but we display it for the same range of the diffusion-based operators to facilitate comparison.

Additionally, to test the influence of the triangulation on the operators, we repeated the experiment with non-uniform sampling, as follows:

·        Points were uniformly sampled on latitude-longitude coordinates, generating a much denser sampling on the poles than on the Equator (Fig. 6, left). Results are displayed in Fig. 7, left. We can see that the error of the cotangent operator is now similar to the one of the diffusion operator. The use of Voronoi areas produces the best results in all cases.

·        Points were randomly sampled on the surface of the sphere, generating a very uneven triangulation (Fig. 6, right). Results are displayed in Fig. 7, right. We can see that now the cotangent operator is worse than the diffusion-based approach but there is also an increase in the size of the neighborhood needed for optimal performance. The use of Voronoi areas has little impact in performance.

In summary, the results we have so far suggest that we should use the Voronoi-area correction for the diffusion operator, as it can bring an advantage for some cases (both in terms of accuracy and size of neighborhood).

 Figure 6: Examples of spheres with non-uniform triangulations. Left: points sampled uniformly in the azimuth-elevation domain. Right: points sampled randomly (with uniform distribution for each coordinate) on the surface of the sphere. Figure 7: Normalized errors of Δ f on spheres with non-uniform triangulation generated by: a) uniformly sampling on the azimuth-elevation domain (Left) or b) random sampling with uniform distribution on each coordinate (Right). Three discrete operatorsare displayed: cotangent cotangent (dash-dot lines), diffusion (dashed lines) and diffusion + Voronoi (solid lines).

Laplacian spectrum and facial symmetry

We investigated the behavior of the Laplacian spectrum for variable degrees of asymmetry of the human face. To this end, we started by constructing a perfectly symmetric facial template obtained by averaging of 100 facial scans (both the original and left-right mirrored version). Averaging was possible thanks to the use of Least Squares Conformal Maps (LSCM) to have homologous representations of all input surfaces (details here).

Having a symmetric template as described above makes it easy to introduce artificial distortions to the facial surface. As an example, we produced an enlargement of the left side of the face by scaling the distance of all left-side points to the symmetry axis (mid-sagittal plane) by a given percentage. Fig. 8 shows the distortion resulting from such a scaling by a 20% factor.

 Fig. 8- Fully symmetric average artificially made asymmetric by expanding the left side by 20%. The color coding provides a qualitative indication of the expansion magnitude with respect to the original

The re-sampling in the 2D domain was designed to produce points that correspond to anatomically symmetric pairs. Recall that, in general, the human face is not strictly symmetric and, in anatomical terms, the hypothetical axis of bilateral symmetry is not contained in a plane. However, most research in facial symmetry relies on the simplification of taking the mid-sagittal plane as the plane of bilateral symmetry, which is then used to measure the deviations of one side with respect to the other [Claes 2011]. This is related to the direct use of 3D coordinates for analysis, while the facial surface can be considered a 2D manifold embedded in 3D. In contrast, the spectrum of the Laplacian relates to the structure of the manifold itself and can be used for analysis without considering reflection over the mid-sagittal plane. So far, however, computation of this spectrum is dependent on the extent of facial region that is captured, which is extremely variable. Hence, the need and utility of the mesh normalization that we perform based on LSCM.

The spectrum was computed by performing eigen-decomposition of the Laplacian operator of each of the two templates mentioned above. The results in this report correspond to the cotangent-based Laplacian (partial inspection of the results with the diffusion-based Laplacian suggests that results are not too different).

An important problem to compute the spectrum is that the, given a discrete Laplacian operator L, we cannot guarantee that its eigenvectors u will be real (i.e. in the general case they are complex). This happens because L is not strictly symmetric, as it should be if it were a faithful representation of its continuous equivalent. There are two main solutions to this problem [Zhang 2010]: i) force the operator to be symmetric, for example by using = (L+LT)/2; ii) decompose the Laplacian into a symmetric matrix S and a mass matrix B, which simplifies obtaining real eigenvectors, especially if B is diagonal. Evidently, the second solution is more principled. The idea is to express the Laplacian as:

Both the cotangent and diffusion approaches allow to easily separate L into a symmetric matrix S and a diagonal matrix B-1. Doing some algebra it is possible to show that an eigen-decomposition of L can be obtained from the eigendecomposition of B-1/2SB-1/2. Specifically, the eigenvalues of B-1/2SB-1/2 are also eigenvalues of L and an eigenvector v of B-1/2SB-1/2 can be transformed into an eigenvector u of L by:

Fig. 9 show the first 14 non-zero eigenvectors for the symmetric template and three manipulated versions of it with increasing degree of asymmetry. To highlight the patterns resulting from the eigenvectors, we have not plotted the values of the eigenvectors v themselves, but a non-linear scaled version |u|1/4, where the exponential is applied element-wise. The effect of this scaling is twofold: i) if highlights the zero-crossings (the dark-blue lines); ii) it simplifies the perception of symmetry by converting anti-symmetric values to symmetric ones.

There are two main observations: i) progressively adding asymmetry resulted in also in a progressive modification of the spectral components; ii) there were some shifts in the positions of certain components, e.g. 9 with 10 and 11 with 12 in Fig. 3, presumably due to very close eigenvalues, as discussed in the work by Ovsjanikov et al [Ovsjanikov 2008].

The above is encouraging regarding the potential to measure symmetry in the spectral domain, but it has to be put in context of the related work, especially since before applying the spectral decomposition we have to put the surfaces in dense correspondence. To the best of our knowledge, most efforts so far have focused on the detection of the plane of bilateral symmetry or on the alignment of a given shape with the mirrored (reflected) version of itself, which is an over-simplification of the problem. We believe that this spectral-based approach could help achieving a more principled analysis.

 Figure 9: Non-zero eigenvectors 1 to 7 (5 left-most columns) and 8 to 14 (5 right-most columns) of the Laplacian, sorted in ascending order from top to bottom, for a symmetric face which is artificially made asymmetric by expanding its left side; the percentage varies from left to right as follows: 0% (fully symmetric), 2%, 6%, 12%, 20%.

References

[Atmosukarto 2010] I. Atmosukarto, K. Wilamowska, C. Heike and L.G. Shapiro. 3D object classification using salient point patterns with application to craniofacial research. Pattern Recognition, 43(1):1502–1517, 2010.

[Belkin 2008] M. Bekin, J. Sun and Y. Wang. Discrete Laplace operator on mesh surfaces. Proc. ACM Symp, on computational geometry, pp 278-287, 2008

[Bronstein 2011] M.M. Bronstein and A.M. Bronstein. Shape recognition with spectral distances. IEEE Transactions on Pattern Analysis and Machine Intelligence, 33(5):1065–1071, 2011.

[Claes 2011] P. Claes, M. Walters, D. Vandermeulen and J.G. Clemen. Spatially-dense 3D facial asymmetry assessment in both typical and disordered growth. Journal of Anatomy, 219(4):444–455, 2011.

[Hormann 2007] K. Hormann, B. Levy and A. Sheffer. Mesh Parameterization: Theory and Practice. ACM SIGGRAPH course notes, 2007.

[Hammond 2008] P. Hammond, C Forster-Gibson, A.E. Chudley, et al. Face–brain asymmetry in autism spectrum disorders. Molecular Psychiatry, 13, 614–623, 2008.

[Khairy 2008] K. Khairy and J. Howard, Spherical harmonics-based parametric deconvolution of 3-D surface images using bending energy minimization. Medical Image Analysis, 12(2): 217–227, 2008.

[Levy 2006] B. Lévy. Laplace-Beltrami eigenfunctions. Towards and algorithm that ‘understands’ geometry. Proc IEEE Conf Shape Modeling and Applications, 2006.

[Li 2006] H. Li and R. Hartley. New 3D Fourier descriptors for genus-zero mesh objects. In Proc. Asian Conference on Computer Vision. LNCS vol. 3851, pp 734–743, 2006.

[Meyer 2002] M. Meyer, M. Desbrun, P. Schröder and A.H. Barr. Discrete differential-geometry operators for triangulated 2-manifolds. Proc Visualization and Mathematics III, pp. 35–57, 2002.

[Ovsjanikov 2008] M. Ovsjanikov, J. Sun, and L. Guibas. Global intrinsic symmetries of shapes. Computer Graphics Forum, 27(5):1341–1348, 2008.

[Qiu 2008] A. Qiu, L. Younes and M.I. Miller. Intrinsic and extrinsic analysis in computational anatomy. Neuroimage, 39(4):1803–1814, 2008.

[Reuter 2006] M. Reuter, F.E. Wolter and N. Peinecke. Laplace-Beltrami spectra as ‘Shape DNA’ of surfaces and solids. Computer Aided Design, 38(4):342-366, 2006.

[Reuter 2009] M. Reuter, S. Biasotti, D. Giorgi, et al. Discrete Laplace-Beltrami Operators for Shape Analysis and Segmentation. Computers & Graphics, 33(3): 381-390, 2009.

[Ruggeri 2010] M.R. Ruggeri, G. Patane, M. Spagnuolo and D. Saupe. Spectral-driven isometry-invariant matching of 3D shapes. International Journal of Computer Vision, 89(2-3):248–265, 2010.

[Taubin 1995] G. Taubin. A signal processing approach to fair surface design. Proc SIGGRAPH, pp. 351-358, 1995.

[Vallet 2008] B. Vallet. Function bases on manifolds. PhD Thesis, Institut National Polytechnique de Lorraine, France, 2008.

[Wang 2007] S. Wang, Y. Wang, M. Jin et al.  Conformal Geometry and Its Applications on 3D Shape Matching, Recognition, and Stitching. IEEE Transactions on Pattern Analysis and Machine Intelligence, 29(7):1–12, 2007.

[Zhang 2010] H. Zhang, O. van Kaick and R. Dyer. Spectral mesh processing. Computer Graphics Forum, 29(6):1865–1894, 2010.