A WAVELET FRAME APPROACH FOR REMOVAL OF MIXED GAUSSIAN AND IMPULSE NOISE ON SURFACES

,


1.
Introduction. Triangulated surface models are widely used in many fields such as computer graphics, reverse engineering, architectural design and terrain modelling. Compared with other surfaces, such as parametric and implicit surfaces, mesh surfaces are quite easy to obtain and operate. Additionally, they can approximate surfaces with arbitrary topology and geometry [2]. Such surfaces are usually acquired by 3D measurement technologies such as digital scanning devices which produce depth maps and ship-based sonar. Due to physical scanning processes and algorithm errors, measured mesh surface models often contain kinds of noise [24]. The noise may severely affect the usability of mesh models, and it is often desirable to remove noise to obtain high-quality surfaces before further processing. The main challenge in this task is to reduce noise and artifacts from an observed mesh surface, while preserving key features such as sharp edges and corners.
The mathematical model for mesh surface denoising can be stated as follows. We denote a triangular mesh by M = (V, E, F), where V = {V (1), · · · , V (n)} is the vertex set, E = {E(1), · · · , E(e)}, E(i) ∈ V × V is the edge set, and F = {F (1), · · · , F (f )}, F (i) ∈ V × V × V is the triangular face set. Here, let V (i) = (V 1 (i), V 2 (i), V 3 (i)) T be the (x, y, z)-coordinates of the vertex V (i). We are interested in computation of a piecewise smooth function S : M → R d on a triangle mesh. The function S can flexibly be chosen to describe vertex positions, texture coordinates, vertex displacements, etc. For a discretization on the triangle mesh, we formulate S by the vector of sample values at the n mesh vertices (1) S : M → [S(V (1)), · · · , S(V (n))].
Let S(V (i)) = S(V (i)) + i N i be the observed noisy surface where i is a random variable and N i := N (V (i)) ∈ R 3 is the normal of the mesh for each vertex V (i). The task of surface denoising is to approximate the true underlying surface S from S. Surface denoising is a fundamental problem in geometry processing and many attempts have been made in the past decade. Based on nonlinear diffusion equations and the theory of geometric evolution problems, some successful models used for image restoration have been extended to denoising surfaces (see e.g. [5,6,7,11,16,27,28] and the references therein). For instance, in [6] an important geometric model for surface denoising was proposed using diffusion and curvature flow. Moreover, authors in [11] established a variational model analogous to the ROF model [23] in the context of geometry processing. Recently, based on the sparsity of the gradient of the face normal field, a variational model in [31] was proposed for mesh denoising. Other approaches for surface denoising include the classical Laplacian smoothing method [2], bilateral filtering [13,17], normal filtering and vertex updating [26]. More recently, based on a tight wavelet frame representation of surfaces, an analysis based surface denoising model was proposed in [9].
When the perturbation in (1) is additive white Gaussian noise, i.e. i ∼ N (0, σ 2 ), it is mostly considered in the above literature for its good characterization of system noise. However, non-Gaussian type noise are also encountered in many real observations due to noisy sensors and channel transmission [15,24]. Moreover, although in image processing, a vast amount of literature is devoted to non-Gaussian noise problems, such as impulse noise, Poisson noise and mixed Poisson-Gaussian noise (see [15,18,25,29] and many references therein), there are few studies in surface denoising with non-Gaussian noise.
In this paper we study how to effectively remove mixed Gaussian and impulse noise in surface functions. These noises are caused by errors in data transmission or faulty memory location in hardware. We assume that either additive Gaussian white noise or impulse noise takes the form of randomly occurring displacement of the vertices along the normal direction. Mathematically speaking, let S be an observed corrupted surface, i.e.
additive Gaussian noise, and 2 N (V (i)), V (i) ∈ V 2 denotes the i.i.d. random-valued impulse noise. The subset V 2 denotes the region where the information of S is missing. It is assumed to be unknown with each element being drawn from the whole region V = {V (i), · · · , V (n)} by Bernoulli trial with a given probability 0 ≤ p ≤ 1. We denote r := |V 2 |/|V| as the level of impulse noise. If V 2 is empty, there is no impulse noise, then the problem in (2) becomes the surface denoising problem for removing Gaussian noise. If the region V 2 is known in advance, this can be considered as a mesh repair problem. If V 2 is a proper subset of V, the observation S of a surface function S is said to be corrupted by the mixed Gaussian and impulse noise. Under the assumption that the region V 2 is unknown in advance, the goal of this paper is to approximate S from the observed surface function S.
The challenge of this problem is to preserve key features of surface functions such as sharp edges and corners, while removing noise and spurious information simultaneously. For image processing problems, the variational model (3) min u∈R m×n F (u) + λΓ(u) has been successfully applied to image denoising, inpainting, deconvolution, etc (see e.g. [3,4,10,15,18,23,29,30]). Here, F (u) is a fidelity term which keeps the image u close enough to the observed image, Γ(u) is a regularization term for modeling a priori knowledge on unknown images, and λ is a positive parameter balancing the two terms. The fidelity term F (u) takes different forms according to different noise statistics and is designed based on the noise characteristic. For example, it is well known that the least square fidelity is used for additive white Gaussian noise. Less well known is the Kullback-Leibler (KL)-divergence fidelity for Poisson noise and 1 -norm fidelity for impulse noise [15,18]. The regularization term Γ(u) is designed based on a priori assumption on original image. One of the assumptions commonly used is the sparsity of the underlying solution in some transformed domain. Such transforms can be discrete gradient used in total variation [23], wavelet tight frame transform [3,4,10,18], local cosine transforms, etc.
We here extend the idea of (3) to surface denoising. To remove mixed Gaussian and impulse noise on surfaces, we choose the fidelity term , and λ is a parameter. This fidelity term was first introduced in [15] for image restoration with mixed or unknown noises. We choose the regularization term Γ(S) = WS 1 := n i=1 |WS(V (i))|, where W is the discrete wavelet frame transform on surfaces, see Section 2.1 for details. The basic idea behind the regularization term is to make use of the interaction between the wavelet frame transform and the 1 -norm. It is well known that the wavelet coefficient sequence of a signal, which is sampled from a piecewise smooth function, is sparse. Furthermore, because of the 1 -norm, the regularization term WS 1 gives preference to a solution S whose wavelet coefficient sequence is sparse, and to keep the edges and corners of surface functions. The extension of image restoration model to surface denoising is not trivial because of the nonlinear nature of the surfaces and the corresponding algorithm [9,31].
Combining the fidelity term and the regularization term together, we have the following wavelet frame based model (4) is an ordinary least squares problem with an 1 -fidelity term and an 1 -regularization term. It means that this problem cannot be solved straight forward by solving only one system of equations. However, if we write the two 1 -terms in a matrix form, (4) becomes a least squares problem with an 1 -fidelity term. Then, iterative solvers like the augmented Lagrangian method (ALM) can be applied to solve (4).
The rest of this paper is organized as follows. In Section 2.1 we introduce the wavelet frame transform on surfaces. In Section 2.2 we derive a wavelet frame based model and the corresponding algorithm to remove the mixed noise on surfaces. In Section 3 we present some numerical examples and discussions. In the last section we give a conclusion.
2.1. Wavelet frame transform on surfaces. Before presenting our proposed regularization for surface denoising, we briefly review a few facts of discrete tight wavelet frame decomposition and reconstruction. The interesting readers should consult [9,10,19] for a more detailed survey. Note that in the discrete setting, an image is a 2D array that can be understood as a vector living on a discrete grid. Then the discrete wavelet tight frame decomposition can be represented as a matrix multiplication W, where W is derived by all the refinement masks (filters) of a wavelet tight frame system. Correspondingly, by the unitary extension principle [22], the reconstruction operator W T can be defined similarly by the transpose of refinement masks and we have W T W = I. In the implementation, these two matrix multiplications are done by using the fast tensor product tight wavelet frame decomposition and reconstruction algorithms instead, which are essentially just the convolution of images by a set of filters.
Noting that the order of the coefficients of each filter in M is harmonious with the neighborhood of the vertex V (k) which forms the column of the data-matrix. We define the (1-level) wavelet frame decomposition on a mesh surface as Then, for numerical simplicity, we define R as Thus, we have RWS = S. We omit discussions on high-level wavelet frame transform on surfaces here and the interested reader is referred to [8,9] for details.

Model and algorithm.
For a given triangular mesh M and a noisy function S defined on M, our task is to remove the mixed Gaussian and impulse noise from S. Here, the function S can be chosen to describe vertex positions, vertex displacements, a piecewise smooth function on a two dimensional smooth manifold, etc. Motivated by the impressive performance of wavelet frame image restoration methods (see e.g. [3,4,10,18,19,25,30]), we also use the wavelet frame coefficients to build a new characterization of functions on mesh surfaces with key features. To remove the mixed noise from S, we consider the following variational model where λ and ν are two positive parameters and W is the wavelet frame transform on surfaces. The first two terms of (6) are the fitting terms with The last term of (6) is a regularization term with which was defined by (5). Since the optimization model (6) is a least squares problem with an 1 -fitting term and an 1 -regularization term, this problem cannot be solved straight forward by solving only one system of equations. To apply certain fast iterative solvers for this 1 / 2 -problem, we first replace (6) with the following matrix form: I is an identity matrix, e and e are vectors of ones, and | · | denotes a column vector of absolute values. Next, we reformulate (7) by introducing a new variable Z: In the following, we apply the augmented Lagrangian method (ALM) [21] to solve (8). First, let us define the augmented Lagrangian function of (8) associated with a given parameter σ > 0: Then, the ALM for solving (8) is summarized as follows: Here, we choose the calculation error ε = 10 −4 . The main task of ALM iteration is to solve the inner subproblem min S,Z L σ k (S, Z; y k ) (9) = min Inverse Problems and Imaging Volume 11, No. 5 (2017), 783-798
First, for fixed S, we consider the minimization of L σ k (S, Z; y k ) in terms of Z. Let Then, the minimization of L σ k (S, Z; y k ) is equivalent to , which has a closed form solution given by: Here, for β > 0, T β is the soft-threshold operator with t βi (η i ) := sgn(η i ) · max{|η i | − β i , 0}. Note that T β (η) is composed of two parts: one is thresholding on the surface function values, T e (σ k (S − S) + y k 1 ); the other is thresholding on the wavelet frame coefficients, T νê (σ k (WS) + y k 2 ). Here, y k := (y k 1 , y k 2 ) T . Substitute (10) into (9), one can check that (11) min where φ ε (t) is the Huber function defined by Therefore, to solve the subproblem (9), we need to find the optimal S of the following unconstrained convex optimization problem (12) min Following [14], we apply the accelerated proximal gradient (APG) method [10,20] to solve (12), where the gradient of H is given by Note that A has full column rank and H(S) is a strictly convex function. Hence, the above minimization has a unique solution.
The APG algorithm is summarized as follows:

end for
Here, we choose the number of iterations p = 25, and the step length of APG algorithm τ = 1. In conclusion, the ALM-APG algorithm for solving (6) is summarized in Algorithm 3.

Algorithm 3 ALM-APG Algorithm
Numerical results and discussions. In this section, we present numerical experiments and compare our method with certain existing methods to verify the effectiveness of our method. The synthetic data are generated by clean functions on meshes with additive Gaussian noise and impulse noise. The meshes are normalized before test. All the examples are tested on a laptop with Intel Core i3 and 6 GB RAM and all models are rendered using flat shading.
We provide three kinds of numerical experiments to test the performance of our method for suppressing Gaussian and impulse noise, and then compare its performance with recently proposed algorithm in [31] and other variational models. The methods are abbreviated as Zhang et al.'s method [31], 2 -variational model, 1 -variational model and our ( 2 + 1 )-variational model in the following. Denoising results of the algorithm in [31] are kindly provided to us by the authors.
There are two common types of impulse noise: salt-and-pepper impulse noise and random-valued impulse noise. Since a vertex coordinate contaminated by randomvalued impulse noise is not as distinctively an outlier as that contaminated by the salt-and-pepper noise, it is more difficult to detect. Thus we only consider the random-valued impulse noise, which is defined as follows: with probability r, the function value S(V (i)) on vertex V (i) is altered to be an uniform random value in the interval between the minimum and maximum value of the surface. With probability 1 − r, the function value S(V (i)) on V (i) is disturbed by additive Gaussian noise N (µ, σ 2 ). The mixed noise level of following examples are given by three parameters r, µ and σ.

Numerical examples.
In this subsection, we show the experimental processes and results on surface denoising of our method. In the first experiment, two noisy surfaces, named 'elephant' and 'bunny', are simulated by displacement of the vertices along the normal direction (see FIGURE 2). We first add mixed Gaussian and impulse noise to clean surfaces, then apply our approach to remove noise. Denoising results are illustrated in FIGURE 2. In the second experiment, we test the performance of surface function denoising. First, we choose manifold Ω 1 to be the unit sphere in R 3 and define a function on Ω 1 as S 1 (x, y, z) := 1 + x 8 + e 2y 3 + e 2z 2 + 10xyz. As in [1], we visualize this function as a kind of offset surface to the manifold, i.e., we consider the new surface as where N ξ denotes the unit outer normal to the manifold Ω 1 at ξ, see  In the last experiment, we consider a graph G, and the vertices set V of G are sampled from a unit sphere. The functions f G : V → R are generated by mapping two images, 'Slope' and 'Eric Cartman', onto the graph G (see FIGURE 4). Let the observed graph data be perturbed by the mixed noise along the normal direction. The denoising results are presented in FIGURE 5.  [31]. It keeps consistent in all parameters of noise level between the two methods.
In   method can preserve most sharp edges and corners better (see the zoomed view in FIGURE 6). In addition, we can also keep the shallow edge well. In  The quality of the denoised result S(V ) can also be measured in terms of signalto-noise ratio (SNR), which is defined by   4. Conclusion. In this paper, we proposed a tight wavelet frame based model to simultaneously suppress mixed Gaussian and impulse noise on surfaces, under the assumption that the impulse noise region is unknown. We considered both surfaces and functions defined on surfaces. Based on topological structure of surfaces, we first defined wavelet frame transform on surface functions. Then we applied the ALM-APG algorithm to solve the proposed model. In addition, we implemented three kinds of experiments to verify the practicability and effectiveness of our method.
In the end, we compare our approach with Zhang et al.'s method [31] and other variational models. Results showed the efficacy of our method.