Smooth Local Interpolation of Surfaces Using Normal Vectors

This paper proposes a simple surface interpolation attaining tangent-plane continuity. It is a natural extension of the local quadratic C0 interpolator developed by the author 2005 in one of his works, which has already been applied successfully to diverse engineering problems. The methodology presented in this paper inherits most of the advantages possessed by the C0 scheme. That is, i The algorithm is efficient and completely local requiring only the position vectors and normals given at the nodes of a patch, and hence it is suitable for parallel processing. ii It converges rapidly to the given surface with the increase in the number of nodes. iii Singular points apexes, sharp edges, etc. and nonmanifolds can be treated quite easily. iv Because of the minimization criteria assigned to the surface coefficients, it is rather robust and amenable to computational analyses. Validity and effectiveness of the proposed technique are demonstrated through numerical examples.


Introduction
There is a significant gap between the requirements on geometric models in the CAD and computational science communities.In the former, surface descriptions using few patches with large degree of freedom and high level of continuity e.g., NURBS, Bézier, Gregory are considered to be desirable.On the other hand, models for numerical simulation e.g., the finite element method are almost always represented by fine meshes, where the major interests are in fast, stable, and accurate recovery of the surface information lost during the process of discretization.With this as background, the author proposed an interpolation scheme suitable for such analyses 1 .It has already been applied successfully to engineering problems including A high-precision machining 2 , Other differences of the method from traditional approaches smooth local interpolators in particular are discussed extensively in 1 , and hence are not repeated here.
Another merit of the C 0 patch is that it can easily be extended to establish tangentplane continuity, as reported in this paper.The proposed interpolation is named G 1 Nagata patch here for convenience of reference.The term does not imply that the interpolation always generates G 1 surfaces; it can reproduce singular points apexes, sharp edges, etc. or may yield cusps similar to the C 0 patch.The new algorithm involves correction using a simple rational function retaining the boundary of the C 0 patch.Hence it is no longer quadratic, but it inherits all other desirable features from the C 0 patch, such as complete locality and capability of handling multiple normals.Similar to the C 0 patch, the present algorithm assigns minimization criteria inherent in the generalized inverse to the surface coefficients, and hence is rather robust.This is notable since accuracy and stability are not compatible with smoothness for most surface interpolators presently available.
This paper is organized as follows.The next section formulates the problem, briefly reviews the C 0 Nagata patch, identifies the conditions on its correction for the tangentplane continuity, and gives a simple solution satisfying the requirement, thus deriving the interpolation scheme in Section 2.2.5.Section 3 assesses the convergence and accuracy of the interpolation through numerical examples, followed by Section 4 summarizing the results with a scope for future work.

Formulation
Assume that position vectors and normals are given at the vertices of a triangular mesh with general topology.The objective here is to recover the curvature and smoothness of the surface through interpolating each patch independently.The process consists of the following three steps: A to replace each edge of a patch with a curve orthogonal to the normals given at the endpoints of the edge, B to determine a parametric quadratic polynomial patch C 0 Nagata patch reproducing the modified boundary to keep the C 0 continuity, C to correct the patch preserving the boundary in order to obtain tangent-plane continuity.
Algorithms for A and B were already proposed by the author in 1 .However, they are essential for the present approach and hence are summarized in the next section for completeness.The C 0 patch is then upgraded to enable local smooth interpolation.

The C 0 Patch
For the time being, assume a smoothsurface.Consider a curve segment on the surface as illustrated in Figure 1.Its endpoints P 0 and P 1 have the position vectors x 0 , x 1 and the unit normals n 0 , n 1 In this paper, a bold typeface is used for vectors column matrices , and bold symbols in square brackets denote matrices. .They are supposed to be known as input.
The step A mentioned at the beginning of Section 2 approximates the curve segment using the quadratic function to reproduce the position and normal vectors at the endpoints P 0 , P 1 .Here ξ ∈ 0, 1 denotes the parameter for the curve, and is the vector connecting the endpoints.The coefficient c in 2.1 introduces curvature to the segment, and hence is named the curvature parameter.It is determined from the boundary conditions at P 0 , P 1 as

2.3
Here 0 is the zero vector.For general surfaces which may have multiple normals at P 0 , P 1 , the curvature parameter can be extended as where n , {d} are constructed by 32 of 1 from the position and normal vectors at P 0 , P 1 .The generalized inverse of the former is explicitly written as where O is the zero matrix, and a k , D k are coefficients defined by 27b and 29 of 1 , respectively.Note that D 0 , D 1 , D 2 are dimensionless and should be treated as zeros when their absolute values are below a small positive threshold ε, for avoiding numerical instability.The interpolator 2.1 is at most quadratic, and hence its approximation capability is rather limited.If the magnitude of the curvature parameter c is too large, failure in the aforementioned singularity detection or insufficient resolution of the original data is suspected.This situation can be handled as follows.Note that the distance between P 0 and P 1 is |d| Figure 1 and the curvature parameter c also has the dimension of length.If the latter increases, the interpolated curve deviates from the line segment P 0 P 1 .Therefore, can be used as a criterion for inadequate interpolation.Here δ c is a dimensionless parameter and hence can be a constant regardless of the geometric scale.If the above holds, the next expression on the right-hand side of 2.5 is chosen for evaluating the generalized inverse irrespective of the decision based on the threshold ε, and the curvature parameter c of 2.4 is recalculated.If it still satisfies 2.6 , the same process is repeated.In the worst case, the generalized inverse 2.4 becomes the zero matrix resulting in c 0 making the curve segment linear.Since the case is obviously against 2.6 , this iteration always terminates and is safe.Similar treatment should be taken also for 2.3 .With the above algorithm as a basic tool, the triangular patch in Figure 2 a is interpolated.Its geometry can be described by where x denotes the position vector to a point on the patch, and η, ζ are its parameters defined within the domain of Figure 2 b .In Figure 2 a , the position vectors to the vertices v 1 , v 2 , v 3 : x 00 ≡ x 0, 0 , and the unit normals at the three points n 00 ≡ n 0, 0 , 2.9 are assumed to be known.The C 0 patch approximates the surface using the quadratic polynomial

2.10
The edges of the patch are first interpolated in the form of 2.1 as

2.12
Here the curvature parameter c is evaluated by 2.3 for smooth surfaces or by 2.4 for general cases.The parametric representation 2.10 of the patch is uniquely obtained from the boundary curves of 2.11 as

2.13
In the subsequent discussion, symmetric expressions play an essential role.For this purpose, barycentric coordinates are introduced to replace the surface parameters η, ζ.In the parameter space of Figure 2 b , the vertices of the triangular patch have the following position vectors:

2.14
Hence the parameter vector for a point on the patch can be described as where β 0 , β 1 , β 2 are the barycentric coordinates for η.They are in one-to-one correspondence with the surface parameters provided that

2.15b
Equations 2.15a -2.15b can be solved for the barycentric coordinates as

2.16
Equations 2.16 and 2.7 show that the following holds on the patch:

2.17
Let V i be the vertex for the parameter vector η i and E i its opposite edge.Equations 2.15a -2.15b give their representation in the barycentric coordinates as

2.18b
Here i c denotes the set of all cyclic permutations of the indices 0, 1, and 2; that is, Introducing the barycentric coordinate vector equations 2.18a -2.18b can be rewritten as where 1 i i 0, 1, 2 are unit vectors defined by As 2.21b implies, the direction along the edge E i in the space of barycentric coordinates is given by Application of 2.15b and 2.16 to 2.13 yields the position vector in the physical space as the following homogeneous polynomial:

2.24
Here are constant vectors.Equation 2.24 is equivalent to the quadratic Bernstein form where i {i 0 , i 1 , i 2 } T is a multi-index, whose valid values form the set I 2,2 , and denote the Bernstein basis functions.Substituting the above into 2.26 and comparing the result with 2.24 yield the coefficients

2.28
which are the quadratic Bézier control points for the C 0 patch.Symmetry of 2.24 greatly simplifies the formulation.For instance, a position vector on any edge E i can be expressed through 2.18b as

The G 1 Patch
Despite its extreme simplicity, the aforementioned interpolation has high geometric accuracy.
Hence it has already found many engineering applications as mentioned in Section 1.However, it generally gives discontinuous normals across patches, even if the original surface is smooth.In this study, parallel normals are considered to be identical.This may not be desirable when visual quality needs to be respected or the tangent-plane continuity is essential for the analysis.To overcome the drawback, the G 1 patch has been designed.It is obtained through modifying the C 0 patch to conform with the following conditions, which are necessary for the complete locality of the algorithm.
A The C 0 and G 1 patches have the same boundary.
B The normals of the G 1 patch on its edge depend only on the position and normal vectors at the endpoints of the edge.
After mathematical preliminaries, the normal of B is defined in Section 2.2.2.The G 1 patch is then determined from the boundary conditions.

Tangent Vector and Cross-Boundary Derivative along the Edges
For the subsequent formulation, directional derivative see, e.g., 6, page 86

2.31
Consider a triangular patch described by x β .Selecting the vector e i of 2.23 as e in 2.30 , the tangent vector along the edge E i can be written as

2.32a
Replacing the above e i with the vector e j for the next edge leads to the cross-boundary derivative Note that t i • and s i • defined above are linear operators generating functions from their arguments.Vector product of 2.32a and 2.32b gives the normal of the patch along the edge E i .
For the C 0 patch of 2.24 , its partial derivative with respect to the barycentric coordinate β i is Substituting this into 2.32a -2.32b yields the tangent vector and the cross-boundary derivative as

2.34b
where the following constant vectors are introduced:

2.35b
The right-hand sides of 2.35a are obtained through 2.25b .Equations 2.34a -2.34b and 2.18a -2.18b reveal that −T ij , S ij represent the tangent vector and the cross-boundary derivative, respectively, at the vertex V j on the edge E i .Similarly, the tangent vectors of the edges E j and E k {i, j, k} ∈ i c meeting at the vertex V i are T ji and −T ki , respectively.Their normalized vector product gives the unit normal of the C 0 patch 2.24 at the vertex V i .By definition, it is perpendicular to the tangent vector 2.35a and the cross-boundary derivative 2.35b at the same location, that is,

2.37b
For ordinary cases where the C 0 patch reproduces the normals of the original surface at the vertices, Figure 2 a indicates

2.38
and hence it is unnecessary to evaluate 2.36 in practice.Moreover, if the surface has the unique normal n i at the vertex V i , it naturally becomes parallel to the one for the adjacent patch.In the subsequent discussion, such smooth surfaces are assumed to focus on the tangent-plane continuity.

Normals on the Edges
The most important feature of the C 0 patch summarized in Section 1 is complete locality; that is, positional continuity across shared edges is automatically assured without requiring the coefficients of the adjacent patches.This is an outcome of the fact that each edge curve is constructed only from the position and normal vectors at its endpoints.If the same information can also determine the tangent plane on the patch boundary, smoothness can be achieved retaining the locality.
According to the definition of 2.36 , the unit normals at the endpoints V j and V k of the edge E i are n j and n k , respectively.These normals tend to be parallel when the edge becomes shorter, if the original surface is smooth.Therefore, the convergence is not lost by an assumption that the normal of the G 1 patch at an arbitrary point on the edge E i is a linear combination of n j and n k as where f and g are unknown functions.The above normal must be orthogonal to the tangent vector of 2.34a , and hence the following needs to be satisfied everywhere on the edge E i :

2.40
Here 2.37a is employed, and the constants are defined.Equation 2.40 can always be established if the unknown functions are selected as Here a i is an arbitrary constant, but the choice normalizes 2.42a and hence is recommended.Based on 2.42a -2.42b and 2.39 , 43 is adopted for the normal vector of the G 1 patch on the edge E i with the following coefficients:

2.44
The special treatment of the case c i d i 0 prevents the normal from vanishing.The corresponding coefficients c i d i 1 make the normal of 2.43 equivalent to 2.39 with f β j , g β k .This satisfies the orthogonality condition of 2.40 and hence causes no problem.
If only one of c i and d i is zero, the normal of 2.43 may vanish at an endpoint of the edge.However, this does not imply difficulty because the normals at the vertices are predefined and hence unchanged after the interpolation; that is, the orthogonality is not necessary to be specified at the endpoints.Thus it is concluded that the proposed normal vector is adequate as long as the C 0 patch matches the boundary conditions.

Conditions on the Correction
The next step is to find a patch whose normal on its boundary is 2.43 .Note that the equation depends only on the edge curve and the normals at its endpoints.Hence the normals automatically coincide on an edge shared by adjacent patches, and tangent-plane continuity is attained for ordinary cases where the original surface is smooth and its normals at the vertices are reproduced by the C 0 patch.
The patch in question is represented by the following function of the barycentric coordinate vector β:

2.45
On the right-hand side, x is given by 2.24 describing the C 0 patch.Since the patch is generally not orthogonal to the normal of 2.43 , the correction term g is added.As mentioned in Section 2.2, this modification is carried out without changing the boundary curves.Therefore, the correction g needs to vanish on the edges:

2.46
This preserves the tangent vector t i x of 2.34a , which is already perpendicular to N i because of 2.40 .If the cross-boundary derivative the tangent-plane continuity is achieved.The first term on the right-hand side can be simplified through 2.34b , 2.37b , and 2.43 as where the following constants are defined:

2.49
The correction g should be chosen to fulfill

2.54
The above derivative is subject to the condition of 2.48a ; that is, Multiplying the above by its denominator and substituting 2.48b into the result yields which is necessary for the tangent-plane continuity, together with 2.51 demanding that γ vanish on the edges.Due to the boundary condition and 2.18b , γ 0 must hold when at least one of the barycentric coordinates β 0 , β 1 , β 2 becomes zero.Therefore, γ can be written in the following form where G is a polynomial:

2.56b
Once the polynomials γ and ω satisfying 2.56a -2.56b are determined, the parametric representation x g of the G 1 patch is obtained through 2.50 and 2.45 .Among an infinite number of such candidates, a simple function is designed in the subsequent consideration.A monomial involved in the polynomial γ of 2.56b has the partial derivative which gives the cross-boundary derivative through 2.32b as

2.58
Assuming that 2.56b is at most quartic, the degree of G on its right-hand side is one or zero.This gives where γ 0 , γ 1 , and γ 2 are unknown coefficients.The cross-boundary derivative of the above polynomial is obtained through 2.58 as Substituting this and 2.43 into the condition 2.56a -2.56b results in whose right-hand side is quartic or lower.Hence ω on the left-hand side is at most quadratic and can be described similarly to 2.24 as where W i and w i are constants.This leads to

2.63
Substituting the above into 2.61 and comparing the coefficients arrives at the condition which is required to reproduce the normal of 2.43 .Advancing the suffixes of 2.64a by one gives C j W k c j n T k γ k , whose ratio to 2.64b is

2.65
Unfortunately, the above equation cannot hold in general since the values of both sides are predefined by 2.44 and 2.49 .This fact makes it difficult to realize 2.59 based on the hypothesis that γ is quartic or lower.Hence, a simple quintic polynomial is considered for γ in the next section.

Example of the Correction
Because of 2.58 , a term in γ needs to be linear with respect to the barycentric coordinate β i to have nonzero cross-boundary derivatives on the edge E i .If the term is quadratic or higher in β j and β k {i, j, k} ∈ i c , it does not affect the cross-boundary derivatives on the other edges.This observation motivates the choice of the quintic polynomial

2.66
where γ 0 , γ 1 , and γ 2 are constant unknown vectors.It is obvious that the above conforms with 2.56b .Equation 2.58 provides the cross-boundary derivative of 2.66 as

2.67
Note that the above result has the coefficient γ i only.Hence the corresponding term can exclusively adjust the normal on the edge E i .Substitution of 2.67 and 2.43 reduces the condition 2.56a to

2.68
Here the right-hand side consists of two scalar terms, and the unknown vector γ i already has three components.Therefore, the polynomial ω on the left-hand side is arbitrary, provided that its structure is consistent with the other side.If 2.62 is assumed again, 2.63 and 2.18b yield

2.69
Substituting the above into 2.68 and comparing the coefficients lead to the following conditions:

2.71
The value of C i in 2.70 is predefined by 2.49 and hence has no degree of freedom.Therefore, the following is generally required:

2.72a
Consequently, only the terms with the coefficient w i remain on the left-hand sides of 2.71 .The coefficient is in proportion to the magnitude of γ i on the right-hand side, and hence can have any nonzero value without loss of generality.The simplest choice is

2.72b
Substitution of 2.72a -2.72b into 2.62 determines the denominator polynomial as and reduces 2.71 to the simple linear problem

2.74
Since the above equation is underdetermined two conditions for three unknowns , it is generally satisfied by an infinite number of candidates.Among them, there is always a unique minimum-error minimum-norm solution

2.75
In general, minimization criteria assigned to redundant degrees of freedom have positive effects improving robustness and convergence; for example, 2.4 used to determine the curvature parameter is the key to the robustness of the C 0 patch 1 .Hence the above solution is adopted to complete the correction term g through 2.73 , 2.66 , and 2.50 .Formula 8 in 1 explicitly gives the generalized inverse on the right-hand side of 2.75 as

2.76c
Here 2.76b relies on the fact that 2.36 is a unit vector.Note that D 0 , D 1 of 2.76c are dimensionless due to the fact that n j , n k in 2.76b are unit vectors, and that the coefficients c i , d i are normalized by 2.44 and 2.42b , and should be treated as zeros when their absolute values are below a small positive threshold ε, for avoiding numerical instability.
The proposed surface correction does not modify the boundary of the C 0 patch, and hence its geometric representation capability is rather limited.Therefore, too large magnitudes of the coefficient vectors γ 0 , γ 1 , γ 2 imply failure in the aforementioned singularity detection or insufficient resolution of the original data.Since the coefficients have the dimension of length, such difficulties can be detected similarly to 2.6 by the condition where δ g is a dimensionless constant.If the above holds, the next expression on the righthand side of 2.76a is used for evaluating the generalized inverse.
It is important to verify the consistency of the solution 2.75 with 2.74 for tangentplane continuity.If the vectors c i n j , d i n k are linearly independent, the condition is fulfilled irrespective of the values on its left-hand side: this yields D 0 / 0 due to 2.76b , 2.76c and hence corresponds to the first case in 2.76a .The solution always gives γ i by 2.75 satisfying 2.74 .The vectors c i n j , d i n k can be linearly dependent in the following two situations.
A n j n k .
The tangent vector −T ij and the cross-boundary derivative S ij at the vertex V j see 2.34a -2.34b become perpendicular also to the normal n k at the vertex V k .The same is true if the vertices V j and V k are interchanged.Therefore,

2.78
This reduces 2.75 to γ i 0, which satisfies the condition 2.74 .
B n j ∦ n k and at least one of c i , d i is 0.
Assume that c i 0. Then 2.49 yields C i d i n T k S ij and the first row of 2.74 gives The above discussion concludes that tangent-plane continuity can be attained except for the following two pathological cases As mentioned at the end of Section 2.2.1, it is assumed that the C 0 patch reproduces the normals at its vertices.: These do not occur in meshes discretized properly.Assume 2.79a , for instance.Then 2.43 makes N i and n k parallel and hence the whole edge E i is on the tangent plane at the vertex V k .The plane, however, does not contain the cross-boundary derivative S ij at the other vertex V j , due to 2.79a .This implies that the original curve of the edge E i has displacement away from the tangent plane near the vertex V j and cannot be approximated by the simple quadratic function 2.29 .Interchanging V j and V k leads to the case of 2.79b .Hence the edge E i should be subdivided if 2.79a or 2.79b is encountered.
Due to 2.73 , 2.17 , and 2.15b , the denominator polynomial ω is positive on the patch excluding its vertices.Although it vanishes at the vertices see 2.18a , smoothness of the correction term g can be maintained, as will be clarified in the following.Consider an arbitrary path approaching the vertex V i of 2.21a in the space of barycentric coordinates.The position vector to a point on the path can be described as ωh r 0 1 r < p q, i, j, k ∈ i c .

2.84b
The above with r 0 as well as 2.66 and 2.50 gives the limit of the correction term lim which does not depend on the path.Therefore, adopting it as the value of g at the vertex V i to redefine 2.50 as preserves the continuity.This function is also smooth: substitution of e 1 i into 2.52 yields the partial derivative .66 and 2.73

2.87
whose limit vanishes at every vertex because of 2.84a -2.84b with r 0 and the boundedness of 2.81 ; that is, lim

2.88
This and 2.30 tell that every directional derivative also has the limit 0 at the vertices.Furthermore, the limit coincides with the value obtained directly from the definition 2.30 : .66 and 2.50 0 ∵ 2.84b with r 1 ,

2.89
where e is an arbitrary vector.This clearly indicates that the correction g never changes the normals at the vertices.Equation 2.89 concludes that the following extended definition of the directional derivative becomes continuous at every vertex:

2.90
This assures the smoothness of the correction term g of 2.86 .The directional derivative on the right-hand side of 2.90 can be evaluated through substituting 2.87 into 2.30 .The partial derivatives with respect to the surface parameters η, ζ can also be obtained similarly through 2.31 .

Algorithm
Finally, the procedure for determining the G 1 patch is summarized below.
A Give the position and normal vectors at the three vertices of the patch.
B Apply the algorithm described in Section 2.1 to determine the coefficients of 2.13 for the C 0 patch and obtain the vectors P i ,

Numerical Examples
Figure 3 shows a typical example of the proposed interpolation.Figure 3 a is a triangular mesh of an ellipsoid with semiaxes of lengths 1, 2, and 3.It has only six vertices and hence is essentially an octahedron.Application of the C 0 patch to the mesh with the normals to the original surface at the vertices gives Figure 3 b .It is closer to the ellipsoid, but its tangent plane is not continuous.This problem is then solved by the G 1 patch algorithm of Section 2.2.5 yielding a completely smooth surface in Figure 3 c .Next, a more general mesh of Figure 4 a is tried.It is obtained through discretization of a solid model and has the position and normal vectors of the original surface at its vertices.The mesh, however, lacks all the other geometric information.The G 1 patch interpolates it as Figure 4 b .Note that the sharp features as well as smoothness of the original model are reproduced properly.This is due to the capability of the proposed approach inherited from the C 0 patch which enables handling an arbitrary number of normals at a vertex.As already mentioned in Section 1, the present methodology has been developed mainly for computational science and engineering.Hence quantitative accuracy has the primary significance rather than visual effect.It is assessed in the following steps.First, meshes for a unit sphere, a right circular cone, and a cylinder with both unit base radius and height, and a torus whose major and minor radii are 2 and 1, respectively, are generated with diverse resolutions.They are processed by the proposed scheme, and geometric errors of the interpolation are evaluated.The result is shown in Figure 5.Both axes are in logarithmic scale.In the figure, the error δ max is the maximum distance of the interpolated surface from the exact geometry, computed using a sufficient number of sample points regardless of the mesh resolution: more than 44 × 10 4 points for the cone; more than 88 × 10 4 points for the cylinder and the sphere; more than 13 × 10 4 for the torus, uniformly distributed on the parameter planes of the patches with all their boundaries involving the singular points, i.e., the apex of the cone and the sharp edges of the discal bases .The mesh scale Δ is defined as the square root of the average patch area.The solid and hollow symbols represent the results by the G 1 and C 0 patches, respectively, both showing decrease in the error δ max with the mesh scale Δ, that is, convergence to the exact surface.The order of convergence for the G 1 patch varies between 2.2 and 4, depending on the types of the surfaces.Its decrease may be caused by the existence of isolated singular points apexes or saddle points.It is notable that the proposed algorithm maintains its stability even for extremely coarse meshes.This robustness may be attributed to the minimization criteria introduced in 2.4 and 2.75 .
For most of the cases in Figure 5, the G 1 and C 0 patches attain the same accuracy.This implies that the maximum errors tend to occur on the patch boundaries shared by both.The discrepancies between the two methods observed for the cone are considered to be due to the existence of the apex.Since normals at singular points have arbitrariness, their definition can change the convergence rate.The performance of the interpolator can naturally be affected also by the mesh subdivision scheme to control the resolution.
Regular points on a surface can be classified into elliptic, parabolic, and hyperbolic points.Since the four models used in the error analysis cover all these types as well as typical singular features, the present technique is expected to have similar effectiveness for general geometries.

Concluding Remarks
The C 0 Nagata patch has been extended to establish tangent-plane continuity, retaining its complete locality and other merits.The proposed interpolation is able to stably obtain the same level of accuracy as that of the C 0 patch, which is already sufficient even for high-precision engineering, as illustrated in Section 1.Nevertheless, there still exist problems demanding smoothness preserved in discretized models.The G 1 patch will find applications in such areas.They include optical and contact analyses, where artificial normal discontinuities may give rise to unwanted jumps in refraction angles and friction forces, respectively.

Journal of Applied Mathematics
As can be imagined from the varying nonintegral order of convergence observed in Figure 5, its theoretical estimation is difficult.Since the algorithm involves the generalized inverse, whose definition changes according to the rank of the matrix, the discussion would be fairly complex and hence deserves future study.
Similar to other interpolators, the performance of the present technique can be affected by the mesh density distribution.In this context, further research efforts will be directed also to adaptive mesh refinement with assured upper bound of the geometric errors due to the approximation by the proposed algorithm.

Figure 1 :
Figure 1: Interpolation of a curve segment on a surface.

Figure 2 :
Figure 2: Interpolation of a triangular patch using the normals at its vertices.
.30 becomes handy.Here e is an arbitrary nonzero vector along the direction of the differentiation.Recalling 2.16 and 2.20 , the partial derivatives with respect to the surface parameters η and ζ are represented by the above operator as e ζ f, e ζ ≡ {0, −1, 1} T .

Figure 3 :
Figure 3: Application of the interpolation to a mesh of a smooth surface.

Figure 4 : 1 δ
Figure 4: Application of the interpolation to a mesh of a surface with sharp features.

Figure 5 :
Figure 5: Interpolation error of the Nagata patches as a function of the mesh scale.

c d 2 , n 10 , n 11 , d 3 ≡ x 11 − x 00 , c 3 ≡ c d 3 , n 00 , n 11 .
1 i he, e ≡ {e 0 , e 1 , e 2 } T , 2.80where h is the distance from the vertex V i β 1 i in the space.The unit vector e is variable and arbitrary as long as the above β satisfies 2.15b ; that is, If one of c i , d i is zero and C i / 0, the tangent-plane continuity on the edge E i is not achieved.Then the edge E i should be subdivided since the case implies lack of mesh resolution see the comment on 2.79a -2.79b .E Determine the unknown vectors γ i i 0, 1, 2 through 2.75 and 2.76a -2.76c .Note that D k should be treated as zero if |D k | < ε or 2.77 holds for positive dimensionless thresholds ε and δ g .F Equations 2.73 , 2.66 , and 2.50 then give the correction g.Adding it to the function x of 2.24 for the C 0 patch yields the parametric representation x g of 2.45 for the G 1 patch.