# An Improved Building Reconstruction Algorithm Based on Manhattan World Assumption and Line-Restricted Hypothetical Plane Fitting

Authors: Xiaoguo Zhang; Guo Wang; Ye Gao; Huiqing Wang; Qing Wang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

An improved patch-based multiview stereo (PMVS) algorithm based on Manhattan world assumption and the line-restricted hypothetical plane fitting method according to buildings’ spatial characteristics is proposed. Different from the original PMVS algorithm, our approach generates seed points purely from 3D line segments instead of using those feature points. First, 3D line segments are extracted using the existing Line3D++ algorithm, and the 3D line segment clustering criterion of buildings is established based on Manhattan world assumption. Next, by using the normal direction obtained using the result of 3D line segment clustering, we propose a multihypothetical plane fitting algorithm based on the mean shift method. Then, through subdividing on the triangle mesh constructed based on the building hypothetical plane model, semidense point cloud can be quickly obtained, and it is used as seed points of the PMVS pipeline instead of the sparse and noisy seed points generated by PMVS itself. After that, dense point cloud can be obtained through the existing PMVS expansion pipeline. Finally, unit and integration experiments are designed; the test results show that the proposed algorithm is 15%∼23% faster than the original PMWS in running time, and at the same time, the reconstruction quality of buildings is improved as well by successfully removing many noise points in the buildings.


---

## Body

## 1. Introduction

Accurate spatiotemporal information is the fundamental guarantee of the implementation of smart cities. The digital building model is one of the core elements of spatiotemporal information, and its automatic 3D reconstruction technology has received significant attentions in the field of geomatics and computer vision [1, 2]. Although several kinds of data sources could be used to reconstruct buildings such as images and LiDAR, image-based building reconstruction has become a research hotspot for it is able to provide very fine building models with relatively low cost [2–6].Among those existing image-based building reconstruction algorithms, structure from motion (SfM) and multiview stereo (MVS) are able to reconstruct 3D point clouds from multiple building images. As the most frequently used MVS algorithm, patch-based multiview stereo (PMVS) proposed by Furukawa is able to reconstruct the dense point cloud building model with very good performance [3, 7, 8].The advantages of PMVS algorithms include the following: (1) no need of bounding volume prior such as convex hull and bounding box; (2) being able to reconstruct both single object and crowded scenes; and (3) having high accuracy and efficiency [8]. However, as many researchers mentioned [4, 5, 9, 10], the main disadvantages of PMVS are the following: (1) the obtained normal information of the reconstructed points is not well consistent with its local geometry, and the problem becomes more obvious when the image capturing configuration is downward shooting or upward shooting; (2) the algorithm itself is with high time and memory complexities, which makes its time and memory consumption unacceptable when processing large-scale high-resolution images; and (3) it sometimes fails to reconstruct the complete point cloud of objects with poor-textured regions and non-Lambertian planes.Recently, some studies were conducted to improve the building reconstruction technology. Shi et al. [9] proposed an improved strategy based on a patch adjusting trick through using the scene geometric information and a multiresolution expanding approach, which improves the normal vector precision and surface smoothness. Ying et al. [11] defined matching strategy constraints in the PMVS feature matching process; as a result, the precision and integrity of the reconstructed point cloud are improved. Wu et al. [12] used the tensor to improve the PMVS algorithm, through fully using the surface geometry information, and the normal vector of PMVS reconstruction point cloud is effectively improved. Fei et al. [5] presented a novel object-based MVS algorithm using structural similarity (SSIM) index matching cost in a coarse to fine workflow, and the efficiency can be greatly improved when the overlap between images is large. Aiming at large UAV image processing in narrow baseline cases, Xiao et al. [4] proposed an image grouping and self-adaptive patch-based multiview stereo matching algorithm aided by DEM information, making the processing efficiency significantly better than that of the PMVS algorithm. For generating accurate large-scale 3D models of buildings employing point clouds derived from UAV-based photogrammetry, a new two-level segmentation approach based on efficient RANSAC shape detection was developed to segment potential facades and roofs of the buildings and extract their footprints [6]. Li et al. [13] used the Canny boundary detection operator to improve the PMVS algorithm and obtained a better reconstruction result. Considering the difficulties in the reconstruction of textureless surface under Manhattan world assumption, Furukawa et al. [14] proposed a method using the point cloud reconstructed by PMVS to extract the hypothetical plane and using the MRF model to reconstruct the depth map of each photo, and this method provides a better performance in poor-texture condition. Zhu and Leow [10] presented an alternative approach for reconstructing textured mesh surfaces from point cloud recovered by using the patch-based MVS method, and their test results showed that the reconstructed 3D models are sufficiently accurate and realistic for 3D visualization in various applications. Li et al. [15] defined a texture probabilistic grammar and developed an algorithm to obtain the 3D information in a 2D scene by training the texture probabilistic grammar from the prebuilt model library, and the algorithm was reported to be superior to the traditional reconstruction method based on point clouds.The above approaches are mostly improved through enhancing the precision and the normal vectors of point clouds, while few studies have been performed on the improvement of efficiency, and the result in the research by Xiao et al. [4] could accelerate the process time, however, it depends on DEM data, which could not be directly used in building reconstruction. To adapt to the massive reconstruction data and rapid data updating of smart cities, the reconstruction efficiency of the PMVS algorithm remains to be improved [9].Basically, the PMVS algorithm uses weak constraints to perform image feature matching and reconstructs the sparse point cloud [3], which causes local optimal results and the reconstructed sparse point cloud to be contaminated by a lot of noise points. When such noise-rich sparse points are used as the seed points for the subsequent patch expansion and optimization process, meaningless expansions could be executed, and the algorithm efficiency will be downgraded since the patch expansion and optimization process is the main time-consuming process in the PMVS algorithm. Therefore, theoretically, the PMVS algorithm could be improved through enhancing the accuracy and density of the seed points.In this study, an improved PMVS algorithm effectively using the spatial characteristics of buildings is proposed, which generates seed points for the PMVS pipeline directly using 3D line segment information. In Section3, a 3D line segment clustering algorithm based on the Manhattan world assumption is presented, which is able to determine the directions of three coordinate axes through clustering line segments into two groups, axis directional and nonaxis directional, quickly and accurately. Then, the hypothetical plane fitting algorithm based on 3D line sets is presented to obtain the building hypothetical plane model. After that, in Section 4, an initial triangle mesh is constructed using the point cloud sampled from 3D line segments for each hypothetical plane, and semidense seed points can be obtained through triangle mesh subdivision with a penalty mechanism; consequently, a dense point cloud model of the building is obtained through existing PMVS expansion and filtering pipeline. Finally, in Section 5, tests are designed to verify the performance of our algorithm, and the experiments show that the proposed algorithm is able to accelerate the building reconstruction speed of the PMVS, and our algorithm successfully removed many noise points that do not belong to buildings.


## 2. The Overall Working Flow of the Improved PMVS Algorithm

The overall working flow of the improved PMVS algorithm proposed in this study is described in Figure1.Figure 1
The overall working flow chart of the improved PMVS algorithm based on spatial plane geometric characteristics of buildings.In this method, the 3D line segments of the building are first identified and extracted using Line3D++ algorithms [16]; after that, the line segments are clustered using its direction information, and then, the hypothetical plane model of the building is fitted. After obtaining the hypothetical planes, the triangle meshes of the building are generated, and an accurate semidense point cloud can be quickly obtained by subdividing the triangle meshes with a penalty mechanism. Finally, the semidense point cloud is processed as the seed points of PMVS, and the dense point cloud model of the building could be constructed through the existing PMVS expansion and filtering pipeline.


## 3. Hypothetical Plane Fitting Algorithm Based on the 3D Line Segment Model of Buildings

The hypothetical planes are the potential planar faces in building. Regarding the data source, the existing building hypothetical plane fitting algorithm can be divided into two categories: (1) based on dense points [14], and the disadvantage is that the acquisition of dense points is a time-consuming process and (2) based on sparse points and sparse 3D line segments [17], and the corresponding issue of such algorithms is that the sparse point cloud reconstructed from SfM usually contains much noise, which may reduce the precision of the fitted hypothetical plane. In addition, SfM depends on feature point detection, and if there is textureless region or region with non-Lambertian reflectance in the building surfaces, it will be very difficult to obtain reliable sparse or dense points.Real-world buildings generally consist of several mutually perpendicular or parallel faces being rectangles or rectangle-like styles. Considering the above spatial geometric features, this paper proposes a hypothetical plane fitting algorithm based on the 3D line segment model of buildings. First, the Line3D++ algorithm by Hofer et al. is applied to reconstruct the 3D line segment model [16]; then, a 3D line segment clustering criterion under the Manhattan world assumption is proposed, and after that, the spatial multihypothesis planes are fitted using the mean shift method. The algorithm has higher reliability and efficiency in the complex reconstruction environment.


### 3.1. 3D Line Segment Clustering Algorithm Based on Manhattan World Assumption

Manhattan world assumption holds that the normal vectors of the planes of a building are usually parallel to one coordinate axis of the same Cartesian coordinate system [18], and it is widely used in building reconstruction and its postprocess [19–22]. A Manhattan world coordinate system is shown in Figure 2.Figure 2
A schematic diagram of the Manhattan world Cartesian coordinate system.Different from those previous research studies, we are using Manhattan world assumption to obtained seed points for the follow-up PMVS expansion and filtering pipeline. We present a 3D line segment clustering algorithm to quickly and accurately extract those axis-aligned and nonaxis-aligned line segment clusters, together with the three axis directions of the Cartesian coordinate system as well. As shown in Figure3, the main steps of the algorithm are as follows:Step 1. Perform preliminary clustering of the 3D line segments of the building using line direction information (refer to Algorithm 1 in the Supplementary Section), and then, the initial cluster set of the 3D line segments are obtained.Step 2. According to the perpendicular relationship among different line segment clusters, all cluster trios, which are perpendicular with each other, are extracted from the initial 3D line segment clusters, and a cost function is defined to find the optimal axis-aligned cluster trio, which has the greatest possibility of being parallel with the Cartesian coordinate axes (refer to Algorithm 2 in the Supplementary Section).Step 3. Extract the nonaxis-aligned clusters using the information of the optimal axis-aligned cluster trio.Step 4. Reclaim those missed line segments during the above clustering process by performing the clustering algorithm using the included angle information between the missed line segment and the cluster direction.Step 5. Filter all clustered line segments and line segment clusters according to the length of line segments and the size of line segment clusters. As a result, those line segments with their length being smaller than certain threshold will be treated as noise and be removed.Figure 3
The working flow of the 3D line segment clustering algorithm based on Manhattan world assumption.The detailed description of Step 1 is as follows. In the process of preliminary clustering of 3D line segments, we defineU as the unclassified line segment set and C as the classified line segment cluster set. First, put all the line segments into U in descending order according to their length. For each element in U, check whether it belongs to any existing cluster according to its included angle and length, and if it does not belong to any existing cluster, generate a new cluster if its length is greater than the threshold value, or otherwise, put it back to U waiting for the next loop operation. The threshold value is the average length of all the line segments in the 3D line segment model, and constrained by that threshold value, the length of the first line segment of the newly added line segment clustering can be avoided to be too short, which ensures the accuracy of the newly added clustering, and is conducive to improving the robustness of the clustering algorithm. Repeat the above operations until U is empty or the number of cycles reaches the upper limit to obtain the initial 3D line segment cluster set C.The detailed description of Step 2 is as follows. In the extraction process of the three axis-aligned line segment clusters, first, the main direction of each initial 3D line segment cluster is calculated; then, the cluster pair whose main directions are perpendicular to each other are identified, and for each perpendicular cluster pair, find the cluster in the remaining cluster set that forms a Cartesian coordinate system; as a result, a candidate axis cluster trio is formed. For each group of candidate axis cluster trio, the cost functionW1 shown in (1) is established to evaluate its credibility, and finally, the optimal axis-aligned cluster trio O0, which is with the greatest W1 value, is obtained from all the candidate axis cluster trios:(1)W1=∑a=1n1La∗maxDk1⟶·la⟶,Dk2⟶·la⟶,Dk3⟶·la⟶∑a=1n1La,where La represents a line segment in set C, Ia⟶ represents the direction of La, La represents the length of La, n1 represents the total number of line segments in set C, and Dk1⟶, Dk2⟶, and Dk3⟶, respectively, represent the main directions of the three axis-aligned line segment clusters in the kth candidate axis cluster trio. The cost function W1 actually tells the number of line segments being parallel to the three axes of the Cartesian coordinate system of the current 3D line segment cluster trio.The extraction process of nonaxis-aligned line segment cluster follows the determination process of the optimal axis cluster trioO0. For each line segment cluster, calculate the included angles between its main direction and the directions of three axis-aligned clustering in O0; if all three angles are greater than a threshold value (the angle is negatively correlated with the direction consistency of line segments in the axis-aligned line segment cluster, and the recommended threshold value is 20°, which is obtained by experiment), the line segment cluster is determined to be a nonaxis-aligned one and is added to the nonaxis-aligned line segment cluster set N. After second clustering and filtering procedure, the final line segment cluster set B can be obtained.


### 3.2. A Multihypothetical Plane Fitting Algorithm Based on the Mean Shift Method

As aforementioned, real-world buildings generally consist of several mutually perpendicular or parallel planes; if we successfully get all those planes’ model, the building could be represented using aB-rep model. To fit planes on a building, we use a spatial multihypothetical plane fitting algorithm based on the mean shift algorithm [23]. First, we generate the sparse point cloud through equidistance point sampling from the 3D line segments, then determine the normal vector of the hypothetical plane according to the main direction of the line segment cluster, and after that, we use the mean shift algorithm to cluster the point cloud and get the spatial multiplane point cloud clusters along the direction of each normal vector. The spatial point cloud cluster of each plane represents the hypothetical plane that needs to be fitted.


#### 3.2.1. Fitting of the Axis-Aligned Plane

The axis-aligned direction is defined to be parallel to one of the three coordinate axes in the Manhattan world Cartesian coordinate system. Based on the Manhattan world assumption, the normal vectors of most building planes are parallel to one of the three coordinate axes; so, the most hypothetical planes of buildings can be extracted by performing spatial point cloud clustering along the three axis-aligned directions. The specific algorithm flow is shown in Figure4, and the steps are as follows:Step 1. Generate point cloud: generate sparse point cloud V1 through equidistance sampling on each line segment in axis-aligned clusters; equidistance sampling means dividing each line segment in axis-aligned clusters O0 into several parts at equal distance d to get all end points of the subline segment, and d is related to the expected size of point cloud for plane fitting.Step 2. Transform coordinates: after obtaining the sparse point cloud V1, for the 3D coordinate vb⟶ of each 3D point, use the formula xb=vb⟶⋅Dk⟶ to convert the coordinates of the point from 3D coordinate to 1D, where Dk⟶ indicates the axis direction of the plane to be fitted to.Step 3. Fit planes based on the mean shift method: the mean shift clustering method is performed on the 1D point cloud obtained by coordinate transformation of V1, and multiple initial point cloud clusters are gotten. Each initial point cloud cluster represents an initial hypothetical plane.Step 4. Filter hypothetical planes: a cost function W2 is defined in (2) to evaluate the credibility of the initial hypothetical plane Pi, and those low-credibility ones are filtered out:(2)W2=n2+∑a=1n2LaLmin,wheren2 is the number of 3D line segments contained in Pi, La represents the length of the ath 3D line segment in Pi, and Lmin represents the length of the shortest line segment in O0. The cost function considers both the number and the length of the line segment in the initial hypothetical plane Pi, and more 3D line segments and longer segments usually mean the higher probability that the initial hypothetical plane represents the real plane. Generally, Pi would be adopted if the following condition is satisfied:(3)W2i>β1×W2max,whereW2i is the W2 value of Pi, W2max is the maximum value of W2 of all initial hypothetical planes, β1 is related to the quantity and precision specified for extracting hypothetical plane, and the recommended value for β1 is 50%.Step 5. Repeat steps 2 to 4 on all axis directions, and the final axis-aligned hypothetical plane fitting result is obtained.Figure 4
The working flow of the plane fitting algorithm along with axis-aligned directions.


#### 3.2.2. Nonaxis-Aligned Plane Fitting

Unlike plane fitting in axis-aligned direction, it is difficult to determine the normal vector of those nonaxis-aligned planes. As shown in Figure5, for the same nonaxis-aligned line segment cluster, there are multiple possibilities of the plane clustering direction (or the normal vector direction of the hypothetical plane). From the information provided by the existing line segment clusters, we can determine the direction of the plane clustering and find the most possible hypothetical plane distribution. Suppose H is the set of nonaxis-aligned line segment, and F is the final set of hypothetical planes. The mean shift algorithm is selected in the paper to conduct plane fitting due to its three advantages: (1) it does not need any prior information about the characteristics of the data model and can be applied to clustering of any feature space; (2) it only need one parameter (bandwidth of the kernel function) as input; and (3) it adopts statistical features and has strong robustness to noise. Figure 6 is the process sequence diagram, and the detailed description (refer to Algorithm 3 in the Supplementary Section) is as follows:Step 1. Initialize the set H of line segments to be clustered: take out each member from set N, the nonaxis-aligned line segment cluster, and add all its line segments to H.Step 2. Generate the sparse point cloud V2: generate sparse point cloud V2 from a 3D line segment in H by means of equidistance point sampling.Step 3. Determine the direction of plane fitting: randomly pick two line segments Lp and Lq, compute the normal vector n⟶ of the corresponding plane S, and decide whether n⟶ is the valid plane clustering direction by counting the number of line segments embedded in plane S in H. If the number is greater than given threshold (the recommended value is 3), n⟶ will be determined to be the effective plane clustering direction, and the subsequent plane fitting operation will be performed; otherwise, two new line segments will be randomly picked from H again, and repeat aforementioned operation.Step 4. Fit planes based on the mean shift method: if n⟶ is the valid clustering direction, perform mean shift clustering on point cloud V2 to obtain the point cloud cluster set Q, which will be used to fit the corresponding hypothetical planes.Step 5. Filter hypothetical planes: for each hypothetical plane in Q, if its line segment number is β2 times greater than total line segment number in H, it will be treated as valid and be added to F, and all corresponding line segments will be removed from H, where the recommended value for β2 is 50%. In addition, those left invalid hypothetical planes in Q will be kept first and processed later.Step 6. Repeat steps 2 to 5 until H is empty or the cycle number reaches the upper limit.Step 7. Find the optimal hypothetical plane cluster: suppose after r cycles of operations from steps 2 to 5, there are r hypothetical plane clusters in total, and Qj represents the jth hypothetical plane. We define a cost function W3 shown in (4) to evaluate the credibility of each hypothetical plane cluster:(4)W3=∑i=1mW2im,whereW2i is the W2 value of the ith hypothetical plane in Qj, and m is the total hypothetical plane number in Qj. From the r hypothetical plane clusters, select the one with the largest W3 as the optimal hypothetical plane cluster and name it as Q′.Step 8. Filter hypothetical planes: filter Q′ using cost function W2 defined in (2), and add Q′ to the set F.Step 9. Traverse each member in the nonaxis-aligned line segment cluster set N, repeat steps 1 to 8 to get the final nonaxis-aligned hypothetical plane fitting result.Figure 5
The ambiguity of the nonaxis-aligned plane fitting, the normal direction.Figure 6
The working flow of the nonaxis-aligned plane fitting algorithm.Additionally, the criterion is that the optimal hypothetical plane cluster is the one with the least number of hypothetical planes, and the maximum sum of hypothetical planes credibility is used to solve the problem of direction ambiguity of hypothetical plane fitting in the nonaxis-aligned direction. Now, the building hypothetical plane model is obtained and can provide reliable building’s plane information for the subsequent triangle mesh subdivision algorithm.


## 3.1. 3D Line Segment Clustering Algorithm Based on Manhattan World Assumption

Manhattan world assumption holds that the normal vectors of the planes of a building are usually parallel to one coordinate axis of the same Cartesian coordinate system [18], and it is widely used in building reconstruction and its postprocess [19–22]. A Manhattan world coordinate system is shown in Figure 2.Figure 2
A schematic diagram of the Manhattan world Cartesian coordinate system.Different from those previous research studies, we are using Manhattan world assumption to obtained seed points for the follow-up PMVS expansion and filtering pipeline. We present a 3D line segment clustering algorithm to quickly and accurately extract those axis-aligned and nonaxis-aligned line segment clusters, together with the three axis directions of the Cartesian coordinate system as well. As shown in Figure3, the main steps of the algorithm are as follows:Step 1. Perform preliminary clustering of the 3D line segments of the building using line direction information (refer to Algorithm 1 in the Supplementary Section), and then, the initial cluster set of the 3D line segments are obtained.Step 2. According to the perpendicular relationship among different line segment clusters, all cluster trios, which are perpendicular with each other, are extracted from the initial 3D line segment clusters, and a cost function is defined to find the optimal axis-aligned cluster trio, which has the greatest possibility of being parallel with the Cartesian coordinate axes (refer to Algorithm 2 in the Supplementary Section).Step 3. Extract the nonaxis-aligned clusters using the information of the optimal axis-aligned cluster trio.Step 4. Reclaim those missed line segments during the above clustering process by performing the clustering algorithm using the included angle information between the missed line segment and the cluster direction.Step 5. Filter all clustered line segments and line segment clusters according to the length of line segments and the size of line segment clusters. As a result, those line segments with their length being smaller than certain threshold will be treated as noise and be removed.Figure 3
The working flow of the 3D line segment clustering algorithm based on Manhattan world assumption.The detailed description of Step 1 is as follows. In the process of preliminary clustering of 3D line segments, we defineU as the unclassified line segment set and C as the classified line segment cluster set. First, put all the line segments into U in descending order according to their length. For each element in U, check whether it belongs to any existing cluster according to its included angle and length, and if it does not belong to any existing cluster, generate a new cluster if its length is greater than the threshold value, or otherwise, put it back to U waiting for the next loop operation. The threshold value is the average length of all the line segments in the 3D line segment model, and constrained by that threshold value, the length of the first line segment of the newly added line segment clustering can be avoided to be too short, which ensures the accuracy of the newly added clustering, and is conducive to improving the robustness of the clustering algorithm. Repeat the above operations until U is empty or the number of cycles reaches the upper limit to obtain the initial 3D line segment cluster set C.The detailed description of Step 2 is as follows. In the extraction process of the three axis-aligned line segment clusters, first, the main direction of each initial 3D line segment cluster is calculated; then, the cluster pair whose main directions are perpendicular to each other are identified, and for each perpendicular cluster pair, find the cluster in the remaining cluster set that forms a Cartesian coordinate system; as a result, a candidate axis cluster trio is formed. For each group of candidate axis cluster trio, the cost functionW1 shown in (1) is established to evaluate its credibility, and finally, the optimal axis-aligned cluster trio O0, which is with the greatest W1 value, is obtained from all the candidate axis cluster trios:(1)W1=∑a=1n1La∗maxDk1⟶·la⟶,Dk2⟶·la⟶,Dk3⟶·la⟶∑a=1n1La,where La represents a line segment in set C, Ia⟶ represents the direction of La, La represents the length of La, n1 represents the total number of line segments in set C, and Dk1⟶, Dk2⟶, and Dk3⟶, respectively, represent the main directions of the three axis-aligned line segment clusters in the kth candidate axis cluster trio. The cost function W1 actually tells the number of line segments being parallel to the three axes of the Cartesian coordinate system of the current 3D line segment cluster trio.The extraction process of nonaxis-aligned line segment cluster follows the determination process of the optimal axis cluster trioO0. For each line segment cluster, calculate the included angles between its main direction and the directions of three axis-aligned clustering in O0; if all three angles are greater than a threshold value (the angle is negatively correlated with the direction consistency of line segments in the axis-aligned line segment cluster, and the recommended threshold value is 20°, which is obtained by experiment), the line segment cluster is determined to be a nonaxis-aligned one and is added to the nonaxis-aligned line segment cluster set N. After second clustering and filtering procedure, the final line segment cluster set B can be obtained.


## 3.2. A Multihypothetical Plane Fitting Algorithm Based on the Mean Shift Method

As aforementioned, real-world buildings generally consist of several mutually perpendicular or parallel planes; if we successfully get all those planes’ model, the building could be represented using aB-rep model. To fit planes on a building, we use a spatial multihypothetical plane fitting algorithm based on the mean shift algorithm [23]. First, we generate the sparse point cloud through equidistance point sampling from the 3D line segments, then determine the normal vector of the hypothetical plane according to the main direction of the line segment cluster, and after that, we use the mean shift algorithm to cluster the point cloud and get the spatial multiplane point cloud clusters along the direction of each normal vector. The spatial point cloud cluster of each plane represents the hypothetical plane that needs to be fitted.


### 3.2.1. Fitting of the Axis-Aligned Plane

The axis-aligned direction is defined to be parallel to one of the three coordinate axes in the Manhattan world Cartesian coordinate system. Based on the Manhattan world assumption, the normal vectors of most building planes are parallel to one of the three coordinate axes; so, the most hypothetical planes of buildings can be extracted by performing spatial point cloud clustering along the three axis-aligned directions. The specific algorithm flow is shown in Figure4, and the steps are as follows:Step 1. Generate point cloud: generate sparse point cloud V1 through equidistance sampling on each line segment in axis-aligned clusters; equidistance sampling means dividing each line segment in axis-aligned clusters O0 into several parts at equal distance d to get all end points of the subline segment, and d is related to the expected size of point cloud for plane fitting.Step 2. Transform coordinates: after obtaining the sparse point cloud V1, for the 3D coordinate vb⟶ of each 3D point, use the formula xb=vb⟶⋅Dk⟶ to convert the coordinates of the point from 3D coordinate to 1D, where Dk⟶ indicates the axis direction of the plane to be fitted to.Step 3. Fit planes based on the mean shift method: the mean shift clustering method is performed on the 1D point cloud obtained by coordinate transformation of V1, and multiple initial point cloud clusters are gotten. Each initial point cloud cluster represents an initial hypothetical plane.Step 4. Filter hypothetical planes: a cost function W2 is defined in (2) to evaluate the credibility of the initial hypothetical plane Pi, and those low-credibility ones are filtered out:(2)W2=n2+∑a=1n2LaLmin,wheren2 is the number of 3D line segments contained in Pi, La represents the length of the ath 3D line segment in Pi, and Lmin represents the length of the shortest line segment in O0. The cost function considers both the number and the length of the line segment in the initial hypothetical plane Pi, and more 3D line segments and longer segments usually mean the higher probability that the initial hypothetical plane represents the real plane. Generally, Pi would be adopted if the following condition is satisfied:(3)W2i>β1×W2max,whereW2i is the W2 value of Pi, W2max is the maximum value of W2 of all initial hypothetical planes, β1 is related to the quantity and precision specified for extracting hypothetical plane, and the recommended value for β1 is 50%.Step 5. Repeat steps 2 to 4 on all axis directions, and the final axis-aligned hypothetical plane fitting result is obtained.Figure 4
The working flow of the plane fitting algorithm along with axis-aligned directions.


### 3.2.2. Nonaxis-Aligned Plane Fitting

Unlike plane fitting in axis-aligned direction, it is difficult to determine the normal vector of those nonaxis-aligned planes. As shown in Figure5, for the same nonaxis-aligned line segment cluster, there are multiple possibilities of the plane clustering direction (or the normal vector direction of the hypothetical plane). From the information provided by the existing line segment clusters, we can determine the direction of the plane clustering and find the most possible hypothetical plane distribution. Suppose H is the set of nonaxis-aligned line segment, and F is the final set of hypothetical planes. The mean shift algorithm is selected in the paper to conduct plane fitting due to its three advantages: (1) it does not need any prior information about the characteristics of the data model and can be applied to clustering of any feature space; (2) it only need one parameter (bandwidth of the kernel function) as input; and (3) it adopts statistical features and has strong robustness to noise. Figure 6 is the process sequence diagram, and the detailed description (refer to Algorithm 3 in the Supplementary Section) is as follows:Step 1. Initialize the set H of line segments to be clustered: take out each member from set N, the nonaxis-aligned line segment cluster, and add all its line segments to H.Step 2. Generate the sparse point cloud V2: generate sparse point cloud V2 from a 3D line segment in H by means of equidistance point sampling.Step 3. Determine the direction of plane fitting: randomly pick two line segments Lp and Lq, compute the normal vector n⟶ of the corresponding plane S, and decide whether n⟶ is the valid plane clustering direction by counting the number of line segments embedded in plane S in H. If the number is greater than given threshold (the recommended value is 3), n⟶ will be determined to be the effective plane clustering direction, and the subsequent plane fitting operation will be performed; otherwise, two new line segments will be randomly picked from H again, and repeat aforementioned operation.Step 4. Fit planes based on the mean shift method: if n⟶ is the valid clustering direction, perform mean shift clustering on point cloud V2 to obtain the point cloud cluster set Q, which will be used to fit the corresponding hypothetical planes.Step 5. Filter hypothetical planes: for each hypothetical plane in Q, if its line segment number is β2 times greater than total line segment number in H, it will be treated as valid and be added to F, and all corresponding line segments will be removed from H, where the recommended value for β2 is 50%. In addition, those left invalid hypothetical planes in Q will be kept first and processed later.Step 6. Repeat steps 2 to 5 until H is empty or the cycle number reaches the upper limit.Step 7. Find the optimal hypothetical plane cluster: suppose after r cycles of operations from steps 2 to 5, there are r hypothetical plane clusters in total, and Qj represents the jth hypothetical plane. We define a cost function W3 shown in (4) to evaluate the credibility of each hypothetical plane cluster:(4)W3=∑i=1mW2im,whereW2i is the W2 value of the ith hypothetical plane in Qj, and m is the total hypothetical plane number in Qj. From the r hypothetical plane clusters, select the one with the largest W3 as the optimal hypothetical plane cluster and name it as Q′.Step 8. Filter hypothetical planes: filter Q′ using cost function W2 defined in (2), and add Q′ to the set F.Step 9. Traverse each member in the nonaxis-aligned line segment cluster set N, repeat steps 1 to 8 to get the final nonaxis-aligned hypothetical plane fitting result.Figure 5
The ambiguity of the nonaxis-aligned plane fitting, the normal direction.Figure 6
The working flow of the nonaxis-aligned plane fitting algorithm.Additionally, the criterion is that the optimal hypothetical plane cluster is the one with the least number of hypothetical planes, and the maximum sum of hypothetical planes credibility is used to solve the problem of direction ambiguity of hypothetical plane fitting in the nonaxis-aligned direction. Now, the building hypothetical plane model is obtained and can provide reliable building’s plane information for the subsequent triangle mesh subdivision algorithm.


## 3.2.1. Fitting of the Axis-Aligned Plane

The axis-aligned direction is defined to be parallel to one of the three coordinate axes in the Manhattan world Cartesian coordinate system. Based on the Manhattan world assumption, the normal vectors of most building planes are parallel to one of the three coordinate axes; so, the most hypothetical planes of buildings can be extracted by performing spatial point cloud clustering along the three axis-aligned directions. The specific algorithm flow is shown in Figure4, and the steps are as follows:Step 1. Generate point cloud: generate sparse point cloud V1 through equidistance sampling on each line segment in axis-aligned clusters; equidistance sampling means dividing each line segment in axis-aligned clusters O0 into several parts at equal distance d to get all end points of the subline segment, and d is related to the expected size of point cloud for plane fitting.Step 2. Transform coordinates: after obtaining the sparse point cloud V1, for the 3D coordinate vb⟶ of each 3D point, use the formula xb=vb⟶⋅Dk⟶ to convert the coordinates of the point from 3D coordinate to 1D, where Dk⟶ indicates the axis direction of the plane to be fitted to.Step 3. Fit planes based on the mean shift method: the mean shift clustering method is performed on the 1D point cloud obtained by coordinate transformation of V1, and multiple initial point cloud clusters are gotten. Each initial point cloud cluster represents an initial hypothetical plane.Step 4. Filter hypothetical planes: a cost function W2 is defined in (2) to evaluate the credibility of the initial hypothetical plane Pi, and those low-credibility ones are filtered out:(2)W2=n2+∑a=1n2LaLmin,wheren2 is the number of 3D line segments contained in Pi, La represents the length of the ath 3D line segment in Pi, and Lmin represents the length of the shortest line segment in O0. The cost function considers both the number and the length of the line segment in the initial hypothetical plane Pi, and more 3D line segments and longer segments usually mean the higher probability that the initial hypothetical plane represents the real plane. Generally, Pi would be adopted if the following condition is satisfied:(3)W2i>β1×W2max,whereW2i is the W2 value of Pi, W2max is the maximum value of W2 of all initial hypothetical planes, β1 is related to the quantity and precision specified for extracting hypothetical plane, and the recommended value for β1 is 50%.Step 5. Repeat steps 2 to 4 on all axis directions, and the final axis-aligned hypothetical plane fitting result is obtained.Figure 4
The working flow of the plane fitting algorithm along with axis-aligned directions.


## 3.2.2. Nonaxis-Aligned Plane Fitting

Unlike plane fitting in axis-aligned direction, it is difficult to determine the normal vector of those nonaxis-aligned planes. As shown in Figure5, for the same nonaxis-aligned line segment cluster, there are multiple possibilities of the plane clustering direction (or the normal vector direction of the hypothetical plane). From the information provided by the existing line segment clusters, we can determine the direction of the plane clustering and find the most possible hypothetical plane distribution. Suppose H is the set of nonaxis-aligned line segment, and F is the final set of hypothetical planes. The mean shift algorithm is selected in the paper to conduct plane fitting due to its three advantages: (1) it does not need any prior information about the characteristics of the data model and can be applied to clustering of any feature space; (2) it only need one parameter (bandwidth of the kernel function) as input; and (3) it adopts statistical features and has strong robustness to noise. Figure 6 is the process sequence diagram, and the detailed description (refer to Algorithm 3 in the Supplementary Section) is as follows:Step 1. Initialize the set H of line segments to be clustered: take out each member from set N, the nonaxis-aligned line segment cluster, and add all its line segments to H.Step 2. Generate the sparse point cloud V2: generate sparse point cloud V2 from a 3D line segment in H by means of equidistance point sampling.Step 3. Determine the direction of plane fitting: randomly pick two line segments Lp and Lq, compute the normal vector n⟶ of the corresponding plane S, and decide whether n⟶ is the valid plane clustering direction by counting the number of line segments embedded in plane S in H. If the number is greater than given threshold (the recommended value is 3), n⟶ will be determined to be the effective plane clustering direction, and the subsequent plane fitting operation will be performed; otherwise, two new line segments will be randomly picked from H again, and repeat aforementioned operation.Step 4. Fit planes based on the mean shift method: if n⟶ is the valid clustering direction, perform mean shift clustering on point cloud V2 to obtain the point cloud cluster set Q, which will be used to fit the corresponding hypothetical planes.Step 5. Filter hypothetical planes: for each hypothetical plane in Q, if its line segment number is β2 times greater than total line segment number in H, it will be treated as valid and be added to F, and all corresponding line segments will be removed from H, where the recommended value for β2 is 50%. In addition, those left invalid hypothetical planes in Q will be kept first and processed later.Step 6. Repeat steps 2 to 5 until H is empty or the cycle number reaches the upper limit.Step 7. Find the optimal hypothetical plane cluster: suppose after r cycles of operations from steps 2 to 5, there are r hypothetical plane clusters in total, and Qj represents the jth hypothetical plane. We define a cost function W3 shown in (4) to evaluate the credibility of each hypothetical plane cluster:(4)W3=∑i=1mW2im,whereW2i is the W2 value of the ith hypothetical plane in Qj, and m is the total hypothetical plane number in Qj. From the r hypothetical plane clusters, select the one with the largest W3 as the optimal hypothetical plane cluster and name it as Q′.Step 8. Filter hypothetical planes: filter Q′ using cost function W2 defined in (2), and add Q′ to the set F.Step 9. Traverse each member in the nonaxis-aligned line segment cluster set N, repeat steps 1 to 8 to get the final nonaxis-aligned hypothetical plane fitting result.Figure 5
The ambiguity of the nonaxis-aligned plane fitting, the normal direction.Figure 6
The working flow of the nonaxis-aligned plane fitting algorithm.Additionally, the criterion is that the optimal hypothetical plane cluster is the one with the least number of hypothetical planes, and the maximum sum of hypothetical planes credibility is used to solve the problem of direction ambiguity of hypothetical plane fitting in the nonaxis-aligned direction. Now, the building hypothetical plane model is obtained and can provide reliable building’s plane information for the subsequent triangle mesh subdivision algorithm.


## 4. An Improved Expanding Algorithm Based on Triangle Mesh Subdivision

Considering the building’s spatial and geometry characteristics, we improve the expanding mechanism of the PMVS algorithm based on triangle mesh subdivision. By using the information of the hypothetical plane model of the building, semidense point cloud is obtained from spatial subdivision of the triangle meshes of the building and is used to replace sparse point cloud generated by the original PMVS algorithm working as the expanding seed points in the PMVS pipeline, so that the efficiency of the PMVS algorithm can be improved.


### 4.1. Working Flow

The point cloud cluster in the hypothetical plane obtained from the above hypothetical plane fitting algorithm is used as the input data, and the working flow of the improved algorithm based on triangle mesh subdivision is shown in Figure7. The detailed steps are as follows:Step 1. Construct the initial Delaunay triangle meshesStep 2. Subdivide and refine the initial modelStep 3. Use the Nelder–Mead nonlinear optimization algorithm to optimize and filter the semidense point cloud through patch optimization and filteringStep 4. Reconstruct the triangle meshes of the point cloud and update the normal informationStep 5. Repeat steps 2 to 4 within n1 times, and a recommended value for n1 is 3Step 6. Set the optimized semidense point cloud as the seed points of the PMVS algorithm, repeat the expansion and filtering process defined in the original PMVS pipeline for n2 times, and the final dense point cloud is obtained, where a recommended value for n2 is 2Figure 7
Working flow of the improved PMVS algorithm based on spatial subdivision triangular mesh.


### 4.2. Building the Initial Triangle Meshes

To implement subdivision, the initial triangle mesh should be obtained first. For each hypothetical plane, the triangle mesh is reconstructed using the GPU-based 2D constrained Delaunay triangulation algorithm [24] with line segments derived from the 3D line segment model as constraint. And mesh information needs to be initialized after the reconstruction, including the normal vector, visibility information, and subdividing penalty value. The normal vector is the average unit normal vector of all triangles planes, which contain this vertex. The visibility information is a set of images in which the vertex is visible, which is introduced by PMVS [3]. The subdividing penalty values are initialized to 0, which will be described in detail in Section 4.5.2.


### 4.3. Subdivision of Triangle Meshes

In the subdivision process, the area of the triangle patch is used to determine whether it needs subdivision. If its area is bigger than a certain threshold, the center point of the triangle patch is added, and the original triangle patch is subdivided into three triangle patches. The process could be described in Algorithm 4 in the Supplementary Section. And the detailed steps are as follows:Step 1. Given the m1th of the triangle tm1 in the initial triangle mesh Ta, the area threshold can be computed as(5)s0=s¯3loop−f,wheres¯ represents the average triangle area of the triangle mesh, loop represents the cycle number of subdivision, and f represents the subdivision penalty value of tm1, which can be computed as(6)f=fm10+fm11+fm123,wherefm10, fm11, and fm12, respectively, represent subdivision penalty values of three vertices in tm1.Two setsUT and UT′ are used to represent the set of triangles to be subdivided and the set of temporary triangles generated during subdivision process, respectively, and compute the area sm1 of tm1; if sm1>s0, add tm1 to UT, or otherwise, repeat step 1 on the next triangle of Ta.Step 2. For the m2th of the triangle tm2 in the set UT, compute the center point vc of tm2, and add the center point vc to the point cloud Pa as a new expanded vertex. Update the normal vector of vc and the visibility information as follows: set the normal vector of vc to the normal vector of plane tm2; and set the visibility information of vc as vc=V0∪V1∪V2, where V0, V1, and V2 represent the visibility information of the three vertices of tm2, respectively.Step 3. After inserting point vc to the point cloud pa as the new expanded vertex, vc and the original three vertices of the triangle t(m2) will form three new triangles, and then add the three new triangles to the set UT′.Step 4. Repeat steps 2 to 4 for each triangle in the set UT, and exchange the set UT and UT′ after each cycle.Step 5. Repeat the step 4 until set UT is empty.Step 6. Repeat steps 1 to 5 for the next triangle until all triangles in Ta are traversed.


### 4.4. Patch Optimization and Filtering

By using the aforementioned triangle mesh subdivision algorithm, the initial triangle meshes are greatly refined. However, such subdivision is executed indiscriminately, and consequently noise vertices could be introduced. Therefore, additional optimization and filtering operations are needed to remove those noise points.In the original PMVS pipeline, three optimization parameters are used in the subdivision optimization process [3]: patch position cp⟶ and two Euler angles related to the normal vector of patch np⟶. This paper fully uses the geometric information provided by the triangle mesh of the hypothetical planes and reduces the parameters in the original nonlinear optimization of PMVS to just one parameter: the patch position cp⟶. And the normal vector of the patch remains unchanged during the optimization process. The patch position of the original PMVS algorithm is limited to the line which connects the center of the patch to the camera’s optical center corresponding to the reference photo. In our algorithm, the patch is restricted to move along the normal direction during the optimization process.Since the normal vector of the patch is regarded as a known parameter in this process, the number of optimization parameters is reduced, and the optimization speed can be therefore improved.


### 4.5. Triangle Mesh Updating Based on Subdivision Penalty Mechanism

#### 4.5.1. Update of Vertex Normal Vector

The point cloud of a given hypothetical plane has changed a lot after patch optimization and filtering. It is necessary to reconstruct the triangle mesh and update the normal vector information of all vertexes in the new triangle mesh using the average of all normal vectors of those triangles containing the vertex.


#### 4.5.2. Subdivision Penalty Mechanism

If we subdivide every triangle in the mesh using the subdivision strategy mentioned in Section 4.3 but without a penalty mechanism, each region of the hypothetical plane will almost have the same point cloud density. When there is a hole in the hypothetical plane, such an algorithm will not only downgrade the processing efficiency but also bring additional noise points.To overcome the problem, a penalty mechanism is introduced to determine whether to perform subdivision on certain triangle patch guided by the results of patch optimization and filtering to improve both the processing efficiency and the accuracy of the generated point cloud.SupposeTb is used to represent the new triangle mesh corresponding to the initial triangle mesh Ta after triangle mesh subdivision, patch optimization, and filtering. The process of the subdivision penalty mechanism could be described in Algorithm 5 in the Supplementary Section. And the detailed steps are as follows:Step 1. Compute the subdivision success rate r of all triangles in the initial triangle mesh Ta.The subdivision success rater is used to decide whether a triangle should be punished. For the ith of the triangle ti in Ta, the subdivision success rate is computed using(7)ri=nsuviventotal,wherentotal is the total number of vertices obtained after the subdivision of the triangle ti, and nsuvive indicates the number of vertices that the triangle ti contains after the improved patch optimization and filtering process. The subdivision success rate ri reflects the probability that ti is in a hole region: generally, the smaller ri is, the more probable ti is in a hole region. As a result, all vertices contained in ti should be punished, which makes them have low probability to be subdivided.Step 2. For the jth vertex vj in the new triangle mesh Tb, use fvj′ to represent its penalty value, and use fvj to represent its penalty in mesh Ta if it does not come from subdivision. If vj comes from subdivision, find the parent triangle tp of vj, and then go to Step 3; otherwise, fvj′=fvj, the penalty value of vj will not update, and then go to Step 4.Step 3. If the subdivision success rate of tp (the parent triangle of vj) satisfies the following conditions:(8)rtp≤β3,then thevj needs to be penalized, and we need modify the penalty value fvj′ of vj using (9); otherwise, just directly update it using (10):(9)fvj′=fu+ftp0+ftp1+ftp23,(10)fvj′=ftp0+ftp1+ftp23,whereftp0, ftp1, and ftp2, respectively, represent the subdivision penalty values of the three vertices of the parent triangle tp, fu is a preset penalty value, β3 is a preset subdivision success rate threshold, and the recommended value for β3 is 0.5.Step 4. Repeat steps 2 to 3 for each mesh vertex in Tb, and finally complete the update of the mesh vertex subdivision penalty value of Tb.


## 4.1. Working Flow

The point cloud cluster in the hypothetical plane obtained from the above hypothetical plane fitting algorithm is used as the input data, and the working flow of the improved algorithm based on triangle mesh subdivision is shown in Figure7. The detailed steps are as follows:Step 1. Construct the initial Delaunay triangle meshesStep 2. Subdivide and refine the initial modelStep 3. Use the Nelder–Mead nonlinear optimization algorithm to optimize and filter the semidense point cloud through patch optimization and filteringStep 4. Reconstruct the triangle meshes of the point cloud and update the normal informationStep 5. Repeat steps 2 to 4 within n1 times, and a recommended value for n1 is 3Step 6. Set the optimized semidense point cloud as the seed points of the PMVS algorithm, repeat the expansion and filtering process defined in the original PMVS pipeline for n2 times, and the final dense point cloud is obtained, where a recommended value for n2 is 2Figure 7
Working flow of the improved PMVS algorithm based on spatial subdivision triangular mesh.


## 4.2. Building the Initial Triangle Meshes

To implement subdivision, the initial triangle mesh should be obtained first. For each hypothetical plane, the triangle mesh is reconstructed using the GPU-based 2D constrained Delaunay triangulation algorithm [24] with line segments derived from the 3D line segment model as constraint. And mesh information needs to be initialized after the reconstruction, including the normal vector, visibility information, and subdividing penalty value. The normal vector is the average unit normal vector of all triangles planes, which contain this vertex. The visibility information is a set of images in which the vertex is visible, which is introduced by PMVS [3]. The subdividing penalty values are initialized to 0, which will be described in detail in Section 4.5.2.


## 4.3. Subdivision of Triangle Meshes

In the subdivision process, the area of the triangle patch is used to determine whether it needs subdivision. If its area is bigger than a certain threshold, the center point of the triangle patch is added, and the original triangle patch is subdivided into three triangle patches. The process could be described in Algorithm 4 in the Supplementary Section. And the detailed steps are as follows:Step 1. Given the m1th of the triangle tm1 in the initial triangle mesh Ta, the area threshold can be computed as(5)s0=s¯3loop−f,wheres¯ represents the average triangle area of the triangle mesh, loop represents the cycle number of subdivision, and f represents the subdivision penalty value of tm1, which can be computed as(6)f=fm10+fm11+fm123,wherefm10, fm11, and fm12, respectively, represent subdivision penalty values of three vertices in tm1.Two setsUT and UT′ are used to represent the set of triangles to be subdivided and the set of temporary triangles generated during subdivision process, respectively, and compute the area sm1 of tm1; if sm1>s0, add tm1 to UT, or otherwise, repeat step 1 on the next triangle of Ta.Step 2. For the m2th of the triangle tm2 in the set UT, compute the center point vc of tm2, and add the center point vc to the point cloud Pa as a new expanded vertex. Update the normal vector of vc and the visibility information as follows: set the normal vector of vc to the normal vector of plane tm2; and set the visibility information of vc as vc=V0∪V1∪V2, where V0, V1, and V2 represent the visibility information of the three vertices of tm2, respectively.Step 3. After inserting point vc to the point cloud pa as the new expanded vertex, vc and the original three vertices of the triangle t(m2) will form three new triangles, and then add the three new triangles to the set UT′.Step 4. Repeat steps 2 to 4 for each triangle in the set UT, and exchange the set UT and UT′ after each cycle.Step 5. Repeat the step 4 until set UT is empty.Step 6. Repeat steps 1 to 5 for the next triangle until all triangles in Ta are traversed.


## 4.4. Patch Optimization and Filtering

By using the aforementioned triangle mesh subdivision algorithm, the initial triangle meshes are greatly refined. However, such subdivision is executed indiscriminately, and consequently noise vertices could be introduced. Therefore, additional optimization and filtering operations are needed to remove those noise points.In the original PMVS pipeline, three optimization parameters are used in the subdivision optimization process [3]: patch position cp⟶ and two Euler angles related to the normal vector of patch np⟶. This paper fully uses the geometric information provided by the triangle mesh of the hypothetical planes and reduces the parameters in the original nonlinear optimization of PMVS to just one parameter: the patch position cp⟶. And the normal vector of the patch remains unchanged during the optimization process. The patch position of the original PMVS algorithm is limited to the line which connects the center of the patch to the camera’s optical center corresponding to the reference photo. In our algorithm, the patch is restricted to move along the normal direction during the optimization process.Since the normal vector of the patch is regarded as a known parameter in this process, the number of optimization parameters is reduced, and the optimization speed can be therefore improved.


## 4.5. Triangle Mesh Updating Based on Subdivision Penalty Mechanism

### 4.5.1. Update of Vertex Normal Vector

The point cloud of a given hypothetical plane has changed a lot after patch optimization and filtering. It is necessary to reconstruct the triangle mesh and update the normal vector information of all vertexes in the new triangle mesh using the average of all normal vectors of those triangles containing the vertex.


### 4.5.2. Subdivision Penalty Mechanism

If we subdivide every triangle in the mesh using the subdivision strategy mentioned in Section 4.3 but without a penalty mechanism, each region of the hypothetical plane will almost have the same point cloud density. When there is a hole in the hypothetical plane, such an algorithm will not only downgrade the processing efficiency but also bring additional noise points.To overcome the problem, a penalty mechanism is introduced to determine whether to perform subdivision on certain triangle patch guided by the results of patch optimization and filtering to improve both the processing efficiency and the accuracy of the generated point cloud.SupposeTb is used to represent the new triangle mesh corresponding to the initial triangle mesh Ta after triangle mesh subdivision, patch optimization, and filtering. The process of the subdivision penalty mechanism could be described in Algorithm 5 in the Supplementary Section. And the detailed steps are as follows:Step 1. Compute the subdivision success rate r of all triangles in the initial triangle mesh Ta.The subdivision success rater is used to decide whether a triangle should be punished. For the ith of the triangle ti in Ta, the subdivision success rate is computed using(7)ri=nsuviventotal,wherentotal is the total number of vertices obtained after the subdivision of the triangle ti, and nsuvive indicates the number of vertices that the triangle ti contains after the improved patch optimization and filtering process. The subdivision success rate ri reflects the probability that ti is in a hole region: generally, the smaller ri is, the more probable ti is in a hole region. As a result, all vertices contained in ti should be punished, which makes them have low probability to be subdivided.Step 2. For the jth vertex vj in the new triangle mesh Tb, use fvj′ to represent its penalty value, and use fvj to represent its penalty in mesh Ta if it does not come from subdivision. If vj comes from subdivision, find the parent triangle tp of vj, and then go to Step 3; otherwise, fvj′=fvj, the penalty value of vj will not update, and then go to Step 4.Step 3. If the subdivision success rate of tp (the parent triangle of vj) satisfies the following conditions:(8)rtp≤β3,then thevj needs to be penalized, and we need modify the penalty value fvj′ of vj using (9); otherwise, just directly update it using (10):(9)fvj′=fu+ftp0+ftp1+ftp23,(10)fvj′=ftp0+ftp1+ftp23,whereftp0, ftp1, and ftp2, respectively, represent the subdivision penalty values of the three vertices of the parent triangle tp, fu is a preset penalty value, β3 is a preset subdivision success rate threshold, and the recommended value for β3 is 0.5.Step 4. Repeat steps 2 to 3 for each mesh vertex in Tb, and finally complete the update of the mesh vertex subdivision penalty value of Tb.


## 4.5.1. Update of Vertex Normal Vector

The point cloud of a given hypothetical plane has changed a lot after patch optimization and filtering. It is necessary to reconstruct the triangle mesh and update the normal vector information of all vertexes in the new triangle mesh using the average of all normal vectors of those triangles containing the vertex.


## 4.5.2. Subdivision Penalty Mechanism

If we subdivide every triangle in the mesh using the subdivision strategy mentioned in Section 4.3 but without a penalty mechanism, each region of the hypothetical plane will almost have the same point cloud density. When there is a hole in the hypothetical plane, such an algorithm will not only downgrade the processing efficiency but also bring additional noise points.To overcome the problem, a penalty mechanism is introduced to determine whether to perform subdivision on certain triangle patch guided by the results of patch optimization and filtering to improve both the processing efficiency and the accuracy of the generated point cloud.SupposeTb is used to represent the new triangle mesh corresponding to the initial triangle mesh Ta after triangle mesh subdivision, patch optimization, and filtering. The process of the subdivision penalty mechanism could be described in Algorithm 5 in the Supplementary Section. And the detailed steps are as follows:Step 1. Compute the subdivision success rate r of all triangles in the initial triangle mesh Ta.The subdivision success rater is used to decide whether a triangle should be punished. For the ith of the triangle ti in Ta, the subdivision success rate is computed using(7)ri=nsuviventotal,wherentotal is the total number of vertices obtained after the subdivision of the triangle ti, and nsuvive indicates the number of vertices that the triangle ti contains after the improved patch optimization and filtering process. The subdivision success rate ri reflects the probability that ti is in a hole region: generally, the smaller ri is, the more probable ti is in a hole region. As a result, all vertices contained in ti should be punished, which makes them have low probability to be subdivided.Step 2. For the jth vertex vj in the new triangle mesh Tb, use fvj′ to represent its penalty value, and use fvj to represent its penalty in mesh Ta if it does not come from subdivision. If vj comes from subdivision, find the parent triangle tp of vj, and then go to Step 3; otherwise, fvj′=fvj, the penalty value of vj will not update, and then go to Step 4.Step 3. If the subdivision success rate of tp (the parent triangle of vj) satisfies the following conditions:(8)rtp≤β3,then thevj needs to be penalized, and we need modify the penalty value fvj′ of vj using (9); otherwise, just directly update it using (10):(9)fvj′=fu+ftp0+ftp1+ftp23,(10)fvj′=ftp0+ftp1+ftp23,whereftp0, ftp1, and ftp2, respectively, represent the subdivision penalty values of the three vertices of the parent triangle tp, fu is a preset penalty value, β3 is a preset subdivision success rate threshold, and the recommended value for β3 is 0.5.Step 4. Repeat steps 2 to 3 for each mesh vertex in Tb, and finally complete the update of the mesh vertex subdivision penalty value of Tb.


## 5. Experiment and Analysis

### 5.1. Test Datasets and Test Environment

To verify the performance of the algorithm presented in this paper, four sets of building sequence photos are used, which are named house1, house2, hall, and coffee shack, respectively. The detailed information such as the photo numbers and resolutions are shown in Table1, and the previews of the four datasets are shown in Figure 8.Table 1
Numbers of photo and resolution in experimental dataset1.
Dataset namePhoto numberResolutionHouse1652424 × 1351House2263071 × 2301Hall611534 × 1008Coffee shack301086 × 7151House1 is from https://github.com/zxg519/3D_reconstruction_dataset/, house2 is from https://demuc.de/colmap/datasets/, hall is from http://mirror.openstreetmap.nl/sources/pmvs-2/data/hall/, and coffee shack is from http://idav.ucdavis.edu/∼hkim/mvs/dataset/.Figure 8
Sequence image previews of the experimental datasets.The hardware configuration is as follows: CPU: Intel Core i5-6500 @ 3.2 GHz quad core; GPU: NVIDIA GeForce GTX 750, 2 GB RAM; and memory: 16G.The experiment carried out below can be divided into three parts. The first part is designed to verify the effectiveness of the hypothetical plane fitting algorithm, and the second part is designed to verify the effectiveness of triangle mesh reconstruction and subdivision. The third part runs the whole algorithm flow and compares the time consumption and the quality of point cloud reconstruction with the original PMVS algorithm.


### 5.2. Hypothetical Plane Fitting

#### 5.2.1. 3D Line Segment Clustering Algorithm Based on Manhattan World Assumption

Before the experiment, the COLMAP’s SfM algorithm [25] is used to perform camera calibration and sparse point cloud reconstruction for each experimental dataset; thus, the camera intrinsic and extrinsic parameters, initial feature point matching information, and sparse point cloud information are obtained. After that, Line3D++ [16], a 3D line segment model reconstruction algorithm, is used to obtain all 3D line segment information for each dataset.These 3D line segment models are used as the input data, and clustering operations are performed using a 3D line segment clustering algorithm based on Manhattan world assumption. The 3D line segment clustering results of the four datasets are shown in Figure9 (different types of line segments are represented using different colours in the figure), and the figure shows that the algorithm successfully extracts the main directions of the building line segment model, which includes the main axis-aligned and nonaxis-aligned directions. It could be found that(1)
It performs effectively on the nonaxis-aligned direction and works robustly on clustering of short line segments. As shown in Figure10, the nonaxis line segments in the datasets house1 and hall are still clustered correctly, although they are very short (about 1 m in Figure 10(a) and 0.2 m in Figure 10(b)), and the roof on the left side of the dataset house2 is not extracted just because the detected segments are too sparse.(2)
It works robustly when there is strong noise caused by natural landscapes. Most of the natural landscapes do not follow the Manhattan world hypothesis, and its reconstructed line segments can easily be filtered out because it is difficult to find their category in line segment clustering. Therefore, besides line segment clustering functionality, it also has the filtering function of line segment noise brought by natural landscape or nonartificial objects. In the dataset coffee shack, the trees introduce a large number of noise line segments to the Line3D++ reconstruction result, as shown in Figure11(a), while after using the line segment clustering algorithm proposed in this paper, most of these noise segments are filtered out, as shown in Figure 11(b).Figure 9
Experimental results of the 3D line segment clustering algorithm based on Manhattan world assumption, and different colours used to represent different line segment clusters. (a) Results of dataset house1; (b) results of dataset house2; (c) results of dataset hall; and (d) results of dataset coffee shack.
(a)(b)(c)(d)Figure 10
Results of the 3D line segment clustering algorithm in nonaxis-aligned directions, and different clusters are shown using different colours. (a) Results of house1 and (b) results of hall.
(a)(b)Figure 11
The comparison results of the 3D line segment model of coffee shack before and after line segment clustering show that the algorithm can filter line segment noise in a certain degree on the line segment noise generated by natural landscape. (a) Before clustering; and (b) after clustering.
(a)(b)


#### 5.2.2. Multihypothesis Plane Fitting Based on the Mean Shift Method

After obtaining the 3D line segment clusters of each dataset, the multihypothesis plane fitting algorithm based on the mean shift method is used to fit the point cloud sampled from the line segment model. The results of plane fitting are shown in Figures12 and 13. In the two figures, the hypothetical plane models of each dataset are presented in a point cloud clustering manner. Each point cloud cluster represents a hypothetical plane using different colours.Figure 12
The fitting results of hypothetical planes of different buildings where different hypothetical planes are represented using different colours (I). (a) House1 in the top view; (b) house1 in the oblique view; (c) house2 in the top view; and (d) house2 in the top view.
(a)(b)(c)(d)Figure 13
The fitting results of hypothetical planes of different buildings where different hypothetical planes are represented using different colours (II). (a) Hall in the front view; (b) hall in the top view; (c) coffee shack in the oblique view; and (d) coffee shack in the side view.
(a)(b)(c)(d)The fitting results show that the main planes of each building in the four datasets (including the main facades and the planes where the roof is located) are correctly extracted.The extraction completeness information of hypothetical planes is shown in Table2, where the completeness rate refers to the ratio of the number of final extracted hypothetical planes to the number of actual hypothetical planes. The actual hypothetical plane refers to the plane determined manually from the 3D line segment model of the dataset (tiny planes are not counted). Table 2 shows that our algorithm maintains a high completeness rate of hypothetical plane extraction for each dataset. Specially, the root causes of the five cases less than 1.00 are as follows: (1) there are too few reconstructed 3D line segments to successfully recover the plane information in the missed hypothetical plane; and (2) the quantitative area of the missed plane is too small or the line segments is too short.Table 2
Extraction completeness of hypothetical planes.
House1House2HallCoffee shackExtracted (X-axis)4263Actual (X-axis)4264Extracted (Y-axis)21172Actual (Y-axis)22202Extracted (Z-axis)0001Actual (Z-axis)0021Extracted (nonaxis)2422Actual (nonaxis)2423Although the algorithm may miss some hypothetical planes with small areas, the main hypothetical planes are correctly extracted, which will show little effect on the subsequent triangle mesh subdivision.The consumed time information of the test is shown in Figure14. The relationship between the consumed time and the total number of 3D line segments shows a linear-like relationship, and the time on line segment clustering is much less than that of plane fitting, which means the time is mainly consumed on the mean shift clustering method.Figure 14
Running time of the 3D line segment clustering algorithm and hypothetical plane fitting algorithm.


### 5.3. Triangle Mesh Construction and Subdivision

#### 5.3.1. Triangle Mesh Construction

In this test case, triangle mesh models are reconstructed using the algorithm described in Section 4.2 on the hypothetical planes model of the four datasets. Since the constrained Delaunay triangulation algorithm used in the paper is limited by the line segment in the 3D line segment model, the reconstructed triangle mesh preserves the geometric features of the building well. As shown in Figure15, the window, doors, and facades can be clearly identified. In addition, a “fake mesh region” phenomenon caused by noise line segments from natural landscape appeared in the results of the dataset coffee shack, which will not affect the final reconstruction results because it will be filtered by the subsequent subdivision penalty strategy.Figure 15
Results of the initial triangular meshes, and different hypothetical planes are presented with different colours. (a) House1; (b) house2; (c) hall; and (d) coffee shack.
(a)(b)(c)(d)


#### 5.3.2. Triangle Mesh Subdivision Algorithm

In this case, the algorithm described in Section 4.3 is used to perform multiround subdivision of the triangle meshes on the dataset house1 with different loops. The test result of a wall is shown in Figure16. The figure shows that the spatial subdivision algorithm can uniformly and densely subdivide the triangle meshes with the normal vectors of points transiting smoothly. There is no erroneous normal vector being contrary to that of the hypothetical plane, which proves the validity of updating the normal vector in the spatial subdivision algorithm.Figure 16
Results of different levels of triangular mesh subdivision of a wall of dataset house1. (a) The initial triangle meshes constructed; (b) normal vectors of vertexes on the initial triangle meshes; (c) the triangle meshes constructed when loop = 1; (d) normal vectors of the triangle meshes when loop = 1; (e) the triangle meshes constructed when loop = 2; and (f) normal vectors of the triangle meshes when loop = 2.
(a)(b)(c)(d)(e)(f)Figure17 shows the results of the subdivision test with a penalty mechanism. The subdivision algorithm with a penalty mechanism is able to treat those triangles discriminately, which effectively avoids the invalid subdivision of the meshes in the real cavity region, and improves the efficiency and accuracy. The subdivided meshes are able to faithfully present the original real plane: sparse in the cavity region and dense in the noncavity region.Figure 17
The subdivision results on two different hypothetical planes with a penalty mechanism (dataset hall), and invalid subdivision is effectively avoided. (a) Mesh of plane 1 after the 1st round; (b) mesh of plane 2 after the 1st round; (c) mesh of plane 1 after the 2nd round; (d) mesh of plane 2 after the 2nd round; (e) mesh of plane 1 after the 3rd round; and (f) mesh of plane 2 after the 3rd round.
(a)(b)(c)(d)(e)(f)


#### 5.3.3. Improved PMVS Patch Optimization Algorithm

In this test case, the improved PMVS patch optimization algorithm in Section 4.4 and the original PMVS patch optimization algorithm are used to optimize the same input point cloud produced by triangle mesh subdivision with different loop numbers separately. The comparison of the consumed time and the success rate are shown in Figures18 and 19.Figure 18
Time expense of two patching optimization algorithms: original PMVS patching optimization algorithm and our patching optimization algorithm.Figure 19
Time expense of two patching optimization algorithms: original PMVS patching optimization algorithm and our patching optimization algorithm.Figure18 shows that the time consumption of the improved patch optimization algorithm is significantly lower, and the running efficiency is increased by 65% on an average comparing with the original algorithm. Figure 20 shows the visual comparison between test results of the original and the improved patch optimization algorithms on point cloud from different viewports, and the processed point cloud is obtained after two rounds of triangle mesh subdivision. As shown in Figure 20, the improved algorithm has the same optimization ability as the original algorithm, and there is no clear visual difference on the two groups of optimized point clouds.Figure 20
Comparison experimental results between the original PMVS algorithm and our improved algorithm, both almost share the same visual output. (a) House1 model from the original PMVS, in side view; (b) house1 model from our algorithm, in side view; (c) house1 model from the original PMVS, in overlook view; and (d) house1 model from our algorithm, in overlook view.
(a)(b)(c)(d)From all above speaking, the comparison shows that the optimization effect of the improved optimization algorithm and the original PMVS optimization algorithm are visually similar. Although in some small area, the improved optimization algorithm could lead to partially missing points because of the failure of optimizing, and the consumed time of the improved optimization algorithm is only 35% of that of the original algorithm.


### 5.4. Integrated Experiment and Results Analysis

In the integrated experiment, algorithms including the hypothetical plane fitting algorithm based on the 3D line segment model, the mesh subdivision algorithm with the penalty mechanism, and the original PMVS expansion algorithm are used to reconstruct dense point clouds on the four datasets, respectively, and the setting parameters in the reconstructing process such as iteration number are shown in Table3.Table 3
Setting parameters of the test1.
House1House2HallCoffee shackn12222Loop times(2, 1)1(2, 1)(3, 2)(5, 2)n212221(2, 1) means loop1 = 2, loop2 = 1, and the same hereinafter.To compare our algorithm with the original PMVS algorithm, the original PMVS algorithm (actually PMVS2) is also run to reconstruct dense point clouds on the four datasets. The original PMVS algorithm has 3 iterations, and its parameters are all default settings.To visually compare the final results, the seed point cloud reconstructed by the original PMVS algorithm and the improved algorithm are shown in Figures21 and 22, respectively. From the two figures, we can find that the seed point cloud of the original PMVS algorithm is sparse and with a lot of noise, but the seed point cloud reconstructed by the improved algorithm is very clean and highly accurate.Figure 21
Seed point cloud reconstructed by the original PMVS algorithm, sparse, and with a lot of noise points. (a) House1; (b) house2; (c) hall; and (d) coffee shack.
(a)(b)(c)(d)Figure 22
Seed point cloud (semidense point cloud) reconstructed by the improved PMVS algorithm, which are much better than those of Figure19. (a) House1; (b) house2; (c) hall; and (d) coffee shack.
(a)(b)(c)(d)The reconstructed point clouds of the original PMVS algorithm and our improved algorithm are compared in Figure23. As shown in Figure 23, the reconstruction results of the improved algorithm are similar to that of the original PMVS algorithm or even no difference in the main part of the building, while they have some points missing in the nonmain part and nonbuilding area, which have been marked with a red circle in the figure.Figure 23
The comparison of reconstruction results of the original PMVS algorithm and our algorithm, in dense point cloud style. (a) House1 from the original PMVS; (b) house1 from our algorithm; (c) house2 from the original PMVS; (d) house2 from our algorithm; (e) hall from the original PMVS; (f) hall from our algorithm; (g) coffee shack from the original PVMS; and (h) coffee shack from our algorithm.
(a)(b)(c)(d)(e)(f)(g)(h)Considering the algorithm efficiency performance, the comparing test results of running time and point cloud density between the original PMVS algorithm and the improved algorithm are shown in Figures24 and 25, and it can be found that the improved algorithm is 15%∼23% faster than the original PMVS algorithm, while the point cloud losing rate is kept within −10.10%∼3.73%, which demonstrates that the improved algorithm efficiency is effectively improved. As known, the losing rate of point cloud refers to the reduction rate of the number of reconstruction points of the optimization algorithm comparing with the original algorithm. In this test, the improved algorithm shows the best performances on the dataset house1 with 10.10% more points and 23.32% less consumed time than the original PMVS algorithm, as shown in Figures 24 and 25, which proves that the improved algorithm is very suitable for buildings with simple geometric features and structures.Figure 24
Comparison of running time between the original PMVS algorithm and the improved PMVS algorithm, and our algorithm consumes less time on all test dataset.Figure 25
Comparison of the number of reconstructed point cloud between the original PMVS algorithm and the improved PMVS algorithm.In addition, as shown in Figure23, the point cloud reconstructed by the original PMVS algorithm contains many noise points (shown in the red oval box in Figure 23(a)), which proves that the improved algorithm not only improves the efficiency but also improves the accuracy of the reconstructed model by removing a large amount of noise points.It could be found that the final reconstructed point clouds lose partial points, and the root cause is that these three buildings have some special structures which are not included in the hypothetical plane model of buildings. And these structures are more independent of the building’s main body. When there are few seed points appearing on these structures, it is difficult to realize the reconstruction of these structures in the expansion iteration process of PMVS. For example, the two chimneys and column towers on the building’s roof in the dataset house2 and the white wall on the left side of the main cabin in the dataset coffee shack were not reconstructed by the algorithm. However, these parts are not the building’s main body, and the original PWMS algorithm also failed to reconstruct them and failed to build a consolidated building model.Overall speaking, from the aforementioned experimental results, we can find that the improved PMVS algorithm achieves similar point cloud reconstruction accuracy with the original PMVS algorithm, but our algorithm is 15% to ∼23% faster than the original PMVS algorithm in running time, and our algorithm successfully removed many noise points in the buildings.


## 5.1. Test Datasets and Test Environment

To verify the performance of the algorithm presented in this paper, four sets of building sequence photos are used, which are named house1, house2, hall, and coffee shack, respectively. The detailed information such as the photo numbers and resolutions are shown in Table1, and the previews of the four datasets are shown in Figure 8.Table 1
Numbers of photo and resolution in experimental dataset1.
Dataset namePhoto numberResolutionHouse1652424 × 1351House2263071 × 2301Hall611534 × 1008Coffee shack301086 × 7151House1 is from https://github.com/zxg519/3D_reconstruction_dataset/, house2 is from https://demuc.de/colmap/datasets/, hall is from http://mirror.openstreetmap.nl/sources/pmvs-2/data/hall/, and coffee shack is from http://idav.ucdavis.edu/∼hkim/mvs/dataset/.Figure 8
Sequence image previews of the experimental datasets.The hardware configuration is as follows: CPU: Intel Core i5-6500 @ 3.2 GHz quad core; GPU: NVIDIA GeForce GTX 750, 2 GB RAM; and memory: 16G.The experiment carried out below can be divided into three parts. The first part is designed to verify the effectiveness of the hypothetical plane fitting algorithm, and the second part is designed to verify the effectiveness of triangle mesh reconstruction and subdivision. The third part runs the whole algorithm flow and compares the time consumption and the quality of point cloud reconstruction with the original PMVS algorithm.


## 5.2. Hypothetical Plane Fitting

### 5.2.1. 3D Line Segment Clustering Algorithm Based on Manhattan World Assumption

Before the experiment, the COLMAP’s SfM algorithm [25] is used to perform camera calibration and sparse point cloud reconstruction for each experimental dataset; thus, the camera intrinsic and extrinsic parameters, initial feature point matching information, and sparse point cloud information are obtained. After that, Line3D++ [16], a 3D line segment model reconstruction algorithm, is used to obtain all 3D line segment information for each dataset.These 3D line segment models are used as the input data, and clustering operations are performed using a 3D line segment clustering algorithm based on Manhattan world assumption. The 3D line segment clustering results of the four datasets are shown in Figure9 (different types of line segments are represented using different colours in the figure), and the figure shows that the algorithm successfully extracts the main directions of the building line segment model, which includes the main axis-aligned and nonaxis-aligned directions. It could be found that(1)
It performs effectively on the nonaxis-aligned direction and works robustly on clustering of short line segments. As shown in Figure10, the nonaxis line segments in the datasets house1 and hall are still clustered correctly, although they are very short (about 1 m in Figure 10(a) and 0.2 m in Figure 10(b)), and the roof on the left side of the dataset house2 is not extracted just because the detected segments are too sparse.(2)
It works robustly when there is strong noise caused by natural landscapes. Most of the natural landscapes do not follow the Manhattan world hypothesis, and its reconstructed line segments can easily be filtered out because it is difficult to find their category in line segment clustering. Therefore, besides line segment clustering functionality, it also has the filtering function of line segment noise brought by natural landscape or nonartificial objects. In the dataset coffee shack, the trees introduce a large number of noise line segments to the Line3D++ reconstruction result, as shown in Figure11(a), while after using the line segment clustering algorithm proposed in this paper, most of these noise segments are filtered out, as shown in Figure 11(b).Figure 9
Experimental results of the 3D line segment clustering algorithm based on Manhattan world assumption, and different colours used to represent different line segment clusters. (a) Results of dataset house1; (b) results of dataset house2; (c) results of dataset hall; and (d) results of dataset coffee shack.
(a)(b)(c)(d)Figure 10
Results of the 3D line segment clustering algorithm in nonaxis-aligned directions, and different clusters are shown using different colours. (a) Results of house1 and (b) results of hall.
(a)(b)Figure 11
The comparison results of the 3D line segment model of coffee shack before and after line segment clustering show that the algorithm can filter line segment noise in a certain degree on the line segment noise generated by natural landscape. (a) Before clustering; and (b) after clustering.
(a)(b)


### 5.2.2. Multihypothesis Plane Fitting Based on the Mean Shift Method

After obtaining the 3D line segment clusters of each dataset, the multihypothesis plane fitting algorithm based on the mean shift method is used to fit the point cloud sampled from the line segment model. The results of plane fitting are shown in Figures12 and 13. In the two figures, the hypothetical plane models of each dataset are presented in a point cloud clustering manner. Each point cloud cluster represents a hypothetical plane using different colours.Figure 12
The fitting results of hypothetical planes of different buildings where different hypothetical planes are represented using different colours (I). (a) House1 in the top view; (b) house1 in the oblique view; (c) house2 in the top view; and (d) house2 in the top view.
(a)(b)(c)(d)Figure 13
The fitting results of hypothetical planes of different buildings where different hypothetical planes are represented using different colours (II). (a) Hall in the front view; (b) hall in the top view; (c) coffee shack in the oblique view; and (d) coffee shack in the side view.
(a)(b)(c)(d)The fitting results show that the main planes of each building in the four datasets (including the main facades and the planes where the roof is located) are correctly extracted.The extraction completeness information of hypothetical planes is shown in Table2, where the completeness rate refers to the ratio of the number of final extracted hypothetical planes to the number of actual hypothetical planes. The actual hypothetical plane refers to the plane determined manually from the 3D line segment model of the dataset (tiny planes are not counted). Table 2 shows that our algorithm maintains a high completeness rate of hypothetical plane extraction for each dataset. Specially, the root causes of the five cases less than 1.00 are as follows: (1) there are too few reconstructed 3D line segments to successfully recover the plane information in the missed hypothetical plane; and (2) the quantitative area of the missed plane is too small or the line segments is too short.Table 2
Extraction completeness of hypothetical planes.
House1House2HallCoffee shackExtracted (X-axis)4263Actual (X-axis)4264Extracted (Y-axis)21172Actual (Y-axis)22202Extracted (Z-axis)0001Actual (Z-axis)0021Extracted (nonaxis)2422Actual (nonaxis)2423Although the algorithm may miss some hypothetical planes with small areas, the main hypothetical planes are correctly extracted, which will show little effect on the subsequent triangle mesh subdivision.The consumed time information of the test is shown in Figure14. The relationship between the consumed time and the total number of 3D line segments shows a linear-like relationship, and the time on line segment clustering is much less than that of plane fitting, which means the time is mainly consumed on the mean shift clustering method.Figure 14
Running time of the 3D line segment clustering algorithm and hypothetical plane fitting algorithm.


## 5.2.1. 3D Line Segment Clustering Algorithm Based on Manhattan World Assumption

Before the experiment, the COLMAP’s SfM algorithm [25] is used to perform camera calibration and sparse point cloud reconstruction for each experimental dataset; thus, the camera intrinsic and extrinsic parameters, initial feature point matching information, and sparse point cloud information are obtained. After that, Line3D++ [16], a 3D line segment model reconstruction algorithm, is used to obtain all 3D line segment information for each dataset.These 3D line segment models are used as the input data, and clustering operations are performed using a 3D line segment clustering algorithm based on Manhattan world assumption. The 3D line segment clustering results of the four datasets are shown in Figure9 (different types of line segments are represented using different colours in the figure), and the figure shows that the algorithm successfully extracts the main directions of the building line segment model, which includes the main axis-aligned and nonaxis-aligned directions. It could be found that(1)
It performs effectively on the nonaxis-aligned direction and works robustly on clustering of short line segments. As shown in Figure10, the nonaxis line segments in the datasets house1 and hall are still clustered correctly, although they are very short (about 1 m in Figure 10(a) and 0.2 m in Figure 10(b)), and the roof on the left side of the dataset house2 is not extracted just because the detected segments are too sparse.(2)
It works robustly when there is strong noise caused by natural landscapes. Most of the natural landscapes do not follow the Manhattan world hypothesis, and its reconstructed line segments can easily be filtered out because it is difficult to find their category in line segment clustering. Therefore, besides line segment clustering functionality, it also has the filtering function of line segment noise brought by natural landscape or nonartificial objects. In the dataset coffee shack, the trees introduce a large number of noise line segments to the Line3D++ reconstruction result, as shown in Figure11(a), while after using the line segment clustering algorithm proposed in this paper, most of these noise segments are filtered out, as shown in Figure 11(b).Figure 9
Experimental results of the 3D line segment clustering algorithm based on Manhattan world assumption, and different colours used to represent different line segment clusters. (a) Results of dataset house1; (b) results of dataset house2; (c) results of dataset hall; and (d) results of dataset coffee shack.
(a)(b)(c)(d)Figure 10
Results of the 3D line segment clustering algorithm in nonaxis-aligned directions, and different clusters are shown using different colours. (a) Results of house1 and (b) results of hall.
(a)(b)Figure 11
The comparison results of the 3D line segment model of coffee shack before and after line segment clustering show that the algorithm can filter line segment noise in a certain degree on the line segment noise generated by natural landscape. (a) Before clustering; and (b) after clustering.
(a)(b)


## 5.2.2. Multihypothesis Plane Fitting Based on the Mean Shift Method

After obtaining the 3D line segment clusters of each dataset, the multihypothesis plane fitting algorithm based on the mean shift method is used to fit the point cloud sampled from the line segment model. The results of plane fitting are shown in Figures12 and 13. In the two figures, the hypothetical plane models of each dataset are presented in a point cloud clustering manner. Each point cloud cluster represents a hypothetical plane using different colours.Figure 12
The fitting results of hypothetical planes of different buildings where different hypothetical planes are represented using different colours (I). (a) House1 in the top view; (b) house1 in the oblique view; (c) house2 in the top view; and (d) house2 in the top view.
(a)(b)(c)(d)Figure 13
The fitting results of hypothetical planes of different buildings where different hypothetical planes are represented using different colours (II). (a) Hall in the front view; (b) hall in the top view; (c) coffee shack in the oblique view; and (d) coffee shack in the side view.
(a)(b)(c)(d)The fitting results show that the main planes of each building in the four datasets (including the main facades and the planes where the roof is located) are correctly extracted.The extraction completeness information of hypothetical planes is shown in Table2, where the completeness rate refers to the ratio of the number of final extracted hypothetical planes to the number of actual hypothetical planes. The actual hypothetical plane refers to the plane determined manually from the 3D line segment model of the dataset (tiny planes are not counted). Table 2 shows that our algorithm maintains a high completeness rate of hypothetical plane extraction for each dataset. Specially, the root causes of the five cases less than 1.00 are as follows: (1) there are too few reconstructed 3D line segments to successfully recover the plane information in the missed hypothetical plane; and (2) the quantitative area of the missed plane is too small or the line segments is too short.Table 2
Extraction completeness of hypothetical planes.
House1House2HallCoffee shackExtracted (X-axis)4263Actual (X-axis)4264Extracted (Y-axis)21172Actual (Y-axis)22202Extracted (Z-axis)0001Actual (Z-axis)0021Extracted (nonaxis)2422Actual (nonaxis)2423Although the algorithm may miss some hypothetical planes with small areas, the main hypothetical planes are correctly extracted, which will show little effect on the subsequent triangle mesh subdivision.The consumed time information of the test is shown in Figure14. The relationship between the consumed time and the total number of 3D line segments shows a linear-like relationship, and the time on line segment clustering is much less than that of plane fitting, which means the time is mainly consumed on the mean shift clustering method.Figure 14
Running time of the 3D line segment clustering algorithm and hypothetical plane fitting algorithm.


## 5.3. Triangle Mesh Construction and Subdivision

### 5.3.1. Triangle Mesh Construction

In this test case, triangle mesh models are reconstructed using the algorithm described in Section 4.2 on the hypothetical planes model of the four datasets. Since the constrained Delaunay triangulation algorithm used in the paper is limited by the line segment in the 3D line segment model, the reconstructed triangle mesh preserves the geometric features of the building well. As shown in Figure15, the window, doors, and facades can be clearly identified. In addition, a “fake mesh region” phenomenon caused by noise line segments from natural landscape appeared in the results of the dataset coffee shack, which will not affect the final reconstruction results because it will be filtered by the subsequent subdivision penalty strategy.Figure 15
Results of the initial triangular meshes, and different hypothetical planes are presented with different colours. (a) House1; (b) house2; (c) hall; and (d) coffee shack.
(a)(b)(c)(d)


### 5.3.2. Triangle Mesh Subdivision Algorithm

In this case, the algorithm described in Section 4.3 is used to perform multiround subdivision of the triangle meshes on the dataset house1 with different loops. The test result of a wall is shown in Figure16. The figure shows that the spatial subdivision algorithm can uniformly and densely subdivide the triangle meshes with the normal vectors of points transiting smoothly. There is no erroneous normal vector being contrary to that of the hypothetical plane, which proves the validity of updating the normal vector in the spatial subdivision algorithm.Figure 16
Results of different levels of triangular mesh subdivision of a wall of dataset house1. (a) The initial triangle meshes constructed; (b) normal vectors of vertexes on the initial triangle meshes; (c) the triangle meshes constructed when loop = 1; (d) normal vectors of the triangle meshes when loop = 1; (e) the triangle meshes constructed when loop = 2; and (f) normal vectors of the triangle meshes when loop = 2.
(a)(b)(c)(d)(e)(f)Figure17 shows the results of the subdivision test with a penalty mechanism. The subdivision algorithm with a penalty mechanism is able to treat those triangles discriminately, which effectively avoids the invalid subdivision of the meshes in the real cavity region, and improves the efficiency and accuracy. The subdivided meshes are able to faithfully present the original real plane: sparse in the cavity region and dense in the noncavity region.Figure 17
The subdivision results on two different hypothetical planes with a penalty mechanism (dataset hall), and invalid subdivision is effectively avoided. (a) Mesh of plane 1 after the 1st round; (b) mesh of plane 2 after the 1st round; (c) mesh of plane 1 after the 2nd round; (d) mesh of plane 2 after the 2nd round; (e) mesh of plane 1 after the 3rd round; and (f) mesh of plane 2 after the 3rd round.
(a)(b)(c)(d)(e)(f)


### 5.3.3. Improved PMVS Patch Optimization Algorithm

In this test case, the improved PMVS patch optimization algorithm in Section 4.4 and the original PMVS patch optimization algorithm are used to optimize the same input point cloud produced by triangle mesh subdivision with different loop numbers separately. The comparison of the consumed time and the success rate are shown in Figures18 and 19.Figure 18
Time expense of two patching optimization algorithms: original PMVS patching optimization algorithm and our patching optimization algorithm.Figure 19
Time expense of two patching optimization algorithms: original PMVS patching optimization algorithm and our patching optimization algorithm.Figure18 shows that the time consumption of the improved patch optimization algorithm is significantly lower, and the running efficiency is increased by 65% on an average comparing with the original algorithm. Figure 20 shows the visual comparison between test results of the original and the improved patch optimization algorithms on point cloud from different viewports, and the processed point cloud is obtained after two rounds of triangle mesh subdivision. As shown in Figure 20, the improved algorithm has the same optimization ability as the original algorithm, and there is no clear visual difference on the two groups of optimized point clouds.Figure 20
Comparison experimental results between the original PMVS algorithm and our improved algorithm, both almost share the same visual output. (a) House1 model from the original PMVS, in side view; (b) house1 model from our algorithm, in side view; (c) house1 model from the original PMVS, in overlook view; and (d) house1 model from our algorithm, in overlook view.
(a)(b)(c)(d)From all above speaking, the comparison shows that the optimization effect of the improved optimization algorithm and the original PMVS optimization algorithm are visually similar. Although in some small area, the improved optimization algorithm could lead to partially missing points because of the failure of optimizing, and the consumed time of the improved optimization algorithm is only 35% of that of the original algorithm.


## 5.3.1. Triangle Mesh Construction

In this test case, triangle mesh models are reconstructed using the algorithm described in Section 4.2 on the hypothetical planes model of the four datasets. Since the constrained Delaunay triangulation algorithm used in the paper is limited by the line segment in the 3D line segment model, the reconstructed triangle mesh preserves the geometric features of the building well. As shown in Figure15, the window, doors, and facades can be clearly identified. In addition, a “fake mesh region” phenomenon caused by noise line segments from natural landscape appeared in the results of the dataset coffee shack, which will not affect the final reconstruction results because it will be filtered by the subsequent subdivision penalty strategy.Figure 15
Results of the initial triangular meshes, and different hypothetical planes are presented with different colours. (a) House1; (b) house2; (c) hall; and (d) coffee shack.
(a)(b)(c)(d)


## 5.3.2. Triangle Mesh Subdivision Algorithm

In this case, the algorithm described in Section 4.3 is used to perform multiround subdivision of the triangle meshes on the dataset house1 with different loops. The test result of a wall is shown in Figure16. The figure shows that the spatial subdivision algorithm can uniformly and densely subdivide the triangle meshes with the normal vectors of points transiting smoothly. There is no erroneous normal vector being contrary to that of the hypothetical plane, which proves the validity of updating the normal vector in the spatial subdivision algorithm.Figure 16
Results of different levels of triangular mesh subdivision of a wall of dataset house1. (a) The initial triangle meshes constructed; (b) normal vectors of vertexes on the initial triangle meshes; (c) the triangle meshes constructed when loop = 1; (d) normal vectors of the triangle meshes when loop = 1; (e) the triangle meshes constructed when loop = 2; and (f) normal vectors of the triangle meshes when loop = 2.
(a)(b)(c)(d)(e)(f)Figure17 shows the results of the subdivision test with a penalty mechanism. The subdivision algorithm with a penalty mechanism is able to treat those triangles discriminately, which effectively avoids the invalid subdivision of the meshes in the real cavity region, and improves the efficiency and accuracy. The subdivided meshes are able to faithfully present the original real plane: sparse in the cavity region and dense in the noncavity region.Figure 17
The subdivision results on two different hypothetical planes with a penalty mechanism (dataset hall), and invalid subdivision is effectively avoided. (a) Mesh of plane 1 after the 1st round; (b) mesh of plane 2 after the 1st round; (c) mesh of plane 1 after the 2nd round; (d) mesh of plane 2 after the 2nd round; (e) mesh of plane 1 after the 3rd round; and (f) mesh of plane 2 after the 3rd round.
(a)(b)(c)(d)(e)(f)


## 5.3.3. Improved PMVS Patch Optimization Algorithm

In this test case, the improved PMVS patch optimization algorithm in Section 4.4 and the original PMVS patch optimization algorithm are used to optimize the same input point cloud produced by triangle mesh subdivision with different loop numbers separately. The comparison of the consumed time and the success rate are shown in Figures18 and 19.Figure 18
Time expense of two patching optimization algorithms: original PMVS patching optimization algorithm and our patching optimization algorithm.Figure 19
Time expense of two patching optimization algorithms: original PMVS patching optimization algorithm and our patching optimization algorithm.Figure18 shows that the time consumption of the improved patch optimization algorithm is significantly lower, and the running efficiency is increased by 65% on an average comparing with the original algorithm. Figure 20 shows the visual comparison between test results of the original and the improved patch optimization algorithms on point cloud from different viewports, and the processed point cloud is obtained after two rounds of triangle mesh subdivision. As shown in Figure 20, the improved algorithm has the same optimization ability as the original algorithm, and there is no clear visual difference on the two groups of optimized point clouds.Figure 20
Comparison experimental results between the original PMVS algorithm and our improved algorithm, both almost share the same visual output. (a) House1 model from the original PMVS, in side view; (b) house1 model from our algorithm, in side view; (c) house1 model from the original PMVS, in overlook view; and (d) house1 model from our algorithm, in overlook view.
(a)(b)(c)(d)From all above speaking, the comparison shows that the optimization effect of the improved optimization algorithm and the original PMVS optimization algorithm are visually similar. Although in some small area, the improved optimization algorithm could lead to partially missing points because of the failure of optimizing, and the consumed time of the improved optimization algorithm is only 35% of that of the original algorithm.


## 5.4. Integrated Experiment and Results Analysis

In the integrated experiment, algorithms including the hypothetical plane fitting algorithm based on the 3D line segment model, the mesh subdivision algorithm with the penalty mechanism, and the original PMVS expansion algorithm are used to reconstruct dense point clouds on the four datasets, respectively, and the setting parameters in the reconstructing process such as iteration number are shown in Table3.Table 3
Setting parameters of the test1.
House1House2HallCoffee shackn12222Loop times(2, 1)1(2, 1)(3, 2)(5, 2)n212221(2, 1) means loop1 = 2, loop2 = 1, and the same hereinafter.To compare our algorithm with the original PMVS algorithm, the original PMVS algorithm (actually PMVS2) is also run to reconstruct dense point clouds on the four datasets. The original PMVS algorithm has 3 iterations, and its parameters are all default settings.To visually compare the final results, the seed point cloud reconstructed by the original PMVS algorithm and the improved algorithm are shown in Figures21 and 22, respectively. From the two figures, we can find that the seed point cloud of the original PMVS algorithm is sparse and with a lot of noise, but the seed point cloud reconstructed by the improved algorithm is very clean and highly accurate.Figure 21
Seed point cloud reconstructed by the original PMVS algorithm, sparse, and with a lot of noise points. (a) House1; (b) house2; (c) hall; and (d) coffee shack.
(a)(b)(c)(d)Figure 22
Seed point cloud (semidense point cloud) reconstructed by the improved PMVS algorithm, which are much better than those of Figure19. (a) House1; (b) house2; (c) hall; and (d) coffee shack.
(a)(b)(c)(d)The reconstructed point clouds of the original PMVS algorithm and our improved algorithm are compared in Figure23. As shown in Figure 23, the reconstruction results of the improved algorithm are similar to that of the original PMVS algorithm or even no difference in the main part of the building, while they have some points missing in the nonmain part and nonbuilding area, which have been marked with a red circle in the figure.Figure 23
The comparison of reconstruction results of the original PMVS algorithm and our algorithm, in dense point cloud style. (a) House1 from the original PMVS; (b) house1 from our algorithm; (c) house2 from the original PMVS; (d) house2 from our algorithm; (e) hall from the original PMVS; (f) hall from our algorithm; (g) coffee shack from the original PVMS; and (h) coffee shack from our algorithm.
(a)(b)(c)(d)(e)(f)(g)(h)Considering the algorithm efficiency performance, the comparing test results of running time and point cloud density between the original PMVS algorithm and the improved algorithm are shown in Figures24 and 25, and it can be found that the improved algorithm is 15%∼23% faster than the original PMVS algorithm, while the point cloud losing rate is kept within −10.10%∼3.73%, which demonstrates that the improved algorithm efficiency is effectively improved. As known, the losing rate of point cloud refers to the reduction rate of the number of reconstruction points of the optimization algorithm comparing with the original algorithm. In this test, the improved algorithm shows the best performances on the dataset house1 with 10.10% more points and 23.32% less consumed time than the original PMVS algorithm, as shown in Figures 24 and 25, which proves that the improved algorithm is very suitable for buildings with simple geometric features and structures.Figure 24
Comparison of running time between the original PMVS algorithm and the improved PMVS algorithm, and our algorithm consumes less time on all test dataset.Figure 25
Comparison of the number of reconstructed point cloud between the original PMVS algorithm and the improved PMVS algorithm.In addition, as shown in Figure23, the point cloud reconstructed by the original PMVS algorithm contains many noise points (shown in the red oval box in Figure 23(a)), which proves that the improved algorithm not only improves the efficiency but also improves the accuracy of the reconstructed model by removing a large amount of noise points.It could be found that the final reconstructed point clouds lose partial points, and the root cause is that these three buildings have some special structures which are not included in the hypothetical plane model of buildings. And these structures are more independent of the building’s main body. When there are few seed points appearing on these structures, it is difficult to realize the reconstruction of these structures in the expansion iteration process of PMVS. For example, the two chimneys and column towers on the building’s roof in the dataset house2 and the white wall on the left side of the main cabin in the dataset coffee shack were not reconstructed by the algorithm. However, these parts are not the building’s main body, and the original PWMS algorithm also failed to reconstruct them and failed to build a consolidated building model.Overall speaking, from the aforementioned experimental results, we can find that the improved PMVS algorithm achieves similar point cloud reconstruction accuracy with the original PMVS algorithm, but our algorithm is 15% to ∼23% faster than the original PMVS algorithm in running time, and our algorithm successfully removed many noise points in the buildings.


## 6. Conclusions

For accelerating building reconstruction, this paper presents a novel method, which generates seed points for the PMVS expansion pipeline using the 3D line segment information. Two algorithms are proposed in the study: a hypothetical plane fitting algorithm using the 3D line segment information of the building and a triangle mesh subdivision and patch optimization algorithm with a penalty mechanism, and the first algorithm works as the input of the second one:(1)
Based on Manhattan world assumption, the hypothetical planes fitting algorithm clusters the 3D line segments by using greedy thought to extract the coordinate directions and which later is used to guide the subsequent procedures including the clustering of point cloud based on the mean shift method and the fitting of hypothetical planes(2)
The triangle mesh subdivision and patch optimization algorithm generates the semidense point cloud using hypothetical plane information with a penalty mechanism, and then, the semidense point cloud is sent to the original PMVS expansion pipeline as the seed points(3)
Experiments show that the improved approach not only accelerates the speed of building reconstruction but also provides visually better reconstruction results and more accurate performance through effectively removing noise points from the building modelThe future work will focus on how to effectively recognize those small hypothetical planes in buildings and recover those small features without being affected by noises, and technologies such as deep learning would be used to recognize the building structures and guide the 3D reconstruction process.



---

*Source: 9267854-2020-09-10.xml*