EE 396: Lecture 12 Ganesh Sundaramoorthi March 26, 2011

1

Level Set Methods

Last lecture, we derived an algorithm for minimizing the functional E({ai , Ri }N i=1 )

= − log

p({ai , Ri }N i=1 |I)




Z N X 1 (I(x) − ai )2 dx + αL(∂Ri ). = 2σn2 Ri

(1)

i=1

The algorithm is based on an alternative minimization scheme where we first fix the regions Ri and solve for the optimal constants ai , and then fixed the constants ai and then perturb the boundaries of ∂Ri (and hence the regions Ri ) in the steepest (i.e., gradient direction) so as to decrease the energy. The gradient descent for E with respect to Ri results in a curve evolution describing the motion of the boundary of Ri in the form : ∂ c(p, t) = f (c(p, t))N (p, t) + ακ(p, t)N (p, t). ∂t

(2)

The straightforward way to implement this evolution numerically is to discretize the curve, and then discretize the equation above : ck+1 (ski ) − ck (ski ) = f (ck (si ))N k (ski ) + αckss (ski ) ∆t

(3)

where ck (ski+1 ) − ck (ski−1 ) , (J is rotation by 90) |ck (si+1 ) − ck (si−1 )| (ck (ski+1 ) − ck (ski ))/∆ski − (ck (ski ) − ck (ski−1 ))/∆ski−1

Nsk (ski ) = Jcks (ski ) = J ckss (si ) =

1 k 2 (∆si

ck (si ) = ith sample point

+ ∆ski−1 )

(4) (5) (6)

∆ski = |ck (ski+1 ) − ck (ski )|.

(7)

Note that as the curve grows or shrinks, we would need to resample the curve so that points on the curve do not become to close or too far from each other. We would also devote considerable programming to deal with the region represented by c merging or splitting. While this is certainly possible, it is quite cumbersome using this representation of the curve. This motivates using a technique known as level set methods [1], which is a way of implementing a curve/surface evolution. The method is not just restricted to curve evolutions in the form (2), indeed, it applies to many other curve/surface evolutions provided they satisfy some important properties that we discuss below. 1

1.1

Implicit Representation of a Curve

First we recall an implicit representation of a simple closed curve : c : S 1 → R2 . An implicit representation of c is a function Ψ : R2 → R such that {x ∈ R2 : Ψ(x) = 0} is the curve c, that is, the level set of Ψ is zero. For example, the implicit representation of a circle :   R cos (2πp) c(p) = , p ∈ [0, 1] (8) R sin (2πp) is Ψ(x1 , x2 ) = x21 + x22 − R2 .

(9)

More generally, a canonical way of constructing an implicit representation of a generic simple connected closed curve c is via the signed distance function. We first define the distance function of c as dc : Ω → R as dc (x) = inf |x − c(p)|, (10) p∈[0,1]

i.e., each point x ∈ Ω, the domain of dc is assigned the minimum distance to the curve c. Now the signed distance function is ( −dc (x) x ∈ int(c) Ψ(x) = (11) dc (x) x∈ / int(c) where int(c) denotes the region enclosed by the curve c. Note that Ψ(x) = 0 when x in on the curve c and nowhere else. Thus, Ψ is an implicit representation of c. The signed distance function is usually used in level set methods, which we describe below, for reasons of numerical stability, but indeed any implicit representation may be used in level set methods. Although it seems costly to compute the signed distance function, in the next lecture, we will describe very quickk algorithms to compute the signed distance function.

1.2

Curve evolution as an implicit function evolution

In level set methods, instead of evolving the curve directly, one evolves an implicit representation of the curve with the property that the zero level set of the implicit function always corresponds to the curve of interest. In otherwords, the evoution of the zero level set of the implicit function should correspond to the evolution of the curve in level set methods. Mathematically, we define an implicit function that changes in time t : Ψ : R+ × R2 → R where the arguments of Ψ are (t, x) ∈ R+ × R2 . We want that the zero level set of Ψ at each time t corresponds to the curve c at time t, i.e., Ψ(t, c(t, p)) = 0, for all p ∈ [0, 1]

(12)

where the curve evolves according to ∂ c(p, t) = F (p, t) (13) ∂t where F is some deformation field of the curve. Now given these two conditions, we calculate the evolution of Ψ by taking the time derivative of (12) : d ∂ Ψ(t, c(t, p)) = Ψt (t, c(t, p)) + ∇Ψ(t, c(t, p)) · c(p, t) dt ∂t = Ψt (t, c(t, p)) + ∇Ψ(t, c(t, p)) · F (p, t),

0=

2

(14) (15)

where Ψt is the derivative with respect to the first argument of Ψ and ∇Ψ is the derivative w respect to the second argument. We rewrite the above equation as Ψt (t, x) = −∇Ψ(t, x) · F (p, t), x = c(t, p) for some p

(16)

and thus we know how the implicit function evolves for points x ∈ R2 that are on the curve c. Note that F does not have meaning outside the curve c, and thus, we need to specify how to define F for points x ∈ R2 that are not on the curve. In general, this is done through extensions, which we will briefly address later. Now, it seems that we are expending much more effort implementing the evolution of the implicit representation than the curve directly; indeed instead of evolving a one-dimesional curve, we are now evolving a 2-dimensional function Ψ. However, if done correctly, one can implement the implicit evolution with nearly the same computational complexity as directly discretizing the curve evolution. The advantages of using the level set method are • resampling the curve when the curve grows/shrinks is taken care of automatically • topological changes of the curve, e.g., when the curve splits into two or more curves or when two or more curves merge together are handled naturally without any additional programming effort the disadvantage of the method is that not all curve evolutions may be implemented via level set methods, indeed, the curve evolution must satisfy a Comparison Principle. That is, if two curves one inside the other are evolving according to the curve evolution, then (13) the curves must not collide. If they were colliding, then level sets of Ψ for different levels would be colliding and then |∇Ψ| → ∞.

1.3

Level Set Implementation of Curvature Flow

We show how to find the corresponding level set evolution of curvature flow: ct = css = κN.

(17)

Let us first express the normal vector of the curve in terms of the implicit representation. Note that differentiating (12) with respect to the arclength parameter s, we have that 0=

∂ ∂ Ψ(t, c(t, s)) = ∇Ψ(t, c(t, s)) · c(t, s) = ∇Ψ(t, c(t, s)) · T (t, s) ∂s ∂s

(18)

where T is the tangent vector to the curve. Therefore, we see that ∇Ψ is perpendicular to the tangent of the curve and thus ∇Ψ(t, c(t, s)) N (t, s) = (19) |∇Ψ(t, c(t, s))| is normal to the curve. If we use the convention that the value of Ψ is negative inside the curve and positive outside (as in the case of the signed distance function), then the above is the outward normal vector. Now the level set evoluton becomes ∇Ψ(t, x) )) = κ(t, x)|∇Ψ(t, x)|, |∇Ψ(t, x)| (20) where the evolution is defined for points x on the curve c. We can express the curvature of the curve in terms of the level set function Ψ (and indeed we have already done it in Lecture 6). Recall that   d ∇Ψ(c(s)) ∇Ψ(c(s)) J = div κ(s) = N (s) · css (s) = N (s) · . (21) ds |∇Ψ(c(s))| |∇Ψ(c(s))| Ψt (t, x) = −∇Ψ(t, x) · (κ(t, x)N (t, x)) = −∇Ψ(t, x) · (κ(t, x) · (−

3

Computing the expression above, we find that κ(x) =

Ψx1 x1 (x)Ψ2x2 (x) − 2Ψx1 x2 Ψx1 Ψx2 (x) + Ψx2 x2 (x)Ψ2x1 (x) . |∇Ψ(x)|3

(22)

Note the above has meaning even for points x not on the curve; indeed recall from Lecture 6 that the expression above is the curvature of the level set of Ψ that passes through the point x. Thus, there is a natural extension of κN beyond the curve c. Thus, the level set evolution for curvature flow is   ∇Ψ(t, x) Ψt (t, x) = |∇Ψ(t, x)| div . (23) |∇Ψ(t, x)| Note that curvature flow does satisfy a comparison principle (as discussed in Lecture 10-11), and thus level sets at different levels will not collide and thus the gradient of Ψ does not blow up. In the equation above, each level set of Ψ (including the zero level set) are simulataneously evolving by curvature flow. To discretize the above equations, one uses central differences for all the spatial differences and a forward difference for the time derivative.

1.4

Level Set Implementation of Region-Based Term

We now consider level set implementation of the curve evolution : ct (t, p) = f (c(t, p))N (t, p)

(24)

where f : Ω → R. Using that Ψt = −∇Ψ · ct , we see that Ψt (t, x) = f (x)|∇Ψ(t, x)|, x = c(t, p) for some p.

(25)

Note that f really only has meaning on the curve c (although it is defined on all Ω), and thus we would need to extend f to points x that are not on the curve. A typical way to do this is to define f at x to be f (y) where y is the closest point on the curve c to x 1 . We call the extension of f , fext . The level set equation above behaves very differently than the level set equation for curvature flow, and indeed, if you were to use central differences for the equation, the equation would be unstable. Thus, it is necessary to analyze the behavior of the PDE to understand the appropriate discretizing scheme. To gain intuition, let us look at the one-dimensional analog of (25) : ( fext (x)Ψx (t, x) Ψx (t, x) > 0 Ψt (t, x) = fext (x)|Ψx (t, x)| = (26) fext (x)Ψx (t, x) Ψx (t, x) < 0 where x ∈ R. Let us also suppose for the moment that fext does not depend on x, and let −b = fext . The equation Ψt (t, x) = −bΨx (t, x) (27) is called a transport equation, and thus, (26) is a coupled pair of transport equations. Note that the previous transport equation can be written Ψt (t, x) + bΨx (t, x) = 0 or

∂ Ψ(t, x) = 0, where v = (1, b); ∂v

1

(28)

In level set methods, one usually does not evolve x for each x ∈ Ω, indeed, in practice, one only evolves Ψ is a narrowband around the curve c (since ultimately we are only interested in how the zero level set evolves). This reduces computational cost significantly. Typically the narrowband is about 4 pixels around the zero level set. As long as the narrowband is thin, the extension defined above is well defined in the sense that there is a unique point y on the curve closest to x.

4

∂ denotes the directional derivative in the direction v. This last expression says the function Ψ is note that ∂v constant in the direction v in the (t, x) space. Thus, along lines that are parallel to v the function remains constant. Such curves where the solution of a PDE are constant are called characteristics of the PDE. Suppose that we have an initial condition for Ψ at time t = 0, that is, Ψ(0, x) = Ψ0 (x) where Ψ0 : R → R. Then we may propogate this information along lines parallel to v to obtain Ψ(t, x) for any t > 0 and x ∈ R. We simply follow a line parallel to v starting from (t, x) and follow the line back to time zero and the corresponding spatial location y, and then Ψ(t, x) = Ψ(0, y) = Ψ0 (y). More precisely, the line through (t, x) and parallel to v is l(s) = (t, x) + sv = (t + s, x + bs) where s ∈ R, and we seek to find the s when l(s) = (0, y). Thus, we see that s = −t and y = x + bs = x − tb. Hence, we see that

Ψ(t, x) = Ψ(0, x − tb) = Ψ0 (x − tb).

(29)

We see that information (i.e., the initial condition Ψ0 ) flows in the direction v, and thus, when we discretize the transport equation, we must make sure that we use only information that is in the opposite direction of the information flow (See Diagram in Class). Hence the discretization of the transport PDE Ψt = −bΨx is ( −bDx+ Ψ(t, x) b < 0 Ψ(t + ∆t, x) − Ψ(t, x) = . (30) ∆t −bDx− Ψ(t, x) b > 0 where Ψ(t, x + ∆x) − Ψ(t, x) ∆x Ψ(t, x) − Ψ(t, x − ∆x) Dx− Ψ(t, x) = ∆x Dx+ Ψ(t, x) =

(31) (32)

are forward and backward differences respectively. If one were to perform Von-Neumann analysis then the stable step size is ∆t 0 ≤ |b| ≤ 1. (33) ∆x In the case of the equation ( bΨx (t, x) Ψx (t, x) > 0 Ψt (t, x) = b|Ψx (t, x)| = , (34) −bΨx (t, x) Ψx (t, x) < 0 the discretization scheme is   bDx+ Ψ(t, x) b > 0, Dx+ Ψ(t, x) > 0    b < 0, Dx− Ψ(t, x) > 0 Ψ(t + ∆t, x) − Ψ(t, x) bDx− Ψ(t, x) =  ∆t −bDx− Ψ(t, x) b > 0, Dx− Ψ(t, x) < 0    −bD+ Ψ(t, x) b < 0, D+ Ψ(t, x) < 0 x x ( max (Dx+ Ψ(t, x), 0) − min (Dx− Ψ(t, x), 0) b > 0 =b . max (Dx− Ψ(t, x), 0) − min (Dx+ Ψ(t, x), 0) b < 0

(35)

(36)

In class, we have analyzed this discretization scheme with examples Ψ0 (x) = |x| and Ψ0 (x) = −|x| the former produces a shock and the latter produces a rarefaction fan. 5

Getting back to the equation of interest : Ψt (t, x) = fext (x)|∇Ψ(t, x)| = fext (x)

q Ψ2x1 (x) + Ψ2x2 (x),

(37)

and applying the reasoning of the one-dimensional equation, we arrive at the following discretizing scheme : v u   u max2 (Dx+ Ψ(t, x), 0) + min2 (Dx− Ψ(t, x), 0)  1 1  t fext (x) > 0   2 (D + Ψ(t, x), 0) + min2 (D − Ψ(t, x), 0)  + max x2 x2 Ψ(t + ∆t, x) − Ψ(t, x) = fext (x) v . u  ∆t  u min2 (Dx+ Ψ(t, x), 0) + max2 (Dx− Ψ(t, x), 0)  1 1  t  fext (x) < 0   + min2 (Dx+2 Ψ(t, x), 0) + max2 (Dx−2 Ψ(t, x), 0) (38) For additional information see, for example, [2].

References [1] S. Osher and J.A. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. Journal of computational physics, 79(1):12–49, 1988. [2] J.A. Sethian. Level set methods and fast marching methods. 1999.

6