Sage Journals: Discover world-class research

Abstract

The k-nearest neighborhoods (kNN) of feature points of complex surface model are usually isotropic, which may lead to sharp feature blurring during data processing, such as noise removal and surface reconstruction. To address this issue, a new method was proposed to search the anisotropic neighborhood for point cloud with sharp feature. Constructing KD tree and calculating kNN for point cloud data, the principal component analysis method was employed to detect feature points and estimate normal vectors of points. Moreover, improved bilateral normal filter was used to refine the normal vector of feature point to obtain more accurate normal vector. The isotropic kNN of feature point were segmented by mapping the kNN into Gaussian sphere to form different data-clusters, with the hierarchical clustering method used to separate the data in Gaussian sphere into different clusters. The optimal anisotropic neighborhoods of feature point corresponded to the cluster data with the maximum point number. To validate the effectiveness of our method, the anisotropic neighbors are applied to point data processing, such as normal estimation and point cloud denoising. Experimental results demonstrate that the proposed algorithm in the work is more time-consuming, but provides a more accurate result for point cloud processing by comparing with other kNN searching methods. The anisotropic neighborhood searched by our method can be used to normal estimation, denoising, surface fitting and reconstruction et al. for point cloud with sharp feature, and our method can provide more accurate result comparing with isotropic neighborhood.

Keywords

k nearest neighbors neighborhood searching hierarchical clustering normal estimation

Introduction

In the field of reverse engineering application, generally, two methods including contact and non-contact scanning devices are used to acquire 3D point cloud of objects,¹ but the collected data is 3D point cloud commonly. It can realize surface reconstruction and product remanufacturing by processing the point cloud. However, there is not any topological structure for scattered point cloud only containing 3D coordinate information.

Most point cloud processing methods depend on the neighborhood structure of point. Efficiency and accuracy of neighborhood searching method directly affect the speed and result of point cloud processing including normal estimation,² point cloud registration,³ denoising⁴ and surface reconstruction.⁵ For example, estimating normal usually needs to search k nearest neighborhood (kNN) of each point for building the local tangent plane; for most of point cloud denoising and reconstruction algorithms, the normal vector and the kNN are the indispensable input parameter. Thus, neighborhood searching is of great significance to point cloud processing and its application.

The most commonly used neighborhood in point cloud processing is kNN. It is defined as follows: Given a point cloud p including n points in IR^d and a positive integer $k \leq n - 1$ , compute the kNN of each point in p. More formally, let p = {p₁,p₂,…,p_n} be a point cloud in IR^d, where $d \leq 3$ ; for each $p_{i} \in P$ , let $Nb_k (p_{i})$ be the k points in P, closest to $p_{i}$ .

In recent year, many kNN searching methods have been put forward. The most typical method is to directly calculate the Euclidean distance between current query point $p_{i}$ and any other points one by one, sorting these points in ascending order according to the distance. Then the closest k points are chosen for kNN. This method, although time-consuming and poor-efficiency, is easy to be implemented. To address this problem, many neighborhood searching methods are proposed to improve the computational efficiency. Those algorithms fall into two categories broadly: data reorganization (DR) algorithm and spatial partition (SP) algorithm.⁶

DR^6–12 algorithm divides point cloud into multi-layered subspaces that are employed to construct tree nodes based upon different splitting rules. In these algorithms, the point cloud dataset is regard as multi-children whose tree nodes are different than each other because of splitting rules. Thus, splitting rules can decide searching efficiency. For instance, KD _ tree⁸ is a widely used nearest neighborhood searching method which is a space-partitioning data structure for organizing points in a k-dimensional space. It employs tree nodes to store space range of partition data. Other tree data structure, such as, Cell _ tree and VP _ tree are also employed for nearest neighborhood searching. Cell _ tree⁹ searches nearest neighbors by dividing the data space into cubes of equal size. VP _ tree¹¹ builds a binary tree using a spherical with distance from the selected advantageous location that is different from Cell _ tree building tree employing cubic cell. The KD _tree is most widely used kNN searching method for point cloud.

SP^13–19 algorithm is also called 3D grid partition method. It partitions the bounding box of point cloud into lots of cells. When searching the nearest neighborhood of a point, SP algorithm first locates the cell of a current query point, calculates the distance of each point in the cell to the current query point, and sorts the point in ascending order by distance. If points in the cell are enough, and the k-th shortest distance is smaller than the distance between the current query point, then stop searching. Otherwise, it continues to search more cells until meet conditions. SP algorithm splits the point cloud dataset into a series of cells, and searching the nearest neighborhood in the neighbor cells, thus reducing the searching scope. It is not only searching nearest neighborhood fast, but also with satisfactory accuracy.

Although there are a lot of kNN searching algorithms, most of them mainly focus on how to improve the computational efficiency without considering the accuracy of the neighborhood for point cloud with sharp feature. When the point cloud model is smooth everywhere, kNN can be used to describe local surface information; however, when the point cloud model contains sharp feature, the kNN of feature point are isotropic, which involve points locating at more than one surface patch.²⁰ Point data processing based on the isotropic neighbors may lead to feature blurring or large processing error. To better preserve the original geometry during data processing, such as normal estimation, denoising and surface reconstruction, the neighbors of point should be anisotropic, located at one continuous surface patch. For feature preserving surface mesh denoising, Wang et al.²¹ searched the anisotropic neighborhoods for each vertex by constructing a weighted dual graph, over which biquadratic Bezier surface patches are fitted and projected. Similarly, Zhu et al.²² segmented the isotropic neighborhood to some anisotropic sub-neighborhoods via local spectral clustering. Although these methods can calculate the anisotropic neighborhood for mesh surface, they could not be directly used for point cloud for neighborhood searching. If the neighbors of points searching by SP and DR are directly applied to point cloud with sharp feature for processing, such as normal estimation and point denoising, it will lead to serious processing error. The motivation of the work was to search appropriate neighborhood differing from the most of algorithms which are reducing searching scope and computational time. Thus, an anisotropic neighborhood searching method for point cloud with sharp feature was proposed. In this work, the KD tree was used to search the kNN, and the isotropic kNN are segmented by mapping the neighbors to Gaussian sphere to form a series of sub-neighbors which is anisotropic.

The rest of the work is organized as follows. Section 2 gives the details of our novel algorithm for neighborhood searching. The experimental results and discussions are presented in Section 3, and the conclusions are drawn in Section 4.

Our method

Normal estimation and feature point detection

The normals of points are estimated and used to map the isotropic neighbors into Gaussian image. Many normal estimation methods can be classified into two classes: local surface fitting (LSF) and Delaunay-based methods.²³ LSF²⁴ is commonly used method to estimated normal vector for point cloud. It is suppose that the surface of point cloud is smooth everywhere. As mentioned in our previous work,²⁰ the kNN of point are usually employed for fitting the least square tangent plane and calculating the covariance matrix of local surface. The normal vector of a point is the eigenvector corresponding to the smallest eigenvalue of the covariance matrix.²⁰ LSF method is also called the PCA method. The estimated normals based on PCA method are accurate in flat regions, but normals of feature points are not accurate enough since neighbors of feature point are isotropic. To accurately estimate the normal vector, many algorithms have been reported,^2,23,25–28 but these algorithm are more time-consuming comparing to the PCA method. Compromising efficiency and accuracy, the PCA method is used for normal estimation for point cloud.

The PCA method for normal estimation and feature detection is as follows: For a point set P={p₁,p₂,…,p_n} where n is the total number of a point cloud. Constructing KD tree for point cloud data P, the kNN of point p_i is represented as $Nb_k (p_{i})$ . Let u be the centroid of $Nb_k (p_{i})$ , and the covariance matrix $C_{3 \times 3}$ can be defined as

C_{3 \times 3} = \frac{1}{k} \sum_{i = 1}^{k} (p_{i} - u) (p_{i} - u)^{T} = \frac{1}{k} {[\begin{matrix} p_{1} - u \\ ⋮ \\ p_{k} - u \end{matrix}]}^{T} [\begin{matrix} p_{1} - u \\ ⋮ \\ p_{k} - u \end{matrix}]

(1)

The covariance matrix $C_{3 \times 3}$ can be resolved into three eigenvectors that arranged in descending order of their eigenvalues $λ_{2} \geq λ_{1} \geq λ_{0}$ , where the eigenvalues $λ_{2}$ , $λ_{1}$ and $λ_{0}$ corresponds to the eigenvectors v₂, v₁ and v₀, respectively. v₀ approximates the surface normal n_i for p_i.

The least eigenvalue $λ_{0}$ is to define the surface variation along the surface normal.²⁰ The surface variation of point p_i is described as

ω_{i} = \frac{λ_{0}}{λ_{0} + λ_{1} + λ_{2}}

(2)

The weight $ω_{i}$ measures the confidence of p_i belonging to a feature. If $ω_{i} = 0$ , $p_{i}$ is planar point; the value of $ω_{i}$ is larger; the local surface where p_i located is sharper. If $ω_{i}$ is bigger than a given value Thresh, then p_i can be seen as sharp area points that containing feature point locating at the intersection area of two or more surface patches. The kNN of planar point is located at a continuous surface patch, which is anisotropic. However, when a point is close to feature area, its kNN is from at least two surface patches, which is isotropic. Therefore, the feature area should be detected and kNN of feature point is segmented to obtain anisotropic neighborhood.

Figure 1(a) is a point cloud model with sharp feature, with the red line representing the intersection of two or more surface patch. Point closes to the intersection areas can be seen as feature points. It can be known that kNN of the feature point located at more than two surface patches. The feature point detection result is shown in Figure 1(b), the red points represent feature points. After the feature point cloud detection, the point cloud can be classified to feature points and non-feature points. The kNN of non-feature points are anisotropic and the kNN of feature points are isotropic. The estimated normal vector based on PCA method is shown in Figure 1(c) with the red line representing normal vector. Since the directions of normal vectors are not consistent, they are adjusted to the same direction for displaying. It can be seen from Figure 1(c), normal vector of feature point in intersection area is not perpendicular to the local planar surfaces because kNN used to fit a local plane are from at least two surface patches. In this case, the estimated normal vectors are inaccuracy, and normal vectors of feature points are regarded as smoothed.

Figure 1.

Feature points detection and normal estimating results for point cloud model based on PCA method; (a) point cloud model; (b) feature point cloud estimation result; (c) normalized normal vectors.

The PCA method can be used to estimate normals of points for smooth surface rapidly, that means normals of points are accurate in flat area instead of feature area. To obtain more accurate normal, the normals need to be refined. Lee et al.²⁹ extend the bilateral filtering to refine the normal of triangular faces of the mesh. Given a triangular face f_i with a unit normal n_i and a centroid $c_{i}$ , the refined normal n_i of f_i is defined as:

n_{i}^{'} = \frac{\sum_{j \in N (i)} ω_{d} (‖ c_{i} - c_{j} ‖) W_{n_{1}} [n_{i} \cdot (n_{i} - n_{j})] n_{j}}{\sum_{j \in N (i)} ω_{d} (‖ c_{i} - c_{j} ‖) W_{n_{1}} [n_{i} \cdot (n_{i} - n_{j})]}

(3)

Where $N (i) = {j : ‖ c_{i} - c_{j} ‖ < ρ = 2 σ_{d}}$ is the neighboring face set of $n_{i}$ , and the corresponding unit normal of the neighboring face f_j in $N (i)$ is $n_{j}$ .

Similarly, in our previous work,²³ we proposed an improved bilateral normal filtering to refine the normal of point cloud as follows:

n_{i}^{'} = \frac{\sum_{p_{j} \in Nb_k_{1} (p_{i})} ω_{d} (‖ p_{j} - p_{i} ‖) W_{n_{2}} (‖ n_{j} - n_{i} ‖) n_{j}}{\sum_{p_{j} \in Nb_k_{1} (p_{i})} ω_{d} (‖ p_{j} - p_{i} ‖) W_{n_{2}} (‖ n_{j} - n_{i} ‖)}

(4)

where $Nb_k_{1} (p_{i})$ represent the k₁-nearest neighbors for point $p_{i}$ (the value of k₁ ranges from 16 to 32); $n_{i}$ and $n_{j}$ are the normal vectors of $p_{i}$ and $p_{j}$ respectively, and normals are normalized and adjusted to the same direction; equations (3) and (4) both consider the same Gaussian function $ω_{d} (d)$ , and the difference lies on the definition of $W_{n_{1}}$ and $W_{n_{2}}$ in equations (3) and (4). Therefore, we analyze $W_{n_{1}}$ and $W_{n_{2}}$ , that is,

{\begin{matrix} W_{n_{1}} (n) = \exp (- {| n_{i} \cdot (n_{i} - n_{j}) |}^{2} / 2 σ_{n}^{2}) \\ W_{n_{2}} (n) = \exp (- {(‖ n_{j} - n_{i} ‖)}^{2} / 2 σ_{n}^{2}) \end{matrix}

(5)

$W_{n_{1}}$ and $W_{n_{2}}$ are the Gaussian functions of the normal vector in equations (3) and (4) respectively, and the $δ_{n}$ is the standard deviation. Consider the unit normal vector $n_{i}$ as the y-axis direction under xoy coordinate system, and the distribution of $W_{n_{1}}$ and $W_{n_{2}}$ in term of $n_{j}$ are shown in Figure 2. From Figure 3 it can be know that, $W_{n_{1}}$ is greater than $W_{n_{2}}$ with the same $n_{j}$ , which means that $n_{j}$ under $W_{n_{1}}$ has more weights on the refining result $n_{i}^{'}$ of $n_{i}$ than that under $W_{n_{2}}$ . As a result, equation (3) is more sensitive to noise and feature point compared with equation (4). In addition, the weight that is, $W_{n_{1}}$ and $W_{n_{2}}$ both decrease as the normal vector difference of $n_{i}$ and $n_{j}$ increases, therefore, a larger the deviation between $n_{j}$ and $n_{j}$ , the less effect on the refining result.

Figure 2.

The distribution of three different weights.

Figure 3.

Schematic of point cloud and normal vector; (a) point cloud in two surface patches; (b) normal vector before refining.

Figure 3 is the schematic of point cloud and normal vector before refining in two surface patches, where the red point represents the feature point. The kNN of feature points are located at two surface patches, and the blue points in surface patch S₂ can be regard as noise or outlier for point in surface S₁ when normal vector refining. Although the kNN of feature point locate at two surface patches at least, most of neighbors of feature points are non-feature points. Normal vectors of non-feature points are similar to each other because some of neighbor points are located at one surface patch. When neighbor points are located at different surface patches, normals of these points are deviated from each other; the surface is sharper; the angle of normal vector between different points is larger.

Even though equation (4) decreases the effect of noise point, it still exists with a certain value. Consequently, more iteration times are required to refine the normal vector. To tackle this problem, we optimize equation (4) as follows:

n_{i}^{'} = \frac{\sum_{p_{j} \in Nb_k_{1} (p_{i})} ω_{d} (‖ p_{j} - p_{i} ‖) W_{n_{3}} (n_{j}, n_{i}) n_{j}}{\sum_{p_{j} \in Nb_k_{1} (p_{i})} ω_{d} (‖ p_{j} - p_{i} ‖) W_{n_{3}} (n_{j}, n_{i})}

(6)

where $W_{n_{3}} (n)$ is represented by:

W_{n_{3}} (n_{j}, n_{i}) = {\begin{matrix} 0 if (n_{j}, n_{i}) n_{j} \geq u_{1} \\ {[(n_{j}, n_{i}) n_{j} - u_{1}]}^{2} otherwise \end{matrix}

(7)

where $u_{1} = \frac{\sqrt{\sum_{p_{j} \in Nb_k_{1} (p_{i})} {| (n_{j}, n_{i}) n_{j} |}^{2}}}{k_{1}}$ .

Essentially, we truncate the normal vectors if the difference between them and $n_{i}$ are greater than the average normal vector difference $u_{1}$ . From Figure 2, the weight $W_{n_{3}}$ is always smaller than $W_{n_{1}} (n)$ and $W_{n_{2}} (n)$ with the same $n_{j}$ . Thus, our method is less sensitive to noise and outlier.

In order to improve the computational efficiency, only normals of feature points are refined. Figure 4 is the normal refine results for two different point cloud model. It can be seen that the refined normal vector is nearly perpendicular to the corresponding local surface. For comparison, the refined normal based on our method and our previous work²³ is shown in Figure 5. Three iterations are performed for both methods when normal refining. The detailed comparison result can be seen the black rectangular area. It can be seen that our method can obtain more precision results for normal refining with the same iteration.

Figure 4.

Normal refining result for point cloud; (a), (b) are the point cloud models; (c), (d) are normal estimation results based on PCA method; (e), (f) are normal refining results.

Figure 5.

The comparison for normal refining result based on our method and our previous work.²³ (a) our previous work²³; (b) our method.

Anisotropic neighborhood seaching

To get the anisotropic neighborhood, the isotropic neighborhood is divided into several sub- neighborhoods by mapping the isotropic kNN into Gaussian image.

Let $T = T (u, v) \subset R^{3}$ be a regular surface with consistent and normalized normal vectors in a three dimensions space. The map $G : T \to T^{2}$ , and $T^{2}$ represents a unit normal of a point p ( $p \in T$ ) in the unit sphere,²⁰ where

T^{2} = {(x, y, z) \subset R^{3} | x^{2} + y^{2} + z^{2} = 1}

The resulting map $G : T \to T^{2}$ is defined as the Gaussian map of the surface T. The Gaussian image can be denoted by $G (T) = {n (p) | p \in T}$ . The unit sphere $T^{2}$ can be called the Gaussian sphere.³⁰

The Gaussian image has very important property: If a connected surface is a plane, the Gaussian image of the point on that surface is a point; Conversely, if the Gaussian image of a connected surface is a point, all the points on that surface are planar points.²⁰ However, there is always error between estimated normal and ground-truth normal. Therefore, the Gaussian image of point cloud dataset in a planar surface display a cluster where data is in dense distribution. Figure 6(a) shows the sketch map of Gaussian image of planar surface

Figure 6.

Schmatic plot of Gaussian map; (a) Gaussian map of planar surface; (b) Gaussian map of dihedral surface patch.

If the estimated normal vector is approximated to the ground-truth normal, for a point locating at the surface edge area, its kNN points may be from two surface, that the kNN points in Gaussian map displaying two dense clusters; And for a point locating at the surface corner area, its Gaussian map displays several clusters, and the number of clusters equals to that of surface patches. Figure 6(b) shows the schematic plot of Gaussian image of dihedral surface, where the red and black clusters correspond to points from two different surfaces patches. Therefore, the clustering algorithm is necessary to separate the data into appropriate cluster on Gaussian sphere.

To segment the isotropic neighborhood of feature points p_i, a larger neighbors $Nb_k_{2} (p_{i})$ are selected. For any point p_i, its neighbors with size k₂ are represented as p_ij; $n_{i}$ is the normalized normal of point p_i; $n_{ij}$ the neighborhood normal. All normal vectors must be normalized and adjusted to the same direction for Gaussian mapping. Then, the Gaussian map of a point p_i can be denoted by $G (p_{i}) = {n_{ij}}$ , j=0,1,…,k₂-1. The clustering algorithm is employed to divide the dataset in Gaussian sphere to appropriate sub-cluster. There are many clustering algorithms have been report in literature, but most of them require knowing the number of clusters to be produced in advance. However, the number of surface patches contained in the isotropic kNN of feature point is unknown, so is the number of clusters in the Gaussian map.

Mean shift algorithm and hierarchical agglomerative algorithm are non-parameter clustering algorithms to cluster data set without specifying the number of cluster, widely used in data mining and machine vision area. However, it is found that the mean shift algorithm usually classifies one cluster data with too many sub-clusters, some of which may just contain one or two points. The hierarchical agglomerative (”bottom-up”) clustering method³¹ was employed to classify data in the work. It begins with each point as a separate cluster, and merges them successively into larger clusters. For the hierarchical agglomerative clustering algorithm, it is an important step to determine the similarity between two clusters.

There are three classic standards, which are single-linkage, complete-linkage and average-linkage measures, to define the distance between two clusters. The average-linkage standard defining the mean distance D_c between elements of each cluster tends to merge two clusters with small difference. It is used to define the distance between two clusters in the work.

The data set $G (p_{i})$ is the input data for hierarchical agglomerative algorithm. Then, the mean distance D_c between elements of each cluster can be defined as

D_{c} (C_{1}, C_{2}) = \frac{1}{| C_{1} | | C_{2} |} \sum_{n_{ix} \in C_{1}} \sum_{n_{iy} \in C_{2}} d (n_{ix}, n_{iy})

(8)

where $| C_{1} |$ and $| C_{2} |$ are the numbers of point data corresponding to clusters C₁ and C₂, respectively; $d (n_{ix}, n_{iy})$ is the Euclidean distance between two points $n_{ix}$ and $n_{iy}$ in the Gaussian map, where x, y = 0,1,…,k₂-1.

d (n_{i}, n_{j}) = 2 (1 - 〈 n_{ix}, n_{iy} 〉)

(9)

where $〈 n_{ix}, n_{iy} 〉$ represents the inner product of $n_{ix}, n_{iy}$ , when the angle between $n_{ix}, n_{iy}$ is less than T. That is, when the distance of two points in Gaussian image is less than $2 (1 - \cos θ)$ , two points are merged into the one cluster. The largest cluster corresponding to the points are selected as the optimal neighbor for the feature point.

Results and discussion

Parameter selection

There are five parameters influencing the neighborhood searching result, requiring the user to prepare in the algorithm. They are Thresh, t, k, k₂, and $θ$ , respectively.

Feature point detection threshold Thresh requires the user to prepare for selecting the feature point where neighbors should be segmented, and normal needs to be modified. The sharpness threshold is related to the “sharpness” of the features of the point cloud. If the surface of point cloud contains very sharp features, it should be a small value; otherwise, it is assigned with a relatively high value. Based on the experimental, the value of Thresh can be 0.001 to 0.1.

t is the iteration time to refining the normal vector. We have the setting in our implementation: t: [1–5].

k is the size of nearest neighborhood, representing the KD tree neighborhood of point data for user to define. After establishing KD tree for point cloud, the kNN of the point cloud are selected by user. When a point is non-feature, the kNN is the final neighbor.

k ₂ is the size of nearest neighborhood of feature point, with neighbors mapped into Gaussian image. Usually, the k₂-nearest neighbors are isotropic, segmented into several sub-neighbors that are anisotropic. The value of k₂ is associated with k. Since the feature point is close to the area of intersection of two or more surface patches, we need to consider the surface complexity to set the value of k₂, thus making the neighbor size of feature point close to the neighbor of regular point as much as possible,.

For example, if the feature points are the intersection of two or more surface patches, it can be set as k₂ = 2k; if the feature points are the intersection of three or more surface patches, then it can set as k₂=3k. It is assumed k₂ = 3k, and k = 16 in default in this work.

Parameter $θ$ is used in hierarchical clustering algorithm that denotes the angle threshold of two normal vectors. That is, when the distance of two points in Gaussian is less than $2 (1 - \cos θ)$ , two points are merged into the one cluster. When points are located at one continuous surface patch, the normal vector of point between each other is nearly parallel to each other in local area, and the angle between normal vectors is small. However, when the distance between two points in Gaussian sphere is large, these points are from two different surface patches.

The sharper surface patch leads to the larger angle between normal vectors of points in different surface patch. $θ$ is an important threshold, if $θ$ is set to a big value, the segmented neighbors is still isotropic; if set to a small value, the isotropic neighbors are segmented into many anisotropic sub-neighbors. Previous researches^25,26,32 have measured the accuracy of estimated normal vector with angle threshold $θ = 10 °$ . In other words, when the angle between the estimated normal vector and ground-truth normal vector is less than $10 °$ , the estimated normal is regard as accurate, and the corresponding point is a flat point; otherwise, the point is non-flat or bad. The threshold $θ$ is similar to normal accuracy measurement, so we set $θ = 10 °$ according to the suggestion in the research.^25,26,32 Namely when the distance between two points in Gaussian image is less than 0.173, two points or clusters are merged into one large cluster.

Applications of point cloud neighborhood

In order to verify the validity of our algorithm, the anisotropic neighbors of point cloud are applied to point data processing. We compare our proposed method with the state-of-the-art neighborhood searching algorithms, such as KD tree¹⁴ and 3D grid¹³ kNN searching method. These algorithms are implemented using Microsoft vision studio 2017 C++, and OpenGL. Since it is hard to visualize neighbors for every point, the neighbors of point are applied for data processing, for example, normal estimation based on local surface fitting, and point data denoising. In the comparison, all the parameters are same except for neighbors.

Normal estimation based on local surface fitting

Normal vector is an important property of a 3D surface. Many effective point cloud processing algorithms depend on the accurate normal as input to produce a precise processing result. The normal estimation can be applied in point cloud denoising, simplification, registration and surface reconstruction. LSF method estimates the normal vector by fitting the least square tangent plane using the local neighborhoods. If neighborhoods of point are anisotropic, the estimated normal is nearly perpendicular to the local tangent plane which means normal are accurate. LSF method is used to estimated normal vector with the KD tree, 3D grid and our anisotropic neighbors as input.

Figures 7 and 8 shows the visualization results of normal vectors based on different neighborhoods as input. Octahedral and fandisk model consists of many continuous surface patches, and feature boundaries can be observed clearly. Normals are normalized and adjusted to the consistent direction for displaying. It can be known that the normal vector of point is nearly perpendicular to the local surface in flat area for 3D grid and KD tree neighbor as input. However, the normal vector of feature point is smoothed, deviated from the local surface because the neighbors of feature point are from at least two surface patches. The estimated normal vectors based on our neighbors are accurate in flat area and feature points because neighbors of feature points are anisotropic, only located at one continuous surface patch.

Figure 7.

Normal vector estimation results for octahedral model based on different neighbors as input: (a) Triangulated model of octahedral point cloud; (b) KD tree neighbor; (c) 3D grid neighbor; (d) Our neighbor.

Figure 8.

Normal vector estimation results for findisk model based on different neighbors as input: (a) Triangulated model of fandisk point cloud; (b) KD tree neighbor; (c) 3D grid neighbor; (d) Our neighbor.

Point cloud denoising

With the development of 3D scanner, it is very convenient to acquire 3D point cloud. However, the collected point cloud data inevitably contains noise due to a variety of internal and external factors interference. Denoising is necessary step to improve the data quality for further data processing. The main problem of point cloud denoising is to remove noise while preserving sharp feature.

Bilateral filtering algorithm³³ can be used to denoise for point cloud and mesh surface while preserving the sharp feature; however, it needs the accurate normal and neighbors as input. Bilateral filtering denosing for point cloud is given as follows: Let $p \in P$ be a point of input point data set P, the bilateral filter denoising result for p is $p^{'}$ .

p^{'} = p + α n

(10)

where $α$ is the filter coefficient.

α = \frac{\sum_{p_{j} \in Nb_K (p)} W_{c} (‖ p - p_{j} ‖) W_{s} (‖ 〈 p - p_{j}, n 〉 ‖) 〈 p - p_{j}, n 〉}{\sum_{p_{j} \in Nb_K (p)} W_{c} (‖ p - p_{j} ‖) W_{s} (‖ 〈 p - p_{j}, n 〉 ‖)}

(11)

Where n is the normal vector of point p; $〈 p - p_{j}, n 〉$ is the inner product of $p - p_{j}$ and n; $Nb_K (p)$ is the K neighbors of point p, with $p_{j} \in Nb_K (p)$ ; $W_{c}$ and $W_{s}$ are weighting functions, where $W_{c} (x) = \exp (- x^{2} / 2 σ_{c}^{2})$ ; $W_{c} (y) = \exp (- y^{2} / 2 σ_{s}^{2})$ ; $σ_{c} = σ_{s} = 0.2$ .

To highlight the influence of neighborhood on denoising results, normal of point is estimated by method introduced in the Literature.²⁶ It can accurately estimate normal vector and robust to noise. Two point cloud which contain sharp feature are employed to test denoising results.

Figures 9 and 10 show the denoising results based on 3D grid, KD tree and our neighbors as input $Nb_K (p)$ . In order to clearly observe the denoising results, the point cloud models are triangulated for displaying. The deviation including maximum error (ME), average error (AE) between denoising model and noise-free model are measured using Geomagic Studio software.³⁴ The values of ME and AE for mounting-foot and gear model are shown in Table 1, and the deviation results are shown in Figures 9 and 10. Figures 9 and 10 and Table 1 show that our neighbors as input for denoising can preserve sharp feature well and generate the least ME and AE values.

Figure 9.

Bilateral filter denoising results for mounting-foot model with KD tree, 3D grid and our neighbors as input; (a) mounting-foot point cloud model with noise; (b) triangulation of Figure 9(a); (c) triangulation of mounting-foot model with noise-free; (d)–(f) Denoising and deviation results based on KD tree, 3D grid and our method, respectively.

Figure 10.

Bilateral filter denoising results for gear model with KD tree, 3D grid and our neighbors as input; (a) gear point cloud model with noise; (b) Triangulation of Figure 10(a); (c) Triangulation of gear model with noise-free; (d)–(f) Denoising and deviation results based on KD tree, 3D grid and our method, respectively.

Table 1.

Comparison of error (mm).

Point model	Error	KD tree	3D grid	Our method
Mount-foot	AE	0.025	0.020	0.007
Mount-foot	ME	2.415	2.430	1.362
Gear	AE	0.002	0.002	0.001
Gear	ME	0.159	0.153	0.072

Computational time analysis

Our proposed method needs to fit a local tangent plane to estimate normal vector and detect feature points, and normal vectors of these points are iteratively refined. After hierarchical clustering in Gaussian sphere, isotropic neighbors are classified into several anisotropic sub-clusters besides constructing KD tree. Therefore, it is more computational than the neighborhood searching method by 3D grid and KD tree. Test computations are implemented on a PC with a 2.5 GHz Intel Core processor and 16 GB RAM. Table 2 shows the results of computational time for different point cloud. Different computers take different amounts of time to execute the algorithm, thus, the computational time in Table 2 just for reference. According to the computational complexity and Table 2, it is obvious that our method is more time-consuming that several times larger than that required by the method of 3D grid and KD tree.

Table 2.

Comparison of computational time.

Pointmodel	Number ofpoint cloud	Computationaltime (s)
		KD tree	3D grid	Ourmethod
Fandisk	11230	0.21	0.19	0.41
Octahedron	32219	0.58	0.51	1.06
Mount-foot	110985	1.92	1.85	4.07
Gear	105309	1.88	1.78	3.89

Conclusion

The neighborhood searching of point data is very import and a basic work for point cloud processing. A large amount of point data processing algorithms need neighbors as input to produce precise processing result, such as normal estimation, denoising, and reconstruction. Since the kNN of feature points is isotropic that may contain points belonging to more than two surface patches, sharp feature of point cloud model are over-smoothed with these neighbors as input for point processing. An anisotropic neighborhood searching method was proposed in the work.

Our method divides a big isotropic kNN of feature points into several anisotropic sub-neighborhoods by mapping the big isotropic kNN neighbors into a Gaussian sphere. The points on Gaussian image formed several sub-clusters corresponding to the anisotropic sub-neighbors. The agglomerative hierarchical clustering algorithm was used to separate the data in Gaussian image into appropriate cluster. The biggest sub-neighborhood corresponding to the biggest sub-cluster was selected as the optimal anisotropic neighbors. Since a large kNN of feature point was segmented into a series of sub-neighbors, the number of neighbors of point is not same everywhere. Neighbors of points are applied to estimate normal and bilateral filtering denoising.

K nearest neighborhood usually is isotropic for point cloud with sharp feature, different from the other k nearest neighborhood search methods, this paper emphasizes on the anisotropic neighborhood searching instead of computational efficiency. Point cloud normal estimation and denoising based on 3D grid, KD tree and our neighbors as input, our method can obtain more precision result for point cloud with sharp feature. The anisotropic neighborhood searching based on our method can be used for point cloud normal estimation, denoising, surface fitting, and reconstruction, et al. It is important to point out that the method in the work is more computational and complex than DR and SP neighborhood searching method. The future work will improve the computational efficiency and reduce complexity for searing anisotropic neighborhood.

Footnotes

Acknowledgements

The authors would like to thank Jiangxi Province Education Department, Jiangxi provincial department of science and technology, National Natural Science Foundation of China and our previous work²⁰ to support this work.

Author contributions

Xiaocui Yuan proposed the idea of this algorithm, and wrote programs to implementation our algorithm and other correlation algorithm. Baoling Liu designed the framework of this paper. Yongli Ma was responsible for data collection and modified paper.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is support by Jiangxi Province Education Department (Grant No. GJJ161122, GJJ161104, GJJ16110), Jiangxi Provincial Department of Science and Technology (Grant No. 20171BAB206037) and National Natural Science Foundation of China (Grant No. 61903176, 62001202).

ORCID iD

Xiaocui Yuan

References

Arvanitis

Lalos

Moustakas

. Denoising of dynamic 3D meshes via low-rank spectral analysis. Comput Graph 2019; 82: 140–151.

Cao

Chen

Zhang

, et al. Normal estimation via shifted neighborhood for point cloud. J Comput Appl Math 2018; 329: 57–67.

Wang

. Iterative closest point registration for fast point feature histogram features of a volume density optimization algorithm. Meas Contr 2020; 53(1–2): 29–39.

Xing

, et al. A dynamic and adaptive scheme for feature-preserving mesh denoising. Graph Model 2020; 110.

Gao

Deng

Qin

, et al. Point cloud and 3-d surface reconstruction using cylindrical millimeter-wave holography. IEEE Trans Instrum Meas 2019; 68: 4765–4778.

Ding

, et al. A new extracting algorithm of k nearest neighbors searching for point clouds. Pattern Recognit Lett 2014; 49: 162–170.

Yen

Shih

, et al. Applying multiple kd-trees in high dimensional nearest neighbor searching. Int J Circ Syst Signal Process 2010, 4(4): 153–160.

Bentley

. K-D trees for semidynamic point sets. In: Proceedings of the sixth annual symposium on computational geometry, Berkley, CA, May 1990, pp.187–197. New York, NY: Association for Computing Machinery.

Clarkson

. Fast algorithms for the all nearest neighbors problem. In: 24th annual symposium on foundations of computer science, Tucson, AZ, 7–9 November 1983, pp.226–232. New York: IEEE.

10.

Arya

Mount

Netanyahu

, et al. An optimal algorithm for approximate nearest neighbor searching. In: SODA ’94: proceedings of the fifth annual ACM-SIAM symposium on Discrete algorithms, Arlington, VA, USA, 23–25 January 1994, pp.573–582. New York: ACM.

11.

Yianilos

. Data structures and algorithms for nearest neighbor search in general metric spaces. In: SODA ’93: proceedings of the fourth ACM-SIAM symposium on discrete algorithms, Austin, TX, 25–27 January 1993, pp.311–321. Philadelphia, PA: Society for Industrial and Applied Mathematics.

12.

Behley

Steinhage

Cremers

. Efficient radius neighbor search in three-dimensional point clouds. In: IEEE international conference on robotics and automation, Seattle, WA, 26–30 May 2015, pp.3625–3630. New York: IEEE.

13.

Yang

Lin

Zhang

, et al. A fast algorithm for finding k-nearest neighbors of large-scale scattered point cloud. Int J Image Graph 2016; 18(4): 399–406.

14.

Piegl

Tiller

. Algorithm for finding all k nearest neighbors. Comput Aided Des 2002; 34: 167–172.

15.

Wei

Zhang

Zhou

. Spatial sphere algorithm for searching k-nearest neighbors of massive scattered points. Acta Aeronaut. Astronaut. Sin 2006; 27: 944–948.

16.

Ping

. Novel algorithm for k-nearest neighbors of massive data based on spacial partition. J South China Univ Technol 2007; 35: 65–69.

17.

Xiong

. Algorithm for finding ~k~-nearest neighbors of scattered points in three dimensions. J Comput Aided Des Comput Graph 2004; 16: 909–908.

18.

Fang

Zhao

, et al. Algorithm for finding k-nearest neighbors based on spatial sub-cubes and dynamic sphere. Geomatics Inf. Sci. Wuhan Univ 2011; 36: 358–362.

19.

Min

Yun

. Study on point cloud management strategy based on Octree-like index. Laser & Optoelectronics Progress, 2020, http://kns.cnki.net/kcms/detail/31.1690.TN.20200114.0735.004.html

20.

Yuan

Peng

, et al. A novel method of normal estimation for 3D surface reconstruction. In: ASME 2015 international design engineering technical conferences and computers and information in engineering conference, 2–5 August 2015, Boston, MA, p.9. New York: American Society of Mechanical Engineers.

21.

Wang

Zhang

. A cascaded approach for feature-preserving surface mesh denoising. Comput Aided Des 2012; 44: 597–610.

22.

Lei

Wei

, et al. Coarse-to-fine normal filtering for feature-preserving mesh denoising based on isotropic subneighborhoods. Comput Graph Forum 2013; 32: 371–380.

23.

Yuan

Lu-Shen

Chen

. Normal estimation of scattered point cloud with sharp feature. Optic Precis Eng 2016; 34(10): 2581–2588.

24.

Hoppe

DeRose

Duchamp

, et al. Surface reconstruction from unorganized points. ACM SIGGRAPH Comput Graph 1992; 26(2): 71–78.

25.

Zhang

Cao

Liu

, et al. Point cloud normal estimation via low-rank subspace clustering. Comput Graph 2013; 37: 697–706.

26.

Wang

Feng

H-Y

Delorme

F-É

, et al. An adaptive normal estimation method for scanned point clouds with sharp features. Comput Aided Des 2013; 45: 1333–1348.

27.

Liu

Yuan

Zhao

. 3D point cloud denoising and normal estimation for 3D surface reconstruction. In: IEEE international conference on robotics and biomimetics, Zhuhai, 6–9 December 2015, pp.820–825. New York: IEEE.

28.

Sun

Wang

. Deep feature-preserving normal estimation for point cloud filtering. Comput Aided Des 2020; 125: 102860.

29.

Lee

Wang

. Feature-preserving mesh denoising via bilateral normal filtering. In: Proceedings of the ninth international conference on computer aided design and computer graphics, Hong Kong, 7–10 December 2005, pp.275–280. New York: IEEE.

30.

Wang

Hao

Ning

, et al. Automatic segmentation of urban point clouds based on the Gaussian map. Photogramm Rec 2013; 28: 342–361.

31.

Hastie

Tibshirani

Friedman

, et al. The elements of statistical learning: data mining, inference and prediction. Math Intel 2005; 27: 83–85.

32.

Schnabel

Klein

, et al. Robust normal estimation for point clouds with sharp features. Comput Graph 2010; 34: 94–106.

33.

Fleishman

Drori

Cohen-Or

. Bilateral mesh denoising. ACM Trans Graph 2003; 22(3): 950–953.

34.

https://www.3dsystems.com/press-releases/geomagic/announces-studio-2013.

Anisotropic neighborhood searching for point cloud with sharp feature

Abstract

Keywords

Introduction

Our method

Normal estimation and feature point detection

Anisotropic neighborhood seaching

Results and discussion

Parameter selection

Applications of point cloud neighborhood

Normal estimation based on local surface fitting

Point cloud denoising

Computational time analysis

Conclusion

Footnotes

Acknowledgements

Author contributions

Declaration of Conflicting Interests

Funding

ORCID iD

References