# Smooth Local Interpolation of Surfaces Using Normal Vectors

Authors: Takashi Nagata
Journal: Journal of Applied Mathematics
 (2010)

---

## Abstract

This paper proposes a simple surface interpolation attaining tangent-plane continuity. It is a natural extension of the local quadraticC0 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.


---

## Body

## 1. 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],(B)
simulation of elastoplastic mechanics  [3, 4],(C)
ray tracing of optical devices  [5].All those applications prohibited the usage of traditional sophisticated surface descriptions, due to severe tolerance as well as geometrical and physical complexity of the systems. In (A) the method readily yielded positional accuracy of around10-13 m, which could not be attained by other interpolators showing poor or no convergence. Less than 10% of the linear surface patches turned out to be necessary to obtain the sufficient precision in (B). The technique reduced the number of the patches required for the ray tracing in (C) by a factor of 109. The interpolator is now referred to as C0 Nagata patch by its users, and hence this paper follows this convention, though the name of the author is frequently dropped for simplicity. The characteristics of the C0 patch listed below have enabled the aforementioned significant performance improvements.   (i)
The algorithm is completely local, requiring only the position vectors and normals given at the vertices of each patch. Therefore, it is suitable for parallel processing.(ii)
The algorithm is simple, computationally inexpensive, and hence amenable to various geometrical and physical evaluations.(iii)
Because the formulation accounts for discontinuity (multiplicity) of normals, sharp edges and singular points as well as non-manifolds can be treated quite easily.(iv)
It has the minimum degree (two) of interpolation necessary for representation of the curvature. This property is desirable especially for ray tracing, contact problems, and so forth, which involve implicitization and inversion, since closed-form solutions may be obtained.(v)
The interpolation assures theC0 continuity, and converges to the original surface rapidly with the increase in the number of nodes, even in the presence of the singular features. Hence error in the normals can be sufficiently small using rather few patches. This implies asymptotic smoothness.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 theC0 patch is that it can easily be extended to establish tangent-plane continuity, as reported in this paper. The proposed interpolation is named G1 Nagata patch here for convenience of reference.(The term does not imply that the interpolation always generates G1 surfaces; it can reproduce singular points (apexes, sharp edges, etc.) or may yield cusps similar to the C0 patch.) The new algorithm involves correction using a simple rational function retaining the boundary of the C0 patch. Hence it is no longer quadratic, but it inherits all other desirable features from the C0 patch, such as complete locality and capability of handling multiple normals. Similar to the C0 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 theC0 Nagata patch, identifies the conditions on its correction for the tangent-plane 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.


## 2. 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 (C0 Nagata patch) reproducing the modified boundary to keep the C0 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 C0 patch is then upgraded to enable local smooth interpolation.


### 2.1. TheC0 Patch

For the time being, assume a smoothsurface. Consider a curve segment on the surface as illustrated in Figure1. Its endpoints P0 and P1 have the position vectors x0, x1  and the unit normals n0, n1  (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.Figure 1
Interpolation of a curve segment on a surface.The step (A) mentioned at the beginning of Section2 approximates the curve segment using the quadratic function  (2.1)x(ξ)≡x0+(d-c)ξ+cξ2(0≤ξ≤1)
to reproduce the position and normal vectors at the endpoints P0, P1. Here ξ∈[0,1] denotes the parameter for the curve, and  (2.2)d≡x1-x0
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 P0, P1 as
(2.3)c(d,n0,n1)={[n0,n1]1-c2[1-c-c1]{n0Td-n1Td}(c≡n0Tn1≠±1),0(otherwise).
Here 0 is the zero vector. For general surfaces which may have multiple normals at P0, P1, the curvature parameter can be extended as(2.4)c=[n]+{d},
where [n], {d} are constructed by (32) of [1] from the position and normal vectors at P0, P1. The generalized inverse of the former is explicitly written as
(2.5)[n]+={1D(0)[a](0)[n]T(D(0)≠0),1D(1)[a](1)[n]T(D(0)=0,D(1)≠0),1D(2)[n]T(D(0)=0,D(1)=0,D(2)≠0),[O](otherwise),
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 P0 and P1 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 P0P1. Therefore,(2.6)δc2(dTd)<cTc
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 Figure2(a)  is interpolated. Its geometry can be described by(2.7){x(η,ζ)∣0≤ζ≤η≤1},
where x denotes the position vector to a point on the patch, and η, ζ are its parameters defined within the domain of Figure 2(b).Interpolation of a triangular patch using the normals at its vertices.
(a)(b)In Figure2(a), the position vectors to the vertices v1, v2, v3:    (2.8)x00≡x(0,0),x10≡x(1,0),x11≡x(1,1)
and the unit normals at the three points(2.9)n00≡n(0,0),n10≡n(1,0),n11≡n(1,1)
are assumed to be known. The C0 patch approximates the surface using the quadratic polynomial  (2.10)x(η,ζ)=c00+c10η+c01ζ+c11ηζ+c20η2+c02ζ2.
The edges of the patch are first interpolated in the form of (2.1) as    (2.11)x(η,0)=x00+(d1-c1)η+c1η2,x(1,ζ)=x10+(d2-c2)ζ+c2ζ2,x(η,η)=x00+(d3-c3)η+c3η2
whose coefficients are given by    (2.12)d1≡x10-x00,c1≡c(d1,n00,n10),d2≡x11-x10,c2≡c(d2,n10,n11),d3≡x11-x00,c3≡c(d3,n00,n11).
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)x(η,ζ)=x00(1-η)+x10(η-ζ)+x11ζ-c1(1-η)(η-ζ)-c2(η-ζ)ζ-c3(1-η)ζ.
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)η0≡{00},η1≡{10},η2≡{11}.
Hence the parameter vector for a point on the patch can be described as(2.15a)η={ηζ}=β0η0+β1η1+β2η2={β1+β2β2},
where β0, β1, β2 are the barycentric coordinates for η. They are in one-to-one correspondence with the surface parameters provided that
(2.15b)β0+β1+β2=1.Equations (2.15a)-(2.15b) can be solved for the barycentric coordinates as  (2.16)β0=1-η,β1=η-ζ,β2=ζ.
Equations (2.16) and (2.7) show that the following holds on the patch:  (2.17)0≤β0,β1,β2.
Let Vi be the vertex for the parameter vector ηi and Ei its opposite edge. Equations (2.15a)-(2.15b) give their representation in the barycentric coordinates as(2.18a)Vi:βi=1,βj=βk=0,(2.18b)Ei:βi=0,βj+βk=1({i,j,k}∈ic).Here ic denotes the set of all cyclic permutations of the indices 0, 1, and 2; that is, (2.19)ic≡{{0,1,2},{1,2,0},{2,0,1}}.
Introducing the barycentric coordinate vector
(2.20)β≡{β0,β1,β2}T,
equations (2.18a)-(2.18b) can be rewritten as(2.21a)Vi:β=1i,(2.21b)Ei:β=1jβj+1k(1-βj)=eiβj+1k({i,j,k}∈ic),where 1i(i=0,1,2) are unit vectors defined by(2.22)10≡{100},11≡{010},12≡{001}.
As (2.21b) implies, the direction along the edge Ei in the space of barycentric coordinates is given by(2.23)ei≡1j-1k({i,j,k}∈ic).
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)x(β)=(x00β0+x10β1+x11β2)(β0+β1+β2)-c1β0β1-c2β1β2-c3β2β0=∑i=02Piβi2+∑{i,j,k}∈icQiβjβk.
Here(2.25a)P0≡x00,P1≡x10,P2≡x11,(2.25b)Q0≡x10+x11-c2=P1+P2-C0,Q1≡x11+x00-c3=P2+P0-C1,Q2≡x00+x10-c1=P0+P1-C2,(2.25c)C0≡c2,C1≡c3,C2≡c1are constant vectors. Equation (2.24) is equivalent to the quadratic Bernstein form(2.26)x(β)=∑i∈I2,2Pi0i1i2ℬ2,i(β),
where i={i0,i1,i2}T is a multi-index, whose valid values form the set I2,2, and  (2.27)ℬ2,i(β)={β02(i={2,0,0}T),β12(i={0,2,0}T),β22(i={0,0,2}T),2β0β1(i={1,1,0}T),2β1β2(i={0,1,1}T),2β2β0(i={1,0,1}T)
denote the Bernstein basis functions. Substituting the above into (2.26) and comparing the result with (2.24) yield the coefficients(2.28)P200≡P0,P020≡P1,P002≡P2,P011≡Q02,P101≡Q12,P110≡Q22,
which are the quadratic Bézier control points for the C0 patch.Symmetry of (2.24) greatly simplifies the formulation. For instance, a position vector on any edge Ei can be expressed through (2.18b) as  (2.29)Ei:x=Pjβj2+Pkβk2+Qiβjβk({i,j,k}∈ic).


### 2.2. TheG1 Patch

Despite its extreme simplicity, the aforementioned interpolation has high geometric accuracy. Hence it has already found many engineering applications as mentioned in Section1. 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 G1 patch has been designed. It is obtained through modifying the C0 patch to conform with the following conditions, which are necessary for the complete locality of the algorithm.   (A)
TheC0 and G1 patches have the same boundary.(B)
The normals of theG1 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 Section2.2.2. The G1 patch is then determined from the boundary conditions.


#### 2.2.1. Tangent Vector and Cross-Boundary Derivative along the Edges

For the subsequent formulation, directional derivative (see, e.g., [6, page 86])(2.30)Def(β)≡limh→0f(β+he)-f(β)h=∂f∂βe
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    (2.31)∂f∂η=∂f∂β∂β∂η=∂f∂βeη=Deηf,eη≡{-1,1,0}T,∂f∂ζ=∂f∂β∂β∂ζ=∂f∂βeζ=Deζf,eζ≡{0,-1,1}T.
Consider a triangular patch described by x(β). Selecting the vector ei of (2.23) as e in (2.30), the tangent vector along the edge Ei can be written as(2.32a)ti[x]≡Deix(β)|Ei=(∂x∂βei)|βi=0=(∂x∂βj-∂x∂βk)|βi=0({i,j,k}∈ic).
Replacing the above ei with the vector ej for the next edge leads to the cross-boundary derivative(2.32b)si[x]≡Dejx(β)|Ei=(∂x∂βej)|βi=0=(∂x∂βk-∂x∂βi)|βi=0({i,j,k}∈ic).Note that ti[·] and si[·] 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 Ei.For theC0 patch of (2.24), its partial derivative with respect to the barycentric coordinate βi is   (2.33)∂x∂βi=2Piβi+Qjβk+Qkβj({i,j,k}∈ic).
Substituting this into (2.32a)-(2.32b) yields the tangent vector and the cross-boundary derivative as(2.34a)ti[x]=((2Pjβj+Qkβi+Qiβk)-(2Pkβk+Qiβj+Qjβi))|βi=0=-Tijβj+Tikβk({i,j,k}∈ic),(2.34b)si[x]=(2Pkβk+Qiβj+Qjβi)-(2Piβi+Qjβk+Qkβj)|βi=0=Sijβj+Sikβk({i,j,k}∈ic),where the following constant vectors are introduced:(2.35a)Tij≡Qi-2Pj=Pk-Pj-Ci,Tik≡Qi-2Pk=Pj-Pk-Ci({i,j,k}∈ic),(2.35b)Sij≡Qi-Qk=Tij-Tkj,Sik≡2Pk-Qj=-Tjk({i,j,k}∈ic).The right-hand sides of (2.35a) are obtained through (2.25b). Equations (2.34a)-(2.34b) and (2.18a)-(2.18b) reveal that -Tij, Sij represent the tangent vector and the cross-boundary derivative, respectively, at the vertex Vj on the edge Ei. Similarly, the tangent vectors of the edges Ej and Ek({i,j,k}∈ic) meeting at the vertex Vi are Tji and -Tki, respectively. Their normalized vector product(2.36)n̂i≡Tki×Tji|Tki×Tji|({i,j,k}∈ic)
gives the unit normal of the C0 patch  (2.24)  at the vertex Vi. 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.37a)n̂iTTji=n̂iTTki=0,(2.37b)n̂iTSji=n̂iTSki=0({i,j,k}∈ic).For ordinary cases where the C0 patch reproduces the normals of the original surface at the vertices, Figure 2(a) indicates  (2.38)n00=n̂0,n10=n̂1,n11=n̂2,
and hence it is unnecessary to evaluate (2.36) in practice. Moreover, if the surface has the unique normal n̂i at the vertex Vi, 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.


#### 2.2.2. Normals on the Edges

The most important feature of theC0 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 Vj and Vk of the edge Ei 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 G1 patch at an arbitrary point on the edge Ei is a linear combination of n̂j and n̂k as  (2.39)N̂i≡fn̂j+gn̂k({i,j,k}∈ic),
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 Ei:  (2.40)0=N̂iTti[x]=(fn̂j+gn̂k)T(-Tijβj+Tikβk)=f(diβk)-g(ciβj).
Here (2.37a) is employed, and the constants  (2.41)ci≡n̂kTTij,di≡n̂jTTik({i,j,k}∈ic)
are defined. Equation (2.40) can always be established if the unknown functions are selected as(2.42a)f=aiciβj,g=aidiβk({i,j,k}∈ic).
Here ai is an arbitrary constant, but the choice  (2.42b)ai={1(ci=di=0),1|ci|+|di|(otherwise)normalizes (2.42a) and hence is recommended. Based on (2.42a)-(2.42b) and (2.39), (2.43)N̂i≡ĉin̂jβj+d̂in̂kβk({i,j,k}∈ic)
is adopted for the normal vector of the G1 patch on the edge Ei with the following coefficients:(2.44){ĉi=d̂i=1(ci=di=0),ĉi=aici,d̂i=aidi(otherwise).
The special treatment of the case ci=di=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 ofci and di 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 C0 patch matches the boundary conditions.


#### 2.2.3. 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 C0 patch.The patch in question is represented by the following function of the barycentric coordinate vectorβ: (2.45)xg(β)=x(β)+g(β).
On the right-hand side, x is given by (2.24) describing the C0 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)g=0onE0,E1,E2.
This preserves the tangent vector ti[x] of (2.34a), which is already perpendicular to N̂i because of (2.40). If the cross-boundary derivative(2.47)si[xg]=si[x]+si[g]
is also normal to N̂i, that is, (2.48a)0=N̂iTsi[xg]=N̂iTsi[x]+N̂iTsi[g],
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  (2.48b)N̂iTsi[x]=(ĉin̂jβj+d̂in̂kβk)T(Sijβj+Sikβk)=Ĉiβjβk({i,j,k}∈ic),where the following constants are defined: (2.49)Ĉi≡ĉin̂jTSik+d̂in̂kTSij({i,j,k}∈ic).
The correction g should be chosen to fulfill (2.46) and (2.48a)-(2.48b). For sufficient generality, an arbitrary rational function (2.50)g(β)=γ(β)ω(β)
is assumed. Here γ, ω are polynomials, and the former is a vector function vanishing on the edges due to (2.46); that is,(2.51)γ=0onE0,E1,E2.
A directional derivative of (2.50) can be evaluated as  (2.52)De(γ(β)ω(β))=Deγ(β)ω-Deω(β)ω2γ,
and hence the cross-boundary derivative (2.32b) for the correction g becomes  (2.53)si[g]=Dejg(β)|βi=0=(Dejγ(β)ω-Dejω(β)ω2γ)|βi=0=(Dejγ(β)ω)|βi=0,
where (2.18b) and (2.51) are employed. The above can be simplified further through applying (2.32b) as  (2.54)si[g]=si[γ]ω|βi=0({i,j,k}∈ic).
The above derivative is subject to the condition of (2.48a); that is,(2.55)0=N̂iT(si[x]+si[γ]ω|βi=0).
Multiplying the above by its denominator and substituting (2.48b) into the result yields(2.56a)0=Ĉiβjβkω|βi=0+N̂iTsi[γ],
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)γ=β0β1β2G(β).Once the polynomials γ and ω satisfying (2.56a)-(2.56b) are determined, the parametric representation xg of the G1 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  (2.57)∂∂βiβ0α0β1α1β2α2=αiβiαi-1βjαjβkαk({i,j,k}∈ic,0<α0,α1,α2),
which gives the cross-boundary derivative through (2.32b) as  (2.58)si[β0α0β1α1β2α2]=(αkβkαk-1βiαiβjαj-αiβiαi-1βjαjβkαk)|βi=0={0(1<αi),-βjαjβkαk(αi=1)({i,j,k}∈ic,0<α0,α1,α2).
Assuming that (2.56b) is at most quartic, the degree of G on its right-hand side is one or zero. This gives(2.59)γ=β0β1β2(γ0β0+γ1β1+γ2β2)=γiβi2βjβk+γjβiβj2βk+γkβiβjβk2({i,j,k}∈ic),
where γ0, γ1, and γ2 are unknown coefficients. The cross-boundary derivative of the above polynomial is obtained through (2.58) as  (2.60)si[γ]=-γjβj2βk-γkβjβk2({i,j,k}∈ic).
Substituting this and (2.43) into the condition (2.56a)-(2.56b) results in  (2.61)Ĉiβjβkω|βi=0=(ĉin̂jβj+d̂in̂kβk)T(γjβj2βk+γkβjβk2)=(ĉin̂jTγj)βj3βk+(ĉin̂jTγk+d̂in̂kTγj)βj2βk2+(d̂in̂kTγk)βjβk3({i,j,k}∈ic)
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  (2.62)ω=∑i=02Wiβi2+∑{i,j,k}∈icwiβjβk,
where Wi and wi are constants. This leads to  (2.63)ω|βi=0=Wjβj2+Wkβk2+wiβjβk({i,j,k}∈ic).
Substituting the above into (2.61) and comparing the coefficients arrives at the condition(2.64a)ĈiWj=ĉin̂jTγj,(2.64b)ĈiWk=d̂in̂kTγk,(2.64c)Ĉiwi=ĉin̂jTγk+d̂in̂kTγj({i,j,k}∈ic),which is required to reproduce the normal of (2.43). Advancing the suffixes of (2.64a) by one gives ĈjWk=ĉjn̂kTγk, whose ratio to (2.64b) is(2.65)ĈjĈi=ĉjd̂i({i,j,k}∈ic).
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.


#### 2.2.4. 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 Ei. If the term is quadratic or higher in βj and βk({i,j,k}∈ic), it does not affect the cross-boundary derivatives on the other edges. This observation motivates the choice of the quintic polynomial(2.66)γ=γiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βk({i,j,k}∈ic)=β0β1β2(γ0β1β2+γ1β2β0+γ2β0β1),
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)si[γ]=-γiβj2βk2({i,j,k}∈ic).
Note that the above result has the coefficient γi only. Hence the corresponding term can exclusively adjust the normal on the edge Ei. Substitution of (2.67) and (2.43) reduces the condition (2.56a) to  (2.68)Ĉiβjβkω|βi=0=(ĉin̂jβj+d̂in̂kβk)Tγiβj2βk2({i,j,k}∈ic).
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)ω|βi=0=(Wjβj2+Wkβk2+wiβjβk)(βj+βk)=Wjβj3+(Wj+wi)βj2βk+(Wk+wi)βjβk2+Wkβk3({i,j,k}∈ic).
Substituting the above into (2.68) and comparing the coefficients lead to the following conditions:(2.70)ĈiWj=ĈiWk=0,(2.71)Ĉi(Wj+wi)=ĉin̂jTγi,Ĉi(Wk+wi)=d̂in̂kTγi({i,j,k}∈ic).
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)Wi=0(i=0,1,2).
Consequently, only the terms with the coefficient wi 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)wi=1(i=0,1,2).Substitution of (2.72a)-(2.72b) into (2.62) determines the denominator polynomial as(2.73)ω=β0β1+β1β2+β2β0,
and reduces (2.71) to the simple linear problem  (2.74){ĈiĈi}=[ĉin̂jTd̂in̂kT]γi({i,j,k}∈ic).
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)γi=[ĉin̂jTd̂in̂kT]+{ĈiĈi}({i,j,k}∈ic).
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 C0 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.76a)[ĉin̂jTd̂in̂kT]+=([ĉin̂j,d̂in̂k]+)T={1D(0)[ĉin̂j,d̂in̂k][a22-a12-a12a11](D(0)≠0),1D(1)[ĉin̂j,d̂in̂k](D(0)=0,D(1)≠0),[O](otherwise),(2.76b)a11≡ĉi2,a22≡d̂i2,a12≡ĉid̂in̂jTn̂k({i,j,k}∈ic),(2.76c)D(0)=a11a22-a122,D(1)=a11+a22.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 theC0 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(2.77)δg2(diTdi)<γiTγi(i=0,1,2),
where δg is a dimensionless constant. If the above holds, the next expression on the right-hand 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 tangent-plane continuity. If the vectors ĉin̂j, d̂in̂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̂in̂j, d̂in̂k  can be linearly dependent in the following two situations. (A)
n̂j∥n̂k.The tangent vector-Tij and the cross-boundary derivative Sij at the vertex Vj (see (2.34a)-(2.34b)) become perpendicular also to the normal n̂k at the vertex Vk. The same is true if the vertices Vj and Vk are interchanged. Therefore,
(2.78)0=n̂kTTij=n̂kTSij=n̂jTTik=n̂jTSik.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 thatĉi=0. Then (2.49) yields Ĉi=d̂in̂kTSij and the first row of (2.74) gives Ĉi=0. Since d̂i is nonzero, (2.74) can hold only if n̂kTSij=0. Similarly, n̂jTSik=0 is required when d̂i=0.The above discussion concludes that tangent-plane continuity can be attained except for the following two pathological cases (As mentioned at the end of Section2.2.1, it is assumed that the C0 patch reproduces the normals at its vertices.):(2.79a)ĉi=0,n̂kTSij≠0,(2.79b)d̂i=0,n̂kTSij≠0.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 Ei is on the tangent plane at the vertex Vk. The plane, however, does not contain the cross-boundary derivative Sij at the other vertex Vj, due to (2.79a). This implies that the original curve of the edge Ei has displacement away from the tangent plane near the vertex Vj and cannot be approximated by the simple quadratic function (2.29). Interchanging Vj and Vk leads to the case of (2.79b). Hence the edge Ei 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 Vi of (2.21a) in the space of barycentric coordinates. The position vector to a point on the path can be described as(2.80)β=1i+he,e≡{e0,e1,e2}T,
where h is the distance from the vertex Vi(β=1i) in the space. The unit vector e is variable and arbitrary as long as the above β satisfies (2.15b); that is,(2.81)|e|=1,e0+e1+e2=0.
Substituting (2.80) into (2.73) yields  (2.82)ω=(1+hei)(hej)+(hej)(hek)+(hek)(1+hei)=(ej+ek)h+O(h2)({i,j,k}∈ic),
which reveals that the denominator polynomial is first order with respect to h. Similarly by (2.80),
(2.83)βjpβkq=ejpekqhp+q
holds. If the ratio(2.84a)βjpβkqωhr=ejpekqhp+q-r-1(ej+ek)+O(h)
has a positive exponent of h in the numerator, it vanishes as it approaches the vertex β=1i for h→0 irrespective of the path; that is  (note that ej+ek never vanishes in the denominator of (2.84a) due to (2.81)), (2.84b)limβ→1iβjpβkqωhr=0(1+r<p+q,{i,j,k}∈ic).The above with r=0 as well as (2.66) and (2.50) gives the limit of the correction term   (2.85)limβ→1ig=limβ→1iγiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βkω=0,
which does not depend on the path. Therefore, adopting it as the value of g at the vertex Vi to redefine (2.50) as   (2.86)g(β)={0(β0=1orβ1=1orβ2=1),γ(β)ω(β)(otherwise)
preserves the continuity. This function is also smooth: substitution of e=1i into (2.52) yields the partial derivative(2.87)∂∂βi(γω)=1ω∂γ∂βi-1ω2∂ω∂βiγ=1ω2((βiβj+βjβk+βkβi)(γiβj2βk2+2γjβiβjβk2+2γkβiβj2βk)-(βj+βk)(γiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βk))(∵(2.66)and(2.73))=γiω2βj2βk2(βiβj+βjβk+βkβi-βi(βj+βk))+γjω2βiβjβk2(2βiβj+2βjβk+2βkβi-βi(βj+βk))+γkω2βiβj2βk(2βiβj+2βjβk+2βkβi-βi(βj+βk))=γiω2βj3βk3+γjω2βiβjβk2(ω+βjβk)+γkω2βiβj2βk(ω+βjβk)(∵(2.73))=1ω2(βj3βk3γi+βiβjβk(ω+βjβk)(βkγj+βjγk))({i,j,k}∈ic),
whose limit vanishes at every vertex because of (2.84a)-(2.84b) with r=0 and the boundedness of (2.81); that is,
(2.88)limβ→1i∂∂βi(γω)=0(i=0,1,2).
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):   (2.89)Deg(1i)=limh→0g(1i+he)-g(1i)h=limh→0g(1i+he)h(∵(2.86))=limh→0γiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βkωh(∵(2.66)and(2.50))=0(∵(2.84b)withr=1),
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)Deg(β)={0(β0=1orβ1=1orβ2=1),De(γ(β)ω(β))(otherwise).
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).


#### 2.2.5. Algorithm

Finally, the procedure for determining theG1 patch is summarized below.(A)
Give the position and normal vectors at the three vertices of the patch.(B)
Apply the algorithm described in Section2.1 to determine the coefficients of (2.13) for the C0 patch and obtain the vectors Pi,Qi(i=0,1,2) through (2.25a)–(2.25c).(C)
Compute the constantsTij,Tik,Sij,Sik({i,j,k}∈ic) of (2.35a)-(2.35b) and the unit normals n̂i(i=0,1,2) of (2.36).(D)
After calculatingĉi, d̂i through (2.44), (2.42b), and (2.41), evaluate (2.49) to get Ĉi(i=0,1,2). If one of ĉi, d̂i is zero and Ĉi≠0, the tangent-plane continuity on the edge Ei is not achieved. Then the edge Ei 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 C0 patch yields the parametric representation xg of (2.45) for the G1 patch.


## 2.1. TheC0 Patch

For the time being, assume a smoothsurface. Consider a curve segment on the surface as illustrated in Figure1. Its endpoints P0 and P1 have the position vectors x0, x1  and the unit normals n0, n1  (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.Figure 1
Interpolation of a curve segment on a surface.The step (A) mentioned at the beginning of Section2 approximates the curve segment using the quadratic function  (2.1)x(ξ)≡x0+(d-c)ξ+cξ2(0≤ξ≤1)
to reproduce the position and normal vectors at the endpoints P0, P1. Here ξ∈[0,1] denotes the parameter for the curve, and  (2.2)d≡x1-x0
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 P0, P1 as
(2.3)c(d,n0,n1)={[n0,n1]1-c2[1-c-c1]{n0Td-n1Td}(c≡n0Tn1≠±1),0(otherwise).
Here 0 is the zero vector. For general surfaces which may have multiple normals at P0, P1, the curvature parameter can be extended as(2.4)c=[n]+{d},
where [n], {d} are constructed by (32) of [1] from the position and normal vectors at P0, P1. The generalized inverse of the former is explicitly written as
(2.5)[n]+={1D(0)[a](0)[n]T(D(0)≠0),1D(1)[a](1)[n]T(D(0)=0,D(1)≠0),1D(2)[n]T(D(0)=0,D(1)=0,D(2)≠0),[O](otherwise),
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 P0 and P1 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 P0P1. Therefore,(2.6)δc2(dTd)<cTc
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 Figure2(a)  is interpolated. Its geometry can be described by(2.7){x(η,ζ)∣0≤ζ≤η≤1},
where x denotes the position vector to a point on the patch, and η, ζ are its parameters defined within the domain of Figure 2(b).Interpolation of a triangular patch using the normals at its vertices.
(a)(b)In Figure2(a), the position vectors to the vertices v1, v2, v3:    (2.8)x00≡x(0,0),x10≡x(1,0),x11≡x(1,1)
and the unit normals at the three points(2.9)n00≡n(0,0),n10≡n(1,0),n11≡n(1,1)
are assumed to be known. The C0 patch approximates the surface using the quadratic polynomial  (2.10)x(η,ζ)=c00+c10η+c01ζ+c11ηζ+c20η2+c02ζ2.
The edges of the patch are first interpolated in the form of (2.1) as    (2.11)x(η,0)=x00+(d1-c1)η+c1η2,x(1,ζ)=x10+(d2-c2)ζ+c2ζ2,x(η,η)=x00+(d3-c3)η+c3η2
whose coefficients are given by    (2.12)d1≡x10-x00,c1≡c(d1,n00,n10),d2≡x11-x10,c2≡c(d2,n10,n11),d3≡x11-x00,c3≡c(d3,n00,n11).
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)x(η,ζ)=x00(1-η)+x10(η-ζ)+x11ζ-c1(1-η)(η-ζ)-c2(η-ζ)ζ-c3(1-η)ζ.
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)η0≡{00},η1≡{10},η2≡{11}.
Hence the parameter vector for a point on the patch can be described as(2.15a)η={ηζ}=β0η0+β1η1+β2η2={β1+β2β2},
where β0, β1, β2 are the barycentric coordinates for η. They are in one-to-one correspondence with the surface parameters provided that
(2.15b)β0+β1+β2=1.Equations (2.15a)-(2.15b) can be solved for the barycentric coordinates as  (2.16)β0=1-η,β1=η-ζ,β2=ζ.
Equations (2.16) and (2.7) show that the following holds on the patch:  (2.17)0≤β0,β1,β2.
Let Vi be the vertex for the parameter vector ηi and Ei its opposite edge. Equations (2.15a)-(2.15b) give their representation in the barycentric coordinates as(2.18a)Vi:βi=1,βj=βk=0,(2.18b)Ei:βi=0,βj+βk=1({i,j,k}∈ic).Here ic denotes the set of all cyclic permutations of the indices 0, 1, and 2; that is, (2.19)ic≡{{0,1,2},{1,2,0},{2,0,1}}.
Introducing the barycentric coordinate vector
(2.20)β≡{β0,β1,β2}T,
equations (2.18a)-(2.18b) can be rewritten as(2.21a)Vi:β=1i,(2.21b)Ei:β=1jβj+1k(1-βj)=eiβj+1k({i,j,k}∈ic),where 1i(i=0,1,2) are unit vectors defined by(2.22)10≡{100},11≡{010},12≡{001}.
As (2.21b) implies, the direction along the edge Ei in the space of barycentric coordinates is given by(2.23)ei≡1j-1k({i,j,k}∈ic).
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)x(β)=(x00β0+x10β1+x11β2)(β0+β1+β2)-c1β0β1-c2β1β2-c3β2β0=∑i=02Piβi2+∑{i,j,k}∈icQiβjβk.
Here(2.25a)P0≡x00,P1≡x10,P2≡x11,(2.25b)Q0≡x10+x11-c2=P1+P2-C0,Q1≡x11+x00-c3=P2+P0-C1,Q2≡x00+x10-c1=P0+P1-C2,(2.25c)C0≡c2,C1≡c3,C2≡c1are constant vectors. Equation (2.24) is equivalent to the quadratic Bernstein form(2.26)x(β)=∑i∈I2,2Pi0i1i2ℬ2,i(β),
where i={i0,i1,i2}T is a multi-index, whose valid values form the set I2,2, and  (2.27)ℬ2,i(β)={β02(i={2,0,0}T),β12(i={0,2,0}T),β22(i={0,0,2}T),2β0β1(i={1,1,0}T),2β1β2(i={0,1,1}T),2β2β0(i={1,0,1}T)
denote the Bernstein basis functions. Substituting the above into (2.26) and comparing the result with (2.24) yield the coefficients(2.28)P200≡P0,P020≡P1,P002≡P2,P011≡Q02,P101≡Q12,P110≡Q22,
which are the quadratic Bézier control points for the C0 patch.Symmetry of (2.24) greatly simplifies the formulation. For instance, a position vector on any edge Ei can be expressed through (2.18b) as  (2.29)Ei:x=Pjβj2+Pkβk2+Qiβjβk({i,j,k}∈ic).


## 2.2. TheG1 Patch

Despite its extreme simplicity, the aforementioned interpolation has high geometric accuracy. Hence it has already found many engineering applications as mentioned in Section1. 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 G1 patch has been designed. It is obtained through modifying the C0 patch to conform with the following conditions, which are necessary for the complete locality of the algorithm.   (A)
TheC0 and G1 patches have the same boundary.(B)
The normals of theG1 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 Section2.2.2. The G1 patch is then determined from the boundary conditions.


### 2.2.1. Tangent Vector and Cross-Boundary Derivative along the Edges

For the subsequent formulation, directional derivative (see, e.g., [6, page 86])(2.30)Def(β)≡limh→0f(β+he)-f(β)h=∂f∂βe
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    (2.31)∂f∂η=∂f∂β∂β∂η=∂f∂βeη=Deηf,eη≡{-1,1,0}T,∂f∂ζ=∂f∂β∂β∂ζ=∂f∂βeζ=Deζf,eζ≡{0,-1,1}T.
Consider a triangular patch described by x(β). Selecting the vector ei of (2.23) as e in (2.30), the tangent vector along the edge Ei can be written as(2.32a)ti[x]≡Deix(β)|Ei=(∂x∂βei)|βi=0=(∂x∂βj-∂x∂βk)|βi=0({i,j,k}∈ic).
Replacing the above ei with the vector ej for the next edge leads to the cross-boundary derivative(2.32b)si[x]≡Dejx(β)|Ei=(∂x∂βej)|βi=0=(∂x∂βk-∂x∂βi)|βi=0({i,j,k}∈ic).Note that ti[·] and si[·] 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 Ei.For theC0 patch of (2.24), its partial derivative with respect to the barycentric coordinate βi is   (2.33)∂x∂βi=2Piβi+Qjβk+Qkβj({i,j,k}∈ic).
Substituting this into (2.32a)-(2.32b) yields the tangent vector and the cross-boundary derivative as(2.34a)ti[x]=((2Pjβj+Qkβi+Qiβk)-(2Pkβk+Qiβj+Qjβi))|βi=0=-Tijβj+Tikβk({i,j,k}∈ic),(2.34b)si[x]=(2Pkβk+Qiβj+Qjβi)-(2Piβi+Qjβk+Qkβj)|βi=0=Sijβj+Sikβk({i,j,k}∈ic),where the following constant vectors are introduced:(2.35a)Tij≡Qi-2Pj=Pk-Pj-Ci,Tik≡Qi-2Pk=Pj-Pk-Ci({i,j,k}∈ic),(2.35b)Sij≡Qi-Qk=Tij-Tkj,Sik≡2Pk-Qj=-Tjk({i,j,k}∈ic).The right-hand sides of (2.35a) are obtained through (2.25b). Equations (2.34a)-(2.34b) and (2.18a)-(2.18b) reveal that -Tij, Sij represent the tangent vector and the cross-boundary derivative, respectively, at the vertex Vj on the edge Ei. Similarly, the tangent vectors of the edges Ej and Ek({i,j,k}∈ic) meeting at the vertex Vi are Tji and -Tki, respectively. Their normalized vector product(2.36)n̂i≡Tki×Tji|Tki×Tji|({i,j,k}∈ic)
gives the unit normal of the C0 patch  (2.24)  at the vertex Vi. 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.37a)n̂iTTji=n̂iTTki=0,(2.37b)n̂iTSji=n̂iTSki=0({i,j,k}∈ic).For ordinary cases where the C0 patch reproduces the normals of the original surface at the vertices, Figure 2(a) indicates  (2.38)n00=n̂0,n10=n̂1,n11=n̂2,
and hence it is unnecessary to evaluate (2.36) in practice. Moreover, if the surface has the unique normal n̂i at the vertex Vi, 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.


### 2.2.2. Normals on the Edges

The most important feature of theC0 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 Vj and Vk of the edge Ei 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 G1 patch at an arbitrary point on the edge Ei is a linear combination of n̂j and n̂k as  (2.39)N̂i≡fn̂j+gn̂k({i,j,k}∈ic),
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 Ei:  (2.40)0=N̂iTti[x]=(fn̂j+gn̂k)T(-Tijβj+Tikβk)=f(diβk)-g(ciβj).
Here (2.37a) is employed, and the constants  (2.41)ci≡n̂kTTij,di≡n̂jTTik({i,j,k}∈ic)
are defined. Equation (2.40) can always be established if the unknown functions are selected as(2.42a)f=aiciβj,g=aidiβk({i,j,k}∈ic).
Here ai is an arbitrary constant, but the choice  (2.42b)ai={1(ci=di=0),1|ci|+|di|(otherwise)normalizes (2.42a) and hence is recommended. Based on (2.42a)-(2.42b) and (2.39), (2.43)N̂i≡ĉin̂jβj+d̂in̂kβk({i,j,k}∈ic)
is adopted for the normal vector of the G1 patch on the edge Ei with the following coefficients:(2.44){ĉi=d̂i=1(ci=di=0),ĉi=aici,d̂i=aidi(otherwise).
The special treatment of the case ci=di=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 ofci and di 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 C0 patch matches the boundary conditions.


### 2.2.3. 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 C0 patch.The patch in question is represented by the following function of the barycentric coordinate vectorβ: (2.45)xg(β)=x(β)+g(β).
On the right-hand side, x is given by (2.24) describing the C0 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)g=0onE0,E1,E2.
This preserves the tangent vector ti[x] of (2.34a), which is already perpendicular to N̂i because of (2.40). If the cross-boundary derivative(2.47)si[xg]=si[x]+si[g]
is also normal to N̂i, that is, (2.48a)0=N̂iTsi[xg]=N̂iTsi[x]+N̂iTsi[g],
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  (2.48b)N̂iTsi[x]=(ĉin̂jβj+d̂in̂kβk)T(Sijβj+Sikβk)=Ĉiβjβk({i,j,k}∈ic),where the following constants are defined: (2.49)Ĉi≡ĉin̂jTSik+d̂in̂kTSij({i,j,k}∈ic).
The correction g should be chosen to fulfill (2.46) and (2.48a)-(2.48b). For sufficient generality, an arbitrary rational function (2.50)g(β)=γ(β)ω(β)
is assumed. Here γ, ω are polynomials, and the former is a vector function vanishing on the edges due to (2.46); that is,(2.51)γ=0onE0,E1,E2.
A directional derivative of (2.50) can be evaluated as  (2.52)De(γ(β)ω(β))=Deγ(β)ω-Deω(β)ω2γ,
and hence the cross-boundary derivative (2.32b) for the correction g becomes  (2.53)si[g]=Dejg(β)|βi=0=(Dejγ(β)ω-Dejω(β)ω2γ)|βi=0=(Dejγ(β)ω)|βi=0,
where (2.18b) and (2.51) are employed. The above can be simplified further through applying (2.32b) as  (2.54)si[g]=si[γ]ω|βi=0({i,j,k}∈ic).
The above derivative is subject to the condition of (2.48a); that is,(2.55)0=N̂iT(si[x]+si[γ]ω|βi=0).
Multiplying the above by its denominator and substituting (2.48b) into the result yields(2.56a)0=Ĉiβjβkω|βi=0+N̂iTsi[γ],
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)γ=β0β1β2G(β).Once the polynomials γ and ω satisfying (2.56a)-(2.56b) are determined, the parametric representation xg of the G1 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  (2.57)∂∂βiβ0α0β1α1β2α2=αiβiαi-1βjαjβkαk({i,j,k}∈ic,0<α0,α1,α2),
which gives the cross-boundary derivative through (2.32b) as  (2.58)si[β0α0β1α1β2α2]=(αkβkαk-1βiαiβjαj-αiβiαi-1βjαjβkαk)|βi=0={0(1<αi),-βjαjβkαk(αi=1)({i,j,k}∈ic,0<α0,α1,α2).
Assuming that (2.56b) is at most quartic, the degree of G on its right-hand side is one or zero. This gives(2.59)γ=β0β1β2(γ0β0+γ1β1+γ2β2)=γiβi2βjβk+γjβiβj2βk+γkβiβjβk2({i,j,k}∈ic),
where γ0, γ1, and γ2 are unknown coefficients. The cross-boundary derivative of the above polynomial is obtained through (2.58) as  (2.60)si[γ]=-γjβj2βk-γkβjβk2({i,j,k}∈ic).
Substituting this and (2.43) into the condition (2.56a)-(2.56b) results in  (2.61)Ĉiβjβkω|βi=0=(ĉin̂jβj+d̂in̂kβk)T(γjβj2βk+γkβjβk2)=(ĉin̂jTγj)βj3βk+(ĉin̂jTγk+d̂in̂kTγj)βj2βk2+(d̂in̂kTγk)βjβk3({i,j,k}∈ic)
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  (2.62)ω=∑i=02Wiβi2+∑{i,j,k}∈icwiβjβk,
where Wi and wi are constants. This leads to  (2.63)ω|βi=0=Wjβj2+Wkβk2+wiβjβk({i,j,k}∈ic).
Substituting the above into (2.61) and comparing the coefficients arrives at the condition(2.64a)ĈiWj=ĉin̂jTγj,(2.64b)ĈiWk=d̂in̂kTγk,(2.64c)Ĉiwi=ĉin̂jTγk+d̂in̂kTγj({i,j,k}∈ic),which is required to reproduce the normal of (2.43). Advancing the suffixes of (2.64a) by one gives ĈjWk=ĉjn̂kTγk, whose ratio to (2.64b) is(2.65)ĈjĈi=ĉjd̂i({i,j,k}∈ic).
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.


### 2.2.4. 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 Ei. If the term is quadratic or higher in βj and βk({i,j,k}∈ic), it does not affect the cross-boundary derivatives on the other edges. This observation motivates the choice of the quintic polynomial(2.66)γ=γiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βk({i,j,k}∈ic)=β0β1β2(γ0β1β2+γ1β2β0+γ2β0β1),
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)si[γ]=-γiβj2βk2({i,j,k}∈ic).
Note that the above result has the coefficient γi only. Hence the corresponding term can exclusively adjust the normal on the edge Ei. Substitution of (2.67) and (2.43) reduces the condition (2.56a) to  (2.68)Ĉiβjβkω|βi=0=(ĉin̂jβj+d̂in̂kβk)Tγiβj2βk2({i,j,k}∈ic).
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)ω|βi=0=(Wjβj2+Wkβk2+wiβjβk)(βj+βk)=Wjβj3+(Wj+wi)βj2βk+(Wk+wi)βjβk2+Wkβk3({i,j,k}∈ic).
Substituting the above into (2.68) and comparing the coefficients lead to the following conditions:(2.70)ĈiWj=ĈiWk=0,(2.71)Ĉi(Wj+wi)=ĉin̂jTγi,Ĉi(Wk+wi)=d̂in̂kTγi({i,j,k}∈ic).
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)Wi=0(i=0,1,2).
Consequently, only the terms with the coefficient wi 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)wi=1(i=0,1,2).Substitution of (2.72a)-(2.72b) into (2.62) determines the denominator polynomial as(2.73)ω=β0β1+β1β2+β2β0,
and reduces (2.71) to the simple linear problem  (2.74){ĈiĈi}=[ĉin̂jTd̂in̂kT]γi({i,j,k}∈ic).
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)γi=[ĉin̂jTd̂in̂kT]+{ĈiĈi}({i,j,k}∈ic).
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 C0 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.76a)[ĉin̂jTd̂in̂kT]+=([ĉin̂j,d̂in̂k]+)T={1D(0)[ĉin̂j,d̂in̂k][a22-a12-a12a11](D(0)≠0),1D(1)[ĉin̂j,d̂in̂k](D(0)=0,D(1)≠0),[O](otherwise),(2.76b)a11≡ĉi2,a22≡d̂i2,a12≡ĉid̂in̂jTn̂k({i,j,k}∈ic),(2.76c)D(0)=a11a22-a122,D(1)=a11+a22.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 theC0 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(2.77)δg2(diTdi)<γiTγi(i=0,1,2),
where δg is a dimensionless constant. If the above holds, the next expression on the right-hand 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 tangent-plane continuity. If the vectors ĉin̂j, d̂in̂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̂in̂j, d̂in̂k  can be linearly dependent in the following two situations. (A)
n̂j∥n̂k.The tangent vector-Tij and the cross-boundary derivative Sij at the vertex Vj (see (2.34a)-(2.34b)) become perpendicular also to the normal n̂k at the vertex Vk. The same is true if the vertices Vj and Vk are interchanged. Therefore,
(2.78)0=n̂kTTij=n̂kTSij=n̂jTTik=n̂jTSik.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 thatĉi=0. Then (2.49) yields Ĉi=d̂in̂kTSij and the first row of (2.74) gives Ĉi=0. Since d̂i is nonzero, (2.74) can hold only if n̂kTSij=0. Similarly, n̂jTSik=0 is required when d̂i=0.The above discussion concludes that tangent-plane continuity can be attained except for the following two pathological cases (As mentioned at the end of Section2.2.1, it is assumed that the C0 patch reproduces the normals at its vertices.):(2.79a)ĉi=0,n̂kTSij≠0,(2.79b)d̂i=0,n̂kTSij≠0.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 Ei is on the tangent plane at the vertex Vk. The plane, however, does not contain the cross-boundary derivative Sij at the other vertex Vj, due to (2.79a). This implies that the original curve of the edge Ei has displacement away from the tangent plane near the vertex Vj and cannot be approximated by the simple quadratic function (2.29). Interchanging Vj and Vk leads to the case of (2.79b). Hence the edge Ei 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 Vi of (2.21a) in the space of barycentric coordinates. The position vector to a point on the path can be described as(2.80)β=1i+he,e≡{e0,e1,e2}T,
where h is the distance from the vertex Vi(β=1i) in the space. The unit vector e is variable and arbitrary as long as the above β satisfies (2.15b); that is,(2.81)|e|=1,e0+e1+e2=0.
Substituting (2.80) into (2.73) yields  (2.82)ω=(1+hei)(hej)+(hej)(hek)+(hek)(1+hei)=(ej+ek)h+O(h2)({i,j,k}∈ic),
which reveals that the denominator polynomial is first order with respect to h. Similarly by (2.80),
(2.83)βjpβkq=ejpekqhp+q
holds. If the ratio(2.84a)βjpβkqωhr=ejpekqhp+q-r-1(ej+ek)+O(h)
has a positive exponent of h in the numerator, it vanishes as it approaches the vertex β=1i for h→0 irrespective of the path; that is  (note that ej+ek never vanishes in the denominator of (2.84a) due to (2.81)), (2.84b)limβ→1iβjpβkqωhr=0(1+r<p+q,{i,j,k}∈ic).The above with r=0 as well as (2.66) and (2.50) gives the limit of the correction term   (2.85)limβ→1ig=limβ→1iγiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βkω=0,
which does not depend on the path. Therefore, adopting it as the value of g at the vertex Vi to redefine (2.50) as   (2.86)g(β)={0(β0=1orβ1=1orβ2=1),γ(β)ω(β)(otherwise)
preserves the continuity. This function is also smooth: substitution of e=1i into (2.52) yields the partial derivative(2.87)∂∂βi(γω)=1ω∂γ∂βi-1ω2∂ω∂βiγ=1ω2((βiβj+βjβk+βkβi)(γiβj2βk2+2γjβiβjβk2+2γkβiβj2βk)-(βj+βk)(γiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βk))(∵(2.66)and(2.73))=γiω2βj2βk2(βiβj+βjβk+βkβi-βi(βj+βk))+γjω2βiβjβk2(2βiβj+2βjβk+2βkβi-βi(βj+βk))+γkω2βiβj2βk(2βiβj+2βjβk+2βkβi-βi(βj+βk))=γiω2βj3βk3+γjω2βiβjβk2(ω+βjβk)+γkω2βiβj2βk(ω+βjβk)(∵(2.73))=1ω2(βj3βk3γi+βiβjβk(ω+βjβk)(βkγj+βjγk))({i,j,k}∈ic),
whose limit vanishes at every vertex because of (2.84a)-(2.84b) with r=0 and the boundedness of (2.81); that is,
(2.88)limβ→1i∂∂βi(γω)=0(i=0,1,2).
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):   (2.89)Deg(1i)=limh→0g(1i+he)-g(1i)h=limh→0g(1i+he)h(∵(2.86))=limh→0γiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βkωh(∵(2.66)and(2.50))=0(∵(2.84b)withr=1),
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)Deg(β)={0(β0=1orβ1=1orβ2=1),De(γ(β)ω(β))(otherwise).
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).


### 2.2.5. Algorithm

Finally, the procedure for determining theG1 patch is summarized below.(A)
Give the position and normal vectors at the three vertices of the patch.(B)
Apply the algorithm described in Section2.1 to determine the coefficients of (2.13) for the C0 patch and obtain the vectors Pi,Qi(i=0,1,2) through (2.25a)–(2.25c).(C)
Compute the constantsTij,Tik,Sij,Sik({i,j,k}∈ic) of (2.35a)-(2.35b) and the unit normals n̂i(i=0,1,2) of (2.36).(D)
After calculatingĉi, d̂i through (2.44), (2.42b), and (2.41), evaluate (2.49) to get Ĉi(i=0,1,2). If one of ĉi, d̂i is zero and Ĉi≠0, the tangent-plane continuity on the edge Ei is not achieved. Then the edge Ei 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 C0 patch yields the parametric representation xg of (2.45) for the G1 patch.


## 2.2.1. Tangent Vector and Cross-Boundary Derivative along the Edges

For the subsequent formulation, directional derivative (see, e.g., [6, page 86])(2.30)Def(β)≡limh→0f(β+he)-f(β)h=∂f∂βe
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    (2.31)∂f∂η=∂f∂β∂β∂η=∂f∂βeη=Deηf,eη≡{-1,1,0}T,∂f∂ζ=∂f∂β∂β∂ζ=∂f∂βeζ=Deζf,eζ≡{0,-1,1}T.
Consider a triangular patch described by x(β). Selecting the vector ei of (2.23) as e in (2.30), the tangent vector along the edge Ei can be written as(2.32a)ti[x]≡Deix(β)|Ei=(∂x∂βei)|βi=0=(∂x∂βj-∂x∂βk)|βi=0({i,j,k}∈ic).
Replacing the above ei with the vector ej for the next edge leads to the cross-boundary derivative(2.32b)si[x]≡Dejx(β)|Ei=(∂x∂βej)|βi=0=(∂x∂βk-∂x∂βi)|βi=0({i,j,k}∈ic).Note that ti[·] and si[·] 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 Ei.For theC0 patch of (2.24), its partial derivative with respect to the barycentric coordinate βi is   (2.33)∂x∂βi=2Piβi+Qjβk+Qkβj({i,j,k}∈ic).
Substituting this into (2.32a)-(2.32b) yields the tangent vector and the cross-boundary derivative as(2.34a)ti[x]=((2Pjβj+Qkβi+Qiβk)-(2Pkβk+Qiβj+Qjβi))|βi=0=-Tijβj+Tikβk({i,j,k}∈ic),(2.34b)si[x]=(2Pkβk+Qiβj+Qjβi)-(2Piβi+Qjβk+Qkβj)|βi=0=Sijβj+Sikβk({i,j,k}∈ic),where the following constant vectors are introduced:(2.35a)Tij≡Qi-2Pj=Pk-Pj-Ci,Tik≡Qi-2Pk=Pj-Pk-Ci({i,j,k}∈ic),(2.35b)Sij≡Qi-Qk=Tij-Tkj,Sik≡2Pk-Qj=-Tjk({i,j,k}∈ic).The right-hand sides of (2.35a) are obtained through (2.25b). Equations (2.34a)-(2.34b) and (2.18a)-(2.18b) reveal that -Tij, Sij represent the tangent vector and the cross-boundary derivative, respectively, at the vertex Vj on the edge Ei. Similarly, the tangent vectors of the edges Ej and Ek({i,j,k}∈ic) meeting at the vertex Vi are Tji and -Tki, respectively. Their normalized vector product(2.36)n̂i≡Tki×Tji|Tki×Tji|({i,j,k}∈ic)
gives the unit normal of the C0 patch  (2.24)  at the vertex Vi. 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.37a)n̂iTTji=n̂iTTki=0,(2.37b)n̂iTSji=n̂iTSki=0({i,j,k}∈ic).For ordinary cases where the C0 patch reproduces the normals of the original surface at the vertices, Figure 2(a) indicates  (2.38)n00=n̂0,n10=n̂1,n11=n̂2,
and hence it is unnecessary to evaluate (2.36) in practice. Moreover, if the surface has the unique normal n̂i at the vertex Vi, 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.


## 2.2.2. Normals on the Edges

The most important feature of theC0 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 Vj and Vk of the edge Ei 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 G1 patch at an arbitrary point on the edge Ei is a linear combination of n̂j and n̂k as  (2.39)N̂i≡fn̂j+gn̂k({i,j,k}∈ic),
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 Ei:  (2.40)0=N̂iTti[x]=(fn̂j+gn̂k)T(-Tijβj+Tikβk)=f(diβk)-g(ciβj).
Here (2.37a) is employed, and the constants  (2.41)ci≡n̂kTTij,di≡n̂jTTik({i,j,k}∈ic)
are defined. Equation (2.40) can always be established if the unknown functions are selected as(2.42a)f=aiciβj,g=aidiβk({i,j,k}∈ic).
Here ai is an arbitrary constant, but the choice  (2.42b)ai={1(ci=di=0),1|ci|+|di|(otherwise)normalizes (2.42a) and hence is recommended. Based on (2.42a)-(2.42b) and (2.39), (2.43)N̂i≡ĉin̂jβj+d̂in̂kβk({i,j,k}∈ic)
is adopted for the normal vector of the G1 patch on the edge Ei with the following coefficients:(2.44){ĉi=d̂i=1(ci=di=0),ĉi=aici,d̂i=aidi(otherwise).
The special treatment of the case ci=di=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 ofci and di 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 C0 patch matches the boundary conditions.


## 2.2.3. 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 C0 patch.The patch in question is represented by the following function of the barycentric coordinate vectorβ: (2.45)xg(β)=x(β)+g(β).
On the right-hand side, x is given by (2.24) describing the C0 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)g=0onE0,E1,E2.
This preserves the tangent vector ti[x] of (2.34a), which is already perpendicular to N̂i because of (2.40). If the cross-boundary derivative(2.47)si[xg]=si[x]+si[g]
is also normal to N̂i, that is, (2.48a)0=N̂iTsi[xg]=N̂iTsi[x]+N̂iTsi[g],
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  (2.48b)N̂iTsi[x]=(ĉin̂jβj+d̂in̂kβk)T(Sijβj+Sikβk)=Ĉiβjβk({i,j,k}∈ic),where the following constants are defined: (2.49)Ĉi≡ĉin̂jTSik+d̂in̂kTSij({i,j,k}∈ic).
The correction g should be chosen to fulfill (2.46) and (2.48a)-(2.48b). For sufficient generality, an arbitrary rational function (2.50)g(β)=γ(β)ω(β)
is assumed. Here γ, ω are polynomials, and the former is a vector function vanishing on the edges due to (2.46); that is,(2.51)γ=0onE0,E1,E2.
A directional derivative of (2.50) can be evaluated as  (2.52)De(γ(β)ω(β))=Deγ(β)ω-Deω(β)ω2γ,
and hence the cross-boundary derivative (2.32b) for the correction g becomes  (2.53)si[g]=Dejg(β)|βi=0=(Dejγ(β)ω-Dejω(β)ω2γ)|βi=0=(Dejγ(β)ω)|βi=0,
where (2.18b) and (2.51) are employed. The above can be simplified further through applying (2.32b) as  (2.54)si[g]=si[γ]ω|βi=0({i,j,k}∈ic).
The above derivative is subject to the condition of (2.48a); that is,(2.55)0=N̂iT(si[x]+si[γ]ω|βi=0).
Multiplying the above by its denominator and substituting (2.48b) into the result yields(2.56a)0=Ĉiβjβkω|βi=0+N̂iTsi[γ],
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)γ=β0β1β2G(β).Once the polynomials γ and ω satisfying (2.56a)-(2.56b) are determined, the parametric representation xg of the G1 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  (2.57)∂∂βiβ0α0β1α1β2α2=αiβiαi-1βjαjβkαk({i,j,k}∈ic,0<α0,α1,α2),
which gives the cross-boundary derivative through (2.32b) as  (2.58)si[β0α0β1α1β2α2]=(αkβkαk-1βiαiβjαj-αiβiαi-1βjαjβkαk)|βi=0={0(1<αi),-βjαjβkαk(αi=1)({i,j,k}∈ic,0<α0,α1,α2).
Assuming that (2.56b) is at most quartic, the degree of G on its right-hand side is one or zero. This gives(2.59)γ=β0β1β2(γ0β0+γ1β1+γ2β2)=γiβi2βjβk+γjβiβj2βk+γkβiβjβk2({i,j,k}∈ic),
where γ0, γ1, and γ2 are unknown coefficients. The cross-boundary derivative of the above polynomial is obtained through (2.58) as  (2.60)si[γ]=-γjβj2βk-γkβjβk2({i,j,k}∈ic).
Substituting this and (2.43) into the condition (2.56a)-(2.56b) results in  (2.61)Ĉiβjβkω|βi=0=(ĉin̂jβj+d̂in̂kβk)T(γjβj2βk+γkβjβk2)=(ĉin̂jTγj)βj3βk+(ĉin̂jTγk+d̂in̂kTγj)βj2βk2+(d̂in̂kTγk)βjβk3({i,j,k}∈ic)
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  (2.62)ω=∑i=02Wiβi2+∑{i,j,k}∈icwiβjβk,
where Wi and wi are constants. This leads to  (2.63)ω|βi=0=Wjβj2+Wkβk2+wiβjβk({i,j,k}∈ic).
Substituting the above into (2.61) and comparing the coefficients arrives at the condition(2.64a)ĈiWj=ĉin̂jTγj,(2.64b)ĈiWk=d̂in̂kTγk,(2.64c)Ĉiwi=ĉin̂jTγk+d̂in̂kTγj({i,j,k}∈ic),which is required to reproduce the normal of (2.43). Advancing the suffixes of (2.64a) by one gives ĈjWk=ĉjn̂kTγk, whose ratio to (2.64b) is(2.65)ĈjĈi=ĉjd̂i({i,j,k}∈ic).
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.


## 2.2.4. 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 Ei. If the term is quadratic or higher in βj and βk({i,j,k}∈ic), it does not affect the cross-boundary derivatives on the other edges. This observation motivates the choice of the quintic polynomial(2.66)γ=γiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βk({i,j,k}∈ic)=β0β1β2(γ0β1β2+γ1β2β0+γ2β0β1),
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)si[γ]=-γiβj2βk2({i,j,k}∈ic).
Note that the above result has the coefficient γi only. Hence the corresponding term can exclusively adjust the normal on the edge Ei. Substitution of (2.67) and (2.43) reduces the condition (2.56a) to  (2.68)Ĉiβjβkω|βi=0=(ĉin̂jβj+d̂in̂kβk)Tγiβj2βk2({i,j,k}∈ic).
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)ω|βi=0=(Wjβj2+Wkβk2+wiβjβk)(βj+βk)=Wjβj3+(Wj+wi)βj2βk+(Wk+wi)βjβk2+Wkβk3({i,j,k}∈ic).
Substituting the above into (2.68) and comparing the coefficients lead to the following conditions:(2.70)ĈiWj=ĈiWk=0,(2.71)Ĉi(Wj+wi)=ĉin̂jTγi,Ĉi(Wk+wi)=d̂in̂kTγi({i,j,k}∈ic).
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)Wi=0(i=0,1,2).
Consequently, only the terms with the coefficient wi 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)wi=1(i=0,1,2).Substitution of (2.72a)-(2.72b) into (2.62) determines the denominator polynomial as(2.73)ω=β0β1+β1β2+β2β0,
and reduces (2.71) to the simple linear problem  (2.74){ĈiĈi}=[ĉin̂jTd̂in̂kT]γi({i,j,k}∈ic).
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)γi=[ĉin̂jTd̂in̂kT]+{ĈiĈi}({i,j,k}∈ic).
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 C0 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.76a)[ĉin̂jTd̂in̂kT]+=([ĉin̂j,d̂in̂k]+)T={1D(0)[ĉin̂j,d̂in̂k][a22-a12-a12a11](D(0)≠0),1D(1)[ĉin̂j,d̂in̂k](D(0)=0,D(1)≠0),[O](otherwise),(2.76b)a11≡ĉi2,a22≡d̂i2,a12≡ĉid̂in̂jTn̂k({i,j,k}∈ic),(2.76c)D(0)=a11a22-a122,D(1)=a11+a22.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 theC0 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(2.77)δg2(diTdi)<γiTγi(i=0,1,2),
where δg is a dimensionless constant. If the above holds, the next expression on the right-hand 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 tangent-plane continuity. If the vectors ĉin̂j, d̂in̂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̂in̂j, d̂in̂k  can be linearly dependent in the following two situations. (A)
n̂j∥n̂k.The tangent vector-Tij and the cross-boundary derivative Sij at the vertex Vj (see (2.34a)-(2.34b)) become perpendicular also to the normal n̂k at the vertex Vk. The same is true if the vertices Vj and Vk are interchanged. Therefore,
(2.78)0=n̂kTTij=n̂kTSij=n̂jTTik=n̂jTSik.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 thatĉi=0. Then (2.49) yields Ĉi=d̂in̂kTSij and the first row of (2.74) gives Ĉi=0. Since d̂i is nonzero, (2.74) can hold only if n̂kTSij=0. Similarly, n̂jTSik=0 is required when d̂i=0.The above discussion concludes that tangent-plane continuity can be attained except for the following two pathological cases (As mentioned at the end of Section2.2.1, it is assumed that the C0 patch reproduces the normals at its vertices.):(2.79a)ĉi=0,n̂kTSij≠0,(2.79b)d̂i=0,n̂kTSij≠0.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 Ei is on the tangent plane at the vertex Vk. The plane, however, does not contain the cross-boundary derivative Sij at the other vertex Vj, due to (2.79a). This implies that the original curve of the edge Ei has displacement away from the tangent plane near the vertex Vj and cannot be approximated by the simple quadratic function (2.29). Interchanging Vj and Vk leads to the case of (2.79b). Hence the edge Ei 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 Vi of (2.21a) in the space of barycentric coordinates. The position vector to a point on the path can be described as(2.80)β=1i+he,e≡{e0,e1,e2}T,
where h is the distance from the vertex Vi(β=1i) in the space. The unit vector e is variable and arbitrary as long as the above β satisfies (2.15b); that is,(2.81)|e|=1,e0+e1+e2=0.
Substituting (2.80) into (2.73) yields  (2.82)ω=(1+hei)(hej)+(hej)(hek)+(hek)(1+hei)=(ej+ek)h+O(h2)({i,j,k}∈ic),
which reveals that the denominator polynomial is first order with respect to h. Similarly by (2.80),
(2.83)βjpβkq=ejpekqhp+q
holds. If the ratio(2.84a)βjpβkqωhr=ejpekqhp+q-r-1(ej+ek)+O(h)
has a positive exponent of h in the numerator, it vanishes as it approaches the vertex β=1i for h→0 irrespective of the path; that is  (note that ej+ek never vanishes in the denominator of (2.84a) due to (2.81)), (2.84b)limβ→1iβjpβkqωhr=0(1+r<p+q,{i,j,k}∈ic).The above with r=0 as well as (2.66) and (2.50) gives the limit of the correction term   (2.85)limβ→1ig=limβ→1iγiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βkω=0,
which does not depend on the path. Therefore, adopting it as the value of g at the vertex Vi to redefine (2.50) as   (2.86)g(β)={0(β0=1orβ1=1orβ2=1),γ(β)ω(β)(otherwise)
preserves the continuity. This function is also smooth: substitution of e=1i into (2.52) yields the partial derivative(2.87)∂∂βi(γω)=1ω∂γ∂βi-1ω2∂ω∂βiγ=1ω2((βiβj+βjβk+βkβi)(γiβj2βk2+2γjβiβjβk2+2γkβiβj2βk)-(βj+βk)(γiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βk))(∵(2.66)and(2.73))=γiω2βj2βk2(βiβj+βjβk+βkβi-βi(βj+βk))+γjω2βiβjβk2(2βiβj+2βjβk+2βkβi-βi(βj+βk))+γkω2βiβj2βk(2βiβj+2βjβk+2βkβi-βi(βj+βk))=γiω2βj3βk3+γjω2βiβjβk2(ω+βjβk)+γkω2βiβj2βk(ω+βjβk)(∵(2.73))=1ω2(βj3βk3γi+βiβjβk(ω+βjβk)(βkγj+βjγk))({i,j,k}∈ic),
whose limit vanishes at every vertex because of (2.84a)-(2.84b) with r=0 and the boundedness of (2.81); that is,
(2.88)limβ→1i∂∂βi(γω)=0(i=0,1,2).
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):   (2.89)Deg(1i)=limh→0g(1i+he)-g(1i)h=limh→0g(1i+he)h(∵(2.86))=limh→0γiβiβj2βk2+γjβi2βjβk2+γkβi2βj2βkωh(∵(2.66)and(2.50))=0(∵(2.84b)withr=1),
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)Deg(β)={0(β0=1orβ1=1orβ2=1),De(γ(β)ω(β))(otherwise).
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).


## 2.2.5. Algorithm

Finally, the procedure for determining theG1 patch is summarized below.(A)
Give the position and normal vectors at the three vertices of the patch.(B)
Apply the algorithm described in Section2.1 to determine the coefficients of (2.13) for the C0 patch and obtain the vectors Pi,Qi(i=0,1,2) through (2.25a)–(2.25c).(C)
Compute the constantsTij,Tik,Sij,Sik({i,j,k}∈ic) of (2.35a)-(2.35b) and the unit normals n̂i(i=0,1,2) of (2.36).(D)
After calculatingĉi, d̂i through (2.44), (2.42b), and (2.41), evaluate (2.49) to get Ĉi(i=0,1,2). If one of ĉi, d̂i is zero and Ĉi≠0, the tangent-plane continuity on the edge Ei is not achieved. Then the edge Ei 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 C0 patch yields the parametric representation xg of (2.45) for the G1 patch.


## 3. Numerical Examples

Figure3 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 C0 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 G1 patch algorithm of Section 2.2.5 yielding a completely smooth surface in Figure 3(c).Application of the interpolation to a mesh of a smooth surface.
(a)(b)(c)Next, a more general mesh of Figure4(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 G1 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 C0 patch which enables handling an arbitrary number of normals at a vertex.Application of the interpolation to a mesh of a surface with sharp features.
(a)(b)As already mentioned in Section1, 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×104 points for the cone; more than 88×104 points for the cylinder and the sphere; more than 13×104 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 G1 and C0 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 G1 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).Figure 5
Interpolation error of the Nagata patches as a function of the mesh scale.For most of the cases in Figure5, the G1 and C0 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.


## 4. Concluding Remarks

TheC0 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 C0 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 G1 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.As can be imagined from the varying nonintegral order of convergence observed in Figure5, 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.



---

*Source: 952420-2010-08-01.xml*