Image Processing

Fig. 1

Echocardiographic frame and its recognized Area

In the next part we focus the attention concerning the choice of the filter function for the edge-detection and the related applicability of the speed term formulation in the continuous model problem.

Fig. 2

Curve at isolevel at $$t=0$$

1.1 The Segmentation Problem

The aim of segmentation is to find a partition of an image into its costituent parts. For front evolution the

—

standard model is adopted [1–3](Fig. 2).

A curve in ${\mathbb {R}}^2$ can be represented as the zero-level line of a function in higher dimension. More precisely, let us suppose that there exists a function $u:{\mathbb {R}}^2\,\times \,{\mathbb {R}}^+\rightarrow {\mathbb {R}}$ solution of the initial value problem:

$\begin{aligned} \left\{ \begin{array}{ll} u_t (x,t)-g_I(x) \,\! |\nabla u(x,t) |=0 &{} \qquad x \in \varOmega ,\;t\ge 0\\ u(x,0)=u_0(x) &{} \qquad x \in \varOmega \end{array}\right. \end{aligned}$

(1)

which is the model of eikonal equation for front evolution for segmentation problem, typically involved in the detection for edges of objects contained in an image. We have to note that the evolving function

always remains a function as long as

is smooth. So we choose a different expression for the speed terms in order to limit the loss of image definition of Gaussian smoothing and by new approach with functional minimization to preserve this property. For ours purposes we intend enough smooth a function $C^2(\varOmega )$ , [4–12].

1.2 Edge-Detection

Given an image, in particular an echocardiographic image, we need to automatically detect the internal contour of object for clinical pourposes. We focus ours attention to a ventricular cavity. We consider a closed subset $\varOmega$ of ${\mathbb {R}}^{2},$ having $\partial \varOmega$ as the external edge of the rectangular image and being $I:\varOmega \rightarrow [0,1]$ the brightness intensity. The “edge-detector” (detector of contours) for segmentation of the image is a positive real coefficient, which is dependent of

’s gradient at every point of the curve [2, 13].

In particular, the model is represented by a filter function such that

$g:{\mathbb {R}}^+\longrightarrow {\mathbb {R}}^+$

where

$g(z)=\frac{1}{1+z}$

decreasing with

and

$\begin{aligned} \lim _{z \rightarrow \infty } g(z) = 0 \end{aligned}$

such that $0 \leqslant g$ and

1.3 Speed Choice for Image Processing

The speed term is dependent of the brightness intensity at every pixel and is directed along the outward normal, starting from an initial elliptic profile

centered in $\varOmega .$

Insofar we define the composite function

$\begin{aligned} g_I(|\nabla I(x)|)= g(|\nabla (G_\sigma * I(x)) |)= \dfrac{1}{1+ |\nabla (G_\sigma * I(x)) |}, \end{aligned}$

(2)

where $G_ \sigma$ *

is the convolution of the image

with a regularizing operator $G_ \sigma$ whose algorithmic implementation we are going to explicit later. The operator $G_{\sigma }$ is a filter that allows, by heat equation, to calculate the brightness-intensity gradient of the image in the presence of discontinuous data.

The discretization of the problem is performed by building a rectangular lattice, which is made fine-grained according to image’s pixel definition. The curve is parametrized by means of a Lipschitz function. Therefore, the differential problem to be numerically solved will look as follows:

Considering the built-in function

such that $E = \left\{ x \in {\mathbb {R}}^2: u_0(x)=0\right\}$

$\begin{aligned} \left\{ \begin{array}{ll} u_t - g_I (|\nabla I(x)|) | \nabla u | = 0 &{} \qquad x \in \varOmega \subset {\mathbb {R}}^2 \times [0,T]\\ u(x,0)=u_0(x) &{} \qquad x \in \varOmega \subset {\mathbb {R}}^2\\ \end{array} \right. \end{aligned}$

(3)

where T is the time horizon.

In the numerical tests the algorithm predicts a stop of the evolving curve at a threshold value

for the speed term, such that, if $g_{I}(x)\leqslant th\Rightarrow g_{I}(x)=0,$ then the threshold parameter allows curve evolution to stop in the presence of different gradient values for the regularized image.

2 New Approach to Regularize and Enhance Image

Analysis of images through variational methods finds application in a number of fields such as robotics, elaboration of satellite data, biomedical analysis and many other real-life cases. By segmentation is meant a search for constituent parts of an image, rather than an improving of its quality or characteristics. Our aim is to develop a criterium for enhancing movie frame in two possible ways by functional minimization. First, we adopt the M-S functional and its approximated form proposed by Ambrosio and Tortorelli [14] to regularize data frame for curve evolution method, instead of Gaussian regularization.

We call it “classic functional” since the function space the variational integral converges over contains functions that have the regularity needed for image gradient calculus. Second, we present a numerical scheme, where a time-dependent parameter is inserted, precisely in the second integral part of functional distinguished by a “time gradient”, for enhancing the internal moving parts of the object.

A detailed analytical treatment and a numerical scheme for minimization of the functional, which involves some delicate conjectures and refined mathematical steps, can be found in [15]. In the following section we recall in brief the essential formulation of the model problem used to regularize and enhance image. Other particular on numerical approximation and function space are explained in the next parts and in references. The reader can look up, for a complete review, the book by Morel and Solimini[16].

2.1 Ambrosio-Tortorelli Algorithm (’90)

This section refers to M-S algorithm for the approximation of $$F(u,K)$$

with a sequence $F_{\varepsilon }$ of regular functionals defined on a Sobolev space. We focus our attention on the Ambrosio-Tortorelli approximation, which is among the most used in image analysis. In this particular approach, the set $$K$$

(or $S_{u}$ ) is replaced by an auxiliary variable $$S$$

(a function) which approximates the characteristic function $(1-\chi _{K}).$

$F_{\varepsilon }\left( u, S \right) =\int _\varOmega \left( u-g \right) ^2 dx + \int _{\varOmega } S^2 \mid \nabla u\mid ^2 dx \quad +\int _{\varOmega } \left( \varepsilon \mid \nabla S\mid ^2 + \frac{1}{4 \varepsilon } (1-S)^2 \right) dx$

If $(u_{\varepsilon },S_{\varepsilon })$ minimizes the functional $F_\varepsilon ,$ then the following result holds [14]:

$\mathop {u_\varepsilon \longrightarrow u}\limits ^{\fancyscript{L} ^2} \;\; \text{ ed } \; S_\varepsilon \longrightarrow 1 - \chi _K \; \; \text{ per } \; \varepsilon \longrightarrow 0^+$

We gather from the problem of minimun related to the Ambrosio-Tortorelli functional, the Euler equations, than by an appropriate approximation will give us an algorithm for minimization.

When we apply Euler equation system, if $u_{0}=g$ is chosen, an enough rapid continuation method will give good results even though contours are not well defined.

2.2 Euler Equation of the Approximated Functional

Given

$\int _{\varOmega } \psi (u, \nabla u) dx,$

the associated Euler equation is

$div \left( \nabla _\xi \psi (u,\nabla u) \right) - \left( \frac{\partial \psi }{\partial u}\right) =0.$

Using Neumann boundary conditions, we get:

$F\left( u\right) =\mu \int _\varOmega \left( u-g \right) ^2 + \int _{\varOmega } (S^2+K_\varepsilon ) \mid \nabla u\mid ^2 dx \quad + \alpha \int _{\varOmega } \left( \varepsilon \mid \nabla S\mid ^2 + \frac{1}{4 \varepsilon } (1-S)^2 \right) dx.$

Euler equation system for $$u,$$

considering that the quadratic gradient $2(S^{2}+K_{\varepsilon })\nabla u^{2}$ is transformed into the Laplacian, looks as follows:

$\begin{aligned}\left\{ \begin{array}{cc} div (2(S^2+K_\varepsilon )\nabla u)=\mu (u-g) &{} \frac{\partial u}{\partial \mathbf n }=0 \; \; \text{ in } \;\partial \varOmega \\ \alpha \varepsilon \Delta S= S\mid \nabla u\mid ^2-\frac{\alpha }{4 \varepsilon }(1-S) &{} \frac{\partial S}{\partial \mathbf n }=0 \; \; \text{ in } \;\partial \varOmega \end{array}\right. \end{aligned}$

where the non linear terms $S^{2}$ and $\nabla u$ cause the system to be non-linear and elliptic, with such a structure that, when $$S$$

is known, the first equation gets linear, while, if it is $$u$$

to be known, it is the second equation to be linear. This suggests the adoption of a two-stage iterative scheme. Where we fix the discrete step of convergence that is enough to give back a regularized image, enhanced in its content edge, as we present in the next section.

3 Approximation of the Models

We are going to describe, in this section, the numerical approximation of the models introduced above. In details, our aim is to present the semilagrangian scheme for curve evolution by numerical solution of eikonal equation. The standard mollifier is constructed by a Gaussian operator as a numerical scheme for heat eqution in two dimensions. The last section reports the procedure to build the approximation of the Mumford-Shah functional by Ambrosio-Tortorelli scheme.

3.1 Semi-Lagrangian Scheme for Eikonal Equation

Semi-Lagrangian (SL) schemes try to mimic the continuous behavior by constructing the solution at each grid point by back integration along the characteristic trajectory passing through the point and reconstructing the value at the foot of the trajectory by interpolation.

The numeric dependence of the domain contains its continuous dependence without any additional condition on $\Delta t$ and $\Delta x$ (space lattice is usually defined for finite differences as the pixel resolution of images). This allow a larger time steps than other schemes where the CFL condition has to be imposed for stability guaranty.

The numerical method for eq. (3) following the semilagrangian scheme by [17]. In our implementation we use a threshold value that has been chosen and discussed in a Master thesis [18].

We construct the protocol step of curve evolution by an oriented adaptation to medical application of HJPACK parallel OpenMp Fortran programming language.

3.2 Ambrosio-Tortorelli Approximation of the M-S Functional

The numerical scheme is made by dividing in two coupled parts with $$u_0=g $$

and

At every step we calculate $u_{1}$ for $S_{0}=1$ solving a linear elliptic equation and, this way, we find $S_{1}$ from the second equation; this process is repeated for a fixed number of iterations.

3.2.1 System Discretization

We use a numerical scheme based on explicit finite differences, over the rectangle $\varOmega$ with step $$h$$

; this way we obtain $$(x,y)=(ih,jh)$$

for $0\leqslant i,j\leqslant N.$ Reducing $\varOmega$ to a square of side 1 and taking $h=\frac{1}{N},$ discrete coordinates will become: $u(ih,jh)\cong u_{i,j}$ and $S(ih,jh)\cong S_{i,j}.$

We use an approximated scheme that is enough to enhancing little areas, characteristic of the echographic image. Some scheme for the Ambrosio-Tortorelli segmentation problem [19].

In order to determine the minimum of the functional we adopt the schema given, for a finite element, in [20] to a finite difference meshgrid through the following iterative scheme:

given a maximum number of iteration $$Nit$$

and a tolerance ‘ $\varepsilon$ , then we construct:

for find , by solving:

$\begin{aligned}\left\{ \begin{array}{ll} div (2(S^2_{n-1}+K_\varepsilon )\nabla u_n)=\mu (u_n-g) &{}\; \; \text{ in } \; \varOmega \\ \frac{\partial u_n}{\partial \mathbf n }=0 &{}\; \; \text{ in } \;\partial \varOmega \\ \end{array}\right. \end{aligned}$

and

by solving:

$\begin{aligned}\left\{ \begin{array}{ll} \alpha \varepsilon \Delta S_n= S_n\mid \nabla u_n\mid ^2-\frac{\alpha }{4 \varepsilon }(1-S_n) &{} \; \; \text{ in } \; \varOmega \\ \frac{\partial S_n}{\partial \mathbf n }=0 &{}\; \; \text{ in } \;\partial \varOmega \\ \end{array}\right. \end{aligned}$

stop for .

From minimization theorem 3.1 Proposition 2.1 [20], it respectively follows that:

is a piecewise $$C^2$$

submanifolds of $\mathbb {R}^2$ .

For any

$$n>1$$” src=”/wp-content/uploads/2016/03/A320009_1_En_16_Chapter_IEq68.gif”></SPAN> there exists <SPAN id=IEq69 class=InlineEquation><IMG alt=$$u_n$$ src=

solution of the respective system which satisfy the bounds:

$\Vert u_n\Vert _{L^\infty } \le \Vert g\Vert _{L^\infty }.$

Then we discretize the equations by finite differences.

3.2.2 The Discreet Divergence

In a numerical scheme a function $u \in {\mathbb {R}}^2$ can be approximated by finite difference, its first order variation on the $$x$$

direction is:

$\;\;\frac{\partial u}{\partial x}\cong \frac{u_{i+1,j}-u_{i,j}}{h}$

and for

direction is:

$\;\;\frac{\partial u}{\partial y}\cong \frac{u_{i,j+1}-u_{i,j}}{h}.$

The divergence of a function $$Z(x,y)$$

by second order of centered difference is given by

$div(Z(x,y) \nabla u(x,y))\cong$

$\cong Z_{i+\frac{1}{2},j}(u_{i+1,j}-u_{i,j})-Z_{i-\frac{1}{2},j}(u_{i,j}-u_{i-1,j})$

$+\,Z_{i,j+\frac{1}{2}}(u_{i,j+1}-u_{i,j})-Z_{i,j+\frac{1}{2}}(u_{i,j-1}-u_{i,j}).$

where

$Z_{i+\frac{1}{2},j}=\frac{1}{2} (Z_{i+1,j}+Z_{i,j}) \;\; ;\;\; Z_{i-\frac{1}{2},j}=\frac{1}{2} (Z_{i,j}+Z_{i,-1j})$

$Z_{i,j+\frac{1}{2}}=\frac{1}{2} (Z_{i,j+1}+Z_{i,j}) \;\; ;\;\; Z_{i,j-\frac{1}{2}}=\frac{1}{2} (Z_{i,j}+Z_{i,j-1})$

By applying to the system of Euler equation for the $n^{th}$ approximated item of the sequence of the Ambrosio-Tortorelli functional we obtain: the term $Z(x,y)=(S^2(x,y) + K_\varepsilon )$ become by fixing every direction of the space, with $K_\varepsilon$ neglected as in [19], we get: for x

$\frac{\partial }{\partial x}\left( \left( S^2+K_\varepsilon \right) \frac{\partial u}{\partial x}\right) \cong$

$\frac{1}{2h^2}\left[ (S^2 _{i+1,j} +S^2 _{i,j})(u_{i+1,j}-u_{i,j})+ (S^2 _{i,j} +S^2 _{i-1,j})(u_{i-1,j}-u_{i,j})\right] ,$

for y

$\frac{\partial }{\partial y}\left( \left( S^2+K_\varepsilon \right) \frac{\partial u}{\partial y}\right) \cong$

$\frac{1}{2h^2}\left[ (S^2 _{i,j+1} +S^2 _{i,j})(u_{i,j+1}-u_{i,j})+ (S^2 _{i,j} +S^2 _{i,j-1})(u_{i,j-1}-u_{i,j})\right] .$

In Image Processing the lattice has a spatial density which is equivalent to image resolution, so we use here finite differences with node resolution equal to the spatial grid-step, so the discretized equation will look as follows

$\begin{aligned}\begin{array}{c} (S^2 _{i+1,j} +S^2 _{i,j})(u_{i+1,j}-u_{i,j}) + (S^2 _{i,j} +S^2 _{i-1,j})(u_{i-1,j}-u_{i,j})\\ +\,(S^2 _{i,j+1} +S^2 _{i,j})(u_{i,j+1}-u_{i,j}) + (S^2 _{i,j} +S^2 _{i,j-1})(u_{i,j-1}-u_{i,j})\\ =\,\mu h^2 (u_{i,j}-g_{i,j}). \end{array}\end{aligned}$