13.472J/1.128J/2.158J/16.940J
COMPUTATIONAL GEOMETRY
Lecture 20
Dr. K. H. Ko
Prof. N. M. Patrikalakis
Copyright c≤2003 Massachusetts Institute of Technology
Contents
20 Advanced topics in di?erential geometry 2
20.1 Geodesics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
20.1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
20.1.2 De?nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
20.1.3 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
20.1.4 Two-point boundary value problem . . . . . . . . . . . . . . . . . . . . . . . . . 5
20.1.5 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
20.2 Developable surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
20.2.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
20.2.2 De?nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
20.2.3 Developable surface in terms of B′ezier surface . . . . . . . . . . . . . . . . . . . 12
20.2.4 Development of developable surface (?attening) . . . . . . . . . . . . . . . . . . 13
20.3 Umbilics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
20.3.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
20.3.2 De?nition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
20.3.3 Computation of umbilical points . . . . . . . . . . . . . . . . . . . . . . . . . . 15
20.3.4 Classi?cation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
20.3.5 Characteristic lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
20.4 Parabolic, ridge and sub-parabolic points . . . . . . . . . . . . . . . . . . . . . . . . . 21
20.4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
20.4.2 Focal surfaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
20.4.3 Parabolic points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
20.4.4 Ridge points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
20.4.5 Sub-parabolic points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
Bibliography 25
Reading in the Textbook
? Geodesics : Chapter 10, pp.265 - 291
? Umbilics : Chapter 9, pp.231 - 264
1
Lecture 20
Advanced topics in di?erential
geometry
20.1 Geodesics
In this section we study the computation of shortest path between two points on free-form
surfaces [14, 11].
20.1.1 Motivation
? ship design
? robot motion planning
? terrain navigation
installation of underwater cable
?
20.1.2 De?nition
? t: Unit tangent vector of C at P
n: Unit normal vector of C at P?
N: Unit surface normal vector of S at P?
? u: Unit vector perpendicular to t in the tangent plane de?ned by N × t.
2
t
N
u
n
k
k
g
k
n
S
P
C
Figure 20.1: De?nition of geodesic curvature.
? We can decompose the curvature vector k of C into N component k
n
, which is called
normal curvature vector, and u component k
g
, which is called geodesic curvature vector
k = k
n
+ k
g
= ?ρ
n
N + ρ
g
u (20.1)
Here ρ
n
and ρ
g
are the normal and geodesic curvatures, respectively and de?ned as
follows:
ρ
n
= ?k N (20.2)·
ρ
g
= k u (20.3)·
? Consequently,
dt
ρ
g
=
ds
· (N × t) (20.4)
? Geodesic paths are sometimes de?ned as shortest path between points on a surface,
however this is not always a satisfactory de?nition.
De?nition: Geodesics are curves of zero geodesic curvature [24].
20.1.3 Governing equations
? The unit tangent vector of the curve C on the surface r is given by
dr(u(s), v(s)) du dv
t = = r
u
+ r
v
(20.5)
ds ds ds
? Hence using chain rules
? ?
2
? ?
2
dt du du dv dv d
2
u d
2
v
= r
uu
+ 2r
uv
+ r
vv
+ r
u
+ r
v
(20.6)
ds ds ds ds ds ds
2
ds
2
3
? Consider that the surface normal N has the direction of normal of the geodesic line ±n
n r
u
= 0, n r
v
= 0. (20.7)· ·
dt
? Since kn =
ds
, equation (20.7) can be rewritten as
dt dt
r
u
= 0, r
v
= 0 (20.8)
ds
·
ds
·
? By substituting equation (20.6) into equations (20.8) we obtain
? ?
2
? ?
2
du du dv dv d
2
u d
2
v
(r
uu
· r
u
) + 2(r
uv
· r
u
) + (r
vv
· r
u
) + E + F = 0 (20.9)
ds ds ds
? ?
2
du du dv
ds ds
2
ds
2
? ?
2
dv d
2
u d
2
v
(r
uu
· r
v
) + 2(r
uv
· r
v
) + (r
vv
· r
v
) + F + G = 0 (20.10)
ds ds ds ds ds
2
ds
2
v
By eliminating
d
2
from equation (20.9) using equation (20.10), and eliminating
d
2
u
from
ds
2
ds
2
?
equation (20.10) using equation (20.9) and employing the Christo?el symbols, we obtain
? ? ?
2
d
2
u du
?
2
du dv dv
+ ?
1
+ 2?
1
+ ?
1
= 0 (20.11)
ds
2
11 22
ds
?
d
2
v du
?
2
12
ds ds ds
? ?
2
du dv dv
+ ?
2
+ 2?
2
+ ?
2
= 0 (20.12)
ds
2
11 22
ds
12
ds ds ds
?
jk
(i, j, k = 1, 2) are the Christo?el symbols de?ned as follows: Where ?
i
?
1
GE
u
? 2F F
u
+ F E
v
?
2
2EF
u
? EE
v
+ F E
u
= =
11 11
?
2(EG ? F
2
)
,
2(EG ? F
2
)
1
GE
v
? F G
u
?
2
EG
u
? F E
v
= =
12 12
?
2(EG ? F
2
)
,
2(EG ? F
2
)
(20.13)
1
2GF
v
? GG
u
+ F G
v
?
2
EG
v
? 2F F
v
+ F G
u
= =
22 22
2(EG ? F
2
)
,
2(EG ? F
2
)
? These two second order di?erential equations can be rewritten as a system of four ?rst
order di?erential equations [6].
du
= p (20.14)
ds
dv
= q (20.15)
ds
dp
2 2
= ??
1
? 2?
1
11
p
12
pq ? ?
1
(20.16)
22
q
ds
dq
2 2
= ??
2
? 2?
2
11
p
12
pq ? ?
2
(20.17)
22
q
ds
4
? Euler Lagrange Equations (Calculus of Variations)
We can also ?nd this result by means of the general rules of the calculus of variations.
We want to minimize
?
b
?
b
I = ds = f(u, v, v˙)du (20.18)
a a
where
dv
2
f(u, v, v˙) =
?
E + 2Fv˙ + Gv˙ , v˙ = (20.19)
du
This leads to the condition
πf d πf
(20.20)
πv
?
du πv˙
= 0
from which we can derive the same di?erential equation for geodesics.
20.1.4 Two-point boundary value problem
? We can solve a system of four ?rst order ordinary di?erential equations (20.14) to (20.17)
as
– Initial-value problem (IVP), where all four boundary conditions are given at one
point, or as
– Boundary-value problem (BVP), where four boundary conditions are speci?ed at
two distinct points.
? The ?rst order di?erential equation for a boundary value problem can be written in
vector form as:
dy
= g(s, y), y(A) = (u
A
, v
A
, p
A
, q
A
)
T
, y(B) = (u
B
, v
B
, p
B
, q
B
)
T
(20.21)
ds
where p
A
, p
B
, q
A
and q
B
are unknowns,
y = (u, v, p, q)
T
(20.22)
2
12
pq ? ?
1 2 2
22
q
2
)
T
g = (p, q, ??
1
? 2?
1
22
q , ??
2
? 2?
2
12
pq ? ?
2
(20.23)
11
p
11
p
? There are two commonly used approaches to the numerical solution of BVP.
1. Shooting method: easy to implement but unstable.
2. Relaxation method: more sophisticated but stable.
? Shooting method
– We assume values at s = A, which are not given as boundary conditions at s = A
and compute the solution of the resulting IVP to s = B.
5
?
? ?
(?)
(?)
– The computed values of y(B) will not, in general, agree with the corresponding
boundary condition at s = B.
– Consequently, we need to adjust the initial values and try again.
– The process is repeated until the computed values at the ?nal point agree with the
boundary conditions and referred as shooting method.
– Formulation: Using the ?rst fundamental form, given p
A
we can obtain q
A
from
2
?Fp
A
± F
2
p
A
? G(Ep
2
q
A
? 1)
A
= .
G
Thus we assume p
A
and solve the di?erential equation as IVP using, say Runge-
Kutta method. Here we also have to assume the entire arc length of the geodesic
path s to stop the integration. Thus the unknowns can be considered as p
A
and s.
?
), the di?erence can be given
B
?
, v
B
If we denote the computed value of (u
B
, v
B
) as (u
as (u
This can be done by employing the Newton’s method
T?
)? v
B
B
?
? u , v
B
B
. We need to adjust p
A
and s to make the di?erence zero.
?1
? ?
?
B
?s
?
B
?u?u
? ? ? ? ? ?
? u
B
p
A
p
A
u
?p
A
B
?= ?
B
?s
?
B
?v?v
s s
? v
B
v
Bn+1
?p
A
n
n
where
?
B
?
( )p
A
B
?
( + ? ) ?p p u
A A
B
?
B
?
( )s
B
?
( + ? ) ?s s u
B
πu πuu u
= , =
πP
A
?p
A
πs ?s
?
B
v
?
( )p
A
B
?
( + ? ) ?p p v
A A
B
?
B
?
( )s
B
?
( + ? ) ?s s v
B
πv πv v
= , =
πP
A
?p
A
πs ?s
Relaxation method ?
dy
– The second method is based on a ?nite di?erence approximation to
ds
on a mesh
of points in the interval [A, B].
– This method starts with an initial guess and improves the solution iteratively and
referred as, direct method, relaxation method or ?nite di?erence method.
– The shooting method is often very sensitive to the unknown initial angles at point
A and unless a good initial guess is provided, the integrated path will never reach
the other point B, while the relaxation method starts with two end points ?xed and
relaxes to the true solution and hence it is much more stable.
– Let us consider a mesh of points satisfying A = s
1
< s
2
< . . . < s
m
= B. We
approximate the n ?rst order di?erential equations by the trapezoidal rule [8].
1Y
k
? Y
k?1
= [G
k
+ G
k?1
], k = 2, 3, . . ., m (20.24)
s
k
? s
k?1
2
6
a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0
where the n-vectors Y
k
, G
k
are meant to approximate y(s
k
) and g(s
k
) with bound-
ary conditions
Y
1
= β = (u
A
, v
A
, p
1
, q
1
)
T
, Y
m
= ? = (u
B
, v
B
, p
m
, q
m
)
T
(20.25)
– Y
1
has n
1
= 2 known components, while Y
m
has n
2
= n ? n
1
= 4 ? 2 = 2 known
components.
u2u1 u3 u4 uM?2 uM?1 uM
a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a32
a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a32
Boundary Conditions
Unknowns
a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0
v1 v2 v3 v4 vM?2 vM?1 vM
a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a32
a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0p1 p2 p3 p4 pM?2 pM?1 pM
a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a32
q1 q2 q3 q4
. . . . . . . . .
qM?2 qM?1 qM
a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0a1a0a1a0a2a0a1a0a1a0a1a0
s1 s2
s3 s4
sM?2 sM?1 sM
Figure 20.2: Mesh points.
s
– This discrete approximation will be accurate to the order of h
2
(h = max
k
{s
k
?
k?1
}). Equation (20.24) forms a system of (m? 1)n nonlinear algebraic equations
with (m ? 1)n unknowns.
– Let us refer to equation (20.24) as
1
F
k
= (F
1,k
, F
2,k
, . . ., F
n,k
)
T
=
Y
k
? Y
k?1
[G
k
+ G
k?1
] = 0 (20.26)
s
k
? s
k?1
?
2
– and the boundary conditions (20.25) as
F
F
1
= (F
1,1
, F
2,1
, . . ., F
n
1
,1
)
T
= Y
1
? β = 0, (20.27)
m+1
= (F
1,m+1
, F
2,m+1
, . . ., F
n
2
,m+1
)
T
= Y
m
? ? = 0
– then we have mn nonlinear algebraic equations
F = (F
T
2
, . . .,F
T
1
,F
T
m+1
)
T
= 0 (20.28)
7
Points Tolerance Iterations Geodesic Distance
m κ
N
L M R L M R
101 1.0E-3 10 1 10 1.661 1.865 1.661
Table 20.1: Numerical conditions and results for the computation of the geodesic path between
corner points of the wave-like surface.
– and can be computed by quadratically convergent Newton iteration, if a su?ciently
accurate starting vector Y
(0)
= (Y
T
2
, . . . , Y
T
1
, Y
T
m
)
T
is provided. The Newton iter-
ation scheme is given by
Y
(i+1)
[J
= Y
(i)
+ ?Y
(i)
(20.29)
(i)
]?Y
(i)
= ?F
(i)
(20.30)
where [J
(i)
] is the mn by mn Jacobian matrix of F
(i)
with respect to Y
(i)
.
x
y
z
Figure 20.3: Geodesic paths on the wave-like bicubic B-spline surface between points of two
corners.
20.1.5 Example
Bilinear surface r(u, v) = (u, v, uv).
E
E
E = 1 + v
2
, F = uv, G = 1 + u
2
u
= 0, F
u
= v, G
u
= 2u
v
= 2v, F
v
= u, G
v
= 0
8
x
y
z
Figure 20.4: Convergence of the right geodesic path in Figure 20.3.
?
1
?2(uv) · v + uv · (2v)
= = 0
11
2
2[(1 + v
2
)(1 + u
2
) ? u v
2
]
?
2
= ?
1
22
= 0
22
= ?
2
?
11
1
v
=
12
?
u
2
+ v
2
+ 1
2
u
=
12
u
2
+ v
2
+ 1
Finally the di?erential equations are given by
du
= p
ds
dv
= q
ds
dp ?2u
=
ds u
2
+ v
2
+ 1
pq
dq ?2v
=
ds u
2
+ v
2
+ 1
pq.
9
20.2 Developable surface
Developable surfaces are a special class of surfaces that can be developed or unfolded onto a
plane without stretching or tearing [1, 9, 5, 16, 15].
20.2.1 Motivation
? In shipbuilding, doubly curved plates are manufactured by roller and line heating pro-
cesses, while the singly curved plates (developable surface) are manufactured by roller
only.
? The use of developable surfaces has several advantages such as lower labor costs in
construction, smaller capital investment in equipment, ease of repair and simple tools for
construction.
? In automobile production, body panels, upholstery and window glass are developable
surfaces.
20.2.2 De?nition
? A ruled surface is de?ned as a surface generated by the motion of a straight line referred
as a generator or ruling [24].
? The mathematical representation of a ruled surface is given by
r(u, v) = r
A
(u) + vβ(u) (20.31)
where r(u) is a 3D curve called the directrix or base curve of the ruled surface and β(u)
is a unit vector which gives the direction of the ruling at each point on the directrix see
Figure 20.5.
ruling or
generator
directrix
α(u)
r(u,v)=
+
rA(u)
rA(u) vα(u)
Figure 20.5: De?nition of ruled surface
10
? An alternate expression based on rulings joining corresponding points on two space curves
r
A
(u) and r
B
(u) is given by
r(u, v) = (1 ? v)r
A
(u) + vr
B
(u) (20.32)
Bilinear surface is a ruled surface. ?
r
r
A
(u) = (1 ? u)b
00
+ ub
10
B
(u) = (1 ? u)b
01
+ ub
11
r(u, v) = (1 ? v)r
A
(u) + vr
B
(u)
= (1 ? u)(1 ? v)b
00
(u) + (1 ? v)ub
10
+ v(1 ? u)b
01
+ uvb
11
? A ruled surface is a developable surface if and only if [7]
˙ β) = 0 (20.33)r
A
· (β × ˙
where × and · are the cross and dot products or
(r
A
? r
B
) · (r˙
A
× r˙
B
) = 0 (20.34)
? The following statements are the equivalent necessary and su?cient conditions for a
surface to be developable [19].
1. Gaussian curvature is zero.
2. Geodesics on a developable surface can be mapped onto straight lines in the plane,
and the straight lines in plane can map into geodesics on a developable surface.
3. The normal vectors on a developable surface along the ruling are parallel.
4. Developable surfaces possess the same tangent plane at all points of the same gen-
erator.
? If the direction of the ruling β(u) is constant, the condition for developability is auto-
matically satis?ed since β(u) = 0. This implies that the ruled surface is a cylinder.˙
? If the direction of the ruling β(u) is given by r
A
(u) ? a, then the condition (20.33)
becomes r˙
A
(u) · ((r
A
? a) × r˙
A
(u)) → 0 and hence the surface is developable. In this case
the surface is a cone with apex a.
? Finally, if β(u) is a tangent vector to r
A
(u), then again the condition (20.33) is satis?ed,
and the resulting developable surface is called a tangential developable surface. In this
case the base curve coincides with the so-called the edge of regression or cuspidal edge.
11
20.2.3 Developable surface in terms of B′ezier surface
? Let us de?ne a developable B′ezier surface in terms of two B′ezier curves (directrices) and
rulings between pairs of points from each curve [1], see Figures 20.6.
? We restrict the two directrices (r
A
(t), r
B
(t)) to lie in parallel planes so that r˙
B
(t) =
ζ(t)r˙
A
(t). Here ζ(t) denotes a linear function of t.
? Then the condition for developability becomes (ζ(t) ? 1)r˙
A
× (r
B
? r
A
) · r˙
A
0 and hence→
it is automatically satis?ed.
? Given the control points of the design curve r
A
(t), we want to determine those of r
B
(t).
Example: r
A
(t) quadratic and r
B
(t) cubic
a
1
a
2
b
0
b
1
b
2
b
3
generator
directrix
r
B
directrix
r
A
a
0
Figure 20.6: A B′ezier Developable Surface with Quadratic r
A
(t) and Cubic r
B
(t)
Two B′ezier Curves?
2
r
A
(t) = s a
0
+ 2sta
1
+ t
2
a
2
2
r
B
(t) = s
3
b
0
+ 3s tb
1
+ 3st
2
b
2
+ t
3
b
3
where s = 1 ? t
? Tangent vectors
r˙
A
(t) = 2[s(a
1
? a
0
) + t(a
2
? a
1
)]
2
r˙
B
(t) = 3[s
2
(b
1
? b
0
) + 2st(b
2
? b
1
) + t (b
3
? b
2
)]
Scalar function
?
ζ(t) = ζ
0
(1 ? t) + ζ
1
t = ζ
0
s + ζ
1
t
12
? Substitute r˙
A
(t), r˙
B
(t) and ζ(t) into r˙
B
(t) = ζ(t)r˙
A
(t) and collect terms and equating
the coe?cients of the independent functions s
2
, 2st and t
2
to zero.
? We get three equations
2ζ
0
(a
1
? a
0
) = 3(b
1
? b
0
)
ζ
0
(a
2
? a
1
) + ζ
1
(a
1
? a
0
) = 3(b
2
? b
1
)
2ζ
1
(a
2
? a
1
) = 3(b
3
? b
2
)
? 6 scalar equations with 10 scalar unknowns (b
0
, b
1
, b
2
, b
3
, ζ
0
, ζ
1
)
? We can ?x b
0
and b
3
.
20.2.4 Development of developable surface (?attening)
? In the manufacture of developable surfaces, it is necessary to calculate the plane devel-
opment of these surfaces [7].
? The development is based on the isometric mapping [13].
– The length of any arc on the developable surface is the same as that of on the
developed plane.
– The coe?cients of the ?rst fundamental forms at the corresponding points are the
same.
– Isometric surfaces have the same Gaussian curvature at corresponding points. Cor-
responding curves on those surfaces have the same geodesic curvature at correspond-
ing points.
– Every isometric mapping is conformal (the angle of intersection of every arbitrary
pair of intersecting arcs are preserved).
Procedures ?
– Map one of the directrix isometrically onto the plane by integrating a system of
ordinary di?erential equation from A to C, see Figure 20.7.
d
2
x dy
g
(s) = 0 (20.35)
ds
2
? ρ
ds
d
2
y dx
+ ρ
g
(s) = 0 (20.36)
ds
2
ds
– Compute the angle of two isoparametric lines at A.
r
u
(0, 0) r
v
(0, 0)
cos β = (20.37)
r
u
(0, 0)
·
r
v
(0, 0)| | | |
13
Generator
B
C
D
θ
Directrix
Directrix
Generator
A
Figure 20.7: Developed surface.
– Since we know the length and direction of the generator at A, we may obtain the
point B.
– Integrate a system of ordinary di?erential equation from B to D.
– Connect C and D.
14
?
20.3 Umbilics
20.3.1 Motivation
? Identi?cation of singular points in the principal direction ?eld
? Fingerprints of shapes
? Object matching
? Object recognition
20.3.2 De?nition
An umbilic is a point on a surface where the normal curvatures in all directions are equal and
the principal curvature directions are indeterminate. Locally, a surface around an umbilical
point can be best approximated by a circle whose radius is equal to a radius of curvature at
the umbilical point.
20.3.3 Computation of umbilical points
For NURBS surfaces, umbilical points can be calculated by solving a non-linear system of
equations derived from the de?nition using the Gaussian (K) and the mean (H) curvatures as
follows [24]:
ρ
1,2
(u, v) = H(u, v) ± H
2
(u, v) ? K(u, v). (20.38)
Let W (u, v) = H
2
? K. The principal curvatures, ρ
1,2
, are real valued functions so that W √ 0
must hold. From the de?nition of the umbilical point we have W (u, v) = 0. With these two
conditions combined, we can infer that at an umbilical point, W (u, v) has a global minimum
[17, 18]. Here, we assume that W is at least C
2
smooth. Then, the condition that W has a
global minimum at an umbilic implies that ≡W = 0. Therefore, at an umbilic the following
equations hold [18]:
πW (u, v) πW (u, v)
W (u, v) = 0, = 0, = 0. (20.39)
πu πv
Given a polynomial parametric surface patch such as a rational B′ezier surface patch, we
can set W =
P
N
, where P
N
and P
D
are polynomials in u and v. With the condition W √ 0,
P
P
D
N
√ 0 is assured since P
D
> 0 is always true under the regularity condition of the surface
[24]. The equation W = 0 is equivalent to P
N
= 0. The ?rst derivative of W is
?W
=
?x
i
(
?P
N
?P
D
D
(i = 1, 2), where x
1
= u and x
2
= v, which is reduced to
?W
= (
?P
N
)/P
D
?x
i
P
D
? P
N
?x
i
)/P
2
?x
i
?x
i
using P
N
= 0. Therefore, equations (20.39) are reduced to [18]
πP
N
πP
N
P
N
(u, v) = 0, = 0, = 0. (20.40)
πu πv
To locate isolated umbilical points, a set of equations (20.40) can be solved by the rounded
interval projected polyhedron algorithm [23, 21]. But when the IPP algorithm encounters re-
gions of umbilical points, it slows down dramatically. A di?erent approach, called the adaptive
15
quadtree decomposition which uses the quadtree decomposition and the convex hull properties
of Bernstein polynomials, can be adopted in such a case [12].
20.3.4 Classi?cation
Umbilical points can be isolated or form lines or regions. They are classi?ed into two types
based on their stability with respect to small perturbations:
? generic
? non-generic
Generic umbilical points are stable with respect to small perturbations.
Isolated generic umbilical points are further categorized into three types as shown in Figure
20.8 [2]:
? lemon
? star
? monstar
Star type umbilical points are further classi?ed into the hyperbolic star and the elliptical star
type umbilical points. Several methods are available for the classi?cation of isolated generic
umbilical points. What follows is an introduction to techniques of umbilical point classi?cation.
Lemon Star Monstar
Figure 20.8: Three generic umbilics
Index
The type of isolated generic umbilical points can be determined by the index. The index is the
amount of rotation that a straight line segment tangent to the lines of curvature experiences
16
when rotating in the counterclockwise direction along a small closed path around an umbilic
[18]. The equation for direct calculation of the index is given as follows [18]:
1
n
?
Index =
2?
i=0
??
i
(20.41)
? ? ? ?
1 1 L+?E M +?F
where ??
i
= ?
(i+1) mod n
??
i
, ? ? ? ??
i
? ? and ?
i
= tan
?1
? , or tan
?1
? .
2 2 M +?F N +?G
Figure 20.9 shows the principal direction ?elds around two distinct types of umbilical points.
The index can distinguish the star type umbilical point from the monstar or lemon type umbil-
Figure 20.9: Umbilics with principal direction ?elds:
1
ical point. If the index is ?
2
, then the umbilical point is of the star type, whereas the umbilical
1
point is of the monstar or lemon type if the index is
2
. Being topological, this distinction is
very robust [10].
Complex θ plane method
An umbilical diagram shown in Figure 20.10 [22] is a comprehensive and easy way to distinguish
the type of an isolated generic umbilical point. In order to use this diagram, the local surface
near an umbilical point has to be represented as a height function or the Monge form with
respect to a local coordinate system as follows [18]:
r = (x, y, h(x, y)). (20.42)
The height function h(x, y) is Taylor expanded at the origin of the local coordinate system.
Then we have
h(x, y) =
ρ
(x
2
+ y
2
), (20.43)?
2
+
1
(ax
3
+ 3bx
2
y + 3cxy
2
+ dy
3
) + O(4),
6
where ρ is the normal curvature at the umbilical point. Let us set
2
C(x, y) = ax
3
+ 3bx
2
y + 3cxy + dy
3
, (20.44)
and denote C(x, y) as cubic form. This expression implies that the local structure of a surface
near an umbilical point is dominated by the coe?cients of C(x, y), i.e. by (a, b, c, d), which
17
determine the type of the umbilical point [20, 22]. It is convenient to represent the cubic part
C(x, y) in the complex plane for analysis purposes. If we set ? = x + iy, then the cubic part
C(x, y) becomes
?
3
C(?) = ω?
3
+ 3??
2
? + 3???
2
+ ω? , (20.45)
with
1
ω =
8
[(a ? 3c) + i(d ? 3b)], (20.46)
1
? = [(a + c) + i(b + d)],
8
where ω = 0. We can express (20.45) in a coordinate system rotated about the normal vector ≥
without losing any essential features to make the coe?cient of ?
3
be equal to 1 [22]. Using
? = ω
1
3
?, the equation (20.45) becomes
3
C(? ?
3 2
? ?
2
) = + 3θ? + 3θ?
?
+ ? , (20.47)
where θ = ?ω
?
1
3
ω
?
2
3
. This means that the cubic part C(x, y) is parametrized with respect to
a single complex variable θ [4, 22]. Therefore, the variations of C(x, y) can be mapped onto
the complex plane [4, 20, 22]. When ω = 0, the equation (20.45) is reduced to
?
C(?) = 3??(?? + ??). (20.48)
This equation corresponds to the in?nity in the θ plane [22, 4, 20], which is not considered in
this discussion.
?
?
Depending on the structure of C(x, y) (in turn C(?)), three characteristic lines are deter-
mined as follows [22, 20]:
1
: β ?
1
3
(2e
i?
+ e
?2i?
),?
θ = 1, ? | |
? ?
2
: β ? (2e
i?
+ e
?2i?
),
where ?
1
and ?
2
are maps from β to the complex θ-plane. They divide the complex plane
into sub-regions as shown in Figure 20.10. Each sub-region corresponds to a speci?c type of
an umbilical point. In Figure 20.10, ES means the elliptic star, HS the hyperbolic star, MS
the monstar and L the lemon. If θ falls on a dividing curve, then the corresponding umbilical
point is of non-generic type. The behavior of such an umbilical point can be analyzed with
more higher order terms [18]. Using this diagram, the type of an umbilical point is easily
determined, see [22, 4, 20].
20.3.5 Characteristic lines
: ? ?
1
3
(2e
i?
+ e
?2i?
)?
1
The cubic form (20.47) is parabolic on the deltoid ?
1
[22]. This implies that there are three
root lines of (20.47) but two of them coincide. Inside the deltoid, the cubic form (20.47) is
18
? ?
?
1
Γ
ESHS
HS
HS
MS
MS
MS
L
L
L
|ω|=1
2
Γ
Figure 20.10: The umbilic diagram adapted from [22]
elliptic [22]. It is hyperbolic outside the deltoid. This classi?cation is directly related to the
number of ridge lines passing through an umbilical point, and the existence of extrema of the
principal curvatures near the umbilical point [4, 10, 18]. Here, a ridge point is de?ned as an
extremum point of a line of curvature, and a ridge line is a set of such points [10, 20]. Inside
the deltoid, the number of ridge curves is three, the extrema of principal curvatures exist, and
an umbilical point is of the elliptical star type [4, 10]. On the other hand, outside the deltoid
only one ridge curve passes through an umbilical point, no extremum of a principal curvature
exists, and the umbilical point is of the hyperbolic star type [4, 10].
|?| = 1
On the circle θ = 1, the cubic form (20.47) is right-angled [22]. When the root directions of | |
the Hessian (20.49) are mutually orthogonal, we call that the cubic form (20.44) is right-angled
[22]. Here, the Hessian is de?ned as
?
?
2
C ?
2
C
?
?x?y
?
?
H
e
(x, y) =
1
det
?
?
?x
2
?
2
C ?
. (20.49)
6 ?
?
2
C
?
?x?y ?y
2
This implies that the maximum and minimum lines of curvature are orthogonal at an umbilical
point to form approximately a plain rectangular grid pattern [22]. This circle is related to the
index [10]. Inside the circle, the index is ?
1
2
, and an umbilical point is of the star type.
19
Outside the circle, the index is
1
2
, and the umbilical point is classi?ed as the lemon or monstar
type [10, 18]. On the circle, the index value is zero, and the umbilical point is of non-generic
type. This circle also has a relation with the birth/death of generic umbilical points under the
evolution of a surface [22].
i?
+ e
?2i?
)?
2
: ? ? (2e
Another cubic form, called the Jacobian cubic form is de?ned to explain the deltoid ?
2
as
follows:
2 2 3
U (x, y) = bx
3
? (2c ? a)x y ? (2b ? d)xy ? cy , (20.50)
whose root lines are tangent to the lines of curvature near an umbilical point [10]. On the
deltoid ?
2
, the cubic form (20.50) becomes parabolic [22]. The Jacobian cubic form (20.50) is
related to the number of extrema of the cubic form (20.44) [22]. The cubic form (20.44) can be
represented as C(r, β) in polar coordinates with x = r cos β and y = r sin β, and the expression
of the direction in which the local extrema of the cubic form C(r, β) occur, i.e.
dC(?)
= 0 [18]
d?
is reduced to the Jacobian cubic form (20.50). Inside the deltoid ?
2
, there are three real root
lines of the Jacobian cubic form (20.50) or three directions of the extrema of the cubic form
(20.44), which implies that three lines of curvature converge to an umbilical point [10, 22].
This umbilical point is classi?ed as the star or monstar type [10, 22]. Outside the deltoid ?
2
,
there is one root line of the Jacobian cubic form (20.50), and no extremum of the cubic form
(20.44) exists. An umbilical point of this case is of the lemon type [10, 22].
20
20.4 Parabolic, ridge and sub-parabolic points
20.4.1 Motivation
? Sophisticated classi?cation of umbilics
? Stronger matching conditions than the umbilic matching
? Registration
? Analysis of surface evolution
20.4.2 Focal surfaces
For any point on a surface r, two points f
1
and f
2
are called focal points which are de?ned as
follows:
1 1
f
1
= r + N, f
2
= r + N (20.51)
ρ
1
ρ
2
where N is a normal vector at P and ρ
1
and ρ
2
are the maximum and minimum principal
curvatures, respectively. The focal surface is a set of focal points which can be considered as
the envelope of the normal to the surface. The focal surface will play an important role in
studying ridges, crests and sub-parabolic lines.
Figure 20.11: Focal Surfaces
21
20.4.3 Parabolic points
Contact function
Consider a surface r in Monge form as follows:
2
z = h(x, y) = a
0
x
2
+ a
1
xy + a
2
y + H.O.T, (20.52)
so that the tangent plane at the origin is the x,y-plane. The parabolic points can be identi?ed
2
by studying the intersection between r and the tangent plane. If a
1
= 4a
0
a
2
then the surface
has A
2
contact with the tangent plane.
De?nition
The parabolic points on a regular surface are points where the tangent plane has specially
high ’contact’ with the surface, i.e. A
2
contact [3]. Alternatively, they can be interpreted as
the points which separate the elliptic and hyperbolic regions of the surface, namely, points or
curves where the Gaussian curvature is zero.
20.4.4 Ridge points
Contact functions
The sphere of curvature is a sphere centered at one of the centers of principal curvature and
having a radius equal to the corresponding radius of curvature. Consider a sphere centered at
(0, 0, r) and passing through the origin such that it is tangent to the surface r there. Then the
equation of the sphere is
2 2 2
x + y + (z ? r)
2
= r (20.53)
The contact function is de?ned with substitution of z = h(x, y).
2
g(x, y) = x
2
+ y + (h(x, y) ?r)
2
?r
2
= 0 (20.54)
Expanding g(x, y) as a power series in x and y gives
r
2 2 2
g(x, y) = x
2
(1 ? rρ
1
) + y (1 ? rρ
2
) ? (b
0
x
3
+ 3b
1
x y + 3b
2
xy + b
3
y
3
)
3
r
(c
0
x
4
ρ
2
y
2
)
2
?
12
+ ···) + (
ρ
1
x
2
+ + (20.55)
2 2
···
De?nition
? The ridge point of a smooth surface is a point where one of the sphere of curvature has
more degenerate contact than the usual A
2
contact and the curve of intersection becomes
ρ
1
g(x, y) = y
2
(ρ
1
? ρ
2
) ?
1
(3b
1
x
2
y + 3b
2
xy
2
+ b
3
y
3
)
3
ρ
3
1
(c
0
x
4
?
12
+ ···) +
1
x
4
+
4
···
22
b
? ?
2
1
b
2
b
3
2
x y= (ρ
1
? ρ
2
) y ?
2(ρ
1
? ρ
2
)
2
?
2(ρ
1
? ρ
2
)
xy ?
6(ρ
1
? ρ
2
)
1
1
+ (ρ
1
? ρ
2
)(c
0
? 3ρ
3
1
))x
4
+ (20.56)?
12(ρ
1
? ρ
2
)
(3b
2
···
1
– If r =
1
(or r =
?
2
) and b
0
= 0 (or b
3
= 0), then the contact function (20.55)
?
1
becomes A
3
singularity.
– The coe?cient P
1
= 3b
2
1
) should not be zero.
1
+ (ρ
1
? ρ
2
)(c
0
? 3ρ
3
? A point is a ridge point relative to a principal direction p if and only if the principal
curvature ρ
1
corresponding to p has an extremum along the line of curvature in the
direction p.
? Ridge points are the pre-image of cuspidal edges on a focal surface of a smooth surface.
20.4.5 Sub-parabolic points
? The point on which one principal curvature has an extremum relative to motion along
the other line of curvature is called the sub-parabolic point on the surface.
? It is the pre-image of a parabolic point on the focal surface.
23
Bibliography
[1] G. Aumann. Interpolation with developable B′ezier patches. Computer Aided Geometric
Design, 8(5):409–420, November 1991.
[2] M. V. Berry and J. H. Hannay. Umbilic points on Gaussian random surfaces. Journal of
Physics A., 10(11):1809–1821, 1977.
[3] J. W. Bruce, P. J. Giblin, and F. Tari. Parabolic curves of evolving surfaces. International
Journal of Computer Vision, 17(3):291–306, 1996.
[4] J. W. Bruce, P. J. Giblin, and F. Tari. Ridges, crests and sub-parabolic lines of evolving
surfaces. International Journal of Computer Vision, 18(3):195–210, 1996.
[5] J. S. Chalfant and T. Maekawa. Design for manufacturing using B-spline developable
surfaces. Journal of Ship Research, 42(3):207–215, 1998.
A. Bj¨[6] G. Dahlquist and
?
orck. Numerical Methods. Prentice-Hall, Inc., Englewood Cli?s,
NJ, 1974.
[7] I. D. Faux and M. J. Pratt. Computational Geometry for Design and Manufacture. Ellis
Horwood, Chichester, England, 1981.
[8] J. H. Ferziger. Numerical Methods for Engineering Applications. Wiley, 1981.
[9] W. H. Frey and D. Bindschadler. Computer-aided design of a class of developable B′ezier
surfaces. R&D Publication 8057, General Motors, September 1993.
[10] P. W. (editor) Hallinan. Two and Three Dimensional Patterns of the Face. A.K. Peters,
Ltd., Wellesley, Massachusetts, 1999.
[11] R. Kimmel, A. Amir, and A. M. Bruckstein. Finding shortest paths on surfaces using
level sets propagation. IEEE Transactions on Pattern Analysis and Machine Intelligence,
17(6):635–640, June 1995.
[12] K. H. Ko, T. Maekawa, N. M. Patrikalakis, H. Masuda, and F.-E. Wolter. Shape intrinsic
?ngerprints for free-form object matching. In G. Elber and V. Shapiro, editors, 8th ACM
Symposium on Solid Modeling and Applications, Seattle, WA, June 2003, In press.
[13] E. Kreyszig. Di?erential Geometry. University of Toronto Press, Toronto, 1959.
[14] T. Maekawa. Computation of shortest paths on free-form parametric surfaces. Journal of
Mechanical Design, Transactions of the ASME, 118(4):499–508, December 1996.
[15] T. Maekawa and J. S. Chalfant. Computation of in?ection lines and geodesics on devel-
opable surfaces. Mathematical Engineering in Industry, 7(2):251–267, 1998.
[16] T. Maekawa and J. S. Chalfant. Design and tessellation of B-spline developable surfaces.
Journal of Mechanical Design, Transactions of the ASME, 120(3):453–461, September
1998.
[17] T. Maekawa and N. M. Patrikalakis. Interrogation of di?erential geometry properties for
design and manufacture. The Visual Computer, 10(4):216–237, March 1994.
[18] T. Maekawa, F.-E. Wolter, and N. M. Patrikalakis. Umbilics and lines of curvature for
shape interrogation. Computer Aided Geometric Design, 13(2):133–161, March 1996.
[19] M. J. Mancewicz and W. H. Frey. Developable surfaces: Properties, representations and
methods of design. R&D Publication 7637, General Motors, 1992.
[20] R. Morris. The sub-parabolic lines of a surface. In G. Mullineux, editor, The Mathematics
of Surfaces VI, Proceedings of the 6th IMA Conference on Mathematics of Surfaces VI,
pages 79–102, Oxford, UK, 1996. Clarendon Press.
[21] N. M. Patrikalakis and T. Maekawa. Shape Interrogation for Computer Aided Design and
Manufacturing. Springer-Verlag, Heidelberg, February 2002.
[22] I. R. Porteous. Geometric Di?erentiation for the intelligence of curves and surfaces.
Cambridge University Press, Cambridge, 1994.
[23] E. C. Sherbrooke and N. M. Patrikalakis. Computation of the solutions of nonlinear
polynomial systems. Computer Aided Geometric Design, 10(5):379–405, October 1993.
[24] D. J. Struik. Lectures on Classical Di?erential Geometry. Addison-Wesley, Cambridge,
MA, 1950.