A mesh-free algorithm for ROF model

The total variation (TV) denoising method is a PDE-based technique that preserves the edges well but has undesirable staircase effect in some cases, namely, the translation of smooth regions (ramps) into piecewise constant regions (stairs). This paper introduces a novel mesh-free approach using TV (ROF model) regularization and radial basis function (RBF) for the numerical approximation of TV-based model to remove the additive noise from the measurements. This approach is structured on local collocation and multiquadric radial basis function. These features enable this strategy not only to eliminate noise from images and preserve the edges but also has the advantage to minimize the staircase effect substantially from real and artificial images which cause the image to look blocky. Experimental results demonstrate that the proposed mesh-free approach is robust and performs well in visual improvement as well as peak signal-to-noise ratio compared with the recent partial differential equation (PDE)-based traditional methods.


Introduction
Image denoising is an inverse problem and is a very active area in the fields of image and signal processing that has been studied for the last three decades. Although there exist different types of noise, here we consider only models for removing additive, zero-mean Gaussian noise. This can be modeled as where z is the given noisy image containing the unknown additive noise (Gaussian noise) η and u is the known actual image, all of which are defined on a domain ∈ R 2 . In literature, there are many effective numerical techniques have been utilized to tackle such models connected with image denoising having additive noise, for instance, in [1,4,13,37,44,45].
The TV filtering [38] has proved to be one of the most successful tool in image processing for the solution of variational based partial differential equation (PDE) restoration problems. In this method, it is supposed that the images are defined on a continuous domain, which results in continuous functional. This functional then leads to a Euler-Lagrange equation. The resulting PDEs are then *Correspondence: d2014017@hhu.edu.cn; chenwen@hhu.edu.cn State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, College of Mechanics and Materials, Hohai University, Nanjing, Jiangsu 210098, People's Republic of China discretized by existing classical numerical methods on a regular grid for smooth solutions. For more details about the TV filtering, see [22,38,39]. The first TV-based model for image restoration having additive noise was proposed by Rudin et al. (ROF) [38]. This model yields very satisfactory results for removing image noise while preserving edges, see [7,37]. However, it also processes some unfavorable properties like staircase effect, loss of image contrast and in time computation due to its nonlinearity and non-differentiability [5,32,33,38]. In [38], the authors proposed an artificial time marching method to the associated Euler-Lagrange equation. This strategy is slow due to its strict stability constraints in the time steps. Also, the artificial time marching method computes the approximate solution, not the exact solution. Recently different procedures have been used to overcome this difficulty and hence some good results have been obtained, for instance, see [14, 18, 29-31, 43, 46]. But still, there is space for improvement. So, in this work, we adapt the mesh-free BRF collocation method to reduce these issues.
During the past decade, RBFs have been observed to be active techniques for the interpolation and approximation of multi-variable smooth functions on scattered data sets [3,6,12]. More recently, an increasing attention has been given to the development of mesh-free methods using RBFs for the numerical solution of PDEs. Most PDEs results have concerned steady state problems with smooth solutions. Recently there has been a growing interest in applying RBF methods to time-dependent PDE problems, again to problems with sufficiently smooth solutions. The RBF techniques have more points of interest and have exhibited superior accuracy as compared with traditional numerical strategies, for example, finite difference method (FDM) [20,23,47], finite element method (FEM) [23], finite volume method (FVM) [21,25], and pseudospectral method [24]. Interested readers can refer to [2,8,9,16,17,[26][27][28]40] for more details about the RBF collocation methods.
Global RBF collocation technique is also easy to implement, gives good accuracy and converges exponentially for solving the PDEs. Although, in this strategy, the interpolation matrix is fully populated and ill-conditioned, and thus sensitive to shape parameter. Thus, it is computationally extremely expensive to apply global collocation method to large scale problems. So in literature, there are many domain type collocation techniques, for example, Kansa technique [19,20] etc, is to settle these issues.
The main advantages of the RBFs for interpolating multidimensional scattered data are discussed in [19,20]. In recent decades, meshless methods have been proved to treat scientific and engineering problems efficiently. The mesh-free method based on the collocation method has been dominated and very efficient. Over the last several decades RBFs have been found to be widely successful for the interpolation of scattered data. RBF methods are not tied to a grid and in turn, belong to a category of methods called mesh-free methods. They apply only a cloud of points without any information about nodal connections. It is (conditionally) positive definite [3,36,42], rotationally and translationally invariant. The RBF approximation is an incredibly powerful tool for representing smooth functions in non-trivial geometries since the method is mesh-free and can be spectrally accurate [10]. RBFs interpolations have been used to remove the Gibbs oscillations from the given arbitrary data points [41] and very useful results have been obtained.
Motivated by the applications of TV-regularization in image restoration and RBF collocation methods for the solution of PDEs, we propose a new mesh-free strategy, with some modifications, of the TV (ROF) filter by RBF approximation to accomplish a new algorithm to solve the associated PDE with minimization of ROF model. This strategy is entirely mesh-less and is not only helpful to restore the image efficiently and resolve the edges due to is discontinuous jumps but also to eliminate the staircase effect and preserve the textures during the restoration process. The numerical treatment in this approach is also easy to implement and faster because of its mesh-free properties as compared to the traditional mesh-based numerical methods.
The rest of the paper is organized as follows. In Section 2, some details are provided related to the applications of TV-regularization and its detail use in ROF model for image restoration. This section also contains the shortcoming in ROF model. This section also includes the details of RBFs and its applications in solving PDEs and comparison with traditional methods. Two mesh-based methods, i.e., implicit and Augmented Lagrangian methods utilized for the solution of ROF model are presented in Section 3. This section also contains proposed method, i.e., BRF collocation method (Kansa method) for the solution of the associated PDE with ROF model. Section 4, describes experimental results and discussion, to compare the three methods for ROF model regarding CPU times, the number of iterations, and quality (peak signal to noise ratio (PSNR)) of the restored images. This section also includes the shape parameter analysis on image restoration and comparison of the proposed method with an other recent method. Section 5, shows the tabulated discussions about the sensitivity of parameters of the proposed method. The conclusion is provided in Section 6. And finally, the details for derivatives for our proposed method are given in an Appendix.

Rudin-Osher-Fatemi (ROF) model
The TV regularization is a process in digital image processing for the noise removal and also an important tool in the inverse problem and numerical [15]. It is the only regularizer which is used to preserve the edges and for removing the highly noisy frequency components from the image. It is also a convex one. So the total variation (TV) for an image u : → R 2 is defined as Rudin et. al (ROF) proposed the first model for image restoration from given noisy image having additive noise using TV regularization in [38]. This model achieved some useful restoration results.
The minimization approach for the model (1) using TV method is given as where the first term is the total variation of u, and the second term is data fitting term, respectively, and λ is the parameter. The fitting parameter λ is used for balancing the denoising and smoothing of the denoised image which usually depends on upon the noise level. The corresponding Euler Lagrange equation of the ROF model (2) is given as follows; with ∂u ∂n = 0 on the boundary of = ∂ . The time marching restoration PDE from (3) is given as follows, for the given u(x, y, 0), and also ∂u ∂n = 0 on ∂ . ROF model has the following disadvantages.
• This model yields staircase effect, in restoring the smooth images in applications where edges are not the main features. • This model also generates to the loss of image contrast during the restoration process. • This model also contains the difficulty with the non-differentiability term in the total variation norm.

Radial basis functions approximation
RBFs are mostly multivariate functions, and their values depend only on the distance from the origin, so that φ(x) = φ(r) ∈ R, x ∈ R n , r ∈ R; or alternatively on the distance from a point of a given set {x j } such that is called the radial function. Some commonly used, globally-supported RBFs are shown in Table 1. RBF interpolation of a continuous multivariate function, can be approximated by a linear combination of RBFs, namely, where γ j are unknown coefficients which must be determined. Using the collocation method, one may write: The above linear system of equations can be expressed in the following N × N linear system matrix form

and the RBF interpolation matrix is given by
However, some RBFs are conditionally positive definite functions as listed in Table 1, such as MQ, IMQ, GA, and TPS.
Hence, polynomials are augmented to Eq. (5) to guarantee that the resultant interpolation matrix is invertible. Such a formulation is expressed as follows with constraints Table 1 [k] denotes the nearest integers less than or equal to k, and N the natural number, c a positive constant which is known as the shape parameter, and CPD denotes the m-order conditionally positive definite functions [3,11] Name of RBF Definition CPD order (m) [k/2]+1 Thin plate splines (TPS) where m represents the polynomial space in which the total degree of all polynomials is then m in N variables [36], Then, Eqs. (6) and (7) yields matrix system of (M + N) Moreover, details of positive definite (PD) and conditionality positive definite (CPD) RBFs are discussed in [3,36] and listed in Table 1. For RBFs containing the shape parameter c, such as as in Table 1, small shape parameters produces more accurate results, but also associated with poorly conditioned interpolation matrix [3,36].

Augmented Lagrangian method (M2)
Chunlin Wu [43] first proposed this numerical method in 2009, for the solution of ROF model (2) and hence many good results have been obtained. So for this, he introduced a new variable p for ∇u and then separate the calculation of the non-differentiability term and the fidelity term. So the minimization problem (2) can be written as which is constrained optimization problem. So, the above problem (8) can be solved by Augmented Lagrangian method as where L rof (u, p, χ) is called the augmented Lagrangian functional, χ = (χ1, χ2) t is the Lagrange multiplier, and r is the positive constant. The system of optimality condition is written as follows The augmented Lagrangian process is used to solve (9), which is stated in the given Algorithm 1. To solve the Algorithm 1 :Algorithm for Augmented Lagrangian method for the ROF model where L rof (u, p, χ k ) is defined is functional (9). minimization problem (13), we split (13) into two subproblems.
For fixed p value, (14) and for fixed u value, For the u-subproblem, the optimality condition is given by the following linear equation which allows us to use the Fast Fourier transform F(u) and hence we get The p-subproblem can be formulated as and hence obtained the following closed solution where w = r∇u − χ k . The given Algorithm 2 shows the solution of sub-problem (13). For further details, we refer Algorithm 2 :Algorithm for Augmented Lagrangian method for the ROF model-solution of sub-problem (13).

Proposed method (M3)
In this subsection, our aim is to introduce a new algorithm, using both RBF interpolation and total variation norm to solve ROF model (4) and to reduce the difficulty associated with the total variation norm in ROF method (4). This methodology obviously uses advantages of both BRF interpolation and total variation norm which leads to the good restoration performance regarding restoration, eliminating the staircase effect, sharp edges and textures. Hence using this method, consistent improvement in PSNR values is obtained.
So, for any radial basis function the following equations satisfied, φ(r) = r 2 in R 2 , i.e r = (x, y). For {xc j } Nc j=1 , given Nc centers, defining of radial basis function without polynomial term one may write: ρ j coefficients in RBF is determined via enforcing the interpolation condition a set of points that usually coincides with Nc centers. The interpolation condition at Nc centers results in a Nc × Nc linear system Aρ = z which must be solved for expansion coefficients of ρ.
Nc × 1 matrices. The matrix A is called interpolation matrix or system matrix, and is given by This system matrix A is Nc × Nc square and is always invertible [36] because it is always positive definite matrix [3,34]. Hence we have where ρ is Nc×1 matrix. The interplant is evaluated using (20) at N evaluation points {x i } N i=1 , through forming (N × Nc) evaluation matrix B which is given as The interplant is then evaluated at N points using matrix vector product to produce u as follows.
which gives approximate solution at any point in .
Where u is N × 1 matrix. As in ROF model, the time marching Eq. (4) can be rewritten as Now, combining Eqs. (24) and (23) we get a new restoration PDE, which is shown in the undergoing non-linear system of equations: where L(u) = u 2 x + u 2 y , u x = H x z, u xx = H xx z, u y = H y z, u yy = H yy z, ∂u ∂n = u n = H n z, and z (0) = z. Since the RBF in the Kansa scheme does not necessarily satisfy the governing Eq. (25), so we have more freedom to choose a RBF. The most popular RBF in the Kansa method is the multiquadric (MQ) [20,34], which usually shows spectral accuracy if an appropriate shape parameter c is chosen. Here, the shape parameter c used in RBF is also one of the most important parameters for the smoothness in our method M3. For the optimal value of c, our proposed methodology gives more accurate and smooth results in image denoising having additive noise. In this technique the shape parameter c and fitting parameter λ depend on the size of the image and the noise level in the image.
Algorithm 3 :Algorithm for proposed method M3 RBF 1. Set N = Nc, n total number of pixel points (N shows the image size i.e., N × N), where N and Nc are called the total number of pixel and center pixel points which are used in the RBF approximation process. 2. Find the ρ according to (21) by multiquadric radial basis function (MQ-RBF) by using step (1). 3. Find u according to (23) by MQ-RBF by using steps (1) and (2).

Results and discussion
In this section, some numerical results are provided to demonstrate the performance of our proposed method M3. Obtained results are compared with the results of methods M1 and M2. The test images are "Lena", "Parrot", "Synthetic1", "Synthetic2", and "Cameraman" which are shown in Fig. 1.
In this paper, it is assumed that N = Nc = the size of the image, for our method M3, for the sake of comparison with methods M1 and M2. Here, Multiquadric radial basis function (MQ-RBF) is utilized for the proposed method M3. To quantify the denoised image, the peaksignal-to-noise ratio (PSNR) is considered. This measure has been commonly used and applied to determine the Whereû is the given image, u is the restored image and M × N is the size of the image. Iterations in our algorithm are terminated when the following condition is satisfied.
where indicates the maximum permissible error. Here, it is set to be 10 −5 .
Here we use the Multiquadric (MQ) RBF to test and compare the results of method M3 with methods M1 and M2. For each point (x j , y j ), Multiquadric (MQ) RBF is defined as equation below.
Example 1: In this first example, the three methods, M1, M2, and M3 are applied and tested on natural images "Lena" and "Parrot" with additive noises (Gaussian noise with mean value zero and standard deviation σ ) with noise levels σ = 20, and σ = 18, which are shown in the Figs. 2 and 3, respectively. In all three figures, (a) and (b) are the original and noisy images while (c), (d), and (e) depict the restored images by the three methods M1, M2, and M3, respectively. In each case, we can clearly notice that the visual quality of restoration by proposed method M3 is quite efficient than that of methods M1 and M2. Method M1 can restore the images, but the quality of restored images are not so good. Also, it creates staircase effect, which is an intrinsic defect of the TV regularization in ROF method as discussed in [7,35]. These reconstructed images are shown in the Figs. 2c and 3c, respectively. In M2, the visual quality of the restored images are good as compared with the method M1. But still, it contains the staircase effect due to the TV regularization. These obtained images are shown in the Figs. 2d and 3d, respectively. In our method M3, the visual quality, preservation of edges, and the minimization of staircase effect of the restored images are far better than that of methods M1 and M2 due to the mesh-free characteristics of MQ-RBF approximation used in our method M3. These denoised images are shown in the Figs. 2e and 3e, respectively. In our method M3, the shape parameter c plays a vital role in image denoising. The range of optimal value for c in this example is set to 1.63 ≤ c ≤ 1.70. Moreover, the PSNR values for the two images "Lena" and "Parrot" for three methods M1, M2, and M3 are listed in Table 2. The bigger the PSNR value, the better the denoising performance. It can be seen from Table 2 that the PSNR values of method M3 are greater than that of methods M1 and M2 for the two images, which shows the dynamic restoration performance of the M3 over M1 and M2. The CPU time of computation and number of iterations required for convergence of the three methods M1, M2, and M3 are also listed in Table 3. It can be observed from Table 3 that the number of iterations and CPU time of computation of the M3 is smaller than that of M1 and M2, which shows the faster restoration performance of our mesh-free algorithm of proposed model M3 over the mesh-based algorithms M1 and M2. The best experimental optimal values of parameters of our technique M3 (shape parameter (c) and fitting parameter (λ) ) for the two images "Lena" and "Parrot" are (1.69, 15.5) and (1.65, 14.7), respectively. So, it is evident from this example, that the performance of our mesh-free based method M3 is superior to that of meshbased methods M1 and M2 regarding visual restoration quality ( PSNR), the number of iterations, and CPU time of computation.
Example 2: In this second example, we study how our algorithm M3 deals with the synthetic images "SynIm-age1" and "SynImage2" having the Gaussian noises. These synthetic images are listed in the Figs. 4 and 5, respectively. The noise level for the two artificial images "synIm-age1" and "synImage2" in this example are σ = 18 and In this example, the best experimental optimal values of the parameters used for the two images "SynImage1", and "SynImage2" are for proposed technique M3 (shape parameter (c) and fitting parameter (λ)) are (1.76, 15.3) and (1.73, 13.9), respectively. Again, the range of the optimal value of shape parameter c for our proposed method M3 in this example is set to 1.72 ≤ c ≤ 1.77. The performance of the three methods M1, M2, and M3, concerning restoration (PSNR values), CPU computation time, and number of iterations required for convergence for the two images, i.e., "SynImage1" and "SynImage2", are recorded in Tables 2 and 3, respectively. Again, from Tables 2 and  3, M3 shows adequate performance than M1 and M2 due to the mesh-free applications of MQ-RBF used in our algorithm M3. Example 3: In this example, the three methodologies M1, M2, and M3 are applied and tested on "Cameraman" having different noise levels. Here, we can observe that the visual quality of restoration (PSNR) of our proposed approach M3 is much better than that of M1 and M2 because of mesh-less applications of MQ-RBFs employed in our proposed method M3 especially when the noise variance is significant. These can be seen from Figs. 6, 7, 8, and Table 4, respectively. Example 4: In this fourth example, we evaluate the performance of three algorithms, M1, M2, and M3, in regions that are dominated by discontinuities by the piecewise contact signals which are shown in Fig. 9.
In areas with edges, we observe that M3 restores and enhances these image features in a productive way compared to M1 and M2, which shows the better performance of edge enhancement of M3 over M1 and M2. These results are shown in Fig. 9c, d, and e, respectively.
Example 5: In this example, we compare the three techniques, M1, M2, and M3, for texture perseveration that are given in Fig. 10. Here, we observed that the textured regions were sufficiently reconstructed by using the M3 as compared to M1 and M2 because of the effectiveness to our proposed method M3. These reconstructed results are shown in Fig. 10c, d, and e, respectively.
Example 6: Here, the homogeneity is checked, and loss (or preservation) is examined for the three techniques M1, M2, and M3 while being applied to "Lena". For this purpose, different lines of the original image compared with noisy and restored images that are shown in Figs. 11 and 12. It is clear that the lines restored by proposed method M3 (shown in Figs. 11d and 12d) are far better than what acquired utilizing methods M1 and M2 that are presented in the figures (Figs. 11c and 12c) and (Figs. 11b and 12b), respectively. The blue line is the original image, and red line is the restored image.

Shape parameter analysis
In this subsection, we compare the image restoration (PSNR) by our proposed method M3 for the different values of shape parameter c for real and artificial images.
Here, we can notice that different values of the shape parameter c affect the image restoration quality of real and artificial images shown in Figs. 13 and 14. The PSNR values of the two images are also listed in Table 5 for different values of shape parameters.

Comparison with an other recent model
Dong et.al proposed a new variational model for image restoration in 2015 [18], and some good restoration results have been obtained on real images. The minimization functional of this model is given as were α and λ are the two regularization parameters. This model can be written a decomposition model as z = u + v + w. Here, u, v, and w are taken as discontinuous, piecewise smooth, and noise components, respectively. By alternating minimization technique and unconstrained transformation technique the author gets the following minimization functional.
where d = ∇ NL u. The split Bregman applied for the solution of functional (29), which is given as The associated Euler Lagrange equation of subproblem (30) for u is The Gauss-Seidel method is utilized to find optimal value of u from Eq. (30), which is  where The optimal value of d can be defined as where shrink(x, γ ) = |x| x .max(|x| − γ , 0). For more information, see [18].
Here, we compare our proposed method M3 with a recent method [18] described above. We can notice that the display results in Figs. 15, 16, and 17 and Table 6 that our model M3 performs well regarding visual quality of restoration (PSNR) than method [18] for the same noise variance taken in method [18].

Sensitivity analysis of parameters
To comment briefly on the choice of the shape parameter (c) and fitting parameter (λ) used in the algorithm M3 described above, it is recommended from our experience that all the two parameters c and λ are more complicated to choose. However, its optimal values are adjusted and tuned according to the noise variance, image size, etc. It has been observed that the range of values allowed are c ∈ [ 1.63, 1.78] and λ ∈ [ 13.2, 16.3], for natural and artificial images. It indicates that all the parameters c and λ are more important for improving denoising performance. Similarly, the number of iterations required for convergence are taken to be in the range [ 7,27] for results with improved PSNR. Thus, the availability of information about the uncertainty of the denoising result on the userchosen parameters (by Trial and Error Method) is helpful to avoid incorrect decisions.
For brevity, for Tables 7 and 8 we shall denote by

Conclusions
In this paper, a new TV based mesh-free algorithm for additive noise removal is presented in which TV regularization (ROF model) is employed in conjunction with MQ-RBF approximation. This algorithm is exploited for the solution of non-linear PDE arisen from the minimization of the associated TV functional of ROF model. The proposed algorithm (Kansa method) is mathematically simple and robust compared with the classical meshbased methods and hence provide more optimal results because of mesh-free applications of MQ-RBF associated with his algorithm. This approach is tested on different artificial and real images for additive noise, and the results are compared with the existing methods. Our experimental results have  shown that the quality of the restoration of images, the number of iterations, and the CPU times with the use of the proposed method are quite good, and the proposed algorithm is quite efficient. We have also noticed that the performance of our proposed method is far better than that of the existing methods regarding restoration quality (PSNR), the number of iterations, and CPU times because of the mesh-free properties of RBF used in our algorithm. The choice of shape parameter c also plays a significant role in this algorithm, which affects the image restoration. The shape parameter analysis has also been discussed here. A comparison with another method in this field is provided as well. However, this method produces an unsymmetrical interpolation matrix. Additionally, sometimes, this approach suffers relatively lower accuracy in boundaryadjacent regions. These problems are under intense study and results will be reported in the subsequent paper.

Appendix
The derivatives in our method M3 for Eq. (25) are given as under: Since Eq. (21) is When we evaluate the derivative at N evaluation points {x i } N i=1 and Nc center points {x j } Nc j=1 , then RBF interpolation, we have or u = Bρ, (36) with (N × Nc) evaluation matrix B, i.e., Then, the derivative becomes from (35) is as follows: or Where    Define H = BA −1 , then above equation (39) can be rewritten as The differentiation matrix can be defined as For second derivative, we have Also The differentiation matrix is well-defined since it is known that the system matrix A is invertible. For any sufficiently differentiable RBF, φ[ r(x)], the chain rule gives for the first derivative, with ∂r ∂x i = x i r .
The second derivative is calculated as follows with