Modified Bilateral Filter for Feature Enhancement in Mesh Denoising

Mesh models resulting from scanners are inevitably noisy; hence, removing the noise in scanned meshes becomes an essential task in the services using three-dimensional mesh models. Filtering-based methods are simple but have some constraints in eliminating noise because they degrade the features in the mesh models while removing the noise. In this study, we design a feature enhancement filter that is combined with a conventional denoising filter to remove the noise while enhancing the features. The designed enhancement filter is applied only to the feature areas in the mesh model. Results from experiments on synthetic and natural scanned models validate that the proposed method can restore false features by integrating conventional filtering-based methods, and outperforms other state-of-the-art methods.


I. INTRODUCTION
As the demand for services based on three-dimensional (3D) mesh models increases, high-fidelity scanners have recently been developed. The mesh models resulting from scanners are inevitably noisy, even when high-fidelity scanners are used [1]. Therefore, removing the noise from the scanned meshes becomes an essential task in the services using 3D mesh models. The most challenging task in mesh denoising is to preserve features while noise is removed. In the literature, many algorithms have been introduced to remove undesired noise while preserving mesh features, such as sharp edges and corners. The algorithms are classified into isotropic and anisotropic methods, filtering-based and optimization-based algorithms, and data-driven methods, according to their core techniques [2], [8], [17], [24], [28].
In isotropic methods, surface geometry is not considered when noisy data is removed. Vollmer et al. [2] used Laplacian smoothing to remove noise without considering the preservation of sharp features. Taubin [3] proposed a non-shrinking and two-step smoothing method. Desbrun et al. [4] used mean curvature flow for surface fairing based on a simplified mass matrix. Liu et al. [5] introduced a smoothing approach that preserves the volume on triangular meshes. Nealen et al. [6] and Su et al. [7] proposed global systems of equations for the isotropic smoothing method.
The isotropic approaches do not consider features in 3D mesh models; whereas, recent researches have focused on anisotropic techniques, which remove the noise while preserving the feature data. Perona and Malik [8] presented an anisotropic diffusion algorithm, which is a partial differential equation (PDE)-based denoising method. Similarly, Clarenz et al. [9] proposed a diffusion tensor for the anisotropic diffusion method. Ohtake et al. [10], and Ohtake and E. Belyaev [11] introduced the anisotropic diffusion method using the denoised face normal, which is derived from a weighted sum of neighboring face normals. The level set surface model, based on the anisotropic diffusion of surface normal was proposed by Tasdizen et al. [12] to derive denoised mesh models. Hildebrandt et al. [13] presented the prescribed mean curvature-based surface evolution method to avoid volume shrinkage and preserve features. Anisotropic diffusion-based methods preserve features in 3D mesh models more efficiently than isotropic methods. However, when the 3D mesh models are contaminated with heavy noise, there is a limit to the noise the anisotropic methods can eliminate.
The filtering-based method to remove noise is one of the most studied techniques because it can eliminate noise effectively with low computational complexity. Because early methods directly filter the positions of vertices to remove noise [14], [15], [16], these methods do not preserve features when the noise intensity is heavy. To overcome this problem, recent methods include steps for normal filtering and vertex updating. Zheng et al. [17] used bilateral filters to update normal filtering. Zhang et al. [18] proposed a joint bilateral filter, which is based on a well-designed guidance normal field. Later, Zhang et al. [19] introduced a scaleaware normal filter, which uses static and dynamic guidance. Yadav et al. [20] presented a normal filter in a statistics framework with the isotropic regularization factor. Arvanitis et al. [21] developed a coarse-to-fine framework based on graph spectral processing, and subsequently refined the initial output by iteratively smoothening face normals and vertices. Zhao et al. [22] proposed a normal filter with guidance normals computed by graph-cut method. Wang et al. [23] introduced an implicit global optimization function with a feature-aware trilateral filter for guidance.
Optimization-based methods introduce a hypothesis for an ideal denoised mesh, and formulate mesh optimization as a global energy function whose minimization satisfies this hypothesis. The optimization-based techniques assume features are sparse on a clean mesh. He and Schaefer [24] introduced optimization-based algorithm, which involves the design of an edge operator for arbitrary triangular meshes, and the algorithm maximizes flat regions piecewise to remove noise. The algorithm shows robust performance even in the presence of heavy noise. However, the highly nonconvex nature of optimization makes the calculation more complex; moreover, a piece-wise flat effect is seen to occur in this method. Wang et al. [25] constructed a weight optimization method to solve these problems. Zhang et al. [26] proposed filtered face normals, which use the total variation of face normals to measure sharp feature sparsity. These optimization-based methods are suitable only for simple models because they are based on the hypothesis that the mesh consists of plain regions.
Recently, various studies based on data-driven methods have also been conducted. A data-driven method was proposed by Diebel et al. [27]. They formulated a pairwise normal potential function for learning, which considers application-specific geometry shapes. More recently, Wang et al. [28] used local geometry descriptors as inputs for training neural networks. The network consists of singlehidden layer feed forward networks in a cascaded way. Wang et al. [29] proposed a similar design that uses local geometry descriptors, where a second normal estimation is utilized to recover features lost during the first step. Li et al. [30] constructed a convolution neural network based on a nonlocal similarity approach. Zha et al. [31] also proposed a CNN-based method, including a cascaded structure. The performance of the data-driven methods depends on the quality of the training dataset. Furthermore, the computation cost of the training process for these methods is high.
In addition to the aforementioned classified methods, there are other methods that use a normal voting tensor [32], [34], classify feature areas [32], [33], and use non-local similarity [35], [36]. Fan et al. [32] used the eigenanalysis of a normal voting tensor (NVT) to separate meshes into feature and nonfeature regions, where denoising algorithms are applied to each region differentially. Bian and Tong [33] separated meshes using the k-means clustering algorithm. Yadav et al. [34] used element-based NVT (ENVT), which is a method similar to [32] and [33], but does not include a process of dividing the mesh area into smaller regions. In this algorithm, mesh denoising is performed using binary optimization on the eigenvalues of the ENVT. In the non-local similarity method, Wei et al. [35] and Li et al. [36] proposed low-rank matrix recovery by co-filtering similar patches. Rudin et al. [37] first introduced total variation (TV) regularization for 2D image denoising. Owing to its good edge-preserving property, it has many extensions for 2D image restoration and denoising [38], [39] as well as for mesh denoising [40]. Zhang et al. [40] combined TV and anisotropic Laplacian regularization. Duchamp and Stuetzle [41] proposed an extension of spline smoothing using a finite element method where sharp features are not preserved. By contrast, the surface reconstruction technique [42] based on splines shows impressive performance.
Among the aforementioned conventional algorithms, filtering-based method is one of the simplest algorithms to eliminate noise. However, the method may fail to preserve features during noise removal, because it applies the denoising filter several times to mesh data. To overcome this problem, we design an efficient enhancement filter, which can be jointly used with filtering-based methods. The combined method removes the noise corruption in mesh data, while enhancing the features in the mesh structure.
This paper is organized as follows. We review the related works in Section II. In Section III, we explain the motivation for our algorithm. In Section IV, we propose an efficient algorithm to remove the noise corruption in mesh data while enhancing the mesh features. The performance of the proposed algorithm is evaluated and compared with those of various conventional methods in Section V. Finally, we present our conclusions in Section VI.

II. RELATED WORKS
This section describes filtering-based methods, which are simple although they may degrade feature information during removal of noise in the mesh structure.
Most recent filtering-based methods, [17], [18], [20], and [43], consist of two-step processes, which modify face normals in the first step and then update vertices according to the modified face normals. Sun et al. [43] classified two-step filtering-based methods into one-stage and two-stage techniques. In one-stage techniques, the two steps are coupled as a pair to give a single overall iteration procedure, whereas two separate iteration phases are performed in twostage techniques.

A. NOTATION
In this section, we define some notations that are used in filtering-based methods. = { , , … , } is a set of faces, where is the face of the mesh. = { , , … , } is the set of coordinates of vertices. ℕ = { , , … , } is the set of normal vectors, where is the normal vector of . The center coordinate and area of are denoted by and , respectively.
denotes the set of edges that constitute the boundary of . We denote as the 1-ring vertex neighborhood of a vertex . Similary, the set of faces that share a vertex is denoted by . For a , the set of 1-ring neighbor faces is defined as follows: The set of neighbor faces for is defined as follows: where is the radius determined by a user or is a fixed value. Fig. 1 shows an example of , Φ , and Ω , where red face is , green faces are included in Φ . Ω consists of blue and green faces.

B. BILATERAL NORMAL FILTER
The bilateral normal filter (BNF) method is a two-stage technique, which updates normals ! times and then updates vertices " ! times. The process of the denoising algorithm using the BNF can be represented as #$%& ! () *&+, -. + #$%& ! 0! 1!, 2 -. . The BNF updates the normal vector using the weighted average of noisy neighbor normal vectors as follows: n n n n (5) Function 7 8 gives small weight when the distance between the centroids (i.e. the geometric centers) of and neighboring faces is large. Function 7 9 gives more weight as the similarity between and neighboring normals increases.

C. GUIDED NORMAL FILTER
The guided normal filter (GNF) method is a modified onestage technique. The process by which the noise removal algorithm uses the GNF method can be represented as The updated normal vector is derived as follows: where ; is a guided normal vector of and is designed to be robust to noise. For a given face , a patch is selected among a set of candidate patches that contain the face. The selected patch has the most consistent normal directions. The average normal of the chosen patch is used as the guidance normal ; for the face . The GNF preserves features better than BNF because ; is more robust to noise than [18]. Notably, Ω is used in GNF whereas faces in Φ are considered in BNF.

D. ROBUST AND HIGH-FIDELITY FILTER
The robust and high-fidelity filter (RHF) method is a onestage technique that has one iteration parameter. The process of the RHF method can be represented as #$%& ! () *&+, + #$%& ! 0! 1!, -. . The RHF updates normal vectors in a similar manner as BNF: where Tukey's bi-weight function is used instead of Gaussian function 7 9 as follows.
Tukey's bi-weight function preserves the features better than Gaussian function because it does not allow diffusion across sharp features. The RHF method removes noise components along sharp features using Tukey's bi-weight function.

E. VERTEX POSITION UPDATING
The two-step techniques, including BNF and GNF, update the coordinates of the vertices using the method proposed by Sun et al. [43]. They constructed an error function using the fact that must be perpendicular to the updated normal @ . Therefore, to minimize the error function, the gradient descent method is used as follows.
RHF also updates the coordinate of vertices using this concept; however, this involves an additional term, AB . The term B is a combination of the differential coordinate and tangential component. It makes the algorithm robust against high-intensity noise. Because term B has a similar effect to isotropic smoothing, using the same A value in all iterations causes a shrinkage problem. To overcome this problem, the value of A decreases as the iterations increase. iterations respectively. From Figs. 2c, 2d, and 2e, we observe that features are not preserved during all iterations. This is because some false features are regarded as features. Fig. 2f is the colored version of Fig. 2e, where the faces, including the false features and correct normals, are painted red and blue, respectively. BNF, GNF, and RHF methods give high weights to neighboring normal, which are similar to the current normal in (3), (6), and (7), respectively. Thus, when algorithms are applied to one of the normals in the red area of Fig. 2f, high weights are assigned to only the neighboring normals in the red area. As a result, even if filtering is performed several times, the noisy normals cannot be corrected. This problem is due to using the distance between and C in the weighting functions of (3), (6), and (7).

III. MOTIVATION
In this paper, we propose an efficient filtering process for normal vectors to solve the aforementioned problem. The basic idea of the proposed algorithm is to recover features using the normals in the blue areas when filtering the normals in the red areas in Fig. 2f. The proposed algorithm consists of denoising and feature enhancement modules, where the denoising part is implemented by a conventional algorithm, such as BNF, GNF, and RHF, based on (3), (6), and (7), respectively. The enhancement module is a form of the filterbased technique that can be used to fix noisy features. When the proposed enhancement filter is applied to the normal vector of the current face , we assign the high-weighted values to the normal vectors of the neighbor faces, which have the feature property. Details of the proposed algorithm and the enhancement module are described in Section IV.

A. DEFINITION OF SET D
We define a direction vector from current face to one of the neighboring faces E as follows: As shown in Fig. 3, F E is a vector from the center of face to that of face E . According to the colors used in Fig. 1, is marked in red, faces in the set G are green, and faces of blue and green regions are in the set H . Based on these notations, we define the set I E for each E in G as follows:  Fig. 4 shows an example of set I which includes , K , and L where directions of F K and F L are closer to the direction F than other F M ∀ M ∈ G , * ≠ 100). We construct the sets I E s for all one-ring neighbor faces E s. For example, we can construct I , I K , I S , I T , I TU , I KK , I U , I , I L , I SS , I SU , and I S in Fig. 4.

B. DEFINITION OF SET V
For each face E in the one-ring neighbor G , we construct the set Γ E ,which consists of indices of neighboring faces satisfying the following condition. 1
Γ E is a set of indices of faces X s, which are in one-ring neighborhood and whose F X s are within an angle Y from F E . If Y is set to a small value, the subsequent processes are noise-sensitive. By contrast, if Y is set to a large value, the computational complexity of the algorithm increases substantially. It means that the selection of Y has a trade-off relationship between noise sensitivity and computational complexity. Y is set to 60° in our algorithm after experimental tests. We construct the sets Γ E s for all one-ring neighbor faces E s. Based on the sets Γ E s, we construct the sets \ E s for all one-ring neighbor faces E s as follows. .

C. FEATURE AND NON-FEATURE REGIONS
To determine whether the current face is in a feature region or not, b (b&cd& ) is calculated as follows. (14) where |\ E | is the total number of faces in \ E . In the feature region, the averaged dot product value of the current normal vector and the normal vectors f , belonging to \ E depends on the geometric characteristic of \ E . If set \ E contains normal vectors that are similar to the normal vector of the current face such as \ SS in Fig. 5, the averaged dot product value is approximated to 1. By contrast, if \ E contains normal vectors whose directions are different from the normal vector of the current face such as \ KK or \ TU in Fig. 5, the averaged dot product value decreases to −1. As the averaged angle between and f , increases, the averaged dot product decreases to −1 through 0 from 1.
In the feature region, as the normal vectors included in \ E have a larger angle with the normal vector of the current face, the averaged dot product value decreases to −1 through 0 from 1. By contrast, in a non-feature region, regardless of the location of \ E , the normal vector of the current face and the normal vectors of the neighboring faces have similar directions, so the averaged dot product values for all \ E s are similar and close to 1.
In the numerator of (14), we calculate the averaged dot product value between the normal vector of the current face and normal vectors f of faces f in set \ E . Among them, we select the minimum value that occurred when the averaged angle between and f s is largest. If is in a feature region, the minimum of the numerator in (14) is a small value (down to −1), because there is at least one \ E whose normal vectors are considerably misaligned to of the current face . By contrast, in the non-feature region, because the normal vectors and f s are mostly aligned, the minimum value is approximately 1.
When the mesh model consists of sphere-shaped parts, although the faces are in the non-feature region where edges are not included, and f s are not aligned, and the minimum value of the numerator of (14) is not approximated to 1. This tendency differs from that of the non-feature regions that do not include sphere-shaped parts. Therefore, in order to classify a sphere-shaped part as a non-feature region, the averaged dot product is normalized by the denominator of (14).
In the denominator of (14), the maximum value is selected among the averaged dot products. In the homogeneous region (i.e., non-feature region) that does not include the sphere-shaped part, most normal vectors f of faces in the set \ E are aligned with . Thus, the maximum value is approximately 1. By contrast, in the feature region, although many faces are misaligned with the current face, there is at least one normal vector of faces in set \ E that is approximately aligned to the normal vector of the current face, because the extremely sharp mesh-like needle-shape rarely occurs. Therefore, in the feature region, the maximum value of the denominator is approximately 1. Note that the denominator of (14) is approximately 1 for the feature as well as non-feature regions (without including the sphere-shaped part). When the non-feature region consists of the sphereshaped parts, the value of the denominator of (14) is approximately equal to that of the numerator. Therefore, the b (b&cd& ) of (14) is approximately 1 for the nonfeature region (regardless of whether sphere-shaped parts are included), whereas it has a much smaller value than 1 for the feature region.
However, we can consider a special case of the feature region where all the neighboring faces of the current face are perpendicular to the current face, although this scenario rarely occurs. In this case, the numerator and the denominator are both 0. In order to avoid the denominator being 0, a very small value of l is added to the denominator in (14). In this special case, the b (b&cd& ) in (14) is 0. Therefore, by comparing the b (b&cd& ) with a threshold value between 0 and 1, the considered regions can be classified into feature or non-feature regions. In our algorithm, l and the threshold m n are set to 10 op and values in the range [0.7, 0.9], respectively. Fig. 6a visualizes the feature regions by using b (b&cd& ) in (14), where the red and blue colors are set to the faces as much as 1 − b (b&cd& ) and b (b&cd& ) of the faces, respectively. In Fig. 6a, to represent it clearly, the assigned color values are clipped in the range of [0, 255]. In this study, when b (b&cd& ) of is less than a threshold m n , we consider that is in a feature region. Using the thresholding, Fig. 6b shows the feature regions with yellow color.

D. TARGET NORMAL VECTOR
In this section, we define the target normal vector r , which is derived from the neighbor faces whose normal vectors are aligned to of the current face . In order to select those neighbor faces, among all \ E s, we choose one whose normal vectors are most aligned to of the current face according to the following equation.
We can explain the meaning of (15) with an example of Fig.  5, where because faces in the orange area are over a feature, the value of (15) corresponding to the orange area is small. On other hand, the pink area has a large value of (15) because all normals of the pink region are similar to . In the example of Fig. 5, the value of s * is set to 88. From s * , the target normal vector r of a in a feature region is derived as follows: where u c is a normalization function of c . 7 8 is a Gaussian function and 4 5 controls the kernel widths of Gaussian function 7 8 .

E. FEATURE ENHANCEMENT FOR NORMAL VECTOR
In our algorithm, the feature-enhanced normal @ of the current face is calculated using the following filter: (17) where 7 8 is a Gaussian function and 7 v is a Tukey's biweight function. 4 5 ′ and 4 6 ′ control the kernel widths of Gaussian functions 7 8 and 7 9 , respectively. The function 7 8 gives high-weighted values to normals that are near the set I E * . For example, if s * is 88 in Fig. 5, high-weighted values are assigned to the normals of faces belonging to I L and I SU because they are close to I E * . By contrast, low weights are given to the normals of the sets far from I E * , such as I TU , I S , and I S . In (17), the center position of the Tukey's bi-weight function is shifted from the current normal vector to the target normal vector, where the target normal vector is computed using (16). The function 7 v enhances features by giving high weights to neighboring normals, which are similar to the target normal vector r of the current face, . Tukey's bi-weight function produces more robust results than the least square and Gaussian-based estimators because it completely cuts off the diffusion at sharp features [20]. As (17) is a filter to enhance the feature, it is more desirable to avoid assigning weights to outliers whose direction is very different from the target normal vector r , so we used Tukey's bi-weight function. Notably, (17) enhances the features; it does not remove the noise. Consequently, the enhancement module based on (17) is combined with a conventional denoising algorithm, such as

F. CRITERION FOR IMPROVEMENT
The proposed technique applies filters to the noisy data iteratively, as many conventional algorithms do, to remove the noise. However, the proposed algorithm performs the filtering for denoising and that for feature enhancement alternatively. When the denoising filter is applied, one of the conventional denoising methods, such as BNF, GNF, and RHF, is used to remove the noise iteratively. In the denoising step, if the improvement effect of the denoising filter is saturated, we apply the feature-enhancing filter of (17). After the features have been enhanced, the denoising process is applied alternatively.
In the end of each iteration, we calculate the improvement factor (IF) to check whether the improvement resulting from the applied filters has been saturated. To calculate the IF in the ! iteration, alteration x of the normal vectors is evaluated as follows.
where y E s are normal vectors obtained after the vertices have been updated by (9) according to the filtered normal vectors @ in (3), (6), (7), and (17). v is the total number of total normal vectors. (18) implies the averaged absolute difference between angles of E s and y E s. To check the IF of the ! iteration effectively, we calculate the moving average of the first derivative of x -. as follows.   Fig. 7 illustrates the flowchart of the proposed algorithm, which consists of a (a) denoising process, (b) feature enhancement, and (c) denoising process without considering feature enhancement. In the denoising part, one of the conventional denoising filters in (3), (6), and (7) is used to remove the noise in the mesh model. Then, the vertices are updated by (9). At the end of each iteration of the denoising part, we evaluate the IF of (19) to check whether the noise is removed sufficiently. If z { -. is less than a threshold m { , the flag of FE is set to 1 and the process is switched to the feature enhancement part (b) in the next iteration. z { -. is the value of z -. calculated with data resulting from part (a). In part (b), we apply the feature enhancement filter of (17) to enhance the degraded features. The feature enhancement is not applied consecutively, because consecutive feature enhancements may amplify the noise. Thus, FE is set to 0 after each application of the feature enhancement. At the end of each iteration of the (b) feature enhancement part, if z |X -. is less than a threshold m |X , we decide that the feature enhancement is saturated sufficiently and no more enhancement is needed by setting NME = 1. Notably, z |X -. is the value of z -. calculated with the data resulting from part (b). After NME is set to 1, only the (c) part is performed until the end of the iteration.

G. CONTRIBUTIONS
Conventional bilateral filter-based methods [17], [18], [20] used Gaussian functions to assign higher weights to the normal vectors whose corresponding faces' centroids are near that of the current face. In addition, a Gaussian function or Tukey's bi-weight function was used to assign higher weights for the neighbor normal vectors, which are similar to the normal vector of the current face. While the conventional methods consider distance as one of the most important parameters, the proposed method considers the direction from the location of the current face to those of the neighboring faces. The proposed method gives higher weights to normal vectors because the locations of neighboring faces are nearer on the line of direction set in (15). This is done so as to give low weights to the neighboring normal vectors across the edge, even though the distance from the current face is low. We proposed the set I E of (11) to define the direction from the location of the current face to those of the neighboring faces. We also defined sets Γ E , \ E of Eqns. (12), (13) to determine the direction in which to assign high weights.
Some conventional mesh denoising methods [25], [34] used covariation matrices and principal component analysis (PCA) to classify meshes into feature and non-feature regions. These classification processes have high computational complexities. By contrast, we proposed a simple method to distinguish between feature regions and non-feature regions, where the parameter b (b&cd& ) of (14) is used. According to the value of (14), the feature enhancement filter can be selectively applied to feature regions.
We also modify Tukey's bi-weight function to assign higher weights for the neighbor normal vectors that are similar to the target normal vector. In (17), the center position of the Tukey's bi-weight function is shifted from the current normal vector to the target normal vector, where the target normal vector is computed using (16). As the target normal vector in the feature region is more robust to noise than the current normal vector, the proposed filter better preserves the feature.

V. EXPERIMENTAL RESULTS
In this section, we evaluate the performance of the proposed algorithm and compare it with those of various conventional methods. As we can see in Fig. 7, the proposed algorithm imports one of the conventional denoising methods, such as the bilateral normal filter (BNF) [17], the guided normal filter (GNF) [18], and the robust and high-fidelity filter (RHF) [20] to remove noise from mesh models. The performance of our proposed algorithm is compared with those of other methods such as + minimization (L0M) [24], the non-local low-rank filter (NLLR) [36], and the cascaded normal regression (CNR) filter [28] as well as BNF, GNF, and RHF. The source codes of GNF, RHF, and NLLR have been kindly provided by their authors, whereas BNF and L0M are implemented by a third party. The trained neural networks of CNR have been kindly provided by their authors [28].
When the conventional methods are integrated with the proposed algorithm, the parameters of the methods are set to harmonize with those of the proposed algorithm. When BNF is used as in [17], it requires two iteration parameters to be set because it is a two-stage method. However, when BNF is integrated with the proposed method, one iteration parameter is required. As RHF is a one-stage method like the proposed method, it can be used in the proposed algorithm without any alteration.

A. PARAMETERS SETTING
The proposed algorithm has some parameters, such as radius in (2), a threshold m n to be compared with b (b&cd& ) of (14), 4 5 } and 4 6 } in (17), a threshold m { to check the effect of the denoising, a threshold m |X for the enhancement.
The radius parameter represents the range of neighborhood of the current face. BNF constitutes neighboring faces with 1-ring neighbor. GNF calculates = 1~× € for a radius to constitute neighboring faces, where 1~ is the average distance among the centroids of all faces, and € is a parameter set by the user. In RHF, is set to 1~× 2. When BNF and RHF are respectively integrated with the proposed algorithm, neighboring faces are constructed in their ways, and when GNF is used, € is set to 2. In (11)- (17) , is set to 1~× 3.5. m n depends on the noise level and is set to values in the range [0.7, 0.9]. From empirical test data, we found that the performance of the proposed algorithm is not sensitive to m n .
The parameters 4 5 and 4 6 in (3), (6), and (7) are set to values used in BNF, GNF, and RHF. The 4 5 } in (17) is set to the average distance 1~ among the centroids of faces as was done by Yadav et al. [20]. 4 6 } is set by the user as in [18], [20]. From empirical test data, we found that m { and m |X are independent of noise levels and characteristics of the mesh models. In the simulations of this section, we set m { and m |X as 1.0 and 0.15, respectively.

B. TEST MODELS
To evaluate the performances of conventional and proposed methods, we used various mesh models, such as Fandisk,  [28]. Fandisk, Block, Sharpsphere and Twelve contain distinct feature areas, whereas Nicolo, Julius and Bunny consist of smooth surfaces. For CAD-generated models, we add impulse noise to Twelve model, and we corrupted other CAD-generated models with Gaussian noise.

C. OBJECTIVE EVALUATIONS
The two most widely used error metrics, mean square angular error (MSAE) and vertex-based surface-to-surface error [44], were used for objective evaluations of the models used in our study. MSAE represents the mean angular difference between the matching normal vectors of the ground-truth and denoised models as follows, where is the normal vector of ground-truth, } is the normal vector of denoised mesh, and |ℕ| is the number of normal vectors. The vertex-based error is a metric that calculates how close the positions of vertices in the denoised model are to those in the ground-truth model. The vertexbased error is defined as follows, (22) where ‚ is the nearest face in the ground-truth mesh to the vertex } of the denoised mesh. The vertex-based error is calculated based on the weighted average of the distance between vertices, where the weighting factors are areas of the surrounding faces of the corresponding vertex. As discussed in [23] and [45], MSAE metric is more correlated to visual quality than vertex-based metric.
In Table 1, the performance of the proposed algorithm incorporating BNF is compared with that of BNF. In simulations with the CAD models, the parameters of BNF were set to the optimal values provided in [18]. The parameters (σ 6 , ! , " ! ) of BNF for Cone and Pyramid models and the parameters (4 6 } , m n , ! ) of the proposed algorithm for all models were set to the optimal values, which provide the best performances in the empirical tests.
From Table 1, we can observe that the proposed method outperforms BNF for all models in respect to MSAE; the proposed algorithm reduces MSAE significantly in the models Fandisk, Block, Sharpsphere, Twelve, and Pyramid, which contain many distinct features. As for Nicolo, Julius,   10) and Bunny models, which have few features and consists of smoothed edges, the proposed algorithm slightly improve the performance in respect to MSAE compared to BNF.
In the " … column, the proposed algorithm outperforms BNF in all models except for Nicolo and Cone models. The performance of the proposed algorithm is approximately equal to that of BNF for Nicolo and Cone models. In the models Fandisk, Block, Sharpsphere, Twelve, and Pyramid, which contain many distinct features, " … is significantly improved.
As we observe in Table 1, the proposed algorithm provides large improvement for models having distinct features, because the main idea of the proposed algorithm is to enhance the features while the defects in the plain regions are removed. Table 2 represents the performances of GNF and the proposed algorithm incorporating GNF for various models. When GNF is implemented, the parameters (€, σ 6 , ! , " ! ) for the CAD models are set to the optimal values provided by [18]. The parameters of GNF for Cone and Pyramid models and the parameters (4 6 } , m n , ! , " ! ) of the proposed algorithm for all models are set to the optimal values, which provide the best performances in the empirical tests.
From MSAEs of BNF and GNF in Table 1 and 2, we observe that MSAEs of GNF are much lower than those of  BNF for the distinct feature models Fandisk, Block, and Twelve. It means that GNF preserves features better than BNF.
In Table 2, in the MSAE column, the proposed algorithm outperforms the GNF method for all models. We observe that the tendencies of MSAE improvement of the proposed algorithm are similar in Tables 1 and 2. From the results for Sharpsphere including the curved sharp edges in Table 2, we observe that MSAE and " … of the proposed algorithm are much smaller than those of GNF. This implies that the proposed method preserves the curved sharp edges better than the GNF method. In the column of " … in Table 2, the proposed algorithm is slightly inferior to GNF in some models, such as Fandisk, Nicolo, Julius, Bunny, and Cone.
The comparison between performances of RHF and the proposed algorithm incorporating RHF is summarized in Table 3. The parameters (σ 6 , A, ! ) of the RHF method for Fandisk and Nicolo models are set to the optimal values provided by [20]. The parameters of RHF for other models and the parameters ( 4 6 } , A , m n , ! ) of the proposed algorithm for all models are set to the optimal values, which provide the best performances in the empirical tests. The tendency of the results shown in Table 3 is very similar to that of Table 2, except for results of the Sharpsphere model.
In Table 4, the performances of the proposed algorithms respectively integrated with BNF, GNF, and RHF methods are compared with those of L0M, NLLR, and CNR, where the best and second-best performances are shown in red and blue colors, respectively. As we observe from Table 4, one of the proposed algorithms shows the best performances of MSAE and " … for all CAD models, such as Fandisk, Block, Sharpsphere, Twelve, Nicolo, Julius, and Bunny. As for Cone and Pyramid models, the proposed algorithm has the second best MSAE and the best " … . Notably, Cone and Pyramid models consist of overall flat regions, except for few feature regions over small areas. This property of those models is matched to the hypothesis of the L0M method. It is the reason why L0M has the best MSAE values for the Cone and Pyramid models.
In Table 5, the time required to complete the processes of various algorithms is compared. BNF is the fastest method because it only considers the distance between the current face and neighboring faces and the difference between the current normal vector and neighboring normal vectors. However, it results in degraded meshes when the models include distinct features, as shown in Figs. 8b and 9b. In the case of OURS+BNF, although the time consumed increases slightly, enhanced results are obtained, as shown in Figs. 8c and 9c. The GNF method needs the longest computation times because it requires many iterations for updating vertices, where the total number of iterations of updating vertices is " ! × ! . In OURS+GNF, the number of iterations of updating vertices is substantially decreased because the proposed filter effectively removes noise with a small " ! while preserving features satisfactorily. For example, in the Block model, the GNF method has a ! value of 40 and a " ! value of 30, resulting in 1200 iterations of updating vertices. By contrast, because the OURS+GNF method has a ! value of 70 and a " ! value of 2, only 140 iterations are performed. Thus, the consumed times of OURS+GNF are smaller than those of GNF. When RHF and OURS+RHF are compared, OURS+RHF requires slightly higher computational times than RHF, because the sets I E , Γ E , \ E of (11), (12), (13) and the target normal vector should be calculated. However, OURS+RHF produces better results than RHF, as shown in Figs. 8g and 9g, in a reasonable time, because the feature enhancement filter of the proposed algorithm is no longer applied when the feature enhancement is saturated according to (19).

D. SUBJECTIVE EVALUATIONS
In this section, the subjective performances of various algorithms are evaluated for a variety of mesh models. In Figs. 8,9,10, and 11, we can visually compare the meshes denoised by BNF, OURS+BNF, GNF, OURS+GNF, RHF, and OURS+RHF.
Figs. 8 and 9 show the results for Fandisk, Block, Sharpsphere, and Twelve models, which have distinct feature areas. In the Fandisk model, challenging regions, including corners, are indicated by two red windows. As shown in the windows, the proposed method ensures that the features are well reconstructed. The region in the red window of the Block model consists of triangular faces whose density is very high. Whereas most conventional methods do not recover features in this region, the proposed method efficiently reconstructs the features in the challenging areas. In the enlarged window of the Sharpsphere model, we observe that the BNF method fails to recover the features. Conversely, the GNF method reconstructs its features to some extent, but does not recover its curved edge. The proposed algorithm, integrated with the BNF or GNF method can perfectly reconstruct the features and effectively enhance the curved edges. In columns (f) and (g) of Fig. 9, both the RHF method and the proposed algorithm incorporating RHF show excellent reconstruction and enhancement of features of the Sharpsphere model. In the Twelve model, one of the edges is enlarged in the red window. Among the results (b), (d), and (f) obtained from the BNF, GNF, and RHF methods for the Twelve model, respectively, the GNF method results in the best visual quality, but the left end of the edge is distorted in the result. Conversely, all denoised results by the proposed algorithm show perfectly reconstructed edges. Note that the BNF method is not effective for impulse noise, but the proposed algorithm preserves the features corrupted from impulse noise. Fig. 10 shows the denoised meshes of the Nicolo, Julius, and Bunny models, which have sparse features. In the Nicolo model, there is no significant difference between visual qualities of meshes resulting from BNF and the proposed algorithm incorporating BNF. When GNF method is applied on the noisy Nicolo model, some wrong features are generated in the nose area of the resulting mesh. Conversely, the proposed algorithm using GNF provides the enhanced features in the nose area. In the test of RHF for Nicolo model, it blurs the model excessively, whereas the proposed   area of a set of neighboring faces is also very small, the proposed algorithms (OURS+BNF and OURS+RHF) do not recover the noisy features, which are located actually in the near neighbor but out of the set H . When the Pyramid model is used as a test data, we observe from results of the enlarged red windows that the conventional methods of BNF, GNF, and RHF blur out corners and edges, whereas the proposed algorithm, integrated with those conventional methods, reconstructs the features efficiently. In Figs  In Fig. 12, the Fandisk models denoised by L0M, NLLR, and CNR methods have much more degraded features in the corner regions than those by the proposed algorithms. When the Block model is used, the L0M method reconstructs the features in the corner regions better than the NLLR and CNR methods. However, L0M does not reconstruct the curved area of the model. Conversely, the proposed method respectively incorporating the BNF and RHF methods, enhances the features in the corner areas efficiently. In the test on the Block model, OURS+GNF show the best performance in the corner areas. In the simulations with the Sharpsphere model, which consists of curved surfaces and edges, the proposed method outperforms the other methods. The results for the Twelve model have a similar trend as those of the Sharpsphere model. In Figs. 12 and 13, we observe that the L0M method has difficulty in making the curved surface into a set of small flat areas, the NLLR method blurs the edge, and the CNR method does not recover the details of the features. Overall, the results using the proposed method have better performances than each of the L0M, NLLR, and CNR methods. In Fig. 14, all results denoised by the L0M method have a set of small flat areas, which degrade the visual quality significantly. Regarding removing the noise, NLLR shows the worst performance among the compared algorithms. OURS+BNF and CNR show the best performances. However, CNR does not recover the details of hair above the ear in the Julius model. OURS+GNF and OURS+RHF enhance the features efficiently.
For the Cone and Pyramid models of Fig. 15, L0M method shows good performance. In the results from NLLR, the edges are slightly blurred and the noise in the flat area is not completely removed. The CNR method removes the noise in the flat area effectively, but the edge feature is considerably blurred. OURS+GNF and OURS+RHF efficiently removes the noise in the smooth surfaces while simultaneously enhancing the edge.

E. CONVERGENCE OF THE PROPOSED ALGORITHM
Figs. 16 and 17 show the convergence patterns of z -. of (19) for various mesh models, such as Fandisk, Block, Nicolo, and Julius, when the proposed algorithm, respectively incorporating the BNF, GNF, and RHF methods, is applied to remove the noise. In these simulations, we corrupted the test models with low-to high-levels of noise. The x-axis and y-axis in these figures represent the iteration number and z -. , respectively. The graphs in the figures show that z -. converges as the iteration number increases. It implies that the change of the normal vectors decreases and z converges to zero regardless of the characteristics and noise level of the models, as the iteration number increases.

F. ROBUSTNESS OF THE PROPOSED ALGORITHM
This section explains the variations in the performance of the proposed algorithm with changes in the value of m n , which is set by the user. Fig. 18 shows the MSAEs of the denoised models when the proposed algorithm is applied to the noisy Fandisk, Block, Sharpsphere, Twelve, Nicolo, Julius, and Bunny models with various values of m n . Fig. 18 shows that the performance of the proposed algorithm, respectively integrating BNF, GNF, and RHF method, is not affected by changes in the values of m n .

VI. CONCLUSION
This paper proposes a method to enhance the features in noisy mesh models. We found that the conventional filteringbased methods update the current normal by giving high weights to neighboring normals, which are most similar to the current normal. Due to this procedure, the normals of the neighboring faces having the wrong dominant features affect the update the normal of the current face. This constrains the performances of the filtering-based methods. Conversely, the proposed algorithm drives the target normal vectors and enhances the degraded features based on the target normal vector, while the noise is removed by integrating one of the conventional filtering-based methods. From various simulations, we showed that the proposed algorithm is efficient at both removing the noises and enhancing the features.