## Semiimplicit CoVolume Scheme

Get Instant Access

We present our method in discretization of Eq. (11.8), although we always use its e-regularization (11.2) with a specific e > 0. The notation is simpler in the case of (11.8) and it will be clear where regularization appears in the numerical scheme.

First we choose a uniform discrete time step t and a variance a of the smoothing kernel Ga. Then we replace time derivative in (11.8) by backward difference. The nonlinear terms of the equation are treated from the previous time step while the linear ones are considered on the current time level, this means semi-implicitness of the time discretization. In the last decade, semi-implicit schemes have become a powerful tool in image processing, we refer e.g. to [3,4,25-27,33,37,51,57,58].

Semi-implicit in time discretization. Let t and a be fixed numbers, I0 be a given image, and u0 be a given initial segmentation function. Then, for

Figure 11.15: Segmentation of the circle with a big gap (Fig. 11.12 (middle)) using e = 1 (top), e = 10-2 (middle), and e = 10-5 (bottom). For bigger missing part a smaller e is desirable. In the left column we see how close to the edges the isolines are accumulating and closing the gap, and in the right we see how steep the segmentation function is along the gap (color slide).

Figure 11.15: Segmentation of the circle with a big gap (Fig. 11.12 (middle)) using e = 1 (top), e = 10-2 (middle), and e = 10-5 (bottom). For bigger missing part a smaller e is desirable. In the left column we see how close to the edges the isolines are accumulating and closing the gap, and in the right we see how steep the segmentation function is along the gap (color slide).

Figure 11.16: Isolines of the segmentation function in the segmentation of the noisy circle (Fig. 11.12 (right)) are shown in time steps 0, 50, 100, and 200. Since the gap is not so big we have chosen s = 10-1 (color slide).

n = 1,..., N, we look for a function un, solution of the equation,

A digital image is given on a structure of pixels with rectangular shape, in general (red rectangles in Fig. 11.18). Since discrete values of 10 are given in pixels and they influence the model, we will relate spatially discrete approximations of the segmentation function u also to image pixels, more precisely, to their centers (red points in Fig. 11.18). In every discrete time step of the method (11.11), we have to evaluate gradient of the segmentation function at the previous step |Vun-1|. For that goal, it is reasonable to put a triangulation (dashed black lines in Fig. 11.18) inside the pixel structure and take a piecewise linear approximation of the segmentation function on this triangulation. Such an approach will give a constant value of the gradient per triangle, allowing simple and clear construction of fully discrete system of equations. This is the main feature of the co-volume [25,56] and finite element [13-15] methods in solving mean curvature flow in the level set formulation.

Figure 11.17: The graph of the segmentation function and its histograms in time steps 100 and 200 for the same experiment as presented in Fig. 11.16. The histograms give a practical advise to shorten the segmentation process in case of noisy images. For a noisy image, the formation of completely piecewise flat subjective surface takes longer time. However, the gaps in histogram of the segmentation function are developed soon. It allows to take any level inside these gaps and to visualize the corresponding level line to get desirable segmentation result (color slide).

Figure 11.17: The graph of the segmentation function and its histograms in time steps 100 and 200 for the same experiment as presented in Fig. 11.16. The histograms give a practical advise to shorten the segmentation process in case of noisy images. For a noisy image, the formation of completely piecewise flat subjective surface takes longer time. However, the gaps in histogram of the segmentation function are developed soon. It allows to take any level inside these gaps and to visualize the corresponding level line to get desirable segmentation result (color slide).

As can be seen in Fig. 11.18, in our method the centers of pixels are connected by a new rectangular mesh and every new rectangle is splitted into four triangles. The centers of pixels will be called degree of freedom (DF) nodes. By this procedure we also get further nodes (at crossing of red lines in Fig. 11.18) which, however, will not represent degrees of freedom. We will call them non-degree of freedom (NDF) nodes. Let a function u be given by discrete values in the pixel centers, i.e. in DF nodes. Then in additional NDF nodes we take the average value of the neighboring DF nodal values. By such defined values in NDF nodes, a piecewise linear approximation uh of u on the triangulation can be built. Let us note that we restrict further considerations in this chapter only to this type of grids. For triangulation Th, given by the previous construction, we construct a co-volume (dual) mesh. We modify a basic

 ' 1 \ s 1 * r 1 S 1 \ / \ 1 ' V ' f ---_ _ _ _ _ _ N 1

Figure 11.18: The image pixels (red solid lines) corresponding to co-volume mesh. Triangulation (black dashed lines) for the co-volume method with degree of freedom nodes (red round points) corresponding to centers of pixels (color slide).

approach given in [25, 56] in such a way that our co-volume mesh will consist of cells p associated only with DF nodes p of Th, say p = 1,...,M. Since there will be one-to-one correspondence between co-volumes and DF nodes, without any confusion, we use the same notation for them. In this way we have excluded the boundary nodes (due to Dirichlet boundary data) and NDF nodes.

For each DF node p of Th, let Cp denote the set of all DF nodes q connected to the node p by an edge. This edge will be denoted by apq and its length by hpq. Then every co-volume p is bounded by the lines (co-edges) epq that bisect and are perpendicular to the edges apq, q e Cp. By this construction, the co-volume mesh corresponds exactly to the pixel structure of the image inside the computational domain ^ where the segmentation is provided. We denote by Epq the set of triangles having apq as an edge. In a situation depicted in Fig. 11.18, every Epq consists of two triangles. For each T e £pq let cpq be the length of epq that is in T, i.e., cpq = m(epqi the portion of epq that is in T, i.e., cpq = m(epqPT), where m is a measure in

R . Let Np be the set of triangles that have DF node p as a vertex. Let uh be a piecewise linear function on triangulation Th. We will denote a constant value of \Vuh\ on T e Th by |VuT | and define regularized gradients by

We will use the notation up = uh(xp), where xp is the coordinate of the node p of triangulation Th.

With these notations, we are ready to derive co-volume spatial discretization. As is usual in finite volume methods [20,34,44], we integrate (11.11) over every co-volume p, i = 1,..., M. We get f 1 un - un-1 f ( 0 Vun \

For the right-hand side of (11.13), using divergence theorem we get f ( 0 Vun \ f g0 dun V-l g0-r )dx = / ,--ds

Xf ncCI

So we have an integral formulation of (11.11)

expressing a "local mass balance" property of the scheme. Now the exact "fluxes" on the right-hand side and "capacity function" on the left-hand side (see e.g. [34]) will be approximated numerically using piecewise linear reconstruction of un-1 on triangulation Th. If we denote g° approximation of g0 on a triangle T e Th, then for the approximation of the right-hand side of (11.14), we get yiy cTvq| , (11.15)

h\Tk.n IVuni hr, ' 1 j and the left-hand side of (11.14) is approximated by un - up 1

where m(p) is a measure in lRd of co-volume p and either

The averaging of the gradients fll.17) has been used in [25, 56], and the approximation (11.18) is new and we have found it very useful regarding good convergence properties in solving the linear systems (see below) iteratively for e ^ 1. Regularizations of both the approximations of the capacity function are as follows: either

Now we can define coefficients, where the e-regularization is taken into account, namely, bn— = Mpmfp), (11.21)

hpq Te£„ |V UT !e which together with (11.15) and (11.16) give the following.

Fully-discrete semi-implicit co-volume scheme. Let u°p, p = 1,..., M, be given discrete initial values of the segmentation function. Then, for n = 1,..., N we look for up, p = 1,..., M, satisfying bnpU + t E ann \unp - up = bn- up-1. (11.23)

qeCp

Theorem. There exists a unique solution (W^, ..., uM) of the scheme (11.23) for any t > 0, e > 0 and for every n = 1,..., N. Moreover, for any t > 0, e > 0 the following stability estimate holds min up < min up < max up < max u 1 < n < N. (11.24)

Proof. The system (11.23) can be rewritten in the form

Applying Dirichlet boundary conditions, it gives the system of linear equations with a matrix, the off diagonal elements of which are symmetric and negative. Diagonal elements are positive and dominate the sum of absolute values of the nondiagonal elements in every row. Thus, the matrix of the system is symmetric and diagonally dominant M-matrix which imply that it always has a unique solution. The M-matrix property gives us the minimum-maximum principle, which can be seen by the following simple trick. We may temporarily rewrite (11.23) in the equivalent form

bP q eCp and let max(w™,..., uM) be achieved in the node p. Then the second term

< uWp on the left-hand side is non-negative and thus max(up,..., uM) = up < up 1 <

max(un 1,...,uM x). In the same way we can prove the relation for minimum and together we have min unp 1 < min unp < max unp < max unp , 1 < n < N, (11.27)

p p p p p p p p which by recursion imply the desired stability estimate (11.24). □

So far, we have said nothing about evaluation of gT included in coefficients (11.22). Since image is piecewise constant on pixels, we may replace the convolution by the weighted average to get := Ga * I0 (see e.g. [37]) and then relate discrete values of I® to pixel centers. Then, as above, we may construct its piecewise linear representation on triangulation and in such way we get constant value of Von every triangle T e Th. Another possibility is to solve numerically a linear heat equation for time t corresponding to variance a with initial datum given by I0 (see e.g. [3]). The convolution represents a preliminary smoothing of the data. It is also a theoretical tool to have bounded gradients and thus a strictly positive weighting coefficient g0. In practice, the evaluation of gradients on discrete grid (e.g., on triangulation described above) always gives bounded values. So, working on discrete grid, one can also avoid the convolution, especially if preliminary denoising is not needed or not desirable. Then it is possible to work directly with gradients of piecewise linear representation of 10 in the evaluation of g^.

Our co-volume scheme in this paper is designed for the specific mesh (see Fig. 11.18) given by the rectangular pixel structure of 2D image. For simplicity of implementation and for the reader's convenience, we will write the co-volume scheme in a "finite-difference notation." As is usual for 2D rectangular grids, we associate co-volume p and its corresponding center (DF node) with a couple (i, j), i will represent the vertical direction and j the horizontal direction. If ^ is a rectangular subdomain of the image domain where n1 and n2 are number of pixels in the vertical and horizontal directions, respectively, then i = 1,..., mi, j = 1,..., m2, mi < ni — 2, m2 < n2 — 2 and M = m1m2. Similarly, the unknown value un p is associated with un i j. For every co-volume p, the set Np consists of eight triangles (see Fig. 11.18). In every discrete time step n = 1,..., N, and for every i = 1,..., m1, j = 1,..., m2, we compute absolute value of gradient on these eight triangles denoted by Gk j, k = 1,..., 8. For that goal, using discrete values of u from the previous time step, we use the following expressions (we omit upper index n — 1 on u):

Gi,j

1,j2

1,j2

In the same way, but only in the beginning of the algorithm, we compute values Gaj, k = 1,..., 8, changing u by in the previous expressions, where is a smoothed image as explained in the paragraph above. Then for every i = 1,..., m1, j = 1,..., m2 we construct (north, west, south, and east)

coefficients n

to define diagonal coefficients

Ci,j = ni,j + wi,j + Si, j + j + mi, jh • If we define right-hand sides at the nth discrete time step by ri,j = jh , then for DF node corresponding to couple (i, j) we get the equation

Ci,ju i,j ni,ju i+1,j wi,ju i,j-1 ju i-1,j ei,ju i,j+1 = ri,j• (11.28)

Collecting these equations for all DF nodes and taking into account Dirichlet boundary conditions, we get the linear system to be solved.

We solve this system by the so-called SOR (successive over relaxation) iterative method, which is a modification of the basic Gauss-Seidel algorithm (see e.g. [46]). At the nth discrete time step we start the iterations by setting U3 = un-X, i = 1,--., mi, 3 = 1,..., m2. Then in every iteration l = 1,... and forevery i = 1,..., m1,3 = 1,..., m 2, we use the following two-step procedure:

Y = (Si,3Ui-1,3 + wi,3ui, 3-1 + ei,3ui, 3+1 + ni,3Ui+1,3 + ri,3 )'Ci,3

We define squared L2 norm of residuum at current iteration by

R — (ci,jui,j — ni,jUi+1,j — wi,jui,j-1 — Si,jUi-1,j — ei,jui,j+1 — ri,j) •

The iterative process is stopped if R® < TOL R0). Since the computing of residuum is time consuming itself, we check it, e.g., after every ten iterations. The relaxation parameter m is chosen by a user to improve convergence rate of the method; we have very good experience with m — 185 for this type of problems. Of course, the number of iterations depends on the chosen precision TOL, length of time step t , and a value of the regularization parameter e also plays a role. If one wants to weaken this dependence, more sophisticated approaches can be recommended (see e.g. [25,35,46] and paragraph below) but their implementation needs more programming effort. The semi-implicit co-volume method as presented above can be implemented in tens of lines.

We also outline shortly further approaches for solving the linear systems given in every discrete time step by (11.23). The system matrix has known (penta-diagonal) structure and moreover it is symmetric and diagonally dominant M-matrix. One could apply direct methods as Gaussian elimination, but this approach would lead to an immense storage requirements and computational effort. On the contrary, iterative methods can be applied in a very efficient way. In the previous paragraph we have already presented one of the most popular iterative methods, namely SOR. This method does not need additional storage, the matrix elements are used only to multiply the old solution values and convergence can be guaranteed for our special structure and properties of the system matrix. However, if the convergence is slow due to condition number of the system matrix (which increases with number of unknowns and for increasing t and decreasing e), faster iterative methods can be used. For example, the preconditioned conjugate gradient methods allow fast convergence, although they need more storage. If the storage requirements are reduced, then they can be very efficient and robust [25,35]. For details of implementation of the efficient preconditioned iterative solvers for co-volume level set method, we refer to [25], cf. also [51]. Also an alternative direct approach based on operating splitting schemes can be recommended [57,58].

In the next section, comparing CPU times, we will show that semi-implicit scheme is much more efficient and robust than explicit scheme for this type of problems. The explicit scheme combined with finite differences in space is usually based on formulations such as (11.9) [7-9,30,31,48-50] where all derivatives are expanded to get curvature and advection terms. Then, e.g., Eq. (11.2) for e = 1 is written in the form

0(1 + ux2)uX1 — 2uX1 ux uX1 x + (1 + 'uXl)uX X 0 0

where us means partial derivative of a function u with respect to a variable s and xi and X2 are spatial coordinates in the plane. In this form, it is not clear (reader may try) which terms to take from previous and which on the current time level, having in mind the unconditional stability of the method. Fully implicit time stepping would lead to a difficult nonlinear system solution, so the explicit approach is the one straightforwardly utilizable. In spite of that, the basic formulation (11.2) leads naturally to convenient semi-implicit time discretization.

Let us recall the usual criterion on numerical schemes for solving partial differential equations: numerical domain of dependence should contain physical domain of dependence. In diffusion processes, in spite of advection, a value of solution at any point is influenced by any other value of solution in a computational domain. This is naturally fulfilled by the semi-implicit scheme. We solve linear system of equations at every time step which, at every discrete point, takes into account contribution of all other discrete values in computational domain.