Restoring medical images with combined noise base on the nonlinear inverse scale space method

1.   Introduction

During image acquisition, storage, processing, and transmission, the sensor is easily affected by various factors. Besides, external factors may cause interference, resulting in some noise during image transmission. The presence of noise results in a blurred image and loss of intricate details from the original image, complicating any subsequent image analysis. The aim is to remove the effects of noise on the quality of the image, a standard problem in image restoration. Image denoising and restoration are some of the major image problems. The mathematical model builds on prior knowledge to degrade the image, and then uses the reverse process of image degradation to restore the image ^[1].

The additive noise model has been extensively studied by scholars. Most literature studies additive Gaussian noise, assuming that a given original image $u$ is destroyed by additive noise $\eta\sim N(\mu, \sigma^2)$. The task is to recover $u$ from the observed image $f$:

Many effective methods have been developed to deal with this problem, and well-known techniques include wavelet approaches ^[2], learning-based methods ^[3,4,5], and variational approaches ^[6]. One of the learning-based approaches is the one that scholars are now committed to pursuing. In the literature ^[7], the authors proposed a new medical image denoising method to denoise images with Gaussian and impulse noise using end-to-end deep learning. Learning-based denoising methods can generally be used to remove specific types of noise. Variational-based methods are also important methods used for image denoising. Based on ^[6], the researchers proposed some improved methods: high-order total variation ^[8], fractional-order total variation ^[9], and Weberized total variation ^[10,11,12].

In addition, there is another kind of noise called multiplicative noise. In a multiplicative noise model, $f$ represents a noisy image, the original image $u$ is destroyed by multiplicative noise $\eta$:

There are many possible distributions of multiplicative noise. Assume that the multiplicative noise in this paper follows a Gamma distribution with mean one with its probability density function given by

where $1/L$ is the variance. Generally, $L$ takes an integer coefficient. The smaller the $L$ is, the more severe the original image is destroyed by noise. $\Gamma(\cdot)$ is a Gamma function, $\Gamma(L) = (L-1)!$. Recently, there is a lot of interest in multiplicative noise removal. Different from additive noise, the distribution of multiplicative noise is usually non-Gaussian, and it has a multiplicative relationship with the original image. Due to its characteristics, it is more difficult and demanding to remove than additive noise. At present, scientists propose many techniques for eliminating multiplicative noise, such as the Lee method ^[13], scale threshold method ^[14], and variational methods ^{[15,16,17,18,19,20,21,22]}.

However, the true observation image has been contaminated by many different types of noise. The above-mentioned denoising method can only remove a certain noise, or the effect of removing other types of noise is not obvious. In this paper we will focus on removing combined additive and multiplicative noise. Most learning-based denoising methods can remove a specific type (additive or multiplicative) of noise in most cases; whereas, variational-based methods have been studied for combined additive and multiplicative noise. There are variational models available for restoring images affected by combined noises. In ^[23], the author provided restoration of images corrupted with additive, multiplicative, and mixed noise by a new total least squares-based image denoising algorithm. In ^[24], the authors developed a variational model and an efficient nonlinear multigrid approach via a robust fixed-point smoother. Other scholars in ^[25] proposed a variational model combining a fractional-order total variation-norm with the fields of experts-image prior ^[26], employing a fast alternating minimization algorithm, Newton method, and a non-convex optimization iPiano algorithm ^[27].

In ^[28], the authors derived a scale space method for inverse problems, which usually starts from a noisy image and then gradually smooths the image. Different from other methods, the inverse scale space method starts from $u(x, 0) = c_0$ and gradually approximates the noise image $f$ as time grows. The inverse scale space method is related to the regularization theory, especially Tikhonov regularization ^[29]. In ^[28], the inverse scale space method worked effectively on quadratic regularization functions, which can derive a linear evolution equation, but it cannot derive good results for total variational functionals. Then, the author in ^[30,31] derived a general method of the inverse scale space as a limit of the iterative refining process. With the new method, the nonlinear inverse scale space method is implemented on the total variation model. Later, Shi and Osher extended the nonlinear inverse scale space method to the total variational model of multiplicative noise in ^[17].

In fact, real images often contain more than one type of noise, but may contain a mixture of additive and multiplicative noise, making the processing of mixed noise more challenging. The above scholars have only denoised specific types of noise using inverse scale space methods. Motivated by these works, we proposed a new variational model to restore the image damaged by combined noise and generalized the nonlinear inverse scale method on our model. This paper is dedicated to recovering images corrupted by a mixture of additive noise and multiplicative noise. The method provides a unified framework for mixed additive and multiplicative noise.

This paper is organized as follows. In Section 2, we introduce some basic mathematical knowledge and review some classic variational models for removing additive noise or multiplicative noise. In Section 3, we propose our new restoration model and prove that the minimum value exists. In Section 4, the relaxed inverse scale space method is introduced for the solution of the proposed model. Some numeric results of the experiments are presented in Section 5 and the paper concludes with Section 6.

2.   Preliminary knowledge

2.1.   Basic knowledge of bounded variation function space

Image denoising methods based on variational methods and partial differential equations (PDEs) have been widely studied because there is a complete set of mathematical theoretical systems and calculation methods.

The basic theory of bounded variation (BV) function space is the theoretical basis of image processing based on the variational method and PDEs. The total variation (TV) model is a classic model of BV function space theory applied to image denoising. Next, we introduce the theory of bounded variation function space.

Definition 2.1 (Total variation). Let $\Omega \subset R^n$ be a bounded open subset and $u \in L^1_{loc} (\Omega)$ be a locally integrable function. The total variation of $u$ is $TV(u) = \left\{sup{\int _{\Omega}udiv\varphi dx; \varphi = (\varphi_1, \varphi_2, \cdots, \varphi_n) \in C^1_0(\Omega)^n, |\varphi|_{L^\infty (\Omega)} \le 1}\right\}$.\\ If $TV(u) < \infty$ is satisfied, then $u$ is called a bounded variation function, denoted as $|u|_{BV(\Omega)}$. The space composed of all bounded variation functions in $\Omega$ is denoted as $BV(\Omega)$.

Definition 2.2. With a grayscale image $u$ as a two-dimensional image, let it be the Sobolev space $W^{(1, 1)}(\Omega)$, and then $\nabla u \in L^1(\Omega)$ and

The total variation of the two-dimensional image $u$ can be defined as

where $\nabla u = (\frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}) = (u_x, u_y)$, denoting the derivative of a two-dimensional gray image in regard to variables $x$ and $y$, is called the gradient of the image, and $|\nabla u| = \sqrt{u_x^2+u_y^2}$ denotes the modulus of the image gradient.

Definition 2.3. (Bounded variation function space). The space of the bounded variation function can be defined as

Under the definition of natural norm $\parallel u\parallel_{BV} = \ \parallel u \parallel_{L^1} +\ TV(u)$, $BV(\Omega)$ forms a Banach space.

Remark 2.1. The function of bounded variation $\parallel u\parallel_{BV}$ is convex in the $BV(\Omega)$ space, but not strictly convex.

2.2.   Review of some related models

2.2.1.   Total variation-based additive mode

● ROF model

In 1992, Rudin, Osher, and Fatemi used the space of a function of bounded variation establishing a total variation regularized image denoising model (ROF) in ^[6]. The authors constructed the following energy functional:

which meets the following two constraints:

The first formula represents that the mean value of noise is $0$, and the second formula represents that the variance of noise is $\sigma^2$. By introducing the Lagrange multiplier, (2.1) and (2.2) are transformed into a new energy functional without restriction.

where $\lambda$ is a positive parameter, which plays an important balance between smoothing and denoising, $u$ denotes the image to be restored, and $f$ denotes the noisy image.

2.2.2.   Total variation-based multiplicative model

Among the many denoising methods for multiplicative noise, variational models have attracted the attention of many scholars. In this section, we review several classic total variation models for removing multiplicative noise.

● RLO model

In 2003, the first variational model was proposed by Rudin, Lions, and Osher in ^[15]. Assume that the noise in (1.2), the mean noise is $1$, with variance $\sigma^2$. In other words, the noise satisfies (2.3).

Use total variation $\int_{\Omega}\mid \nabla u \mid dxdy$ as the regular term. The model (RLO) is expressed as follows:

where $\lambda$ and $\mu$ are two Lagrange multipliers.

● AA model

In 2008, based on the maximum a posterior (MAP) estimator, Aubert and Aujol in ^[16] derived a non-convex variational model (AA) (2.4) and transformed the solution of the AA model into solving its corresponding Euler-Lagrange equation. This model is suitable for the multiplicative noise of the Gamma distribution with a mean value of $1$.

where $\lambda$ is parameter.

Although the AA model is non-convex, it still proves the existence of its minimum and gives a sufficient condition to ensure the uniqueness of the solution. However, since the AA model is non-convex, its solution may not be the optimal solution, and the solution depends on the initial value selection.

● I-divergence model

In 2009, Steidl and Teuber in ^[18] proposed a convex variational model using I-divergence as the fidelity term to remove multiplicative noise.

Definition 2.4. Let $u \in R^n_{+}$ and $f \in R^n_{+}$. The I-divergence of $f$ from $u$ is defined in ^[32] by

With the property $I(f \parallel u) \ge 0$, if and only if $u = f$, the equal sign holds. Ignore the constant terms, the data fidelity term is simplified to $\int_{\Omega} (u-f\log u) dx$, and the denoising model is as follows:

where $\lambda$ is parameter.

Although I-divergence is generally used as a fidelity term to remove Poisson noise, the author applied Douglas-Rachford splitting techniques and alternating direction methods of multipliers, and proved that it can also effectively remove Gamma noise.

3.   The proposed model and analysis

In this section, we propose a novel variational model combining additive noise (AN) and multiplicative noise (MN). This model is more common in real life, and its solving process is also more challenging. We begin with an observation image contaminated by AN and MN, with the assumption that AN and MN are independent or unrelated. In ^[24,33], the following model is proposed for formulating the image:

where $k_0$ and $k_1$ are constants and $\eta$ represents Gaussian noise with variance $1$.

We propose the following convex model:

where $\alpha_1$ and $\alpha_2$ are parameters that can be adjusted to achieve a balance between the additive and multiplicative parts of the fidelity.

● When $\alpha_1 = 0$, this reduces to the model in ^[6].

● When $\alpha_2 = 0$, this reduces to the model in ^[18].

For the proof of the existence of the model solution, the boundary condition of the solution must be proved ^[34].

Lemma 3.1. (Boundedness). Let $f \in L^2(\Omega)$ and $\underset{\Omega}{\inf} f > 0$, where $\Omega \subset R^2$. If the optimization problem Eq (3.1) has a solution $\hat{u}$, then $\underset{\Omega}{\inf} f \le \hat{u}\le \underset{\Omega}{\sup} f$.

Proof. Define $\alpha = \underset{\Omega}{\inf} f$, $\beta = \underset{\Omega}{\sup} f$. Notice that $u > f$, $u-f\log u$, and $\left(u-f\right)^2$ strictly monotonically increase. Then

Besides, Lemma 1 in ^[35] states that there is

Combining (3.2)–(3.4),

If $u \le \beta$, the equation is valid.

Since $\hat{u}$ is the minimum solution of Eq (3.1), the equation holds when $u = \hat{u}$, thus $\hat{u} \le \beta$.

It is also the case that

Thus $\hat{u}\ge \alpha$.

So there is $\alpha \le \hat{u} \le \beta$, that is, $\underset{\Omega}{\inf} f \le \hat{u}\le \underset{\Omega}{\sup} f$. □

Theorem 3.1. (Existence). Let $f \in L^2(\Omega)$ and $\underset{\Omega}{\inf} f > 0$, where $\Omega \subset R^2$, and then there exists at least one solution in the solution space $BV(\Omega)$ by Eq (3.1).

Proof. Define $\alpha = \underset{\Omega}{\inf} f$, $\beta = \underset{\Omega}{\sup} f$. Because there is $u \equiv \beta \in BV(\Omega)$, the solution space is not empty ^[10]. Consider the optimization problem's minimization sequence ${u_n}\in BV(\Omega)$. According to Lemma 3.1, we have $\alpha \le u_n \le \beta$.

Notice that $\left\{u_n\right\}$ is a bounded sequence in $BV(\Omega)$, $\Omega$ is also bounded, and then

Moreover, by the definition of a minimization sequence, we know there is a constant $C > 0$ and $E(u_n) \le C$. When $u_n = f$, $\int_{\Omega}(u_n-f \log u_n)dxdy$ reaches a positive minimum $\int_{\Omega}(f-f \log f)dxdy$ and $\int_{\Omega}(u_n-f)^2dxdy$ reaches $0$. Since $TV(u_n) < \infty$, we have

From (3.7) and (3.8), there will be evidence that $\left \{ u_n \right \}$ is consistently bounded in $BV(\Omega)$ space with respect to $n$. Because of the compactness of $BV(\Omega)$ space, there are subcolumns in $\left \{ u_n \right \}$ (this remains as $\left \{ u_n \right \}$ for simplicity) and function $u$ in $BV(\Omega)$ space such that $\left \{ u_n \right \}$ converges strongly to $u$ in $L^1(\Omega)$ space. More, assume that $u_n(x, y) \to u(x, y)$, a.e. $\forall x, y \in \Omega$.

The Lebesgue control convergence theorem allows us to deduce

According to the lower semi-continuity of the function, we have the following formula:

Since $\left \{ u_n \right \}$ is a minimization sequence, it is the smallest solution to the optimization problem Eq (3.1), and then Eq (3.1) has at least one solution in its solution space.

The objective function $E(u)$ in the optimization problem Eq (3.1) is convex, and the uniqueness of the solution thus follows directly from the convexity of the objective function. □

4.   Numerical methods

In this section, we will introduce the basic principles and solution process of the inverse scale space method and apply this method to our proposed model. Some scholars have used the inverse scale space method to remove additive noise, and some scholars have used this method to remove multiplicative noise. However, no one has used this method to remove combined noise on our model. Encouraged by this, we have made some improvements on the inverse scale space method to remove combined noise.

4.1.   Inverse scale space method

In the following we generalize nonlinear inverse scale space flow ^[30,31] to the proposed model (3.1), a relaxation method using two evolution equations to approximate the real evolution equations.

The method relies on a convex variance problem:

where $J(u)$ is the regularization term, and $G(u, f)$ is the data fidelity term. We assume that $f$ has a zero mean: $\int_{\Omega}fdxdy = 0$.

Definition 4.1 (Subdifferential). Let $f$ be a proper lower-semi-continuous convex function on $R^m$, where $R^m$ is $m-$dimensional Euclidean space. The subdifferential of $f$ is defined for $y \in R^m$ by

$\partial f: = \left\{ y \in R^m and f(z) \ge f(x)+ \left \langle y, z-x \right \rangle, \forall z \in R^m \right\}$.

Definition 4.2 (Bregman distance). Assuming $J: X \to R$ is a convex function, where $X \in R^m$, $u, v \in X$, $p \in \partial J(v)$, then the Bregman distance between $u$ and $v$ is defined as $D^p_J(u, v)$:

From the definition of the subdifferential Bregman distance $D^p_J(u, v) \ge 0$, the Bregman distance does not satisfy symmetry, $D^p_J(u, v) \ne D^p_J(v, u)$. Assume $p \in \partial J(u) \cap \partial J(v)$, and then

Applying Bregman's regularization method, (4.1) can be rewritten as

According to the generalized form of Bregman distance, $D(u, u_{k-1}) = J(u)-J(u_{k-1})-\left \langle u-u_{k-1}, p_{k-1} \right \rangle$, where $p_{k-1} \in \partial J(u_{k-1})$, so (4.2) can be rewritten as

Omitting the constant term (the term about $u_{k-1}$) irrelevant to minimization yields

The Euler-Lagrange equation of (4.3) is

where $\partial _u G$ denotes the variation of $G$ with respect to $u$.

We are led to the relation:

Let $\lambda = \triangle t$, $p(u_k) = p(u_k(t))$, and

When $\triangle t \to 0$, we arrive at a continuous process

The initial conditions for (4.6) are $u \mid_{t = 0} = c_0$, $p \mid_{t = 0} = 0$, $\int_{\Omega} \partial _u G(c_0, f) = 0$.

The flow begins with $u(0) = c_0$ and it converges back to the image $f$ as $t \to \infty$, i.e., $\underset{t \to \infty}{\lim} u(t) = f$. By (4.6) the image $u(t)$ flows from the smoothest possible image $(u(0) = c_0)$ to the noisy image $f$. We are trying to use the flow to denoise the image.

Since the relationship between $p$ and $u$ is difficult to calculate in the case of non-linearity, (4.6) cannot be directly calculated, so we propose a relaxation method to approximate the direct flow.

The first-order optimality condition and $v_k = p_k/\lambda$ yield the following formula from (4.5):

Therefore,

and a system of equations about $u_k$(4.3) and $v_k$(4.7) can be expressed as

where $u_0 = c_0$, $v_0 = 0$, $k = 1, 2, 3, \cdots$, $\int_{\Omega} \partial _u G(c_0, f) = 0$.

It can also be rewritten as

After converging to a minimizer $u_k$ in the first formula of (4.8), it is easy to compute $\partial_u G(u_k, f)$ and update $v_k$. However, the amount of calculation in this process is quite high.

The standard way of solving $u_k$ and $v_k$ is first to evolve a steepest descent flow for $u_k$, having a fixed $v_k$, based on the Euler-Lagrange equations of (4.8):

The corresponding diffusion equation is

The update for $v_k$ in the second formula of (4.8) can be viewed as an iterative descent in $v_{k-1}$ for minimizing $G(u_k, f)$. This is an indirect minimization, which affects $u_k$ by its coupling with $v_k$. We write the solution for $v_k$ in the following manner:

$\tau$ is a time variable similar to $t$. This extends the definition of the sequence $v_k$ to a continuous formulation.

We propose two continuous flows $u(t)$ and $v(t)$ to approximate the sequences $u_k$ and $v_k$. Let $\tau = m t$, hence we obtain the following:

Replacing $v_k$ in (4.9) by $v(t)$ and $u_k$ in (4.10) by $u(t)$ yield the relaxed inverse scale space flow:

with $p \in \partial J(u)$ and initial conditions $u \mid_{t = 0}\ = \ c_0, v \mid_{t = 0}\ = \ 0$, $\int_{\Omega} \partial _u G(c_0, f) = \ 0$.

4.2.   Analysis of our model

4.2.1.   Inverse scale space flow of our model

Consider our model

In ^[17,20], the authors used the inverse scale space method to solve the denoising model containing only one fidelity term. For the two fidelity terms in our paper, we do the following processing to apply the inverse scale space method theory.

The greatest common divisor is the greatest common divisor shared by two or more integers, and the parameters of the fidelity term of the model proposed in this paper may be small, thus we generalize the concept of the greatest common divisor as follows. We extend the concept of the greatest common divisor from integers to decimals. Define $\lambda = d(\alpha_1, \alpha_2)$,

where $\lambda$ denotes the greatest common divisor of $\alpha_1$ and $\alpha_2$, $\mu_1$ and $\mu_2$ are integers, and $\alpha_1$ and $\alpha_2$ can be integers or decimals. Below are two examples of $\lambda$.

When $\alpha_1 = 1$, $\alpha_2 = 2$, then $\lambda = d(1, 2) = 1$. $\lambda$ is an integer.

When $\alpha_1 = 0.36$, $\alpha_2 = 0.45$, then $\lambda = d(0.36, 0.45) = 0.09$. $\lambda$ is a decimal.

In the model of this paper, the function $G(u, f)$ is composed of two terms and the coefficients of the fidelity term are two different decimals, whereas the traditional inverse scale space method only applies to one fidelity term, so we promote the concept of the greatest common divisor to combine $d(\alpha_1, \alpha_2)$ with the traditional inverse scale space method.

We will mark the last two terms of $E(u)$ as $g(u, f)$. Let $g(u, f) = \lambda G(u, f)$, and therefore,

Thus, it is possible to obtain:

Then, we get the following inverse scale space model from (4.6):

We also have

from (4.10).

We can also obtain the corresponding relaxed inverse scale space flow of the model proposed in this paper:

Its initial conditions are $u_0 = c_0$, $p_0 = 0$, $\int_{\Omega} \partial _u G(c_0, f) = 0$.

The procedure of solving (4.12) is summarized as follows.

4.2.2.   The convergence behavior

It is obvious that these equations $(\frac{\partial u}{\partial t} = 0, \frac{\partial v}{\partial t} = 0)$ are the equilibrium: $u = f, v = \frac{p(f)}{\lambda}$. We still try to show that for any $f \in BV$, the solutions converge to this steady state.

The following energy function is defined to analyze the convergence:

where $q \in \partial J(f), u = u(t), v = v(t)$.

Proposition 4.1. Let $u(0), v(0)$ be an initial value such that $e(0) < \infty$. Let $u(t), v(t)$ be the solution of (4.12) with $\lambda > 0, \alpha > 0$. Then the energy $e(t)$ decreases monotonically. Moreover, when $t \to \infty$, then $e(t)\to 0$.

Proof. We calculate the time derivative of the energy functional to obtain

In steady state, $u-f = 0$, this is equivalent to $\mu_1 \left(\frac{f}{u}-1\right)+\mu_2(f-u) = 0$. We replace $\mu_1 \left(\frac{f}{u}-1\right)+\mu_2(f-u)$ by $f-u$.

According to the definition of the Bregman distance, we know $D(u, f) = J(u)-J(f)-\left \langle u-f, \partial J(f) \right \rangle$ and $D(f, u) = J(f)-J(u)-\left \langle f-u, \partial J(u) \right \rangle$. Therefore, $D(u, f)+D(f, u) = \left \langle u-f, \partial J(u)-\partial J(f) \right \rangle$.

Finally, we get

which implies that the energy $e(t)$ decreases monotonically. Because of $e(t) \ge 0$ and $\frac{de(t)}{dt} \le 0$, we get $\underset{t \to \infty}{\lim}e(t) = 0$. □

The dual function of $u$ is $E^{*}(p): = \underset{u}{\sup} \left\{\left \langle u, p \right \rangle -E(u)\right\}$, and then $p \in \partial_u E(u)$ and $u \in \partial_p E^{*}(p)$ are equivalent. If we can compute the dual function $E^{*}$, we can get an explicit relationship for $u(p)$. Based on the above conditions, we can get the following estimator.

First, the derivatives of the fitted function and $u(t)$ in the time direction are calculated separately:

Because $u > 0$, $u$ is bounded, $\mu_1 > 0$, and $\mu_2 > 0$, we might as well set $\frac{\mu_1}{u}+\mu_2 = a > 0$, and then

where $G^{*} = \partial_{pp} E^{*}$ is the Hessian matrix of the dual function.

If $E^{*}$ is strictly convex, then there is a constant $b > 0$ that makes

where $\forall \varphi, \omega \in BV^{*}$. Therefore combine (4.13), (4.15), and (4.16), and we get

Using Gronwall's inequality, there is

if $t > s$. Thus, $t \to \infty$, $u \overset{L^2}{\rightarrow } f$.

In the estimate of $L^2$ above, $f$ can be either a clean image or an image containing noise. Assume $f$ is a clean image and $E(f) < \infty$, and we can get the Bregman distance error estimator

If $q(t)\in L^2 \cap \partial E(f)$, we can get estimators for $p(t)$ and $q(t)$,

Integrate the above formula:

Since $p(0) = 0, \parallel p(t)-q(t) \parallel ^2_{L^2}$, $D(u(s), f)$ are all non-negative and $D(f, u(s)) \ge D(f, u(t))$ $(s \le t)$, $D(f, u(t)) \le \frac{\parallel q(t) \parallel^2_{L^2}}{2at}$. Thus, $t \to \infty$, $u \overset{D}{\rightarrow }f$.

When $f$ is a noisy image, $w$ is a clean image, and $E(w) < \infty$, there may be $E(f) = \infty$. In this case, there are the following propositions:

Proposition 4.2. Based on the above conditions, if $\parallel f-u(t) \parallel_{L^2} \ > \ \sigma$ and $\parallel f-g \parallel_{L^2} \ \le \ \sigma$, and $\sigma$ is the standard deviation of the noise, then $D(g, u(t))$ is monotonically decreasing in time.

Proof.

When $\parallel f-u(t) \parallel_{L^2} \ > \ \parallel f-g \parallel_{L^2}$, the right-hand side of the inequality is negative, and $D(g, u(t))$ is monotonically decreasing in time. □

4.2.3.   Initial conditions

For the initialization of the inverse scale space, $u(0) = c_0$, we use the fact that

which means that $\int_{\Omega} \partial_u G(u(t), f) = 0$ is time invariant for all time.

As $t \to 0$, the expression vanishes, leading to

Combining $\partial_u G(u, f) = \mu_1 \left(1-\frac{f}{u}\right) +\mu_2 (u-f)$, we obtain

which means that

so

The solution of this Eq (4.17) for $c_0$ is: $c_{01} = f, c_{02} = -\frac{\mu_1}{\mu_2} < 0$.

Thus the initial condition of ISS flow is given: $c_0 = f$.

4.2.4.   Second order in time formulation

Our relaxed model can be written in a general form as

The first equation in (4.18) can be written in $u$ with a second-order derivative in the time domain.

Replacing $\partial_t v$ in the above formula with the second equation in (4.18), we can obtain

The above formula can be rewritten as

with the initial conditions $u\mid_{t = 0}\ = \ c_0, u_t \mid_{t = 0}\ = \ -\lambda \partial _u G(0, f)$.

As in ^[30,31] we approximate $p(u) = -\Delta u$, which leads to

We fix $u$ as a constant $u_0$ and $f$ as a constant $f_0$. Rewriting (4.19) in the frequency domain $\xi$ by taking the Fourier transform, we obtain the characteristic equation

with two solutions

Both roots exist for all frequencies if

and we know

Thus, we should require the parameters $\alpha$ and $\lambda$ meet

to ensure monotone evolution of $u$ toward $f$.

5.   Numerical results

In this section, some numerical experiments on medical images show the performance of our proposed model, expanding the application field of the model. The peak signal-to-noise ratio (PSNR) and the structural similarity (SSIM) are currently the most widely used methods for evaluating image quality. The larger the value of PSNR is, the closer the denoised image is to the original image, which means that the effect of image restoration is better.

First, we introduce the mean square error (MSE).

where $f_{ij}$ and $u_{ij}$ respectively represent the pixel value of the restored image and the original image at $(i, j)$, and $M$ and $N$ are the number of rows and columns of the image.

The specific calculation formula of the PSNR is as follows:

where $L$ denotes the maximum gray level of the image. For 8-bit grayscale images, $L$ is equal to $255$.

The expression for SSIM is:

where $\mu_f$ and $\mu_u$ are the means of $f$ and $u$, respectively; $\sigma_f^2$, $\sigma_u^2$ are the variances of $f$ and $u$, respectively, $\sigma_{fu}$ is the covariance of $f$ and $u$, and $c_1$ and $c_2$ are constants used to maintain stability.

5.1.   Analysis of medical image denoising performance with ISSM

5.1.1.   Different figure and noise level test

As shown in Figure 1, our numerical experiments use six real medical images, which are called "Adenocarcinoma", "Chest", "Brain", "CXR", "Knee", and "Foot".

First, a more thorough analysis is given. We have compared the denoising effect of different sizes of the same image at the same noise level $(25, 0.1)$ using the ISSM denoising method of this paper and the results are shown in Table 1. According to Table 1, we can find that as the image size increases, the PSNR value after denoising also increases. As the image size is larger, the number of pixels contained in the image is more and the observed image is also clearer. The result is that image processing takes longer. Based on the above considerations, we chose an image of size 512*512 for the experiment in this paper.

In the first experiment, three images with a size of $512*512$ were destroyed by different degrees of noise: $(k_0, k_1) = (10, 0), (25, 0), (25, 0.1), (25, 0.2), (50, 0)$. The PSNR values of the image contaminated by noise and the PSNR values after denoising are shown in Table 2. After a lot of experiments, in order to make the noise reduction effect better, we chose the parameters $dt = 0.03, \; \beta = 0.01, \alpha_1 = 0.01, \alpha_2 = 0.005$, and the number of iterations equal to $800$. It can be seen that the PSNR values of the three images have been improved after denoising by our model.

Figure 2 shows the PSNR values of the six medical CT images tested in Figure 1 at different noise levels, from which we can see that the denoising effect of these six images varies at the same noise level. Among them, under the noise level of $(10, 0)$ and $(25, 0)$, image "Knee'' has the highest PSNR value after denoising and the best recovery effect; under the noise levels of $(25, 0.1)$, the PSNR of images "Adenocarcinoma'', "Brain'', "Knee'', "Foot'', are basically about the same, but image "Brain'' is slightly higher; under the noise level of $(25, 0.2)$, the PSNR of image "Brain'' is the highest and the denoising effect is the best; under the noise level of $(50, 0)$, image "Adenocarcinoma" has the highest PSNR value, which is not distinctly high compared with other tested images.

We plot the PSNR and running time (CPU) of the tested images after denoising at different noise levels as bar charts (as shown in Figure 3), where 3(a) is the PSNR after denoising and 3(b) is the running time (CPU). The dark blue, red, yellow, light blue, green, and purple bars in the figure represent images "Adenocarcinoma", "Chest", "Brain", "CXR", "Knee", and "Foot", respectively. From 3(a), it can be found that the denoising ability of different images is different for different noise levels. From 3(b), it can be seen that the denoising runtime for the "Chest'' image is longer for all the different noise levels, and the runtime for the other five images is basically the same. We can also see from 3(b) that there is not much difference in the running time when the noise level is relatively small; however, when the denoising level is large, the difference in the running time is more obvious, and the running time of the X-ray (1(d)–1(f)) is more similar compared to that of the CT image (1(a)–1(c)).

5.1.2.   Different perspectives

In the second experiment, we illustrate from all aspects that our model has good performance. We take $(k_0, k_1) = (25, 0.1)$ to denoise the image in Figure 1. The choice of parameters remains $dt = 0.03, \beta = 0.01, \alpha_1 = 0.01, \alpha_2 = 0.005$.

Figure 4 shows that our proposed model can remove the mixture of additive and multiplicative noise to some extent.

In Figure 5, the histogram of the image shows the distribution of different gray levels in the total pixels. The horizontal axis represents the value of the pixel, and the vertical axis represents the number of times that the corresponding pixel value appears. It can be seen that the histogram of the image destroyed by noises is completely different from the original image, but the histogram of the denoised image is very close to the original image.

In Figure 6, we show the $250$th line of the tested images. The original image, noisy image, and denoised image are represented by blue, purple, and yellow solid lines, respectively. From these figures, we can see that when the original image is contaminated by noise, the one-dimensional image oscillates greatly. The yellow solid line is highly consistent with the blue, which represents that the denoising effect of our proposed model is very effective.

In Figure 7, we look at the same areas of the tested image and mark them with a red rectangle. The three-dimensional surface map of the original image is relatively flat, and that of the noisy image is rugged and accompanied by many sharp protrusions. It can be clearly seen that the three-dimensional surface map is highly consistent with the original image, which concludes that our proposed model has a good denoising effect.

5.1.3.   $\alpha_1$ and $\alpha_2$-dependent test

In the third experiment, we first selected several different sets of parameter values and conducted experiments, and the results are shown in the following Table 3. According to Table 3, we can find that different parameter values affect the denoising effect (PSNR value and SSIM value), and therefore, to further find the suitable parameter values, we utilize a grid search for them.

We examine the impact of changing the values of $\alpha_1 \in [0.001, 0.05]$ and $\alpha_2 \in [0.001, 0.05]$ on the proposed variational method. Our method is executed with the parameters $dt = 0.03, \beta = 0.01$, the number of iterations $N = 700$, and a noise level equal to $(25, 0.1)$. As clearly seen in Figure 8(a), we set $\alpha_2 = 0.005$, according to the PSNR value, and it can be seen that when $\alpha_1 \in [0.001, 0.015]$, the effect of removing this type of noise is excellent. From 8(b), we set $\alpha_1 = 0.003$, and it can be seen that when $\alpha_2 \in [0.001, 0.01]$, the effect of removing this type of noise is excellent. In Figure 8(c), we show a three-dimensional plot of the changes in the PSNR values corresponding to changes in the regularization parameters $\alpha_1$ and $\alpha_2$. The darker the red, the higher the PSNR value, and the red areas are those with higher PSNR values, which we find are consistent with the conclusions reflected in 8(a) and 8(b). For more details, see Figure 8.

From Figure 8, we can further see that although the size of parameters $\alpha_1$ and $\alpha_2$ do not have a significant impact on the PSNR value after image denoising as a whole, there are individual parameter values that will make the PSNR value drop dramatically, or even become negative. In order to make the model's denoising effect more robust and accurate, we need to discuss the parameters.

5.1.4.   $N$-dependence test

In the fourth experiment, in order to decide the value of the $N$, we study the relation between the PSNR and the number of the iterations $N$. The Figure 9 shows the relation between the PSNR and the number of iterations $N$ on "Adenocarcinoma", "Brain", and "Knee" images corrupted by the combined additive and multiplicative noise with noise intensity $(k_0, k_1) = (25, 0.1)$. Figure 9 shows that the curves of the relationship between the PSNR and number of iterations are slightly different from other images. In 9(a), we can see that the curve first shows an upward trend and then decreases, but finally reaches a steady state; in 9(b) and 9(c), we can see that the curve shows an upward trend and finally reaches a steady state. Therefore, for images "Adenocarcinoma" and "Brain", the optimal number of iterations $N$ for removing the mixed noise is in the interval $(1000, 1250)$, and for image "Knee", the optimal number of iterations $N$ for removing the mixed noise is around $1500$.

5.2.   ISSM and GDE denoising performance analysis

5.2.1.   Different algorithms

In the following, we still test the six images of Figure 1. A noise of $(k_0, k_1) = (25, 0.1)$ is added to the images and denoised with gradient descent (GDE) and the inverse scale space algorithm (ISSM), which shows the PSNR and running time (CPU) after denoising as shown in Table 4 where higher PSNR values and smaller CPUs are in bold. In this experiment, the parameter values are tuned as $dt = 0.03, \beta = 0.01, (\alpha_1, \alpha_2) = (0.01, 0.005)$, and the number of iterations is equal to $500$. In Table 4, we can find that the PSNR and SSIM values of the ISSM method are slightly higher than those of the GDE, and for images 1(a), 1(c), 1(f), the running time of the GDE is slightly longer than that of the ISSM, whereas for images 1(b), 1(d), 1(e), the running time of the GDE is about $2$ seconds longer than that of the ISSM. In general the ISSM is a little bit shorter in terms of running time than the GDE. Therefore the method in this paper is better than the GDE in terms of PSNR, SSIM, and runtime.

We used the "Adenocarcinoma" image to compare the ISSM method with the GDE method in detail. As shown in Figure 10, 10(a) shows that the GDE denoising method yields a greater number of low-pixel values than the ISSM method. 10(b) shows that the ISSM method and the GDE method are essentially identical at line $370$th of the image after denoising. However, in certain instances, the ISSM method outperforms the GDE method. As shown in 10(c), the 3D surface maps of the images denoised by the ISSM method and the GDE method are basically the same. Thus, Figure 10 as a whole reflects that the inverse scale space method (ISSM) is superior to the gradient descent method (GDE).

5.2.2.   Testing of other types of medical images

In this subsection of the experiment, we tested brain tumour MRI and fetal head ultrasound images, as shown in Figure 11. The original image, the noisy image, the ISSM denoised image, and the GDE denoised image are shown from left to right in Figure 12, where the noise level of the noisy image is $(25, 0.1)$. As can be seen from the figure, visually we can see very little difference with the naked eye, so we analyzed and compared them again from a quantitative point of view, as shown in Table 5.

In Table 5, we can clearly see that the PSNR and SSIM values obtained by our method in this paper are higher than those of GDE, which indicates that our method can effectively denoise another two types of images, MRI and ultrasound, besides X-ray and CT; and according to the comparison of the last two columns of the running time, it can be seen that our method has a slight improvement in comparison with GDE, which further illustrates the advantages of our method for the adaptability of medical image denoising.

5.2.3.   Other noise type denoising

In this section of the experiment, we add salt and pepper noise with a noise level of $0.02$ to the image in Figure 11 and denoise it with the method of this paper and GDE, and the results are shown in Figure 13 and Table 6. From the figure, we can find that both this paper's method and GDE can effectively remove part of the salt and pepper noise, and ISSM is better than the GDE denoising effect, but both ISSM and the GDE method cannot effectively retain the texture details. From Table 6, it can be seen that the PSNR value after denoising is not high, and there is not much difference between the two methods.

6.   Conclusions

In this paper, we combine the data fidelity terms from the ROF model and the I-divergence model to restore medical images with combined noises. We have proved the existence and uniqueness of the minimum solution of the model. The inverse scale space method has been used to solve the non-linear PDE resulting from minimizing the proposed functional. Experimental results demonstrate that the proposed model can remove combined noise on medical images from different perspectives. The inverse scale space method on this model improves PSNR better and has faster runtime than traditional methods on medical images. Parameter sensitivity is also discussed.

Medical imaging is often affected by noise, resulting in image quality degradation, while the model and method proposed in this paper can effectively denoise medical images, thus enabling doctors to more accurately identify lesions and reduce the possibility of misdiagnosis and omission of diagnosis, thus improving diagnostic accuracy and enhancing lesion detection. It can also help patients to reduce the number of repeated examinations and enhance the patient experience. Learning-based methods are currently more popular research. Combining the modeling approach of this paper with deep learning as well as further exploring real medical images with noise are directions for future work.

Author contributions

Conceptualization, Chenwei Li and Donghong Zhao; methodology, Chenwei Li and Donghong Zhao; software, Chenwei Li; formal analysis, Chenwei Li and Donghong Zhao; writing—original draft preparation, Chenwei Li; writing—review and editing, Chenwei Li and Donghong Zhao; visualization, Chenwei Li. All authors have read and approved the final version of the manuscript for publication.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors thank the reviewers for their helpful comments and valuable suggestions. This research is funded by the National Natural Science Foundation of China (grant number 12371481) and Youth Teaching Talents Training Program of USTB (2016JXGGRC-002).

Conflict of interest

All authors declare no conflicts of interest.