## Numerical Solutions

The numerical solution for region force diffusion is discussed in detail in Section 10.5.1, but the detailed numerical solutions for RAGS level set representation are only presented in Appendix A as they are not critical to the understanding of the concepts underlying RAGS. In fact, the whole of this section can be skipped without loss of continuity.

10.5.1 Numerical Solutions for Region Force Diffusion for RAGS

Initially, a mesh grid needs to be selected, with final accuracy directly dependent on its resolution. However, due to the nature of a digital image, the grid resolution is constrained to the pixel level. It was shown in Section 10.4.2 that the steady solution of (10.10) can be achieved by computing the equilibrium state of (10.12). The initial state of the region force vector field R is given by the gradient of the region map R. Simple central differences can be used to approximate VR, resulting in vectors that are then diffused. Let Ax and Ay be the grid spacing, At be the time step, and i, j, and n represent the spatial position and time. The partial derivative of time can be approximated by forward difference as yK+1 _ yU

The spatial partial derivatives can be solved using central differences approximation given by

V2u _ ui+1,j + ui,j+1 + ui-1,j + ui^.j—1 ~ 4ui,j (10 24)

AxAy

The solutions to partial derivatives of v(x, y, t) are similar to those of u(x, y, t). The weighting functions given in (10.11) can be easily computed. Thus, substituting the partial derivatives into (10.12) gives the following iterative solution:

where,

A — ââ (Vn+i'j+j+j+uti'j - 4uïj) - q (0ij - Rx'ij)

and a = AAXAj(vi+ij + j + j + <-ij - Kj) - q(Oij(vlj - Ryij)], where Rxij and Ryij are partial derivatives of R. They can also be approximated by central differences as

The convergence is guaranteed with the time step restriction of (10.13).

10.5.2 Numerical Solution for the Level Set Implementation of RAGS

As in the numerical solution for vector diffusion, a computational grid is required. Once the grid is chosen, the initial level sets 0(x, t) = 0 can be defined with the property that the zero level set corresponds to the initial contours of the snake. The signed-distance transform can be used to build the initial level sets. A brute-force Euclidean distance transform would be computationally infeasi-ble. Practically, accuracy is required only near the initial contours, and discrete values based on grid distance can suffice further away. A positive sign is given to the points outside the contours, and a negative sign is applied to the points inside.

As shown in (10.17), the snake evolves according to four forces that can be categorized into three types based on the nature of their motions:

1. The first motion is of a smoothing and collapsing nature with speed proportional to its curvature as shown in Fig. 10.1. It can be numerically approximated using central differences, because the curvature is only dependent on the contour. It is independent of time and spatial position.

2. The second is expanding or shrinking with a spatially constant speed, characterized by ag( ) in the normal direction of the curve. However, when the constant term exists, the normals can collide with each other while evolving. Thus shocks, or corners, will form and once a shock has developed, some information will be lost as it evolves. This means that shocks cause irreversibility; information cannot be recovered by tracing 'backwards' in time. Generally, no new information can be created while evolving, which is referred to as an entropy condition. Central difference approximation cannot be used to approximate the gradient in this case, as it suffers from shocks where the entropy condition is invoked. An upwind scheme can be used, as an entropy-satisfying scheme, that engages information upwind of the direction of its propagation. In other words, in order to achieve a stable numerical scheme, the numerical domain of dependence should contain the mathematical domain of dependence. Thus, in order to approximate the gradient of the constant term, it is important to first know which way the speed function points, and whether it is negative or positive. Then we can choose proper backward or forward difference approximations.

3. The third type of motion in (10.17) is contributed by the underlying static velocity field, the direction and strength of which are based on spatial position. It is independent of the shape and position of the snake. The motion of contours under this velocity field can be numerically approximated through upwind schemes by checking the sign of each component of the velocity field and constructing one-sided upwind differences in the appropriate direction. For a positive speed component, backward difference approximation is used, otherwise forward difference approximation should be applied.

By using these approximation schemes, (10.17) can be numerically implemented. The detailed numerical solutions for RAGS are presented in Appendix A. For general numerical solution to level sets, including concepts such as entropy condition and upwind scheme, the interested reader is referred to works by Sethian [16,17] and by Osher et al. [18].