AN ORIGINAL ALGORITHM FOR BIM GENERATION FROM INDOOR SURVEY POINT CLOUDS

: Nowadays, it is essential to ﬁnd new strategies, able to perform the ﬁrst step of the scan-to-BIM process, by retrieving the geometrical information contained in point clouds that are so easily collected through laser scanners and range cameras. This paper presents a new algorithm for the automatic extraction of the layout and the height of a small indoor environment from its point cloud. In particular, the algorithm was tested on a point cloud of 600000 vertices, selected from the dataset of the ISPRS benchmark on indoor modelling. The preliminary results are encouraging: the 3D shape (layout and height) of the investigated room is effectively reconstructed.


INTRODUCTION
Building Information Modelling (BIM) offers several benefits in the design, planning, and construction of new buildings as well as in the management and maintenance of the existing ones.Frequently, however, BIMs are not available for older constructions, and, in some cases, even blueprints are missing.Consequently, it is essential to find new strategies to efficiently generate BIMs for the already existing and used buildings (Xiong et al., 2013).At the same time, the recent advances in surveying technology (e.g.laser scanners and range cameras) make relatively easy to collect dense point clouds of indoor environments, even if the current processing solutions are not yet ripe to analyse the massive datasets produced.Therefore, it is increasingly necessary to develop automatic and efficient algorithms able to rapidly extract at least the geometric BIM features from indoor point clouds, such as layouts and room heights.
In this paper, the algorithm presented in (Capocchiano et al., 2017) is extended and enhanced.Originally designed to extract layouts from the 3D models of rooms acquired by means of a low-cost range camera, the algorithm is now able to process also very dense point clouds provided by professional laser scanners.Specifically, the implementation was refined in order to: • increase the computational efficiency of the algorithm; • retrieve also the height of the scanned room (supposed constant) along with the layout.
In particular, the new version was implemetned in Python and tested on a point cloud of 600000 vertices (Figure 1), selected from the dataset of the ISPRS benchmark on indoor modelling (Khoshelham et al., 2017).
This paper is organized as follows: Section 2 describes the workflow of the algorithm in detail and Section 3 presents and discusses the obtained results.Finally, in Section 4, some conclusions are drawn and future prospects are outlined.* Corresponding author.

THE ALGORITHM
The algorithm retrieves the layout and the heigth of an indoor environment (e.g. a room) starting from its 3D point cloud, which have to include the ceiling, the upper portion of the walls and a small portion of the floor (Figure 1).If the mesh is also available, as is common when working with mobile devices equipped with range cameras, it could be useful to carry out a preliminary statistical analysis of the length of the polygon edges (for the selection of the triangulation parameters, as it will be described later).However, from this point on, only the information contained in the vertices of the point cloud are considered by the algorithm.
In particular, the user has to set the values of different parameters belonging to the following three categories: • parameters depending on the available computing power (e.g.maximum number of executions for a loop); • parameters depending on the features of the specific sensor whereby the point cloud was collected; • parameters strictly related to the specific survey, such as the number of the sides of the room.
The reference frame associated with the point cloud is assumed arbitrary (even though professional scanners include sensors capable to detect the vertical direction by means of gravity), thus the equation of the plane modeling the ceiling is identified through the RANdom SAmple Consensus (RANSAC) algorithm (Fischler and Bolles, 1981).RANSAC is iterated multiple times: each iteration takes as input the inlier set of the previous one (Ravanelli et al., 2015) and works with a threshold reduced by 1 cm with respect to the previous iteration.The parameters required are: • the starting threshold, depending on the sensor accuracy and the real ceiling flatness; • the stop criteria.In point clouds collected by low-cost range cameras, the ceiling may present a variable slope: an excessively small threshold would force RANSAC to restrict the inlier set to a portion of the ceiling with uniform slope, i.e. invalidating the plane estimation.To prevent this issue, conditions on the minimum threshold, the maximum angle (the angle obtained by inverting the scalar product of the two normal unit vectors to the planes in question and is comprised between 0 • and 180 • ) between two consecutive estimated planes and the maximum reduction factor (the relationship between the number of inliers of the iteration N, and the number of inliers of the iteration N -1 must be greater than a certain value, typically between 0.8 and 0.9) for the number of inliers, should be considered.
Eventually, the final RANSAC result is refined by estimating the plane parameters only on the inlier set through the Ordinary Least Squares method.After that, a new reference frame with the Z' axis orthogonal to the ceiling plane (pointing outward) is defined (Figure 2): the new coordinates of the vertices are subsequently computed by means of quaternions (Sansò, 1973).This reference frame transformation is useful because in this way the room height can be easily computed as the absolute value of the minimum of Z'(Figure 3).Subsequently, the portion of the point cloud exceeding a specific distance from the plane is excluded from further processing (Figure 3).This distance depends on the features of the environment surveyed, e.g. the presence of windows, false ceilings, furniture, air conditioning and in general every object that can obstruct the scanning process and hide parts of the ceiling-walls border.The remaining points are thus orthogonally projected on the ceiling plane, reducing the original 3-dimensional problem into a 2-dimensional one (Figure 4).
Then, considering that only the boundary of the point cloud is interesting for the perimeter detection, namely for the layout extrac-tion, a cloud reduction is performed to keep an acceptable execution time in the following processing steps.This simplification is virtually mandatory when very dense point clouds collected by professional laser scanners are used and can be summarized in the following steps: 1. the bounding box of the 2D point cloud (from now on P0) is subdivided into a grid of squares s'ij with side ls; 2. the squares containing at least one point are included in the reduced point cloud P1if they share at least one vertex with an empty square or if they are placed at the border of the bounding box; at this stage the final reduced point cloud Pr is equal to P1; 3. the squares containing at least one point previously excluded from Pr and sharing at least one vertex with a square which belongs to Pr form a "monitored area" M1; 4. the bounding box of the full 2D point cloud is subdivided again into squares (s"ij) with side length ls/k1 (k1 ∈ N ).
A second reduced point cloud P2 is assembled under the above-mentioned conditions; 5. if s" ab ⊂ P2, s" ab ⊂ M1 and s" ab ⊂ P1 then the s" ab is added to the final reduced point cloud Pr.If s" ab ⊂ P1 but s" ab ⊂ P2 then s" ab is removed from Pr. Whenever a square s' cd belonging to M1 has a point added to the final reduced point cloud Pr, every square s' ab that share a vertex with such square is added to M1, and the control operated at the beginning of this point is reiterated, until no more points are added to Pr; 6. it is possible to define a new monitored area M2 (same criteria adopted for M1, this time with s"ij squares) together with a new simplified point cloud P3 comprised of squares s"'ij with side length ls/k1k2 (k1, k2 ∈ N ), in order to refine again P1 following the procedure set up in the previous step.
In general, the procedure described in step 6 can ideally take place many times with an arbitrary number of reduced point clouds, but limiting the process to P2, with a proper choice of the parameter k1, allows to obtain satisfactory results.Moreover, it is highly recommended to avoid removing further points from the simplified cloud of the previous step (as explained in step 5) because of the excessive narrowing of the border area that would occur in this case.Regarding the choice of ls, the task is complicated by conflicting needs.A small value of ls assures that all the border points required for the layout extraction are included in the simplified point cloud (Figure 5).In this way, though, the border of a hole placed in the central portion of the point cloud can be erroneously added to the reduced point cloud (Figure 5b).On the other hand, a large value of ls prevents the inclusion of central holes, but at the cost of a lower resolution in the border zone: significant points may be excluded if located near reflex interior angles (in room with concave geometry) and a large number of non-significant points can be included in the s' squares in P1 together with border ones (Figure 6).
The aim of the monitored area is precisely to find a compromise between these two equally inconvenient configurations: a value of ls comparable to the width of the holes ensures the removal of as many central points as possible, whereas the small value of ls/k1 (or even more ls/k1k2) enhances the resolution of Pr only where it is necessary, operating inside P1 and M1 (Figure  boundary can only be consistently defined in relation to a subset of the R2 plane rather than the current set of isolated points.Therefore, a Delaunay Triangulation is applied to the reduced point cloud (Figure 8a), imposing a constraint on the maximum edge length, i.e. the Maximum Edge Limit (MEL).A high value for the MEL can cause an alteration in the reconstruction of the final layout because, particularly in concave room, triangular meshes that do not match any real portion of the room are included in the surface from which the boundary will be computed later (Figure 8b).An excessively small MEL produces instead the opposite problem: only isolated portions of the point cloud characterized by small distances between the points are triangulated, and the final layout will not represent the real one at all (Figure 8b).
The correct choice of the MEL parameter depends on the resolution of the sensor employed for the survey: professional laser scanners provide very dense point cloud that can easily be triangulated in full with a MEL of 2 or 3 cm.Instead, the density of a point cloud provided by a low-cost range camera is sensibly lower, so that MELs up to 15 cm may be required.Furthermore, in case of a point cloud originally provided with meshes, as afore mentioned, a statistical analysis can provide a useful advice for the MEL choice: for instance, setting the MEL equal to the maximum edge length of the original triangulation assures a full Once the final meshes are identified, considering that in a Delaunay Triangulation a segment can only belong to one or two different triangular meshes, it is possible to easily detect the segments of the point cloud boundary because these ones satisfy the first case.These boundary segments form several closed simple polygonal chains (Figure 9), one of which contains all the others and approximates the cloud external border.
However, smaller polygonal chains surround the holes of the point cloud (Figure 10), either the original ones or the biggest one corresponding to the central portion of the cloud removed during the cloud reduction process.The separated polygonal chains are defined as a sequence of boundary segments that share one endpoint.The proximity between an eventual hole and the external region, or between two or more holes, can cause the connection of more than two boundary segments through a point p (Figure 11a).This quite frequent situation may lead to ambiguity or make the external polygonal chain merge with another chain.To prevent this issue, all the vertices pi that are end-point of more than two boundary segments are checked through the following steps: 1. all the segments connected to pi, including the boundary segments (those that we want to identify) and the internal ones, are selected (right panel of Figure 11a);  2. the segments are sorted according to the angle formed with the x axis, counted counterclockwise in the [0, 2π) interval; 3. two boundary segments connected to pi will be in the same polygonal chain only if they are consecutive, or first and last, in the sorted list of the previous step, and if they are not part of the same triangular mesh (Figure 11b); After this check, it is possible to isolate each polygonal chain from the others, and to select the only one meaningful for the last section of the algorithm: the one that has the same bounding box as the original 2D point cloud (Figure 12).
The last step of the algorithm is aimed to approximate the found boundary polygonal chain with a polygon of n sides, where n is equal to the number of sides of the room surveyed.This task requires once again the use of RANSAC, this time employing a linear model.Firstly, RANSAC is recursively applied for n times on the vertices of the polygonal chain, with the outlier set of an iteration serving as the input set for the next one.In each iteration, the RANSAC estimated straight line corresponds to one side of the real room, often the longest available at that point (Figure 13).
If there are two or more sides of the room that are representable as different segments of the same straight line, a check on the inlier set of each RANSAC estimation is needed to avoid that one iteration removes the points corresponding to more than one side (Figure 14): • the maximum distance between a pair of points belonging to the inlier set is computed; one point of this pair is thus selected as a starting point; • the points of the inlier set are sorted according to the distance from the starting point; • the distance from every pair of consecutive points in the sorted list is calculated; • distances exceeding a reasonable maximum distance parameter (for instance 2×MEL) mark a separation within the inlier set: points before and after that distance belong to two (or more) different real side of the room; • one group is selected as inlier set for the current RANSAC iteration, the other group(s) join the outlier set as input for the next iteration.
This check is useful not only for the above-mentioned case but also for rooms with a concave shape, where the straight line approximating one side can intercept small groups of points located in distant parts of the room.These points are false-inliers and are removed from the inlier set in the same way.The choice of the RANSAC threshold parameter depends primarily on the sensor performance in modelling the upper portion of the wall as flat.In the case of a zigzag pattern, typical of low-cost range cameras, a working threshold can be selected only with proper experience.Instead, point clouds collected by laser scanner present a more realistic flat pattern and require smaller thresholds, typically under 5 cm.If the threshold selected is too high, there is the concrete risk that the inlier set includes not only the points belonging to the side of one room, but also few points of the bordering two sides, a circumstance that removes some useful inliers from the estimation of those sides.In order to increase the accuracy, the final RANSAC result is refined by estimating the line parameters only on the inlier set through the Ordinary Least Squares method.
Eventually, the quality of the overall estimation after n iterations is measured by the number of outliers left after the last iteration: for many point clouds it is impossible to bring their number to 0, but the smallest achievable number corresponds in general to the best layout estimation.Hence, the entire estimation (n iterations) described earlier, is repeated indefinitely until the number of final outliers is equal or minor to nestimation/α , where α is a parameter related only to the available computing power.After α estimations, the number of acceptable final outliers increases by 1.Moreover, there are outliers consistently excluded by all RANSAC iterations: they derive from errors in the scanning process, corresponding to objects or physical barriers located in the ceiling-wall border, or to small curved portions of the 3D point cloud (corresponding often to real angles in the real room).To reduce the impact of these points on the algorithm runtime, the points that remain outliers for the totality of the first α estimations are removed from the input set of the subsequent estimations; same goes, if necessary, for steady outliers occurring in the second or third section of α estimation, although this is a rarer circumstance.
Once the final set of inliers is found, a raw layout is available (Figure 15): it consists of one line equation and two endpoints expressing the position of the inlier set, for each side of the room.
Finally, the actual layout is produced extending the line segments along their lines, two by two, up to the point of their intersection (Figure 16).

RESULTS
As mentioned before, the algorithm developed was tested on the 3D point cloud of 600.000 vertices selected from the ISPRS dataset on indoor modelling.The obtained layout is reported in Figure 17: the results are very promising, the algorithm is able to model effectively the 3D shape of the investigated room, at least from a visual inspection.Furthermore, a height of 2.67 m was estimated for the input point cloud.As the computational efficiency, the algorithm average runtime on the input point cloud is, at the moment, five minutes on a Desktop PC with an IntelCore i7 3 GHz CPU and 16 GB of RAM.Both the runtime and the system requirements can be considerably narrowed with informatic expertise: the parallelization of the code and the switch to more efficient and fast languages could optimize the CPU use, while a careful selection of the Data Structure could avoid bottlenecks and an excessive usage of RAM resources.Furthermore, the portability of the software is threatened only by the amount of data that have be processed in the 3D section of the algorithm, even though it is possible to consider a reduction of the cloud density, based for instance on a required minimum distance between every couple of points.

CONCLUSIONS AND FUTURE WORK
In this paper, a new algorithm for the automatic extraction of the layout and the height of a small environment was presented.In particular, the algorithm was tested on a point cloud of 600000 vertices, selected from the ISPRS benchmark on indoor modelling (Khoshelham et al., 2017).
The preliminary results are encouraging: the 3D shape (layout and height) of the investigated room is effectively reconstructed.However it is still necessary to deepen the analysis, by testing the algorithm on more and different environments, performing also a quantitative accuracy assessment, using for instance the evaluation criteria described in (Khoshelham et al., 2018).
Furthermore, at this stage, only room with polygonal layout and flat surfaces for walls and ceiling are taken into consideration, but thanks to the adaptability of RANSAC to different models in both its 2D and 3D variants, curve-shaped room together with room characterized by slanted walls or ceiling, could soon become feasible input.Finally, a possible future development is to release the implementation of the algorithm as a Free and Open Source Software (FOSS).

Figure 3 :
Figure 3: Vertical section of the input point cloud and determination of the room heigth (h = |zmin|).

Figure 7 :
Figure 7: Final reduced point cloud: (a) overview, (b) details. triangulation of the 2D point cloud, even though a smaller MEL corresponding to the value of the statistical mode could produce an equally optimal result.

Figure 8 :
Figure 8: Delaunay Triangulation: (a) detail; (b) influence of the value of the MEL.

Figure 9 :
Figure 9: Boundary polygonal chains for the input point cloud.

Figure 10 :Figure 11 :
Figure 10: Details of the smaller polygonal chains with relative portions of the Delaunay triangulation.

Figure 13 :
Figure 13: RANSAC iterations for identifying the sides of the room.

Figure 14 :
Figure 14: Example of an uncorrect identification of the side.

Figure 16 :
Figure 16: Inlier locations along the border of the polygonal chain.In black the outliers.