Chapter 3 Surface Classification Methods

Among the goals of this thesis is to devise an algorithm that automatically detects parts of a discretized surface, originating from a simulation workflow, that benefits from a higher resolution of the simulation domain. This requires mathematically analyzing the surfaces so that the regions of interest can be efficiently detected.

In chapter 2 three different surface representations and how to obtain one of their basic geometric properties (i.e., the surface normal) have been discussed. Surfaces have more geometric properties than just the surface normal vector (e.g., the tangent plane). Another property of orientable surfaces, which will be discussed later, is the surface curvature. The surface curvature describes how the direction of the normal vector changes when it is moved from one point on the surface to a neighboring point in an infinitesimal area around the point. Thus, the surface curvature can be used to classify parts of a surface.

In this chapter the mathematical concept of surface curvature on continuous surfaces is introduced. This is followed by a discussion on how these concepts can be discretized for the previously discussed surface representations. For implicit surfaces several methods for calculating surface curvatures are reviewed. The last section in this chapter discusses how discrete surfaces can be analyzed without prior knowledge of surface curvatures.

3.1 Curvature of Continuous Surfaces

The constructions and definitions presented within this section are extracted and simplified from [25]. Only the results needed for further investigations in this work are discussed and proofs are omitted. The goal of this section is to give a rigorous understanding of how the curvatures of a continuous surface are constructed and which properties of the surface they express.

3.1.1 2D Surface (Curves)

First 2D surfaces are discussed. In differential geometry 2D surfaces are referred to as parametrized curves and are defined as follows [25]:

For example consider the map \(\gamma (t) = (t^3, t^2)\), the parametrized curve created by this expression is shown in Figure 3.1.

This curve has a point \(\mathbf {p} = (0,0)\) for which \(\gamma '(\mathbf {p})=0\), these points are called singular points.

For the further investigations in this section it is important that a parametrized curve \(\gamma (t)\) does not have any singular points. Thus, the definition of a parametrized curve is extended to a regular parametrized curve.

Excluding curves with singular points allows for the definition of the length of a curve. The arc length of a regular parametrized curve is defined as follows:

A regular curve that fulfills \(\|\gamma '(t)\| = 1\) for all \(t \in I\) is called a curve parametrized by arc length. For such curves the arc length is given by \(s(t) = t - t_0\).

The Tangent Line and Curvature of a Curve Parametrized by Arc Length

In each point \(\gamma (t)\) on a curve parametrized by arc length there exists a straight line, defined by \(\gamma '(t)\) which has a length of \(\|\gamma '(t)\| = 1\). This vector passes through the point \(\gamma (t)\) and is called the tangent line (\(\gamma '(t)\)). In singular points it is not possible to define the tangent line and so, it is not possible to discuss geometric properties of such curves.

Consider a curve parametrized by arc length \(\gamma \) which in a point \(\gamma (t)\) fulfills \(\gamma ''(t) \neq 0\), in this case a unit vector \(\vec {n}(t)\) is well-defined by \(\gamma ''(t) = \|\gamma ''(t)\|\vec {n}(t)\). Furthermore, differentiating \(\gamma '(t) \cdot \gamma '(t) = 1\) yields \(\gamma ''(t) \cdot \gamma '(t) = 0\), which proves that \(\gamma ''(t)\) and \(\gamma '(t)\) are orthogonal. Therefore, the vector \(\vec {n}(t)\) is normal to \(\gamma '(t)\) and is called the normal vector in \(\gamma (t)\).

Since \(\gamma '(t)\) has length \(1\), the norm of the second derivative \(|\gamma ''(t)|\) measures the rate of change of the angles of the neighboring tangent planes at the tangent plane in \(\gamma (t)\). In other words the curve parametrized by arc length \(\gamma (t)\) pulls away from the tangent line in \(\gamma (t)\) at a rate determined by \(\|\gamma ''(t)\|\).

A straight line \(\gamma (t) = \vec {u} t + \vec {v}\) with constant vectors \(\vec {u}, \vec {v}\) has a curvature of \(0\). Inversely if \(\|\gamma ''(t)\| = 0\) is integrated the curve \(\gamma (t)\) has to be a straight line. Thus, the curvature of a curve can be used as a metric to distinguish between parts of a curve that are flat and parts that are not, e.g., curved.

Furthermore, it can be shown that each regular parametrized curve can be expressed as a curve parametrized by arc length  [25]. Thus, the concept of the curvature can be expanded to any regular parametrized curve. The curvature of the regular parametrized curve \(\gamma (t)\) is given by the curvature of the reparametrized (i.e., parametrized by arc length) curve \(\gamma ^*(t)\).

3.1.2 3D Surface

In this section a similar construction to the one presented in the previous section is given for 3D surfaces. However, defining geometric properties like the tangent plane or the curvature for 3D surfaces (which are also referred to as regular surfaces) requires additional considerations. The first major difference is that a regular surface is not defined as a function but as a set [25].

To further elaborate the meaning of Definition 3.1.5 consider Figure 3.2.

A regular surface is a subset of the plane that gets deformed in \(\mathbb {R}^3\). The deformation is handled in such a way that the resulting surface has no sharp points, edges or self-intersections. This definition allows for the definition of a tangent plane in each point on the surface and ensures that the surface is smooth enough to apply methods from calculus. The map \(\mathbf {x}\) is called a parametrization of the surface \(S\). A surface \(S\) may have several parametrizations, however, certain properties of the surface do not depend on the parametrization but only on the set \(S\).

Definition 3.1.5 is a technical definition of a regular surface, which is useful for the theoretical investigation of surfaces and their properties. In this work, the geometric properties of discretized surfaces are investigated, thus, two less technical descriptions of regular surfaces are given which result from Definition 3.1.5.

Proposition 3.1.6 can be interpreted in a more descriptive way. The set of points in \(\mathbb {R}^3\) that fulfill an implicitly defined function \(f(x,y,z) = a\), which has a non-vanishing gradient in each point of the image, describes a regular surface. Consider, for example, the following function \(f(x,y,z) = x^2 + y^2 + z^2 -1\), all points except for \(\mathbf {p} = (0,0,0)\) have a non-vanishing gradient, however, \(f^{-1}(0)\) is not in the preimage of \(f\), thus it is a regular surface (i.e., a sphere with radius \(1\)).

Proposition 3.1.7 guarantees that in a neighborhood around each point on a regular surface \(S\) the surface can be expressed as the graph of a differentiable function.

The Tangent Plane and Curvatures of a Regular Surface

Consider a regular surface \(S\) and a point on the surface \(\mathbf {p} \in S\), and a regular parametrized curve \(\gamma :(-\epsilon , \epsilon ) \rightarrow S\) with \(\gamma (0) = \mathbf {p}\). Since \(\gamma \) is a regular parametrized curve the tangent line \(\gamma '(0)\) exists. Let \(\mathbf {x}:U \subset \mathbb {R}^2 \rightarrow S\) be a parametrization of \(S\) and let \(\mathbf {q} \in U\). Then the set of all tangent lines of regular parametrized curves on the surface \(S\) passing through \(\mathbf {p}\) span a 2D subspace. The plane \(d\mathbf {x}_q\), that passes through \(\mathbf {x}(\mathbf {q}) = \mathbf {p}\), does not depend on the parametrization and is called the tangent plane to \(S\) denoted by (\(T_{\mathbf {p}}(S)\)) at the point \(\mathbf {p}\). The basis vectors that span \(T_{\mathbf {p}}(S)\) are written as \(\mathbf {x}_u\) and \(\mathbf {x}_v\) and are determined by the parametrization. By choosing a parametrization \(\mathbf {x}\), a unit normal vector in the point \(\mathbf {p} \in \mathbf {x}(U)\) can be defined as follows

\begin{align} \vec {n}(\mathbf {p}) = \frac {\mathbf {x}_u \times \mathbf {x}_v}{\|\mathbf {x}_u \times \mathbf {x}_v\|}(\mathbf {p}), \end{align} which is called the normal vector of the surface \(S\) in \(\mathbf {p}\) (see Figure 3.3).

It is not always possible to continuously define a unique normal vector \(\vec {n}(\mathbf {p})\) on every point of a regular surface, see, for example, the Möbius strip (see Figure 3.4).

When a normal vector and a closed curve that traces over the Möbius strip is chosen, and the normal is moved along this curve it ends up in the same point, however, the normal now points in the opposite direction. This observation can also be described with the help of parametrizations. On each point of the Möbius strip there exist different parametrizations of the surface, these parametrizations describe the same tangent plane. However, the normal vectors defined by the parametrizations point in opposite directions (see Figure 3.4) and it is not possible to continuously change from one parametrization to the other.

For the following discussions it is essential that the discussed surfaces allow to continuously define a unique normal vector in each surface point. To that end a parametrization of a neighborhood of each point on a regular surface is fixed, the set of all these parametrizations is called a family of coordinate neighborhood. If this family of coordinate neighborhood describes a positive movement around each point of the surface it is called an orientation of the surface which is formally defined as follows:

Intuitively, Definition 3.1.8 expresses that on an orientable surface the change in the direction of the normal vector is continuous when switching between two neighboring parametrizations.

If a regular surface \(S\) is defined by a differentiable function \(f:U \subset \mathbb {R}^3 \rightarrow \mathbb {R}\) as follows \(S = \{ (x, y, z) \in \mathbb {R}^3: f(x, y, z) = a\}\) and \(f^{-1}(a) \neq 0\) then \(S\) is orientable.

It should be mentioned that an orientation can locally be defined on every regular surface, however, orientability is a global property that concerns the entire surface.

To gain further insights into the geometry of a surface the orientation (i.e., \(\mathcal {N}\)) defined on the entire surface is further investigated.

Figure 3.5 visualizes how the Gauss map operates on the hyperbolic paraboloid. Another easily visualized example of the Gauss map is the image of the torus which simply is the entire surface of the unit sphere.

The derivative of the Gauss map \(d\mathcal {N}_{\mathbf {p}}\) in a point \(\mathbf {p}\) maps the tangent plane \(T_{\mathbf {p}}(S)\) into itself. Now consider a curve \(\gamma (t)\) in \(S\) with \(\gamma (0)=\mathbf {p}\), when the normal vector of the surface \(\vec {n}_S (\mathbf {p})\) is restricted to the curve \(\gamma (t)\) the map \(d\mathcal {N}_{\mathbf {p}}(\gamma (0)) = \vec {n}_S' (\mathbf {p})\) results in a tangent vector in \(T_{\mathbf {p}}(S)\). So \(d\mathcal {N}_{\mathbf {p}}\) measures how fast \(\vec {n}_S (\mathbf {p})\) pulls away from the surface in the curve \(\gamma (0)\). This is a similar construction to the one discussed for regular parametrized curves with the difference that the curvature is measured by a linear map and not by the norm of a vector. The construction is independent of the chosen regular parametrized curve, so the following definition can be given:

This definition results in a curvature value dependent on the chosen curve on the surface, which leads to the question if these curves can be classified in some way. To get further insights into this question the Gauss map \(\mathcal {N}\) of a surface \(S\) is further investigated. The negative differential \(-d\mathcal {N}\) of the Gauss map (sometimes referred to as the Weingarten map or shape operator) is a self-adjoint linear operator (\(II\)). Thus, its Eigenvectors are orthogonal, and the Eigenvalues are real [25]. The matrix \(II\) is also known as the second fundamental form of the surface \(S\). Recall, that \(-d\mathcal {N}_{\mathbf {p}}\) measures the normal curvature for each regular parametrized curve passing through \(\mathbf {p}\). Since the map \(-d\mathcal {N}_{\mathbf {p}}\) is a self-adjoint linear operator the maximum and minimum normal curvatures in each point can be calculated, which leads to the following definition:

The corresponding Eigenvectors to the principal curvatures are called the principal directions of \(S\) in \(\mathbf {p}\). Since the map \(-d\mathcal {N}\) is a linear map the trace and the determinant of the map can be calculated.

The mean and Gaussian curvature are used to define two special classes of surfaces. Surfaces with a constant Gaussian curvature of \(0\) are called developable surfaces, these surfaces can be "rolled" into the plane, or in other words folded out of paper. A cylinder is an example of a developable surface and the sphere is an example of a non-developable surface.

Surfaces with a constant mean curvature of \(0\) are called minimal surfaces, if these surfaces are bounded they minimize the Dirichlet energy in each point [65]. The plane is a trivial example of a minimal surface, where an example for a non-trivial minimal surface is the helicoid. Minimal surfaces play an important role in the remainder of this work.

Until this point the curvatures of a surface have been derived from the definition of the Gauss map without any explicit reference to a local coordinate system (i.e., a parametrization). For the calculation of the surface curvatures of discrete surfaces a description of the Gauss map in a local coordinate system \(\{\mathbf {x}_u, \mathbf {x}_v\}\) is required and is given by [25]

\begin{align} \label {eq:WeingartenMapCoords} d\mathcal {N} = \begin{pmatrix} a_{11} & a_{12}\\ a_{21} & a_{22} \end {pmatrix} := -\underbrace {\begin{pmatrix} e & f\\ f & g \end {pmatrix}}_{II} \underbrace {\begin{pmatrix} E & F\\ F & G \end {pmatrix}^{-1}}_{I}, \end{align} with

\begin{align} \begin{array}{ccl} E = \mathbf {x}_u \cdot \mathbf {x}_u & \; & -f = \vec {n}_v \cdot \mathbf {x}_u = \vec {n}_u \cdot \mathbf {x}_v\\ F = \mathbf {x}_u \cdot \mathbf {x}_v & \; & -e = \vec {n}_u \cdot \mathbf {x}_u\\ G = \mathbf {x}_v \cdot \mathbf {x}_v & \; & -g = \vec {n}_v \cdot \mathbf {x}_v \end {array}, \end{align} \(I\) (see Equation 3.2) is also referred to as the first fundamental form. The mean and Gaussian curvature can thus be expressed in any basis as

\begin{align} \label {eq:GaussCurveDet} K &= \frac {\det (II)}{\det (I)},\\ \label {eq:MeanCurveDet} H &= \frac { \mathrm {trace} (I \: \mathrm {adj}(II))}{2\det (I)} \end{align} where \(\text {adj}(M)\) stands for the adjugate matrix.

Some methods of calculating the curvatures of discretized surfaces further down in this chapter only calculate the mean and the Gaussian curvature. The principal curvatures can be calculated from the mean and Gaussian curvature by solving the following equation

\begin{align} \kappa _{1,2} = H \pm \sqrt {H^2-K}. \end{align}

3.2 Curvature Calculation for Point Clouds

When estimating the surface curvatures of a point cloud similar obstacles, such as the approximation of the surface normal discussed in Section 2.1.1, have to be overcome. Since there is no information about the spatial coherence of the points in a point cloud, a suitable neighborhood has to be estimated first. After a suitable neighborhood has been identified, recall Proposition 3.1.7, which states that in a local neighborhood a surface can be expressed by a differentiable function \(f(x,y)\), which is also called the height function [66]. The goal now is to estimate \(f(x,y)\) with the points in the previously estimated neighborhood. \(f(x,y)\) can be expressed by a Taylor series of order \(n\)

\begin{align*} f(x,y) = J_{\vec {B},n}(x,y) + \mathcal {O}(\| (x,y) \|^{n+1}), \end{align*} \(J_{\vec {B},n}(x,y)\) is called a jet of order \(n\) or an \(n\)-jet. To obtain an \(n\)-jet of the height function defined by the points in a neighborhood around a point \(\mathbf {p}\) it has to be approximated. The \(n\)-jet is approximated in the least square sense which leads to the following optimization problem

\begin{align*} \min \left ( \sum _{i=1}^N (J_{\vec {B},n}(x_i,y_i) - f(x_i,y_i))^2 \right ). \end{align*} This optimization problem results in a vector \(\vec {B}\) that holds the first \(n\) coefficients of the Taylor polynomial associated with the \(n\)-jet.

\begin{align} \label {taylor_jet} J_{\vec {B},n} = B_{0,0} + B_{1,0}x +B_{0,1}y + \frac {1}{2}(B_{2,0}x^2 + 2 B_{1,1}xy + B_{0,2}y^2)+ \dots . \end{align} To calculate the curvatures of a surface requires at least a \(2\)-jet. The coefficients of the Taylor polynomial can then be used to calculate the Weingarten map (see Equation 3.2) with [66]

\begin{align} \begin{array}{ccl} E = 1 + B^2_{1,0} & \; & -f = \frac {2B_{0,2}}{\sqrt {1+B^2_{1,0}+B^2_{0,1}}}\\ F = B_{1,0}B_{0,1} & \; & -e = \frac {2B_{1,1}}{\sqrt {1+B^2_{1,0}+B^2_{0,1}}}\\ G = 1 + B^2_{0,1} & \; & -g = \frac {2B_{0,2}}{\sqrt {1+B^2_{1,0}+B^2_{0,1}}} \end {array}. \end{align} The here presented discussion on 3D surfaces can be similarly adapted for 2D surfaces [66].

3.3 Curvature Calculation for Surface Meshes

Estimating the curvatures of a 3D surface mesh requires less computational effort than for point clouds, since the neighborhood of a vertex and information about the geometry is already present in the discretization. The general idea for calculating the curvatures of a surface mesh is to associate an infinitesimal area around a vertex with the value of an operator, further referred to as spatial averaging. Two different operators are required to calculate the Gaussian and mean curvatures. When the area around a vertex converges towards zero the value of the quotient between these operators and the area converges towards the desired curvature value [67]. Thus, first an appropriate area around each vertex on a surface mesh has to be determined. An appropriate area around a given vertex is given by the so-called Voronoi region [67]. Furthermore, each Voronoi cell (i.e., the Voronoi area each triangle contributes to the Voronoi region) minimizes the error introduced by the spatial averaging and the discretization [68].

Before the calculation of the Voronoi region is presented another concept has to be introduced. The set of all vertices incident to a vertex \(\mathbf {x}_i\) is referred to as the 1-ring neighborhood of \(\mathbf {x}_i\) (\(N_1(\mathbf {x}_i)\)) (see Figure 3.6a).

Thus, the Voronoi region of a vertex \(\mathbf {x}_i\) in a triangle mesh is calculated as follows

\begin{align} \label {voronoi_region} A_{\text {Voronoi}}= \frac {1}{8}\sum _{j\in N_1(\mathbf {x}_i)} (\cot \alpha _{ij} +\cot \beta _{ij}) \|\mathbf {x}_i - \mathbf {x}_j \|^2, \end{align} see Figure 3.6b. This definition only holds for triangle meshes with non-obtuse triangles, which can in general not be guaranteed. The operators discussed further down in this section also result in valid curvature values for 1-ring neighborhoods which consist only of obtuse triangles. Therefore, the calculation of the area has to be adapted (see Algorithm 1) [67]. The algorithm checks each triangle in the 1-ring neighborhood, if a triangle is obtuse a fraction of its area is added to \(A_{\text {mixed}}\), otherwise the area of the Voronoi cell is added. \(A_{\text {mixed}}\) denotes the surface area of the mixed region around a vertex, the set of points making up the mixed region is referred to as \(A_M\).

Mean Curvature

As previously hinted, surfaces with a constant mean curvature of \(0\) describe a special class of surfaces: minimal surfaces. Minimal surfaces minimize the surface area of a surface, which lead to an ample body of literature [69, 70]. One of the investigated operators is the mean curvature normal operator which establishes a direct relation between the mean curvature flow and the minimization of the surface area [71]

\begin{align} \label {mean_curve_average} \lim _{\text {diam}(A_{\mathbf {p}})\rightarrow 0} \frac {\nabla A_{\mathbf {p}}}{A_{\mathbf {p}}} = 2H_{\mathbf {p}}\vec {n}_{\mathbf {p}}, \end{align} where \(\nabla A_{\mathbf {p}}\) stands for the gradient, and \(A_{\mathbf {p}}\) for the area around the point \(\mathbf {p}\). The mean curvature normal operator, also known as the Laplace-Beltrami operator, is defined as follows:

The integral of the Laplace-Beltrami operator over the area \(A_{\text {mixed}}\) in a vertex \(\mathbf {x}_i\) on a surface mesh can be discretized as follows [67]

\begin{align*} \iint \limits _{A_M} \textbf {H}(\mathbf {x}_i) dA = \frac {1}{2} \sum _{j \in N_1(\mathbf {x}_i)} (\cot \alpha _{ij} +\cot \beta _{ij}) (\mathbf {x}_i - \mathbf {x}_j). \end{align*} After inserting the discretized value of the Laplace-Beltrami operator and \(A_{\text {mixed}}\) into Equation 3.10. The mean curvature of a vertex \(\mathbf {x}_i\) on a surface mesh can be expressed as

\begin{align} H_{\mathbf {x}_i} = \frac {\|\sum _{j \in N_1(\mathbf {x}_i)} (\cot \alpha _{ij} +\cot \beta _{ij}) (\mathbf {x}_i - \mathbf {x}_j)\|}{4A_{\text {mixed}}}. \end{align}

Gaussian Curvature

The Gaussian curvature, like the mean curvature, can also be expressed as a fraction of an infinitesimal area around a point as follows [72]

\begin{align*} \lim _{\text {diam}(A_{\mathbf {p}})\rightarrow 0} \frac {A_{\mathbf {p}}^G}{A_{\mathbf {p}}} = K. \end{align*} To elaborate the meaning of \(A_{\mathbf {p}}^G\), recall the definition of the Gauss map \(\mathcal {N}\), it maps all normal vectors from a surface into the unit sphere. Consider a small surface area around a point \(\mathbf {p}\) on an orientable surface, then the image of that area under \(\mathcal {N}\) describes a small area on the surface of the unit sphere (e.g., \(A_{\mathbf {p}}^G\)). Instead of directly calculating the integral over the image of the Gauss map the Gauss-Bonnet theorem is used [67]

\begin{align*} \iint \limits _{A_M}K \, dA + \sum _{i=0}^k \varepsilon _i = 2 \pi . \end{align*} This is a simplified version of the theorem since the surface area of the sphere is investigated which has a geodesic curvature of \(0\) [25]. The \(\varepsilon _i\) stand for the external angles of the boundary of the Voronoi region. For Voronoi regions it is easy to see that the angle at the vertex \(\mathbf {x}_i\) is \(\theta _i = \varepsilon _i\), because both edges of the Voronoi cell are perpendicular to the edges of the triangle and, therefore, \(\theta _i + \alpha _i = \pi \) (see Figure 3.7). Thus, the Gaussian curvature of a vertex \(\mathbf {x}_i\) on a surface mesh can be expressed as

\begin{align} K_{\mathbf {x}_i} = \frac {\left (2\pi - \sum _{j=1}^{f} \theta _j\right )}{A_{\text {mixed}}}, \end{align} where \(f\) stands for the number of all triangles in \(N(\mathbf {x}_i)\).

3.4 Curvature Calculation for Implicit Surfaces

The computationally least expensive way to estimate the surface curvatures of a discrete surface considered in this work is the estimation for implicit surfaces, since, derivatives can be approximated directly from the discretization.

In this section three different approaches of calculating the surface curvatures of an implicitly defined surface are presented: The general formula for implicit surfaces, mean curvature estimation by calculating the variation of the normal vector, and mean curvature estimation through the shape operator. These three approaches are further analyzed in Chapter 5.

3.4.1 General Formula

The goal is to express Equations 3.4 and 3.5 with the derivatives of an implicit function \(\phi \). This is achieved in the following way, the determinants and traces are transformed such that they are expressed with respect to a local parametrized \(\mathbf {x}\) and the normal vector \(\vec {n}\) defined by the chosen parametrization [73]

\begin{align} K &= \frac {(\mathbf {x}_u \times \mathbf {x}_v) \cdot (\vec {n}_u \times \vec {n}_v)}{ \| \mathbf {x}_u \times \mathbf {x}_v \|^2}, \\ H &= \frac {(\mathbf {x}_u \times \mathbf {x}_v) \cdot ((\mathbf {x}_v \times \vec {n}_u) - (\mathbf {x}_u \times \vec {n}_v))}{ 2 \| \mathbf {x}_u \times \mathbf {x}_v \|^2}. \end{align} Afterwards, the terms of the above equations are transformed such that they are expressed as the derivatives of the implicit function. This is achieved with the help of matrix identities, which transform the above expressions into [73]

\begin{align} \label {eq:GeneralFormulaContinousGauss} K &= \frac {\nabla \phi \: \text {adj}(\text {Hess}(\phi )) \:\nabla \phi ^T }{ \| \nabla \phi \|^4}, \\ \label {eq:GeneralFormulaContinousMean} H &= \frac {\nabla \phi \: \text {Hess}(\phi ) \: \nabla \phi ^T - \| \nabla \phi \|^2 \: \text {trace}(\text {Hess}(\phi )) } { 2 \| \nabla \phi \|^3}, \end{align} where \(\text {Hess}(\phi )\) stands for the Hessian matrix of \(\phi \). Furthermore, the mean curvature of an implicit surface can be expressed as

\begin{align} \label {eq:MeanCurveDivNormal} H = \nabla \cdot \frac {\nabla \phi }{\| \nabla \phi \|}, \end{align} which can be seen by expanding Equation 3.16 and Equation 3.17 and comparing the resulting terms.

Now that an expression for the curvatures of an implicit surface is available it has to be discussed how these formulas can be used on discretized implicit surfaces. Considering a discretized implicit surface as defined in Section 2.3 (e.g., a level-set function), the derivatives of a level-set function can be approximated by central differences as follows

\begin{align} D_x(\phi _{i,j,k}) &\approx \frac {\phi _{i+1,j,k} - \phi _{i-1,j,k}}{2\Delta x}, \label {eq:lowerOrderX}\\ D_{xx}(\phi _{i,j,k}) &\approx \frac {\phi _{i+1,j,k} - 2\phi _{i,j,k} + \phi _{i-1,j,k}}{\Delta x^2}, \label {eq:lowerOrderXX}\\ D_{xy}(\phi _{i,j,k}) &\approx \frac {\phi _{i+1,j+1,k} - \phi _{i-1,j+1,k} - \phi _{i+1,j-1,k } + \phi _{i-1,j-1,k}}{{4\Delta x \Delta y}}, \label {eq:lowerOrderXY} \end{align} where \(\phi _{i,j,k}\) stands for the \(\phi \)-value of the level-set function at the grid point \((i,j,k)\), if the subscript indices are omitted it is assumed that they are \((i,j,k)\). Using these approximations and an expanded form of Equation 3.16, the discretized mean curvature of a level-set function can be approximated as follows

\begin{align} \label {eq:GeneralFormula} H = \frac {\splitdfrac {D_x(\phi )^2(D_{yy}(\phi ) + D_{zz}(\phi ))} {\splitfrac {+D_y(\phi )^2(D_{xx}(\phi )+ D_{zz}(\phi ))} {\splitfrac {+ D_z(\phi )^2(D_{xx}(\phi ) + D_{yy}(\phi ))} {\splitfrac {- 2 D_x(\phi )D_y(\phi )D_{xy}(\phi )} {\splitfrac {-2 D_y(\phi )D_z(\phi )D_{yz}(\phi )} {-2 D_z(\phi )D_x(\phi )D_{xz}(\phi )}}}}}} {2\| \nabla \phi \|^3}. \end{align} This method of calculating the mean curvature will be referred to as the General Formula method. The Gaussian curvature can be calculated similarly by using an expanded form of Equation 3.15

\begin{align} \label {eq:gaussianCurvature} K = \frac {\splitdfrac {D_x(\phi )^2(D_{yy}(\phi )D_{zz}(\phi ) - D_{yz}(\phi )^2)} {\splitdfrac {+ D_y(\phi )^2(D_{xx}(\phi )D_{zz}(\phi ) - D_{xz}(\phi )^2)} {\splitdfrac {+ D_z(\phi )^2(D_{xx}(\phi )D_{yy}(\phi ) - D_{xy}(\phi )^2)} {\splitdfrac {+2D_x(\phi )D_y(\phi )(D_{xz}(\phi )D_{yz}(\phi ) - D_{xy}(\phi )D_{zz}(\phi ))} {\splitdfrac {+2D_x(\phi )D_z(\phi )(D_{xy}(\phi )D_{yz}(\phi ) - D_{xz}(\phi )D_{yy}(\phi ))} {+2D_y(\phi )D_z(\phi )(D_{xy}(\phi )D_{xz}(\phi ) - D_{yz}(\phi )D_{xx}(\phi ))}}}}} } {(D_x(\phi )^2 + D_y(\phi )^2 + D_z(\phi )^2)^2}. \end{align} Calculating the Gaussian curvature requires the same derivatives as calculating the mean curvature.

The finite difference approximations used for calculating the surface curvatures require a certain amount of grid points to be calculated. The set of all grid points required for the calculation of a set of finite differences is referred to as a finite difference stencil, or just stencil (\(\eta (\mathbf {x})\)). A finite difference stencil is always considered with respect to the central grid point \(\mathbf {x} = (i,j,k)\). Two finite difference stencils are of particular importance for this work, the \(7\) point star stencil (\(\eta _S(\mathbf {x})\)), and the \(19\) point plane stencil (\(\eta _P(\mathbf {x})\)). Figure 3.8 shows an illustration of these stencils. Calculating the surface curvature with the General Formula requires a plane stencil.

Lastly it should be mentioned that the maximal curvature a level-set function can describe is bound by \(\pm 1/\Delta x\) [15], since the smallest circle radius which can be represented with a given resolution is \(\Delta x\).

3.4.2 Shape Operator

Recalling Equation 3.17, the formula states that the mean curvature of an implicit surface can be expressed as the gradient of the normal vector. This expression can be expanded and leads to the following expression [75]

\begin{align} \label {eq:CurveWithHess} \nabla \cdot \frac {\nabla \phi }{\| \nabla \phi \|} = \nabla \vec {n} = -\frac {1}{\| \nabla \phi \|} (I - \vec {n}\vec {n}^T) \text {Hess}(\phi ), \end{align} where \(I\) stands for the Identity matrix.

Assuming that the level-set function fulfills the Eikonal equation (see Equation 2.8) with \(F(\mathbf {x}) =1\), then the term \(1 / \| \phi \|\) can be ignored. The matrix \((I - \vec {n}\vec {n}^T)\) projects onto the tangent plane of the implicit surface. Restricting \(\nabla \vec {n}\) to the tangent plane shows that it is symmetric [75]. Thus, there exists an orthonormal basis \([\vec {p_1}, \vec {p_2}]\) in \(\mathbb {R}^2\) of the tangent plane, this basis can be expanded to an orthonormal basis in \(\mathbb {R}^3\) as follows \([\vec {p_1}, \vec {p_2}, \vec {n}]\). When the gradient of the normal vector is expressed in this basis it results in the following matrix

\begin{align} \nabla \vec {n} = \begin{pmatrix} \kappa _1 & 0 & \sigma _1 \\ 0 & \kappa _2 & \sigma _2 \\ 0 & 0 & 0 \end {pmatrix}. \end{align} Here, \(\kappa _1\) and \(\kappa _2\) stand for the principal curvatures with the principal directions \(\vec {p_1}\) and \(\vec {p_2}\), \(\sigma _1\) and \(\sigma _2\) stand for the so-called flow line curvatures [75]. So the mean curvature of the surface can be expressed as \(\text {trace}(\nabla \vec {n})/2\). Furthermore, the trace of a matrix is independent of the chosen basis, thus the mean curvature of a level-set function \(\phi \) with \(\| \nabla \phi \| = 1\) can be expressed as

\begin{align} \label {eq:shapeOperator} H = \frac {\text {trace}(\text {Hess}(\phi _{i,j,k}))}{2} = \frac {D_{xx}(\phi _{i,j,k}) + D_{yy}(\phi _{i,j,k}) + D_{zz}(\phi _{i,j,k})}{2}. \end{align} This method is referred to as the Shape Operator method. In Section 4.1.2 a method is discussed that describes how a level-set function can be constructed that fulfills \(\| \nabla \phi \| = 1\). Calculating the mean curvature of the zero level-set using the Shape Operator method requires a star stencil \(\eta _S\). It should be mentioned that the quality of the mean curvature calculated with this method on a discretized surface is highly dependent on the quality of the discretization, which is discussed further in Section 5.1.3.

3.4.3 Variation of Normal

The third method calculates the mean curvature of the level-set function by again using Equation 3.17 which is approximated through the Euler-Lagrange derivative [76, 77]. In addition to the first-order central differences, this approximation requires the first-order forward and backward differences which can be approximated by

\begin{align} \label {eq:FDForward} D^+_x(\phi _{i,j,k}) \approx \frac {\phi _{i+1,j,k} - \phi _{i,j,k}}{\Delta x}, \end{align}

\begin{align} \label {eq:FDBackward} D^-_x(\phi _{i,j,k}) \approx \frac {\phi _{i,j,k} - \phi _{i-1,j,k} }{\Delta x}. \end{align} In turn, the mean curvature can be approximated as [77]

\begin{multline} H = \nabla \cdot \frac {\nabla \phi }{\| \nabla \phi \|} \approx \frac {1}{2\Delta x}\bigg (\bigg ( \frac {D^+_x(\phi )}{\|\vec {g}^+_x \|} - \frac {D^-_x(\phi )}{\|\vec {g}^-_x \|} \bigg ) + \\ \bigg ( \frac {D^+_y(\phi )}{\|\vec {g}^+_y \|} - \frac {D^-_y(\phi )}{\|\vec {g}^-_y \|} \bigg ) + \bigg ( \frac {D^+_z(\phi )}{\|\vec {g}^+_z \|} - \frac {D^-_z(\phi )}{\|\vec {g}^-_z \|} \bigg )\bigg ), \end{multline} where \(\vec {g}^{\pm }_x\) is defined as

\begin{align} \begin{split} \vec {g}^{\pm }_x := (D^{\pm }_x(\phi _{i,j,k}), &\frac {1}{2}(D_y(\phi _{i\,j,k}) + D_y(\phi _{i\pm 1,j,k})),\\ &\frac {1}{2}(D_z(\phi _{i\pm 1,j,k}) + D_z(\phi _{i,j,k}))), \end {split} \end{align}

\begin{align} \begin{split} \vec {g}^{\pm }_y := (D^{\pm }_y(\phi _{i,j,k}), &\frac {1}{2}(D_x(\phi _{i,j\pm 1,k}) + D_x(\phi _{i,j,k})),\\ &\frac {1}{2}(D_z(\phi _{i,j\pm 1,k}) + D_z(\phi _{i,j,k}))), \end {split} \end{align}

\begin{align} \begin{split} \vec {g}^{\pm }_z := (D^{\pm }_z(\phi _{i,j,k})), &\frac {1}{2}(D_x(\phi _{i,j,k\pm 1}) + D_x(\phi _{i,j,k})),\\ &\frac {1}{2}(D_y(\phi _{i,j,k\pm 1}) + D_y(\phi _{i,j,k})). \end {split} \end{align} To calculate all the required finite differences a plane stencil \(\eta _P\) is required, this method is referred to as the Variation of Normal method. This method only calculates the mean curvature, thus, the Gaussian curvature has to be estimated separately by using Equation 3.22.

3.4.4 Curvature of 2D Implicit Surfaces (Curves)

An analogous construction to the one presented in Section 3.4.1 can be done for regular curves (i.e., 2D surfaces) [73]. Recall that regular curves only have one type of curvature which can be approximated for a 2D level-set function with the formula

\begin{equation} \label {curvature} \kappa = \frac {D_y(\phi )^2 D_{xx}(\phi ) - 2 D_x(\phi ) D_y(\phi ) D_{xy}(\phi ) + D_x(\phi )^2 D_{yy}(\phi )}{\| \nabla \phi \|^3}. \end{equation}

3.5 Intuitive Approach to Surface Classification

When considering the construction of the normal curvature presented in Section 3.1.2 one can see that the normal curvatures are defined through the curvature of a parametrized curve and the cosine of the angle between the surface normal and the normal of the considered parametrized curve. It is tempting to drop the formal mathematics behind curvature calculation and try to directly use the angle between the normal vectors of surface points in a neighborhood on a discretized surface for surface classification. However, this approach has several drawbacks which are discussed in the remainder of this section. The most striking drawback is, that this approach can only estimate angles with points of the discretization without any knowledge of the underlying parametrization of the surface.

Point Clouds

Consider a point cloud and a neighborhood \(M_i(\mathbf {x}_i)\) of a point \(\mathbf {x}_i\). The relative spatial position of the points in the neighborhood with respect to the investigated point \(\mathbf {x}_i\) is not known. So, it is possible that two points in the neighborhood describe the same curve on the surface and that one point is farther away from \(\mathbf {x}_i\). In this scenario it is not known which of the two points should be chosen for the classification of the surface which can lead to wrong classifications. Therefore, this approach is not suited for point clouds.

Surface Meshes

The triangles of a surface mesh are flat, thus, they have a mean and Gaussian curvature of 0. A similar argument can be used for edges, since the surface can only bend in one direction over an edge. Consequently, using the angle between normal vectors can of course only be used on vertices on a surface mesh, but there is no straightforward way to calculate the surface normal of a vertex. Which again, as with point clouds, makes this approach unsuited for surface meshes.

Level-Set Functions (Implicit Surfaces)

Finally, considering level-set function, some merit can be found in applying the normal vector angle approach. All points on a regular grid have the same distance from each other, so one could check if the points that lie in a box stencil (see Figure 3.9) are surface points (i.e., the level-set function changes sign from the point \(\mathbf {x}_i\) to this neighboring point). However, depending on the position of the zero level-set more points in the stencil are required, i.e., \(21\) in 2D and \(81\) in 3D, which is a lot more than calculating the derivatives for curvature calculation. Further drawbacks of this method on level-set functions are discussed in Section 5.1.4.