Speckle Noise Removal by Energy Models with New Regularization Setting

In this paper, we introduce two novel total variation models to deal with speckle noise in ultrasound image in order to retain the fine details more effectively and to improve the speed of energy diffusion during the process. Firstly, two new convex functions are introduced as regularization term in the adaptive total variation model, and then, the diffusion performances of Hypersurface Total Variation (HYPTV) model and Logarithmic Total Variation (LOGTV) model are analyzed mathematically through the physical characteristics of local coordinates. We have shown that the larger positive parameter in the model is set, the greater energy diffusion speed appears to be, but it will cause the image to be too smooth that required adequate attention. Numerical experimental results show that our proposed LOGTV model for speckle noise removal is superior to traditional models, not only in visual effect but also in quantitative measures.


Introduction
With the development of digital image technology, a large number of digital images are transmitted and compressed through various channels. However, image corruption usually is unavoidable during transmission and storage, and the resultant noise quite often seriously affects the visual effect of the image. Clearly, high-quality images are desirable in many areas, such as in medical imaging and pattern recognition. Thus, image denoising is a critical step in image processing and computer version, which plays an important role in various applied areas, especially in medical imaging, video processing, and remote sensing. On the other hand, it is also an important preprocessing process for other image processing that relies on subsequent processing.
In 1992, Rudin et al. [10] proposed a denoising model based on the total variation: where Ð Ω jDujdx = sup f Ð Ω u div ðφÞ | φ ∈ C 1 c ðΩ, ℝ n Þ, kφk ∞ ≤ 1g represents the TV regularization term, u : Ω ⟶ ℝ is an original image that is without noise, u 0 = u + n is noisy image, and n represents the Gaussian random noise with mean zero and standard deviation σ. λ 1 > 0 represents the regularization parameter which can be used to balance the fidelity terms and regularized terms in TV model, and jDuj represents the L 1 norm of the image gradient.
It is well known that medical ultrasonic images may have many noises of speckle, which will bring a significant problem in terms of the quality of ultrasonic images and cover up the lesions of certain important tissues. Further, it brings great difficulties to the actuate diagnosis and identification of certain specific diseases and can create the potential risk of missed diagnosis and misdiagnosis. Thus, it is very desirable to eliminate the speckle noise in ultrasound image and simultaneously retain the important features in practices. As mentioned in article [3], the speckle noise in medical ultrasonic images can be formulated in the following form: where f is a noisy image, and n represents the Gaussian random noise with zeros mean and standard deviation σ.
In this paper, we focus on the image denoising form by using variation method, where the noisy image is ultrasound speckle noise. Based on the model (2) and the characteristics of the Gaussian distribution, Krissian et al. in the article [16] derived a date fidelity term: Within the variational framework, the data fidelity function is derived from the degradation model (2). One of the technical approaches to solve the variational model is the regularization technique, which minimizes the cost function to obtain stable and accurate solutions. In general, the image denoising variation method is to consider: where TVðuÞ is a regularization term that represents a prior information about the object to be restored, and F ðu, f Þ is a fidelity term to ensure that the restoration u is not far from the original observation f . λ > 0 represents the regularization parameter which can balance the fidelity term and the regularization term.
In [11], motivated by the classical ROF model [10], the authors proposed a convex variational model for removing the speckle noise in ultrasound image. The convex variational model involves the TV regularization term and convex fidelity term (see Equation (5)): where Ð Ω jDujdx represents the TV regularization term, and Du represents the directional gradient of u. The existence and uniqueness of the solution of model (5) is proved in [11]. In this paper, we call the model in (5) the "JIN's model." In [17], the authors proposed a well-balanced speckle noise reduction (WBSN) model that can detect edges.
Although TV regularization is effective for image denoising, it also leads to some staircase effects that is undesirable. In order to solve this problem, many methods based on improved TV regularization are proposed, such as highorder TV regularization [18][19][20], several hybrid TV regularization [21], the improved infimal convolution [22,23], nonlocal TV model [24,25], fractional order TV model [26], and anisotropic TV model [27,28]. Fractional theory [29,30], wavelet [31], and statistical information [32][33][34] are also employed to deal with intensity inhomogeneity or noise. Although these denoising methods can reduce the staircase effects in the restoration of additive noise images, however, there are more staircase effects that appeared in the restoration of speckle noise images. In this paper, we introduce HYPTV model and LOGTV model to reduce speckle noise and staircase effects in ultrasound images in an effective way.
In numerical algorithms, most of the energy function minimization problems can be transformed to an Euler-Lagrange equation and then be solved by using the finite difference method. However, the choice of adequate regularization terms is critical in terms of solution accuracy. Moreover, when solving the Euler-Lagrange equation, the energy diffusion form of the noise image in different regions is required to handle differently based on the physical characteristics of local coordinates, since the diffusion velocity of different parameters to the energy function is different in the two models. The new proposed HYPTV model and the LOGTV model, as shown in this paper, not only can preserve the edges of the restored images well when restoring the ultrasonic image with speckle noise but also better reduce the staircase effect generated during the recovery process.
The rest of this paper is as follows: in Section 2, we review some background knowledge. In Section 3, we propose two new models based on variation; meanwhile, we not only analyze diffusion performance of the proposed models but also give the corresponding numerical algorithms. Section 4 shows five different experiments and results. The paper ends with concluding remarks in Section 5.

Preliminary
is extremized only if h satisfies the partial differential equation:

Journal of Function Spaces
For model (5), the corresponding Euler-Lagrange equation is as follow: where ∇ and div, respectively, represent gradient operators and divergence operators. Using gradient descent method, we can get the model as follows: where n ! is the unit out normal vector of ∂Ω.
Without losing generality, in the following, we consider the grayscale images as M × N matrices.

2.
Let u ∈ U = C 2 2 ðΩ, RÞ, g = ðg 1 , g 2 Þ ∈G = C 2 2 ðΩ, R 2 Þ, the gradient operators on the space U and the divergence operators on the space G are defined as: div : are the first-order forward and backward discrete derivation operators in the x-direction and y-direction, respectively, which are defined as: Applying the above gradient operators and divergence operators to model (9), we can obtain the equivalent minimization problem.
The following facts are easily checked: Theorem 4. If functions ϕ 1 and ϕ 2 are convex and have the same domain definition, then ϕ = ϕ 1 + ϕ 2 is also convex.
then ϕ is convex.

The Condition of TV Regularization Term.
Although TV regularization is very effective in image restoration, it usually generates some staircase effects. Thus, it is suggested in literature to use general variational methods, i.e., to consider: where φðsÞ represents a potential function, and the case φðsÞ = s leads to the total variation regularization term. In the literature [35], the author Costanzino chooses φðsÞ = s 2 that leads to the well-known harmonic model. In practice, we prefer good smoothing in some domain where the intensity of variations is relatively weak. This can be achieved by requiring a function φðsÞ to satisfy the following conditions: Near the edge of the image, the intensity of variations is strong. If we would like to preserve the edge, then the function φðsÞ should satisfy the following conditions: With the conditions (16) and (17), the function φðsÞ is convex and nondecreasing function, such as: There functions are convex and nondecreasing on s ≥ 0, as shown in Figure 1.
In this paper, we will use two new functions φ 1 ðsÞ = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 + α 1 s 2 p and φ 2 ðsÞ = s ln ðα 2 + sÞ, which appear to be quite effective for image processing, in particular, for ultrasound image denoising. Obviously, the functions φ 1 ðsÞ and φ 2 ðsÞ are convex and nondecreasing.

The Proposed Restoration Model
In this section, we propose adaptive total variation model for image restoration. We use the finite difference method to solve the Euler-Lagrange equation directly, and then find the minimum value of the energy function.

The Adaptive Total Variation Model.
Apply the selected function to model (15), we propose adaptive total variation model, where Ð Ω φ i ð|∇u | Þdx is a regularization term, and i = 1, 2; Ð Ω ððf − uÞ 2 /uÞdx is a fidelity term; λ represents the regularization parameter which can balance fidelity term and regularization term.
Firstly, the energy functional EðuÞ is convex, which guarantees the existence of the minimal solution of the model (19). Theorem 6. The energy functional EðuÞ is convex. That is to say, for any u 1 , u 2 ∈ Ω, and t ∈ ½0, 1, we have: where Proof. Firstly, the function φ i ðsÞ is convex; according to Definition 3, for any u 1 , u 2 ∈ Ω, and t ∈ ½0, 1, we have Therefore, This proof is established. Secondly, the uniqueness of the minimum solution of the model (19) can also be proved. Proof. According to Theorem 6, we have If u 1 ≠ u 2 , then the above assumption gives a contradiction that u 1 is not a minimize solution.

Diffusion Performance.
In this subsection, we mainly analyze the diffusion performance and diffusion speed of the energy function of HYPTV model and LOGTV model.

Diffusion
Performance of HYPTV Model. Firstly, we use the finite difference method to solve the HYPTV model, associated with the potential function φ 1 ðsÞ = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 + α 1 s 2 p . From Definition 2, we can obtain the corresponding Euler-Lagrange equation HYPTV model that as follows: Using gradient descent method, Equation (25) can be transformed to: where n ! is the unit out normal vector of ∂Ω.
In order to analyze the diffusion performance, local image coordinate system ξ-η is established. As shown in Figure 2, the η-axis represents the direction parallel to the image gradient at the pixel level, and the ξ-axis is the corresponding vertical direction.
According to Figure 2, we can know: Journal of Function Spaces So, Equation (26) can be rewritten as: where u ηη = u 2 x u xx + 2u x u y u xy + u 2 y u yy ∇u j j 2 : The ψ 1 1 ð|∇u | Þ and ψ 1 2 ð|∇u | Þ are control functions of the diffusion along the ξ-direction and η-direction, respectively. Now, we consider the diffusion of image restoration. Some test images are shown in Figure 3.  (19) is isotropic. In other words, the energy diffusion rate along direction ξ and direction η is very close in the process of image restoration in the smooth region. And the rate of energy diffusion is obviously positively correlated with the parameter α 1 .
(2) Sharp area. When |∇u | ⟶∞, we obtain lim |∇u|→∞ ðψ 1 2 ð|∇ u | Þ/ψ 1 1 ð|∇u|ÞÞ = 0. This shows that the diffusion form of the energy Equation (19) is anisotropic. In other words, the energy diffusion rate in ξ-direction in Equation (28)  One can see that the larger the parameter α 1 is set, the smaller the limit becomes. And the rate of energy diffusion is obviously positively correlated with the parameter α 1 .

Diffusion
Performance of LOGTV Model. Secondly, we use the finite difference method to solve the LOGTV model, which the potential function is φ 2 ðsÞ = s ln ðα 2 + sÞ.
Using gradient descent method, Equation (30) can be transformed to: where n ! is the unit out normal vector of ∂Ω.
Hence, Equation (31) can be rewritten as: where ψ 2 1 ð|∇u | Þ and ψ 2 2 ð|∇u | Þ are control functions of the diffusion along the ξ-direction and η-direction, respectively. Now, we consider the diffusion of image restoration.  (19) is anisotropic. However, the gradient of noise image is relatively large, so in the smooth region, whatever the value of α 2 , it has little effect on the model.  (19) is anisotropic. In other words, the energy diffusion rate in ξ-direction in Equation (28) is much larger than that in the η-direction in the sharp region. But the gradient |∇u | does not exceed 255, so lim |∇u|→255 ðψ 2 ð|∇u | Þ/ ψ 1 ð|∇u|ÞÞ = ð510α 2 + 255 2 Þ/ð225α 2 + 255 2 + ðα 2 + 255Þ 2 * ln ðα 2 + 255ÞÞ. One can see that the larger the parameter α 2 is set, the smaller the limit becomes. And the rate of energy diffusion is clearly positively correlated with the parameter α 2 .
3.3. Numerical Implementation. We will describe the corresponding numerical algorithm in this section. Firstly, the HYPTV model can be solved by discretization as follows: where T 1 ðj∇u k jÞ = α 1 / ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 + α 1 j∇u k j 2 q , and Δt represents time step. Furthermore, the iterative formula can approximate as: for i = 1, ⋯, M; j = 1, ⋯, N, and M × N represent the size of the image. Here: where m½a, b = ððsign a + sign bÞ/2Þ ⋅ min ð½jaj, jbjÞ, and δ > 0 is a positive parameter that is close to zero. With boundary conditions:

10
Journal of Function Spaces Now, note Equation (25), the two sides are multiplied by ðf − uÞu/f + u, and then, the integral on the domain Ω can be obtained: Because the Gaussian noise n have mean 0 and variance σ 2 , we can obtain: In the process of iteration, we always use the previous solution to calculate the next solution. The optimization algorithm for HYPTV model is given in the following (Algorithm 1).
Secondly, the LOGTV model can be solve by discretization as follows: where T 2 ðj∇u k jÞ=ððln ðα 2 + j∇u k jÞÞ/j∇u k jÞ+ ð1/ðα 2 + j∇u k jÞÞ. Furthermore, the iterative formula can approximate as: Similar to the HYPTV model, the optimization algorithm for LOGTV model is given in the following (Algorithm 2).

Experimental Results
In this section, we present numerical results to demonstrate the effectiveness of the HYPTV and LOGTV model in image restoration. Firstly, to evaluate the quality of restored images, we use the peak signal-to-noise ratio (PSNR) value and the structure similarity (SSIM) index, which are defined as follows: where u ∈ ℝ m×n is the clean image, and u ∈ ℝ m×n is the restored image. μ a is the average of a, σ a is the standard deviation of a, and c 1 and c 2 are some constants for stability. Secondly, we calculated the noise deviation reduction (NDR) at each iteration as a convergence condition;  Journal of Function Spaces where u k represents the results of the kth iterations, respectively. Finally, the stopping condition (NDR) for the HYPTV and LOGTV models is as follows: In the numerical experiment, we will use the noise image as the initial value, that is, f = u 0 . Moreover, the gray values of all original images are in range [0, 255].

Denoising Effect of Different Parameters of HYPTV Model and LOGTV Model.
In this example, we use different parameter α i ði = 1, 2Þ values in the algorithm to test the effect of HYPTV and LOGTV models. The test images is shown in Figure 3(a), and the noise levels are σ = 2 and σ = 3. Figure 4 shows the different PSNR values and iteration numbers when different β values are used in the HYPTV model algorithm. Figure 5 shows the different PSNR values and iteration numbers when different β values are used in the LOGTV model algorithm. Firstly, we can also see that PSNR is the largest in α i ði = 1, 2Þ = 15. Secondly, with the increase of parameters, the speed of image restoration is faster. Based on the above analysis, when α i ði = 1, 2Þ = 15, the denoising performance of HPYTV and LOGTV models is close to the best. Therefore, in the following experiment, we choose α i ði = 1, 2Þ = 15 in HPYTV and LOGTV models.

Denoising Effect of the HYPTV and LOGTV Model.
In this subsection, we use the algorithm to test the effect of the HYPTV and LOGTV model on image denoising. In Figures 6 and 7, we show the original images, the noise images, and the restored images by HYPTV and LOGTV models.   13 Journal of Function Spaces images can be got by using the HYPTV model and the LOGTV model. "Noise image PSNR" is the peak signalto-noise ratio of noisy images and original images. "Denoising image PSNR" is the peak signal-to-noise ratio of restored images and original images. "Iter" is the number of iterations of the algorithm. From the results, it is obvious that the HYPTV model and the LOGTV model are fairly effective in reducing the speckle noise and edge-preserving.

Reduction of Staircase.
In this subsection, we test the reduction of the staircase effect in HYPTV model and LOGTV model by image in Figure 8. Figure 8(a) shows original images ("house"). Figure 8(b) shows detailed images of Figure 8(a), respectively. Figure 8(c) shows noise images which noise standard deviation σ = 3. Figure 8(d) shows detailed images of Figure 8(c). The original images contain a lot of details information, such as textures edges and in homogeneous regions. Figure 9 displays the restoration results of the noisy "house" image, respectively. Figures 9(a)-9(e) was restored ROF model [10], ATV model [36], JIN's model [11]  14 Journal of Function Spaces are shown in Figure 3("house," "peppers," "boat," "pirate," and "bird"), with two sizes of 512 × 512 and three sizes of 256 × 256. Figure 10 shows some noise images with different standard deviation. Figures 11-15 display the restoration results for images ("peppers," "boat," "house," "pirate," and "bird") through ROF model, ATV model, JIN's model, HYPTV model, and LOGTV model. The noise versions of "peppers", "boat" and "house", and "pirate" and "bird" are obtained by model (2) with standard deviations 2, 3, and 4, respectively. In addition, the detailed images of the restored images are also displayed.  Table 2, the noise standard deviation σ = 2. We can observe that although the four models can effectively remove the noises while preserving the edges and details, the restored images by HYPTV model and LOGTV model have better visual effect with less staircase effects than by the ROF model, ATV model, and JIN's model. The noise standard deviation σ = 3, 4. The visual effect of restored images by the ROF model and ATV model is particularly poor, but JIN's model, HYPTV model, and LOGTV model can effectively remove the noises. Finally, Table 2 shows that LOGTV model has higher PSNR and SSIM values than other four models. This means that our proposed LOGTV model is available in reducing the speckle noise in some images.

Denoising
Results of Real Ultrasound Images. In this subsection, we test some real ultrasound images. Figure 16 shows the experimental results of real ultrasound images by applying JIN's model, HYPTV model, and LOGTV model. Table 3 shows the different iteration for the different test images by using the JIN's model, HYPTV model, and LOGTV model. We find that LOGTV model is much effective than JIN's model and HYPTV model in obtaining the satisfactory restored images.

Concluding Remarks
In this paper, we propose a new speckle noise restoration model based on adaptive TV method. Two new convex func-tions are introduced as the TV regularization term. By analyzing the diffusion performance of the proposed two models, one can see that the LOGTV model has faster diffusion speed than the HYPTV model. Moreover, we introduced