A Novel Method for Density Analysis of Repaired Point Cloud with Holes Based on Image Data

Abstract

Repairing point cloud holes has become an important problem in the research of 3D laser point cloud data, which ensures the integrity and improves the precision of point cloud data. However, for the point cloud data with non-characteristic holes, the boundary data of point cloud holes cannot be used for repairing. Therefore, this paper introduces photogrammetry technology and analyzes the density of the image point cloud data with the highest precision. The 3D laser point cloud data are first formed into hole data with sharp features. The image data are calculated into six density image point cloud data. Next, the barycenterization Bursa model is used to fine-register the two types of data and to delete the overlapping regions. Then, the cross-section is used to evaluate the precision of the combined point cloud data to get the optimal density. A three-dimensional model is constructed for this data and the original point cloud data, respectively and the surface area method and the deviation method are used to compare them. The experimental results show that the ratio of the areas is less than 0.5%, and the maximum standard deviation is 0.0036 m and the minimum is 0.0015 m.

1. Introduction

3D laser scanners are widely adopted for the scanning of images. In the process of capturing the data while scanning, data corresponding to a partial point cloud may be missed due to the limitations of the instrument and also due to the environment. The point cloud with a hole cannot truly reflect the shape of the target; thus, it must be repaired.

Researchers have proposed many methods for repairing the point cloud with a hole. Ju [1] presented a robust method for repairing arbitrary polygon models. The method is guaranteed to produce a closed surface that partitions the space into disjoint internal and external volumes. Given any model represented as a polygon soup, an inside/outside volume using an octree grid is constructed, and the surface by contouring is reconstructed. Qiu et al. [2] established a triangle patch based on the point cloud around the hole. Subsequently, two types of isoparametric curves and their intersections are obtained by partitioned choosing of two arbitrary parameters within the equation of triangle patch, and the hole data are repaired by selecting the intersection points. This method relies on the ambient point cloud data and is suitable for curves with relatively small curvature. Bischoff et al. [3] proposed a repair method based on octree. In this method, the morphological operation is applied to establish the spatial topological relationship of the original data. The geometrical relationship and spatial morphology of the point cloud hole are reconstructed based on the topological relationship to achieve hole repair. Xin et al. [4] recognize the boundary of the hole based on the nature of the boundary edge of the adjacency triangle. An initial fill on the hole is accomplished using the wavefront method and the relationship of the included angle of the triangle vertex. The mesh of the hole is refined based on the curvature standard, and finally, an adjustment of geometry on the mesh vertex of the repaired hole is performed to make a natural transition with the surrounding grid. The experimental results indicate that this algorithm is simple, stable, and can repair different types of holes. Quinsat et al. [5] proposed a method to take the a priori knowledge of the numerical model as the nominal mesh. After identifying the digitized holes and calculating the differences between the nominal mesh and the point cloud, the nominal mesh is deformed. This deformation is determined by minimizing the deformation energy of the mesh. Centin et al. [6] proposed a method to repair the hole using Poisson surface reconstruction. An input mesh M along with its boundaries is given. An implicit function is derived by sampling a set of directional points from M, which are used to calculate Poisson surface reconstruction. An improved Delaunay refinement process is then introduced to generate hole patches with seamless transitions and no self-intersections. The hole patches are finally merged with the input mesh by robustly matching the input and by completing the boundary rings, which are stitched together to produce the final output. Li et al. [7] proposed a hole repair algorithm based on the Poisson equation. The predicted surface is fitted by solving the Poisson equation, which is triangulated and stitched seamlessly with the original hole. After that, the direction of the newly formed triangular surface is adjusted according to the normal vector information of the hole boundary region to achieve the effect of feature enhancement. Lin et al. [8] presented a novel feature-preserving hole-filling algorithm. The experimental data are divided into the featured holes and the nonfeatured holes. The spline guided tensor voting is proposed to restore the feature curves. The plane-guided tensor voting is proposed to restore the nonfeatured holes. Geng et al. [9] proposed a way to repair the holes in the terracotta warriors, provided that the missing parts of the model are stored in a database. The boundary of the hole is identified based on which fragment model with a roughly equal area is found in the database. Then, under the orthogonal constraint of double sparse representation, the optimal fragmentation model is predicted according to the vertex position error and edge smoothing error of triangular mesh. According to the registration function, the matching degree set of feature point pairs between the repaired model and the optimal fragmentary model is determined. The second-order umbrella operator is used to smooth the boundary of the model. Wang et al. [10] proposed a method based on the GA-BP neural network for the automatic repair of point cloud holes. The holes are identified and the interpolation points are selected by the method of equal step growth in the hole tone. Later, the interpolation points are taken as the input data of the GA-BP neural network model, and the predicted values are calculated to complete the repair of point cloud data. This method has high automation. Gai et al. [11] proposed a fitting approach to fill the holes based on structure from motion. After extracting the hole boundary by the fringe projection with a two-dimensional phase, the registration of the SFM point cloud and the fringe projection point cloud is carried out, and supplementary points are extracted. The holes are then filled based on a radial basis function on the point cloud added with the supplementary points. Zheng et al. [12] proposed a method of repairing a drill bit surface in the process of laser cladding a robot repairing drill bit. The 3D model of the worn drill bit is sliced by point cloud and the cubic B spline curve is used to fit the point cloud model. Finally, this method is used to select the actual machining points of the manipulator. Fan et al. [13] discussed a robust gap-filling method for extracting power lines from ground laser scanning data. A hierarchical clustering method is used to repair the gap based on the neighborhood relationship of the candidate nodes of the powerline.

Point cloud density is one of the main factors affecting precision. Guo et al. [14] pointed out that precision and robustness of feature extraction are greatly influenced by the point density. Additionally, an alternative way to improve the precision in low density point clouds is to increase the point density. Mat Zam et al. [15] obtained and registered the point cloud data of four different resolutions in the landslide area. The higher the resolution, the higher the density of the point cloud data. The experimental results show that the precision of ultra-high-resolution point cloud data is the highest, indicating that the higher the density of point cloud data, the higher the precision. Wang et al. [16] studied the optimum point cloud density for different scale DEM products under different terrain conditions. The experimental results show that when the terrain is undulating, it is necessary to increase the point cloud density appropriately to improve the precision. Zhang et al. [17] studied the effect of point cloud density on the precision of single leaf area index (LAI). With the same voxel size, the inversion value of the leaf area index of single wood increased with the increase of point cloud density. Du et al. [18] studied the effect of the density of the point cloud data on earthwork calculations. Su et al. [19] discussed the linear relationship between bamboo canopy volume and point cloud density. Zhou et al. [20] obtained the point cloud data of the mining subsidence area and found that higher point cloud resolution and stable subsidence prediction parameters resulted in a smaller fitting error and dense point cloud data, which also led to a higher modeling precision.

The state-of-the-art methods available to repair the hole in the point cloud utilize the spatial geometry relationship between the hole and its surrounding point cloud. However, for some point cloud data with sharp holes, there is no correlation between the missing data and the data around the holes, as shown in Figure 1b. The repairing result of using the curvature of the data around the holes is shown in Figure 1c. This repair method has low precision. Therefore, photogrammetry technology is introduced to repair the sharp hole based on image point cloud data. First, the point cloud is partly deleted manually, and these points are treated as missing point cloud data. Then, the image is resolved into image point cloud data to register with each other. Finally, due to the difference in the acquisition methods of these two kinds of data, the data precision is also different. The influence of density on the repair precision is investigated using image point cloud data for improving the precision. Optimal density is selected, whose precision results are compared to those obtained using original point cloud data. This method of repairing point cloud holes is not affected by the size of the study area, and it is convenient, fast, and has high precision.

The organization of this paper is as follows. Section “Materials and Methods” focuses on data preparation, from data collection to data processing. Six kinds of data are obtained, which are the combination of the missing point cloud data (M-point) and the image point cloud data (I-point) of six different densities. Section “Repair precision assessment on I-point” presents the selection of the optimal point cloud density (O-i-point) by the method of cross-section, and its precision is analyzed by the surface area and deviation method. Section “Discussion” presents the reason for affecting precision and the rationality of data results. Section “Conclusions” contains the concluding remarks and future work. The data processing flow chart is shown in Figure 2.

2. Materials and Methods

2.1. Design of Research Case

The surfaces of most of the buildings are composed of plane, smooth, and non-smooth curved structures. Therefore, these are the three main objects considered in this research work.

2.1.1. Plane Surface

The selected data for a plane surface were the intersection part between the horizontal plane and the vertical plane of the buildings. The research object was a statue shown in Figure 3. The outlined part in the red rectangular box was considered the research object. When the instrument was to be set at a low altitude, the sharp holes tended to occur within the outlined intersection part, leading to an incomplete point cloud.

2.1.2. Curved Surface

For the curved surface, smooth surface and non-smooth surface were chosen. The selected object was a landscape stone (LS), as shown in Figure 4. The outlined curved surfaces within the blue rectangular boxes were the smooth surface and the non-smooth surface, respectively, and were named S-LS and NS-LS.

2.2. The Laser Point Cloud Data of a 3D Object

2.2.1. Data Collection

A RIEGL-1000 3D laser scanner was used to capture the point cloud data. The scanning distance was about 5 m with a scanning interval of 0.002 m. After data preprocessing, such as de-noising, splicing, and compression of the point cloud data, the original point cloud data (O-point) were obtained. Here are the detailed preprocessing steps.

• De-noising on the point cloud

In this paper, the neighborhood average method was used for denoising. The k-D Tree algorithm [21] was used to establish the relationship between point cloud data. The neighborhood of any point in the point cloud data was then calculated. The average distance between this point and all point cloud data in k-neighborhood was finally calculated and compared with the threshold value. If the average distance was greater than the threshold value, it was considered a noise point and should be deleted; if the average distance was less than the threshold value, it was considered a non-noise point and should be retained. This method was implemented in MATLAB 2016 software. This method is suitable for point cloud data with relatively scattered noise points and has the characteristics of fast operational speed and obvious denoising effect, which is in line with the data characteristics of this paper.

2.

Splicing on the point cloud

Among the three kinds of data, only the statue was scanned in two stations, so point cloud data need to be spliced.

Point cloud splicing is mainly divided into two steps. First, coarse stitching: at least four pairs of points were manually selected from the two-point cloud data and completed initial registration according to the rigid transformation formula. Second, fine registration: the classical ICP algorithm [22] was used to complete the precise registration of the point cloud. This method is implemented in Riscan Pro software. The final error was shown to be 0.0026 m.

3.

Compression on the point cloud

In this paper, the minimum distance method [23] was used to compress point cloud data. The minimum distance method is simple, easy to implement, and fast in operation. This method is implemented in MATLAB 2016.

2.2.2. Calculation of Density

The density of a 3D point cloud data is defined as the number of points within a unit area. The greater the number of points, the higher the density of the point cloud. The density can also be interpreted as the average distance between a point and other adjacent points. The shorter this distance, the denser the point cloud. In this paper, these two methods were used to calculate the density of the point cloud.

Calculation of Density Based on Distance

Based on k-nearest algorithm, m adjacent points $q_{ij}\left(j=1,2,3\dots m\right)$ around a point $p_i\left(i=1,2,3\dots n\right)$ were selected successively.

The average distance between every point and adjacent points was calculated according to Equation (1) [24]. Then, according to Equation (2), the density of point cloud data was calculated by averaging the density of all of the points.

$$\overline{d_i}=\frac{1}{m}{\displaystyle \sum }_{j=1}^{m}D\left(p_i,q_{ij}\right)$$

(1)

In Equation (1), p_i is the point cloud; q_ij is the adjacent point of p_i; m is the number of adjacent points; ${\overline{d}}_i$ is the average distance of each point from its adjacent points, which is the density at this point; and $D\left(p_i,q_{ij}\right)$ denotes the distance between point p_i and its adjacent points q_ij.

$$\overline{d}=\frac{1}{n}{\displaystyle \sum }_{i=1}^n\overline{d_i}$$

(2)

In Equation (2), $\overline{d}$ is the average density of all the points.

Calculation of Density Based on the Mesh Grid

Uniform mesh grids are constructed on the point cloud, and the density is interpreted as the average of the number of points within one mesh grid [24]. The interval of mesh grids is denoted as $\overline{d}$ or multiples of $\overline{d}$, and the mesh grid is constructed. The number of solid mesh grids and the number of point clouds within every solid mesh grid were recorded. The average value was calculated and rounded to obtain the density of the point cloud.

Based on the above methods, the densities of the O-point were calculated and are listed in

Table 1. The experimental data’s point cloud density.

Classification	Density-Based on Distance (m)	Density-Based on Mesh Grid (Pieces)
statue	0.0029	48
S-LS	0.0019	219
NS-LS	0.0017	165

.

2.2.3. Generation of Hole

As per the experimental requirement, some points were manually deleted from the original data (O-point) to obtain the missing point cloud data (M-point), as shown in

Table 2. Schematic of O-point and M-point.

Classification	Statue	S-LS	NS-LS
O-point
M-point

.

2.3. Image Data

The camera used in this paper is a non-measurement type Nikon D600 model with 2426 million active pixels and the highest resolution of 6016 × 4016. The image data are captured based on the requirement of photos.

The structure-from-motion (SFM) algorithm and the patch-based multiple view system (PMVS) algorithm based on the path model were utilized to resolve the image data into image point cloud data (I-point), as shown in Figure 5. Here are the detailed steps.

• Feature extraction and matching [25]: SIFT algorithm was used to extract the image features, and image matching was carried out. The basic matrix was calculated by the RANSAC algorithm, and the wrong matching pairs were eliminated;
• Zhang Zhengyou camera calibration [26]: The internal parameter matrix of the camera was calibrated by the Zhang Zhengyou camera calibration method;
• Camera attitude calculation;
• Point cloud data generation by triangulation method [27];
• Bundle adjustment: A bundle adjustment algorithm was introduced to optimize the camera attitude and 3D coordinates;
• Dense matching based on patch model [24]: PMVS was used for image matching, diffusion, and filtering to generate dense point clouds.

The whole process was implemented in PhotoScan Pro software.

Because the precision of the M-point was higher than that of the I-point, resolved from the image, to ensure the precision of repair, the density of the I-point was increased by a factor of 1, 1.2, 1.5, 2, 2.5, and 3 compared to the density of M-point.

2.4. Registration of Two Kinds of Data

The coarse registration was accomplished by manual selection of points. To improve the repair precision, fine registration must be carried out. The Bruse model was one of the seven parameter coordinate transformations, which is suitable for small angles, and so coarse registration must be carried out before using it [28]. When it was applied in a small range, there was a strong correlation between translation parameters and rotation and scale parameters, which led to an ill-conditioned solution model [29,30]. Moreover, the rotation center was located at the origin of coordinates, so the point cloud data may not have been evenly distributed in all quadrants, which would reduce the precision. Based on these two reasons, the origin of the coordinate system was redefined as the barycentric position, and the translation parameter was eliminated to improve precision, which is the barycenterization Bursa model [31].

Therefore, the barycenterization Bursa model was adopted to achieve fine registration between M-point and the above six different density data.

Finally, the adjacent points of the I-point within the M-point can be found using the k-nearest algorithm. I-point-1, I-point-1.2, I-point-1.5, I-point-2, I-point-2.5, and I-point-3 can be obtained by deleting the point cloud data within the superposition section, as shown in

Table 3. The registration results of two kinds of data.

Classification	I-Point-1	I-Point-1.2	I-Point-1.5	I-Point-2	I-Point-2.5	I-Point-3
statue
S-LS
NS-LS

.

3. Repair Precision Assessment on I-Point

3.1. Optimal Density Based on the Method of Cross-Section

The method of cross-section selects two arbitrary planes to truncate the original 3D point cloud data, and the intersection lines and points that are formed are projected onto a plane for comparison. In most cases, intersection points are more commonly used than intersection lines [32]. In this paper, two planes were used to truncate the point cloud data, and the point cloud data between the two planes were projected onto one plane for comparison. The selected position in which the holes are to be repaired is shown in

Table 4. The truncated point cloud data by the method of cross-section.

Classification	Statue	S-LS	NS-LS
The position truncated plane using the method of cross-section
The truncated display of point cloud data

.

It is difficult to directly observe the difference between I-point and O-point using different densities in the selected cross-section. Therefore, in this paper, curve fitting was performed on the selected cross-section to quantitatively analyze the effect of density to make an optimal density choice. For convenience, the cross-section data of O-point is called section-origination (S-origination), and the cross-section data of I-point is called section-1 (S-1), section-1.2 (S-1.2), section-1.5 (S-1.5), section-2 (S-2), section-2.5 (S-2.5), and section-3 (S-3).

3.1.1. Quantitative Analysis of Point Cloud Cross-Section of the Statue

As shown in Table 4, the selected cross-section of the statue can be approximately viewed as a combination of two straight lines, i.e., the point cloud data is fitted as a piecewise function, indicated in Equation (3).

$$f\left(x\right)=\{\begin{array}{c}\begin{array}{cc}{k}_{1}x+b_1& x<a\end{array}\\ \begin{array}{cc}{k}_{2}x+b_2& x\ge a\end{array}\end{array}$$

(3)

In Equation (3), k₁, k₂, b₁, b₂ are fitted coefficients.

The parameters of the piecewise function obtained by fitting the point cloud data are shown in

Table 5. The parameters of fitted piecewise function for the statue.

	Fitted Coefficient k1	Fitted Coefficient b1	Fitted Coefficient k2	Fitted Coefficient b2
S-origination	50.0512	15.9863	0.0187	−4.6980
S-1	49.9438	15.9408	0.0308	−4.6936
S-1.2	47.1638	14.7943	0.0352	−4.6919
S-1.5	48.9306	15.5231	0.0367	−4.6913
S-2	45.5217	14.1155	0.0340	−4.6924
S-2.5	46.6630	14.5862	0.0399	−4.6905
S-3	46.8562	14.6660	0.0370	−4.6915

. Figure 6 shows the linear fitted plots of S-origination and S-2.

In Figure 5, we select points in the x-axis and calculate the corresponding y coordinates of the fitted lines. Then, a difference between the calculated y values of various densities and the y values in the original data are obtained, and an average is calculated. For the vertical line, a single point is selected on the x-axis, and for the horizontal line, multiple points taken on the x-axis at regular intervals are selected. The y values obtained from curve fitting are shown in

Table 6. Comparison of data points of statue.

	x (m)	−0.4142	−0.4130	−0.3980	−0.3830	−0.3680	−0.3530	−0.3380	Difference (mm)
y (m)
S-origination		−4.6878	−4.7057	−4.7054	−4.7051	−4.7048	−4.7046	−4.7043
S-1		−4.6889	−4.7063	−4.7058	−4.7054	−4.7049	−4.7045	−4.7040	−0.3
S-1.2		−4.6871	−4.7064	−4.7059	−4.7053	−4.7048	−4.7043	−4.7038	0.0
S-1.5		−4.6881	−4.7064	−4.7059	−4.7053	−4.7048	−4.7042	−4.7037	−0.1
S-2		−4.6877	−4.7064	−4.7059	−4.7054	−4.7049	−4.7044	−4.7039	−0.1
S-2.5		−4.6884	−4.7069	−4.7063	−4.7057	−4.7051	−4.7045	−4.7039	−0.4
S-3		−4.6883	−4.7068	−4.7063	−4.7057	−4.7051	−4.7046	−4.7040	−0.4

. The results shown in Table 6 indicate that the average difference in data with various densities is small. With the increase in density, the difference values first decrease and then increase. The minimum difference is achieved for S-1.2.

3.1.2. Quantitative Analysis of Point Cloud Cross-Section of S-LS

Table 4 indicates that the selected cross-section of S-LS can be approximately regarded as a combination of several curves. By comparing the fitting precision of degrees 1–9 of the polynomials, the appropriate degree is selected. The sum of squares due to error (SSE), root mean squared error (RMSE), and coefficient of determination (R-square) are selected to evaluate the precision.

The fitting results are shown in Figure 7, and the precision analysis is shown in

Table 7. Precision analysis of the different curve fitting polynomials used for the S-LS cross-section.

Degree of the Polynomial	SSE	R-Square	RMSE
1	0.6232	0.5691	0.0088
2	0.5888	0.5928	0.0085
3	0.4413	0.6948	0.0074
4	0.0903	0.9375	0.0033
5	0.0856	0.9408	0.0033
6	0.0199	0.9862	0.0016
7	0.0151	0.9896	0.0014
8	0.0095	0.9934	0.0011
9	0.0075	0.9948	0.0010

.

Based on the principle of simple calculation and high precision, the polynomial fitting of degree six is utilized as given by Equation (4).

$$f\left(x\right)=p_1\times x^6+p_2\times x^5+p_3\times x^4+p_4\times x^3+p_5\times x^2+p_6\times x^1+p_7$$

(4)

In Equation (4), p₁~p₇ are fitted coefficients.

The parameters of the function obtained by fitting the point cloud data are shown in

Table 8. The parameters of curve fitted function of S-LS.

	Fitted Coefficient p1	Fitted Coefficient p2	Fitted Coefficient p3	Fitted Coefficient p4	Fitted Coefficient p5	Fitted Coefficient p6	Fitted Coefficient p7
S-origination	17.3606	96.8805	218.3503	253.2991	158.7827	50.8180	−1.8516
S-1	18.7270	103.4032	230.5510	264.5488	164.0353	51.9402	−1.7779
S-1.2	18.7302	103.4309	230.6409	264.6915	164.1535	51.9887	−1.7702
S-1.5	18.9817	104.8192	233.7631	268.3472	166.5003	52.7711	−1.6645
S-2	16.8416	93.1438	207.9368	238.7717	148.0528	46.8403	−2.4315
S-2.5	14.9288	82.6928	184.4851	211.1082	129.9656	40.6306	−3.3048
S-3	14.9171	82.6350	184.3714	210.9968	129.9096	40.6175	−3.3057

.

Figure 8 Shows the fitted curve of S-origination and S-2.

In Figure 4, points on the x-axis are selected at equal intervals and the y values of fitted curves are calculated as shown in

Table 9. Comparison of data points of the S-LS.

	x (m)	−1.3	−1.15	−1	−0.85	−0.7	−0.55	Difference (mm)
y (m)
S-origination		−8.3545	−8.3477	−8.3557	−8.3424	−8.3166	−8.3273
S-1		−8.3544	−8.3469	−8.3567	−8.3462	−8.3195	−8.3275	−1.0
S-1.2		−8.3544	−8.3469	−8.3568	−8.3462	−8.3195	−8.3276	−1.0
S-1.5		−8.3545	−8.3469	−8.3568	−8.3462	−8.3194	−8.3276	−1.0
S-2		−8.3530	−8.3470	−8.3561	−8.3456	−8.3202	−8.3274	−0.7
S-2.5		−8.3540	−8.3478	−8.3570	−8.3459	−8.3201	−8.3271	−1.1
S-3		−8.3540	−8.3478	−8.3569	−8.3459	−8.3201	−8.3271	−1.1

. The results indicate that the average difference in data with various densities is small. It can be observed that with the increase in density, the average difference decreases first, and then it increases. The minimum difference is achieved for S-2.

3.1.3. Quantitative Analysis of Point Cloud Cross-Section of NS-LS

Table 4 indicates the existence of sharp holes in the selected cross-section of NS-LS, which can also be approximately regarded as a combination of several curves. As in Section 3.1.2, the curve fitting results using polynomials of different degrees for the cross-section of NS-LS are shown in Figure 9, and the precision analysis is shown in

Table 10. Precision analysis of the different curve fitting polynomials used for the NS-LS cross-section.

Degree of the Polynomial	SSE	R-Square	RMSE
1	16.5800	0.0848	0.0255
2	6.1830	0.6587	0.0156
3	6.1340	0.6614	0.0155
4	3.0960	0.8291	0.0110
5	1.8100	0.9001	0.0084
6	1.7700	0.9023	0.0083
7	0.9385	0.9482	0.0061
8	0.6440	0.9645	0.0050
9	0.6379	0.9648	0.0050

.

Based on the principle of simple calculation and high precision, the polynomial fitting with nine degrees is utilized as given by Equation (5).

$$f\left(x\right)=p_1\times x^9+p_2\times x^8+p_3\times x^7+p_4\times x^6+p_5\times x^5+p_6\times x^4+p_7\times x^3+p_8\times x^2+p_9\times x^1+p_{10}$$

(5)

In Equation (5), p₁~p₁₀ are fitted coefficients.

The parameters of the function obtained by fitting the point cloud data are shown in

Table 11. The parameters of curve fitted function of NS-LS.

	Fitted Coefficient p1	Fitted Coefficient p2	Fitted Coefficient p3	Fitted Coefficient p4	Fitted Coefficient p5
	Fitted Coefficient p6	Fitted Coefficient p7	Fitted Coefficient p8	Fitted Coefficient p9	Fitted Coefficient p10
S-origination	−8.5265 × 102	−2.8326 × 104	−4.1685 × 105	−3.5655 × 106	−1.9528 × 107
	−7.0996 × 107	−1.7124 × 108	−2.6411 × 108	−2.3617 × 108	−9.3216 × 107
S-1	−8.3842 × 102	−2.7887 × 104	−4.1092 × 105	−3.5198 × 106	−1.9308 × 107
	−7.0318 × 107	−1.6994 × 108	−2.6268 × 108	−2.3549 × 108	−9.3219 × 107
S-1.2	−8.4677 × 102	−2.8183 × 104	−4.1558 × 105	−3.5626 × 106	−1.9560 × 107
	−7.1308 × 107	−1.7253 × 108	−2.6704 × 108	−2.3976 × 108	−9.5073 × 107
S-1.5	−8.4549 × 102	−2.8146 × 104	−4.1514 × 105	−3.5597 × 106	−1.9550 × 107
	−7.1295 × 107	−1.7256 × 108	−2.6719 × 108	−2.4000 × 108	−9.5222 × 107
S-2	−7.9051 × 102	−2.6244 × 104	−3.8591 × 105	−3.2978 × 106	−1.8043 × 107
	−6.5516 × 107	−1.5780 × 108	−2.4296 × 108	−2.1683 × 108	−8.5375 × 107
S-2.5	−7.5988 × 102	−2.5183 × 104	−3.6959 × 105	−3.1515 × 106	−1.7200 × 107
	−6.2282 × 107	−1.4954 × 108	−2.2939 × 108	−2.0384 × 108	−7.9854 × 107
S-3	−7.5027 × 102	−2.4852 × 104	−3.6452 × 105	−3.1063 × 106	−1.6941 × 107
	−6.1294 × 107	−1.4702 × 108	−2.2528 × 108	−1.9993 × 108	−7.8199 × 107

.

Figure 10 shows the fitted curve of S-origination and S-2.

In Figure 10, points on the x-axis are selected at equal intervals and the y values of the fitted curves are calculated as shown in

Table 12. Comparison of data points of NS-LS.

	x (m)	−4.25	−4.1	−3.95	−3.8	−3.65	Difference (mm)
y (m)
S-origination		8.1422	8.1022	8.0749	8.0394	8.0799
S-1		8.1433	8.1023	8.0748	8.0395	8.0810	0.4
S-1.2		8.1435	8.1025	8.0749	8.0396	8.0810	0.5
S-1.5		8.1432	8.1024	8.0748	8.0396	8.0810	0.4
S-2		8.1428	8.1023	8.0747	8.0396	8.0809	0.3
S-2.5		8.1430	8.1025	8.0748	8.0398	8.0807	0.4
S-3		8.1427	8.1025	8.0748	8.0399	8.0808	0.4

. The results shown in Table 12 indicate that the average difference of data with various densities is small, and the minimum difference is achieved for S-2.

Based on the above analysis, the maximum and minimum average values of every density are small. The optimal density for the three kinds of data is 1.2 times of M-point for the statue, 2 times of M-point for S-LS and NS-LS Importantly, the large amount of point cloud data leads to a slow processing speed, and the large difference in the density of point cloud data also affects the 3D modeling process. The average difference between the selected optimal density and S-1.2 data is relatively small with a small amount of data. Taking all of the above factors into account, the I-point with 1.2 times of M-point is taken as the optimal image point cloud data (O-i-point). The O-i-point and M-point are combined to act as the combining point cloud data (C-point), as illustrated in the third column “I-point-1.2” of Table 3.

3.2. Precision Assessment

The reliability of the proposed method in repairing the point cloud hole using image data has to be discussed, mainly for the C-point data. The commonly used assessment methods include surface area method, volume method, deviation evaluation, cross-section method, etc. For this paper, the experimental data are available for a plane but not for volume, and so the volume method is not applicable for precision assessment. The surface area method and deviation method are utilized to perform precision assessment because the cross-section method has been used to select O-i-point. For the convenience of expression, the constructed model based on O-point is referred to as the original model (O-model), and the constructed model based on C-point is referred to as the repair model (R-model), as shown in

Table 13. Comparison of O-model and R-model.

Classification	Statue	S-LS	NS-LS
O-model
R-model

.

3.2.1. Surface Area Assessment

The areas of both the O-model and R-model are calculated and compared to make a judgment on whether the proposed method reduces the detailed surface characteristics of original complete point cloud data [33]. The area of a triangle can be calculated by Heron’s formula [32], and is given by Equation (6).

$$S_{\Delta ABC}=\sqrt{L\cdot \left(L-a\right)\cdot \left(L-b\right)\cdot \left(L-c\right)}$$

(6)

In Equation (6), a, b, c are three side lengths of the triangle, and L denotes half of the perimeter of the triangle.$S_{\Delta ABC}$ is the area of the triangle. A, B, C are the three vertices of a triangle.

In software Geomagic Studio 2013, the surface areas calculated using the two models are shown in

Table 14. The surface area of the model of experimental data.

	Model Classification	Area (m2)	Difference of Area (m2)	Area Ratio
statue	O-model	1.0721	0.0082	0.476%
	R-model	1.0803
S-LS	O-model	2.3564	0.0085	0.361%
	R-model	2.3649
NS-LS	O-model	2.2167	0.0094	0.424%
	R-model	2.2261

.

3.2.2. Deviation Assessment

The deviation method analyzes the spatial difference between the tested model and reference model by calculating their shortest distance, that is the distance between the nearest triangular of the two models [34]. In this paper, the test model is the R-model and the reference model is the O-model. This method is implemented in the software Geomagic Studio 2013. The results of the deviation analysis for all of the experimental data are shown in Figure 11. The deviation analysis values are shown in

Table 15. Deviation analysis values of each experimental data.

	Maximum Value (mm)		Minimum Value(mm)		Average Value(mm)		Standard Deviation Value (mm)
	Positive Deviation	Negative Deviation	Positive Deviation	Negative Deviation	Positive Deviation	Negative Deviation
statue	3.7	−3.4	0.1	−0.0	0.6	−1.4	1.5
S-LS	10.0	−10.0	0.2	−0.1	2.2	−1.3	3.1
NS-LS	11.0	−11.0	0.4	−0.2	3.2	−1.1	3.6

.

In Figure 11, the green color denotes zero deviation, the blue color denotes negative deviation, and the red color denotes a positive deviation. The darker the color, the higher the deviation.

In Table 15, the standard deviation for the statue is 1.5 mm. In Figure 11a, the vertical plane is mostly light blue in color, and the horizontal plane is light yellow, almost green, which indicates a relatively small deviation. The horizontal plane is less deviated than the vertical plane. It may be because of errors in the acquisition of the image data.

In Table 15, the standard deviation for S-LS is 3.1 mm. In Figure 11b, the deviation is relatively small because the center of the hole is light in color. The farther to the boundary, the darker the color, indicating the larger the deviation. Additionally, the left side is positive deviation, the right side is negative deviation. This phenomenon may be due to the error of point cloud registration. Although the point cloud should be as uniform as possible in manual selection points, human error cannot be avoided.

In Table 15, the standard deviation for NS-LS is 3.6 mm. In Figure 11c, the deviation is relatively small because the boundary of the hole is a light color, and the deviation is large because the non-smooth part at the center of the hole is a dark color. Additionally, the boundary side is positive deviation, the center side is negative deviation. It shows that the precision of image data is less than that of 3D laser scanning for sharp feature points.

On the whole, the statue has the smallest standard deviation, followed by S-LS and NS-LS. It shows that the precision of image data is highest for the plane and lowest for the non-smooth surface.

4. Discussion

In the previous method of repairing sharp holes in 3D point cloud data with image data, it was not considered whether the density of image data would affect the repair precision. Through analysis in this paper, it was found that the precision of 1.2 times of the corresponding M-point was higher. Here are the factors that affect the precision.

• Precision of point cloud data acquisition

The precision of the point cloud will change with the distance [35]. The farther the distance is, the smaller the precision is. The Riegl VZ-1000 3D laser scanning imaging system was adopted in this paper. The scanning precision of the system within 100 m is 5 mm, the scanning distance in this paper is about 5 m, and the precision of point cloud reaches 5 mm, which meets the requirements of general applications.

2.

Precision of image data acquisition

The non-measurement camera used to obtain image data should have a fairly high resolution. The Photoscan software requires a camera with a minimum of 50,000 pixels. Nikon D600, with an effective pixel of 24.26 million, was selected in this paper.

3.

Resolution precision of image data

The solution of image data was to calculate two-dimensional image data into three-dimensional point cloud data, and the calculated point cloud data were used for hole repair. Therefore, the solution precision of image data is directly related to the repair precision.

4.

Splicing precision of two kinds of data

The splicing precision of point cloud data and image point cloud data directly affects repair precision [36]. The splicing precision is generally measured by the average distance between the nearest points, and the splicing precision in this paper reached the sub-millimeter level. Compared with the literature, the precision is higher [37].

5.

Precision relationship between point cloud data and image point cloud data

According to existing research, the higher the density of the point cloud, the higher the precision [14,15,16,17,18,19,20]. Since the precision of laser point cloud data in this paper was slightly higher than that of image point cloud data, it is reasonable to conclude that the density of the optimal image point cloud data is slightly higher than that of the point cloud.

If the position of the target object and the surrounding environment change, the precision of the two data sources will change. In general, the image data are observed at close range. It can be considered that the precision of image data is stable, and the precision of point cloud data will decrease with the increase of the distance between the station location and the target object. Therefore, with the increase of distance, the precision of image data remains unchanged, while the precision of point cloud data gradually decreases. In this way, the precision of the two kinds of point cloud data will gradually approach. Therefore, the optimal density of the image point cloud data calculated in this paper may change.

6.

Fusion of hole edge data

The data at both ends of the hole edge position are from two data sources, respectively, and the precision and density may be different. Therefore, the fusion of edge data will appear abrupt, which may cause the result of no smooth transition of edge data, thus, affecting the repair precision [38].

5. Conclusions

In this paper, holes in laser scanner point cloud data were repaired based on photogrammetry technology. The precision of the 3D point cloud was slightly larger than that of the image data. To improve the precision of repair, several I-points with different densities were generated and compared to select the optimal density using the cross-section method. Experimental results indicated that the density of the selected optimal I-point is 1.2 times the corresponding M-point by taking the number of point clouds into account. The O-i-point and M-point were combined for comparison with the O-point using the surface area method and the deviation method. The experimental results indicated that the area ratio was less than 0.5% when the surface area method was used, and the maximum and minimum standard deviations were, respectively, 0.0036 m and 0.0015 m when the deviation method was used.

When evaluating the precision, this paper did not evaluate the boundary part after combining two types of data. It remains to be studied whether the edge part is smooth or not. Additionally, the scope of experimental data in this paper is small, and whether the large range of data is suitable for the research results in this paper remains to be studied.