A High-Order Scheme for Image Segmentation via a modified Level-Set method

TL;DR

This paper introduces a novel high-order level-set scheme for image segmentation that enhances stability and accuracy by modifying the velocity function and employing an adaptive filtered scheme, validated through synthetic and real image tests.

Contribution

It presents the first 2D adaptive filtered scheme for level-set image segmentation, combining high-order accuracy with stability and efficiency.

Findings

Improved segmentation accuracy on synthetic and real images.

Enhanced stability of the level-set evolution.

Efficient implementation of the high-order scheme.

Abstract

In this paper we propose a high-order accurate scheme for image segmentation based on the level-set method. In this approach, the curve evolution is described as the 0-level set of a representation function but we modify the velocity that drives the curve to the boundary of the object in order to obtain a new velocity with additional properties that are extremely useful to develop a more stable high-order approximation with a small additional cost. The approximation scheme proposed here is the first 2D version of an adaptive "filtered" scheme recently introduced and analyzed by the authors in 1D. This approach is interesting since the implementation of the filtered scheme is rather efficient and easy. The scheme combines two building blocks (a monotone scheme and a high-order scheme) via a filter function and smoothness indicators that allow to detect the regularity of the approximate solution adapting the scheme in an automatic way. Some numerical tests on synthetic and real images confirm the accuracy of the proposed method and the advantages given by the new velocity.

Figures

Tables

Table 1. Table 1: Summary of the initial conditions considered in the numerical tests.

	$c$	$\tilde{c}$
Expansion case	$u_0$ as Paraboloid (59)	Datum 1: $u_0$ as Paraboloid (59)
(Case a)		Datum 2: $u_0$ as Distance from ${\Gamma }_0^1$ (61)
Shrinking case	$u_0$ as Truncated pyramid
(Case b)	(or tent) (60)	$u_0$ as Distance from ${\Gamma }_0^2$ (62)

Table 2. Table 2: Test 1a. Errors and number of iterations using L ∞ superscript 𝐿 L^{\infty} norm, Datum 1 1 1 and the parameters μ = 2 𝜇 2 \mu=2 , K r ​ e ​ g = 0 subscript 𝐾 𝑟 𝑒 𝑔 0 K_{reg}=0 , t ​ o ​ l = 0.0005 𝑡 𝑜 𝑙 0.0005 tol=0.0005 , varying the number of nodes. Best results are in bold.

$c$	Monotone			AF-LW			WENO
$\#Nodes$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$
$102$	$\mathrm{𝟖𝟒}$	$0.1025$	$0.2321$	$213$	$X$	$X$	$217$	$X$	$X$
$202$	$\mathrm{𝟏𝟓𝟐}$	$0.0526$	$0.1172$	$394$	$X$	$X$	$399$	$X$	$X$
$402$	$\mathrm{𝟐𝟖𝟖}$	$0.0265$	$0.0593$	$745$	$X$	$X$	$669$	$X$	$X$

Table 3. Table 3: Test 1a. Errors and number of iterations using L ∞ superscript 𝐿 L^{\infty} norm, Datum 1 1 1 and the parameters μ = 2 𝜇 2 \mu=2 , K r ​ e ​ g = 0 subscript 𝐾 𝑟 𝑒 𝑔 0 K_{reg}=0 , t ​ o ​ l = 0.0005 𝑡 𝑜 𝑙 0.0005 tol=0.0005 , varying the number of nodes. Best results are in bold.

$\tilde{c}$	Monotone			AF-LW			WENO
$\#Nodes$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$
$102$	$50$	$0.0748$	$0.1694$	$48$	$0.0693$	$0.1568$	$\mathrm{𝟒𝟕}$	$0.0886$	$0.2008$
$202$	$100$	$0.0427$	$0.0950$	$\mathrm{𝟗𝟔}$	$0.0363$	$0.0808$	$\mathrm{𝟗𝟔}$	$0.0469$	$0.1045$
$402$	$199$	$0.0208$	$0.0466$	$\mathrm{𝟏𝟗𝟔}$	$0.0203$	$0.0454$	$\mathrm{𝟏𝟗𝟔}$	$0.0240$	$0.0537$

Table 4. Table 4: Test 1a. CPU times in seconds related to the Tab. 3 .

$\#Nodes$	Monotone	AF-LW	WENO
102	$0.19$	$0.64$	$4.75$
202	$1.35$	$5.61$	$68.38$
402	$9.77$	$40.20$	$1028.51$

Table 5. Table 5: Test 1a. Errors and number of iterations varying # ​ N ​ o ​ d ​ e ​ s # 𝑁 𝑜 𝑑 𝑒 𝑠 \#Nodes ( L ∞ superscript 𝐿 L^{\infty} norm). Best results are in bold.

$\tilde{c}$	Monotone			Datum $1$			Datum $2$
$\#Nodes$	$tol$	$\mu$	$K_{reg}$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$
$102$	$0.001$	$2$	$0$	$\mathrm{𝟓𝟎}$	$0.0748$	$0.1694$	$\mathrm{𝟓𝟎}$	$0.0776$	$0.1757$
$202$	$0.0005$	$2$	$0$	$\mathrm{𝟏𝟎𝟎}$	$0.0427$	$0.0950$	$\mathrm{𝟏𝟎𝟎}$	$0.0420$	$0.0935$
$402$	$0.00025$	$2$	$0$	$\mathrm{𝟏𝟗𝟗}$	$0.0208$	$0.0466$	$200$	$0.0199$	$0.0446$

Table 6. Table 6: Test 1a. Errors and number of iterations varying # ​ N ​ o ​ d ​ e ​ s # 𝑁 𝑜 𝑑 𝑒 𝑠 \#Nodes ( L ∞ superscript 𝐿 L^{\infty} norm). Best results are in bold.

$\tilde{c}$	AF-LW			Datum $1$			Datum $2$
$\#Nodes$	$tol$	$\mu$	$K_{reg}$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$
$102$	$0.001$	$2$	$0$	$48$	$0.0693$	$0.1568$	$\mathrm{𝟒𝟕}$	$0.0658$	$0.1490$
$202$	$0.0005$	$2$	$0$	$\mathrm{𝟗𝟔}$	$0.0363$	$0.0808$	$\mathrm{𝟗𝟔}$	$0.0363$	$0.0808$
$402$	$0.00025$	$2$	$0$	$196$	$0.0203$	$0.0454$	$\mathrm{𝟏𝟗𝟓}$	$0.0201$	$0.0450$

Table 7. Table 7: Test 1a. Errors and number of iterations varying # ​ N ​ o ​ d ​ e ​ s # 𝑁 𝑜 𝑑 𝑒 𝑠 \#Nodes ( L ∞ superscript 𝐿 L^{\infty} norm). Best results are in bold.

$\tilde{c}$	WENO			Datum $1$			Datum $2$
$\#Nodes$	$tol$	$\mu$	$K_{reg}$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$
$102$	$0.001$	$2$	$0$	$\mathrm{𝟒𝟕}$	$0.0886$	$0.2008$	$\mathrm{𝟒𝟕}$	$0.0886$	$0.2008$
$202$	$0.0005$	$2$	$0$	$\mathrm{𝟗𝟔}$	$0.0469$	$0.1045$	$\mathrm{𝟗𝟔}$	$0.0469$	$0.1045$
$402$	$0.00025$	$2$	$0$	$\mathrm{𝟏𝟗𝟔}$	$0.0240$	$0.0537$	$\mathrm{𝟏𝟗𝟔}$	$0.0242$	$0.0541$

Table 8. Table 8: Test 2a. CPU times in seconds related to the brain tests (Figs. 5 - 6 ).

Figure	Mon.(left)	Mon.(right)	AF-LW	WENO
5	$12.95$	$16.54$	$71.18$	$1424.18$
6	$9.63$	$20.31$	$91.41$	$1318.59$

Table 9. Table 9: Test 2b. Errors and number of iterations using L 1 superscript 𝐿 1 L^{1} norm, c ~ ~ 𝑐 \widetilde{c} and the parameters t ​ o ​ l = 0.00005 𝑡 𝑜 𝑙 0.00005 tol=0.00005 , μ = 5 𝜇 5 \mu=5 , K r ​ e ​ g = 3 subscript 𝐾 𝑟 𝑒 𝑔 3 K_{reg}=3 , varying the image size. Best results are in bold.

	Monotone			AF-LW			WENO
Image size	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$
$170\times 170$	$\mathrm{𝟖𝟑}$	$0.0465$	$0.2436$	$87$	$0.0439$	$0.2300$	$87$	$0.0448$	$0.2348$
$340\times 340$	$\mathrm{𝟐𝟐𝟖}$	$0.0118$	$0.2476$	$262$	$0.0078$	$0.1628$	$264$	$0.0079$	$0.1648$

Table 10. Table 10: Test 2b. CPU times in seconds related to the brain tests (Fig. 7 ).

Image size	Monotone	AF-LW	WENO
$170\times 170$	$0.80$	$3.75$	$44.10$
$340\times 340$	$8.39$	$46.07$	$933.90$

Table 11. Table 11: Test 3a. Errors and number of iterations using L 1 superscript 𝐿 1 L^{1} norm, Datum 1 1 1 , velocity c ~ ~ 𝑐 \widetilde{c} , and the parameters μ = 2 𝜇 2 \mu=2 , K r ​ e ​ g = 5 subscript 𝐾 𝑟 𝑒 𝑔 5 K_{reg}=5 , varying the tolerance. Best results are in bold.

$\tilde{c}$	Monotone			AF-LW			WENO
tol	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$
$0.00005$	$433$	$0.0794$	$0.6424$	$\mathrm{𝟑𝟗𝟔}$	$0.0672$	$0.6424$	$432$	$0.0591$	$0.4788$
$0.000025$	$\mathrm{𝟓𝟎𝟖}$	$0.0599$	$0.4848$	$544$	$0.0438$	$0.3544$	$523$	$0.0462$	$0.3744$

Table 12. Table 12: Test 3a. Errors and number of iterations using L 1 superscript 𝐿 1 L^{1} norm, Datum 2 2 2 , velocity c ~ ~ 𝑐 \widetilde{c} , and the parameters μ = 2 𝜇 2 \mu=2 , K r ​ e ​ g = 5 subscript 𝐾 𝑟 𝑒 𝑔 5 K_{reg}=5 , varying the tolerance. Best results are in bold.

$\tilde{c}$	Monotone			AF-LW			WENO
tol	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$
$0.00008$	$\mathrm{𝟑𝟖𝟐}$	$0.1031$	$0.8344$	$451$	$0.0640$	$0.5180$	$411$	$0.0684$	$0.5536$
$0.00004$	$\mathrm{𝟒𝟓𝟑}$	$0.0737$	$0.5964$	$542$	$0.0501$	$0.4052$	$472$	$0.0565$	$0.4572$

Table 13. Table 13: Test 3a. CPU times in seconds related to the chess horse tests (Tabs. 11 - 12 ).

Datum	tol	Monotone	AF-LW	WENO
$1$	$0.00005$	$7.26$	$26.22$	$378.07$
$1$	$0.000025$	$8.78$	$42.23$	$462.55$
$2$	$0.00008$	$6.51$	$37.15$	$377.47$
$2$	$0.00004$	$8.02$	$44.16$	$443.66$

Table 14. Table 14: Test 3b. Errors and number of iterations using L 1 superscript 𝐿 1 L^{1} norm, velocity c ~ ~ 𝑐 \widetilde{c} and the parameters t ​ o ​ l = 0.00005 𝑡 𝑜 𝑙 0.00005 tol=0.00005 , μ = 2 𝜇 2 \mu=2 , K r ​ e ​ g = 3 subscript 𝐾 𝑟 𝑒 𝑔 3 K_{reg}=3 . Best results are in bold.

Monotone			AF-LW			WENO
$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$
$\mathrm{𝟏𝟖𝟕}$	$0.0303$	$0.2452$	$197$	$0.0273$	$0.2212$	$196$	$0.0262$	$0.2120$

Table 15. Table 15: Test 3b. CPU times in seconds related to the horse test (Tab. 14 ).

Monotone	AF-LW	WENO
$3.14$	$14.10$	$195.13$

Table 16. Table 16: Test 4b. Errors and number of iterations with velocity c ~ ~ 𝑐 \widetilde{c} , using L 1 superscript 𝐿 1 L^{1} norm in the stopping criterion, and the parameters t ​ o ​ l = 0.0001 𝑡 𝑜 𝑙 0.0001 tol=0.0001 , μ = 2 𝜇 2 \mu=2 , and K r ​ e ​ g = 2 subscript 𝐾 𝑟 𝑒 𝑔 2 K_{reg}=2 . Best results are in bold.

Monotone			AF-LW			WENO
$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$	$N_i$	$P$-$Er{r}_{rel}$	$P$-$Er{r}_1$
$312$	$0.0097$	$0.0316$	$\mathrm{𝟐𝟗𝟗}$	$0.0074$	$0.0240$	$306$	$0.0064$	$0.0208$

Table 17. Table 17: Test 4b. CPU times in seconds related to the grains test (Tab. 16 ).

Monotone	AF-LW	WENO
$8.71$	$41.26$	$676.91$

Table 18. Table 18: Test 6a. CPU times in seconds related to the Pneumonia tests (Figs. 18 - 19 ).

Figure	Monotone	AF-LW	WENO
18	$1.80$	$8.05$	$83.83$
19	$5.01$	$28.91$	$202.373$

Equations

$${\Gamma }_t:=\{(x,y):v(t,x,y)=0\}.$$

$$\left\{\begin{array}[]{ll}v_{t}+c(t,x,y)|\nabla v|=0,&(t,x,y)\in\mathbb{(}0,T)\times\mathbb{R}^{2},\\ v(0,x,y)=v_{0}(x,y),&(x,y)\in\mathbb{R}^{2},\end{array}\right.$$

$$\left\{\begin{array}[]{ll}v_{0}(x,y)<0,&(x,y)\in\Omega_{0},\\ v_{0}(x,y)=0,&(x,y)\in\Gamma_{0},\\ v_{0}(x,y)>0,&(x,y)\in\mathbb{R}^{2}\setminus\overline{\Omega}_{0},\end{array}\right.$$

$$c_1(x,y)=\frac{1}{\left(1+|\nabla (G\ast I(x,y)){|}^{\mu }\right)},\mu \ge 1,$$

$$c_2(x,y)=1-\frac{|\nabla (G\ast I(x,y))|-m}{M-m},$$

$$v_{0}(x,y)=\left\{\begin{array}[]{ll}-dist\left\{(x,y),\Gamma_{0}\right\}\hbox{ in }\Omega_{0}\\ dist\left\{(x,y),\Gamma_{0}\right\}\hbox{ outside }\Omega_{0}\end{array}\right.$$

$$(x_0,y_0)=(x_c,y_c)-v(t,x_c,y_c)\frac{\nabla v(t,x_c,y_c)}{|\nabla v(t,x_c,y_c)|}.$$

$$\tilde{c}(x,y,v,v_x,v_y)=c\left(x-v\frac{v_x}{|\nabla v|},y-v\frac{v_y}{|\nabla v|}\right),$$

$$v_0(x,y)={\overline{v}}_0\left(\sqrt{x^2+y^2}\right).$$

$$d(C):={\overline{v}}_0^{-1}(C)-r_0,\text{with }{\overline{v}}_0^{-1}(C)\ge 0,$$

$$\tilde{c}(x,y,v,v_x,v_y)=c\left(x-d(v)\frac{v_x}{|\nabla v|},y-d(v)\frac{v_y}{|\nabla v|}\right).$$

$$H(z,v,p)=\tilde{c}(z,v,p)|p|.$$

$$\left\{\begin{array}[]{l}\dot{z}(s)=\nabla_{p}H\\ \dot{v}(s)=\nabla_{p}H\cdot p-H\\ \dot{p}(s)=-\nabla H-H_{v}p,\end{array}\right.$$

$$\frac{\partial H}{\partial p_1}=$$

$$=$$

$$\nabla_{p}H=\left(\begin{array}[]{l}\frac{d(v)p_{2}}{|p|^{2}}\left(p_{1}\frac{\partial c}{\partial\zeta}-p_{2}\frac{\partial c}{\partial\xi}\right)+\widetilde{c}\frac{p_{1}}{|p|}\\ \frac{d(v)p_{1}}{|p|^{2}}\left(p_{2}\frac{\partial c}{\partial\xi}-p_{1}\frac{\partial c}{\partial\zeta}\right)+\widetilde{c}\frac{p_{2}}{|p|}\\ \end{array}\right)\hbox{ that implies }\nabla_{p}H\cdot p=\widetilde{c}(z,v,p)|p|.$$

$$\left\{\begin{array}[]{l}\dot{z}(s)=\nabla_{p}H\\ \dot{v}(s)=\widetilde{c}(z,v,p)|p|-\widetilde{c}(z,v,p)|p|=0\\ \dot{p}(s)=-\nabla\widetilde{c}(z,v,p)|p|+d^{\prime}(v)\nabla\widetilde{c}(z,v,p)|p|^{2}=\nabla\widetilde{c}(z,v,p)|p|(d^{\prime}(v)|p|-1),\end{array}\right.$$

$$d^{\prime }(v)=|p{|}^{-1},$$

$$\left\{\begin{array}[]{l}\dot{z}(s)=\nabla_{p}H,\\ \dot{v}(s)=0,\\ \dot{p}(s)=0,\end{array}\right.$$

$$d^{\prime }(v)=|\nabla v_0{|}^{-1},$$

$$d^{\prime }(v)=\frac{d}{dv}\left({\overline{v}}_0^{-1}(v)\right)=\frac{1}{{\overline{v}}_0^{\prime }(w)},$$

$$|\nabla v_0(x,y)|$$

$$=\left|\left(\frac{{\overline{v}}_0^{\prime }\left(\sqrt{x^2+y^2}\right)x}{\sqrt{x^2+y^2}},\frac{{\overline{v}}_0^{\prime }\left(\sqrt{x^2+y^2}\right)y}{\sqrt{x^2+y^2}}\right)\right|$$

$$=\frac{{\overline{v}}_0^{\prime }\left(\sqrt{x^2+y^2}\right)}{\sqrt{x^2+y^2}}\sqrt{x^2+y^2}={\overline{v}}_0^{\prime }\left(\sqrt{x^2+y^2}\right),$$

$$z^0(0)=z^{\ell }(0)-d(v_0(z^{\ell }))\frac{\nabla v_0(z^{\ell })}{|\nabla v_0(z^{\ell })|},$$

$$\dot{\underline{z}}(s)=$$

$$=$$

$$\dot{\underline{z}}(s)=$$

$$=$$

$$\frac{d}{ds}v(s,\underline{z}(s))=v_s+v_x\dot{\underline{z}}(s)=v_s+v_x{\dot{z}}^{\ell }(s)=\frac{d}{ds}v(s,z^{\ell }(s))=0,$$

$$v_t+c(\underline{z})|\nabla v|=0,(t,x,y)\in (t_n,t_{n+1})\times {\mathbb{R}}^2,$$

$$|c(\underline{z})-c({\underline{z}}^{\prime })|$$

Full text

**Keywords: ** Image segmentation, level-set method, Hamilton-Jacobi equations, filtered scheme, smoothness indicators.

AMS subject classifications:

68U10, 35F21, 35Q68, 65M06, 65M25

1 Introduction

The Level-Set (LS) method has been introduced by Osher and Sethian in the 1980s [35, 31] and then used to deal with several applications, e.g. fronts propagation, computer vision, computational fluids dynamics (see the monographs by Sethian [36] and by Osher and Fedkiw [30] for several interesting examples). This method is nowadays very popular for its simplicity and for its capability to deal with topological changes. In fact, the main advantage of the LS method is the possibility to easily describe time-varying objects, follow shapes that change topology, for example when a shape splits in two, develops holes, or the reverse of these operations. For the image segmentation problem the application of the LS method is based on the evolution of a curve according to a normal velocity based on the gray-levels of the image, typically the curve is described as the 0-level set of a representation function (or level-set function).

In $\mathbb{R}^{2}$ the LS method corresponds to define an initial closed curve $\Gamma_{0}$ using an auxiliary function $v_{0}$ which has to change sign on $\Gamma_{0}$. The evolution of that curve at time $t$ is denoted by $\Gamma_{t}$ and is represented by the 0-level set of a function $v$, i.e.

[TABLE]

This function $v$ is the unique viscosity solution of the following evolutive nonlinear equation of Hamilton-Jacobi type

[TABLE]

where $\nabla v:=(v_{x},v_{y})$ denotes the spatial gradient of $v$. Usually, the velocity $c(t,x,y)$ does not change sign during the evolution and the orientation depends on the type of evolution (outward for an expansion and inward for a shrinking). Typically $v_{0}$ must be a proper representation of the initial front $\Gamma_{0}$, satisfying

[TABLE]

where $\Omega_{0}$ is the region delimited by $\Gamma_{0}$ (or the reverse inequalities).

The LS method can handle velocities also depending on physical quantities in order to describe several phenomena. Typical examples are:

1. a)

$c(t,x,y)$, isotropic growth with time varying velocity 2. b)

$c(t,x,y,\eta)$, anisotropic growth, dependent on normal direction 3. c)

$c(t,x,y,k)$, Mean Curvature Motion, with $k(t,x,y)$ mean curvature to the front at time $t$.

The literature on the LS method is huge as well as the range of applications where the method has been successfully applied. We refer the interested reader to [13, 28] for the isotropic case, to [25] for the anisotropic growth and to [14, 10] for the curvature case (see also the monographs [36, 30] and the references therein). Depending on the choice of the velocity, the front evolution will be described by first or second order partial differential equations, e.g. case c) gives rise to a degenerate second order equation. Here we limit ourself to first order problems (2) and velocities which depend only on time and position and we consider a velocity sign constant in time since this is the more relevant case for image segmentation (a change in sign is relevant for dislocations in material science and has been studied in [9, 12]).

To set our contribution into perspective, let us mention that the segmentation problem has been solved by various techniques which mainly rely on two different approaches: variational methods and active contour methods. For the first approach the interested reader can look at [8, 7] and the references therein. For the link between the two classes of methods, see [14]. To have an idea of other segmentation methods, we refer the reader to the surveys [40, 11] and to the papers [16, 38] for the so called “balloon model” first introduced in [15]. This method is based on the introduction of a potential giving a driving force to the segmentation process and typically leads to second order partial differential equations that are at present outside the reach of the filtered scheme proposed in this paper. For high-order Runge-Kutta methods in the framework of image segmentation we refer to [34]. As already said, we will apply the level set method based on (2) looking for an accurate numerical method. High-order methods have been proposed for (2) and most of them are based on non oscillatory local interpolation techniques that allow to avoid spurious oscillations around discontinuities of the solution and/or jumps in the gradient. These techniques were originally developed for conservation laws (see the seminal paper [21] and the references therein), the research activity on essentially non oscillatory (ENO) methods has been rather effective and a number of improvements have been proposed in e.g. [27, 24, 3] so that now ENO and Weighted ENO (WENO) techniques are rather popular in many applications (see [39] for a recent survey). We should also mention that later on these techniques were successfully applied to the numerical solution of Hamilton-Jacobi equations [32] opening the way to other applications. The reader interested in image processing will find an application to image compression in [1, 2] and an application to image segmentation in [37]. However, a general convergence theorem for ENO/WENO schemes is still missing and their application is a rather delicate issue. These limitations have motivated further investigations and a new class of high-order methods for evolutive Hamilton-Jacobi equations have been proposed looking for a different approach based on *“filtered schemes” * introduced in this framework by Lions and Souganidis [26] (it is important to note that the name does not refer to filtering a noise as is common in the imaging community but rather to the presence of a filter function as we will explain later). The class of filtered schemes is based on a simple coupling of a monotone scheme with a high-order scheme. Monotone schemes are convergent to the weak (viscosity) solution but they are known to be at most first order accurate, whereas high-order schemes gives a higher accuracy but in general are not stable. The crucial point is the coupling between the two schemes which is obtained via a filter function that selects which scheme has to be applied at a node of the grid in order to guarantee (under appropriate assumptions) a global convergence. The construction of these schemes is rather simple as explained by Oberman and Salvador [29] because one can couple various numerical methods and leave the filter function deciding the switch between the two schemes. A general convergence result has been proved by Bokanowski, Falcone and Sahu in [5] and recently improved by Falcone, Paolucci and Tozza [18] with an adaptive and automatic choice of the parameter governing the switch. Note that the adaptation of the parameter depends on some regularity indicators in every cell, these indicators are computed at each iteration and this guarantees convergence. Some contributions to extend filtered schemes to second order problems can be found in [20] for the stationary case and in [6] for the evolutive case and financial applications.

Our contribution

The contribution of this paper is twofold: from a theoretical point of view, we propose a modified velocity for the segmentation problem in the level-set approach. This new velocity is important (and necessary) for the numerical approximation of high-order schemes, like the filtered scheme considered here and its properties will allow to avoid the re-initialization procedure used e.g. in [28]. We also improve the accuracy of the method by applying an adaptive high-order filtered scheme to the segmentation problem. This requires a 2D extension of the scheme proposed in 1D in [19], for which convergence has been proved under rather general assumptions in [18], thanks to the definition of new full 2D smoothness coefficients (see Sect. 3). We will show that this choice is competitive with respect to other high-order schemes, as the WENO scheme, in terms of accuracy and computational cost, analyzing several experiments. This contribution is more efficient in terms of CPU time and is interesting from a theoretical point of view since a precise convergence result for WENO scheme is still missing.

Paper organization

The paper is organized as follows: In Sect. 2, we give the idea behind the new definition of the velocity function, with details on its construction and explaining why it is so important to introduce it when applying high-order approximation schemes. In Sect. 3, we briefly present the Adaptive Filtered (AF) scheme and recall some of its basic elements. In Sect. 4, we give some information on the implementation and we give a sketch of the Algorithm 1 for the solution of the segmentation problem via the new AF scheme. Finally, the performances of this new method for the segmentation problem are illustrated in Sect. 5 where we compare it with a classical monotone scheme (first order accurate) and another high-order scheme (the WENO scheme) on a series of virtual and real images presenting a detailed error analysis of the numerical experiments.

2 Image segmentation via a modified LS method

The boundaries of one (or more) object(s) inside a given image are characterized by an abrupt change of the intensity values $I(x,y)$ of the image, so that the magnitude of $|\nabla I(x,y)|$ can be used as an indication of the edges. Let us assume that $\nabla I$ exists at least almost everywhere, the definition of the velocity $c(x,y)$ in (2) plays a crucial role and must be defined in a proper way in order to guarantee that the curve evolution is close to [math] when the front is close to an edge, this will stop the evolution. The normal velocity will be positive or negative depending on the case, expanding or shrinking, respectively, and if the velocity is just given in term of $I$ it will not change sign during the evolution, as it could happen when the velocity depends on other geometrical properties of the curve (e.g. the curvature). We will focus on velocities which ignore the curvature since the numerical approach we present here is well adapted to first order evolutive Hamilton-Jacobi equations (the extension to second order problems goes beyond the scopes of this paper). Note that, in order to reduce the noise, several methods have been proposed in the literature. A very simple one is to use the convolution with a Gaussian kernel. This can be obtained by evolving the original function $I$ of the gray levels according to the heat equation for a short time interval (see Sect. 4 for more details). Due to the regularizing effect of the heat equation this also guarantees that $\nabla I$ exists.

Several definitions of the velocity function $c(x,y)$ have been proposed in literature. A typical example is

[TABLE]

where $\mu$ is used to give more weight to the changes in the gradient, if necessary. In [13] the authors proposed that velocity with $\mu=2$, and in [28] with $\mu=1$. According to this definition, the velocity takes values in $[0,1]$ and has values that are close to zero at points where the image gradient is high and equal to 1 where $I$ is constant.

Another possible choice has been proposed in [28] and has the form

[TABLE]

where $m$ and $M$ are respectively the minimum and the maximum values of $|\nabla(G\ast I(x,y))|$. This velocity has similar properties with respect to (4) and takes values in $[0,1]$ but is close to [math] if the magnitude of the image gradient is close to its maximal value, and equal to $1$ otherwise.

It is clear that both definitions have the desired properties, but with slightly different features. More precisely, in the first case the velocity depends more heavily on the changes in the magnitude of the gradient, allowing for an easier detection of the edges but possibly producing false edges inside the object (e.g. when specular effects are present in the image). The velocity (5) is smoother inside the objects, being less dependent on the relative changes in the gradient, but might present some problems in the detection of all the edges if at least one of those is “more marked”.

2.1 Extension of the velocity function

The edge-stopping function which is defined choosing one of the above mentioned velocities, has a physical meaning only on the front $\Gamma_{t}$ since it was designed precisely to force the [math]-level set to stop close to the edges. As it has been observed in [28], its meaning does not come from the geometry of $v$ but only from the configuration of the front $\Gamma_{t}$. Using one of the classical velocity functions introduced before, as it will be clarified by the numerical tests in Sect. 5, high-order schemes produce unstable results since the numerical approximation of $v$ can start to produce spurious oscillations near the edges where the front should stop. This problem can be solved by adding a limiter as in [5] but here we present a different technique that avoids the use of a limiter and adapts automatically the scheme according to the regularity of the solution in order to produce more accurate results. To this end, we need to extend the image-based velocity function $c(x,y)$ to all the level sets of the representation $v$ in order to give a physical meaning to the speed used in the whole domain.

Following some of the ideas discussed in [28], we want to extend in a simple way the velocity to the whole domain. Our approach exploits the choice of the initial condition $v_{0}$ (which is free) and allows to avoid all the heavy computations required by the numerical approach solution proposed by the authors in [28]. The modification has an interesting interpretation in terms of the method of characteristics, as we will see later in this section, and allows for stable numerical results. Thus, recalling their approach, the first property that the velocity has to satisfy is:

Property 2.1**.**

Any external (image-based) speed function that is used in the equation of motion written for the function $v$ should not cause the level sets to collide and cross each other during the evolutionary process.

To present the main idea, let us consider the signed distance to the initial [math]-level set as the initial representation function. This is the classical choice

[TABLE]

where $\Omega_{0}$ is the internal region delimited by $\Gamma_{0}$. Therefore, with this choice we can define the velocity extension as follows:

Property 2.2**.**

The value of the speed function $c(x,y)$ at a point $P$ lying on a level set $\{v=C\}$ is exactly the value of $c(x,y)$ at a point $Q$, such that the point $Q$ is a distance $C$ away from $P$ and lies on the level set $\{v=0\}$.

Note that the point $Q$ is uniquely determined whenever the normal direction in $P$ is well defined. In fact, $Q=P-c(x,y)\eta(P)$, where $\eta$ is the outgoing normal and this will provide a definition also when the level sets are non convex. To develop a formal argument we will assume that the normal is sufficiently regular.

In order to compute the point $Q$ on the [math]-level set associated to each point $P$ of any level set, the authors introduce in [28] on pages 162-164 a neighborhood of the 0-level set called Narrow-Band. They fix the width $\delta$ of the band and the related number of iterations $l$ needed for the [math]-level set to reach the boundary of the narrow-band (this is done by the procedure explained at page 164 of [20]). This procedure forces to modify the representation function since the solution is updated only inside the band and a reinitialization step becomes strictly necessary in order to restore the meaning of the distance function. Quoting from their paper, such procedure either requires at least $O(N^{3})$ computations for a grid of $N$ points in each direction or requires to compute the solution of an associated stationary eikonal equation. In order to reduce the computational complexity, the method proposed here is based on a direct assignment of the associated point on the [math]-level set, which requires only $O(N^{2})$ operations to compute the modified velocity at each iteration, and greatly simplifies the problem. Moreover, it allows for the use of representation functions (i.e. initial conditions) even more regular than the signed distance function. The idea is straightforward and is based on the fact that the evolution is oriented in the normal direction to the front. If the reciprocal position of the level sets is also known (that is why we must choose wisely the initial condition) and we make all the points in the normal direction to the [math]-level set evolve according to the same law it is reasonable to expect that all such points will keep their relative distance unchanged as time flows.

To illustrate and motivate our modification, let us still consider the distance to $\Gamma_{0}$ (6) as initial condition and let us consider the shrinking case as example. Then, by construction, all the $C$-level sets are at a distance $C$ from the [math]-level set, as stated by Property 2.2. Hence, if we consider a generic point $(x_{c},y_{c})$ on a $C$-level set, then it is reasonable to assume that the closest point on $\Gamma_{0}$ should be

[TABLE]

Therefore, it seems natural to define the extended velocity $\widetilde{c}(x,y)$ as

[TABLE]

which coincides with $c(x,y)$ on the [math]-level set, as it is needed. The same approach can be applied as long as the initial distance between the level sets is known. In that case, if we want higher regularity to the evolving surface, which would be preferable in the case of high-order schemes such as those we use in the numerical tests, we can define an appropriate initial condition, for example, by simply rotating a regular function in one space dimension. More precisely, let us consider a regular function $\overline{v}_{0}:\mathbb{R}^{+}\to\mathbb{R}$ such that $\overline{v}_{0}(r_{0})=0$, where $r_{0}$ is the radius of the initial circle $\Gamma_{0}$ (e.g. the right branch of a parabola centered in the origin), and let us define $v_{0}(x,y)$ rotating its profile, that is

[TABLE]

Then, it is clear that the C-level sets of $v_{0}$ are located at a distance

[TABLE]

from the [math]-level set and, according to our previous remarks, they should keep this property as time evolves. Consequently, also in this case we can define

[TABLE]

More details on the function $d(v)$ will be given in Sect. 2.2. For simplicity, in the last construction we assumed the representation function to be centered in the origin, but it is straightforward to extend the same procedure to more general situations. Note also that if we have only one object to be segmented (or we are considering the shrinking from the frame of the picture) we can always use a representation function centered in the origin since we can choose freely the domain of integration, given by the pixels of the image.

2.2 Motivations of the new velocity function

The modification of the velocity $c(x,y)$ into $\widetilde{c}(x,y,v,v_{x},v_{y})$ defined in (11) with $d(v)=0$ if $v=0$, is to follow the evolution of the [math]-level set and then to define the evolution on the other level sets accordingly. This allows to give a geometrical interpretation of the new velocity and to establish some properties that will also guarantee existence and uniqueness for the first order evolutive problem. The new velocity can be seen as a characteristic based velocity. As a first step let us analyze the characteristics of the equation, in these computations we assume that we always have the necessary regularity. In particular, we assume $v\in C^{2}(\Omega)$ (or at least $C^{2}$ in space and $C^{1}$ in time) and $c(x,y)\in C^{1}(\Omega)$. We introduce the notations of the vectors $z:=(x,y)$ and $p:=(p_{1},p_{2})=(v_{x},v_{y})$ that will be used only in this section. Using these notations, the Hamiltonian in our case can be written as

[TABLE]

Let us introduce the method of characteristics, writing the usual system

[TABLE]

where $\dot{f}$, $f=z,v,p$ denotes the derivative with respect to the variable $s$, $\nabla_{p}H$ the gradient with respect to $p$, $H_{v}$ the partial derivative with respect to the (scalar) value $v$, and $\nabla H$ the usual spatial gradient of $H$ already introduced in Sect. 1. In our case, defining for brevity the point $(\xi,\zeta):=(x-d(v)\frac{p_{1}}{|p|},y-d(v)\frac{p_{2}}{|p|})$, we obtain

[TABLE]

and analogously $\frac{\partial H}{\partial p_{2}}$, so that we obtain

[TABLE]

Therefore, the system (13) becomes in our case the following:

[TABLE]

where $d^{\prime}(v)$ denotes the derivative of the distance $d$ with respect to the value $v$.

Now, choosing $d^{\prime}(v)$ such that

[TABLE]

we have the final system

[TABLE]

which states that, as long as the function $\widetilde{c}$ remains smooth enough ($\frac{\partial c}{\partial\xi}\approx 0$ and $\frac{\partial c}{\partial\zeta}\approx 0$), the characteristics are basically directed in the normal direction and along them both the height and the gradient are preserved. Looking at the third relation of (18) and at the choice (17), since $p(s)\equiv p(0)=\nabla v_{0}$ along the characteristics, we can choose simply

[TABLE]

which is the trivial case with the function $d(v)=v$ and also for $d(v)$ given by the previous definition (10). In fact, using the inverse function theorem, we have

[TABLE]

with $w$ such that $\overline{v}_{0}(w)=v$. Moreover, recalling the definition (9), we can compute

[TABLE]

and then it is enough to consider $(x,y)$ such that $w=\sqrt{x^{2}+y^{2}}$. From a numerical point of view, we compute the function $d(v)$ analytically, by using (19) that exploits our knowledge of the initial condition $v_{0}$, e.g. given by (9). Thanks to the previous computations, we reached a good understanding of the nature of the evolution given by (2)-(11), but we still have not justified the main motivation that led us to define (11), that is to make all the level sets of $v$ evolve according to the same law. More precisely, we have to show that, if we consider the evolution of two points on the same characteristic but on two different level sets, say the [math]-level set $z^{0}(s)$ and a generic level set $z^{\ell}(s)$, then their relative distance (along the characteristic) does not change during the evolution. This fact would imply that, if we choose the level sets of $v_{0}$ to be such that

[TABLE]

then the points $\underline{z}(s):=z^{\ell}(s)-d(v(z^{\ell}))\frac{p(z^{\ell})}{|p(z^{l})|}$ are always on the [math]-level set of $v$. In order to prove this last statement, let us proceed by a simple differentiation, dropping the dependence on $z^{\ell}$ for brevity,

[TABLE]

Recalling the relations in the system (18) with respect to $\dot{v}(s)$ and $\dot{p}(s)$, we can write

[TABLE]

This calculation shows that the points $\underline{z}(s)$ and $z^{\ell}(s)$ evolve according to the same law along characteristics. Note that if (22) holds then $\underline{z}(s)\equiv z^{0}(s)$ till the characteristics do not cross. In fact, computing the total derivative with respect to $s$, we have

[TABLE]

since the points of $z^{\ell}(s)$ are on the same level set, and so are those of $\underline{z}(s)$, as we wanted. This directly implies that the points $\underline{z}:=\left(z-d(v)\frac{\nabla v}{|\nabla v|}\right)$ are on the [math]-level set of $v$ as long as the gradient is preserved.

In order to reduce the computational cost and to simplify the implementation, as will be better explained in Sect. 4, we decided to split the computation of the Hamiltonian into two steps. First we compute the points $\underline{z}$ and then we update the numerical solution of the full problem (2) with the velocity (11) via the following simplified problem with isotropic velocity

[TABLE]

where, as usual, $t_{n}=t_{0}+n\Delta t$, and $\Delta t$ is the time step. Using this procedure the dependencies on the gradient and on the value of $v$ are frozen at every iteration, leaving just an explicit dependence on the variable $t$. This is why we consider velocities depending only on space variables when defining the numerical schemes.

Let the original velocity $c(x,y)$ be Lipschitz continuous (this is the classical assumption) and let us denote by $L_{c}$ its Lipschitz constant. We conclude showing that our modified velocity $\tilde{c}$ is still Lipschitz continuous. This point is important to guarantee existence and uniqueness of the characteristics and of the viscosity solution for the evolutive problem (2) driven by the new velocity. Let $z=(x,y)$ and $z^{\prime}=(x^{\prime},y^{\prime})$ be two points in the plane and $\underline{z}$ and $\underline{z}^{\prime}$ the corresponding points used in the definition of the new velocity $\tilde{c}$. Provided that the solution and the normal vector are Lipschitz continuous, we have

[TABLE]

where $\eta$ represents the normal unitary vector at the point and $C_{1}$, $C_{2}$ are two appropriate constants (remember that the distance from the level set stays bounded during the evolution).

3 The Adaptive Filtered Scheme

In this section we will introduce and illustrate the AF scheme we will use to approximate the viscosity solution of the problem (2). It is important to note that the name does not refer to filtering a noise as is common in the imaging community but rather to the presence of a filter function as we will explain later in this section. For more details on the AF scheme, see [19]. We assume that the Hamiltonian $H$ and the initial data $v_{0}$ are Lipschitz continuous functions in order to ensure the existence and uniqueness of the viscosity solution [17]. For a detailed presentation of uniqueness and existence results for viscosity solutions, we refer the reader to [17] and [4].

Now, let us define a uniform grid in space $(x_{j},y_{i})=(j\Delta x,i\Delta y)$, $j$,$i\in\mathbb{Z}$, and in time $t_{n}=t_{0}+n\Delta t$, $n\in[0,N_{T}]$, with $(N_{T}-1)\Delta t<T\leq N_{T}\Delta t$. Then, we compute the numerical approximation $u^{n+1}_{i,j}=u(t_{n+1},x_{j},y_{i})$ with the simple formula

[TABLE]

where $S^{M}$ and $S^{A}$ are respectively the monotone and the high-order scheme dependent on both space variables, $F:\mathbb{R}\to\mathbb{R}$ is the filter function needed to switch between the two schemes, $\varepsilon^{n}$ is the switching parameter at time $t_{n}$, and $\phi_{i,j}^{n}$ is the smoothness indicator function at the node $(x_{j},y_{i})$ and time $t_{n}$, based on the 2D-smoothness coefficients defined in [33] and briefly recalled later on in this section. The AF scheme here introduced is convergent, as proven in [18]. In the sequel the gradient components will be denoted by the usual notation $(p,q)$ and $p^{+}$, $p^{-}$ will be right and left discrete derivatives with respect to $x$ (similar notations apply to $q$ that denotes the discrete partial derivative with respect to $y$).

The two schemes composing the AF scheme can be freely chosen, provided that they satisfy the following assumptions:

Assumptions on $S^{M}$: The scheme is consistent, monotone and can be written in differenced form

[TABLE]

for a Lipschitz continuous function $h^{M}(x,y,p^{-},p^{+},q^{-},q^{+})$, with $D_{x}^{\pm}u^{n}_{i,j}:=\pm\frac{u^{n}_{i,j\pm 1}-u^{n}_{i,j}}{\Delta x}$ and $D_{y}^{\pm}u^{n}_{i,j}:=\pm\frac{u^{n}_{i\pm 1,j}-u^{n}_{i,j}}{\Delta y}$.

Assumptions on $S^{A}$: The scheme has a high-order consistency and can be written in differenced form

[TABLE]

for a Lipschitz continuous function $h^{A}(x,y,p^{-},p^{+},q^{-},q^{+})$ (in short), with

[TABLE]

Example 3.1**.**

As examples of monotone schemes in differenced form satisfying the hypotheses stated before, we can consider the simple numerical hamiltonian

[TABLE]

for the eikonal equation

[TABLE]

or, for more general equations also depending on the space variables, we can use the 2D-version of the local Lax-Friedrichs hamiltonian

[TABLE]

with

[TABLE]

where $I(a,b):=[\min(a,b),\max(a,b)]$. This scheme is monotone under the restrictions $\frac{\Delta t}{\Delta x}\cdot\alpha_{x}+\frac{\Delta t}{\Delta y}\cdot\alpha_{y}\leq 1$.

Example 3.2**.**

An example of numerical hamiltonian $h^{A}$ satisfying the assumptions required is the Lax-Wendroff hamiltonian

[TABLE]

where $D^{\pm}_{x}u$, $D_{x}u$, $D^{2}_{x}u$ are, respectively, the usual one-sided and centered one-dimensional finite difference approximations of the first and second derivative in the $x$-direction (analogously for the $y$-direction), whereas for the mixed derivative we use

[TABLE]

Note that the derivatives of $H$ can be computed either analytically or by some second order numerical approximation. In particular, to compute the derivative $H_{x}$, we can simply use

[TABLE]

and analogously for $H_{y}$.

For more details on the construction of $S^{M}$ and $S^{A}$ and other examples of possible numerical hamiltonians, see [33, 18].

In our approach, in order to couple the two schemes, we need to define three key quantities:

The filter function F, which must satisfy

1. (a)

$F(r)\approx r$ for $|r|\leq 1$ so that if $|S^{A}-S^{M}|\leq\Delta t\varepsilon^{n}$ and $\phi_{i,j}^{n}=1\Rightarrow S^{AF}\approx S^{A}$, 2. (b)

$F(r)=0$ for $|r|>1$ so that if $|S^{A}-S^{M}|>\Delta t\varepsilon^{n}$ or $\phi^{n}_{i,j}=0\Rightarrow S^{AF}=S^{M}$.

Several choices for $F$ are possible, different for regularity properties. In this paper, we will consider the discontinuous filter already used in [5] and defined as follows:

[TABLE]

which is clearly discontinuous at $r=-1,1$ and satisfies trivially the two required properties. 2. 2.

If we want the scheme (28) to switch to the high-order scheme when some regularity is detected, we have to choose $\varepsilon^{n}$ such that

[TABLE]

in the region of regularity at time $t_{n}$, that is

[TABLE]

Proceeding by Taylor expansion for the monotone and the high-order Hamiltonians, by (39) we arrive to a lower bound for $\varepsilon^{n}$. The simplest numerical approximation of that lower bound is the following

[TABLE]

in which we have used the usual notation for the gradient, i.e. $(p,q):=(v_{x},v_{y})$ and

[TABLE]

The definition of $\widetilde{h}^{M}_{p^{-}},\widetilde{h}^{M}_{q^{+}},\widetilde{h}^{M}_{q^{-}}$ follows from (42) in an analogous way. All the derivatives of $H$ are computed at $(x,y,D_{x}u^{n},D_{y}u^{n})$ and the finite difference approximations around the point $(i,j)$, using $K>\frac{1}{2}$. See [33] for more details. 3. 3.

For the definition of a function $\phi$, needed to detect the region $\mathcal{R}^{n}$, we require

[TABLE]

with $I_{i,j}:=[x_{j-1},x_{j+1}]\times[y_{i-1},y_{i+1}]$. In order to proceed with the construction, we split the cell $I_{i,j}$ into four subcells, denoted by the superscript ‘$\vartheta_{1}\vartheta_{2}$’, for $\vartheta_{1},\vartheta_{2}=+,-$, according to the shift with respect to the center $(x_{j},y_{i})$. Then, extending the classical WENO approach proposed in [23] to multiple spatial dimensions, we measure the regularity of the solution inside each subcell by computing the smoothness coefficients as rescaled $L^{2}$ norms of the Lagrange polynomial $P_{k}^{\vartheta_{1},\vartheta_{2}}(x,y)$ interpolating the values of $u^{n}$ on the considered stencil, that is

[TABLE]

where $|\vartheta|$ denotes the number of ‘$-$’ in $(\vartheta_{1},\vartheta_{2})$, for $k=0,1$. Note that we have dropped the dependence on the time step $t^{n}$ for brevity and we have used the shorter notation $u_{[t,s]}$ to denote the multivariate undivided difference of $u$ of order $t$ in $x$ and $s$ in $y$. The previous formula can be used to obtain all the needed quantities as long as the following ordered stencils are used to compute the undivided differences

$\bullet$ $\mathcal{S}_{0}^{--}=\{x_{j-1},x_{j},x_{j+1}\}\times\{y_{i-1},y_{i},y_{i+1}\}$, $\mathcal{S}_{1}^{--}=\{x_{j},x_{j-1},x_{j-2}\}\times\{y_{i},y_{i-1},y_{i-2}\}$;

$\bullet$ $\mathcal{S}_{0}^{+-}=\{x_{j+1},x_{j},x_{j-1}\}\times\{y_{i-1},y_{i},y_{i+1}\}$, $\mathcal{S}_{1}^{+-}=\{x_{j},x_{j+1},x_{j+2}\}\times\{y_{i},y_{i-1},y_{i-2}\}$;

$\bullet$ $\mathcal{S}_{0}^{++}=\{x_{j+1},x_{j},x_{j-1}\}\times\{y_{i+1},y_{i},y_{i-1}\}$, $\mathcal{S}_{1}^{++}=\{x_{j},x_{j+1},x_{j+2}\}\times\{y_{i},y_{i+1},y_{i+2}\}$;

$\bullet$ $\mathcal{S}_{0}^{-+}=\{x_{j-1},x_{j},x_{j+1}\}\times\{y_{i+1},y_{i},y_{i-1}\}$, $\mathcal{S}_{1}^{-+}=\{x_{j},x_{j-1},x_{j-2}\}\times\{y_{i},y_{i+1},y_{i+2}\}$.

Since these coefficients are such that

$\bullet$ $\beta_{k}=O(\Delta^{2})$, with $\Delta:=\max\{\Delta x,\Delta y\}$, if the solution is smooth in $\mathcal{S}_{k}$;

$\bullet$ $\beta_{k}=O(1)$ if there is a singularity in $\mathcal{S}_{k}$,

according to the usual WENO procedure we weight the obtained information and focus on the ‘inner’ stencil, denoted by the subscript ‘[math]’, by computing

[TABLE]

where $\alpha^{\vartheta_{1}\vartheta_{2}}_{k}=\frac{1}{(\beta^{\vartheta_{1}\vartheta_{2}}_{k}+\sigma_{\Delta})^{2}}$, with $\sigma_{\Delta}=\Delta x^{2}+\Delta y^{2}$, which represents the measure of smoothness of the solution in each subcell. Once we have computed the four indicators, we couple the information by defining

[TABLE]

which, as can be shown by exploiting the properties of the coefficients $\beta_{k}$, is such that

[TABLE]

where $\rho$ is a point of discontinuity in the gradient. At this point, in order to reduce the amplitude of the oscillations around the optimal value $\frac{1}{2}$ in regular regions (the $O(\Delta)$ term), we use the mapping first introduced in [22] to propose a modification of the original WENO procedure, called M-WENO, that is

[TABLE]

which, using Taylor expansion around $\frac{1}{2}$, directly gives $\omega^{*}=\frac{1}{2}+O(\Delta^{3})$ when the solution is regular in $I_{i,j}$. Finally, in order to define our function $\phi$, it is enough to take

[TABLE]

with $M<\frac{1}{2}$ (e.g. $M=0.1$ in the tests reported in Sect. 5), a number that can also depend on $\Delta x$ and $\Delta y$ . For more details and other possible constructions for the definition of $\phi$, we refer the interested reader to [33, 18].

Remark 3.3**.**

Choosing $\varepsilon^{n}\equiv\varepsilon\Delta x$, with $\varepsilon>0$ and $\phi_{i,j}^{n}\equiv 1$, we get the Filtered Scheme of [5], so here we are generalizing that approach to exploit more carefully the local regularity of the solution at every time $t^{n}$ and cell $I_{i,j}$.

4 Numerical implementation of the modified LS method

Before illustrating the numerical tests, let us first give some comments on the numerical schemes composing the AF scheme adopted for the tests in Sect. 5. The main issue concerning the local Lax-Friedrichs and the Lax-Wendroff schemes defined by (3.1) and (3.2), respectively, is the need to compute the one-directional velocities $H_{p}$ and $H_{q}$ which depend also on $\frac{\partial c}{\partial\xi}$ and $\frac{\partial c}{\partial\zeta}$, as visible in (15). Moreover, in order to implement the local Lax-Friedrichs scheme we should be able to compute the maximum of $|H_{p}|$ (resp. $|H_{q}|$) uniformly with respect to $p$ (resp. $q$), which is a very intricate matter due to the (possible) low regularity of $\widetilde{c}$. In fact, if we focus on the usual behavior of $c(x,y)$ in the proximity of a relevant edge, we can expect the derivatives $\frac{\partial c}{\partial\xi}$ and $\frac{\partial c}{\partial\zeta}$ to be really big. This is not surprising since the front decelerates rapidly in the neighborhood of an edge. In addition, in order to solve the full model (2)-(11), we should take into account also the remaining dependence of $H(x,y,v,\cdot,\cdot)$ when deriving the second-order Lax-Wendroff scheme and, clearly, the formula to compute the threshold $\varepsilon^{n}$. Finally, concerning the Courant-Friedrichs-Lewy (CFL) condition, it is necessary to compute $\max\{|H_{p}|,|H_{q}|\}$ with the full formula (15). Consequently, $\lambda$ could be excessively small due to the low regularity of $\widetilde{c}$. In this latter case, we would clearly need an adaptive mesh refinement technique to reduce the computational cost.

In order to avoid most of these complications in the numerical implementation, we choose to approximate the solution of the simplified problem (26), adjusting the velocity $\widetilde{c}$ according to (11) at each time step. Using the simplified problem (26), we can use the simple relation

[TABLE]

avoiding to take the maximums over all the possible values of $p$ and $q$, which is instead required for the resolution of the full problem (2)-(11) with anisotropic velocity $\widetilde{c}$. Analogous comment holds for $H_{q}$. Lastly, from the numerical point of view, the use of this simplification brings another fundamental consequence: when we apply the numerical schemes to solve (26), we are considering, formally, a problem with bounded velocities $\max\{|H_{p}|,|H_{q}|\}\leq 1$. This implies that we can choose the following CFL condition:

[TABLE]

using the relation (50), which is a less restrictive condition with respect to the original one coming from the full problem (2)-(11).

In the following, we will use the same notations introduced in Sect. 3, except for the number of time steps $N_{T}$, which will be replaced by the total number of iterations $N_{i}$ used by the scheme, since now we are looking for an asymptotic solution (in some stationary sense). The maximum number of iterations, which is fixed at the beginning of the procedure, will be denoted by $N_{\max}$.

Let us give some details on the precise numerical implementation, commenting the main procedures involved in the (sketched) Algorithm 1.

Let us set the parameters of the simulation, which are the power $\mu$ in (4), the number of iterations $K_{reg}$ of the heat equation for the Gaussian filter, the tolerance $tol>0$ of the stopping criterion, the amplitude of the pixels $(\Delta x,\Delta y)$ and, subsequently, the time step $\Delta t$ according to the CFL condition (51).

Then, at each iteration $n=0,\dots,N_{i}$, which has to be interpreted in the sense “until convergence” (note that $N_{i}$ is not known a priori, but depends on the stopping criterion described in Step 3 and can be equal to $N_{max}$ in case of not convergence of the scheme), we repeat the following steps.

Step 1. For $i=0,\cdots N_{y},j=0,\cdots,N_{x}$, with $(N_{x}+1)\times(N_{y}+1)$ the size of the input image, we precompute the matrix $\widetilde{c}(x_{j},y_{i},u_{i,j},D_{x}u_{i,j},D_{y}u_{i,j})$ at the beginning of each iteration using central finite difference approximations for the first order derivatives $D_{x}u_{i,j}$ and $D_{y}u_{i,j}$. Note that the quantities only depend on $(i,j)$ also through $u$. Clearly, this method is valid only as long as the representation function $u$ remains smooth at all the level sets, and should be justified in the case of singular edges (although we will not pursue this precise matter). Moreover, in general the point

[TABLE]

is not a point of the grid $(x_{j},y_{i})$.

To reconstruct the correct value (or at least a reasonable approximation) there exist different possible implementations. For example, a simple bilinear reconstruction from the neighboring values

[TABLE]

where we have used the notation

[TABLE]

with the other cases following an analogous definition. Another possibility, which we have used in the numerical tests since it seems to give nicer results in terms of the shape of the approximate representation $u$, consists in taking as $\left(x_{j_{u}},y_{i_{u}}\right)$ the point such that

[TABLE]

Note that this construction is well defined only if $|\nabla u_{i,j}|\not=0$. Therefore, we define the updated velocity matrix as

[TABLE]

and we use $\widetilde{c}^{n}$ as an isotropic velocity in the next step.

Step 2. We approximate the problem (26) using the AF scheme (28), with the local Lax-Friedrichs scheme (3.1) as $S^{M}$ and the Lax-Wendroff scheme (3.2) as $S^{A}$. From now on, we will refer to this AF implemented scheme with the acronym AF-LW.

We add homogeneous Neumann boundary conditions to the problem (26) in all our experiments in order to not alter the average intensity of the image.

Step 3. In this last step, we describe how approximating the front $\Gamma_{t}$. Since $\Gamma_{t}$ is a curve, it is composed by points that are not all grid points belonging to our mesh. Hence, in order to approximate the position of the front at each time step $t$, we consider a neighborhood $\theta_{\delta}$ of the front $\Gamma_{t}$ of radius $\delta=\max\{\Delta x,\Delta y\}$. In this way, stopping the evolution as soon as the front ceases to move will be equivalent to require that the neighborhood $\theta_{\delta}$ ceases to move. In order to apply that procedure, at each iteration $n$ we store the values of the points $(x_{j},y_{i})$ such that $u_{i,j}^{n}$ changes sign in a matrix $F_{i,j}^{n}$ (we use the closest points on the grid, that are $(i,j\pm 1)$ and $(i\pm 1,j)$), and set $F^{n}_{i,j}=0$ otherwise. In this way we automatically store the disposition of the front with an error of order $\delta=\max\{\Delta x,\Delta y\}$, i.e. we approximate the position of the front, as desired.

We will continue to do that at each iteration until the matrix $F$ at two consecutive iterations will be “close enough”. For this reason, we consider two stopping rules:

[TABLE]

where $\tau>0$ is the prescribed tolerance a priori chosen,

[TABLE]

where now a dependence on the discretization parameters appears.

In our implementation we use one of the two stopping rules above introduced combined with a condition on the number of allowed iterations (i.e. $n<N_{max}$).

5 Numerical simulations

In this section we present a series of numerical experiments on both synthetic and real images, comparing the results obtained by the AF scheme with those obtained by the simple monotone scheme and a high-order scheme which uses the Total Variation Diminishing (TVD) Runge-Kutta (RK) of third order in time and the WENO scheme of second/third order in space. More in details, for the WENO scheme we use the same efficient implementation suggested in [23], Remark 1 on page 2130, which we extended to the 2D case, adding the improvement presented in [2], which consists in choosing $\sigma=\Delta x^{2}$ instead of $\sigma=10^{-8}$ that is the value used in the original paper [23] by Jiang and Peng. Note that we used this scheme instead of the third to fifth order scheme used in the numerical tests of [23] for comparison reasons, otherwise it would be not comparable by order. The first aim is to show the possible improvements of the modified model with respect to the classical formulation. In fact, after extensive numerical simulations, we noticed that the classical model is not well defined when using high-order schemes, since they can produce heavy oscillations as soon as the Lipschitz constant of the representation function becomes too big. This effect causes the stopping rule to be practically ineffective (independently of the norm used) in most cases particularly when using the AF scheme or the WENO scheme, whereas the simple monotone scheme seems to give always stable results. Note that, when the singularity develops, the representation function becomes more and more vertical as time flows.

The numerical tests illustrated in this section will compare the results also varying the initial datum or varying the norm for the stopping rule defined in (57) and (58). Moreover, for synthetic images, we also vary the space steps $\Delta x$ and $\Delta y$, that we consider equal to each other ($\Delta x=\Delta y$), and, therefore, the total number of pixels.

Now, let us specify the initial condition used in each case. When the velocity is defined by the classical model, in the expansion case (Case a) we use the paraboloid

[TABLE]

where $r>0$ is the radius of the initial circle and $\frac{1}{2}r^{2}$ a value chosen in order to cut the surface from above (therefore we have a flat surface at the numerical boundary),whereas in the shrinking case (Case b) we use the truncated pyramid (with a square or rectangular base depending on the image frame), that is

[TABLE]

where $[a_{x},b_{x}]\times[a_{y},b_{y}]$ is the frame of the image, $-0.2$ is the value at which we truncate the pyramid and $2$ is the steepness of the faces of the surface (the smaller the value, the less steep the faces are). By this choice we use a slightly more regular front with respect to the discontinuous representation that simply changes value crossing the frame of the image, still being able to keep the whole surface outside the region occupied by the objects to segment.

For the modified velocity $\widetilde{c}$, in the expansion case (Case a) we consider two different initial data: the paraboloid as in (59) (Datum 1) or the following signed distance function (Datum 2)

[TABLE]

where $\Gamma^{1}_{0}$ is the usual circle centered in $(0,0)$ with radius $r=0.5$ unless otherwise stated. Instead, for the shrinking case (Case b) we only use the signed distance function

[TABLE]

with $\Gamma^{2}_{0}$ representing the frame of the image. We summarize the different initial conditions in Tab. 1.

In order to give a quantitative evaluation of the performances in addition to the qualitative analysis, in the tables we compare the results in terms of number of iterations $N_{i}$ and relative error in pixels, defined as

[TABLE]

where $P_{ex}$ and $P_{a}$ are the number of pixels inside the exact and approximated boundaries of the object(s), respectively. Note that we can compute the “exact” object only if the background is really smooth (in the synthetic cases it is always uniform), because we usually use a comparison with a “threshold” for the values of $I(x,y)$ in order to select the regions occupied by the object (exact object). Whereas, for the approximated object we will count the pixels for which $u_{i,j}\leq 0$. Moreover, we measure the error also with a closely related quantity, that is

[TABLE]

in order to show some dependence on the discretization parameters.

If the schemes do not converge in the fixed maximum number of allowed iterations $N_{\max}$, we will put a “$-$” inside the tables, in place of $N_{i}$. For all our tests, we will set $N_{\max}=2000$. Moreover, in case the front does not stop correctly on the boundary of the object, thus giving an unstable and unusable result, we will put an “X” in correspondence of the errors column. For each test, we specify all the values of the parameters involved ($\mu$, $K_{reg}$, $tol$ and $\Delta x=\Delta y$ or $\#Nodes$), the norm used in the stopping rule and the chosen velocity function. For all the numerical tests presented in this paper, we use CFL number $\lambda=\max\left\{\frac{\Delta t}{\Delta x},\frac{\Delta t}{\Delta y}\right\}=\frac{1}{2}$, $K=1$ in the formula (2) for the computation of $\varepsilon$, and the velocity function $c_{1}$ defined in (4) for the classical model, referring to it simply as the classical $c$. All the numerical tests have been implemented in language C++, with plots and computation of the errors in MATLAB. The computer used for the simulations is a Notebook Asus F556U Intel Core i7-6500U with speed of 2.59 GHz and 12 GB of RAM.

5.1 Synthetic tests

Let us begin by a simple synthetic example. The main goal here is to compare the behavior of the three schemes with respect to the use of the classical velocity $c$ and the modified one $\widetilde{c}$, and the performances varying the number of grid nodes.

Test 1. Rhombus

For this first test, we perform the simulations only in the case of an expansion, since the results do not vary much in the shrinking case.

The rhombus considered is defined by the equation

[TABLE]

that produces a final front visible in Fig. 1 for each scheme using the velocity $\widetilde{c}$. As clearly visible, this rhombus presents some heavily marked corners, which causes some serious troubles when using the filtered scheme or the WENO scheme for the solution of the classical model. In fact, looking at Tab. 2, we can see that the two high-order schemes achieve convergence, in the sense of the iterative stopping rule (57), but after having lost track of the boundary (the front overcomes the edge of the object and then keeps expanding). This happens also varying the number of grid nodes (from $102$ to $402$). On the contrary, the monotone scheme converges, the two errors reported decrease by refining, still using the same parameters $\mu,K_{reg}$ and same tolerance $tol$. The difficulties of the high-order schemes with the classical velocity $c$ are cleary visible looking at Fig. 2, in which the contour plots of the representation function obtained by the three schemes are visible: the evolution of the level-sets for the monotone scheme expands and finally coincides at the final time (last row on the left) with the boundary of the object. Instead, for the two high-order schemes we can see different times in order to show well when the oscillations on the left and right edges begin to increase causing the loss of the rhombus boundary in the approximation.

For the modified velocity $\widetilde{c}$, the behavior of the schemes is different, see the contour plots in Fig. 3. Note that the yellow level-set on the last row, particularly for the WENO scheme, is not the 0-level set. The effects of the modified velocity on the evolution is clear comparing Figs. 2 and 3. In the first case all the level sets expand towards the boundary of the rhombus, eventually collapsing onto each other. This gives rise to a discontinuity around the front, which causes instabilities for high-order schemes. On the contrary, using $\tilde{c}$ these schemes are stable, since the level sets remain equally spaced around the [math]-level set and the representation function is still Lipschitz continuous. Looking at Tab. 3, we can observe that the monotone scheme converges in a lower number of iterations $N_{i}$ and with lower errors with respect to the correspondent results obtained by using the classical velocity and reported in Tab. 2. Note that we used the same parameters, tolerance and initial datum (the paraboloid, Datum 1) for both velocities to make a fair comparison. The AF-LW scheme converges and give always better results in terms of both errors and number of iterations with respect to the monotone scheme with $c$ or $\widetilde{c}$. Also the WENO scheme converges in a lower number of iterations with respect to the monotone one, even if with greater errors. This is due to the fact that the WENO scheme needs to recompute the velocity $\widetilde{c}$, which makes the errors accumulate (since we use Runge-Kutta of third order, we need to recompute $\widetilde{c}$ three times). Hence, all the three schemes benefit of the new definition of velocity and the AF-LW scheme with $\widetilde{c}$ gives the best performances. Regarding the speed of the schemes, looking at Tab. 4 which contains the CPU times in seconds related to the three schemes needed to obtain the results reported in Tab. 3, the AF-LW scheme needs a longer CPU time with respect to the monotone one, as expected. However, the AF-LW scheme is faster than the WENO scheme and its best performances are reached in a short time, at most about 40 seconds for the last refinement of the grid, compared to the $1028.51$ seconds necessary to the WENO scheme.

A more correct analysis can be made by decreasing the tolerance $tol$ together with the number of nodes $\#Nodes$. In Tabs. 5, 6, and 7 we analyze the behavior of the schemes with respect to the initial data. Looking at Tab. 5 we can note that the monotone scheme is not influenced by the change of the initial datum in terms of number of iterations. In fact, comparing the two columns related to $N_{i}$, only one iteration in the last row is different. With respect to the errors, small changes are visible, with a small improvement using Datum 2 for the last two rows. The AF-LW scheme (Tab. 6) has better performances with lower (or equal) errors and lower (or equal) $N_{i}$ when using Datum 2. For the WENO scheme (see Tab. 7) no changes are visible in terms of $N_{i}$ for all the refinements and both the initial data. In terms of errors, only a small change in the last row is visible, with a slight preference for Datum 1. For all the schemes, the errors decrease when we refine the grid as expected. Comparing the three tables, Tabs. 5, 6, and 7, we can note that the AF-LW scheme get always better results in terms of number of iterations and errors with both initial data with respect to the monotone scheme, expect for the last refinement ($\#Nodes=402$) with Datum 2, in which the monotone scheme seems to be a little bit better, due to round off errors. The WENO scheme get always the greatest errors using both intial data, for the reasons explained before, in a number of iterations which differ from those of the AF scheme for at most one iteration in some cases. This simple synthetic example illustrates very well the limits of a high-order scheme like the WENO one, when some singularities occur, and how the proposed AF scheme is able to deal with them.

5.2 Real tests

In this section we consider real images also coming from biomedical applications. Thanks to the error formulas (63) and (64), we can give a sort of quantitative evaluation of the performances of the schemes in terms of “pixels error” (i.e. the number of pixels composing the area of the object to be segmented).

Test 2. Brain ($340\times 340$ pixels)

The first real test focuses on a biomedical image of a human brain. For this test, we approximate the relevant “external” boundary of the brain via a front expansion (Case a) or a shrinking (Case b), so that we start from inside or outside as visible in Fig. 4.

Looking at Figs. 5 and 6, it is worth to note that the two high-order schemes recognize better the boundary of the object starting from two different initial datum using the same tolerance $tol=0.00001$. The differences and the best resolutions are clearly visible looking at the central part close to the bottom of the brain figures. In fact, the monotone scheme stops too early (see the top-left pictures in both Figs. 5 and 6). In order to get better results, a smaller tolerance parameter is necessary for the monotone scheme, see the top-right pictures in both the considered figures. In that case the monotone scheme can increase its resolution even if with a higher number of iterations and in any case not more accurate than the high-order schemes. Comparing the AF and the WENO schemes, both obtain more accurate results with a lower number of iterations with respect to the monotone scheme. More in details, with initial Datum 1, AF scheme converges in a lower number of iterations, with Datum 2 the WENO scheme uses less iterations to get convergence. In any case, the CPU times are really different and the AF scheme is always much more faster than the WENO scheme (aroud 15-20 times faster), as visible looking at Tab. 8, in which the CPU times in seconds related to the simulations of the three schemes visible in Figs. 5-6 are reported. Clearly, the AF-LW scheme needs more CPU time with respect to the monotone scheme, this is mainly due to the computation of the smoothness indicators, but is still fast and competitive since it needs only one minute and half in the worst case (to get a better result). In Case a, hence, the qualitative evaluation is enough to understand which scheme provides the better results in a reasonable CPU time.

For the shrinking case, a quantitative error evaluation is needed in addition to the qualitative one. Looking at Fig. 7, we can see that all the three schemes recognized the desired more external boundary. Analyzing the errors reported in Tab. 9, we observe that the AF-LW scheme produces lower values in both errors $P$-$Err_{rel}$ and $P$-$Err_{1}$ with respect to the other schemes for both the resolutions considered (input image size $170\times 170$ and $340\times 340$ pixels). The errors in the second row have been computed using the mask visible in Fig. 4 on the right and are related to Fig. 7. In that test of the brain, case b, the WENO scheme obtains errors closer to those of the AF scheme, lower with respect to the monotone scheme for both the considered image sizes. Comparing the CPU times with the two different resolutions, we note that with half size and less than half time we can obtain better accuracy with respect to the $P$-$Err_{1}$ error using the AF-LW scheme instead of the monotone one. Comparing only the two high-order schemes, again the CPU times of the WENO scheme are about $11$ or $20$ times greater than those necessary to the AF scheme.

Test 3. Horse Chess ($184\times 256$ pixels)

We choose the image of the horse piece in the game of chess visible in Fig. 8 and we approximate its boundary from inside (Case a), varying the initial datum, and outside (Case b), always using the modified model with velocity $\widetilde{c}$.

Starting our analysis of the results from the expansion case, looking at Fig. 9 we can note some differences between the schemes around the mouth of the horse, at the top at the beginning of the horse’s mane, and at the bottom left, i.e. at the end of the horse’s mane, where the monotone scheme seems to stop too early. Looking at Fig. 10, we note in addition other differences and worse performances of the monotone scheme, due probably to the specularities inside the image. Regarding the errors reported in Tab. 11, the WENO scheme get lower errors with Datum 1 and $tol=0.00005$, even if with a greater number of iterations with respect to the AF scheme, whereas with $tol=0.000025$ we obtain the inverse situation, i.e. the AF scheme get lower errors but with more iterations with respect to the WENO scheme. Looking at Tab. 12, we can note that using the initial Datum 2 the AF scheme provides always the best performances in terms of both errors $P$-$Err_{rel}$ and $P$-$Err_{1}$. Regarding the CPU times, Tab. 13 shows that the WENO scheme is really slow compared to the other schemes, whereas the AF scheme always converges in less than one minute.

Analyzing the results in the shrinking case, small differences can be noted looking at Fig. 11, particularly between the two high-order schemes. For the monotone scheme, we note that the final front stops just before the boundary of the object. So, also in this case, the only qualitative analysis is not enough. A quantitative analysis is shown in Tab. 14, in which both the high-order schemes are more accurate than the monotone scheme as expected. In this case WENO is the most accurate scheme, having lower values in both errors, even if with a CPU time around 14 times greater than that of the AF scheme (see Tab. 15).

In order to show the effectiveness of our implementation of the modified velocity, in Fig. 12 we collected the final representations obtained by the monotone scheme.

Although in the first case some new fronts arise due to the specularities inside the image, we can still recognize that the gradient of the initial condition is preserved and that the level sets do not collide during the evolution. In particular, in the middle image the distance function to the [math]-level set is still clearly visible. The third image is related to the results obtained by using the classical velocity $c$ in Case b, with the same tolerance used with $\widetilde{c}$ and visible in Fig. 11. We note that with the classical model oscillations appear throughout the whole horse, even if focusing on the [math]-level set the scheme performs well. A better behavior of all the level-sets can be obtained but using a greater tolerance $tol$, e.g. $tol=0.0005$.

Test 4. Grains ($300\times 300$ pixels)

In this test we consider the shrinking front in presence of multiple separate objects to be segmented. The final fronts obtained by the three schemes are visible in Fig. 13, in which almost no differences are visible. Looking at the errors and number of iterations reported in Tab. 16, we can note that the two high-order schemes are more accurate than the first order monotone one, with a number of iterations $N_{i}$ very close between the three schemes, the lowest number is obtained by the AF scheme. In terms of errors, the WENO scheme is the most accurate, even if with a really great CPU time as ever (see Tab. 17). We report in Fig. 14 the contour plots of the initial datum and the final representation obtained by the AF scheme. It is rather interesting to see that, although the representation function has a rather complex evolution and splits in various parts, the final solution is still a distance function to the [math]-level set, now composed by different peaks each relative to a single grain.

Test 5. Geometric shapes ($640\times 480$ pixels)

In some of the previous real tests, we have seen that with a careful tuning of the parameters in the stopping rule, the monotone scheme can get comparable results with respect to those obtained by the high-order schemes, at least from a visible point of view. This is not always possible, especially in more critical situations, e.g. when the difference between the background and the objects we want to segment is less marked, as shown in this test. In Fig. 15 we reported the initial front in red for the shrinking case. The more difficult object to detect will be the ellipse, due to the lighter gray levels closer to the white background.

Looking at the final fronts in Fig. 16, obtained by varying the tolerance in the stopping criterion, we can note that the AF scheme recognizes very well the boundaries of all the objects, differently from the monotone one. In fact, if we stop too early the schemes, with $tol=0.0005$, we note that the monotone scheme still has to conclude the recognition of the “spray” shape below on the left. But if we adopt a smaller tolerance, e.g. $tol=0.0001$ or $tol=0.00005$, in order to give “more time” to the monotone scheme to achieve all the boundaries, the scheme improves the detection of the spray shape, but loses the boundary of the ellipse.

Higher-order approximation of the evolution can clearly give more stability, especially when the contrast is not satisfactory as in the present situation.

This numerical test clearly shows the difficulties that the monotone scheme can encounter. It needs a really complicated manual tuning of the parameters, and nonetheless not always gives good and reliable results (as in this case), whereas the AF scheme overcomes this difficulty thanks to its properties.

Test 6. Pneumonia ($191\times 150$ pixels)

Finally, we conclude our numerical tests considering a Pneumonia image in the expansion case starting from an initial datum here composed by four equal paraboloids (each defined as Datum 1) or cones (defined as Datum 2), placed as visible in Fig. 17, with [math]-level sets composed by circles of radius $r=0.125$.

In this test the behavior of the three schemes is confirmed in terms of CPU time, looking at Tab. 18, and in terms of accuracy, as visible in Figs. 18 and 19, where the final fronts obtained by the high-order schemes recognize better the object starting from the two different initial data, Datum 1 and Datum 2, respectively. In fact, in the first figure, Fig. 18, it is evident looking at the lung on the left; in the second case, Fig. 19, it is clear especially looking at the right lung, mostly for visualizing the better performances of the AF scheme. Moreover, comparing in vertical the two figures, focusing on each scheme with different initial datum, we can note that all the three schemes seems to prefer the distance function (Datum 2).

The good behavior of the AF scheme and of the modified model in the case of a merging fronts is confirmed by the contour plots of the final representations, shown in Fig. 20. Regarding the WENO scheme (last column), some oscillations arise at some point in the evolution and increase as time flows (see the bottom-left part of the final representation using Datum 2). This is why a slightly greater tolerance has been used in that case to obtain the reported result.

6 Conclusions and perspectives

In this work we have proposed a new velocity function for the level-set method in order to improve image segmentation and we have extended to 2D an adaptive filtered scheme originally developed in 1D [18]. We have shown that the use of the new velocity function $\widetilde{c}$ allows to get more accurate results stabilizing high-order schemes as the AF or the WENO schemes, for which the classical velocity introduces some instabilities destroying the convergence. Moreover, the new velocity function can also be applied to the simple monotone scheme getting better results even in the first order approximation and does not require the re-initialization of the front in order to keep an accurate tracking of the 0-level set. From the numerical point of view, the AF scheme is based on two building blocks: a monotone scheme and a high-order scheme. The filter function allows to couple the two schemes in a rather simple way and easily switches from one scheme to the other according to new smoothness indicators. In terms of CPU time, the AF scheme is less expensive than the WENO scheme and offers a good option to improve the accuracy of the monotone scheme. It is interesting to note that these two changes are rather effective for the segmentation of synthetic and real images according to many simulations, some of them are presented in Sect. 5. A qualitative analysis of the resulting segmentations is illustrated by the pictures of the last section whereas a more accurate comparison of the schemes is based on a quantitative analysis of the pixel errors. Moreover, the sensitivity of the AF scheme seems to be rather low with respect to the presence of noise and only few steps of a linear filter are required to obtain the necessary regularization of the gray levels of the input image $I$. Although in this paper the image segmentation is obtained just using the classical first order equation, a possible extension to second order problems can be considered, this extension is motivated by the inclusion of curvature terms in the evolutive equation as done in the literature. The analysis of the adaptive filtered scheme to second order non linear equations goes beyond the scopes of this paper and will be object of a further investigation.

Acknowledgments

We would like to thank the National Group INdAM-GNCS for the financial support given to this research.