Mixed geometry information regularization for image multiplicative denoising

Among traditional variational denoising methods, regularization models based on geometric information exhibit superior performance. Models with TV regularization or minimal surface regularization have the advantage of edge protection, but these models bring bad experimental effects such as loss of image contrast [23]. At the same time, the curvature-based model can effectively maintain contrast, but the edge protection ability is reduced and the computational overhead increases several times [27, 41]. We can find that models based on single geometric information enjoy some advantages but also have disadvantages. Then a natural idea is that we can design a mixed geometric model that can combine their advantages and compensate for each other’s shortcomings.

Suppose the image are degraded by the multiplicative gamma noise pollution level $L=4$. As displayed in Fig. 1, during the multiplicative denoising process, the minimal surface model has lower image contrast than the model based on curvature. In terms of visual effects, curvature-based model can protect image details better but over-smooths the edges. To further improve their denoising performance, we use a mixed geometric information model. As shown in the results in Fig. 1, the mixed geometric information model has the best edge structure protection capabilities and comprehensive results.

Unlike additive noise, the extent of corruption caused by multiplicative noise is linked to the signal strength of its clean image. Based on the above analysis, we design gray level indicator for adaptive denoising. Regions with higher original signal strength indicate a greater degree of pollution, necessitating stronger denoising in such regions. However, multiplicative noise severely destroys almost all information in the clean image, making it challenging to directly obtain a corrupted image. Therefore, to achieve our objective, we use a modified gray level indicator. For the image $u(x,y)$, the gray level indicator $\alpha$ can be expressed as

where $G_{\sigma}$ is the two-dimensional convolution kernel with standard deviation $\sigma$, $M=\sup_{x\in\Omega}\left(G_{\sigma}*u\left(x\right)\right)$. $\alpha$ utilizes a Gaussian convolution kernel to pre-smooth the noisy image and approximate the clean image. This operation makes the indicator less sensitive to noise of scale smaller than $\sigma$. Normalization allows the indicator to be regarded as a weight function. The higher the signal strength, the more attention it warrants, and active minimization of the surface should be prioritized. As depicted in the Fig. 2, the indicator effectively characterizes the pollution degree of the noisy image and accurately captures the spatial variation of the image.

Assuming that the gray image $u(x,y)$ is defined on the domain $\Omega$, and is considered as a regular surface in three-dimensional space. Encouraged by the success of mean curvature in handling additive noise, we incorporate the mean curvature of surfaces into the our proposed model. We reconsider the image surface as the zero level set of the level set function $\phi$ [42], $\phi(x,y,z)=u(x,y)-z$. Therefore, its unit outer normal is

Further, the mean curvature of a surface is defined as the divergence of the unit normal

Clearly, the smoothness of the image can be effectively represented by the curvature. Through such a transformation with the level set method, the value of $\phi$ can be used to rapidly and accurately calculate the mean curvature using methods such as the finite difference method.

Considering that the mean curvature regularization of the surface is used in the model, we regularize the surface area rather than directly limiting the length of the level set curve. This also has the advantage that it avoids the singularity caused by division by zero. Specifically, constraining all contours and regularizing the elastic energy of the level set function, we get the following regularization term

Further, by combining the gray level indicator and fidelity term designed for the multiplicative noise process, we propose a model $E\left(u\right)$ that mixes area and curvature geometric information

where $\alpha(u)$ is the gray level indicator, and $b$ and $\lambda$ are positive constants utilized to balance the fidelity term and the regularization term. To better address the characteristics of multiplicative denoising and achieve the regional adaptive effect, we employ the gray level indicator $\alpha$ from [12]. To further illustrate the rationality of the model, we will further introduce the degradation model of the model, each of which is a well-known classic denoising model.

In order to further illustrate the advantages of mixing area terms with curvature terms, we conducted the following experiments. In particular, the free parameters and adaptive part of the model can degenerate the model into many simpler forms despite the high complexity of the model. (1) $\alpha(u)\neq 0$ and is a constant, $b\neq 0$, degenerate into regularization resembling Euler’s elastica; (2) $\alpha(u)\neq 0$ and is a constant, $b=0$, degenerate into minimal surface regularization; (3) $\alpha(u)$ is not a constant, $b=0$, degenerate into adaptive minimal surface regularization; (4) $\alpha(u)$ is not a constant, $b\neq 0$, it is a multiplicative deduction based on the mixed geometry information of the surface. We next analyze various visual denoising effects of the model.

In the experiment, the noise parameter $L=4$ was initially selected. The traditional minimal surface model exhibits a conspicuous phenomenon of fragmentation, which leads to more false boundaries in the result, which is not conducive to the subsequent analysis. Upon employing the gray level indicator, Fig. 3d illustrates a significant reduction in the staircasing effects of the minimal surface model, and Fig. 3e shows a smoother result. However, a “fish scale phenomeno” emerges, accompanied by noticeable jagged edges at the shape’s periphery. In contrast to previous models, our proposed model effectively restores image details and successfully avoids the staircasing effects. The amalgamation of area and curvature geometric information yields mutual benefits, surpassing the denoising capabilities of each individual geometric feature in multiplicative denoising.

We proceed by taking time as the evolution parameter for advancing the fourth-order equation to facilitate the explicit update of $u$. Initially, we assume the image size is $M\times N$. Let $\tau$ represent the time step and define the grid as follows

Generally $\Delta x=1$ and $\Delta y=1$. Denoting $u_{i,j}$ as the value of $u$ at $(x_{i},y_{j})$, various approximations of the derivative in the $x$ direction and the $y$ direction are listed

We exclusively provide the discrete scheme for the second term in (8), with a similar discretization approach applicable to the remaining terms. The second term is reformulated as

where the forms of $V_{1}$ and $V_{2}$ are

Furthermore, we list an approximation of $\nabla\cdot V$

Therefore, we need to calculate $\left.V_{1}\right|_{i+\frac{1}{2},j}$, $\left.V_{1}\right|_{i-\frac{1}{2},j}$, $\left.V_{2}\right|_{i,j+\frac{1}{2}}$, $\left.V_{2}\right|_{i,j-\frac{1}{2}}$. Taking $\left.V_{1}\right|_{i+\frac{1}{2},j}$ as an example, the discretization of others is omitted for simplicity. Let $\Psi:=\kappa\sqrt{1+|\nabla u|^{2}}$, and presented below is the approximate scheme for $V_{1}$ at $\left(i+\frac{1}{2},j\right)$.

where the minmod function is used to increase the anisotropy of the equation

The approximation of $\nabla u$ and $\kappa$ at $\left(x_{i},y_{j}\right)$ can be expressed as

Let $W=\left(W_{1},W_{2}\right):=\frac{\alpha(u)+b\kappa^{2}}{\sqrt{1+|\nabla u|^{2}}}\nabla u$ and $u_{i,j}^{n}$ represents the $n$th update result of $u_{i,j}$, then the evolution process can be expressed as

where $i=1\cdots M$, $j=1\cdots N$. Discretization of Neumann boundary conditions is given by

The explicit scheme evolves in a small time step size due to the constraints imposed by CFL conditions. This restriction diminishes computational efficiency making iteration in this scheme unfeasible.

While the theoretical foundation of model (7) effectively preserves the texture features, the incorporation of high-order nonlinear terms introduces substantial challenges in devising efficient numerical algorithms. Addressing the numerical instability associated with the explicit finite difference method and aiming to retain disorder in diffusion directions, we employ the additive operator splitting (AOS) algorithm as proposed in [43] for the numerical solution of our model.

To begin, we rewrite the evolution equation (9) as

where

The AOS algorithm decomposes the original diffusion into horizontal and vertical directions diffusion, then obtains $u_{1}^{n+1}$ and $u_{2}^{n+1}$ from $u^{n}$, and finally uses the average value of $u_{1}^{n+1}$ and $u_{2}^{n+1}$ as the iteratively updated $u^{n+1}$. From the form of (23), we can easily calculate the diffusion matrices $D_{x}^{n}$ and $D_{y}^{n}$ in the $x$ and $y$ directions of the size $M\times N$ image, taking the $x$ direction as an example:

where boundary conditions are correspondingly defined as $u_{i,0}^{n}=u_{i,1}^{n},u_{i,N+1}^{n}=u_{i,N}^{n}$, then one can derive each element in the matrix accordingly

The values $g_{i,k+\frac{1}{2}}$ and $g_{i,k-\frac{1}{2}}$ are approximated through forward and backward differences. Subsequently, we adjust the original AOS algorithm and add an additional $F^{n}$ term. Upon acquiring the diffusion matrix, one-dimensional diffusion is independently applied to the rows and columns of $u^{n}$ for $K$ iterations.

where $l\in\left\{x,y\right\}$ and $u^{n+1}$ can be solved in an efficient way.

Following the utilization of the Thomas algorithm for the alternate computation of two variables, the denoising result for iterative updates is derived by averaging the values of $u_{1}^{n+1}$ and $u_{2}^{n+1}$

After the tridiagonal matrices $D_{x}^{n}$ and $D_{y}^{n}$ are obtained through (24), the problem of image denoising is transformed into the problem of solving linear system. The AOS algorithm possesses absolute stability, consistently yielding outstanding experimental results even under large time step sizes. However, in order to maintain the accuracy of the algorithm, it is common to set a smaller time step size. Consequently, in experiments, a delicate balance must be struck between solution accuracy and computational efficiency. The choice of an appropriate time step size involves a tradeoff. Algorithm 1 outlines the AOS algorithm for solving (22).

While the AOS algorithm has been used to successfully solve numerous nonlinear models, it is accompanied by the drawbacks of reduced accuracy and ‘streaking artifacts’ [45, 46]. It is worth noting that in order to solve the problem of energy functional minimization, Shen et al. [40, 47] transformed the gradient flow solution of the original problem by introducing a scalar auxiliary variable, and proposed the scalar auxiliary variable (SAV) method. As a burgeoning numerical algorithm in the field of image processing, SAV represents a general and efficient solution framework with simple numerical scheme, accurate calculation accuracy, and unconditional stability [48, 41, 49].

Redefine the new energy functional (7), which can be rewritten as:

where $L$ is a non-negative symmetric linear operator, typically independent of the function $u$ and $\varepsilon_{1}\left[u\right]$ containing a nonlinear term related to $u$. As all terms in problem (7) are nonlinear operators, we take $\varepsilon_{1}\left[u\right]$ equal to the total energy functional subtract linear term inner product. This choice ensures that $\varepsilon\left[u\right]$ is still equal to the total energy functional, allowing them to be harmonized into a unified form. Furthermore, the $L^{2}$ gradient flow corresponding to the minimization of the functional (30) is

By the chain rule, we get that the derivative of energy with respect to time

Therefore, the gradient descent method effectively maintains the reduction of energy during the evolution process. Obviously, the objective functional possesses a well-defined lower bound, so there exists a positive constant $C\textgreater 0$ such that $\varepsilon_{1}\left[u\right]=E\left(u\right)-\frac{\gamma}{2}\left(u,Lu\right)+C\textgreater 0$. To facilitate analysis, we introduce the scalar auxiliary variable $r=\sqrt{\varepsilon_{1}}$ and present the equivalent form of the gradient flow within the SAV framework

Next, we introduce the first-order scheme of SAV

Although the update of $u$ in (34) may appear intricate, the computational cost primarily focuses on solving a linear system with constant coefficients, and the associated calculation cost is comparatively modest. Moreover, we can efficiently address (34) through the following process:

Substitute (34b) and (34c) into (34a), leading to

where

Multipling the (35) by the $b^{n}$ and integrating, we obtain

Then there holds

Finally, we substitute (38) back into (35) to get $u^{n+1}$, which gives

In summary, the first-order scheme of the SAV algorithm can be summarized as Algorithm 2.

SAV approaches of second-order or higher-order schemes can also be easily constructed, typically exhibiting energy stability. Therefore the selection of the time step size becomes less critical, as it does not necessitate adherence to the CFL condition. In what follows, we introduce a second-order scheme within the SAV framework, employing the Crank-Nicolson scheme.

where the value of $\widetilde{u}^{n+\frac{1}{2}}$ can be approximated using an explicit difference method, as exemplified by the following formulation

Then the iterative update of $u^{n}$ can be derived in a similar process in the first-order scheme. First, substitute (40b) and (40c) into (40a), we get

where

Similarly, multipling (43) by $b^{n}$ and integrating, then there holds

Furthermore, it arrives at

Finally, combining (42) and (45), it leads to the iterative update result of $u^{n+1}$

During the practical evolution process, the energy functional may experience rapid decay within a short time period, followed by slight changes in the remaining periods. Therefore, it becomes imperative to opt for a small step size when the energy undergoes rapid changes and a relatively larger step size during periods of energy stability. To address this, we implement an adaptive time step strategy according to [47],

where $\rho$ is a default safety coefficient, $tol$ is a constant representing reference tolerance and $e$ is the relative error in the iteration, $\tau_{min}$ is the minimum time step size, $\tau_{max}$ is the maximum time step size.

We now give the unconditional energy stability properties of the numerical scheme.

Multiply the three equations of (34) by $\mu^{n+1}$, $\frac{u^{n+1}-u^{n}}{\tau}$, and $\frac{2r^{n+1}}{\tau}$, respectively. Then integrate (34a) and (34b), and sum up them. We readily derive the following theorem on energy dissipation. The proof follows from the proof given in [47].

Similarly, multiply the three equations of scheme (40) by $\mu^{n+\frac{1}{2}}$, $\frac{u^{n+1}-u^{n}}{\tau}$, and $\frac{r^{n+1}+r^{n}}{\tau}$, respectively, integrate the first two equations, and finally sum up all the three equations. Using the linear and symmetric properties of $L$, we get the following theorem.