Numerical Implementation

Equations (11a) and (11b) can be solved by treating u and v as functions of time and solving   The steady-state solution (as $t\rightarrow \infty$) of these linear parabolic equations is the desired solution of the Euler equations (11a) and (11b). Note that these equations are decoupled, and therefore can be solved as separate scalar partial differential equations in u and v. The equations in (12) are known as generalized diffusion equations, and are known to arise in such diverse fields as heat conduction, reactor physics, and fluid flow [30]. For us, they have appeared from our description of desirable properties of external fields for active contours. Diffusion is a natural outcome given the desired ``filling in'' property.

For convenience, we rewrite Equation (12) as follows where

$$\begin{eqnarray*}b(x,y) & = & {f_x}(x,y)^2 + {f_y}(x,y)^2 \\ c^1(x,y) & = & b(x,y)f_x(x,y) \\ c^2(x,y) & = & b(x,y)f_y(x,y) \end{eqnarray*}$$

Any digital image gradient operator (cf. [31]) can be used to calculate f_x and f_y. In the examples shown in this paper, we use simple central differences. The coefficients b(x,y), c¹(x,y), and c²(x,y), can then be computed and fixed for the entire iterative process.

To set up the iterative solution, let the indices i, j, and n correspond to x, y, and t, respectively, and let the spacing between pixels be $\Delta x$ and $\Delta y$ and the time step for each iteration be $\Delta t$. Then the required partial derivatives can be approximated as

$$\begin{eqnarray*}u_t & = & \frac{1}{\Delta t}(u_{i,j}^{n+1} - u_{i,j}^{n})\\ v... ...x\Delta y}(v_{i+1,j} +v_{i,j+1}+v_{i-1,j}+v_{i,j-1} - 4v_{i,j}) \end{eqnarray*}$$

Substituting these approximations into (13) gives our iterative solution to GVF: where

$$r = \frac{\mu\Delta t}{\Delta x\Delta y}$$ (11)

Convergence of the above iterative process is guaranteed by a standard result in the theory of numerical methods (cf. [32]). Provided that b, c¹, and c² are bounded, (14) is stable whenever the Courant-Friedrichs-Lewy step-size restriction $r \le 1/4$ is maintained. Since normally $\Delta x$, $\Delta y$, and $\mu$ are fixed, using the definition of r in (15) we find that the following restriction on the time-step $\Delta t$ must be maintained in order to guarantee convergence of GVF:

 

$$\Delta t \le \frac{\Delta x \Delta y}{4 \mu}$$ (12)

The intuition behind this condition is revealing. First, convergence can be made to be faster on coarser images -- i.e., when $\Delta x$ and $\Delta y$ are larger. Second, when $\mu$ is large and the GVF is expected to be a smoother field, the convergence rate will be slower (since $\Delta t$ must be kept small).