Hole Filling of Single Building Point Cloud Considering Local Similarity among Floors

: In the process of acquiring a point cloud, due to the shielding of the building itself and other ground objects (such as vegetation), the collected point cloud cannot cover the building facades surface uniformly and completely, and can produce several holes of different sizes. The existing hole-ﬁlling methods cannot obtain satisfying repaired results. To solve this problem, this paper proposes a new hole-ﬁlling approach considering the local similarity among ﬂoors. The proposed approach ﬁrst detects holes in the building facade surface and then identiﬁes the building ﬂoors; next, according to the fact that building data are similar in the same position on different ﬂoors, the holes of the building facade point cloud can be ﬁlled using the data of different ﬂoors. The study examines three large buildings to verify the proposed method. The experimental results show that the proposed method outperforms the state-of-the-art ﬁlling methods; the ﬁlling results of the proposed method are closer to the real shape and can be smoothly connected with the original building point cloud. Moreover, the proposed method has greater accuracy when comparing the root mean square error (RMSE).


Introduction
The oblique photogrammetry or airborne LiDAR is the main way to collect threedimensional (3D) data of large-scale urban scenes [1,2]. However, the 3D model of an urban scene constructed via oblique photogrammetry or airborne LiDAR technology comprises a whole and a set of many objects [3]. Although this kind of reconstructed 3D model can satisfy the scene visualization requirements, it is difficult to perform recognition and attribute assignment to target objects. Further, it cannot be applied to spatial analysis. Thus, to obtain a meaningful 3D model of each building, the single building extraction processing should be performed on the buildings, which are the main objects in the 3D model [2]. However, in acquiring the building point cloud, due to the shielding of the building itself and other ground objects (such as vegetation), the collected building point cloud is incomplete and may contain some holes. Especially at the bottom of the building, the point cloud of facade surfaces is incomplete after filtering processing due to the shielding of vegetation and vehicles. For example, Figure 1 shows an example of a building point cloud; there are many holes at the bottom of the building due to the shielding of vegetation. Therefore, these holes need to be filled to obtain the complete point cloud for each single building, which can lay a foundation for the automatic construction of its 3D model. data, some objects, such as windows, do not reflect the laser and appear as holes in LiDAR data. The second type group consists of the holes generated by the shieldin objects, such as vegetation and vehicles; these holes are mainly located at the bottom the building. The former type of hole is unique for the building and does not need to repaired, whereas the latter is a blind area generated by the obstruction of obstacle data acquisition. These holes do not belong to the structure of the building itself; t must be filled to achieve complete modeling of the building. With the commonly used hole-filling methods, data fitting and interpolation are formed on the hole areas to complete the hole filling according to the similarity betw each building hole and its neighborhood data. Such methods mainly include implicit face fitting [4,5] and mesh model repairing [6], which can process holes of various sha and sizes. However, they ignore the original data features of the building, and the paired data are relatively smooth and do not retain the sharp features of the building c pletely.
To solve the problems of existing methods, a new hole-filling approach conside the local similarity among floors for a building point cloud is proposed. The propo method can repair the holes of the building point cloud and lay a foundation for the c plete modeling of a single building. More specifically, our contributions can be sum rized as follows: (1) We propose a hole-filling method considering the local similarity among floors can repair holes automatically without manual intervention in the processing. (2) Experimental results on three different buildings show that the proposed met outperforms the state-of-the-art hole-filling methods.
The remaining part of this paper is organized as follows: Section 2 briefly introdu related works. The principle of the proposed hole-filling method is introduced in Sec 3. Section 4 presents the detailed experimental results and discussion. Finally, conclus and future work are presented in Section 5.

Related Works
The existing hole-filling methods for point clouds are mainly based on surface fit algorithms. Based on the high similarity between the hole and its neighborhood data, existing methods perform interpolation repair on the hole by analyzing the data's spa The holes of the building point cloud can be classified into two groups. The first group consists of the holes in the building itself. For instance, while collecting the LiDAR data, some objects, such as windows, do not reflect the laser and appear as holes in the LiDAR data. The second type group consists of the holes generated by the shielding of objects, such as vegetation and vehicles; these holes are mainly located at the bottom of the building. The former type of hole is unique for the building and does not need to be repaired, whereas the latter is a blind area generated by the obstruction of obstacles in data acquisition. These holes do not belong to the structure of the building itself; they must be filled to achieve complete modeling of the building.
With the commonly used hole-filling methods, data fitting and interpolation are performed on the hole areas to complete the hole filling according to the similarity between each building hole and its neighborhood data. Such methods mainly include implicit surface fitting [4,5] and mesh model repairing [6], which can process holes of various shapes and sizes. However, they ignore the original data features of the building, and the repaired data are relatively smooth and do not retain the sharp features of the building completely.
To solve the problems of existing methods, a new hole-filling approach considering the local similarity among floors for a building point cloud is proposed. The proposed method can repair the holes of the building point cloud and lay a foundation for the complete modeling of a single building. More specifically, our contributions can be summarized as follows: (1) We propose a hole-filling method considering the local similarity among floors that can repair holes automatically without manual intervention in the processing. (2) Experimental results on three different buildings show that the proposed method outperforms the state-of-the-art hole-filling methods.
The remaining part of this paper is organized as follows: Section 2 briefly introduces related works. The principle of the proposed hole-filling method is introduced in Section 3. Section 4 presents the detailed experimental results and discussion. Finally, conclusions and future work are presented in Section 5.

Related Works
The existing hole-filling methods for point clouds are mainly based on surface fitting algorithms. Based on the high similarity between the hole and its neighborhood data, the existing methods perform interpolation repair on the hole by analyzing the data's spatial correlations and using them to fit the surface model where the hole is located [7]. The commonly used surface fitting algorithms include the moving least squares (MLS) method [8], B-spline curve method [9][10][11], radial basis function interpolation (RBF) [5], and deep learning method [12]. Although the obvious features in the missing data area can be completely repaired by using the RBF to fit the hole surface, the constructed implicit surface is smooth, and it is difficult to ensure a smooth connection with the original data [13]. Wang proposed a hole-filling method in which the main idea is to first construct a mesh model for the point cloud, and then to repair the holes with an implicit surface [14]. However, this method was only suitable for small-sized holes. The feature point matching method uses the idea of "divide and conquer" to decompose the complex hole into several simple holes, and then uses the wave-front algorithm to complete the hole filling [15]. However, when the hole's area is large and essential information is lost, the local repair effect of the hole is poor since the important feature points cannot be matched. Wang used the gray theory to complete the repair and refinement of the target surface, which only had an ideal performance for the specific model and lacked broad applicability to other models [16].
Brunton proposed another method to fill the holes of a point cloud [17]; it projected the hole's boundary in the 3D models to the two-dimensional (2D) plane, used the improved Delaunay triangulation algorithm to divide the region of the hole on the 2D plane, and then used the minimized energy function to convert the 2D plane to the mesh surface in 3D space. Unfortunately, this method was limited by the size of the hole boundary and caused some damage to the original model to a certain extent. Based on Brunton's method, Fortes proposed a shape feature determination method [4]. Based on the minimum energy function, it generated hole patches for the original holes. However, this method requires good quality data around the holes. Wang proposed a method based on feature line extraction to fill the holes [18]. Using this method, feature lines were first extracted to obtain the sequence of control points of each hole, and Non-Uniform Rational B-Splines (NURBS) was used to fit and fill the holes. Frias proposed a method of filling the holes based on the superposition and comparison of the two cutwaters [19]. First, two cutwaters are segmented, then the translation-fitting of the whole cutwater over the damaged one is performed. Second, the damaged cutwater can be reconstructed when two cutwaters of the bridge point cloud are correctly aligned. The method can obtain better results, but it requires much manual intervention in the processing; therefore, the method is difficult to apply to buildings in the city.
Although the methods mentioned above can complete the repair of the holes in a building point cloud, they ignore the building's unique geometric and texture features in the filling process. The repair results for buildings with complex textures and geometric features are non-ideal, and the repaired data cannot be smoothly connected with the original data. These existing problems have motivated us to propose a new hole-filling method for a single building point cloud considering the inter-floor characteristics of the building. The proposed method can repair the holes and obtain the complete point cloud of the single building, which lays a foundation for constructing the 3D model of the single building.

Proposed Method
This section gives the main details of the proposed holes filling approach. First, the main idea of the proposed method is introduced. Then, the method of detecting holes in building facades is given. Next, the main processes of identifying building floors are presented. Finally, the hole-filling approach based on the local similarity among floors is proposed.

Main Idea
The geometric and textural features of high-rise buildings present some patterns, especially in China; features of the point cloud in similar parts between different floors of buildings have high similarity [20]. Our hole-filling method for a building point cloud is proposed. The main process is shown in Figure 2. First, the building point cloud is extracted from all point clouds, and the boundary of the hole is identified. Second, the building floors are identified according to the building's height. Third, the similarities Remote Sens. 2022, 14, 1900 4 of 16 between the hole data and the same area of adjacent floors are calculated. Finally, the holes are filled, and a 3D model of the single building is obtained.

Main Idea
The geometric and textural features of high-rise buildings present some patterns, especially in China; features of the point cloud in similar parts between different floors of buildings have high similarity [20]. Our hole-filling method for a building point cloud is proposed. The main process is shown in Figure 2. First, the building point cloud is extracted from all point clouds, and the boundary of the hole is identified. Second, the building floors are identified according to the building's height. Third, the similarities between the hole data and the same area of adjacent floors are calculated. Finally, the holes are filled, and a 3D model of the single building is obtained.

Building Rotation
Since the orientation of the building point cloud in the coordinate system is arbitrary, to better perform the following work and simplify the calculation, it is essential to rotate and translate the building parallel to the coordinate axis. The principal component analysis is a good approach to rotate buildings.
Given the building point cloud = { , , ⋯ , }, its centroid is ̅ and the covariance matrix is COV. Then, the eigenvector of the covariance matrix is calculated, which represents three mutually perpendicular coordinate axes. Therefore, a rectangular spatial coordinate is established, which is the new coordinate system of the building point cloud. Finally, the rotation of the building point cloud is achieved. The ̅ of P can be obtained from Equation (1). The covariance matrix COV of P can be obtained from Equation (2).

Building Rotation
Since the orientation of the building point cloud in the coordinate system is arbitrary, to better perform the following work and simplify the calculation, it is essential to rotate and translate the building parallel to the coordinate axis. The principal component analysis is a good approach to rotate buildings.
Given the building point cloud P = {p 1 , p 2 , · · · , p n }, its centroid is p and the covariance matrix is COV. Then, the eigenvector of the covariance matrix is calculated, which represents three mutually perpendicular coordinate axes. Therefore, a rectangular spatial coordinate is established, which is the new coordinate system of the building point cloud. Finally, the rotation of the building point cloud is achieved. The p of P can be obtained from Equation (1). The covariance matrix COV of P can be obtained from Equation (2).
The covariance matrix is computed; the three eigenvectors 3 , y 3 , z 3 ] respectively, and they correspond to the three mutually perpendicular directions of the building data. Then, the three eigenvectors are unitized and rotated parallel to the coordinate axis. When the point cloud is rotated, the z-axis of the building remains unchanged.
The rotation matrix R is computed using Equation (3), and all original data can be rotated to the new coordinate system. The relationship between the original point p where t x , t y , t z is the translation value of coordinates. To simplify the calculation, the coordinates of the point cloud of the building are shifted near the origin. The coordinate of the point at the left bottom for the building is [0, 0, 0] when the building is shifted.

Hole Boundary Extraction
(1) Searching k nearest neighbor points and constructing the perpendicular plane of a normal vector The k-d tree is a classical method of point cloud query which is a variant of a binary tree in multi-dimensional space [21]. Herein, the k-d tree was used to obtain the k nearest neighbor points (Q{q 1 , q 2 , · · · , q k }) of point p n . On how to choose the better k value, reference [22,23] gave some basic criteria. Based on reference [22,23], k is set as 20 in this paper. Afterward, points p n and Q are fitted into a spatial plane (x-y plane). Next, to compute the normal of p n x p , y p , z p as N(A, B, C), the perpendicular plane of the normal vector (L) can be obtained from Equation (5): Equation (5) is converted into the general form (Equation (6)), where D in Equation (6) can be converted into Equation (7): To judge the uniformity of the distribution of k nearest neighbor points (Q), it is important to project Q on the perpendicular plane of the normal vector (L), which is shown in Figure 3.
Remote Sens. 2022, 14, x FOR PEER REVIEW 5 of 16 where t , t , t is the translation value of coordinates. To simplify the calculation, the coordinates of the point cloud of the building are shifted near the origin. The coordinate of the point at the left bottom for the building is [0, 0, 0] when the building is shifted.

Hole Boundary Extraction
(1) Searching k nearest neighbor points and constructing the perpendicular plane of a normal vector The k-d tree is a classical method of point cloud query which is a variant of a binary tree in multi-dimensional space [21]. Herein, the k-d tree was used to obtain the k nearest neighbor points (Q{ , , ⋯ , }) of point pn. On how to choose the better k value, reference [22,23] gave some basic criteria. Based on reference [22,23], k is set as 20 in this paper. Afterward, points pn and Q are fitted into a spatial plane (x-y plane). Next, to compute the normal of ( , , ) as N( , , ), the perpendicular plane of the normal vector (L) can be obtained from Equation (5): Equation (5) is converted into the general form (Equation (6)), where D in Equation (6) can be converted into Equation (7): To judge the uniformity of the distribution of k nearest neighbor points (Q), it is important to project Q on the perpendicular plane of the normal vector (L), which is shown in Figure 3.  The coordinate of the neighboring point q i is (x q , y q , z q ), and the coordinate of its projection point is q i x q , y q , z q . Thus, the projection point q i can be computed from Equation (8): (2) Identifying hole border points The maximum angle difference is used to judge the uniformity of the distribution of k nearest neighbor points (Q) of point pn [24]. If the neighboring points are not distributed around the point pn, pn is viewed as the border point; otherwise, it is not. For instance, in Figure 4a, the neighboring points of the center point are not distributed around the center point; thus, it is a border point. In Figure 4b, the neighboring points of the center point are uniformly distributed around the center point; thus, it is not a border point. The coordinate of the neighboring point qi is ( , , ), and the coordinate of its projection point is ′ ( , , ). Thus, the projection point ′ can be computed from Equation (8): (2) Identifying hole border points The maximum angle difference is used to judge the uniformity of the distribution of k nearest neighbor points (Q) of point pn [24]. If the neighboring points are not distributed around the point pn, pn is viewed as the border point; otherwise, it is not. For instance, in Figure 4a, the neighboring points of the center point are not distributed around the center point; thus, it is a border point. In Figure 4b, the neighboring points of the center point are uniformly distributed around the center point; thus, it is not a border point.  In Figure 3, the edge ( ) is used as the starting edge, and the angles between each projection point and are computed. Then, these angles are sorted into ascending order: the maximum value (2π) and minimum (0) values are added to the angle set; thus, a new ordered set ( = {0, , , ⋯ , , 2 }) can be obtained. The differences in adjacent objects in are calculated from Equation (9), and the difference set = { , , ⋯ , } can be obtained. If there exists an > threshold , the point is determined to be a border point. The selection of threshold is based on reference [22,23]; in reference [23], they give a range of threshold that does not exceed 2 , and for a building point cloud, the threshold is 2 . Therefore, the threshold is set as 2 in this paper.
(3) Connecting border points The border points are connected using the vector deflection angle (Equation (10)) and boundary coherence [22,23]. For example, the adjacent points of point are and ( Figure 5). If only the distance is considered when the border points are connected, points , , and belong to the same boundary. However, they should belong to different boundaries. Thus, the vector deflection angle is used to sort the border points. The main processes of connecting border points are as follows: First, any border point can be viewed as the starting point, and then the boundary loop can be extracted using the breadth-first search. Next, any unconnected border point can be chosen as the initial point; the above process is performed. Finally, all boundaries can be identified exactly. In Figure 3, the edge (p n q 1 ) is used as the starting edge, and the angles between each projection point and p n q 1 are computed. Then, these angles are sorted into ascending order: the maximum value (2π) and minimum (0) values are added to the angle set; thus, a new ordered set (S = 0, α 1 , α 2 , · · · , α k−1 , 2π ) can be obtained. The differences in adjacent objects in S are calculated from Equation (9), and the difference set M = {M 1 , M 2 , · · · , M k } can be obtained. If there exists an M i > threshold δ, the point p n is determined to be a border point. The selection of threshold δ is based on reference [22,23]; in reference [23], they give a range of threshold δ that does not exceed π 2 , and for a building point cloud, the threshold δ is π 2 . Therefore, the threshold δ is set as π 2 in this paper.
(3) Connecting border points The border points are connected using the vector deflection angle (Equation (10)) and boundary coherence [22,23]. For example, the adjacent points of point N 3 are N 4 and N 5 ( Figure 5). If only the distance is considered when the border points are connected, points N 3 , N 4 , and N 5 belong to the same boundary. However, they should belong to different boundaries. Thus, the vector deflection angle is used to sort the border points. The main processes of connecting border points are as follows: First, any border point can be viewed as the starting point, and then the boundary loop can be extracted using the breadth-first search. Next, any unconnected border point can be chosen as the initial point; the above process is performed. Finally, all boundaries can be identified exactly.
where N i , N i+1 , and N i+2 are adjacent points.

Building Floor Identification
To fill holes by using similar features at the same place on different floors, it is essential to identify the different building floors. In identifying the building floors, the point cloud of the top appendage (such as the attic) and the basement are not considered. For instance, in Figure 1, there is a small attic that is located in the left front corner of the roof. The main processes of identifying the building floors are listed below: Traversing all points, the maximum and minimum elevation ( and ) of the building point cloud can be found.
The 3D point clouds are projected on the spatial grid (the grid size (D) is two to three times the point cloud density [25]); this process is mainly to avoid the effect of false elevation.
The maximum elevation ( ) of the point cloud in each grid unit is determined; the elevation value with the highest occurrence frequency is viewed as the evaluation of the roof elevation of the building (ℎ ).
To ensure that the lowest elevation of the point cloud in each grid unit is the facade root of the building, the point with an elevation that is more than in the grid units is removed.
The lowest elevation ( ) of the point cloud in each grid unit is determined, and the value with the highest occurrence frequency is viewed as the bottom elevation of the building (ℎ ).
The actual number of floors of the building is N, which is set based on visual judgment; thus, the height of each floor is h (Equation (11)). The same area of the building's adjacent floor is used to repair the holes of the building facade point cloud. According to Equation (12), the bounding box of the hole is shifted

Building Floor Identification
To fill holes by using similar features at the same place on different floors, it is essential to identify the different building floors. In identifying the building floors, the point cloud of the top appendage (such as the attic) and the basement are not considered. For instance, in Figure 1, there is a small attic that is located in the left front corner of the roof. The main processes of identifying the building floors are listed below: Traversing all points, the maximum and minimum elevation (H max and H min ) of the building point cloud can be found.
The 3D point clouds are projected on the spatial grid (the grid size (D) is two to three times the point cloud density [25]); this process is mainly to avoid the effect of false elevation.
The maximum elevation (H max ij ) of the point cloud in each grid unit is determined; the elevation value with the highest occurrence frequency is viewed as the evaluation of the roof elevation of the building (h max ).
To ensure that the lowest elevation of the point cloud in each grid unit is the facade root of the building, the point with an elevation that is more than H max +H min 2 in the grid units is removed.
The lowest elevation (H min ij ) of the point cloud in each grid unit is determined, and the value with the highest occurrence frequency is viewed as the bottom elevation of the building (h min ).
The actual number of floors of the building is N, which is set based on visual judgment; thus, the height of each floor is h (Equation (11)).

Filling Holes
Let the top height of each floor of the building with N floors be {h 0 , h 1 , h 2 , · · · , h N }. Given the minimum bounding box of a hole {(x max , y max , z max ), (x min , y min , z min )}, if h i−1 < z max ≤ h i and h j ≤ z min < h j+1 , the hole spans m (m = i − j) floors. The minimum bounding box of a similar area is {(x max , y max , z max + mh), (x min , y min , z min + mh)}. The coordinates of the data in a similar area and the hole are p source (x source , y source , z source ) and p target x target , y target , z target , respectively.
The same area of the building's adjacent floor is used to repair the holes of the building facade point cloud. According to Equation (12), the bounding box of the hole is shifted mh units upward, and the new bounding box of the hole can be obtained. If the new bounding box has complete data, the hole can be directly repaired. Otherwise, the bounding box of the hole continues to be shifted h units upward. The holes can be wonderfully filled by the above processes; moreover, the filled results can contain the texture information as well as geometry information.

Experimental Evaluation and Discussion
Herein, Unmanned Aerial Vehicle (UAV) oblique photography was conducted by "DJI Phantom 4" on the Nanjing Normal University campus; the datasheet of the "DJI Phantom 4" is shown in Table 1. Moreover, Pix4D [26] was used for data processing of the oblique photography images to obtain the whole area point cloud. The building boundary was identified from the DOM, and the building point cloud was extracted. An experiment was performed based on the building point cloud, and the effectiveness and accuracy of the proposed method were evaluated. In the experiments, k was set as 20 and the threshold δ was set as π 2 .  Figure 6 shows an extracted building point cloud, which is a pure point cloud that has filtered out vegetation and obstacles. Some holes that are marked by red circles appear on the facade surface of the building, which are mainly located on the bottom floor of the building near the ground.
Remote Sens. 2022, 14, x FOR PEER REVIEW 8 bounding box has complete data, the hole can be directly repaired. Otherwise, the bou ing box of the hole continues to be shifted ℎ units upward. The holes can be wonderf filled by the above processes; moreover, the filled results can contain the texture in mation as well as geometry information.

Experimental Evaluation and Discussion
Herein, Unmanned Aerial Vehicle (UAV) oblique photography was conducted "DJI Phantom 4" on the Nanjing Normal University campus; the datasheet of the Phantom 4" is shown in Table 1. Moreover, Pix4D [26] was used for data processing o oblique photography images to obtain the whole area point cloud. The building bound was identified from the DOM, and the building point cloud was extracted. An experim was performed based on the building point cloud, and the effectiveness and accurac the proposed method were evaluated. In the experiments, k was set as 20 and the thr old was set as 2 .  Figure 6 shows an extracted building point cloud, which is a pure point cloud has filtered out vegetation and obstacles. Some holes that are marked by red circles app on the facade surface of the building, which are mainly located on the bottom floor of building near the ground. The important part of filling holes is to identify hole boundaries correctly. Th holes in the building were identified using the maximum angle difference. Figure 7 sh the holes in the front facade of the building; the red marks within the yellow circle are identified hole boundaries. As depicted in Figure 7, the proposed method can identify accurate boundary of all holes with different shapes and sizes. Since the holes are in ferent positions, they have complex features. Thus, it is a challenging task to fill the ho The important part of filling holes is to identify hole boundaries correctly. The 13 holes in the building were identified using the maximum angle difference. Figure 7 shows the holes in the front facade of the building; the red marks within the yellow circle are the identified hole boundaries. As depicted in Figure 7, the proposed method can identify the accurate boundary of all holes with different shapes and sizes. Since the holes are in different positions, they have complex features. Thus, it is a challenging task to fill the holes.  Some results of filling the holes are presented in Figures 9 and 10, where the lo and right figures show the filled results of the front and lateral holes, respectively. proposed method can wonderfully fill the holes and preserve their original features. C paring the original data and filled results, it can be seen that the proposed method holes with different shapes and sizes quite well, and the complete building point c can be obtained. As expected, the filled results are the same as the original building.  Some results of filling the holes are presented in Figures 9 and 10, where the l and right figures show the filled results of the front and lateral holes, respectively proposed method can wonderfully fill the holes and preserve their original features. C paring the original data and filled results, it can be seen that the proposed method holes with different shapes and sizes quite well, and the complete building point c can be obtained. As expected, the filled results are the same as the original building. Some results of filling the holes are presented in Figures 9 and 10, where the lower and right figures show the filled results of the front and lateral holes, respectively. The proposed method can wonderfully fill the holes and preserve their original features. Comparing the original data and filled results, it can be seen that the proposed method fills holes with different shapes and sizes quite well, and the complete building point cloud can be obtained. As expected, the filled results are the same as the original building.  To better analyze the performance of the proposed method, it was compared with the MLS [8]. Figure 11 presents the filled results of some holes using the MLS and our method, respectively. Figure 11a,d,g shows the holes before filling; Figure 11b,e,h shows the repaired results using the proposed method. Figure 11c,f,i shows the filled results by the MLS. The point cloud of the repaired holes using the MLS is evenly distributed, but cannot be smoothly transitioned from the original point cloud, which is quite different from the original data. In contrast, the filled result by the proposed method can be smoothly connected with the original point cloud, which is similar to the distribution of the original building point cloud. Especially for the boundary, no discontinuity was  To better analyze the performance of the proposed method, it was compared the MLS [8]. Figure 11 presents the filled results of some holes using the MLS an method, respectively. Figure 11a,d,g shows the holes before filling; Figure 11b,e,h s the repaired results using the proposed method. Figure 11c,f,i shows the filled resu the MLS. The point cloud of the repaired holes using the MLS is evenly distributed cannot be smoothly transitioned from the original point cloud, which is quite diff from the original data. In contrast, the filled result by the proposed method ca smoothly connected with the original point cloud, which is similar to the distributi the original building point cloud. Especially for the boundary, no discontinuity To better analyze the performance of the proposed method, it was compared with the MLS [8]. Figure 11 presents the filled results of some holes using the MLS and our method, respectively. Figure 11a,d,g shows the holes before filling; Figure 11b,e,h shows the repaired results using the proposed method. Figure 11c,f,i shows the filled results by the MLS. The point cloud of the repaired holes using the MLS is evenly distributed, but cannot be smoothly transitioned from the original point cloud, which is quite different from the original data. In contrast, the filled result by the proposed method can be smoothly connected with the original point cloud, which is similar to the distribution of the original building point cloud. Especially for the boundary, no discontinuity was found, and the filled result and the original point cloud can be connected well. Thus, the proposed method is superior to the MLS; the filled holes are more appropriate. found, and the filled result and the original point cloud can be connected well. Thus, the proposed method is superior to the MLS; the filled holes are more appropriate. To better evaluate and analyze the performance of the proposed method, some areas of buildings were artificially selected to dig holes with different features and sizes. Then, the artificial holes were filled using different methods, which are marked by red rectangles in Figure 12. The filled results of artificial holes using different approaches are presented in Figure  13. Figure 13a,b shows the original point cloud of artificial holes 1 and 2, respectively; Figure 13c,d shows the corresponding repaired results using the proposed method; Figure  13e,f shows the corresponding filled results using the MLS. Overall, Figure 13 shows that To better evaluate and analyze the performance of the proposed method, some areas of buildings were artificially selected to dig holes with different features and sizes. Then, the artificial holes were filled using different methods, which are marked by red rectangles in Figure 12.
Remote Sens. 2022, 14, x FOR PEER REVIEW 11 of 16 found, and the filled result and the original point cloud can be connected well. Thus, the proposed method is superior to the MLS; the filled holes are more appropriate. To better evaluate and analyze the performance of the proposed method, some areas of buildings were artificially selected to dig holes with different features and sizes. Then, the artificial holes were filled using different methods, which are marked by red rectangles in Figure 12. The filled results of artificial holes using different approaches are presented in Figure  13. Figure 13a,b shows the original point cloud of artificial holes 1 and 2, respectively; Figure 13c,d shows the corresponding repaired results using the proposed method; Figure  13e,f shows the corresponding filled results using the MLS. Overall, Figure 13 shows that The filled results of artificial holes using different approaches are presented in Figure 13. Figure 13a,b shows the original point cloud of artificial holes 1 and 2, respectively; Figure 13c,d shows the corresponding repaired results using the proposed method; Figure 13e,f shows the corresponding filled results using the MLS. Overall, Figure 13 shows that the filled results using the proposed method are similar to the original building point cloud. Compared with Figure 13a,c, it is shown that the filled result can preserve the original building features and objects' structural information, such as windows, quite well. Compared with Figure 13a,e, it is shown that the filled results using the MLS fit the smooth surface without preserving the details. Thus, the filled results cannot be smoothly connected with the original point cloud. Furthermore, although the artificial holes span two floors, the filled results using the proposed method are ideal.  Figure 13a,c, it is shown that the filled result can preserve the original building features and objects' structural information, such as windows, quite well. Compared with Figure 13a,e, it is shown that the filled results using the MLS fit the smooth surface without preserving the details. Thus, the filled results cannot be smoothly connected with the original point cloud. Furthermore, although the artificial holes span two floors, the filled results using the proposed method are ideal. The filled results were further evaluated and analyzed using root mean square error (RMSE) [8]. Table 2 presents the RMSE of different filled results. The original data were The filled results were further evaluated and analyzed using root mean square error (RMSE) [8]. Table 2 presents the RMSE of different filled results. The original data were viewed as real values to calculate the RMSE of the two methods. The RMSE of the proposed method is less than that of the MLS. Therefore, the accuracy of the proposed method is higher, and the filled results using the proposed method are similar to the original data. Hence, the proposed method can effectively fill the complex holes of the building point cloud and preserve the original features of the building quite well, and the filled results are similar to the original building point cloud. From this comparison, we can see that the proposed method always outperforms the MLS. The 3D model of the single building can be constructed when the complete building point clouds are obtained; the result is shown in Figure 14. The 3D model has a complete surface that does not show any obvious holes. Therefore, the effectiveness and practicability of the proposed method can be verified. viewed as real values to calculate the RMSE of the two methods. The RMSE of the proposed method is less than that of the MLS. Therefore, the accuracy of the proposed method is higher, and the filled results using the proposed method are similar to the original data. Hence, the proposed method can effectively fill the complex holes of the building point cloud and preserve the original features of the building quite well, and the filled results are similar to the original building point cloud. From this comparison, we can see that the proposed method always outperforms the MLS. The 3D model of the single building can be constructed when the complete building point clouds are obtained; the result is shown in Figure 14. The 3D model has a complete surface that does not show any obvious holes. Therefore, the effectiveness and practicability of the proposed method can be verified. To better evaluate the effectiveness of the proposed method, the proposed hole-filling method is applied to two other buildings with different sizes of holes. Figures 15 and  16 show building A and B before filling, respectively. Due to the shielding of vegetation and other ground objects, the collected building point cloud is incomplete and contains various holes, especially at the bottom of the two buildings. As shown in Figure 15, the facade point cloud of building A contains large holes, which are difficult to fill. From Figure 16, we can observe that the first-floor point cloud of building B contains very unique holes.
The proposed method is applied to building A; Figure 17 presents the filled result. We can observe that the different holes can be repaired and the original features and objects' structural information can be preserved well. However, we can also find that the boundaries of large holes in two floors cannot be not filled well; some discontinuities are found. Figure 18 shows the repaired results of building B facade point cloud by the proposed method; it can be seen that the filled results preserve its original features. It can be seen that the proposed method can wonderfully fill the holes and preserve its original features. Moreover, the repaired result can be smoothly connected with the original point cloud, which is similar to the distribution of the original building point cloud. To better evaluate the effectiveness of the proposed method, the proposed hole-filling method is applied to two other buildings with different sizes of holes. Figures 15 and 16 show building A and B before filling, respectively. Due to the shielding of vegetation and other ground objects, the collected building point cloud is incomplete and contains various holes, especially at the bottom of the two buildings. As shown in Figure 15, the facade point cloud of building A contains large holes, which are difficult to fill. From Figure 16, we can observe that the first-floor point cloud of building B contains very unique holes.
The proposed method is applied to building A; Figure 17 presents the filled result. We can observe that the different holes can be repaired and the original features and objects' structural information can be preserved well. However, we can also find that the boundaries of large holes in two floors cannot be not filled well; some discontinuities are found. Figure 18 shows the repaired results of building B facade point cloud by the proposed method; it can be seen that the filled results preserve its original features. It can be seen that the proposed method can wonderfully fill the holes and preserve its original features. Moreover, the repaired result can be smoothly connected with the original point cloud, which is similar to the distribution of the original building point cloud.              Comparative experimental results were examined to evaluate the proposed method using three different datasets. The results show that the filling method proposed in this paper successfully repairs the holes caused by the shielding of the building itself and other ground objects. The filled result by the proposed method can be smoothly connected with the original point cloud, which is similar to the distribution of the original building point cloud. Through comparison, we can see that the proposed method has greater reliability and outperforms the state-of-the-art holes filling method MLS. The proposed method presents an important advance in filling holes; it can repair holes automatically without manual intervention in the processing. Therefore, it has value in the field of 3D geospatial information collection based on photogrammetry.   Although the proposed method shows a major improvement in filling holes, several aspects can be improved. First, the filled results for the large holes in more than two floors show some discontinuities, especially on the boundary. Second, the proposed method could be improved to be applied to atypical buildings, especially opera houses or stadiums.

Conclusions
In the view of the application targets of constructing 3D models of buildings, the holes of the building facade point cloud need to be filled. This paper proposes a novel method for filling holes of the building facade point cloud considering local similarity among floors. More importantly, the proposed method can fill holes automatically without manual intervention in the processing. The experimental results demonstrate that the proposed method can fill the holes caused by the shielding of the building itself and other ground objects (such as vegetation) quite well, and completely preserve the building's original geometric and textural features. Compared with the state-of-the-art method MLS, the proposed method can better preserve the original features of the building facades point cloud and has greater accuracy; moreover, the filled results can be smoothly connected with the original data. Thus, the filled result using the proposed method can better provide the complete building point cloud, and the 3D model of the building can be constructed based on these data. Although the proposed method represents a great advantage in filling holes, some limitations need to be addressed; for example, the proposed method is difficult to apply to floor buildings and atypical buildings.
In the next work, we intend to construct a refinement 3D model of the building, such as indoor. Second, for some buildings with atypical facades, mathematical morphology [27] could be a better method for filling large holes.
Author Contributions: X.M., Y.W. and Y.S. provided the core idea for this study. X.M. and Y.W. implemented the method and carried out the experimental evaluation. X.M., Y.W. and Y.S. wrote the main manuscript. Y.S., K.Z., J.Q. and Y.H. provided comments and suggestions for this paper. All authors have read and agreed to the published version of the manuscript.