Additive Operator Splitting Scheme
for a General Mean Curvature Flow
and Application in Edges Enhancement

Abstract.

Many models that use non-linear partial differential equations (PDEs) have been extensively applied for different tasks in image processing. Among these PDE-based approaches, the mean curvature flow filtering has impressive results, for which feature directions in the image are important.

In this paper, we explore a general model of mean curvature flow, as proposed in [ G.I. Barenblatt, Self-similar intermediate asymptotics for nonlinear degenerate parabolic free-boundary problems that occur in image processing, Proceedings of the National Academy of Sciences of the United States of America (2001)], [ G.I. Barenblatt and J.L. Vazquez, Nonlinear diffusion and image contour enhancement, Interfaces and Free Boundaries (2003)]. The model can be re-arranged to a reaction-diffusion form, facilitating the creation of an unconditionally stable semi-implicit scheme for image filtering. The method employs the Additive Operator Split (AOS) technique. Experiments demonstrated that the modified general model of mean curvature flow is highly effective for reducing noise and has a superior job of preserving edges.

1. Introduction

Partial differential equation (PDE) based techniques have been widely applied for various image processing tasks over the past few decades [30, 4, 32, 7, 14]. Designing numerical schemes that take into consideration accuracy, stability, and computational cost is fundamental for the effective implementation of these schemes.

Local geometric properties of images are usually dominated by two main directions: the direction perpendicular to the edge, called the flow-line, and parallel to the level or isophote, which is aligned with the edge. By selecting the appropriate degree of diffusion in these two directions, several approaches can be achieved [3, 26, 18, 9, 10].

Mean curvature flow is a prominent example where diffusion only occurs along edges. The theoretical properties of this were first studied by Gage, Hamilton and Huisken in the 1980s [16, 17, 15]. It is known that a plane curve moving at a normal speed equal to its curvature will shrink to a point, its shape becoming smoother and circular. In higher dimensions, the phenomena are more complex, with no classification is available yet. In the context of image processing, Osher and Sethian [22, 27] provided a more rigorous view, realization that the iso-intensity contours of an image can be moved under their curvature was achieved. This led to numerous papers where the images were viewed as a set of level contours, then moved under their curvature. Image smoothing by level-set curvature motion [19, 20] preserves edge information by thwarting diffusion in the edge direction. This work showed that, in addition to this basic approach, a natural stopping criterion can also be chosen to prevent over-smoothing a given image.

A general mathematical framework for feature-preserving image smoothing, applicable to gray-level, vector-value (color) images, volumetric images, and movies was achieved by Sochen, Kimmel, and Malladi in [28, 29]. The main idea is to view the image as a two-dimensional manifold embedded in a hybrid spatial-feature space. The authors in [29] showed that many classical geometric flows emerge as special cases in this view, along with a new flow known as the Beltrami flow, which moves a Gray-level image under a scaled mean curvature. By following a different approach, Yezzi in [33] arrived at a similar equation.

As a result of the certain degeneracy of the asymptotic forms of equations mentioned in [20, 29], Barenblatt [5] noted the possibility of constructing a more general class of equations that generalized mean curvature motion and Beltrami flow. An asymptotic treatment in the one-dimensional case of this class of equations has been investigated theoretically in [5, 6]. The authors report that the general mean curvature flow equation forms a sharp step in the vicinity of edges.

The task of smoothing noisy images and enhancing their edges generally involves solving partial differential equations using numerical integration. This process is often the most time-intensive aspect of nonlinear image processing algorithms. Recently, a lot of research efforts have focused on developing numerical methods and techniques for solving equations dependent on curvature in image processing (see, e.g., [24, 25, 21, 34]).

Our objective is to numerically solve the modified general model of mean curvature flow and to confirm the one-dimension theoretical results obtained in [5]. This is achieved through the development of finite difference schemes for the two-dimensional model. We aim to investigate two critical aspects: The phenomenon of edge enhancement within the image, and the model’s efficacy in noise reduction. Using explicit finite difference schemes can be very time-consuming due to the scaling and small time step requirements in order to achieve stability. Therefore, developing fast and unconditionally stable schemes is preferable. The additive operator splitting (AOS) method, introduced by Weickert et al. [31] for the nonlinear diffusion flow, as an unconditionally stable scheme, is effective and efficient. Applying the (AOS) method in the modified general model of mean curvature flow requires writing the model to allow this method’s application, which we will detail in this paper.

The paper is structured as follows: In Section 2, we present the mathematical model associated to a general mean curvature flow for edge enhancement. Section 3 is devoted to the presentation of the numerical processing, particularly splitting the model and spatial discretization are presented, then the proposed numerical scheme for solving the discretized problem are exposed. Section 4 presents the numerical experiments for image processing. Finally, the paper concludes by a short conclusion.

2. A general mean curvature flow

The level set method [22, 27] elegantly handles mean curvature flow. Essentially, it views the curve as the level set $u=0$ of a function $u$ constrained to solve the degenerate diffusion equation:

(1) $${\partial }_{t}u=\left|\nabla u\right|\mathrm{div}\left(\frac{\nabla u}{\left|\nabla u\right|}\right),$$

where $\nabla =({\partial }_x,{\partial }_y)$denotes the gradient operator and $\mathrm{div}$ stands for the divergence operator. To avoid a division by zero, according to [13], a positive parameter is adopted to regularize the term $\left|\nabla u\right|$. Thus, we have:

(2) ${\partial }_{t}u(x,y,t)$ $=$ $\sqrt{1+{\left|\nabla u\right|}^2}\mathrm{div}\left(\frac{\nabla u}{\sqrt{1+{\left|\nabla u\right|}^2}}\right)$

$=$ $\frac{(1+u_y^2)u_{xx}-2u_x{u}_y{u}_{xy}+(1+u_x^2)u_{yy}}{1+u_x^2+u_y^2},$

which called the mean curvature flow equation.

Sochen et al. [29] used a differential-geometric approach with various assumptions regarding the image intensity flux, leading them to derive the following equation for image intensity $u$

(3) $${\partial }_{t}u(x,y,t)=\frac{(1+u_y^2)u_{xx}-2u_x{u}_y{u}_{xy}+(1+u_x^2)u_{yy}}{{\left(1+u_x^2+u_y^2\right)}^2},$$

which called the Beltrami flow (selective mean curvature flow). Here, $x$ and $y$ are the Cartesian coordinates in the image plane, $t$ is time. Thus, according to [20] and [29], image processing is reduced to the solution of the chosen equation under an initial condition $u(x,y,0)=u_0(x,y)$ corresponding to a gray level of the image being processed.

As a result of a certain degeneracy of the asymptotic forms of equations $\left(\text{<a href="\#S2.E2" title="Equation 2 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">2</span></a>}\right)$ and $\left(\text{<a href="\#S2.E3" title="Equation 3 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">3</span></a>}\right)$, Barenblatt in [5], proposed to consider a more general class of equations

(4) $${\partial }_{t}u(x,y,t)=\frac{({\lambda }^2+u_y^2)u_{xx}-2u_x{u}_y{u}_{xy}+({\lambda }^2+u_x^2)u_{yy}}{{\lambda }^{-2\alpha }{\left({\lambda }^2+u_x^2+u_y^2\right)}^{\alpha +1}},$$

where $\alpha$ and $\lambda$ are constant parameters $\left(\lambda \ne 0\right)$. In the other words, it is appropriate to consider a more general flow defined by:

(5) $${\partial }_{t}u={\left(1+{\left|\frac{\nabla u}{\lambda }\right|}^2\right)}^{\frac{1}{2}-\alpha }\mathrm{div}\left(\frac{\nabla u}{\sqrt{1+{\left|\frac{\nabla u}{\lambda }\right|}^2}}\right).$$

Here, the variable $\lambda$ plays the role of a contrast parameter to distinguish between regions to be smoothed and edges to be preserved. For $\left|\nabla u\right|>\lambda$, edges are preserved, signifying the diffusion effect is small; whereas for , the diffusion coefficient has a high amplitude, leading to a strong smoothing effect. In other words, regions where $\left|\nabla u\right|>\lambda$ are considered as edges, and the diffusion process has a low effect. Perona and Malik in [23] have pointed out the problem with the setting of $\lambda$. This presents an interesting problem because a suitable value of the parameter $\lambda$ turns into a general problem of region separation, object extraction, and edge preservation.

Equation $\left(\text{<a href="\#S2.E5" title="Equation 5 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">5</span></a>}\right)$ is regard as a generalization of mean curvature flow and the Beltrami flow. To justify this, consider the special case $\alpha =0$ and $\lambda =1,$ to obtain the mean curvature flow equation $\left(\text{<a href="\#S2.E2" title="Equation 2 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">2</span></a>}\right)$, and we can obtain Beltrami flow equation $\left(\text{<a href="\#S2.E3" title="Equation 3 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">3</span></a>}\right)$ for $\alpha =1$. It can be observed that in the case of $\alpha =\frac{1}{2}$, it is the Perona-Malik diffusion equation with the Charbonier edge indicator function $g(s)=\frac{1}{\sqrt{1+{(s/\lambda )}^2}}$.

The difference of equation $\left(\text{<a href="\#S2.E5" title="Equation 5 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">5</span></a>}\right)$ to mean curvature flow and Beltrami flow is the exponent $\alpha$ on the term ${\left(1+{\left|\nabla u\right|}^2\right)}^{\frac{1}{2}-\alpha }$ on the right-hand-side. This factor is used for the enhancement of edges. Indeed, it controls the speed of the diffusion: if $\left|\nabla u\right|$ has a small mean in the neighborhood of an interior point of a smooth region of an image, and the diffusion is therefore strong. On the other hand, at the position of an edge, $\left|\nabla u\right|$ is large, the diffusion spread is lowered, since ${\left(1+{\left|\nabla u\right|}^2\right)}^{\frac{1}{2}-\alpha }$ is small for $\alpha >\frac{1}{2},$ whereas for $\alpha \le \frac{1}{2}$, the factor becomes much large then $1$, this leads to a higher evolution velocity close to an edge.

The asymptotic treatment of this model, for $\alpha =0$ and $\alpha =1,$ done in the papers [20, 29], shows the enhancement of the intensity contrasts by formation of regions of large intensity gradients, i.e., the normal component of the image intensity gradient becomes quite large. This phenomenon was analyzed, for $\alpha \ge 0,$ in [5, 6], where a further 1-D simplification of the model was proposed in order to focus on the boundary layer where large gradients concentrate. Authors proved that the edge enhancement will take place if any equation of class $\left(\text{<a href="\#S2.E4" title="Equation 4 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">4</span></a>}\right)$ will be used. In [12, 11] the one-dimensional model have been investigate by using the traveling profiles method proposed in [8].

Solving nonlinear PDEs using explicit methods can be very time consuming due to the scaling and small time step requirement. In this paper our goal is to build a fast and reliable method to solve the general mean curvature flow equation (GMCF) $\left(\text{<a href="\#S2.E5" title="Equation 5 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">5</span></a>}\right)$ and it is based on AOS technique by using a semi-implicit finite difference scheme. The mathematical model is achieved by the following additional homogeneous Neumann boundary condition expressing the conservation of energy

$$\frac{\partial u}{\partial n}=0,$$

where $n$ represents the outward vector normal to the boundary of the domain where the image is defined.

3. Additive Operator Splitting scheme for the general mean curvature flow

To numerically solve equation (5), we utilize the Additive Operator Splitting (AOS) technique, originally introduced by Weickert et al. [31], as an efficient, reliable, and unconditionally stable schemes for nonlinear diffusion in image processing. We adapt the AOS technique to the general mean curvature flow (GMCF), which, although distinct from nonlinear diffusion, can be treated similarly.

Splitting the GMCF

Consider the general mean curvature flow compared with the nonlinear diffusion of a gray-level image. Let $u$ be the pixel value. Equation (5) can be rewritten as:

(6) $${\partial }_{t}u=h^{\frac{1}{2}-\alpha }\mathrm{div}\left(\frac{\nabla u}{\sqrt{h}}\right),$$

where

$$h=1+{\left|\frac{\nabla u}{\lambda }\right|}^2=1+{\left(\frac{u_x}{\lambda }\right)}^2+{\left(\frac{u_y}{\lambda }\right)}^2.$$

To simplify Equation (6), we have:

${\partial }_{t}u$ $=\frac{1}{h^{\alpha -1/2}}\left[\frac{1}{\sqrt{h}}\left(u_{xx}+u_{yy}\right)-\frac{1}{2}\left(\frac{u_x{h}_x}{h\sqrt{h}}+\frac{u_y{h}_y}{h\sqrt{h}}\right)\right]$

$=\frac{u_{xx}+u_{yy}}{h^{\alpha }}-\frac{1}{2}\left(\frac{u_x{h}_x}{h^{\alpha +1}}+\frac{u_y{h}_y}{h^{\alpha +1}}\right)$

$=\left(1-\frac{1}{2\alpha }\right)\frac{u_{xx}+u_{yy}}{h^{\alpha }}+\frac{1}{2\alpha }\frac{u_{xx}+u_{yy}}{h^{\alpha }}-\frac{1}{2}\left(\frac{u_x{h}_x}{h^{\alpha +1}}+\frac{u_y{h}_y}{h^{\alpha +1}}\right)$

$=\left(1-\frac{1}{2\alpha }\right)\frac{u_{xx}+u_{yy}}{h^{\alpha }}+\frac{1}{2\alpha }\left\{\frac{u_{xx}+u_{yy}}{h^{\alpha }}-\alpha \left(\frac{u_x{h}_x}{h^{\alpha +1}}+\frac{u_y{h}_y}{h^{\alpha +1}}\right)\right\}$

$=\left(\frac{2\alpha -1}{2\alpha }\right)\frac{u_{xx}+u_{yy}}{h^{\alpha }}+\frac{1}{2\alpha }\left\{\frac{\partial }{\partial x}\left(\frac{u_x}{h^{\alpha }}\right)+\frac{\partial }{\partial y}\left(\frac{u_y}{h^{\alpha }}\right)\right\}.$

We set $g=\frac{1}{h^{\alpha }}$. By reducing the above equation, we obtain a reaction-diffusion form:

(7) $${\partial }_{t}u=\frac{2\alpha -1}{2\alpha }g\mathrm{△}u+\frac{1}{2\alpha }\nabla \left(g\nabla u\right).$$

In this form, the equation is not a pure diffusion equation. It has both an “parabolic” edge-preserving and an “hyperbolic” edge-sharpening terms. In addition, the reaction-diffusion form of equation $\left(\text{<a href="\#S3.E7" title="Equation 7 ‣ Splitting the GMCF ‣ 3. Additive Operator Splitting scheme for the general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">7</span></a>}\right)$ hides the mixed derivative ${\partial }_{xy}u$, there by making it conducive to the AOS approach.

In more concise terms, the equation can be rearranged as:

(8) $${\partial }_{t}u=\left(A_x+A_y\right)u,$$

where $A_x$ and $A_y$ are the following differential operators:

(9) $$\{\begin{array}{c}{A}_x=\frac{1}{2\alpha }\frac{\partial }{\partial x}\left(g\frac{\partial }{\partial x}\right)+\frac{2\alpha -1}{2\alpha }g\frac{{\partial }^2}{\partial x^2}\\ A_y=\frac{1}{2\alpha }\frac{\partial }{\partial y}\left(g\frac{\partial }{\partial y}\right)+\frac{2\alpha -1}{2\alpha }g\frac{{\partial }^2}{\partial y^2}\end{array}$$

Applying the backward difference formula to $\left(\text{<a href="\#S3.E8" title="Equation 8 ‣ Splitting the GMCF ‣ 3. Additive Operator Splitting scheme for the general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">8</span></a>}\right)$ we get,

(10) $$\frac{U^{q+1}-U^q}{dt}=\left(A_x+A_y\right)U^{q+1}.$$

Here, the superscript $q$ corresponds to the current time step, while $q+1$ denotes the next time step. The subscripts $i,j$ indicate the discrete pixel location. $U_{i,j}^q$ are the known values, and $U_{i,j}^{q+1}$ are the unknown values to be determined. Using $U^{q+1}$ on the right side of $\left(\text{<a href="\#S3.E10" title="Equation 10 ‣ Splitting the GMCF ‣ 3. Additive Operator Splitting scheme for the general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">10</span></a>}\right)$, the integration scheme becomes implicit and unconditionally stable,

(11) $$\left[I-dt\left(A_x+A_y\right)\right]U^{q+1}=U^q,$$

where $I$ is the identity matrix. Before proceeding in time, we calculate the values of the edge indicator function $g$, using the known values of $U^q$. Thus, the scheme is only semi-implicit. Although $g$ depends of the gradient of $U$, we treat it like a given function of $(x,y)$, making the governing PDE “quasi-linear”.

Note that equation $\left(\text{<a href="\#S3.E11" title="Equation 11 ‣ Splitting the GMCF ‣ 3. Additive Operator Splitting scheme for the general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">11</span></a>}\right)$ includes a large bandwidth matrix, because all equations related to new pixel values $U^{q+1}$ are coupled. Our aim is to decouple the set $\left(\text{<a href="\#S3.E11" title="Equation 11 ‣ Splitting the GMCF ‣ 3. Additive Operator Splitting scheme for the general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">11</span></a>}\right)$ so that each row and each column of pixels can be handled separately. For this, we rearrange the equation into the following form:

(12) $$U^{q+1}={\left[I-dt\left(A_x+A_y\right)\right]}^{-1}{U}^q={\left[\frac{1}{2}\left(I-2dt{A}_x\right)+\frac{1}{2}\left(I-2dt{A}_y\right)\right]}^{-1}{U}^q$$

Of course, we do not intend to invert the matrix to solve the linear set. This is merely a symbolic form used for further derivation.

Equation $\left(\text{<a href="\#S3.E12" title="Equation 12 ‣ Splitting the GMCF ‣ 3. Additive Operator Splitting scheme for the general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">12</span></a>}\right)$ is interesting form in the sense that the system matrix is decomposed. The problem is that the decomposed system matrix is inside the ${\left[\right]}^{-1}$ operator. Instead, we aim to construct the solution in parts as follows:

(13) $$U^{q+1}=\left({\left[\frac{1}{2}\left(I-2dt{A}_x\right)\right]}^{-1}+{\left[\frac{1}{2}\left(I-2dt{A}_x\right)\right]}^{-1}\right)U^q$$

The problem is that the right hand sides of equations $\left(\text{<a href="\#S3.E12" title="Equation 12 ‣ Splitting the GMCF ‣ 3. Additive Operator Splitting scheme for the general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">12</span></a>}\right)$ and $\left(\text{<a href="\#S3.E13" title="Equation 13 ‣ Splitting the GMCF ‣ 3. Additive Operator Splitting scheme for the general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">13</span></a>}\right)$ are not equal, as can be easily verified. Therefore, we pose the question if there exist a simple variable $m$, when used to multiplay the right hand side of $\left(\text{<a href="\#S3.E13" title="Equation 13 ‣ Splitting the GMCF ‣ 3. Additive Operator Splitting scheme for the general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">13</span></a>}\right)$, would make them equal:

$m\left({\left[\frac{1}{2}\left(I-2dt{A}_x\right)\right]}^{-1}+{\left[\frac{1}{2}\left(I-2dt{A}_x\right)\right]}^{-1}\right)U^q=$

$={\left[\frac{1}{2}\left(I-2dt{A}_x\right)+\frac{1}{2}\left(I-2dt{A}_y\right)\right]}^{-1}{U}^q$

This can be simplified into:

$$m\sum _{i=1}^{i=2}{\left[\frac{1}{2}\left(I-2dt{A}_i\right)\right]}^{-1}{U}^q={\left[\sum _{i=1}^{i=2}\frac{1}{2}\left(I-2dt{A}_i\right)\right]}^{-1}{U}^q$$

Putting $I-2dt{A}_i=B^{-1}$, we find:

$$4m{B}^{-1}=B^{-1}.$$

Thus we have:

$$m=\frac{1}{4}$$

Based on this, to use the additive operator splitting scheme given by equation $\left(\text{<a href="\#S3.E12" title="Equation 12 ‣ Splitting the GMCF ‣ 3. Additive Operator Splitting scheme for the general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">12</span></a>}\right)$, we multiply the right hand side by ${\displaystyle \frac{1}{4}}$, and obtain the following:

$$U^{q+1}=\frac{1}{4}\sum _{i=1}^{i=2}{\left[\frac{1}{2}\left(I-2dt{A}_i\right)\right]}^{-1}{U}^q=\sum _{i=1}^{i=2}{\left[4⁤\frac{1}{2}\left(I-2dt{A}_i\right)\right]}^{-1}{U}^q$$

which is the same as:

$$U^{q+1}={\left(2I-4dt{A}_x\right)}^{-1}{U}^q+{\left(2I-4dt{A}_y\right)}^{-1}{U}^q.$$

Now, the whole idea of this scheme is to bring the equation to a simpler form, allowing us to use efficient block-wise solvers, such as Tri-diagonal matrix algorithm.

Introducing the notations $V={\left(2I-4dt{A}_x\right)}^{-1}{U}^q$ and $W={\left(2I-4dt{A}_y\right)}^{-1}{U}^q$. The solution now includes two components:

(14) $$U^{q+1}=V+W.$$

We finally obtain the equation sets for $V$ and $W$ as follows:

(15) $$\{\begin{array}{c}\left(2I-4dt{A}_x\right)V=U^q\\ \left(2I-4dt{A}_y\right)W=U^q\end{array}$$

The differential operators $A_x$ and $A_y$ are similar, hence we derive here the difference equation for a single row of pixels. Equation for the column of pixels is identical. Consider a row with $N+1$ pixels enumerated from $0$ to $N$.

Discretisation

To discretize the first derivative, we approximate it using a centered finite difference scheme as follows:

$$\begin{array}{c}{\left(\frac{\partial V}{\partial x}\right)}_{i+1/2}=\frac{V_{i+1}-V_i}{dx}+𝒪\left(h_x^2\right)\\ {\left(\frac{\partial V}{\partial x}\right)}_{i-1/2}=\frac{V_i-V_{i-1}}{dx}+𝒪\left(h_x^2\right)\end{array}$$

For the classical finite difference approximation of the second derivative, we utilize the classical three points scheme:

$${\left(\frac{{\partial }^{2}V}{\partial x^2}\right)}_i=\frac{V_{i+1}-2V_i+V_{i-1}}{d{x}^2}+𝒪\left(h_x^2\right).$$

Note that the omitted error terms are of order $h_x^2$.

(16) $$4{\left(A_{x}V\right)}_i=\frac{2}{\alpha }\left[\frac{g_{i+1/2}\left(V_{i+1}-V_i\right)-g_{i-1/2}\left(V_i-V_{i-1}\right)}{h_x^2}+\left(2\alpha -1\right)g_i\frac{V_{i+1}-2V_i+V_{i-1}}{h_x^2}\right]$$

To avoid establishing the values of the edge indicator function $g$ at non-nodal points $i-1/2$ and $i+1/2$, we average the values of the two neighbour nodes:

(17) $$g_{i+1/2}=\frac{g_{i+1}+g_i}{2},\text{ and }{g}_{i-1/2}=\frac{g_{i-1}+g_i}{2}$$

Substituting these averages into the equation (16), we get

(18) $4{\left(A_{x}V\right)}_i$ $=\frac{1}{\alpha h_x^2}[(g_{i-1}+(4\alpha -1)g_i)V_{i-1}-$

$-(g_{i-1}+g_{i+1}+2(4\alpha -1)g_i)V_{i-1}+(g_{i+1}+(4\alpha -1)g_i)V_{i+1}]$

Denote:

$${\beta }_i=\frac{dt}{\alpha h_x^2}\left(g_{i-1}+\left(4\alpha -1\right)g_i\right),\text{ and }{\gamma }_i=\frac{dt}{\alpha h_x^2}\left(g_{i+1}+\left(4\alpha -1\right)g_i\right)$$

From equation (15) and the above expression (18), the finite difference equation becomes

(19) $$-{\beta }_i{V}_{i-1}^{q+1}+(2+{\beta }_i+{\gamma }_i)V_i^{q+1}-{\gamma }_i{V}_{i+1}^{q+1}=U_i^q$$

At each time step, we solve an algebraic system where the matrix $\left(2I-4dt{A}_x\right)$ is an irreducible strictly diagonal dominant matrix and then invertible. Such a system can be efficiently solved using block-wise solvers, such as the Tri-diagonal matrix algorithm. Moreover, by using the Von-Neumann method, we can demonstrate that this semi-implicit scheme is unconditionally stable.

A similar equation holds for $W$ in $y-$dimension for $j=0,1,\mathrm{\dots }M$.

(20) $$-{\beta }_j{V}_{j-1}^{q+1}+(2+{\beta }_j+{\gamma }_j)V_j^{q+1}-{\gamma }_j{V}_{j+1}^{q+1}=U_j^q$$

where

$${\beta }_j=\frac{dt}{\alpha h_y^2}\left(g_{j-1}+\left(4\alpha -1\right)g_j\right),\text{ and }{\gamma }_j=\frac{dt}{\alpha h_y^2}\left(g_{j+1}+\left(4\alpha -1\right)g_j\right)$$

For the first and last nodes, the Neumann boundary conditions hold:

$$V_0=V_1\text{, }{V}_N=V_{N-1}\text{ and }{W}_0=W_1\text{, }{W}_M=W_{M-1}.$$

4. Numerical experiment

In this section, we examine the impact of the model (5) presented in this paper on a diverse range of gray scale images. Additionally, we compare the results obtained with those from the Perona-Malik (for $\alpha =1/2$), and Beltrami flow (for $\alpha =1$) models, highlighting the model’s efficacy in enhancing edges and reducing noise. All numerical experiments were conducted using the AOS method with $h_x=h_y=1$ and $dt=1$.

The images in Fig. 1 demonstrate the effects of diffusion at different values of the contrast parameter $\lambda$. For small values $\lambda$, the image is smoothed while preserving details better than with larger values of $\lambda$. In the first line, even if we see that the PM model has a nice edge enhancement property, we have to take into account the importance of the contrast coefficient $\lambda$, to which the diffusion coefficient is sensitive. Therefore, the problem of determining the value of $\lambda$ is interesting because the appropriate value of $\lambda$ turns into a general problem of region separation, object extraction, and edge preservation. In the second and third lines, it is evident that the images obtained with the Beltrami flow and the general model for $\alpha =2$ are clearer than those obtained with the PM model. Therefore, larger values of the exponent $\alpha$ have done a superior job of preserving edges, even for large values of the contrast parameter $\lambda$.

Although a larger number of iterations were required to achieve an equivalent level of smoothing for larger values of $\alpha$, the improvement in edge preservation can be seen very clearly by looking at the letters in the airplane wing (Fig. 2).

From the images in Fig. 3, one can easily see that the larger exponent $\alpha$ again done a superior job of preserving edges, even for a large number of iterations.

As a conclusion from the experiments presented, it can be said that larger values of $\alpha$ lead to sharper edges, small values lead to smoothed edges, and medium values of $\alpha$ lead to compatibility between the two phenomena. The choice of the value of $\alpha$ remains according to the application and the results to be obtained.

The second simulation concerned the efficacy of the model studied in reducing noise. Experiments were conducted on distorted images with a Gaussian and salt-and-pepper noise. Fig. 4 and Fig. 5 show how noise reduction and edge preservation can be combined using this model. We can observe that it is possible to obtain sharp edges even after a large number of iterations for the large values of $\alpha$.

Figure 1. Comparison of different filtering models. Rows from top to bottom: PM model ($\alpha =0.5$), Beltrami flow ($\alpha =1$), GMCF $\left(\text{<a href="\#S2.E5" title="Equation 5 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">5</span></a>}\right)$ for $\alpha =2$,GMCF $\left(\text{<a href="\#S2.E5" title="Equation 5 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">5</span></a>}\right)$ for $\alpha =3$. Column from Left to Right: $\lambda =1,2,4,6$. Results after 50 iterations.

Figure 2. Comparison of different filtering models after 200 iterations. Rows from top to bottom: PM model ($\alpha =0.5$), Beltrami flow ($\alpha =1$), GMCF $\left(\text{<a href="\#S2.E5" title="Equation 5 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">5</span></a>}\right)$ for $\alpha =2$,GMCF $\left(\text{<a href="\#S2.E5" title="Equation 5 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">5</span></a>}\right)$ for $\alpha =2$. Column from Left to Right: $\lambda =2,4,6$.

Figure 3. Comparison of different filtering models. Rows from top to bottom: PM model ($\alpha =0.5$), Beltrami flow ($\alpha =1$), GMCF $\left(\text{<a href="\#S2.E5" title="Equation 5 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">5</span></a>}\right)$ for $\alpha =2$. Column from Left to Right: Iterations 50,100,200. $\lambda =2$.

Figure 4. Comparison of different filtering models after 50 iterations with $\lambda =4$. From Left to Right: Noisy image, PM model ($\alpha =0.5$), Beltrami flow ($\alpha =1$), GMCF $\left(\text{<a href="\#S2.E5" title="Equation 5 ‣ 2. A general mean curvature flow ‣ Additive Operator Splitting Scheme for a General Mean Curvature Flow and Application in Edges Enhancement" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">5</span></a>}\right)$ for $\alpha =2$.

Figure 5. Comparison of different filtering models after 1000 iterations. First row: (a) original image, (b) noisy image corrupted by Gaussian noise, (c) noisy image corrupted by salt-and-pepper noise. Second and Third rows: Filtred versions of images (b) and (c). First column: PM model ($\alpha =0.5$), Second column: Beltrami flow ($\alpha =1$), Third column: GMCF for $\alpha =1.5$, Forth column: GMCF for $\alpha =2$.

5. Conclusion

In this study, we re-arrange the governing equation for the modified general model of mean curvature flow proposed in [5, 6] into a reaction-diffusion form, where the reaction and diffusion terms are explicitly represented. This reaction-diffusion form enables the development of an unconditionally stable semi-implicit scheme for image filtering. The method is based on the Additive Operator Split (AOS), originally applied by Weickert [31] for the nonlinear diffusion flow. The values of the edge indicator function are used from the previous step in scale, while the pixel values of the next step approximate the flow. This approach leads to a semi-implicit linearized difference scheme. The computational time required for image filtering can be reduced by up to ten times or even more compared to the explicit scheme, depending on the scale step value, with no loss of accuracy. Whereas, using explicit finite difference schemes can be very time-consuming due to the scaling and small time-step requirements in order to achieve stability.

Experiments have demonstrated that the modified general model of mean curvature flow is highly effective in reducing noise and excels in edges preservation for large values of exponent $\alpha$.