Construction and smoothing of triangular

Coons patches with geodesic boundary curves

R. T. Farouki,(b) N. Szafran,(a) L. Biard(a)
(a)
Laboratoire Jean Kuntzmann, Université Joseph Fourier – Grenoble, France
(b)
Department of Mechanical and Aeronautical Engineering,
University of California, Davis, CA 95616, USA

Abstract
Given three regular space curves r1 (t), r2 (t) r3 (t), t ∈ [ 0, 1 ] that define a
curvilinear triangle, we consider the problem of constructing a C 2 triangular
surface patch R(u1 , u2 , u3 ) bounded by these three curves, such that they are
geodesics of the constructed surface. Results from a prior study [6] concerned
with tensor–product patches are adapted to identify constraints on the given
curves for the existence of such geodesic–bounded triangular surface patches.
For curves satisfying these conditions, the patch is constructed by means of a
cubically–blended triangular Coons interpolation scheme. A formulation of
thin–plate spline energy in terms of barycentric coordinates with respect to a
general domain triangle is also derived, and used to optimize the smoothness
of the geodesic–bounded triangular surface patches.
Key words: geodesic curves, surface reconstruction, Coons interpolation

1. Introduction
In a prior study [6], sufficient–and–necessary conditions for the boundaries
of a four–sided analytic surface patch to be geodesic curves of the surface were
identified. A Coons interpolation scheme was then developed [7] to construct
polynomial and rational examples of such geodesic–bounded surface patches,
given boundary curves that satisfy the existence conditions. In the present
paper, these results are extended to the case of geodesic–bounded triangular
surface patches, parameterized in terms of barycentric coordinates. Also, the
familiar Cartesian expression for thin–plate spline energy is transformed to
barycentric coordinates defined on a general triangular domain T . In general,

Preprint submitted to Elsevier July 1, 2009

the barycentric form contains 21 terms, whose coefficients depend explicitly
on the geometry of T . This form is used to optimize the shape of geodesic–
bounded triangular patches with respect to residual free parameters.
The motivation for these studies comes from the problem of constructing
computer representations of free–form surfaces from positional/orientational
measurements obtained with the Morphosense — a flexible ribbon–like device
with embedded microsensors that assumes the shape of a geodesic when laid
on a smooth physical surface [15]. By orienting the Morphosense in different
directions on the surface, it may be divided into a collection of rectangular
and triangular patches. Note that constructing triangular geodesic–bounded
patches is actually a more fundamental problem, since a closed surface of
genus zero cannot be entirely covered by four–sided patches.
The remainder of this paper will be organized as follows. After reviewing
some basic facts concerning differential geometry and barycentric coordinates
on triangular domains in Section 2, the problem of constructing triangular
surface patches with geodesics as boundary curves is described in Section 3,
and the methodology for solving it is presented in Section 4. The properties
of the constructed surface patches are then discussed, in the context of the
Gauss–Bonnet theorem, in Section 5. Section 6 develops the thin–plate spline
energy in terms of barycentric coordinates on general triangle domains, while
Section 7 presents computed examples, with shape optimized in terms of this
energy. Finally, Section 8 summarizes and assesses key results of the paper.

2. Preliminaries

For linearly–independent unit vectors u and v, n

and a unit vector n such that n ⊥ u and n ⊥ v, let
(u, v)n denote the oriented angle between u, v in v
the sense of n. Specifically, the angle A = (u, v)n
is defined (see Fig. 1) by u

sin A = det(u, v, n), cos A = hu, vi . (1) Fig. 1. Angle

measurement.
We assume the reader is familiar with the elementary differential geometry
of curves and surfaces [3, 16] — see also [6] for pertinent background.
For a space curve r(t) let (e(t), n(t), b(t)) be the Serret–Frenet frame, and

2
k(t) and τ (t) the curvature and torsion, defined at non–inflectional points by

r0 (t) r0 (t) × r00 (t)

e(t) = , b(t) = , n(t) = b(t) × e(t) , (2)
||r0 (t)|| ||r0 (t) × r00 (t)||

||r0 (t) × r00 (t)|| det(r0 (t), r00 (t), r000 (t))
k(t) = and τ (t) = . (3)
||r0 (t)||3 ||r0 (t) × r00 (t)||2
For a regular curve, satisfying kr0 (t)k 6= 0 for all t, the tangent e(t) is defined
at every point. At inflection points (where r0 (t), r00 (t) are linearly dependent),
the curvature k(t) is zero, and the principal normal n(t), binormal b(t), and
torsion τ (t) are undefined. Henceforth, we specifically exclude the possibility
of “pathological” inflections — i.e., points where the curvature vanishes, and
the left and right limits n− and n+ of the principal normal are neither parallel
nor anti–parallel, since no solution to the geodesic interpolation problem can
be found in the presence of such points.

The barycentric coordinates u = (u1 , u2 , u3 ) of a point p = (x, y) with respect

to a domain triangle T with vertices pi = (xi , yi ), i = 1, 2, 3 are given by
    
u1 y2 − y3 x3 − x2 x2 y3 − x3 y2 x
u2  = 1  y3 − y1 x1 − x3 x3 y1 − x1 y3  y  , (4)
∆
u3 y1 − y2 x2 − x1 x1 y2 − x2 y1 1

where ¯ ¯
¯ 1 1 1 ¯
¯ ¯
∆ = ¯¯ x1 x2 x3 ¯¯ . (5)
¯ y1 y2 y3 ¯
For distinct points p, q the directional derivative along the vector V = q − p
is defined [4] by
X 3
∂
DV = (ui (q) − ui (p)) .
i=1
∂ui
Hence, the derivatives along directions parallel to each side si of T , defined
by ui = 0, is
∂ ∂
Dsi = − , i = 1, 2, 3 (6)
∂ui−1 ∂ui+1
where all indices are reduced modulo 3 to the set {1, 2, 3}.

3
 p1
u1 = 1

Ds3 Ds2

Dn 3
Dn 2

Dn 1 u1 = 0

p2 p3
Ds1
u3 = 0 u2 = 1 u3 = 1 u2 = 0

Fig. 2. Barycentric coordinates relative to the triangle (p1 , p2 , p3 ).

We then define (see Fig. 2) a derivative Dni in the inward direction along
each side si by

Dni = αi−1 (u) Dsi+1 − αi+1 (u) Dsi−1 , i = 1, 2, 3 (7)

where the three “blending functions” defined [8] by

αi (u) = u2i (3 − 2 ui + 6 ui−1 ui+1 ) , i = 1, 2, 3

satisfy
3
X
αi (u) ≡ 1 .
i=1

Note also that, since αi−1 (u) + αi+1 (u) ≡ 1 along the side ui = 0 of T ,
the derivative Dni coincides with the derivative Dsi+1 at the point pi−1 and
−Dsi−1 at the point pi+1 , along the side si .

3. The three-geodesic interpolation problem

Consider, as shown in Fig. 3, three regular space curves r1 (t), r2 (t), r3 (t)
with t ∈ [ 0, 1 ] such that

4
 ri+1 (1) = pi = ri+2 (0) (8)
p1

for i = 1, 2, 3, where all indices are

reduced modulo 3 (so that 0 → 3, r 3 (t)
r 2 (t)
4 → 1, 5 → 2). These curves
are assumed “sufficiently smooth,”
and their derivatives at the corner
points pi , i = 1, 2, 3 are assumed r 1 (t)
p3
to be linearly independent, so as p2

to define the tangent plane Πi for

the interpolating surface R at pi . Fig. 3. Triangular patch boundaries.

3.1. The geodesic interpolation problem

Our intent is to construct an oriented surface patch that interpolates the
specified curves as its boundaries, in such a manner that they are geodesic
paths of the constructed surface. In terms of the barycentric coordinates u =
(u1 , u2 , u3 ) relative to the triangle with vertices p1 , p2 , p3 the interpolating
surface R(u) must have the following properties.

• Interpolation property:

R(0, 1 − t, t) = r1 (t) ,
R(t, 0, 1 − t) = r2 (t) , t ∈ [ 0, 1 ] . (9)
R( 1 − t, t, 0) = r3 (t) ,

• Geodesic property: the principal normal n( t) at each non–inflectional

point of the curves ri (t) must be normal to the surface R(u).

Let (ei (t), ni (t), bi (t)) denote the Serret–Frenet frame and ki (t), τi (t) the
curvature and torsion of each of the curves ri (t), as defined by (2) and (3).
Let r(t) be the “concatenation” of the three boundary curves ri (t), defined by
r(t) = ri (t − i + 1) for t ∈ [i − 1, i] and i = 1, 2, 3, and consider the principal
normal n(t) of this concatenated curve r(t), defined for t ∈ [ 0, 3 ] − {0, 1, 2}.
This normal n(t) is simply the concatenation of the principal normals ni (t) of
the individual boundary curves ri (t). Note that the domain of the parameter
t is considered to [ 0, 3 ] modulo 3.
Now let N(u) be the unit normal to the interpolating surface R(u), and
let N(t) for t ∈ [ 0, 3 ] be its restriction to the boundary r(t) of the surface

5
patch. For a continuous normal function N(t), we define the crossing angles
Ai at each corner pi by
¡ ¢
Ai = ei+1 (1), ei−1 (0) N(pi ) , i = 1, 2, 3 . (10)

Assuming that the principal normals of the two boundary curves meeting at
each corner agree (modulo sign) with the unit surface normal N, we consider
the values σiL , σiR ∈ {−1, +1} defined by

N(0) = N(p1 ) = σ2R n2 (1) = σ3L n3 (0) ,

N(1) = N(p2 ) = σ3R n3 (1) = σ1L n1 (0) , (11)
N(2) = N(p3 ) = σ1R n1 (1) = σ2L n2 (0) .

p1 N n2
e2
b2
N
b3 = h3
Fig. 4 shows the Frenet frames h2
(ei , ni , bi ) and the Darboux r2 (t)

frames (ei , hi , N), with hi (t) = e3 n3 N

b1 = h1

N(t) × ei (t), along the patch r3 (t)

e1

boundary curves i = 1, 2, 3. r1 (t)

p2 p3
n1

Fig. 4. Vectors along patch boundaries.

3.2. Existence conditions for an interpolating surface
From Proposition 7 in [6], a surface interpolating the boundary curves as
geodesics can be found if and only if they satisfy the following constraints:

(C1) osculating constraint: the principal normals of the boundary curves

that meet at each corner pi must agree modulo sign. This implies that,
at each corner pi , the boundary curves meeting there have osculating
planes orthogonal to the surface tangent plane Πi ;

(C2) global normal orientation constraint: a continuous unit vector

function N(t) must exist, such that N(t) = ± n(t) for all t ∈ R ;

(C3) corner geodesic crossing constraint: for the crossing angles Ai

and values σiL , σiR ∈ {−1, +1} defined by (10) and (11), the boundary

6
 curves ri (t) must satisfy the following curvature and torsion constraints
at each corner pi :

[ σi+1,R ki+1 (1) − σi−1,L ki−1 (0) ] cos Ai

+ [ τi+1 (1) + τi−1 (0) ] sin Ai = 0 . (12)

4. Three–sided geodesic interpolation

Given three space curves free of inflection points, that satisfy conditions
(C1)–(C3), we wish to construct a triangular patch R(u1 , u2 , u3 ) interpolating
the curves as geodesic boundaries. A cubically–blended Coons interpolation
scheme [8] will be invoked to construct the surface (see also [6, 9, 11, 15] for
related prior work on interpolation of geodesic boundary curves).

4.1. Triangular Geodesic Coons interpolation

p1

r2 (t)
For a triangular patch, the T3(t)

cubic Coons interpolation r3 (t)

T2(t)
scheme requires a tangent
T1(t)
vector Ti (t) to be defined
along each of the boundary r1 (t)
curves ri (t) — see Fig. 5. p2 p3

Fig. 5. Tangent vectors Ti (t) along patch

boundary curves.
In order to satisfy the geodesic constraint, the tangent vectors Ti (t) must lie
in the surface tangent plane at each point of ri (t), so these tangent vectors
can, in general, be expressed in the form
£ ¤
Ti (t) = di (t) cos(αi (t)) ei (t) + sin(αi (t)) hi (t) , i = 1, 2, 3 , (13)

where the angle functions αi (t) specify the inclination of the vectors Ti (t)
relative to the boundary curve tangents ei (t). As in [6], we will now show
that the vector fields Ti (t) must satisfy the following two conditions:
• Interpolation of the corner derivatives (see Section 4.2):

Ti (0) = − r0i−1 (1) and Ti (1) = r0i+1 (0) , i = 1, 2, 3 (14)

7
 — i..e, the tangent vectors Ti must coincide with the derivatives of the
two curves ri−1 (t), ri+1 (t) meeting the curve ri (t) at its extremities.
• Twist crossing constraints at the corners (see Section 4.3):
The derivatives of the tangent vectors Ti (t) must satisfy the following
relations, which define the twist vectors Wi at the patch corners.

T0i−1 (0) = − T0i+1 (1) =: Wi , i = 1, 2, 3 . (15)

Consider now the C 1 -triangular Coons patch R(u), as developed in [8]:

3
X
R(u) = αi (u) Ri (u) , (16)
i=1

with

Ri (u) = ri−1 (ui+1 ) + ui−1 [ Ti−1 (ui+1 ) − Ti−1 (0) ]

+ ri+1 (1 − ui−1 ) + ui+1 [ Ti+1 (1 − ui−1 ) − Ti+1 (1) ]
− pi − ui−1 ui+1 Wi . (17)

A straightforward computation yields

R(u)|ui =0 = ri (ui−1 ) , i = 1, 2, 3 , (18)

and
¯
Dni R(u)|ui =0 = [ αi−1 (u) Dsi+1 − αi+1 (u) Dsi−1 ] R(u)¯ui =0
· ¸X ¯
∂ ∂ ∂
3 ¯
¯
= − αi+1 (u) − αi−1 (u) αi (u) Ri (u)¯
∂ui ∂ui+1 ∂ui−1 i=1 ¯
ui =0
= Ti (ui−1 )
− αi+1 (u)|ui =0 [ r0i−1 (1) + Ti (0) + ui−1 (T0i−1 (1) + Wi+1 ) ]
+ αi−1 (u)|ui =0 [ r0i+1 (0) − Ti (1) + ui+1 (T0i+1 (0) − Wi−1 ) ] . (19)

Assuming that the three vector fields Ti (t) satisfy conditions (14)–(15),
the patch R(u) will interpolate the boundary curves ri (t) and the tangent
vectors Ti (t) precisely:

R(u)|ui =0 = ri (u) , Dni R(u)|ui =0 = Ti (u) , i = 1, 2, 3 . (20)

8
which proves – with our hypothesis – that the patch R(u) interpolates the
boundary curves ri (t) as geodesics.
Consider now constraints (14) and (15). We will establish conditions on
the vector fields Ti (t) in order to satisfy these constraints. In particular, we
show that the twist constraints (15) can only be satisfied if the conditions
(C1)–(C3) are satisfied by the curves ri (t).

4.2. Interpolation of corner derivatives

ri-1
(t) ri+1(t)
Ti r'i+1
Ti (0)= -r'i-1 (1)
(1)= (0)

ei+1
(0)

hi
ei
(0)
(0)

ri A i-1
N(pi+1)
(t)

hi-1 N(pi-1) ei
A i+1 (1)
(1)

ei-1 (1)

Fig. 6. Vectors along the curve ri (t).

To satisfy the constraints (14) on the vector field Ti (t), defined by (13), we
introduce the following scalar coefficients diL , diR , αiL , αiR (see Fig. 6) for
i = 1, 2, 3:

diL := di (0) = ||r0i−1 (1)|| , diR := di (1) = ||r0i+1 (0)|| ,

(21)
αiL := αi (0) = π − Ai+1 , αiR := αi (1) = Ai−1 .

Then, as in [6], we consider the functions

αi (t) = αiL H0 (t) + αi1 H1 (t) + αi2 H2 (t) + αiR H3 (t) , (22)

di (t) = diL H0 (t) + di1 H1 (t) + di2 H2 (t) + diR H3 (t) , (23)
where Hi (t) are the cubic Hermite polynomials, and the parameters αi1 , αi2 ,
di1 , di2 will be specified in Section 4.3.

9
4.3. Twist crossing constraints at corners pi

r'i+1 (1)
ei+1
(1)

Ai N(pi )
ei-1 (0)

r'i-1
(0)

hi+1 (1) hi-1

(0)

ri-1(t)
ri+1
(t)

Fig. 7. Vectors at corner pi .

At corner pi , the vectors ei+1 (1), hi+1 (1), ei−1 (0), hi−1 (0) are coplanar and
related (see Fig. 7) as follows:

ei+1 (1) = cos Ai ei−1 (0) − sin Ai hi−1 (0) ,

ei−1 (0) = cos Ai ei+1 (1) + sin Ai hi+1 (1) .

Hence (noting from Section 3 that Ai 6= 0, π) we deduce that

cos Ai ei−1 (0) − ei+1 (1) ei−1 (0) − cos Ai ei+1 (1)
hi−1 (0) = , hi+1 (1) =.
sin Ai sin Ai
(24)
Then using the definitions (13), (22), (23), relations (11), (21) and the Serret-
Frenet relations as in [6], with

bi+1 (1) = −σi+1,R hi+1 (1) , bi−1 (0) = −σi−1,L hi−1 (0) ,

we express each of the derivatives occurring in (15) with respect to the frame
(ei+1 (1), ei−1 (0), N(pi )) as
h cos Ai i
T0i−1 (0) = − d0i−1 (0) + ||r0i+1 (1)|| 0
αi−1 (0) ei+1 (1)
sin Ai
1
0
− ||r0i+1 (1)|| αi−1 (0) ei−1 (0)
sin Ai
h i
+ ||r0i+1 (1)|| ||r0i−1 (0)|| − cos Ai σi−1,L ki−1 (0) + sin Ai τi−1 (0) N(pi ) ,

10
 1
T0i+1 (1) = − ||r0i−1 (0)|| αi+10
(1) ei+1 (1)
sin Ai
h cos Ai i
0 0 0
+ di+1 (1) + ||ri−1 (0)|| αi+1 (1) ei−1 (0)
sin Ai
h i
+ ||r0i+1 (1)|| ||r0i−1 (0)|| cos Ai σi+1,R ki+1 (1) + sin Ai τi+1 (1) N(pi ) .

Thus, the twist crossing constraint (15) at corner pi , together with definitions
(21)–(23), yields the following system of equations


 cos Ai 1

 − di−1,1 + di−1,L αi−1,1 = di+1,R αi+1,2 ,

 sin Ai sin Ai


 1 cos Ai
− di−1,L αi−1,1 = − di+1,2 − di+1,R αi+1,2 , (25)
 sin Ai sin Ai



 − cos Ai σi−1,L ki−1 (0) + sin Ai τi−1 (0) =



 − cos Ai σi+1,R ki+1 (1) − sin Ai τi+1 (1) .

As in [6], identifying components of the twist crossing constraint (15) in

the direction of N(pi ) at corner pi yields the geodesic crossing relation (12)
at that point, which is satisfied through the assumption that the boundary
curves obey conditions (C1)–(C3) in Section 3.2. Thus, multiplying equations
(25) by sin Ai yields the following linear equations in the four unknowns di−1,1 ,
αi−1,1 , di+1,2 , αi+1,2 :
½
− sin Ai di−1,1 + cos Ai di−1,L αi−1,1 − di+1,R αi+1,2 = 0 ,
(26)
− di−1,L αi−1,1 + sin Ai di+1,2 + cos Ai di+1,R αi+1,2 = 0 .

This linear system is of rank 2, and if Si denotes the set of its solutions, we
have dim(Si )=2 (see Section 7 for examples illustrating the influence of the
free parameters in Si , i = 1, 2, 3). Thus, we are able to define a twist vector
Wi at the patch corner pi , and we can now state the following result.

Proposition 1. Given three regular space curves, that satisfy the conditions
(C1)–(C3) specified in Section 3.2, an oriented triangular surface patch R(u)
exists, that interpolates these curves as boundaries in such a way that they are
geodesics of the surface. Conversely, when any of the conditions (C1)–(C3)
is not satisfied, such an interpolating surface can not be constructed.

11
5. Gauss–Bonnet theorem
The well–known Gauss–Bonnet theorem relates the integral of Gaussian
curvature over a region of a surface to the integral of the geodesic curvature
along the boundary of that region. This theorem is of particular relevance
to the geodesic–bounded surface patches constructed herein.
Theorem 1. Global Gauss-Bonnet Theorem [3, 16].
Let R ⊂ S be a regular simple region of an oriented surface S and Γ1 , . . . , Γk
be closed, simple, piecewise–regular curves forming its boundary ∂R. Suppose
that Γ1 , . . . , Γk are all positively oriented, and let A1 , . . . , Ak be the external
angles at the junctures of the curves Γ1 , . . . , Γk , as defined by (10). Then
k Z
X ZZ k
X
ki (s) ds + K dσ + Ai = 2π , (27)
i=1 Γi R i=1

where ki (s) is the geodesic curvature of the curve Γi , s is arc length along it,
and K is the Gaussian curvature of the surface S with area element dσ.
Since the geodesic curvature vanishes along any geodesic curve, it follows
from the above theorem that if we construct a triangular surface patch R(u)
as described in Proposition 1, so that its boundary curves are geodesics, the
Gaussian curvature K of this surface must satisfy
ZZ X3
K dσ = 2π − Ai , (28)
R i=1

where A1 , A2 , A3 are the corner angles at the vertices p1 , p2 , p3 .

Fig. 8. Gauss–Bonnet triangle.

12
 For example, givenPthree regular space curves ri (t) satisfying conditions
(C1)–(C3), such that 3i=1 Ai = 2π, as in Fig. 8, any interpolating surface
R(u) constructed as described in Proposition 1 must satisfy
ZZ
K dσ = 0 , (29)
R

which indicates that the interpolating surface must contain both elliptic and
hyperbolic points (where K > 0 and K < 0) — see Fig. 9.

Fig. 9. Two views of a surface interpolating the Gauss–Bonnet triangle of

Fig. 8 as geodesics. The free parameters dij and αij are set equal to zero.
Consider the problem of smoothing the interpolating surface. From (28),
we see that the integral of the Gaussian curvature will be the same for all
geodesic–bounded interpolating surfaces R(u). Thus, one natural choice is
to minimize the variation of Gaussian curvature over the surface. Another
choice is based on minimizing the thin–plate spline energy, adapted to the
context of triangular patches as described below.

6. Thin–plate spline energy for triangular patches

For a bivariate function f (x, y) the thin–plate spline energy per unit area
is defined by
2 2 2
E = fxx + 2fxy + fyy . (30)

13
This expression is clearly invariant under any translation, since it depends
only on the second derivatives of f . One can easily see that it is also invariant
under any rotation, specified by
· 0¸ · ¸· ¸
x cos θ sin θ x
= .
y0 − sin θ cos θ y

Applying the derivative operators

∂ ∂x0 ∂ ∂y 0 ∂ ∂ ∂
= 0
+ 0
= cos θ 0 − sin θ 0 ,
∂x ∂x ∂x ∂x ∂y ∂x ∂y
∂ ∂x0 ∂ ∂y 0 ∂ ∂ ∂
= + = sin θ 0 + cos θ 0 ,
∂y ∂y ∂x0 ∂y ∂y 0 ∂x ∂y
twice to f , one can verify that

fx20 x0 + 2fx20 y0 + fy20 y0 ≡ fxx

2 2
+ 2fxy 2
+ fyy .

Hence, for Cartesian coordinates, the thin–plate spline energy per unit area
(30) is a “universal” expression.
Consider now the formulation of the thin–plate spline energy in terms of
the barycentric coordinates (4) with respect to a domain triangle T with the
vertices pi = (xi , yi ) for i = 1, 2, 3. To transform derivatives with respect to
(x, y) into derivatives with respect to (u1 , u2 , u3 ) we use the relations

∂ ∂u1 ∂ ∂u2 ∂ ∂u3 ∂

= + +
∂x ∂x ∂u1 ∂x ∂u2 ∂x ∂u3
· ¸
1 ∂ ∂ ∂
= (y2 − y3 ) + (y3 − y1 ) + (y1 − y2 ) ,
∆ ∂u1 ∂u2 ∂u3
∂ ∂u1 ∂ ∂u2 ∂ ∂u3 ∂
= + +
∂y ∂y ∂u1 ∂y ∂u2 ∂y ∂u3
· ¸
1 ∂ ∂ ∂
= (x3 − x2 ) + (x1 − x3 ) + (x2 − x1 ) . (31)
∆ ∂u1 ∂u2 ∂u3

Clearly, these derivative operators depend explicitly on the chosen vertices

of the domain triangle T . Consequently, there is no “universal” barycentric
formulation of the expression (30). It was noted above that (30) is invariant
under transformations of the Cartesian coordinates that correspond to rigid
motions. However, mappings between barycentric coordinates with respect to

14
different domain triangles define general affine transformations, not just rigid
motions. Hence, each domain triangle T has its own associated expression
for the thin–plate spline energy, dependent on the vertices of T .
We denote the second derivatives of f with respect to (u1 , u2 , u3 ) by

∂ 2f
fjk = , 1 ≤ j, k ≤ 3 .
∂uj uk

The expression for E in barycentric coordinates with respect to an arbitrary

domain triangle T contains 21 terms ¡ ¢ — the squares of the 6 second derivatives
f11 , f22 , f33 , f12 , f23 , f31 and the 62 = 15 cross terms arising from their pair–
wise products. The coefficients of these 21 terms will depend on the chosen
vertices of T . Using the partial derivatives (31) and writing xjk = xj − xk
and yjk = yj − xk for 1 ≤ j, k ≤ 3 we obtain

∆2 fxx = 2
y23 f11 + y312 2
f22 + y12 f33
+ 2 (y23 y31 f12 + y31 y12 f23 + y12 y23 f31 ) ,
2
∆ fxy = x32 y23 f11 + x13 y31 f22 + x21 y12 f33
+ (x13 y23 + x32 y31 )f12 + (x21 y31 + x13 y12 )f23 + (x21 y23 + x32 y12 )f31 ,
2
∆ fyy = x232 f11 + x213 f22 + x221 f33
+ 2 (x32 x13 f12 + x13 x21 f23 + x21 x32 f31 ) ,

Defining the sides of the domain triangle T by s1 = p3 − p2 , s2 = p1 − p3 ,

s3 = p2 − p1 , one can then verify that expression (30) becomes

∆4 E = |s1 |4 f11
2
+ |s2 |4 f22
2
+ |s3 |4 f33
2

+ 2 [ |s1 |2 |s2 |2 + (s1 · s2 )2 ] f12

2

+ 2 [ |s2 |2 |s3 |2 + (s2 · s3 )2 ] f23

2

+ 2 [ |s3 |2 |s1 |2 + (s3 · s1 )2 ] f31

2

+ 2 (s1 · s2 )2 f11 f22 + 2 (s2 · s3 )2 f22 f33 + 2 (s3 · s1 )2 f33 f11

+ 4 f11 [ |s1 |2 (s1 · s2 ) f12 + (s1 · s2 )(s3 · s1 ) f23 + |s1 |2 (s3 · s1 ) f31 ]
+ 4 f22 [ |s2 |2 (s1 · s2 ) f12 + |s2 |2 (s2 · s3 ) f23 + (s1 · s2 )(s3 · s1 ) f31 ]
+ 4 f33 [ (s2 · s3 )(s3 · s1 ) f12 + |s3 |2 (s2 · s3 ) f23 + |s3 |2 (s3 · s1 ) f31 ]
+ 4 [ |s2 |2 (s3 · s1 ) + (s1 · s2 )(s2 · s3 ) ] f12 f23
+ 4 [ |s3 |2 (s1 · s2 ) + (s2 · s3 )(s3 · s1 ) ] f23 f31
+ 4 [ |s1 |2 (s2 · s3 ) + (s3 · s1 )(s1 · s2 ) ] f31 f12 . (32)

15
If `1 = |s1 |, `2 = |s2 |, `3 = |s3 | are the lengths of the domain triangle sides,
and θ1 , θ2 , θ3 are its interior angles at the vertices p1 , p2 , p3 so that

s1 · s2 = − `1 `2 cos θ3 , s2 · s3 = − `2 `3 cos θ1 , s3 · s1 = − `3 `1 cos θ2 ,

the expression (32) can also be written as

∆4 E = `41 f11
2
+ `42 f22
2
+ `43 f33
2

+ 2 `21 `22 (1 + cos2 θ3 ) f12

2
+ 2 `22 `23 (1 + cos2 θ1 ) f23
2
+ 2 `23 `21 (1 + cos2 θ2 ) f31
2

+ 2 `21 `22 cos2 θ3 f11 f22 + 2 `22 `23 cos2 θ1 f22 f33 + 2 `23 `21 cos2 θ2 f33 f11
+ 4 `21 f11 (`2 `3 cos θ2 cos θ3 f23 − `1 `2 cos θ3 f12 − `3 `1 cos θ2 f31 )
+ 4 `22 f22 (`3 `1 cos θ3 cos θ1 f31 − `1 `2 cos θ3 f12 − `2 `3 cos θ1 f23 )
+ 4 `23 f33 (`1 `2 cos θ1 cos θ2 f12 − `2 `3 cos θ1 f23 − `3 `1 cos θ2 f31 )
+ 4 `21 `2 `3 (cos θ2 cos θ3 − cos θ1 ) f31 f12
+ 4 `1 `22 `3 (cos θ3 cos θ1 − cos θ2 ) f12 f23
+ 4 `1 `2 `23 (cos θ1 cos θ2 − cos θ3 ) f23 f31 .

An adaptation of the thin–plate spline energy (30) is often invoked in

the problem of “smoothing” a rectangular (tensor–product) surface patch
R(u, v) for (u, v) ∈ [ 0, 1 ] × [ 0, 1 ]. Specifically, the expression

E = |Ruu |2 + 2 |Ruv |2 + |Rvv |2

is employed in lieu of (30). The equivalent expression for a triangular surface

R(u1 , u2 , u3 ) parameterized by barycentric coordinates on a domain triangle
T is obtained by replacing f11 , f12 , etc, by R11 , R12 , etc, with the cross
terms f11 f22 , f12 f23 , etc, replaced by the dot products R11 · R22 , R12 · R23 ,
etc. Here Rjk = ∂ 2 R/∂uj ∂uk for 1 ≤ j, k ≤ 3. The three patch corner points
p1 = R(1, 0, 0), p2 = R(0, 1, 0), p3 = R(0, 0, 1) are a natural choice for the
vertices of the barycentric coordinates domain triangle T .

7. Smoothing methods and computed examples

We seek optimal values for the parameters dij , αij in the set S = ∪3i=1 Si
which provide smooth interpolating surfaces. For this purpose, we propose
to minimize one of the following functionals.

16
Criterion A — minimization of the thin plate spline energy E expressed in
terms of barycentric coordinates relative to the reference domain triangle T ,
as developed Section 6: Z
min E du .
S T

Criterion B: — minimization of the expression

Z ³ X ´
min |D2si sj R(u)|2 du ,
S T 1≤i,j≤3

where D2si sj = Dsi ◦ Dsj and T is the reference domain triangle.

Fig. 10 illustrates the smoothing of geodesic–bounded triangular patches
using these two measures, for the case of the boundary curves shown in Fig. 8.

(1) (2)

(3) Criterion A (4) Criterion B

(1) (2)

17
 (3) Criterion A (4) Criterion B

Fig. 10. Two differents views of the Gauss–Bonnet example of Fig. 8,

before and after smoothing: (1) the three boundary curves; (2) the initial
reconstructed surface with free parameters dij equal to zero; (3) the
reconstructed surface after smoothing according to criterion A; and (4) the
reconstructed surface after smoothing according to criterion B. The
optimization is carried out with respect to the parameters dij — the
parameters αij are evaluated from (26). For criterion A and criterion B,
the smoothness measure is reduced by a factor of approximately 2.

Figs. 11–13 illustrate triangular patches with boundary curves defined by

geodesics on some simple quadric surfaces — a sphere, cylinder, and cone.

Fig. 11. Three geodesics on a sphere and the reconstructed surface.

18
Fig. 12. Three geodesics on a cylinder and the reconstructed surface.

(a) (b)

(c)

19
 Fig. 13. Geodesic triangle on a cone: (a) the cone and the three geodesic
curves; (b) the cone, the three geodesic curves and the interpolating surface;
and (c) the three geodesic curves and the interpolating surface alone.

8. Conclusion
Given three space curves that define a curvilinear triangle, a method has
been presented for constructing a triangular surface patch bounded by these
curves, such that they are geodesic curves of the constructed surface. For the
construction to be feasible, the given boundary curves must satisfy certain a
priori compatibility conditions. In order to optimize the smoothness of the
geodesic–bounded triangular patches, a formulation of the thin–plate spline
energy in terms of barycentric coordinates with respect to a general domain
triangle was also derived. This formulation is more complicated (containing
21 terms) than the familiar Cartesian expression, and its coefficients depend
explicitly on the shape of the domain triangle.
The availability of triangular geodesic–bounded patches, in conjunction
with earlier formulations for rectangular geodesic–bounded patches, allows
one to reconstruct surfaces from general networks of geodesic paths on them,
as physically determined by the Morphosense device.

Acknowledgment. Part of this work was accomplished during the visit

of the third author to the Department of Mechanical & Aeronautical Engi-
neering, University of California, Davis.

References
[1] Coons, S., Surfaces for Computer Aided Design, Technical Report,
M.I.T., 1964 (available as AD 663 504 from National Technical Information
Service, Springfield, VA 22161).
[2] Coons, S., Surface patches and B-spline curves, R. Barnhill and R. Riesen-
feld, editors, Computer Aided Geometric Design, Academic Press, 1974.
[3] Do Carmo, M. P., Differential Geometry of Curves and Surfaces,
Prentice–Hall, 1976.
[4] Farin, G., Triangular Bernstein–Bézier patches, Computer Aided Geo-
metric Design 3, 83–127, 2009.

20
[5] Farin, G., Curves and Surfaces for CAGD, 5th Edition, Academic Press,
2002.

[6] Farouki R. T., N. Szafran, L. Biard, Existence conditions for Coons

patches interpolating geodesic boundary curves, Computer Aided Geo-
metric Design, 2009, in press, doi:10.1016/[Link].2009.01.003

[7] Farouki R. T., N. Szafran, L. Biard, Construction of Bézier surface

patches with Bézier curves as geodesic boundaries, Computer Aided De-
sign, to appear, 2009.

[8] Gregory J. A., P. Charrot, A C 1 triangular interpolation patch for

computer–aided geometric design, Computer Graphics and Image Pro-
cessing 13, 80–87, 1980.

[9] Paluszny, M., Cubic polynomial patches through geodesics, Computer

Aided Design 40, 56-61, 2008.

[10] Peters, J., Local smooth surface interpolation: A classification, Com-

puter Aided Geometric Design 7, 191-195 1990.

[11] Sánchez–Reyes J., R. Dorado, Constrained design of polynomial surfaces

from geodesic curves, Computer Aided Design 40, 49-55, 2008.

[12] Sarraga, R. F., G1 interpolation of generally unrestricted cubic Bézier

curves, Computer Aided Geometric Design 4, 23-39, 1987.

[13] Shirman, L. A., Séquin, C. H., Local surface interpolation with Bézier
patches, Computer Aided Geometric Design 4, 279-295, 1987.

[14] Sprynski, N., Reconstruction de courbes et surfaces à partir de données

tangentielles, Laboratoires CEA/LETI, LJK, Thèse de l’université Joseph
Fourier, Grenoble; 5 Juillet 2007.

[15] Sprynski N., N. Szafran, B. Lacolle, L. Biard, Surface reconstruction via

geodesic interpolation, Computer Aided Design 40, 480-492, 2008.

[16] Struik, D. J., Lectures on Classical Differential Geometry, Dover

(reprint), 1988.

21