Shape offsets via level sets

R Kimmel and A M Bruckstein*

In Appendix A, it is proved that the evolution equation

An algorithm for shape offsetting & presented that & based
on level-set propagation. This algorithm avoids the
topological problems encountered in traditional offsetting { ~ = IN(s, t)
(5)
algorithms, and it deals with curvature singularities
by including an 'entropy condition' in its numerical X(s, 0) = Xo(s)
implementation. postulating that each point on the curve moves in the
direction of the instantaneous normal
shape offsets, prairiefires, numericalalgorithms, Huyffensprinciple
1
N(s, t) = [ - y~(s, t), x,(s, t)] T
(x2 (s, t) + y2~(s,t)) 1/2
Algorithms for shape offsetting are of great importance with velocity 1 provides the same flow as Equations 4
in computer-aided design (CAD), computer-aided manu- almost everywhere (where problems do not occur).
facturing (CAM),the numerical control (NC) of machines, This propagation rule is the so-called 'prairie-fire'
computer graphics and related fields. The need for offset model for shape evolution (see, for example, References
shapes arises in applications such as the numerical 1 and 2). If we could determine the solution of Equations
control of sawing machines or milling machines in the 4 or Equations 5 for all t > 0, we could generate not only
car industry. the offset shape (at t=L), but a whole class of offset
The problem of shape offsetting can be formulated as shapes. Hence, it is important to analyze these equations
follows: given a simple, closed planar curve closely and obtain good numerical algorithms for their
Xo(s) = Ix(s), y(s)] T (1) solution. This was indeed done in the context of
flame-propagation models and shape analysis in com-
where s is an arbitrary curve parameterization, find an puter vision, with an emphasis on the analysis of the
offset curve which is simple and closed (or has simple possible singularities arising on the propagated curves.
and closed components) and is almost everywhere given The following singularities (or 'shocks') may occur3:
by
• If there is a local maximum of the curvature such that
XL(s) = Xo(s) + N(s, 0)L (2) 1/k<L, a curvature-discontinuity 'shock' is formed
Equation 2 represents a curve running 'parallel' to Xo(s), after time t = 1/k (see Figure la), and propagating the
where L is the displacement of the offset curve, and N(s, 0) curve beyond t = 1/k < L via Equation 4 forms a 'cusp'
is the unit normal at the point Xo(S) given by (see Figure lb).
• At places where the original curve has breakpoints and
1 the derivatives are not well defined, problems of
N(s, 0) = (x,2 (s) + y,~ (s)) 1/2 [ - Y s ( S ) , X s ( S ) ] T (3)
determining the intersection of the propagated lines
arise (see Figure lc).
Consider X(s, t) to be a curve continuously changing in • Self intersections of the propagating curve may occur,
time, so that, for all t, X(s, t) = Xo(s) + tN(s, 0). This curve
giving rise to some difficult topological problems (see
evolution can be described differentially by Figure ld).

{~t 't) = 1N(s'0)

X(S,0)--Xo(S)
(4)
HISTORICAL REVIEW
The issue of generating offset curves has often been dealt
with in the CAD literature. Several approaches to the
Electrical Engineering Department, Technion, Haifa 32000, Israel problem have been proposed. In the CAD/CAM appli-
*Computer ScienceDepartment, Technion, Haifa 32000, Israel
Paper received: 9 March 1992. Revised:6 July 1992 cations that we consider, there are often shape boundaries

154 0010-4485/93/030154-09 © 1993 Butterworth-Heinemann Ltd com puter-aided design

 Shape offsets via level sets

curve. The translation of each control point in the normal

direction is given by "f.=L(l+kD). Pham 7 uses the
uniform cubic B-spline to represent curves. She also uses
a version of B-splines with control points on the curve,
and the offset curves are improved by adding control
points ('knots'). Elber and Cohen s measure the distance
error to offset curves (and surfaces), and use it to place
new control points. Farin a provides a recursive procedure
to offset B6zier segments that is similar to the one
proposed by Klass. Meek and Walton 1° observe that,
II b for clothoidal splines, the offset curve remains 'in the
family' (i.e. the offset is also a clothoidal spline).
Wang and Jiang 11 use a vector representation of the
shape comprising circular arcs and straight segments only.
Clockwise directed vectors form the shape boundaries,
and multiplication of consecutive vectors yields vectors
indicating the offset direction. In his recent book x2, Held
attacks the shape-offsetting problem using Voronoi
diagrams. His aim is to design the course of a tool creating
a hole 'pocket' in solid material. He restricts his methods
C to linear and circular segments whose representation can
Figure 1. Problems arise when singularities form in easily be fed into CNC machines. An alternative
propagated curve approach to the offset problem is provided by Saeed et
al. in Reference 13. They propose the use of morpho-
logical methods to formulate the offset operation. Indeed,
offset shapes are closely related to dilated or eroded
that are defined by standard parameterized curves such shapes, as defined in mathematical morphology. Related
as line segments, circular arcs or splines. In such cases, problems are skeleton .finders in computer vision and
one could concentrate on determining evolution equa- morphology 14, fat curves in computer graphics ~5, and
tions for the control points defining such curves. the calculation of Euclidean distance maps from plane
However, the offset curves of splines are not necessary c u r v e s 16.
splines themselves, and straightforward translations of The algorithms proposed so far for the offsetting
spline control points only rarely produce correct offset problem deal with edge-intersection problems, shocks,
curves. Klass 4 approximates the offset of a B,spline cusps and self loops in complex, and rather unnatural,
segment by moving the endpoints of the curve, and ways. In these algorithms, the curve-offsetting stage is
calculating the new tangents based on replacing the followed by a procedure aimed at detecting the afore-
curvature k with a ~that obeys 1/k = l / E - L. He proposes mentioned problems and repairing the offset curve
an iterative procedure for calculating the tangents. A accordingly. In the sequel, we present a way to approxi-
spline approximation of the offset curve is generated, and mate the offset shape using level-set propagation on a
then the distance between the two curves is measured, rectangular grid designed to produce results of the desired ~
and, if this distance deviates considerably from L, the accuracy. This is done by applying an algorithm devised
spline is split into several spline segments and the process by Osher and Sethian ~v for the stable and efficient
goes on. Tiller and Hanson s replace the B-spline propagation of wavefront curves in the plane. This
parametric representation of the shape by a 'rational numerical method generates offset curves according to
B-spline' representation with the control points located Equation 5 and a physically motivated 'entropy con-
on the curve itself. The offset operation is then carried dition', and it inherently avoids the topological problems
out by moving the control points a distance of L on the that required special attention with previous algorithms.
normal direction. These approaches require the develop-
ment of some sophisticated procedures to deal with
the problems of loops, shocks, cusps, self-intersection
singularities etc. Hoschek 3 attacks the problem of finding NEW ALGORITHM
self-intersection points on offset curves using an iterative We propose to generate shape offsets via an ingenious
geometrical algorithm. The algorithm is used to eliminate algorithm invented in fluid dynamics for solving equa-
the tails and loops that arise in the offset curves. It tions of the type of Equations 5. This algorithm translates
requires as input two points that are close enough to the problem of curve evolution into a problem of
each intersection point to guarantee convergence, and it 3D-surface evolution, so that the curves changing
does not solve the so-called island problem (see Figure according to Equations 5 are zero (or level) sets of the
ld). Coquillart 6 suggests a new way of translating the time-varying surface. As the 3D surface evolves, it
rational-B-spline control points to preserve circles. The inherently avoids the generation of curve shocks by
translation is controlled by the local curvature k and the implementing a physically motivated entropy condition.
distance D from the B-spline control points to the given The algorithm that produces the desired results works

volume 25 number 3 march 1993 155

R Kimmel and A M Bruckstein

on a square grid with a resolution determined by the

desired accuracy of the results, and it is based on a
recently discovered efficient numerical implementation of
the surface-evolution equations 17

Huygens principle and entropy condition

According to the Huygens principle18, the solution of the
curve propagation of Equations 5 at time t, X(s,t),
corresponds to the envelope generated by the set of all
the disks of radius t centred on the initial curve Xo(s).
)
Problems occur in curve evolution when the normals to
the initial curve collide or cross and hence the curvature
becomes singular. To obtain the solution according to
Huygens' principle after a singularity develops, an
'entropy condition' should be enforced on the propa-
gating curve. One can regard the curve as the wavefront
of a propagating prairie fire separating two areas: the a
shape interior which is not burnt, and the already-burnt
exterior. The flame propagates in the direction of the
original curve normals (the so-called 'ignition curves').
If two ignition curves collide at some time t*, neither one
should have any effect on the propagating curve at t > t*.
The principle 'what was burnt until t cannot burn beyond
t' 1s is the natural 'entropy condition' of this type of curve
evolution. See Figure 2.
The direct approach to propagating the curve can
be referred to as the 'Lagrangian' formulation, because
the coordinates (s and t) are front-dependent 19. The
Lagrangian formulation is the direct numerical approxi-
mation of Equations 5:
Ox(s, t) = y~(s, t)

f Ot
~y(s,t)_
~, Ot
(x2(s, t) + y2(s, t)) u2
x,(s,t)
(x~(s,t)+ y~(s,t)) u2
Taking the discrete approximation of x~ and y~ as central
(6)

derivatives in place (s) and a forward-derivative approxi- b

mation in time (t) yields a numerical-propagation
scheme. The direct numerical propagation of a curve Figure 2. The principle 'what was burn until t cannot burn
according to Equations 6 is both numerically unstable beyond t' is the natural 'entropy condition' o f our curve
and suffers from topological problems (see References 17 evolution
and 19). To avoid the various problems that occur in [(a) According to Huygens' principle, the front of the evolvingcurve
this approach, such as the need for reparameterization is constructed by the front of all the disks of radius t centered on the
initial curve. (b) If two ignition curves collide at some time t*, neither
in order to keep numerical stability and to solve one should have any effecton the propagating curve at t > t*. Observe
topological problems of self intersections by an external that the dotted circle does not have any effecton the front.]
control procedure, the 'Eulerian formulation', described
below, was developed.

the curves X(s, t) as if propagated by Equations 5, and

Solution via Eulerlan formulation also obey the entropy condition. If dp(x,y,t)=O along
The Eulerian scheme is a recursive procedure which X(s, t), then, by the chain rule, we. have
propagates the curve while inherently implementing the
entropy condition. Introduce a smooth function ~b(x,y, t) ~t~(x, y,t)+ ~xdP(x(s,t),y(s,t),t)x,
that is arbitrarily initialized so that ~b(x,y,0)=0 yields
the curve X(s, 0). Assume that X(s, 0) is a closed curve, + g~ ~b(x(s, t), y(s, t), t)yt = 0
and restrict 4) to be negative in the interior and positive
in the exterior of the level set ~b(x,y, 0 ) = 0.
or
The basic idea is to determine an evolution of the
surface ~b(x,y, t) so that the level sets qb(x, y, t) = 0 provide q~t+ V~bXt(s, t) = 0 (7)

156 computer-aided design

 Shape offsets via level sets

where Some schemes based on this idea, such as the

Lax-Friedrichs and Godonov schemes, are presented in
Reference 17. The simplest flow function from our
implementation point of view is the so-called HJ flow,
is the gradient of the function 4~(x,y, t) at time t, at the for which, for H(u)=f(u2), the numerical flow can be
point (x, y). The scalar velocity of each curve point in its given by using in Equation 14 the function
normal direction is g.j(uT, uT+1)=f((min(uT, 0)) 2 + ((max(uT+ 1,0)) 2) (15)
v = N(s, t)'X,(s, t) (8) and the appropriate (weak) entropy solution of ~b can be
In our case, we need to 'impose' v = 1. The gradient V¢ written as
is always normal to the curve given by ~b(x,y, t) = 0 so that ¢7 +' = dP7- At'g(O_ #PT,O +#PT) (16)
v¢ where
N(s, t) = -
IIV¢ll
D_ ~7- 4'7 - ¢7-
the minus sign indicating the inward direction of Ax
propagation, and hence
and
v¢
v = N'X t = - - - Xt = 1 (9) '/'7+ 1 - ~b7
IIV4,11 n+cr-
Ax
Substituting this into Equation 7 yields the surface-
This is a so-called lst-order scheme. More sophisticated
evolution equation
higher-order schemes are introduced in Reference 17. The
~,-IIV¢ll =0 (10) above, scheme is readily extended to more than one
dimension; for example, for H(u, v) = f ( u 2, V2) (in our case
Sethian called this approach Eulerian, since the coordi-
nates here are the natural physical coordinates (x,y).
U=¢x, v=e~,),
~t)n..+ l _ _
Therefore, if we have a surface ~b propagating according u - ~ bnu - A t g(D-~bu,
o x n x n .
D+¢u, y n
D-~bu, y n
D+~bu) (17)
to Equation 10 with the level set ¢(x, y, 0 ) = 0 coinciding
where
with X(s, 0), then ~b(x,y, t) = 0 produces X(s, t) propagated
according to Equations 5, while inherently solving the g.J = f((min( D~- ¢7j, 0) z + (max(D~.~b~'~,0)) z;
topological problems. To drive a numerical scheme for (min(DL ¢~'j,0))z + (max(D~ ~b~'i,0)) 2) (18)
the surface-propagation equation, we follow Reference
17, and show the connection with Hamilton-Jacobi The result is the following algorithm
methods, weak solutions and conservation laws.
Consider the 1D equation of the type of Equation 10:
Algorithm
¢ , - live II = ¢(x, t),- (¢#)1n = 0 (11)
If we define u - ¢ x and H [ u ] = - ( u 2 ) 1/2, the differenti-
(1) Choose a function ¢(x,y,0) such that
ation of the above with respect to x results in a so-called
Hamilton-Jacobi equation, in a conservation-law form: • ~b(x,y, 0) = 0 provides the initial curve X(s, 0),
• ~b(x,y,0) < 0 in the interior of the initial
u,+ [ H [ u ] ] = = 0 (12)
curve,
The 'weak' solution of the above equation is defined as • ¢(x, y, 0) > 0 in the exterior of the initial
a function u(x, t) that satisfies curve,
• ¢(x, y, 0) is Lipschitz-continuous.
dt u(x, t) dx = H[u(a, t)] - H [ u ( b , t)] (13) The next section discusses the possible ways of
'choosing' the initialization ~b(x,y, 0).
To devise a numerical scheme, define u~'= u(iAx, nat). A
differential scheme of three points is said to be in (2) Propagate ~b on an x - y grid of desired spatial
conservation form if there is a 'flow' function g(u~, us) resolution according to
such that 4,,-IIV¢,ll =o
n+l n n n n
Ui --U i
-
g(Ui,Ui+l)--~(Ui-l,U n)
(14) discretized using any conservation-form numeri-
At Ax cal scheme.
where g(u,u)=H(u) is the consistency condition. A (3) Stop after n = L / A t time steps, and find the
scheme is said to be monotone ]f " ui"+ 1 - F ( u i "- l , u i ,"u i +" l) contour (level set) ~b(x,y, L ) = 0 which is XL(s).
is an increasing monotone function of its three variables.
A basic result in numerical analysis is that a scheme which See further below possible ways of implementing the
is monotone and can be represented in a conservation level-set finder. The result is a weak solution of
form automatically obeys the entropy condition 2°. Equations 5 that obeys the entropy conditions.

volume 25 number 3 march 1993 157

R Kimmel and A M Bruckstein

,~(t2) h=Ax=Ay=C=l, then the values of the 0(x,y,0)

function on the grid varies in the interval [ - 1, 1]. The
values of the open interval ( - 1,1 ) are given only to grid
points at a distance of less than the mesh size from the
curve. This initialization problem is quite simple, and
z~/ y ...... "':::::::':::":'" "":21"-"'"
can be regarded as a problem of finding a very tight offset
neighborhood to the initial curve before topological
problems can even begin to affect the results. See Figure 4
for the results of the offsetting process with such an
initialization.
{x(s,tl),y(s,tl)} Initializing qS(x,y,0) when the shape outline is a
sequence of line segments and circular arcs (the common
NC case (see Reference 12)) is quite simple, and it can
be dealt with as follows:
• Find all the intersection points of the shape outline
with the given grid.
• For each grid point (i,j), define a cell as ~'ij= [(oij,
(oi+ l,jcki+ 1,j+ 1, (bl,j+ 1]. For each grid point (i,j), check
X
whether the boundaries of the cell JVij are intersected
Figure 3. When q5 propagates in time, the function may by the shape boundaries.
stay continuous while the offset curves form two separate
close curves which are no longer connected
i
l.tO~

The algorithm automatically enforces the entropy

condition, and frees us from the need to take care of
lee,
possible topological changes (see Figure 3). The algorithm,
by choosing forward derivatives, readily deals with shock
formation and propagation within the numerical flow. °e,

Initialization 6e,

It is obvious that an initialization of the type

4~(x, y, O) 40

+ d((x, y), X(s, 0)) (x, y) ~ exterior ofX(s, 0)

!
(x, y) e interior ofX(s, 0) I°

(x,y) ~ X(s,0)
(19)
where d is the (minimal) Euclidean distance of the point
from the curve, would be a reasonable initialization, but,
on the other hand, such an initialization is, in fact,
equivalent to solving the problem of offsetting, since we
could produce X(s, d) as the locus of all the points where
(a(x, y, O)=d.
However, to provide a proper solution, we can use the
fact that ck(x,y,O) only needs to be continuous, and
initialize dp(x,y,O) on the x - y grid as follows:
4,(x, y, 0)
min[+d((x,y),X(s,O)),C] (x,y) ~ exterior of

I X(s, 0) eo,

"i
=,~max[-d((x, y),X(s, O)), - C ] (x,y) ~interiorof
X(s,O)
(x, y) ~ X(s, 0)
(20) i
I, :. " ~o ,h ,h
where C is an arbitrary constant. If we choose Figure 4. Offsetting of simple shape on 128 × 128 grid

158 computer-aided design

 Shape offsets via level sets

we can use the gray levels of the image in order to initialize

~b(x, y, 0). For example, if the gray level of the shape is

i black (~-0), the background is white (---1), and the

boundaries pass through the gray level gray ( --- 1/2), then
we can take ~p(x,y,O)=l(x,y)-l/2, as the required
initialization, making direct use of the continuity of the
gray levels in the picture, without any extra calculations
(see Figure 6 for the Postscript halftoned version of a
gray-level test picture, and Figure 7 for the offsets of this
shape).
After initialization has been completed, the function ~b
is propagated according to the above-described algorithm
for n = L/At iteration steps. Finally, a contour finder must
be invoked to produce the resulting contour X(s, L) from
ck(x,y, L)=0. A description of the offsetting algorithm
can be found in Appendix B.

Contour finder
Following Reference 21, a simple contour finder for
X(s,L) can be generated in the following manner: for
each grid point (i,j), use the same cell definition as in
zoo
the previous section X o. Now, if m a x [ X u ] < 0 or
min[Yu] > 0, then the contour X(s,L) does not pass
|o through the cell. Otherwise, find the entrance and exit
points of ~b=0 by linear interpolation; this provides a
line segment of X(s, L) belonging to the contour. The line
60
segments need to be neither ordered nor directed in the
same direction to display the desired contour (see Figure
qo.
8); however, using additional information such as the
knowledge of interior and exterior points, one can
produce any desired representation of the curve, such as
Jo, a planar polygon, or a curve comprising cubic or any
other polynomial arcs.

I, I, '.. J, 11, do

Figure 5. Offsetting polygons

• If there is an intersection, attach a value to the two

grid points on the two sides of the intersection point
so that the interior grid point is given a value within
the interval ( - 1, 01 and the exterior grid point a value
within [0, + 1). These two values are samples of a linear
approximation of the ~) function whose zero is the
intersection point.
• Set values for the rest of the grid points as follows:
assign - 1 to the interior points and + 1 to the exterior
points of the given shape. (Chains of 'dependent' points
may occur. When this happens, attach values to the
dependent points so that the intersection points are
the zeros of the linear ~b approximation). Figure 5
shows an example of polygon offsetting using this
initialization process.
Note that, if we must offset a shape that is provided as Figure 6. Miekey Mouse in gray-level picture (128 x 128
a dark region on a light background in a given picture, pixels)

volume 25 number 3 march 1993 159

R Kimmel and A M Bruckstein

CONCLUDING REMARKS
We have described a new method of generating shape
offsets so that topological problems are inherently
avoided, and shocks, cusps and other singularities are
also readily dealt with in an efficient numerical scheme.
The algorithm works on an image grid with a resolution
chosen according to the desired accuracy. It is easy to
implement the algorithm in parallel using each mesh
point as a small calculating device which communicates
with its four close neighbors. In each iteration, we need
to calculate the values of c~(x,y,t) for those grid points
close to the current contour, and the rest of the grid
points serve as sign holders. This can be exploited to
reduce calculation effort.
In summary, we propose to introduce to the CAD/CAM
field some recent advances in numerical methods for fluid
~, I0 " I, d, J, dynamics. We have shown that wavefront-propagation
methods in fluid dynamics also provide a new approach
to the problem of shape offsetting.

ACKNOWLEDGEMENTS
We would like to thank Professor James Sethian for
O
providing us with his reports on numerical methods of
front evolution. We thank Professor Moshe Israeli,
Professor Ben Kimia, Professor Allen Tannenbaum and
Professor J Rubinstein for their help in introducing us
•to the world of numerical analysis and moving surfaces.
R Kimmel's work was partly supported by the
Ollendorf Fund. A Bruckstein's work was supported in
part by the Fund for the Promotion of Research at the
Technion, Haifa, Israel.

. ,~, ,I, REFERENCES

I Blum, H and Nagel, R N 'Shape description using
Figure 7. Offsetting Mickey Mouse (grid size 128)
symmetric axis features' Pattern Recognition Vol 10
(1978) pp 167-180
2 Calabi, L and Hartnett, W E 'Shape recognition,
Cell {i, j }
prairie fires, convex deficiencies and skeletons' Amer.
Math. Mon. Vol 75 (1968) pp 335-342
3 Ho~ehek, J 'Offset curves in the plane' Comput.-
Aided Des. Vol 17 No 2 (1985)
4 Klan, R 'An offset spline approximation for plane
cubic splines' Comput.-AidedDes. Vol 15 No 5 (1983)
pp 297-299
5 Tiller, W and Hanson, E G 'Offset of two dimensional
profiles' 1EEE Trans. Comput. Graph. & Applic. Vol
4 No 9 (1984) pp36-46
6 Ceqmllart, S 'Computing offsets of B-spline curves'
Comput.-AidedDes. Vo119 No 6 (1987) pp 305-309
7 Pham, B 'Offset approximation of uniform B-splines'
Comput.-AidedDes. Vo120 No 8 (1988) pp 471-474
i~,j.1, ~i.1,#.~
8 Eiher, G and Cohen, E 'Error bounded variable
Figure 8. Contour finder." find a line segment in .At u distance offset operator for free form curves and

160 computer-aided design

 Shape offsets via level sets

surfaces' Int. J. Comput. Geom. & Applic. Vol I No 1 Sethian, J A 'Numerical methods for propagating fronts'
(1991) pp 67-78 in Concus, P and Finn, R (Eds.) Variational Methods for
Free Surface Interfaces Springer-Verlag, USA (1987)
9 Farin, G 'Curvature continuity and offsets of
piecewise conics' A C M Trans. Graph. Vol 8 No 2
(1989) pp 89-99
10 Meek, D S and Walton, D S 'Offset curves of
clothoidal splines' Comput.-Aided Des. Vol 22 No 4 APPENDIX A
(1990) pp 199-200
11 Wang, R and Jiang, W 'An algorithm of the Normal evolution
offset shape' Comput. Graph. Vol 15 No 3 (1991) In this section, we show that, in a curve evolution
pp 435--439 described by
12 Held, M On the Computational Geometry of Pocket
Machining Springer-Verlag, Germany (1991) -- X(s, t) = N(s, t)
~t
13 Saeed, S E O, de Pennington, A and Dodsworth, J R the direction of the normal does not change in time
'Offsetting in geometric modeling' Comput.-Aided prior to shock formation (following References 17, 22
Des. Vol 20 No 1 (1988) pp 67-74 and 23).
14 Riazanoff, S, Cervelle, B and Chorowicz, J 'Para-
metrisable skeletonization of binary and multi-level Proof." Let us first define the metric along the curve and
images' Pattern Recognition Lett. No 11 (1990) the arc-length parameter 5:
pp 25-33
g(s,t)- ~, X =(x,2 +Ys)
2 1/2
15 Yao, C and Rokne, J 'Fat curves' Comput. Graph.
Forum No 10 (1991) pp 237-248
16 Vincent, L 'Exact Euclidean distance function by
chain propagation' IEEE Trans. Comput. V. Pattern
fo
Now we can define the tangent, curvature and the
Recognition (May 1991)
normal in a standard way:
17 Osher, S J and Sethian, J A 'Fronts propagating
with curvature dependent speed: algorithms based r_ a x=l__ax
on Hamilton-Jacobi formulations' J. Comput. Phys. c~ g c~s
Vol 79 (1988) pp 12-49
18 Sethian, J A 'Curvature and the evolution of fronts'
Comm. Math. Phys. Vol 101 (1985) pp 487-499
19 Sethian, J ~A'A review of recent numerical algorithms --T
for hypersurfaces moving with curvature dependent N- t~ _ 1 _~ T
speed' J. Differ. Geom. No 33 (1989) pp 131-161
f~sT kg &
20 Sod, G A Numerical Methods in Fluid Dynamics
Cambridge University Press (1985) It is easily established that the curve evolution induces
21 Sethian, J A and Strain, J 'Crystal growth and
dendritic solidification' J. Comput. Phys. Vol 98 -- T = - kgN
&
(1992)
22 Kimia, B B, Tannenbaum, A and Zucker, S W 'On -- N = kgT
the evolution of curves via a function of curvature. ~s
Part 1: The classical case' J. Math. Anal. & Applic. The evolution of g can be computed as follows:
Vol 163 (1992) pp 438-458
23 Kimia, B B 'Toward a computational theory of -g\N , x
shape' PhD Thesis Dep. Electrical Engineering,
McGill University, Canada (1990)
= /

BIBLIOGRAPHY
Mitchel, A R and Grifliths, D F The Finite Difference = 2(gT, kgT)
Method in Partial Differential Equations John Wiley,
USA (1985) =2g(kg)

volume 25 number 3 march 1993 161

R Kimmel and A M Bruckstein

Hence, APPENDIX B

Offsetting algorithm
dr
(1) Initialize O°=c~(ih,jh, O), where h = A x = A y ,
Then, the evolution of the tangent is as follows:
according to the curve representation in hand.
~b(x, y, 0) = 0 is the implicit representation of the
ot =ot initial curve. See above for ways in which to
determine O(x,y,O). Note that, if shape is
- g' ex+-l£- x defined by a dark region on a white back-
g 2 dS g ~t C3S ground, set c~(x,y,O) equal to the gray-level
image of the object at the camera resolution.
=_g'T+!~N (2) For n = 1 to L/At do
g g c~s
~.id+ 1 = ~.,j + zxt . . .
~ {(mm[(~bl,j - ~bi-, j), O] )2
_ kgT+l_kgT
g g + (max[(~bT+ x,j- ~bTj),0]) 2
=0 + (min[ (~Tj- ckT,j-~ ), 0]) 2
Thus, the normal (to the stationary tangent direction) is n 0 ]) 2 } 1/2
+(max[(d)Tj+l-~b,.j),
also stationary; in other words, there is no change in
(3) Find the level set ~b(x, y, L) = 0 by activating the
the normal direction while the curve evolves, and no
contour-finder procedure as described above on
shocks appear. Therefore, Equations ~1are equivalent to
Equations 5. []

162 computer-aided design