Abstract
This paper proposes a novel method for parametrisation and remeshing incomplete and irregular polygonal meshes. Spherical harmonics basis functions are used for parametrisation. This involves least squares fitting of spherical harmonics basis functions to the surface mesh. Tikhonov regularisation is then used to improve the parametrisation before remeshing the surface. Experiments show that the proposed techniques are effective for parametrising and remeshing polygonal meshes.
1. Introduction
Polygonal meshes are often used to represent the surfaces of graphical models. Often a mesh is irregular or unsatisfactory due to over- or undersampling or the triangularisation technique used. There is thus the need to improve the quality of the mesh by remeshing.
There are two main categories of remeshing techniques: parametrisation methods [1–5] and mesh adaptation strategies [6–8]. Parametrisation methods indirectly improve a given mesh by parametrising it, often through bijective mapping from the original mesh to a 2D planar domain. New sampling can then be performed on this planar domain and, using the inverse of the original mapping function, is mapped back to the surface to form a new mesh. Unlike parametrisation techniques, mesh adaption strategies directly improve meshes using relaxation or modification processes to iteratively make local improvements to the mesh [9], until either some size or quality criterion is reached.
Parametrisation and adaption approaches can also be combined. An example of this is a method that assumes a mesh is formed of parametric patches bounded by parametric curves [10] and iteratively improves the mesh by calculating the quality of a patch. Vertices and therefore faces are inserted into some areas and deleted from others until the global quality criterion is achieved. The refinement of a mesh region is determined by a local error computed at each vertex in a patch, which is the difference between the surface curvatures of the original patch and the new patch.
Another example of the combined approaches is the isometric parametrisations (MIPs) method [11]. Here, the remeshing is an adaption strategy on a mesh parametrised to a planar domain and the final result is mapped back to the domain of the original mesh. The remeshing process involves an umbrella operator that selects a vertex to be improved while fixing all neighbouring vertices and treats edges connecting the centre vertex and the neighbouring vertices as “springs.” Each spring is assigned an energy, and the minimisation of the combined energy of all springs is used to determine the new position of the centre vertex. Essentially, this minimisation can be seen as a push and pull to reach an equilibrium of energy. The process is repeated for all vertices until the mesh reaches a certain specified quality or the difference in movement in the vertex positions between iterations is below a certain threshold.
Such local improvement methods tend to be greedy in their approach, and therefore the global quality of the resulting mesh can suffer [12]. These meshes may also be inefficient when carrying out numerous sampling operations in 3D space [3]. Furthermore, methods that parametrise surfaces to a single domain are often constrained to meshes that are topologically similar to that domain.
One group of methods exploits the fact that some meshes can be transformed into a single topologically planar mesh through a series of cuts [13]. These cuts are stored in a graph and are applied to the original mesh prior to parametrisation. The cuts allow the mesh to be simplified for parametrisation using a planar domain. Creating the cut graph for these meshes is the major challenge of such methods. Meshes of genus-1 or above will require multiple cuts and the particular location of these can be difficult to compute automatically. The particular location of each cut needs to be considered carefully to strike a balance between the goal of low distortion and discontinuities introduced with each cut. Even though methods such as these overcome the topological limitations of mesh parametrisation, they often suffer from visible discontinuities along the seams of the cuts. Where the discontinuities are not tolerated by some applications, changing the base domain is often worthwhile [12].
Another approach uses harmonic mappings to parametrise individual segments, before each segment is remeshed [14]. Although this method can produce high quality meshes, with very little to no discontinuities, and is computationally robust, the partitioning of the mesh can result in unnecessary segmentation. For example, it may segment a straight tube into two segments [14]. This is due to the planar domain used, limiting the topology of the mesh. There is a considerable advantage to using spherical domain for parametrisation, as many meshes can be parametrised to a sphere without the need to cut the mesh [15]. Planar meshes that can be parametrised by a square domain can be mapped to a spherical domain using the process of UV mapping [12].
Spherical parametrisation is not without its problems. Many methods face the issue of not being able to guarantee a bijective (one-to-one) mapping between the sphere and the mesh. Some attempts to overcome this create a geometry image that approximates the mesh [16]. First, the mesh is mapped to the spherical domain, and then an arbitrary polyhedron is spherically parametrised. The 2D geometry image is then created by unfolding the polyhedron and then remeshing it by mapping back to the original surface domain. However, the method demonstrates an inherent limitation of spherical parametrisation. A trade-off is required between stretching of the mesh and conformality, as both cannot be attained simultaneously for highly deformed shapes.
This paper proposes a technique for parametrising 3D meshes for remeshing which uses the theory of spherical harmonics to approximate a continuous surface. Spherical harmonics are a natural basis for representing functions defined over spherical and hemispherical domains. They have been used for many applications in 3D modelling, including face recognition [17], lighting and systems [18, 19], and diffusion imaging [20–22]. Spherical harmonics have been suggested as a parametrisation technique for surfaces [13]. To use spherical harmonics basis, a solution to find the weights of the basis functions is obtained through solving a linear system of equations. To overcome computational challenges in using linear least squares to fit a large mesh using a large number of basis functions, regularisation is investigated for reducing the effects of numerical outliers and computational inefficiency. Once the mesh has been parametrised as a combination of spherical harmonic basis functions, the new mesh defined on a sphere can be remeshed to approximate the original surface. Further to this, the consequences of surface parametrisation as a solution to surface inpainting are addressed.
2. Methods
2.1. Spherical Harmonics
For a sufficiently smooth surface, represented by a function , an infinite series of spherical harmonic basis functions can be used to represent it in the following form [23]:where and are the polar and azimuth angles in a spherical coordinate system. As , this representation becomes an exact description of the surface . The spherical harmonic functions are defined by , with order and degree (), and is some weighting coefficient for . The degree describes the number of basis functions to be computed for each order.
Spherical harmonics are formed from the set of solutions to the 3D Laplace equation, given here as a combination of associated Legendre polynomials, :The associated Legendre polynomials, , are defined as follows:
The order specifies the number of polynomial terms that a harmonic function contains. Relatively smooth functions can be obtained using low order spherical harmonics, with increasing order corresponding to an increased number of frequencies captured by the spherical harmonics. The higher the order, the more precise the approximation. Figure 1 shows the spherical harmonics for even orders up to —note that is the reflection of .
The expression of spherical harmonics uses complex domain functionals in (2). However, for ease of programming, and since many of the desired properties are still present, only the real part of is used, denoted by . This is calculated as follows:whereThe function in (1) may be solved for to calculate the weighting of each basis function and thereby allow an analytical representation of the surface, using the inverse transform known as the spherical harmonics transform:However, this is often not possible to solve. As an alternative, (1) can be written as a series of linear equations in matrix form, for some finite on some discrete sampling of the surface, , for . The system is then solved for coefficients , :where , , and . The unique indexing for each pair is as used in [13].
2.2. Regularisation
The spherical harmonics basis is used to approximate the surface using least squares optimisation. The procedure minimises the sum of the squared residuals to the curve determined by the best fit. This involves solving a set of linear equations, where, for an overdetermined system of equations, , with unknown, the best approximation, is calculated such that the sum of square differences between and is minimised. Often, the solution does not exist, and instead the best approximation of is found, minimising the function . This can be done more simply by solving the equivalent system of equations . For the spherical harmonics system, in (7), the system to approximate the weighting coefficients .
Regularisation is often used to optimise solutions. This is done by introducing additional information, usually in the form of a penalty for complexity, for example, restrictions to ensure smoothness. Tikhonov regularisation is one of the most commonly used methods for regularisation [24, 25]. The additional information introduced is described as the Tikhonov matrix, which is added to the original minimisation term, acting as a smoothness constraint. The Tikhonov matrix, , can be constructed in a number of ways. One such way is to use unimodal regularisation () [26]. Alternatively, the construction may be based on the Laplace-Beltrami smoothing operator [27–29].
Due to the linear least squares allowing outlying points to have a disproportionate effect on the fitting, there is often a degree of noise that presents itself when fitting the spherical harmonics to the surface using higher orders of . The noise is amplified as the order increases, due to the massive increase in the number of calculations performed. As the matrix system gets bigger, a higher number of numerical inaccuracies occur, increasing the number of outlying points generated. Consequently, constructing the Tikhonov matrix based on the Laplace-Beltrami operator is the natural choice in this case, as the regularisation introduced will smooth the harmonic signal produced. The regularisation term is introduced in minimisation term as the following constraint:where is the basis matrix, is the radii vector, is the coefficient vector, is a tuning parameter used to balance the fidelity of the data and the smoothness requirement [30], and is the Tikhonov matrix. Minimising this newly formed functional results in an equivalent solution to the following equation: where and has the following structure:where is the spherical harmonics order associated with the th coefficient.
L2 regularisation is more computationally efficient than L1 regularisation as it has an analytical solution. It also gives nonsparse outputs. Sparsity refers to the number of values in a vector or matrix that are nonzero where a dense matrix or vector has mostly nonzero entries. Conversely, a sparse matrix or vector contains very few nonzero entries. Using an L2-norm penalty gives nonsparse coefficients which is desirable with spherical harmonics whereas using an L1-norm penalty would result in sparse outputs. Therefore, L2-norm penalty has been used instead of L1-norm. The multiplicity of eigenfunctions is beyond the scope of this paper and readers are invited to read [31].
2.3. Spherical Meshing
Typically surface representations as polygonal meshes are described in Cartesian coordinates (). For the surface mesh to be represented by , the mesh must be transformed into spherical polar coordinates () about the origin, where , the distance from the origin to the point on the surface at inclination and azimuth : On calculating the spherical harmonics coefficients and obtaining the approximation of continuous functional , a value for may be obtained for any given and . As a result, it is possible to define a new set of points, or a new mesh, on a unit sphere (or any sphere, disregarding its radius). An example of regularly sampled spherical meshes is shown in Figure 2, which can be used to remesh an approximated surface. The mapping from polar coordinates to Cartesian coordinates is , , and .
(a)
(b)
3. Implementation
3.1. Surface Parametrisation
Given a mesh, with vertices transformed into spherical coordinates treating the geometric centre of all mesh points as the origin, a set of basis functions for angles and is calculated up to the specified order . Each basis function is calculated according to (4) and the values indexed in -by- matrix, where is the number of vertices in the mesh and . As described previously, for the th vertex, the calculated values for each basis function are added to the matrix at , where .
The system of equations is set up as in (7), where , the spherical radius value at angle , and solved using the Laplace-Beltrami Tikhonov regularisation described in Section 2.2, using to control the added smoothness constraint. The resulting vector of coefficients can then be used to weight spherical harmonics in the functional (1), describing the surface which the mesh represents. As the spherical harmonics series has been truncated to a maximum order, , the computed values of (corresponding to radius values in a spherical representation) are only approximations of the original values.
Multiplying the calculated coefficients by the spherical harmonics values for the original set of mesh vertex angles, the residual error for each vertex can be obtained in the form , where are the bases, is the set of approximated coefficients, and is values of for each pair of and .
The newly reapproximated vertices can also be visualised by converting the values from spherical to Cartesian coordinates and retaining the face information of the original mesh. This face information is in the form of a list. This list does not contain the vertex values themselves but the indices for each vertex. Therefore, each line in the list contains three (for a triangular mesh) indices that comprise a single face. The face information remains valid as long as the order of initial vertices and their approximations stay the same.
3.2. Remeshing
Once the coefficients for the mesh have been calculated, they can be used in the remeshing process. New sampling is defined on angular sphere and the spherical angles of the vertices are used to calculate values for the basis functions. Note that these basis functions are still limited by used when generating the set of coefficients to ensure the matrix dimensions agree. A new basis matrix of dimension is formed, where is the number of vertices defined on the sphere. The new values of : are calculated by multiplying out the basis functions and coefficients. The angles and for the new spherical mesh, coupled with approximation , can be converted into a set of Cartesian coordinates describing the remeshed surface of the original mesh. Again, face information defined on the sphere can be applied to these new vertices.
4. Results
4.1. Parametrisation
It is important for a remeshing algorithm to approximate an input surface well. Therefore, there is a need to measure how close the reconstructed surface is to the original. The numerical measures used in [13] to describe the fit are the mean residual and the maximum residual, and consequently these measures are used here. The residual is calculated for each vertex in the original and approximate mesh, as described in Section 3.1, and for two given meshes, the surfaces are approximated and results given. These meshes, both representing a human head, exhibit different levels of features. They are distinguished in this paper as Mannequin and Obrubovka and shown in Figure 3.
(a)
(b)
The residuals of parametrisation of the meshes by spherical harmonics, solved only by linear least squares for different values of , are given in Table 1, with reconstructed meshes shown in Figure 4. As the order of the spherical harmonics increases, the residual errors decrease. This is as expected, since a higher order captures higher frequencies in the surface and is a better representation of the structure.
(a)
(b)
In addition to the quantitative metrics, the resulting meshes are plotted and shaded with the local vertex residuals, shown in Figure 5. Visual analysis shows the performance of spherical harmonics as a parametrisation technique is suitable for a large proportion of a mesh, with the highest errors occurring where the curvature in the mesh is sharpest, notably the nose and brow region. As the maximum order is increased, these effects are reduced, as demonstrated by the results presented in Table 1. The results here agree with the observation that spherical harmonics perform well as a parametrisation technique [13]. The reconstruction in Figure 4 is visually very similar to the original mesh and spherical harmonics representation can capture main surface features from the original mesh well. However, at lower orders, very sharp features are difficult to reconstruct, as can be seen in the ears of the Mannequin mesh at .
(a)
(b)
4.2. Remeshing
Each of the meshes parametrised has been remeshed using a unit icosphere (4-time subdivided regular icosahedron) at increasing orders: , 20, 40, and 60. The calculation of the basis coefficients has been done using linear least squares only and has not been subject to regularisation. There are two factors that influence the number of calculations that must now be performed: the value of and the number of vertices in the new mesh. Increasing the order, , and increasing the resolution of the mesh will significantly increase the number of calculations required to solve the linear system of equations. This increase in the number of calculations is due to the larger number of basis functions that would need to be calculated and utilised in the linear system.
Figure 6 shows the remeshing results for the two models. These show outlying regions affected by high frequency harmonics, which increase in noise as the order increases. At high levels, the original mesh’s structure is virtually unrecognisable in the remeshed figure. The “spikes” are amplified particularly in regions of low detail or resolution in the original mesh. The mesh created through the remeshing process inherits many properties of the given spherical sampling mesh, in this case a pseudouniform triangular grid over the entire surface. Increasing the resolution of the mesh much greater than that of the initial mesh can significantly reduce the interpolating effects of the remeshing process.
4.3. Tikhonov Regularisation
Tikhonov regularisation was introduced to aid solving the discrete spherical harmonics transform, adding a smoothing constraint to the resulting surface representation and reducing computational runtime in solving the linear system [32]. Residuals errors of the results for the original meshes parametrised using Tikhonov regularisation are shown in Table 2 for increased order, , and various values of , the weighting coefficient for the Tikhonov matrix. Note that when , the result is the linear least squares solution. The residuals for the results of varied combinations of and are shown in Figure 8, showing the lowest residuals (both mean and maximum) for high order () and low regularisation ().
Using the same process for remeshing, the regularised coefficients are used to approximate the meshes on finer, more uniform sampling. The effects of the smoothing reduce the impact of the high frequency harmonics, allowing high orders to increase refinement of the mesh, while avoiding numerical errors and “spiky” outliers. The resulting meshes using the icosphere sampling are shown in Figure 7 at different orders, showing the effectiveness of the remeshing technique and the substantial improvement regularisation brings about compared to the results using only linear least squares (Figure 6).
(a)
(b)
(a)
(b)
(c)
(d)
However, the regularisation process, by nature, smoothes the entire mesh surface resulting in some loss of sharp edges and folds that are present in the original mesh. Another factor that can result in the loss of sharp regions is the vertex placement used in the remeshing process. Where vertices are not specifically placed at edge regions, it can oversmooth the surface. As shown in both Table 2 and Figure 8, the smoothing effects of the regularisation parameter decrease the accuracy of the spherical harmonics approximation; however, a compromise must be reached between the smoothing constraint and the numerical anomalies experienced. The combination of high order and very low regularisation coefficient, , gives global minima in the residual plots for the combinations tested. Visually, the remeshed plots show clearly retained structure over the regular grid (Figure 7) without the presence of the “spikes” seen in the nonregularised linear least squares representation.
4.4. Surface Inpainting
When surface meshes are created, they may contain undersampled regions rendering the mesh incomplete with holes. For example, when using 3D range scanners to produce a mesh of an object, there may be holes on the surface of the mesh due to occlusion. In order to use the mesh for processes such as animation, texturing, and simulation, the holes on the surface need to be closed in a preprocessing step. Filling these holes is known as surface inpainting. Current methods for surface inpainting incorporate hole detection schemes that use the local context of the surface on the boundary of the hole in order to make a judgement on how to close the hole. Information gathered by the local context includes the surface curvature and sampling rate. Locating the hole on a mesh, especially a large highly sampled mesh, can be computationally expensive.
Using spherical harmonics, holes on the mesh can be closed without the need to locate the holes first. This inpainting process is a natural consequence of the spherical harmonics remeshing process. This consequence is the new mesh inheriting the properties of the sphere mesh used in the remeshing process and effectively interpolating across the hole. The closure of the hole occurs when the surface is remeshed but this is done after computing the continuous surface description. This means that the coefficients for the surface are calculated when there is missing information for the surface which may result in their incorrect determination. The larger the hole in the mesh, the larger the error in the coefficient values. A sufficiently large hole causes the newly generated mesh to “balloon” where the hole was originally located, similar to the spikes experienced in the nonregularised remeshing. Minor ballooning can be mitigated through increasing the spherical harmonics order and regularisation parameter. Increasing the order of the spherical harmonics allows for higher frequencies to be captured and therefore a closer approximation, of the surface, can be achieved, while regularisation will constrain these frequencies from having too much of an effect where the surface is unknown.
5. Conclusions
This paper proposes a technique for remeshing using spherical harmonics as a parametrisation technique. By representing a mesh through a finite series of spherical basis functions with increasing frequencies, it has been possible to represent the surface over an arbitrary computational grid. The spherical functional allows seamless parametrisation of nonoccluded genus-0 meshes. The results have demonstrated low residual errors between approximated vertices and their original positions, showing a clear positive correlation between accuracy and the order of spherical harmonics basis. The coarse-to-fine hierarchy of spherical harmonics allows for the resolution of the model to be altered with ease. This method could be used as a preprocessing algorithm for processes that require a regular mesh such as animation, texturing, simulation, and modelling.
The approach presented was tested with both quantitative and visual analysis of results. Spherical harmonics perform well as a parametrisation technique, although when the basis equations were solved using only linear least squares a source of error occurred in regions with little detail available in the original mesh. This was observed in the remeshing, with errors increasing with basis order, a result of using higher frequencies to represent the surface. As a result, Tikhonov regularisation as a smoothing measure was investigated and found to be effective at constraining high frequency errors without significantly reducing the approximated detail. The Tikhonov matrix used is based on the Laplace-Beltrami smoothing operator, so for highly weighted regularisation, oversmoothing may occur. The best results were observed with high spherical harmonics order and low (nonzero) levels of regularisation.
The computational complexity of this method is affected by two main variables: the number of vertices, , of the original mesh to be remeshed and the value of . The larger the values of and , the more the basis functions that need to be computed. With more basis functions, there are more computations required in solving the linear system in (9). The complexity of calculating the basis functions becomes .
The technique has some limitations. The first is a limitation regarding features in the original mesh. It is common for vertices to be placed along edges of a surface in order to capture sharpness of the edge. However, remeshing may not account for these features, and the surfaces are by nature rounded by spherical harmonics. Research efforts in the future may be focussed towards feature detection and feature preservation schemes, for example, searching for sharp edges in the mesh and preserving these during the remeshing by placing vertices along the original edges.
Additionally, there is compromise between regularisation parameter and spherical harmonics order, which are inversely proportional in producing the highest accuracy. Despite this, if the weight of the regularisation matrix in solving the spherical harmonics system is too low, the noisy spikes can occur during remeshing. This can be manually selected, or similar effects may be investigated with regularisation in alternative solvers, such as those discussed in [13], which deals with calculating approximations of high resolution meshes.
In its present form, the technique here does not explicitly deal with the presence of occluded regions in its mapping to spherical domain. This is especially important in higher genus surfaces or those with structural complexities, such as the ears of the Stanford bunny. Ongoing research is investigating surface segmentation to remesh individual components. The main challenge here is combining remeshed components in a seamless way.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.