Sharp Feature Detection as a Useful Tool in Smart Manufacturing

: Industry 4.0 comprises a wide spectrum of developmental processes within the management of manufacturing and chain production. Presently, there is a huge effort to automate manufacturing and have automatic control of the production. This intention leads to the increased need for high-quality methods for digitization and object reconstruction, especially in the area of reverse engineering. Commonly used scanning software based on well-known algorithms can correctly process smooth objects. Nevertheless, they are usually not applicable for complex-shaped models with sharp features. The number of the points on the edges is extremely limited due to the principle of laser scanning and sometimes also low scanning resolution. Therefore, a correct edge reconstruction problem occurs. The same problem appears in many other laser scanning applications, i.e., in the representation of the buildings from airborne laser scans for 3D city models. We focus on a method for preservation and reconstruction of sharp features. We provide a detailed description of all three key steps: point cloud segmentation, edge detection, and correct B-spline edge representation. The feature detection algorithm is based on the conventional region-growing method and we derive the optimal input value of curvature threshold using logarithmic least square regression. Subsequent edge representation stands on the iterative algorithm of B-spline approximation where we compute the weighted asymmetric error using the golden ratio. The series of examples indicates that our method gives better or comparable results to other methods.


Introduction
Industry 4.0 is the trend towards automation and data exchange in manufacturing technologies and processes. Industry 4.0 factories have machines which are augmented with wireless connectivity and sensors, connected to a system that can visualize the entire production line and make decisions on its own. Smart manufacturing is a broad category of manufacturing that employs computer-integrated manufacturing, high levels of adaptability and rapid design changes, digital information technology, and more flexible technical workforce training [1,2].
Our work covers the part of automation in manufacturing and reverse engineering for design changes. The widespread availability of cheap commercial depth sensors or multi-camera setups leads to their common usage in manufacturing engineering, especially in rapid product development [3]. Therefore, the point cloud has become the regular type of the representation used e.g., in reverse Usually, 3D scanners are provided with software that creates a physical object in the two-phase setup: reconstructing CAD model from the scanned point data, and subsequent STL model output. However, the implemented methods are usually sufficient only for simple and smooth objects. More complex objects with holes, sharp edges or corners often require prior knowledge of how to edit the scanned data. The data post-processing is time-consuming and can potentially lead to errors. We present an algorithm that automatically detects sharp features in arbitrary point cloud data. Our algorithm is based on a modified region-growing method. The input value to the algorithm is only one parameter-the normal vector threshold θ N (optimal start value and the description is presented in Section 2.1). Region-growing algorithm eliminates the planar points and the remaining points are probably edge points.
Before the final manufacturing process or 3D printing, a correct model is essential for STL generation [5] and subsequent slicing [4,6]. Figure 2 shows the common problematic part (jagged and inaccurate borders) that occurs frequently. To solve this problem, in Section 3.2 we propose the novel edge representation based on B-spline approximation. The asymmetric error evaluation is based on bell function with weighting using golden ratio [7] to improve the edge representation of the models.
Moreover, our B-spline approximation is also sufficient for direct STL slicing from point clouds that extracts the sectional contours (curves) directly without model reconstruction. The overview of these methods is e.g., in [8]. The presented B-spline approximation can be used without limitations in this slicing procedure because it obviously iteratively minimizes the shape-error of the layer curves.

Prior Work
Point cloud segmentation and classification plays a key role in point cloud processing in RE/RP usage. Segmentation is the process of grouping point clouds into multiple homogeneous regions with similar properties whereas classification is the step that labels these regions. Segmentation methods are divided into many categories and the choice of the methods depends on the area of interest, sampling density, data quality and explicit structure of point cloud. We can use the categorization described in [9,10]: region-growing, edge-based methods, model fitting, and hybrid methods.
Most of the methods start with region growing that was purposed by Besl [11] for segmentation of images. The region-growing algorithm was later expanded into 3D in two different versions: unseeded [12] and seeded [13]. Region growing is applied in plenty of different applications, e.g., extracting features and planar surfaces from 3D outcrop point clouds [14], surface segmentation, and edge feature lines extraction from fractured fragments of relics [15] or urban environment modelling [16] and also 3D city modelling (literature survey e.g., in [17]).
Region growing is frequently improved. For example, Demaris et al. [18] detect the closed sharp edges using the region growing with normals given by building a local mesh and least-squares planes through the nearest neighborhood points. The correct edge detection in provided by graph theory (minimum spanning tree). In [16], region-growing step is performed on an octree-based voxelized input point cloud representation to extract major segments.
The edge-based methods are usually used with point clouds transformed into the range images [19]. Using the high-level features (scan line grouping technique) as segmentation primitives instead of individual pixels ensure the high computational speed. However, this edge detection is not applicable for noisy data.
Model fitting methods assume that a surface can be described as a composition of simple, canonical geometric shapes (watertight reconstruction). Elementary methods use RANSAC to robustly find planes, spheres, cylinders, cones, and torii [20]. However, there occurs a problem with noise, any fine-grained details are likely to be treated as noise with the primitive prior, if they are unable to be represented as a union of smaller primitives (a survey is described in [21]). The partial solution is hybrid methods that employ primitives with other prior [22].
Region-growing, edge-based methods and model fitting do not need pre-processing, but there exist other segmentation methods based on pre-processing. These algorithms are capable of reconstructing data incompleteness [23,24] or construct a triangular mesh of given point cloud [25][26][27]. Afterwards, the formed mesh can be processed using e.g., the discrete differential geometry where local extreme of the surface detect the edge lines [25] or algebraic methods based on bivariate polynomials [27]. We must note that some of these methods work with the machine learning methods [28].
After the point cloud segmentation and classification, we need a correct model representation and tessellation (usually STL) with subsequent slicing to extract layer-based additive manufacturing [29]. There exist various approaches that combine different technologies. Traditional STL approach is described e.g., in [30]. B-spline or NURBS surfaces are used frequently, e.g., authors in [31,32] presents direct slicing from NURBS without STL. Also, B-Rep model composed of planes, spheres, cylinders and cones from a 3D mesh is fitting for slicing [33,34]. Some of the authors omit the model reconstruction and make the direct generation of sectional contour. For example, authors in [35] employ the curve skeleton of the model to slice through the surface mesh edges. B-spline curve fitting of the layers where the curve curvature is determined by a circle fitting procedure is presented in [8].
The development of neural networks and deep learning also touches the research in point cloud analysis. Some methods are promising but there is still no one-size-fits-all approach. Reviews of the methods are covered e.g., in [36][37][38].
Present literature contains most of the recently known principles. Main scientific contribution is in their improvements and applications. For example, work [39] (published 2020) uses B-spline representation and normal computation using normal curvature directly computed by partial derivatives of the surface. The omitting of the pre-defined threshold value is the main idea in [40] (published 2019). The authors suggest the composition of spatial FFT-based filtering and boundary detection, which together allow for direct generation of low noise tessellated surfaces from point cloud data. Ref [41] presents the feature sensitive point cloud simplification that is based on insensitive support vector regression.
The advantages of the well-known methods (region-growing or PCA) are obvious-robustness, computed easiness, stable algorithms. Because of that, we chose these classical methods and our main idea was (1) the simplification and automatic setting of input values (described in Section 2.1) (2) find optimal error evaluation in B-spline representation that is suitable for engineering objects. The problem of error computation is described in Section 2.2 and we proposed asymmetric weighting based on golden ration.

Methods and Materials
The common industrial engineering components consist of flat, smooth or curved surfaces. A measurable difference in surface normals or curvature is relatively small on most of the model surfaces. Contrary to that, in the neighborhood of a sharp feature, such as a sharp edge or corner, there is a significant change in the surface normal and curvature. Thus, if we remove all points of the flat, smooth or curved surfaces, the remaining ones are the points on the sharp features. With this idea in mind, we build our method, similarly to [42], on the analysis of the eigenvalues of the covariance matrix for every point's k-nearest neighbor. We use region growing not to find planar surfaces or regular models [43] but to remove the points with low curvature (Section 2.1).
The important input values of region growing are normal and curvature thresholds. Their setting substantially changes the segmentation results, but their choice is unintuitive and difficult especially for beginners. We find the geometric correlation between these two values and we compute the optimal initial value that is applicable in the case of regular engineering objects. This is described in the following subsections. This approach is robust to the noise and is suitable for subsequent B-spline approximation.
The B-spline approximation was chosen due to its good properties and well-known computation. The iterative process ensures sufficient quality. The error computation is based on the asymmetrically weighted distance (asymmetric distance is presented in [44]) but in our case, the weight is described with the golden ratio. The asymmetric measurement replaces the outliers and protects the curve from the zig-zag effect. This approach is described in Section 2.2. Computed B-spline curve can represent edges in point clouds. Moreover, described B-spline approximation is also possible to use in the methods for direct slicing.

Edge Points Detection
The processing of a given point cloud starts with the region-growing segmentation. Surface normal computation is based on Principal Component Analysis (PCA) [45] that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables. The description is in Section 2.1.1.
Edges can be also classified by Gauss mapping [42,46,47]. The point normal vectors are projected on the Gauss sphere and subsequent clustering distinguish the point type (corner, edge, planar). This method is time-consuming so that we use the region growing.

PCA-Normal Estimation
The first step in the PCA algorithm is the computation of the covariance matrix C for each point P in the point cloud.
whereP is the centroid of k-nearest neighborhood points P i , i = 1, 2, . . . k of a point P. Eigenvalues (λ 0 , λ 1 , λ 2 ) and eigenvectors ( v 1 , v 2 , v 3 ) of this matrix defines the covariance ellipsoid. We distinguish between these three cases: These types of covariance ellipsoids define the neighborhood of the given point (Figure 3a-c). In the case of an oblate ellipsoid ((λ 0 ≤ λ 1 ≤ λ 2 ) ∧ (λ 1 ≈ λ 2 )), the point is obviously on an almost planar area and the eigenvector corresponding with the smallest eigenvalue is point normal.

Region Growing-Initial Normal Threshold Value thr θ N
Region-growing segmentation algorithm usually requires input parameters: minimal and maximal cluster size, k-nearest neighbor points, normal threshold thr θ N and curvature threshold thr θ C . Value θ N is the angle between two-point normals and θ C is curvature. They strongly affect the segmentation results and geometrically, they are tightly connected. The setting of the initial values is unintuitive and not user-friendly so that we focus on the optimal automatic estimation of θ N . Figure 4a shows the definition of the curvature θ of the curve l at the point P 1 . In the limit as ∆t → 0, we obtain: In the case of a point cloud, we consider the parameter ∆t as the minimal distance d min of the line between two adjacent points ( Figure 4b). The curvature θ C is computed by Therefore, the region-growing segmentation can be controlled only by one parameter-the normal vector threshold thr θ N . The following part describes the process of the ideal initial value of thr θ N .
We prepared a set of testing point clouds M i , i = 1, 2, . . . , 6 to cover the spectrum of different shapes from simple one to the complex ones ( Figure 5). These point cloud data were acquired by ATOS Compact Scan 2M.We set a testing range of maximal values of normal vector threshold thr θ N = 1, . . . , 15. It means that the segmentation removes all points in clusters where the normal vector threshold is lower than thr θ N . Let n out be a number of the removed points and n in be several points in a point cloud.
We applied the region-growing algorithm on our data-sets with parameters thr θ N and we computed the number of removed points n out . The percentage of remaining points is defined as: We notice that the dependence between thr θ N and perc N has logarithmic progress, the logarithmic regression is the best choice to compute the dependence: The graphs and final results are in Section 3.1.

B-Spline Edge Representation
B-spline fitting belongs to the favorite tools for processing of a set of unorganized, possibly noisy data points in computer graphics, computer vision and CAD/CAM. The main advantage of B-spline approximation is based on the well-known formulation of B-spline [48]: where B i,p (u) is the B-spline basis function of degree p with control points P i . The parameter u is from the nonperiodic and nonuniform knot interval u =< u 0 , u 1 , . . . , u n+p+1 >. We set u 0 = 0 and u n+p+1 = 1. Let {Q k } m k=0 be the set of input points. The fitting algorithm search for control points P i , knot vector u of the spline curve defined in Equation (6) and also the parameters u = { u k } m k=0 that satisfied the equation: and First, we compute the initial curve shape. The authors in [44] work with the initial circle (center of the circle is at the mean of the point cloud, the radius of the maximal distance of a point to center) or an ellipse which the main axes are computed by PCA. We start with frequently used chordal and centripetal method. The chordal method computes the parameters as: where d = ∑ m k=1 |Q k − Q k−1 |. Centripetal method is similar, we only add square root of the line length. Centripetal method is more relevant in case that the data takes sharp turns.
Secondly, we set the knot vector u with condition that every knot span contains at least one u k . The internal knot spans are defined by [49]: and function f loor is the largest integer such i ≤ jd. Subsequently, coordinates of the control points P i are approximated by the standard technique of linear least-squares fitting [48]. We minimize the function: We can describe the solution in matrix form (the system of n − 1 × n − 1 equations) as: where matrix BB represents the B-spline function of degree p evaluated for all u k , k = 1, . . . , n − 1.
Matrix P are wanted control points and vector R is the vector of n − 1 points: and The Equation (11) can be solved by Gaussian elimination because the matrix B T B is positive definite and well-conditioned.
We must deal with the main problem throughout the iteration process: a choice of the distance measure to error determination. The formulation of the problem is simple, we have unorganized data points (in our case possible edge points) with non-uniform distribution with considerable noise. This problem can be formulated as a nonlinear optimization problem. Given input points {Q k } m k=0 , we want to compute control points P i that minimize a general objective function [44]: d 2 (P(u), Q k ) is the distance of point Q k to the curve P(t), f s is regularization term to ensure a smooth curve, λ is weight of f s . In our algorithm, we work with the asymmetric weighted point distance measure. We introduce a novel error term in Equation (15) that is based on golden ration. It ensures that the inner points (on the same side as a normal vector) have lesser influence than the outer. This advantage eliminates the problematic outliers and improves the shape of the curve ( Figure 6). The distance d k is the signed distance of the point Q k and the closest point on the computed B-spline curve. d k is positive if the point Q k is on the same side as normal n k -see Figure 7a. This error term is: k = 0, 1, 2 . . . , m, where σ defines the width of the transition of the weighting function with respect to the signed distance. The value ϕ is the golden ratio 1.61803398875. The objective function has a form: The weighting function w ϕ (d k ) and objective function f ϕ are modelled in Figure 7b. The asymmetric error detects the parts where the curve accuracy is insufficient. The improvement is done by local knot insertion. We know the parameter u k of the nearest curve point C(u k ) for every Q k . Therefore, we insert all knots u k in the insufficient parts and recompute the curve using knot insertion scheme. If u k ∈< u k , u k+1 ), then the new control points P i are: where In the parts, where the error is high, we add the new knot in the value of the nearest curve point. We tested the proposed improvement on the set of point clouds C 1 , C 2 , C 3 , C 4 and C 5 . The test results and thorough comparison are in Section 3.

Results and Discussion
This section contains two main parts: Estimation of normal threshold parameter and B-spline fitting with a novel asymmetric error. In the first part, we made a series of testing objects-from almost planar ones to the complex ones with holes. We proceeded these models with region-growing algorithm and using our proposed method (Section 2.1.2) we estimate the optimal value of the normal threshold. This value can be set as a universal input threshold value. The inexperienced end-users of the edge detection algorithm do not need any complex manual on how to control the results.
The second step contains a B-spline fitting algorithm for the edge representation. The input of this part are points detected in the first step. The asymmetric error with the golden ratio (Equation (15)) is used in the iteration process as described in the previous section. On the set of different curves (smooth, complex) we show the results and we make a comparison with method [44].
Feature detection algorithm was implemented in Microsoft Visual Studio C++ with Point Cloud Library 1.8.0. The B-spline approximation was programmed in MATLAB R2018a. The point cloud data-sets were obtained by ATOS Compact Scan 2M. Testing was provided on the MacBook Pro, 2.6 GHz Intel Core i7, 16 GB 1600 MHz DDR3, NVIDIA GeForce GT 750M 2048 MB, SSD.

Estimation of Normal Threshold thr θ N
The first analyses examines the impact of parameter thr θ N . The set of testing point clouds M i , i ∈ 1, 6 is in Figure 5. Let s in be several points in the given point cloud (Table 1). We applied the region-growing algorithm on our data-sets with the initial parameters thr θ N ∈ N, thr θ N = 1, . . . , 15 and we find the number of removed points n out . The percentage of remaining points is computed by: Table A1 provides the results obtained by the analysis of different values thr θ N . We can see that the value of perc N is tightly related to the complexity of the model. Table A2 (middle column) presents point charts of the values thr θ N and perc N . We can see that independently on the model complexity, it has the logarithmic progress. Therefore, the dependence between thr θ N and perc N can be generally described with logarithmic regression using Equation (5). The final equations and graphs are in Table A2 (last column).
We want to determine the best possible input value thr θ N that is applicable for common point clouds of the engineering models. We start with the average regression equation (the average of the equations in Table A2): The optimal value of perc is set empirically as 1-6%. Therefore, the values of thr θ N are in interval R(1), R(6) . Obviously, we get: The arithmetic mean of R(1) and R(6) is θ opt = 3.86 is set as a default initial value in the proposed feature detection algorithm. In the case of insufficient results, this value can be changed incrementally. Table 2 shows the segmentation results for models M i , i = 1, 2, . . . , 6 using the default initial value θ opt = 3.86. It is obvious that the results are sufficient and this value helps inexperienced users to start segmentation without any prior knowledge of curvature or normal threshold. By applying the incremental change of this value, the results are improved to obtain sufficient quality for actual use. Table 2. The final segmentation using the initial value θ opt = 3.86.

Input
Output Input Output We made a comparison of our method with similar approaches presented in [18,47]. We built a similar set of models for the evaluation (see Figure 8). The initial value of the curvature threshold was set to the pre-defined initial value θ N and this value was incrementally changed to obtain the best results. In Table 3, model in column (a), there is significant difference between the results-the incorrectly detected cylinder parts in results [18]. Our method with thr θ N = 4 provides visibly better results. Models (b), (c) and (d) are similar to models in [47]. We can see that our method provides comparable results. In the case of these models, we had to decrease the value of θ N because these models contain a lot of planar surfaces. Table 3. The comparison of the testing results (upper row: model from [18,47]), bottom row: edge detection with the incremental change of the threshold value).

B-Spline Curve Fitting
This section presents the results of the proposed B-spline approximation described in Section 2.2. We tested the computation on the set of curves C i , i = 1, . . . , 5 (Figure 9). The curves were extracted from point clouds-square curve (200 × 200 mm) from the cubic model, circle (radius 200 mm) from the cylindrical model. Remaining models were generated in Blender studio to prove the versatility of the algorithm. Table 4 gives the numbers of points of these curves.  The set of curve points was tested using presented B-spline curve fitting with adaptive knot insertion. The most difficult part was the optimal setting of the error computation-how to measure the quality of approximation if we do not know the correct shape. The common approach of the error computation works with point distance that offers the robustness. The tangent distance [50] provides the substantially faster convergence. The squared distance with asymmetric weighting [44,51] benefits from both and performs a stable method. We used the asymmetric weighting, but we proposed a novel evaluation of asymmetric error: We set the error threshold to 1 mm. In Figure 10a, we can see the initial B-spline curve (dashed blue line) that is computed by the chordal method and least square method described in Section 2.2. Red points are new B-spline control points and black crosses are the points in the point cloud. Figure 10b shows the distance of the point cloud points to the fitting curve (blue line segments). The value of asymmetric error determines whether we proceed the knot insertion and insert new control points in the inaccurate positions. Figure 11 shows the result of second iteration of B-spline fitting. We implemented the method [44] and proposed computation given by Equation (22). We computed average error of approximation E m and E gold in Table 5 (we set σ = 0.02, 0.09, 0.1). We can see that in the case of σ = 0.09 and σ = 0.1 the error is lower and the shape of the curve is visibly better (see Figure 11).  (c) Figure 11. The comparison of second iteration using method [44] (red line) and proposed method given by Equation (22) Table A3 shows different results with the number of control points and the corresponding asymmetric error. It is obvious that the algorithm is applicable in the case of different structures and work well with corner points.

Conclusions
The main contribution of this article is a method that improves the detection and representation of the edges of scanned engineering parts. This area is challenging because the scanning device (usually a 3D scanner) is frequently not able to cover the sharp edges with a sufficient amount of measured points. Therefore, the used techniques have a problem with the correct representation of the object shape. A manual time-consuming work in some commercial program is needed for fine-tuning of the scanned objects. Our article describes the overall workflow of the automatic edge detection in the arbitrary point cloud and subsequent correct edge representation using splines. We provide a user-friendly solution based on a single control parameter that can work with different shapes of engineering parts as well as other 3D scans generally. This method is also suitable e.g., for the reconstruction of the sharp edges of the buildings based on airborne laser scanning. The density of the measured points is also limiting the ability of the algorithms to reconstruct the precise shape of the walls. Simple adjustment of a single control parameter for the particular type of scanned object can substantially improve the results.
As the first step of our algorithm, we used the standard region-growing segmentation. We improved the part of principal component analysis (PCA) algorithm that is necessary for normal computation. We made a series of tests with variable shapes (from almost planar to complex ones) and using geometry calculus and numerical methods we implemented new variable-normal vector threshold thr θ N . Subsequently, we set this computed optimal value as the start value. Users can change incrementally this threshold in dependence on the shape of the object. Using one interconnected value to control the result is more suitable and very comfortable for end-users that are not able to understand to all algorithm parameters. Region-growing method removes the planar points and the remaining points are the possible edge points. The next step was challenging. How to compute the error of the edge representation when the correct edge points are missing. We chose to work with iterative B-spline method using the asymmetrically weighted point distance. The reasons are discussed in Section 3. We decided to use the golden ratio as a parameter (Section 2.2). It ensures that the inner points (on the same side as a normal vector) have lesser influence than the outer and the result is optically "perfect". This advantage eliminates the problematic outliers and improves the shape of the curve. We made the algorithm resistant to noise and prevented the zig-zag effect.
Our work has many possible extensions. Reverse engineering is tightly connected with 3D printing technology. Printed object correctness depends on the appropriate slices (G-code). The first follow-up can be an adjustment of the B-spline slicing method. Furthermore, the visualization of the scanned objects can be improved. Approximation of the curve points can be added to the triangulation net. Substantial advance can be an incorporation of the T-splines. T-spline enables irregular grid and we can represent sharp edges using different local knot vectors. Certainly, many scientific areas take advantage of huge progress in neural networks and deep learning. Point cloud analysis is not an exception. We follow the research in this area but the usage is now limited by the sensitivity of parameters and most methods work just with small point clouds [38].

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: