Improving Initial Flattening of Convex-Shaped Free-Form Mesh Surface Patches Using a Dynamic Virtual Boundary

This study proposes an efﬁcient algorithm for improving ﬂattening result of triangular mesh surface patches having a convex shape. The proposed approach, based on barycentric mapping technique, incorporates a dynamic virtual boundary, which considerably improves initial mapping result. The dynamic virtual boundary approach is utilized to reduce the distortions for the triangles near the boundary caused by the nature of convex combination technique. Mapping results of the proposed algorithm and the base technique are compared by area and shape accuracy metrics measured for several sample surfaces. The results prove the success of the proposed approach with respect to the base method.


INTRODUCTION
Finding planar counterparts of curved surfaces is a commonly known problem in the field of production (such as ships, toys, clothes, shoes and furniture industries). In the conventional method, during manufacturing, 2D pieces are sheared and assembled to form the final product. Ideally, strain factor should be taken into account in these shearing and assembling processes, because the strain produces an elastic energy that reduces quality and creates material fatigue in the final product. The conventional design process used in these industries is based on the trial-and-error method. The designer sketches 2D pieces on paper and makes a prototype of the desired product. If the result is unsatisfactory, the designer changes the patterns based on his experience and makes another prototype. This prototyping and

RELATED WORK AND MOTIVATION
The process of flattening a complex surface is the phenomenon of projecting a surface of 3D space onto a planar surface according to certain rules [8]. A lot of studies have been done so far, and when examined, the flattening of complex surfaces can be summarized under three basic flattening approaches: geometric flattening, mechanic flattening and the hybrid flattening, which is a combination of the first two [9]. The presented study deals mainly with geometric flattening techniques. In the study by Bodduluri and Ravani [10], freeform surfaces were first divided into simple sub-surfaces and then approached to developable surfaces. In that work, a new technique has been proposed for the geometric design of developable surfaces based on rational Bezier and B-spline curves.
In another study [11], a grid approximation flattening algorithm was proposed. In that study, an approximate local coarse projection relationship was established between a given parametric surface patch and its triangular area by minimizing the projection error function.
Ping presented a flattening algorithm [12] for complex surfaces in 1997. In the study, the metal sheet was divided into a series of strip areas and then approximated to a striped surface. Each striped surface was then divided into triangular meshes and each triangular element was mapped onto plane. From the results obtained, it is seen that the smoothed surface contains overlap and gaps. In the following year, along with some researchers, mainly by the same author, a study [13] was conducted to address the overlaps and gaps arising between stripes.
A deep research was done about interactive parametrization design in 2002 [14]. The novel family of intrinsic parametrizations of surface meshes was introduced in the study. The study produced least-distorted parametrizations by proposing a robust and fast technique. The reported results show that technique is effective in automatically designing optimal maps with or without boundary conditions.
In the early 2000s, two theoretical research findings were released by Floater. In the first one [15], he used Mean Value Theorem in the calculation of harmonic functions simplifying and improving methods for parametrization and morphing. The other study [16] was about using discrete form of Rado-Kneser-Choquet theorem for harmonic maps. Provided that the boundary of the triangulation is homeomorphic to a disk, then such mappings are one-to-one fulfilling a discrete maximum principle, stated in the study.
In another study, an innovative technique that alleviates the drawback of convex combination approach was presented [17]. The technique proposed using a virtual boundary serving to get rid of the high distortions produced by convex combination approach for the triangles near the boundary. The results reported in that study shows that the technique produced a parametrization of 3D triangular mesh surfaces with less distortion than the original approach.
In the work of Yang et al. [18], an optimization algorithm using the rational bilinear re-parametrization was proposed to improve the uniformity and orthogonality of iso-parametric curves for general rational Bézier surfaces. The reported results show that their algorithm generates more uniform and orthogonal iso-parametric curves across the rational Bézier surfaces.
In a study [19], Liu et al. proposed a fast mesh parametrization algorithm based on subdivision technique. Their method based on 4-point interpolatory subdivision approximation is indeed an extension of the chordal parameterization to surface case. Their algorithm reportedly avoids computation of linear system of equations, which makes it computationally efficient than previous methods. The proposed subdivision algorithm which is easy to implement is of high efficiency, especially in case of large number of vertices.
A simple and fast low-stretch mesh parametrization method was proposed in 2004 [20]. The study based on Floater's shape preserving parametrization and took its success one step further. As the parameter domain is a convex-shaped polygon, the optimization process in their proposed method does not generate flipped triangles as reported in that study. The efficiency and speed of the method was tested on several complex surfaces and compared with other techniques in literature.
In this work, an efficient technique has been developed to utilize the barycentric mapping theory and a dynamic virtual boundary technique in order to ensure that the length information of arbitrary 3D mesh surfaces is preserved as much as possible

340
computer systems science & engineering and to find the 2D planar equivalents. The parametrization of the 3D surface can be described as an operation to produce a function that will correspond to a 2D planar surface. This concept plays an important role in geometry processing, since it allows to open the way to convert 3D modeling problems into 2D space, and doing it here is relatively easy. Owing to the advances in range scanning and 3D printing technologies, effective parametrizations are important for large meshes since complex 3D surfaces are widely used everywhere. In this context, an efficient algorithm for improving flattening result of triangular mesh surface patches homeomorphic to a disk is presented in this study. The proposed method incorporates a dynamic virtual boundary to reduce the distortions for the triangles near the boundary caused by the nature of convex combination approach. The proposed algorithm produces good parameterization outputs with relatively low distortions, which can be presented as input to a further energy releasing process by an energy-based approach. The rest of this work is organized as follows: Section 3 gives the theory and concept of barycentric mapping. The proposed initial mapping approach is elaborated in Section 4. Metrics for measuring accuracy of planar developments are given in Section 5. An illustrative example about the proposed approach is given in Section 6. Experimental results and discussion are presented in Section 7. Finally, Section 8 concludes the paper.

BARYCENTRIC MAPPING THEORY
Barycentric mapping is a very efficient mapping method to produce parameterization of triangular surfaces [21]. This approach is based on Tutte's barycentric mapping theory [22] and also provides a basis for surface parameterization and approximation studies [23] of Floater, who is renowned for its successful works in the fields. This theorem states that for a triangular surface homeomorphic to a disk, if boundary nodes are on a convex polygon and internal nodes are a convex combination of their neighbors, then all (u, v) planar coordinates form a valid parameterization [21,23]. A scheme is needed for determining specification of convex combinations which guarantees a planar shape preserving local shapes of the 3D surface. Let S(X,F) be a triangulated surface defined by vertices (X={( and facets F (a matrix storing indices of triangles generated by Delaunay Triangulation). Considering the definitions here, the barycentric mapping method is composed of the following steps.
Step 1: Boundary node coordinates of the polygon on which the mapping will be done is specified. This polygon may be unit circle, square or an arbitrary enclosure resembling the original 3D surface. To summarize, are determined to be on the polygon (which is defined in D ⊂ R 2 ) in an anti-clockwise manner.
if i th and j th nodes are connected by an edge, Step 3: u 1 ,…,u n are defined as the solution of the following linear equations where Solving these equations yields P = P(F, ..,n; j =1,...,N . Here P is an embedding of S(X,F) in R 2 and u 1 , …, u N represent positions of nodes connected to each other by straight lines.
Tutte studied how to draw straight lines of planar graph [24] years ago. He suggested for large scale graphs including triangulated ones that Eq. 3 be run as λ i, j =1/d i for neighboring nodes (i = 1, . . . , n) where d i represents number of neighbors of i th node. Therefore, u i can be regarded as barycenter of neighbors. He has proven that there exists a unique solution, and P is a solution consisting of straight lines [23]. In other words, u 1 ,…,u n are distinct guaranteeing that no two edges intersect except at common end points of triangles [24]. Considering the general case in Eq. 2, the equations can be rewritten as follows in order to show that Eq. 3 has a unique solution [20,22].
The reason for arranging the equation as above is that boundary nodes are fixed on border line of the polygon and then internal nodes are obtained by solving Eqs. 2 and 3. When considering that u i has two separate components u and v, this equals two matrix equations: These two linear system of size n (number of internal nodes) are solved in order to find two column vectors storing positions of internal nodes: (u 1 , . . . , u n ) T and (v 1 , . . . , v n ) T . Right hand side of the equations, b 1 and b 2 store weighted coordinates of neighboring nodes. A is square matrix of size n × n and its elements are: That is to say, the matrix takes the following form.
The existence of unique solution of Eq. 3 is dependent on the non-singularity of the matrix A. In other words, this is obtained by computing the inverse of A. There are a few ways to solve Eq. 3. Sparse iterative and direct methods are the ones which are most effective for large-scale meshes. On the other hand, Gauss-Siedel solvers can be used for reasonable small meshes [21].

PROPOSED INITIAL MAPPING APPROACH
The proposed initial mapping method for computing planar equivalents of 3D surface patches is based on Floater's vol 34 no 6 November 2019 barycentric mapping theory, mathematical foundations of which is given in the previous chapter. The proposed approach effectively improves outputs of the base method by incorporating virtual boundary nodes computed dynamically depending upon the complex surface patch to be flattened. Fig. 1 illustrates virtual boundary nodes created dynamically in order to improve flattening result of 3D surface patch and its planar equivalent.
In the convex combination approach, fixing the real boundary nodes of a surface patch directly on a polygon causes the triangles near the boundary in 2D parametric space to be highly deformed and distorted [17] . So, a set of virtual nodes encompassing real nodes of the 3D surface patch are adaptively created. For the ease of computability, virtual nodes are placed on the trajectory of a circle, whose diameter is determined according to total triangular area of the 3D surface patch to be flattened. Radius of the circle that hosts dynamic virtual nodes is computed according to the formula: where A i is the area of the i th triangle in the facet model, n is the triangle number in the facet model. As the boundary polygon in the barycentric mapping method fixes the boundary nodes of planar equivalent of the 3D surface patch at the beginning, these boundary nodes remain constant throughout the planar development process. Thus, determination of the bounding polygon is critically important for the output of planar development to be of high fidelity to its original 3D counterpart.
Number of virtual boundary nodes to be placed on the trajectory of the circle affects the extent of similarity to the shape of original 3D surface patch. Therefore, number of virtual boundary nodes should be large enough to provide reasonable resolution on the trajectory. It is determined to be proportional to total node number of the 3D surface patch. How the determination of number of virtual boundary nodes affects the initial flattening result is explored for several resolution values later in the study. Boundary nodes are positioned on the trajectory of virtual circle of radius r v at equal angle intervals α v determined by the following equation: The proposed initial mapping approach with a dynamic virtual boundary is presented in Algorithm 1. Globally, the proposed algorithm takes two arguments X and F as input, and produces the planar equivalent P of S(X,F) as output. The output P consists of n nodes in the u-v coordinate system, i.e. planar coordinates.
In the barycentric mapping approach, neighborhood weights are computed by the simple expression λ i, j = 1/d i , which evenly distributes weighting shares inversely proportional to the number of neighbors. In the proposed study, these weights are determined by computing reciprocal of the distances between neighbors and then normalizing the result by a scaling factor. For this purpose, the expression λ i, j = 1/d i , used by Tutte for representation of graphs [19], is replaced by the two expressions given in steps 10 and 11 of Algorithm 1. The Euclidean distance between each neighboring vertex in the 3D triangular mesh is computed by the expression in the step 10. The expression in the step 11 normalizes the result dividing d i j by the total cost of neighbors of the each central vertex. Thus, it is aimed to preserve shape of the triangular mesh surface locally in the parametric space [25], which is substantiated by the test results in the following sections.

ACCURACY MEASUREMENT
Generally, two metrics are used in order to evaluate accuracy of the planar surface being developed. These metrics are area and shape accuracies, details of which are explained under the following subheadings.

Area Accuracy
During the development process, planar surface area may change, and the final relative area difference is computed by the formula [26]: where A is the real area of the original mesh surface patch before the development and A is the area of the corresponding surface patch on the planar surface after the development.
A relatively complex relation of the double integral must be solved in order to calculate the area of a surface defined by S(u, v) [26]. Since the facet model is adopted in order to represent a complex surface in the approach, it is not possible to determine the analytical form of the surface S in this study. Therefore, the area A can approximately be calculated simply by adding up area of each triangle in the facet model using the formula: where A i is the area of the ith triangle and n is the number of triangles.

Shape Accuracy
Shape accuracy represents the total difference rate between the length of arc on the surface mesh and its corresponding edge length [26]. The final relative shape accuracy is computed by the formula: where L is the real length of an arc on the original mesh surface patch before the development and L is the length of the corresponding edge on the plane after the development.
Assuming that the curve segment L is on the curve defined by the formula C = S (u(t), v(t)) on the surface S, computing the segment length requires to solve a complex double integral problem. Similar to the aforementioned problem related to the adopted model, it is not possible to determine the analytic form of the curve C. Therefore, the length L can approximately be calculated by simply adding up lengths of each triangle in the facet model using the formula: vol 34 no 6 November 2019 Algorithm 1 Initial Map With Virtual Boundary Input: A given 3D triangular mesh surface S(X,F) where X contains 3D points (nodes) and F contains vertex indices of triangles. Output: The initial 2D shape P of the surface S. 1: Compute the radius of the virtual boundary circle r v according to the Eq. (8). 2: Determine n v , the number of virtual boundary nodes according to the total number of original surface 3: Compute the step angle of the virtual nodes α v according to the Eq. (9). 4: Append computed virtual nodes to the original mesh surface S(X,F) to obtain the augmented surface S*(X,F). 5: Examine the 3D mesh surface S*(X,F) and determine boundary nodes (n+1,…,N). 6: Sort X(size of which is N x 3)in an order so that the boundary nodes follow internal nodes. 7: Choose u n+1 , …,u N to be K -sided convex polygon D ⊂R 2 in an anticlockwise sequence, corners of which are placed at angular intervals α v on the virtual circle. 8: for i from 1 to n do 9: for j from 1 to N do 10: where L i is the length of the ith triangle and m is the number of triangles.

NUMERICAL EXAMPLE
The proposed initial mapping algorithm has been run on a simple surface patch of 7-vertices (named Surf1) in order for the intermediate results to be tracked easily. For this purpose, the simple surface patch given in Fig. 3 is presented to the algorithm. The cloud of points, defined by X, corresponding to the vertices of the mesh surface patch Surf1 and the topology data F defining its triangular facets are as follows: The number of virtual boundary nodes is arbitrarily chosen to be n v = 8 for this sample surface patch. Augmented surface patch S*, which is obtained by incorporating a virtual boundary polygon is given in Fig. 4.
The augmented surface with virtual boundary vertices, S* represented by X* and F* is given by Herein, the new vertices given in green color in the matrix X*, corresponding to the points in 3D space are the vertices augmented to the original mesh surface patch by using a virtual boundary polygon. Translation offset values added to positions of original vertices are calculated to be rx=ry=r v /2=0.48 for this sample surface patch. When the augmented surface is presented as the input to Algorithm 1, the matrix Λ Λ Λ is computed to be as the following: Size of matrix Λ Λ Λ is determined both by the numbers of original and augmented mesh surface patches. For this sample surface, its size becomes 7×15. First seven vertices, coordinates of which are given in black color, of the matrix given in Eq. 15 correspond to internal nodes. On the other hand, the remaining eight vertices, coordinates of which are given in green color are the virtual boundary nodes computed dynamically for the sample surface. Then, the square matrix A is formed using the elements of the matrix Λ Λ Λ as in the way defined in Eq. 7. Flattening result (in u-v parameter space) for the subject augmented surface is given in Fig. 5a.
Ultimate flattening result of the original mesh surface patch Surf1, which is given in Fig. 5b is obtained by removing the virtual nodes from the flattening result of the augmented surface. At the end of initial planar development process, area and shape accuracies for the subject mesh surface are obtained as E S = 0.7570 and E C = 0.5825 under the chosen conditions.

EXPERIMENTAL RESULTS AND DISCUSSION
The proposed initial flattening algorithm is tested on several arbitrary triangular mesh surface patches of various sizes and only three of them are included in this section. One of the triangular mesh surface patches is Surf1, given in the previous section and the other two are Surf2 and Surf3, seen in Figs. 6 and 7. The cloud of points, defined by X, corresponding to the vertices of the mesh surface patch and the topology data F defining its triangular facets for Surf2 and Surf3 are given in Eq.17 and 18 respectively.       1 27 27 31 9 9 10 1 9 31 16 16

Effect of Virtual Boundary Node Number on Flattening Accuracy
The three sample mesh surfaces mentioned above are given as input to Algorithm 1 and initial flattening results are obtained for several number of virtual boundary nodes. However, only one flattening result is given for each sample surface patch due to the page limitation. Flattening results of the sample surface patches for n v = 64 are given in Figs. 8, 9 and 10.
For all of the sample surface patches, virtual boundary polygons with 4, 6, 8, 12, 16, 32 and 64 nodes have been formed and their mappings to the plane have been realized. Area and shape accuracies measured for the obtained flattening results are tabulated in Table 1.
Area and shape accuracies for the sample surface patch Surf1 are also illustrated in the bar graph in Fig. 11. When the number of virtual boundary nodes is n v = 4, the area accuracy and the shape accuracy is calculated to be E S = 0.8875 and E C = 0.7098, respectively. The accuracies decreased gradually with the increase in the number of virtual boundary nodes and are measured as E S = 0.2038 and E C = 0.4202 for n v = 64. Area and shape accuracies for Surf2 are also illustrated in the bar graph in Fig. 12. When the number of virtual boundary nodes is n v = 4, the area accuracy and the shape accuracy is calculated to be E S = 0.9211 and E C = 0.6974, respectively. The accuracies decreased gradually with the increase in the number of virtual boundary nodes and are measured as E S = 0.4906 and E C = 0.2697 for n v = 64.
Area and shape accuracies for Surf3 are also illustrated in the bar graph in Fig. 13. When the number of virtual boundary nodes is n v = 4, the area accuracy and the shape accuracy is calculated to be E S = 0.8400 and E C = 0.6521, respectively. The accuracies decreased gradually with the increase in the number of virtual boundary nodes and are measured as E S = 0.3351 and E C = 0.2905 for n v = 64. The test results reported here indicates that increasing the number of virtual boundary nodes (based on the size of the surface patches) by a certain extent improves area and shape accuracies of the planar projection obtained in the flattening result.

Effect of Neighborhood Weights on Flattening Accuracy
In the section describing the principle of the barycentric mapping theory, the method of calculating the neighborhood weight values λ i, j for each node has been described in detail. For each node in the basic method proposed by Tutte in the drawing of graph and later adopted by Floater, this value is simply calculated to be inversely proportional to the number of neighbors. However in the proposed method, it is calculated using the inverse of the Euclidean distance between neighboring vertices of a node (with the changes given in the steps 10 and 11 in Algorithm 1). In this way, the local geometry shape in the 3D space is transferred in an efficient way to the parametric space where the planar development is realized. For this purpose, the basic method has been executed for all of the sample surface patches under the same conditions as the proposed method, and the results are given in Tables 2-4. The area and shape accuracies given in Tables 2-4 for the three sample surface patches are also presented visually in the  Figs. 14-16. The number of virtual boundary nodes is increased from 4 to 64 at certain intervals to provide the virtual boundary polygon in various resolutions. In each case, the initial projection results of the basic method and the proposed method for the surface patch appear in the bar graph in terms of area and shape accuracy. It is seen from both Fig. 14 and Table 2 that the difference in area accuracy starts from 0.8903 and decreases to 0.8247 in the basic method, while in the proposed method it starts from 0.8875 and decreases to 0.2038 for this sample surface patch. Similarly, the difference in shape accuracy starts from 0.6355 in the basic method and decreases to 0.5676, while in the proposed method it starts from 0.7098 and decreases to 0.4202.
It is seen from both Fig. 15 and Table 3 that the difference in area accuracy starts from 0.9164 and decreases to 0.5825in the basic method, while in the proposed method it starts from 0.9211 and decreases to 0.4906 for Surf2 surface patch. Similarly, the difference in shape accuracy starts from 0.6632 in the basic method and decreases to 0.3437, while in the proposed method it starts from 0.6974 and decreases to 0.2697.
It is seen from both Fig. 16 and Table 4 that the difference in area accuracy starts from 0.8535 and decreases to 0.6623 in the basic method, while in the proposed method it starts from 0.8400 and decreases to 0.3351 for Surf2 surface patch. Similarly, the difference in shape accuracy starts from 0.6141 in the basic method and decreases to 0.4253, while in the proposed method it starts from 0.6521 and decreases to 0.2905.
The coordinates of the nodes on the u − v parametric space are plotted in Fig. 17, with the projection of the sample surface for n v = 64 by the basic method and the proposed method. All these results show that the proposed method is better than the basic method (for the same iteration number in finding the inverse matrix by iterative methods) regarding flattening accuracies. The proposed method produces planar equivalents with relatively high fidelity to their original surfaces when compared to those of the basic method.

CONCLUSION
This paper proposes an effective approach for flattening convexshaped mesh surface patches. The proposed method exploits barycentric mapping theory and a dynamic virtual boundary approach to improve initial flattening result. Barycentric mapping, also known as convex-combination approach is a base method suitable for realizing mappings of convex-shaped mesh surfaces. A dynamic virtual boundary is efficaciously incorporated so as to reduce the distortions of the triangles near the real boundary of the surface. In the proposed approach, accuracy of base barycentric mapping method is improved by using the inverse of the Euclidean distance between neighboring vertices of a node, which is normalized by the total cost of its neighbors. In this way, the local geometry of the surface in 3D the space is transferred efficiently to the parametric space where the planar development is realized. The experimental results prove that the proposed approach effectively improves flattening results regarding area and shape accuracy. In a future work, initial mapping results obtained may further be enhanced by incorporating an energy-based flattening approach to release strain energy inherent in them.