Automatische puntenwolk naar BIM conversie: van scan naar design

scan naar design

Automatische puntenwolk naar BIM conversie: van

Academiejaar 2019-2020

Master of Science in de industriële wetenschappen: elektronica-ICT Masterproef ingediend tot het behalen van de academische graad van

Begeleiders: dr. ir. Michiel Vlaminck, prof. Hiep Luong

Promotor: prof. dr. ir. Wilfried Philips

Studentennummer: 01609227

De auteur en promotor geven de toelating deze scriptie voor consultatie beschikbaar te stellen en delen ervan te kopiëren voor persoonlijk gebruik. Elk ander gebruik valt onder de beper-kingen van het auteursrecht, in het bijzonder met betrekking tot de verplichting uitdrukkelijk de bron te vermelden bij het aanhalen van resultaten uit deze scriptie.

The author and promoter give the permission to use this thesis for consultation and to copy parts of it for personal use. Every other use is subject to the copyright laws, more specically the source must be extensively specied when using from this thesis.

Gent, Mei 2020

De promotor De begeleiders De auteur

Woord vooraf

Bij het schrijven van deze thesis heb ik veel bijgeleerd over het implementeren van verschillende algoritmes, alsook van het optimaliseren hiervan. Daarnaast heb ik ook veel bijgeleerd over het voeren van een literatuuronderzoek, waardoor enkel wetenschappelijke bronnen gebruikt zouden worden bij het verder uitvoeren van de opdracht. Zelfstandig werken, kritisch denken en creatief op zoek gaan naar een oplossing zijn vaardigheden die ook zijn geperfectioneerd door het werken aan dit onderzoek.

Het resultaat zoals het er nu is, zou ik niet bereikt hebben zonder de hulp van bepaalde mensen. Daarom dank ik graag nog de volgende mensen:

Prof. Hiep Luong, voor het opvolgen van de progressie en het bieden van extra inzichten. Dr. ir. Michiel Vlaminck, ook voor de opvolging en de geboden hulp.

Dr. ir. Filip Rooms, voor de inzichten.

Prof. dr. ir. Wilfried Philips, voor het promoten van de thesis.

Aaron Mestdagh Gent 31 mei 2020

Automated point cloud to BIM conversion: from

scan to design

Aaron Mestdagh

Supervisor(s): prof. dr. ir. Wilfried Philips, dr. ir. Michiel Vlaminck, prof. Hiep Luong

Abstract—This article tries to extend the research on automating the pro-cess of transforming point cloud data to Building Information Models. This is done by first detecting planes in the data. These planes are at first solely represented by their mathematical equation. Next, extra constraints are applied to these planes to more accurately represent a wall. In a third step, holes are introduced that represent door openings and windows. Lastly, the results were evaluated both visually and using more objectively defined metrics.

Keywords— Point clouds, BIM, Building Information Model, plane de-tection, plane segmentation

I. INTRODUCTION

A

When people want to do heavy renovations on a building, an architect might want to plan this all out and he will do this us-ing a BIM. The problem that emerges here is the fact that this technology is quite modern and a lot of existing buildings do not have a BIM available. The process of creating this model how-ever, is very labor intensive. It would be necessary to measure the dimensions of all structures manually, which would be very time consuming.

The solution to this problem is to automate this process by ac-quiring a scan on the building using technology such as LIDAR. The data retrieved by these devices are point clouds and could be transformed to a building information model.

II. STATE-OF-THE-ART

As a start of this project, information about the evolutions that were already made on this topic was collected. The approach that was attempted here, was to start off by firstly transform-ing the data to planes while keeptransform-ing efficiency of the algorithm in mind. So research was done on the different ways to detect planes as well as how the data could be represented in a way that increases execution times as much as possible. Lastly, the actual segmentation of structural elements needed to be done. This was researched as well.

A. Plane detection

There are several ways to perform plane detection in point cloud data. One of the most common and relatively simple ways of doing this is using the RANSAC algorithm [1]. Another way is to use Hough transform [2]. The third way is using cluster-ring. A clustering algorithm [3] can be used to group all of the points in the data to match planes.

B. Datastructures

There were two major types of storing the data considered. One of these structures was an octree [4]. This is a structure that is heavily used in 3D processing environments. It has also been used specifically for point cloud applications as well in [5] and [6].

The next structure that was also considered is k-d tree [7]. This is basically a binary tree applied in three dimensions. It is possible to have this in more dimensions as well. This struc-ture allows a quick search based on location, which is a very common problem in point clouds. This search would otherwise require an iteration over every point.

C. Segmentation

There were different ways of segmenting the structural ele-ments of a BIM. One of the ways of doing this, was to work one floor at a time while projecting points onto twodimensional planes representing the floor [8]. This work was also continued and perfected later [9].

Another way that this segmentation could be done was using the smoothness of the surface measured with point clouds [10]. This method was developed to function in an industrial environ-ment and to detect all types of shapes with a smooth surface.

III. DATA PREPARATION

Before starting the actual segmentation of geometrical shapes, the scanning data must first be prepared to deal with the algo-rithms being used. One of the problems of the data is the fact that it could be noisy. The distribution of the data in 3D space is also not uniform. Both of these problems are caused by ei-ther imperfections on the measuring device, or just because of the way it works. A LIDAR sensor for example will perform a circular scan around itself. This means that areas close to the sensor will be more densily sampled than areas far away from the scanner.

An example of the data that was used in this work on can be seen in Fig. 1.

To solve this problem, an algorithm called ‘Voxel downscal-ing’ was applied to the data. This will divide the entire 3D space in cells, otherwise known as ‘voxels’. Once this discretisation is done, all points in the voxels are averaged to result in just one point for every voxel. This highly reduces the density on the areas close to the scanner while also averaging out the noise in the dataset.

This way, the amount of data also lowers a significant bit. This improves execution times, which can make for a more real time experience.

Fig. 1. Point cloud of an appartment, stitched together from RGB-D images, so every point also has a color.

IV. PLANE DETECTION

As stated earlier in section II-A, there are several ways of de-tecting planes in point clouds. Because this set deals with mul-tiple planes and because a determenistic approach was prefered, clustering was chosen detect the planes all at once.

A. Clustering

The algorithm that was chosen to detect all the planes in the dataset is the clustering algorithm. To do this, multiple clus-tering algorithms were considered, including the agglomerative hierarchical clustering algorithm proposed in [3]. Though DB-SCAN scan, which is short for ‘Density-Based Spatial Cluster-ing of Applications with Noise’, was chosen here. DBSCAN is chosen for its unique properties that work with the data given [11]. The size of the data is quite large, so an algorithm with high efficiency was needed to divide the points in clusters and DBSCAN offers a good efficiency. Other than that, the algo-rithm should also be able to separate clusters of all shapes from each other, which is a feature most clustering algorithms do not provide. Lastly, the data still deals with a certain amount of noise, despite the voxel downscaling. There are also some mis-takes present in the dataset due to reflections in the glass of win-dows for example. DBSCAN is able to deal with both of these problems because it clusters based on the density and these mis-takes are most often low in density. When chosing the parame-ters of DBSCAN, a minimum size of clusparame-ters can also be cho-sen, which eliminates small mistakes like those reflections for example.

B. Normal estimation

Clustering in 3D space would not be sufficient for DBSCAN to work in this case: it needs extra features to be able to fit it in a plane. Just calculating a normal vector for every point in the dataset seemed to be sufficient to describe a point and cluster them using DBSCAN into planes. These normal vectors were determined by getting the neighbouring points calculating a lo-cal plane using the least squares method [12].

This results in 6 values that were used to describe a point, be-ing the 3D location as well as an additional 3 values to describe the normal vector of that point.

C. Post processing

After the clustering is performed, there were still artifacts present in the data. Even when the best parameters were se-lected after experimenting with DBSCAN, there were still small clusters present in the data that do not represent a wall or floor, as well as several big clusters that were indexed separately while they should be grouped together. A last problem was that par-allel planes that were close to each other were also grouped to-gether while they should be separated.

These problems were fixed using four steps. The first step was to remove clusters that do not represent planes. When there was a smooth transition of normal vectors between points, there was a possibility that these would be grouped together despite not being part of a plane. To do this, a plane determined us-ing the least squares method is fitted onto the cluster of points and the average distance of the points to this plane is calculated. Next, small clusters that contained an insignificant amount of points were removed from the data. Thirdly, different groups that were roughly located in the same plane and not too far from each other were grouped together in the same cluster. As a last step, RANSAC was applied on the clusters separately to detect parallel planes that were too close to each other to separate using the clustering algorithm.

An example of how this clustering eventually performed is shown in Fig. 2.

Fig. 2. Result of the clustering algorithm. Every cluster of points is labeled with a different color.

V. SEGMENTATION

After dividing the points in appropriate groups, they were used to further segment this data into actual planes with con-straints, so they can represent the structural elements in a BIM. These elements being walls and floors.

Next, the ground floor was detected by using the largest clus-ter as plane that represents the ground floor. Afclus-ter declus-termining the ground floor, the other most present directions of the walls were determined. Before further segmentation was done, the entire point cloud was transformed using a rotation matrix [13] constructed by using the vectors determined in the previous step to align the building walls with the axes of the Cartesian system.

A. Walls

By aligning the floors and walls with the axes of the system, the different floor levels could be determined quite easily. Walls can now be determined for each floor separately. To do this, the planes are projected as lines onto the floor surfaces [8]. This results in what can be seen in Fig. 3. These lines are supposed to be extended so they represent complete walls that touch each other. This is done by calculating all possible intersection points and filtering down these points until there are no walls intersect-ing with each other. Walls cannot go through another wall, they can only be extended to touch that wall.

Fig. 3. Projection of walls onto a 2D plane, in this case this is the ground floor.

B. Floors

To calculate the dimensions of the floors, the same intersec-tion points that were used to calculate the walls are utilized. These intersection points are used to create a grid of all possi-ble cells on that floor which could be occupied with a value that either represents an occupied floor cell or an empty floor cell. All occupied floor cells were then triangulated to have a cov-ered floor. These occupied cells also get compared to whether or not there are points present in this cell. If there are no points present in the cell, the cell is not marked as part of the floor. Fig. 4 shows the result of this algorithm executed on the ground floor of an appartment.

VI. HOLES

One remark is that these planes are all combinations of rect-angles that do support concave shapes as well. But they do not support having enclosed areas that do not contain any points. These are the so called ‘holes’ in the walls and could represent doors or windows. These holes are detected by creating a grid of every 2D surface of the planes, these grids can then be filled in by values that represent the presence of points in the cell. If there are a lot of cells that have no points, this spot will be recognized as a hole. With the holes taken into account, the representation of the model looks like what is shown in Fig. 5.

Fig. 4. Triangulated ground floor. Red lines mark the wallplanes, yellow marks the scanned points on the floor and the grey represents the occupied in area.

Fig. 5. Result of segmented walls and floors with calculated holes. The colors are the pointcloud clustered into groups to make the comparisson.

VII. EVALUATION

For the models being used, there was no ground truth avail-able. There was therefore no way of testing if the determined dimensions were actually correct. However, there is a way of describing if the segmented parts accurately describe the point cloud that was used to determine these structures. We used two metrics used to determine the accuracy. The first metric repre-sents the fraction of points that are assigned to a usable cluster that are still included in the structure. The second metric rep-resents the area that is actually covered by the points that are assigned to a cluster.

A. Original data

These metrics are first applied to the original data. The exper-iments demonstrated that the algorithm to determine the floor planes are quite accurate. The fraction of points used here could be close to 100%, but there are also cases where this is closer to 80%. The fraction of the area occupied resulted in about 85%. For walls, the results do not go as high: 75% of the points are still represented by a plane in that case and 72% of the area was occupied with points. All of the inaccuracies of these results can be explained by the fact that clutter in the dataset causes holes

of all kinds of shapes. These shapes are not supported by the current algorithm and causes the dimensions of holes to be esti-mated wrongly. It is even possible that there is not supposed to be a hole in that location, but there is no way of telling that with the current method.

B. Robustness

To test if noise, mistakes in data or extra clutter has a big im-pact on the algorithm, several robustness tests have been applied on the data. One test applied an extra random factor to the data: a uniform noise function relocated the points up to a maximum of 3 cm from their original location. In another test we added random noise. Deletion of a portion of data was also tested as well as the addition of random planes to test if it would be rec-ognized as clutter. None of these additions had a significant change in results on the metrics that were discussed in section VII-A. However, when the data gets too noisy, it is possible that extra walls are detected. These walls are the result of planes that become too thick, which results in a single plane being rec-ognized as two separate planes. So it is possible that some of the detected structures change when the parameters turn op too high, but overall, there is little change noticeable.

VIII. CONCLUSION

This article proposed a method to segment point clouds into planes that define the shape of the building. Walls and floors are as accurately represented as possible for the data that was being used. The algorithm proved to be quite resistant against noise and clutter. It should perform well enough for modern scanners. Several shapes are not supported, such as more rounded struc-tures, though it becomes really hard to distinguish between clut-ter and important shapes when every possible shape gets taken into account. A next step in this research would be to export the calculated information so it can be used in CAD software.

REFERENCES

[1] Michael Ying Yang, Wolfgang F¨orstner Plane Detection in Point Cloud Data Department of Photogrammetry, Institute of Geodesy and Geoinfor-mation, University of Bonn, 2010

[2] Dorit Borrmann, Jan Elseberg, Kai Lingemann, Andreas N¨uchter The 3D Hough Transform for Plane Detection in Point Clouds: A Review and a new Accumulator Design 3D Research Center and Springer, 2011

[3] Chen Feng and Yuichi Taguchi and Vineet R. Kamat ast Plane Extraction in Organized Point Clouds Using Agglomerative Hierarchical Clustering Mitsubishi Electric Research Laboratories, 2014

[4] Donald Maegher Geometric Modeling Using Octree Encoding Rensselaer Polytechnic Institute, 1981

[5] Anh-Vu Vo and Linh Truong-Hong and Debra F. Laefer and Michela Bertolotto Octree-based region growing for point cloud segmentation Else-vier, 2015

[6] Yun-Ting Su and James Bethel and Shuowen Hu Octree-based segmenta-tion for terrestrial LiDAR point cloud datain industrial applicasegmenta-tions Else-vier, 2016

[7] Michael Greenspan and Mike Yurick Approximate K-D Tree Search for Efficient ICP Dept. of Electrical and Computer Engineering, School of Computing, Queens University, 2003

[8] H´el`ene Macher and Tania Landes and Pierre Grussenmeyer Point Clouds Segmentation As Base For As-Built BIM Creation Photogrammetry and Geomatics Group, ICube Laboratory, 2015

[9] H´el`ene Macher and Tania Landes and Pierre Grussenmeyer From Point Clouds to Building Information Models: 3D Semi-Automatic Reconstruc-tion of Indoors of Existing Buildings Photogrammetry and Geomatics Group, ICube Laboratory, 2017

[10] T. Rabbani and F. A. van den Heuvel and G. Vosselman Segmentation Of Point Clouds Using Smoothness Constraint Section of Optical and Laser Remote Sensing, TU Delft, 2006

[11] Kamran Khan and Saif Ur Rehman and Kamran Aziz and Simon Fong and S. Sarasvady, DBSCAN: Past, Present and Future Department Of Computer Science SZABIST, 2014

[12] Estimating Surface Normals in Noisy Point Cloud Data Stanford Graphics Laboratory, Stanford University, 2003

[13] James Diebel Representing Attitude: Euler Angles, Unit Quaternions, and Rotation Vectors Stanford University, 2006

Inhoudsopgave

1 Inleiding 1

1.1 Building Information Model . . . 1

1.2 Probleemstelling . . . 2

1.3 Bricsys . . . 4

2 State of the art 5 2.1 Vlakdetectie . . . 5 2.1.1 RANSAC . . . 6 2.1.2 Hough transformatie . . . 10 2.1.3 Clustering . . . 13 2.2 Datastructuur . . . 17 2.2.1 Octree . . . 17 2.2.2 K-d tree . . . 18 2.3 Segmentatie . . . 20 2.3.1 Histogram aanpak . . . 20 2.3.2 Muren . . . 22 2.3.3 Gladheid . . . 22 3 Aanpak 23 3.1 Voorbereiding . . . 23 3.1.1 Exploratie . . . 24 3.1.2 Reductie . . . 26 3.1.3 Synthetische data . . . 31 3.2 Clustering . . . 37 3.2.1 Algoritme . . . 38 3.2.2 Optimalisatie . . . 50 3.2.3 Post processing . . . 55 3.3 Segmentatie . . . 64 3.3.1 Projectie . . . 64 3.3.2 Muurdetectie . . . 70 3.3.3 Vloerdetectie . . . 76 3.4 Gatendetectie . . . 83 viii

Inhoudsopgave ix 4 Evaluatie 86 4.1 Foutenintroductie . . . 86 4.1.1 Uniforme ruis . . . 86 4.1.2 Verwijdering . . . 87 4.1.3 Toevoeging . . . 88 4.2 Synthetische set . . . 90 4.2.1 Standaardmeting . . . 90 4.2.2 Robuustheidstesten . . . 95 4.3 gescande data . . . 99 4.3.1 Standaardmeting . . . 99 4.3.2 Robuustheidstesten . . . 101 5 Conclusie 104 5.1 Bespreking . . . 104 5.2 Verder werk . . . 105 Bibliograe 106

Hoofdstuk 1

Inleiding

Tegenwoordig werken architecten niet langer op papier bij het ontwerpen van een gebouw. Dit gebeurt allemaal digitaal, in CAD modellen. CAD staat voor `Computer Aided Design' en is dus eigenlijk een softwarepakket dat men gebruikt om verscheidene zaken in een computer te ontwerpen. Het onderwerp van deze thesis zal zich focussen op dit laatste. CAD wordt tegenwoordig verkozen boven het gebruik van papier omwille van het feit dat men bijvoorbeeld in 3D het geheel al veel beter en makkelijker kan voorstellen. CAD is daarnaast ook in staat om verschillende zaken automatisch te bepalen, iets wat voor een ontwerper zeker een goede zaak is. Ten slotte is het digitaal werken iets wat de aanpasbaarheid van het model, alsook het verspreiden ervan ten goede komt.

1.1 Building Information Model

Bij dit onderzoek staat BIM centraal, BIM staat afgekort voor `Building Information Model'. Een BIM is in essentie ook een CAD model, maar bevat nog heel wat extra informatie. Een gewoon CAD model bevat vooral structurele en geometrische eigenschappen, daarnaast kan het ook nog materiaaleigenschappen bevatten. Een building information model zal echter alle informatie bevatten die een architect mogelijks zou kunnen nodig hebben voor het construe-ren van een gebouw. Deze zal dus dynamischere eigenschappen bevatten, zoals bijvoorbeeld energieverbruik. Het is nuttig om te weten of een muur bijvoorbeeld een steunmuur is. Een extra toevoeging op CAD is dat BIM een structuur heeft die zodanig ontworpen is dat ver-anderingen in het model zeer eenvoudig en snel gemaakt kunnen worden. Alle objecten zijn daarin volgens de beperkingen van een gebouw aan elkaar gelinkt. Als er dan iets verandert, zal alles wat hiermee te maken heeft ook geüpdatet worden. Visueel komt dit neer op wat bijvoorbeeeld in guur 1.1 te zien is. Dit is hier gewoon een 3D voorstelling van een gebouw. Maar achterliggend is nog veel extra informatie aanwezig.

Hoofdstuk 1. Inleiding 2

Figuur 1.1: 3D model van een constructie (BIM) met twee verschillende voorstellingen. Bron: Cherkaoui (2016)

1.2 Probleemstelling

Het is nu bekend dat architecten gebruik maken van digitale software om BIM's te creëren alvorens men zou beginnen aan het bouwen van een constructie. Een probleem dat zich nu voordoet is het feit dat enorm veel gebouwen al geconstrueerd zijn zonder dat hier eerst een BIM model voor gemaakt is. Dit is relatief nieuwe technologie en dus zijn er nog enorm veel gebouwen die niet voorzien zijn van dergelijk model. Een situatie waarin dit wel heel handig zou zijn, is het geval waarbij men een gebouw zou willen verbouwen. Men wil dan vaak structurele eigenschappen wijzigen en dit is vaak niet doenbaar zonder eerst de juiste overwegingen te maken. Om hier zo goed mogelijk op in te spelen moet er dus een BIM gecreëerd worden van het gebouw. Echter, er zijn heel veel afmetingen nodig om zo'n model te kunnen opstellen, men moet de exacte dimensies van het gebouw weten en invoeren om een 3D model van het gebouw te creëren. Kan dit deel niet geautomatiseerd worden? Het probleem die zich hier dus voordoet is het feit dat architecten nood hebben aan een BIM van bestaande gebouwen, maar dat dit een tijdrovend proces is. Dit tijdrovend proces kan minstens voor een deel geautomatiseerd worden.

Hoofdstuk 1. Inleiding 3 In plaats van alles manueel op te meten, kan er gebruik gemaakt worden van reeds bestaande technologie om een gebouw te `scannen'. Om geometrie van iets op te meten is er nood aan een sensor die diepte kan opmeten. Een welbekend toestel dat hiervoor kan gebruikt worden is een LIDAR sensor. LIDAR is de afkorting van `LIght Detection And Ranging' en zal met behulp van lasers de afstand van een bepaald punt tot de sensor meten. LIDAR werkt op basis van het Time Of Flight principe om afstanden op te meten. Deze zal dus lichtpulsen uitsturen en meten hoe lang het duurt vooraleer dat de reectie de sensor terug bereikt:

d= c· t 2· n

Bovenstaande formule illustreert hoe dit in zijn werk gaat. Hierbij is d de afstand, c de snelheid van het licht, t de tijd vooraleer de lichtpuls wordt opgemeten door de sensor en n de brekingsindex. Een LIDAR sensor gaat dit proces herhalen voor meerdere `pixels'. De sensor kan hierbij aanzien worden als een soort camera die volledig om zich heen kijkt, maar deze meet diepte op in plaats van kleur. Het resultaat is een 360 graden beeld van dieptepunten. Dit kan dan op zijn beurt omgezet worden naar een eenvoudiger Cartesiaans stelsel.

Uiteindelijk zal hiermee een zogenaamde puntenwolk bekomen worden. Dit is simpelweg een verzameling van de 3D punten die door de sensor zijn opgemeten. Een voorbeeld van zo'n puntenwolk is te zien in guur 1.2.

Figuur 1.2: Puntenwolk gecapteerd met een LIDAR sensor. Bron: Brooks (2017)

Hoofdstuk 1. Inleiding 4 In bovenstaande guur zie je het resultaat van één scan. Hierin is duidelijk een circulair patroon herkenbaar. Men kan visueel uitmaken waar de scanner zich bevond bij het maken van deze opmeting. Eén scan biedt vaak echter te weinig informatie. Daarom zijn de scans die gebruikt werden in deze thesis meestal combinaties van meerdere scans die aan elkaar gestitcht zijn.

Merk op dat er zeker meerdere soorten hardware zijn om dergelijke soorten data te verzame-len. De datasets die hier besproken worden zijn wel steeds opgemeten met behulp van een LIDAR sensor. Hierbij werd rekening gehouden met de standaard specicaties van zo'n toe-stel, zoals bijvoorbeeld de resolutie en accuraatheid ervan. Met behulp van bijvoorbeeld een dieptecamera of nog andere manieren om diepte op te meten kan men evengoed puntenwolken creëren.

1.3 Bricsys

Deze thesis werd gemaakt in samenwerking met Bricsys. Bricsys is een bedrijf dat software levert. Meer speciek CAD software, genaamd BricsCAD. BricsCAD is een software die op zijn beurt uit nog extra onderdelen bestaat. Eén van die onderdelen is gefocust op het maken van Building Information Models, namelijk BricsCAD BIM. Dat is dan ook de use case waarop hier gewerkt zal worden. Bricsys heeft reeds onderzoek gevoerd naar het automatiseren van het `Scan to BIM'-proces en hebben daarbij machine learning gebruikt om van ruwe geometrische vormen naar eectieve gebouwstructuren te gaan. De onderliggende onderdelen zijn echter nog niet aanwezig. Dat is dan ook waar de focus van dit onderzoek ligt.

Figuur 1.3: Logo Bricsys

Hoofdstuk 2

State of the art

In dit hoofdstuk wordt er besproken welke manieren er reeds bestaan voor het verwerken van puntenwolken tot BIM's. Op het laagste niveau moeten de punten eerst ingedeeld worden in vlakken. Dit zal namelijk de eerste stap vormen naar het creëren van muren en vloeren. In een latere stap moeten deze vlakken, die op dat punt met slechts heel minimale parameters beschreven zijn, eectieve gebouwstructuren vormen. Onder deze gebouwstructuren worden vloeren en muren met concrete dimensies verstaan.

2.1 Vlakdetectie

Het was eerst de opdracht om de punten zodanig in te delen dat deze behoren tot eenzelfde vlak. Dit kan op verschillende manieren gebeuren, elk met zijn voordelen. Om dit probleem op te lossen werden drie manieren gevonden die het vaakst gebruikt werden, namelijk RANSAC, Hough transformatie en clustering.

Hoofdstuk 2. State of the art 6

2.1.1 RANSAC

Wellicht de manier die het vaakst werd teruggevonden in de literatuur is RANSAC. RANSAC is een afkorting en staat voor 'RANdom SAmple Consensus'. Het is een methode die alge-meen kan gebruikt worden om vanalles te gaan bepalen. Het wordt gebruikt om in feite de parameters van een bepaald model te schatten. Hier zou dit model een vlak kunnen zijn, maar RANSAC wordt bijvoorbeeld ook gebruikt voor het stitchen van data. Dit kunnen puntenwol-ken zijn, maar men kan dit ook breder bekijpuntenwol-ken. Zo kunnen bijvoorbeeld ook foto's aan elkaar gestitcht worden met behulp van RANSAC. Het bestaat uit een iteratief proces, waarbij het aantal iteraties instelbaar zijn. Hoe meer iteraties, des te beter de schatting zal zijn, voordat het geheel convergeert. Bij elke iteratie zal RANSAC een aantal elementen uit de dataset nemen die willekeurig gekozen zijn, vandaar random sample consensus. Met behulp van deze elementen zullen de parameters van het model bepaald worden. Daarna wordt dit model op de rest van de dataset geverieerd en wordt er als het ware een score gegeven van hoe goed dit model past bij de beschikbare dataset. Dit wordt dus een aantal keer herhaald en uiteindelijk worden de parameters met de beste score gekozen.

RANSAC heeft een enorm voordeel in het feit dat het goed met outliers overweg kan. Als dit dus wordt toegepast op een dataset met veel ruis, dan zal het vlak nog steeds gedetecteerd worden. Of misschien nog belangrijker, als er meerdere vormen en vlakken aanwezig zijn, dan zal RANSAC er nog steeds in slagen aparte vlakken te detecteren.

(Yang & Förstner, 2010) beschrijft hoe RANSAC in zijn werk gaat bij het detecteren van vlakken. Hierbij zijn meerdere nog steeds meerdere implementaties mogelijk, maar de een-voudigste en waarschijnlijk ook de snelste implementatie gaat steeds slechts rekening houden met 3 punten om de parameters te bepalen. Om de parameters van een vlak te vinden zijn er minstens 3 punten nodig, hierbij kunnen er meerdere in rekening gebracht worden, maar dit zal dan een licht andere aanpak vereisen. Na het bepalen van de parameters met behulp van de 3 punten, kunnen de overige punten in de dataset vergeleken worden of deze passen in het vlak. De eenvoudigste manier om dit te bepalen is gewoon door de afstand te berekenen tot het vlak. Punten die binnen een bepaalde threshold vallen, behoren tot het vlak. En het vlak waartoe het meeste punten behoren, wordt dan gekozen als het correcte vlak. Daarna kunnen de parameters nogmaals opnieuw geschat worden waarbij alle punten in de consen-sus gebruikt worden. Dit kan dan bijvoorbeeld gebeuren aan de hand van de least squares methode. Algoritme 1 beschrijft in pseudocode hoe dit in zijn werk gaat.

Hoofdstuk 2. State of the art 7 Algoritme 1: RANSAC voor vlakdetectie

bestSupport = 0;

bestPlane(3, 1) = [0, 0, 0]; bestStd = ∞;

planepoints = []; i = 0;

e = 1 - forseeable-support/length(point-list) ; // forseeable − support is het maximum waarschijnlijk aantal punten die tot eenzelfde vlak kan behoren N = round(log(1 - α) / log(1 - (1 − e)3_{)) ; // α is de minimum probabiliteit dat}

er minstens 1 goede observatie gevonden wordt in N iteraties while i ≤ N do

j = pick 3 points randomly among (point-list); pl = pts2plane(j);

dis = dist2plane(pl, point-list); s = nd(abs(dis)≤t);

if length(s) < bestSupport or (length(s) = bestSupport and st < bestStd) then bestSupport = length(s); bestPlane = pl; bestStd = st; planepoints = s; end i = i + 1; end

Merk op dat dit algoritme ook rekening zal houden met de spreiding van de punten. Indien twee combinaties van punten een vlak beschrijven die evenveel punten bevatten die tot het vlak behoren, dan zal de combinatie die het minst spreiding heeft verkozen worden. Dit wordt hier voorgesteld met behulp van standaarddeviatie.

Dit beschrijft hoe er één enkel vlak gedetecteerd kan worden. In de aangeleverde datasets zijn er echter veel meer vlakken aanwezig. Volgens dit algoritme zal het grootste vlak van de dataset bekomen worden in zo goed als alle gevallen. Het aantal iteraties moet immers hoog genoeg gehouden worden. Met het grootste vlak wordt hier het vlak bedoeld waartoe het meest aantal punten behoren. Dit vlak kan nadien uit de dataset verwijderd worden en men kan dan vervolgens op de aangepaste dataset dit algoritme herhalen totdat er geen vlakken meer gevonden worden.

Hoofdstuk 2. State of the art 8 Algoritme 2: Meerdere vlakken detecteren

planes = [ ]; planepoints = [ ]; planesize = ∞;

while planesize > t do

plane, points = ndPlane(point-list) ; // findP lane verwijst hier naar het eerder beschreven RANSAC algoritme

planesize = length(points); if planesize > t then planes.append(plane); planepoints.append(points); end end

Algoritme 2 beschrijft concreet de werking hiervan. Men kan hierbij instellen met behulp van parameter t hoe groot een vlak ten minste moet zijn alvorens het herkend zal worden als vlak.

Figuur 2.1: Detectie van dominante vlakken (in het rood aangeduid) in de dataset (in het groen aangeduid) voor twee opeenvolgende iteraties waarbij telkens het reeds gedetecteerde vlak uit de dataset verwijderd wordt.

Figuur 2.1 illustreert het algoritme. Links is een eerste vlak gedetecteerd, aangeduid in het rood, dit stelt de eerste iteratie van het algoritme voor. Rechts kan opgemerkt worden dat het gevonden vlak uit de eerste iteratie niet langer aanwezig is in de dataset en dat er dus een nieuw vlak gedetecteerd wordt.

Hoofdstuk 2. State of the art 9 RANSAC heeft ook nog enkele varianten die hiervan zijn afgeleid. MLESAC Torr & Zisser-man (1996) is hier een voorbeeld van. MLESAC is de afkorting van `Maximum Likelyhood Estimator SAmple Consensus'. Dit algoritme zal een minder willekeurige keuze maken dan RANSAC door het feit dat dit algoritme me probabiliteiten werkt. Daar waar RANSAC werkt met het aantal inliers om te bepalen of het model het beste is, zal MLESAC de probabiliteit dat dit model het beste model is hanteren. Ook MAPSAC, wat afgekort staat voor `Maximum A Posteriori SAmple Consensus' gebruikt een gelijkaardig principe van probabiliteiten om een betere schatting te maken van de samples die gebruikt zullen worden. MAPSAC is eigenlijk een uitbreiding op MLESAC. Ook zijn er variaties die nog gelijkaardig aan RANSAC werken, maar voor meerdere modellen tegelijk gebruikt worden. Lara-Alvarez et al. (2009) illustreert een variant op RANSAC die met meerdere modellen tegelijk overweg kan.

Hoofdstuk 2. State of the art 10

2.1.2 Hough transformatie

Hough transformatie is een welbekende manier om geparametriseerde vormen te zoeken in een dataset. Daarbij is vooral het vinden van lijnen en cirkels in een 2D setting bekend. Hough transformatie is echter ook een techniek die kan toegepast worden in 3D. Dit wordt dan uiteraard uitgevoerd op een puntenwolk, in 2D komt het vaak voor dat dit uitgevoerd wordt op pixels. Borrmann et al. (2011) beschrijft de werking van Hough transformatie voor vlak detectie.

Daarbij wordt het geheel getransformeerd naar een ander coördinatenstelsel rekening houdende met de parameters van het vlak. Een vlak in 3D ruimte wordt normaal als volgt voorgesteld:

a· x + b · y + c · z = ρ

Deze vorm wordt ook wel de `Hesse normaal' vorm genoemd. Hierbij wordt een vlak voor-gesteld door een punt in het vlak met coördinaten [x, y, z], een normaalvector [a, b, c] die loodrecht staat op het vlak in kwestie en dan nog een parameter ρ die de afstand tot de oorsprong voorstelt. Om de Hough transformatie uit te voeren, wordt deze formule nog iets anders geïnterpreteerd, daarbij wordt de normaalvector in poolcoördinaten omgezet, uitein-delijk wordt dan volgende vergelijking bekomen:

x· cos(θ) · sin(φ) + y · sin(θ) · sin(φ) + z · cos(φ) = ρ

Figuur 2.2: Normaalvector poolcoördinaten, het coördinaat aan het einde van de rode pijl is de normaalvector in het Cartesiaans stelsel. In poolcoördinaten wordt dit φ, θ en ρ. Bron: Nordmann et al. (2018)

Hoofdstuk 2. State of the art 11 Figuur 2.2 toont aan hoe de conversie van een Cartesiaans stelsel naar poolcoördinaten in zijn werk gaat. Wanneer men elk punt nu in deze formule zou invoeren bekomt men niet langer een reeks punten, maar eerder een reeks 3D sinusoïden. Deze sinusoïden kunnen geïnterpreteerd worden als alle mogelijke parameters die een punt door dat vlak zou kunnen hebben. Daarbij variëren immers de twee hoeken van de normaalvectoren die in de bovenstaande formule reeds besproken zijn.

Een snijding tussen 2 van deze sinusoïden stelt een rechte voor die de twee punten hiervan verbindt. Wanneer er 3 van deze graeken elkaar snijden kan er gesproken worden van een snijpunt op de parameters die een vlak beschrijft waartoe alle 3 deze punten behoren. Figuur 2.3 maakt dit duidelijk.

Figuur 2.3: Sinusoïden gevormd door Hough transformatie. De rode, groene en blauwe graek stellen respectievelijk de sinusoïdes voor voor de punten [1, 0, 0], [0, 1, 0], [0, 0, 1]. De zwarte stip illustreert het snijpunt van de 3 graeken die dan θ, φ en ρ van het vlak bepaalt.

Hoofdstuk 2. State of the art 12 Nu kan men inzien dat elke 3 punten in de dataset een vlak kunnen vormen, ongeacht of deze visueel een vlak vormen met punten die een hoge dichtheid bezitten. Er zullen dus heel veel van deze snijpunten aanwezig zijn. Wat hier nu moet gebeuren, is het zoeken van de plekken waar heel veel van deze snijpunten dichtbij elkaar aanwezig zijn. De densiteit van het aantal snijpunten moet hiervoor onderzocht worden. Een manier om dit te doen is het discretiseren van de hough ruimte en zo cel per cel te kijken waar het meeste punten aanwezig zijn. Er kan zo een instelbare drempelwaarde zijn die bepaalt wanneer er genoeg punten zijn zodanig dat dit als een vlak aanzien kan worden. Het discretiseren zorgt er dan wel voor dat de parameters minder exact bepaald worden. Hier zal een extra stap nodig zijn om met de gevonden punten nogmaals de parameters apart te berekenen zodanig dat een goed resultaat bekomen wordt. Een andere manier om de vlakken nu te vinden is een densiteitsgebaseerd clustering algoritme uitvoeren op de berekende snijpunten, deze zal dan alle plaatsen met het meest aantal snijpunten eruit halen en zo kan ook het centrum van de cluster bepaald worden die wellicht een accurate schatting is van de parameters.

Hoofdstuk 2. State of the art 13

2.1.3 Clustering

Een derde optie om de punten te groeperen in te delen die vlakken voorstellen, is een cluste-ringmethode. Daarvan wordt in de literatuur iets minder gevonden, maar Feng et al. (2014) beschrijft toch een manier die gefocust is op het segmenteren van vlakken in een puntenwolk op zeer snelle manier. Deze snelheid en het feit dat het een eerder deterministische aanpak is om de vlakken te vinden, is wat dit algoritme zeker interessant maakt.

AHC

AHC is hier wat voorgesteld wordt, dit staat voor `Agglomerative Hierarchical Clustering'. Deze manier van werken gaat eigenlijk steeds punten gaan combineren en zo bottom-up steeds grotere vlakken bekomen. Dit gaat dus eigenlijk omgekeerd in zijn werk dan hoe RANSAC bijvoorbeeld zou werken, waar een vlak gecreëerd wordt en dan achteraf de punten ermee vergeleken worden.

Figuur 2.4: Principe van Agglomerative Hierarchical Clustering. Bron: Feng et al. (2014)

Hoofdstuk 2. State of the art 14 Figuur 2.4 illustreert hoe dit in zijn werk gaat voor 2D punten. Typisch zal AHC afstand gebruiken als parameter voor het groeperen van datapunten. Daarbij zijn er nog verschillende keuzes van hoe de afstand nu net bepaald kan worden. Figuur 2.5 toont de verschillende mogelijkheden.

Figuur 2.5: Verschillende manieren om afstand te bepalen, 1. single-link 2. centroid 3. complete-link.

Echter, wanneer gewerkt wordt met vlakken kan de gewone afstand tussen clusters niet ge-bruikt worden als metriek, punten kunnen immers ook dichtbij elkaar liggen zonder dat deze tot hetzelfde vlak behoren. De metriek die in deze paper daarom gebruikt wordt is MSE, oftewel `Mean Squared Error'.

Het algoritme zal dus steeds buren samenvoegen en berekenen of de MSE niet te groot wordt. Indien dit wel het geval is, dan zal er niets samengevoegd worden. Dit is het voorgestelde algoritme:

Hoofdstuk 2. State of the art 15 Algoritme 3: Agglomerative Hierarchical Clustering

point-groups = group for every point in point-list; smallest-mse = 0;

while smallest-mse < t do temp-smallest-mse = ∞; for group in point-groups do

neighbours = group.getNeighbours(); for neighbour in neighbours do

mse = calcMSE(neighbour, group); if mse < t then

Merge(neighbour, group); end

if mse < temp-smallest-mse then temp-smallest-mse = mse; end end end smallest-mse = temp-smallest-mse; end

Een probleem dat zich bij dit algoritme voordoet is vooral het feit dat de buren moeten bepaald worden. In de literatuur is dat geen probleem, want de data is daar geordend. Het gaat hier om een dieptebeeld, de punten worden hier dus eigenlijk voorgesteld in een tweedimensionaal rooster, zoals een afbeelding. Buren zijn in dergelijk rooster heel eenvoudig te vinden. Men hoeft in dat geval gewoon de naburige cel aan te spreken.

Wat ook opvallend is in dit onderzoek, is dat de datasets niet heel groot zijn. Er werd gebruik gemaakt van een kinect sensor, die een resolutie heeft van 640x480 pixels, wat toch neerkomt op een relatief kleine dataset. In dat geval is real-time snelheid natuurlijk wat eenvoudiger te garanderen. Dit werd ook niet bewezen in een 360 graden omgeving. De kinect camera heeft een gezichtsveld van 57 graden horizontaal en 43 verticaal, het is dus niet zeker of dit in elke omgeving en met meer data nog steeds zou werken. Het zou ook vereisen dat de data geordend moet worden.

Het principe is wel nog steeds interessant, er zijn namelijk andere clusteringtechnieken die gebruikt kunnen worden die niet op deze limieten stoten. Ook zijn er manieren om het vinden van buren te versnellen indien er geen geordende data beschikbaar is.

Hoofdstuk 2. State of the art 16 Normaalschatting

Indien er andere clusteringtechnieken gebruikt worden die wel met afstandsmetrieken werken, dan is de informatie van enkel de 3D positie niet voldoende om hier een correcte clustering uit te voeren. Daarom werd op zoek gegaan naar extra eigenschappen die dit beter kunnen bepalen. De meest voor de hand liggende, en eigenlijk enige extra eigenschap die bevestigt ofdat een punt in een bepaald vlak ligt, is de normaalvector.

Als men een puntenwolk bekijkt is er immers een richting waarin een punt kan wijzen door rekening te houden met de omliggende punten. Zhao et al. (2019) beschrijft een proces om normaalvectoren te bepalen die ook speciek toegepast is op de data verkregen door LIDAR sensoren, wat hier ook de focus is. Daar wordt beschreven dat normaalvectoren in punten-wolken doorgaans bepaald worden door in de omgeving van elk punt te kijken en telkens een lokaal vlak te vormen die met behulp van de least squares methode bepaald wordt. De nor-maalvector van dit vlak is dan in essentie de nornor-maalvector van dit punt. Figuur 2.6 illustreert het gewenste resultaat.

Figuur 2.6: Ideale geval van normaalschatting. De blauwe lijn stelt het vlak voor, de zwarte punten de eigenlijke data en de rode pijlen zijn hier de richting van de normaalvectoren.

Deze manier van werken is echter minder bestand tegen ruis, die toch wel aanwezig is in opge-meten puntenwolken. Ook wordt er detail verloren aan de scherpe randen van objecten. Deze problemen kunnen gedeeltelijk opgelost worden te werken met gewogen waarden bij de least squares methode en de range van punten die bekeken worden te bewerken. Om een heel gede-tailleerde schatting te verkrijgen van de normaalvectoren zullen de nieuwe voorgestelde metho-des in deze paper wellicht handig zijn. Maar om te functioneren in een clusteringomgeving en met de doeleinden die hier nagestreefd worden, is een basisschatting van de normaalvectoren voldoende, zoals de testen later zullen bewijzen.

Deze least squares methode aan de hand van eigen-analyse wordt ook beschreven in Mitra & Nguyen (2003) en biedt voldoende accuraatheid om verdere clustertechnieken op toe te passen.

Hoofdstuk 2. State of the art 17

2.2 Datastructuur

Er wordt in dit project gewerkt met een grote hoeveelheid aan data. Er zijn dus optimalisaties nodig die het geheel gaan versnellen, anders lijkt een uitkomst die zowel snel als correct is onmogelijk. De beste manier om zo'n optimalisatie te implementeren, is het creëren van een datastructuur die toelaat om het zoeken in de data heel wat te versnellen. In de ruwe data is de voorstelling immers helemaal niet optimaal. Dit is in feite gewoon een lijst van punten. Belangrijke zaken, zoals het bepalen van afstanden, zoeken van naburige punten,.. is een heel intensief proces in een dataset die niet geordend is en verder geen geoptimaliseerde structuur heeft.

2.2.1 Octree

De eerste datastructuur die gevonden werd, die het geheel zou kunnen optimaliseren is de `octree'. Dit is een datastructuur die al een hele tijd in de omloop is en wordt reeds gebruikt om allerhande 3D applicaties te versnellen in executietijd. Maegher (1981) bespreekt de werking van een octree en past dit in verschillende scenario's toe waarin het zijn nut kan bewijzen.

Figuur 2.7: Principe octree Bron: Paul (2015)

Figuur 2.7 stelt deze datastructuur visueel voor. Men zal de 3D ruimte dus indelen in 8 cellen en daarna kunnen al deze cellen elk nogmaals ingedeeld worden in 8 cellen en dit kan zo blijven doorgaan. Men kan in principe kiezen hoeveel lagen er zo aanwezig zijn in deze octree, maar op een bepaald moment heeft verder opdelen weinig zin. Dit is iets dat moet afgetast worden met hoe dicht de data nog is.

Hoofdstuk 2. State of the art 18 Ook in de omgeving van puntenwolken worden octrees toegepast Schnabel & Klein (2006). Het is zelfs reeds gebruikt om speciek aan segmentatie te doen Vo et al. (2015) en Su et al. (2016).

Het is duidelijk dat het opzoeken van punten heel wat sneller zal gaan. Er moet in dit geval niet meer geïtereerd worden over de hele lijst om bijvoorbeeld naburige punten te vinden. In deze implementatie hoeft men enkel dezelfde cel en naburige cellen te controleren. En omwille van de geordende datastructuur zijn de naburige cellen ook zeer snel te verkrijgen.

De twee problemen die zich voordoen in een puntenwolk zijn `Wat zijn de x aantal dichtst-bijzijnde punten?' en `Welke punten liggen op een maximale afstand y van punt z?'. Voor beide van deze problemen zorgt een octree ervoor dat er veel minder punten moeten overlopen worden.

2.2.2 K-d tree

Een andere datastructuur die ook prominent aanwezig is in de 3D wereld om allerhande algoritmes te optimaliseren, is de k-d tree Greenspan & Yurick (2003). K-d tree is afgekort voor `k-dimensiontree'. En deze naam is zo gekozen omdat de k hier variabel is. Er wordt hier, net zoals bij de octree, opgesplitst in subcellen. Maar daar waar de cellen in een octree steeds in 4 even grote cellen opgesplitst worden, wordt er in een k-d tree slechts in 2 cellen opgesplitst, en deze cellen hoeven niet per se even groot te zijn.

Figuur 2.8: Principe k-d tree Bron: Gowers (2020)

Figuur 2.8 geeft weer hoe een k-d tree zijn opdeling doet. De eerste opdeling gebeurt volgens de rode lijn. Merk op dat, aangezien er slechts er in 2 cellen wordt opgedeeld, er steeds volgens 1 enkele as wordt opgedeeld. Vervolgens gebeurt er splitsing volgens de groene vlakken, dit zijn er 2, want er zijn reeds 2 cellen. De derde splitsing is gemaakt door de vier blauwe vlakken. Dit stelt hier een driedimensionale k-d tree voor die 3 keer gesplitst is.

Hoofdstuk 2. State of the art 19 Na 3 keer opsplitsen zijn er bij een k-d tree evenveel cellen als bij een octree. Maar de verdeling van deze cellen zit wel anders in elkaar. De splitsing gebeurt aan de hand van hyperplanes, punten aan de ene kant van de hyperplane behoren tot de ene cel, punten aan de andere kant behoren tot de andere cel. De locatie van deze hyperplanes kunnen op verschillende manieren bepaald worden. Een manier om dit bijvoorbeeld te doen is om de mediaan te nemen van alle punten die in de bovenliggende cel aanwezig zijn. Men kan in principe ook gewoon het midden nemen of een andere metriek.

Een k-d tree is dus eigenlijk een binaire boom, maar dan toegepast op meerdere dimensies. Een k-d tree kan ook bestaan in meer dan 3 dimensies.

Van de k-d tree zijn ook bepaalde uitvoeringstijden voor algoritmes bekend, dit wordt in de `big-O-notation' weergegeven in tabel 2.1.

Tabel 2.1: Uitvoeringstijden van algoritmes bij k-d trees.

Algoritme Best case Worst case Space O(n) O(n) Search O(log(n)) O(n) Insert O(log(n)) O(n) Delete O(log(n)) O(n)

Hoofdstuk 2. State of the art 20

2.3 Segmentatie

2.3.1 Histogram aanpak

Het bepalen van vlakken alleen is onvoldoende om een BIM te bekomen, hier moeten nog muren en vloeren uit gehaald worden. Hiervan zijn verschillende aanpakken mogelijk, zoals de aanpak in Macher et al. (2015). Daarbij worden eerste de vloeren gesegmenteerd door te kijken naar de verdeling van de punten in de z-richting, wat hier als verticale richting wordt aanzien. Daarna kan in slices gekeken worden naar de andere punten. In deze techniek wordt eigenlijk geheel in 2D gewerkt in plaats van in 3D om de complexiteit die geïntroduceerd door deze derde dimensie achterwege te laten en zo een sneller werkend algoritme te bekomen.

Figuur 2.9: Vloersegmentatie aan de hand van puntenverdeling. De pieken stellen de verdiepingen voor. `desk' kan hier mogelijks fout bepaald worden.

Bron: Macher et al. (2015)

Figuur 2.9 geeft een histogram weer die de verdeling in de verticale richting toont. Er zijn duidelijke pieken zichtbaar die horizontale vlakken voorstellen. De grote vlakken stellen dan natuurlijk de vloeren voor.

Een probleem met deze techniek is dat ervan uitgegaan wordt dat de scan gealigneerd staat met de assen in het Cartesiaans stelsel. Of toch ten minste dat de verticale as gealigneerd is met de verticale richting van het gebouw. Dit is zo goed als nooit het geval. Scans worden soms bewust niet rechtopstaand genomen omdat dit een betere verdeling van punten kan opleveren.

Hoofdstuk 2. State of the art 21 Maar belangrijker nog is dat bij het stitching proces, waarbij dus meerdere scans bij elkaar gevoegd worden, er ook vaak voor gezorgd wordt dat het geheel wat gedraaid wordt. Deze rotatie kan waarschijnlijk wel bepaald worden, maar dit zal opnieuw neerkomen op het bepalen van normaalvectoren. En wanneer deze dan eenmaal bepaald zijn, zijn er andere manieren die wellicht makkelijker zullen werken. Of het kan natuurlijk ook zijn dat de oriëntatie van de scanner bekend is, in welk geval het rechtopzetten van de data heel eenvoudig wordt. Wat wel een interessant gegeven is aan deze methode, is het feit dat er veel in 2D gewerkt wordt, op vele plaatsen kunnen projecties gemaakt worden die inderdaad het werk veel minder complex maken en algemeen sneller laten verlopen.

Macher et al. (2017) Volgt een gelijkaardige aanpak. Daar wordt ook eerst segmentatie per verdieping gedaan volgens het histogram principe. Ook wordt er gebruik gemaakt van een 2D-projectie waarbij er gewerkt wordt met een binaire afbeelding. Deze paper werkt verder op de bereikte resultaten van Macher et al. (2015).

Figuur 2.10: Binaire afbeelding van een verdieping. Wit duidt op de aanwezigheid van informatie (punten) op die locatie.

Bron: Macher et al. (2017)

Figuur 2.10 geeft een voorbeeld weer van zo'n binaire afbeelding. Hierin kan de splitsing van verschillende kamers makkelijk afgelezen worden, elke kamer wordt immers volledig door lege cellen afgesloten. Dit gaat er dan natuuurlijk ook vanuit dat alle deuren tijdens het scannen gesloten waren om zo een duidelijke afsluiting te hebben tussen de kamers. Dit is echter een moeilijker gegeven om algemeen bij elke dataset te doen werken. Er zijn andere manieren om nog steeds de kamers hier goed uit af te leiden, maar deze zullen iets ingewikkelder zijn dan de eenvoud die hier aanwezig is.

Hoofdstuk 2. State of the art 22

2.3.2 Muren

Na het verdelen van de punten in verschillende verdiepingen en kamers, wordt er vanaf dit punt verder gezocht naar vlakken die de muren kunnen voorstellen. Dit kan op de manieren die reeds in voorgaande secties zijn voorgesteld. De verdere stap hierin is het eectief extraheren van muren uit deze resultaten. Macher et al. (2017) stelt hierbij voor om bij de indoorscan de verschillende vlakken te combineren op basis van hun afstand van elkaar en of deze al dan niet parallel zijn aan elkaar. De afstand tussen twee vlakken kan hierbij op meerdere manieren bepaald worden. Men kan de dichtst mogelijke afstand tussen de twee vlakken nemen, de verst mogelijke afstand, of een gemiddelde. Daarbij dienen ook nog verschillende drempelwaarden gekozen te worden. Een daarvan is hoeveel de rotatie van de twee vlakken van elkaar mogen verschillen om deze als parallel te aanschouwen, de literatuur stelt deze drempelwaarde op 5 graden. Een andere parameter is de afstand die de vlakken hebben ten opzichte van elkaar. De paper in kwestie kiest deze op maximaal een halve meter van elkaar. Daarnaast moet er nog rekening gehouden worden met het minimum aantal punten die een vlak moet bevatten alvorens deze als een muur wordt aanzien. Maar dit zal natuurlijk afhangen van de resolutie van de scanner en welke pre-processing er aan de puntenwolk gebeurd is.

2.3.3 Gladheid

Een compleet andere aanpak die ook werd gevonden in de literatuur om de puntenwolk te segmenteren, is deze die gebruik maakt van gladheid om allerhande objecten in de puntenwolk te segmenteren Rabbani et al. (2006). Deze is vooral gefocust op een industriële setting waar veel verschillende vormen aanwezig zijn en waar vlakken onvoldoende zijn om de omgeving te beschrijven. Dit is dus zeer voordelig in die zin. Het grote nadeel hiervan is dat clutter heel moeilijk weg te krijgen is, aangezien elk soort object gewoon mee gesegmenteerd wordt. Het algoritme werkt op het principe van region growing. Dit principe is wat vergelijkbaar met clusteringalgoritmes die eerder reeds besproken werden. Het gebruikt ook enkel de nor-maalvectoren als extra informatie in de dataset die ook berekend werd op de reeds besproken manier met plane tting. Hoewel gladheid dus waarschijnlijk zelf minder handig lijkt om hier concreet toe te passen, het gebruikt toch ook region growing, wat al een paar keer voorgekomen is en zijn nut binnen dit domein al bewezen heeft.

Hoofdstuk 3

Aanpak

Na het verzamelen van de nodige informatie om dit project te kunnen starten, werden de meest veelbelovende zaken geïmplementeerd. Dit werd gedaan volgens een iteratief proces dat elke stap evalueert om vervolgens een verbetering te implementeren of over te gaan naar de volgende stap. Concreet werd de eerste stap aanzien als de voorbereiding. Hierbij werd de data goed bestudeerd en voorbereid om de volgende algoritmes op te kunnen toepassen. De volgende stap bestaat uit het omzetten van de data naar een meer bruikbare vorm. De laatste stap bestaat er dan in om de eectieve geometrische vormen te segmenteren. Daarna werden er nog verschillende variaties op de data geïmplementeerd om een evaluatie te kunnen uitvoeren. Deze variaties houden bepaalde vormen van ruis in.

3.1 Voorbereiding

Het voorbereiden van de ontwikkeling begon met het kiezen van het beste platform om in te werken. Al snel werd de keuze van programmeertaal gereduceerd tot C++ of python. In C++ was er sprake van een veelgebruikte bibliotheek genaamd `PCL' ofwel `Point Cloud Library'. Deze bevat veel ingebouwde functionaliteiten voor puntenwolken. Het nadeel hiervan echter is het veel trager zijn in het opzetten van projecten. C++ is ook niet de beste programmeer-taal om snelle ontwikkelingen in voort te zetten, wel is het zeer goed om eciënte code te ontwikkelen.

Om de snelheid van de progressie te bevorderen, werd uiteindelijk gekozen voor python, er is immers een binding gemaakt in python voor PCL. Deze bibliotheek zou dus ook beschikbaar zijn in een veel eenvoudigere omgeving. Later werd dan ook nog overgestapt op open3D als bibliotheek, die ook heel veel functionaliteiten voorzien heeft voor puntenwolken en heel wat beter ondersteund is. Voor PCL is immers niet alle inhoud die in C++ aanwezig is, ook aanwezig in python.

Hoofdstuk 3. Aanpak 24

3.1.1 Exploratie

Er zijn vele verschillende datasets aangeleverd om de algoritmes op te testen. Figuur 3.1 illustreert een voorbeeld van zo'n puntenwolk. Dit is de set waarop gedurende de rest van het project het meest op gewerkt werd, aangezien deze de meest gedenieerde vormen heeft. Muren en vloeren zijn hier heel herkenbaar op deze dataset. Er zijn ook datasets waarop dit heel wat minder duidelijk is. Figuur 3.2 toont bijvoorbeeld een scan die eerder ruizig is en waar niet echt herkenbare vormen in aanwezig zijn. Dit is mogelijks ook een scan die buiten is genomen, en dus niet de binnenkant van een gebouw heeft gescand, wat voor een BIM model eigenlijk niet zal werken.

Figuur 3.1: Puntenwolk van een gang met zijkamers.

Hoofdstuk 3. Aanpak 25 Wanneer men erder gaat op de duidelijkere dataset, zijn hier een paar eigenschappen die werden bekeken alvorens hier enig algoritme op werd toegepast. Dit is wellicht een van de kleinere datasets, maar bevat nog steeds 4.215.405 punten. Dit is toch een hele hoop data die verwerkt moet worden. Er zal dus zeker eciëntie in het achterhoofd gehouden moeten worden, bij zoveel data is het immers makkelijk om een zeer tijdsintensieve berekening te starten.

Een belangrijk gegeven is ook dat er artefacten aanwezig zijn. LIDAR zit immers met een probleem dat de lasers kunnen reecteren in vensters. Dit is dan goed zichtbaar op Figuur 3.1. Aan de rechterzijde van de afbeelding kunnen er punten opgemerkt worden die eigenlijk buiten het gebouw vallen. Het ziet er echter wel naar uit dat deze relatief goed gekarakteriseerd kunnen worden door een lage densiteit in punten en ook door het feit dat de vormen niet echt goed herkenbaar zijn. Figuur 3.3 geeft een beter beeld van hoe deze artefacten eruit zien. Daarop is goed te zien dat de dichtheid van de punten iets minder is dan voor de punten die een correcte meting voorstellen. Met een algoritme dat omkan met ruis zou dit dus oplosbaar moeten zijn.

Figuur 3.3: Artefacten ten gevolge van reecties zichtbaar aan de rechterkant.

Bij deze puntenwolk wordt ervan uitgegaan dat punten een maximale afwijking van 2 tot 3 centimeter hebben van het object waarop deze geprojecteerd zijn. Dit is zo gekozen omdat dit het worst case scenario voorstelt voor standaard LIDAR scanners.

Het laatste gegeven dat werd bekeken, is de oriëntatie van de puntenwolk. In de literatuur werden immers al aannames gedaan dat de scan steeds rechtopstaand gebeurt. Figuur 3.4 toont echter dat dit niet het geval is voor deze dataset. Ofwel is er dus een techniek nodig die deze rotatie gaat berekenen, ofwel is er een ander algoritme vereist om de vlakken te kunnen segmenteren.

Hoofdstuk 3. Aanpak 26

Figuur 3.4: Oriëntatie van de puntenwolk.

3.1.2 Reductie

Zoals eerder vermeld, bevat een van de kleinere datasets reeds 4.215.405 punten. Het ene algoritme kan al beter overweg met grote hoeveelheden data dan andere. RANSAC is typisch een algoritme dat minder last heeft van het feit dat een dataset heel groot wordt, daar is vooral de metriek die de fout berekent bepalend voor de snelheid van het algoritme, maar deze is al-gemeen O(n). Ook Hough-transformatie zal hier lineair mee stijgen. Van clusteringalgoritmes hangt dit af van algoritme tot algoritme. Deze zijn typisch wel optimaliseerbaar met behulp van datastructuren zoals bijvoorbeeld de k-d tree.

O(n) is niet slecht, maar dit loopt toch al op voor een dataset van die grootte, daarom dat geopteerd werd voor het verkleinen van deze dataset. Dit vooral om de uitvoeringstijden dus te kunnen verlagen zodanig dat sneller tests uitgevoerd kunnen worden. Om dit te doen, werd gestoten op de term `voxel downscaling', wat in het Nederlands `voxelreductie' wordt genoemd.

Voxelreductie is een algoritme dat de hele dataset opdeelt in een 3D rooster. Alle punten aanwezig in eenzelfde cel worden gecombineerd zodanig dat er maar 1 punt meer aanwezig is in de cel. Dit combineren gebeurt door het gemiddelde te nemen van de locatie van deze punten. Algoritme 4 beschrijft wat er gebeurt bij voxelreductie.

Hoofdstuk 3. Aanpak 27 Algoritme 4: Voxelreductie

cells(xsize / cellsize, ysize / cellsize, zsize / cellsize); for point in point-list do

cell = ndAppropriateCell(point) ; // Discretiseren van de coördinaten en opzoeken in het rooster

cell.addPoint(point); end

newpoints = [ ]; for cell in cells do

cellsum = [0, 0, 0]; for point in cell do

cellsum = cellsum + point; end

newpoint = cellsum / length(cell); newpoints.append(newpoint); end

Hoofdstuk 3. Aanpak 28 Hier is dus eigenlijk maar 1 parameter belangrijk, en dat is de grootte van de cellen, die ook wel voxel genoemd worden. Hoe groter deze voxels genomen worden, des te meer detail verloren gaat uit de dataset. Dit process kan eigenlijk wat vergeleken worden met een blurlter in beeldverwerking, waarbij ook een gemiddelde genomen wordt. Een voxel is hier dan analoog aan een pixel.

Ruis

Dit downscalen heeft dus het voordeel dat de data gereduceerd wordt. Daardoor worden de volgende algoritmes sneller. Maar daarnaast zijn er nog twee extra voordelen die men kan toekennen aan voxelreductie. Ten eerste zal dit ruis onderdrukken. Dit is dus opnieuw de analogie van een blurlter in beeldverwerking. Daarbij zullen ook ruizige elementen onder-drukt worden door het gemiddelde te nemen van verschillende pixels. De nauwkeurigheid van standaard LIDAR-sensoren ligt ongeveer op 2 cm. Bij minder accurate sensoren kan dit soms oplopen tot 3 cm. Als men dus voxels met een grootte van 6 cm zou nemen, dan gaat deze de ruis normaalgezien optimaal onderdrukken. Figuur 3.5 geeft weer wat het eect is van reduceren op ruis.

Figuur 3.5: Ruisonderdrukking ten gevolge van voxelreductie, rood zijn de opgemeten punten, blauw is de grondwaarheid en groen is het uitgemiddelde resultaat.

Merk op dat bovenstaande guur een ideaal geval weergeeft. Hoe minder een voxel gecentreerd staat op de muur, des te minder goed het resultaat zal zijn. Ondanks dat zal het resultaat wel altijd een punt zijn die dichter op de muur komt te liggen. Ruis wordt dus steeds wel een beetje onderdrukt.

Hoofdstuk 3. Aanpak 29 Uniforme data

Een tweede eect ten gevolge van de opmeting, en ook eventueel ten gevolge van het stitchen, is de oneen verdeling van de punten in de puntenwolk. Een LIDAR-sensor zal met een bepaalde resolutie 360 graden om zich heen gaan scannen. Dit heeft als eect dat objecten dichtbij de scanner veel meer scanpunten toegekend krijgen dan objecten veraf. Dit fenomeen wordt geïllustreerd in Figuur 3.6. Daarin is duidelijk te zien dat de kleinere cirkels die dichtbij de scanner staan dichter opeen staan dan de cirkels verder weg van de scanner.

Figuur 3.6: Dichtheid van de scanpunten. Hoe verder de scanpunten verwijderd zijn van de scanner, des te lager de densiteit van de puntenwolk wordt.

Bron: Dewan et al. (2017)

Deze oneen verdeling van de punten zou in verdere verwerking kunnen zorgen voor problemen. Voxelreductie is ook voor dit probleem een oplossing. Er kunnen immers niet meer punten aanwezig zijn in de dataset dan dat er voxels zijn. Dus ook waar de puntenwolk een heel grote densiteit heeft, zullen er na de voxelreductie evenveel punten aanwezig zijn als op andere objecten.

De optimale parameter in deze zin kan gekozen worden door na te denken waar de puntenwolk de kleinste dichtheid heeft, maar nog steeds informatie bevat die men zeker wil behouden. Dat is dan ook de dichtheid die men in de gehele puntenwolk zou willen hebben.

Hoofdstuk 3. Aanpak 30 Uitvoering

Voor de uitvoering werd uiteindelijk gewerkt met een voxelbreedte van 5 cm. Dit stemt ongeveer overeen met wat aangeraden wordt vanuit het standpunt van ruisonderdrukking. Maar wat nog belangrijker bleek te zijn dan dat argument, was dat een voxelgrootte van 5 cm het meest ideaal bleek te zijn om te werken met de algoritmes die nog volgen. De snelheid op dit niveau is acceptabel, maar de algoritmes werken ook nog steeds op voldoende wijze. Als men de cellen groter zou nemen, dan zijn er accuraatheidproblemen en als men de cellen kleiner zou nemen, dan loopt de uitvoeringstijd enorm op. Visueel is er alvast slechts een heel miniem verschil merkbaar. Dit kan men zien op Figuur 3.7. Men kan daar zien dat de puntenwolk wel degelijk wat minder dicht is, maar de vormen zijn nog steeds zeer duidelijk herkenbaar.

Figuur 3.7: Visueel eect van reductie, links is de originele data, rechts is de data na reductie.

Hoofdstuk 3. Aanpak 31

Tabel 3.1: Dataverlies na voxelreductie

Dataset Voor reductie Na reductie Percentage appartement 29885654 126746 0.424 slaapkamer 16888124 83628 0.495 bestuurskamer 25292488 122423 0.484 gang 1 4215405 321273 7.621 gang 2 1622607 412820 25.44 lobby 33250267 218504 0.657 loft 22572887 146669 0.650 gang 3 2703862 368181 13.62 dow 3063059 2679014 87.46 Hasselt 19615568 16536192 84.30 Schoonmeersen 2033197 1504116 73.98 world ufo 2209719 2022789 91.54 world ufo full 10750846 4129130 38.41

In tabel 3.1 is het opvallend hoeveel variatie er zit in het aandeel punten die verwijderd worden uit de dataset. Bij sommige datasets blijft er amper een halve procent aan data over, terwijl er bij andere meer dan 80% behouden wordt. Dit is vooral te wijten aan het feit dat de datasets van verschillende afkomsten zijn. De ene reeks datasets waren reeds wat opgeschoond, terwijl dit bij de andere datasets niet het geval is. Daarom worden er ook veel minder punten uitgelterd voor sommige sets. De datasets die nog niet opgeschoond zijn, zijn afkomstig van een externe bron. Ook visueel is dit goed op te merken, vaak beschrijven de bestanden van externe afkomst een veel kleinere ruimte, maar toch is het aantal punten in de dataset tot 10 keer groter dan in de pcd bestanden. Hoe deze sets eruit zien, wordt verder besproken bij de evaluatie.

Dit bewijst ook het nut van het implementeren van dergelijke opschoning. Zolang de dichtheid van de gebruikte dataset groter is dan 1 punt per 125 cm³ (5 cm · 5 cm · 5 cm), dan zullen de volgende algoritmes volgens verwachting blijven werken. Dit zorgt er eigenlijk voor dat de rest van het geheel generiek behandeld kan worden.

3.1.3 Synthetische data

Problemen bij echte data

Bij de data exploratie werd tot de conclusie gekomen dat deze vaak onvolledig is. Dit dan in de zin van muren die niet ofwel heel weinig gescand zijn. Er waren vooral ook veel muren die slechts aan 1 kant gescand zijn zodanig dat er geen dikte kan bepaald worden van de muur.

Hoofdstuk 3. Aanpak 32

Figuur 3.8: Een onvolledige scan, de rode lijn stelt de contouren voor die de muur waarschijnlijk zou moeten hebben.

Op Figuur 3.8 kan een voorbeeld aanschouwd worden waarbij de scan niet volledig is. Dit stelt een muur voor uit de aangeleverde dataset en visueel zag deze scan er het meest volledig uit. Dus ook in de meest volledige scans zijn er inaccuraatheden op de muren, vloeren en plafonds. Dit is iets waar in de toekomst rekening mee gehouden zal moeten worden.

Ook is er een heel deel `clutter' aanwezig in de dataset. Clutter wijst eigenlijk op alles in de dataset dat geen gebouwstructuren voorstellen. Dit kunnen alle soorten objecten zijn, stoelen, tafels, banken,... Deze clutter kan er potentieel voor zorgen dat de algoritmes niet goed werken en is het dus interessant om deze te laten werken op data die deze niet bevat.

Hoofdstuk 3. Aanpak 33 Figuur 3.9 toont twee voorbeelden van deze clutter. Er zijn hierbij ook verschillende vormen van clutter op te merken. Op de Figuur is het linkervoorbeeld een verzameling van punten die niet echt een gedenieerde vorm heeft. Dit is dan ook relatief makkelijk te herkennen en onderscheiden van een muur. Het rechtervoorbeeld is echter anders, dat is waarschijnlijk een stuk meubilair dat ook gewoon uit vlakken bestaat. Dit kan zeker problemen vormen aangezien vele van de andere algoritmes steunen op het detecteren of gebruiken van vlakken. Dit stuk clutter kan dus zeker een nadelig eect hebben op het resultaat.

Creatie

Dit zijn problemen die men zeker niet zomaar kan achterwege laten en waar later ook zeker mee zal gehandeld moeten worden. Maar het is wel interessant om ook te kijken naar een ideaal geval. Daar waar alle vlakken volledig opgemeten zijn en waar er geen clutter aanwezig is. Dit type data is wellicht niet te vinden in reële scans. Daarom werd er geopteerd om te werken met synthetische data.

Een tweede zaak die dit geheel ook interessant maakt, is het feit dat men van synthetische data met 100% zekerheid kan zeggen wat de dimensies zijn van de verschillende gebouwstructuren. Er is dus sowieso een grondwaarheid aanwezig over de data waarmee de resultaten vergeleken kunnen worden. Dit is niet het geval voor de aangeleverde data. Daarvoor waren er geen dimensies van het gebouw beschikbaar, dus moet de evaluatie ook op een andere manier gebeuren.

De syntethische data werd op zo'n danige manier opgebouwd dat deze zoveel mogelijk gelijkt op hoe echte data verkregen wordt. Daarmee wordt bedoeld dat de muren eectief als 2 parallele vlakken werden voorgesteld en dat een virtuele LIDAR scanner hier dan zijn projecties op maakt. Het resultaat hiervan is te zien in Figuur 3.10.

Hoofdstuk 3. Aanpak 34

Figuur 3.10: grondwaarheid van de synthetische dataset, vloer en plafond zijn hier niet zichtbaar om de binnenste muren goed te kunnen zien.

Omdat men hier uitgaat van ideale data, worden de kamers ook eenvoudig voorgesteld. Deze zijn allemaal balkvormig. Dit heeft als gevolg dat alle muren, vloeren en plafonds perfect rechthoekig zijn. De virtuele scanner wordt steeds in exact het midden van de kamer ge-plaatst. Deze kan dan met een instelbare resolutie 360 graden in horizontale richting en 180 graden in verticale richting de scans uitvoeren. Het algoritme is ook zo opgesteld dat er ook nog steeds ruis aanwezig kan zijn in de dataset. Op die manier is de data niet zodanig on-realistisch dat alles toegepast hierop niet meer zou werken op echte data. Ook kan het voor extra interessante informatie zorgen wanneer de robuustheid van de algoritmes getest moeten worden. De pseudocode voor deze scan is te vinden in algoritme 5.

Hoofdstuk 3. Aanpak 35 Algoritme 5: Scannen van synthetische data

Function ScanRoom(planes, center): list-points = [ ];

for h < horizontal-resolution do

horizontal-angle = 360 / horizontal-resolution * h; for v < vertical-resolution do

vertical-angle = 180 / vertical-resolution * v - 90;

direction = calculateRayDirection(horizontal-angle, vertical-angle) ; // richting in het Cartesiaans stelsel berekenen op basis van de horizontale en de verticale hoek

point = getClosestIntersection(planes, center, direction) ; // Snijpunten met alle vlakken berekenen en het dichtste punt nemen als resultaat list-points.append(point); v = v + 1; end h = h + 1; end return list-points;

De resolutie van de virtuele scanner werd heel erg hoog gezet, aangezien de stap van voxelre-ductie toch aanwezig is om de data te reduceren tot een werkbare hoeveelheid. Op die manier is ook zeker alle nodige data opgemeten. Op Figuur 3.11 staat het resultaat van deze scan weergegeven.

Hoofdstuk 3. Aanpak 36 Merk op dat hier op deze guren ook ramen en deuren aanwezig zijn. Deze data is dus opgebouwd in de veronderstelling dat deuren openstaan en ramen niet detecteerbaar zijn door de sensor. Dit is echter een aanname die niet overal kan gedaan worden. Deuren staan meestal open op de gebruikte datasets, maar wanneer deze gesloten zijn, is het heel dikwijls onmogelijk het onderscheid tussen deur en muur visueel te maken. Ook voor ramen geldt een gelijkaardige aanname. Wanneer men een sensor gebruikt gebaseerd op het TOF principe, zal glas niet als onzichtbaar aanschouwd worden. Wat wel gebeurt, is dat deze gewoon in de muur verwerkt worden waardoor onderscheid maken ook weer zo goed als onmogelijk wordt, of er zijn reecties zichtbaar, wat eigenlijk foute data is. In de synthetische dataset werd zo gewerkt omwille van de eenvoud, en ook om te kunnen experimenteren met algoritmes die gaten in puntenwolken kunnen herkennen.

Op Figuur 3.11 is de puntenwolk nog heel erg dicht, daar wordt eerst nog voxelreductie op toegepast om deze bruikbaar te maken, het resultaat hiervan is te zien op Figuur 3.12.

Figuur 3.12: Gegenereerde dataset na voxelreductie.

De dikte van de voorstelling van de punten is in Figuur 3.12 ook wat aangepast om het geheel wat beter zichtbaar te maken. Ondanks het feit dat hier aan voxelreductie werd gedaan, zijn er nog steeds cirkelvormige guren zichtbaar, vooral op de vloer en het plafond. Dit betekent echter niet dat de punten daar dichter opeen gepakt zitten, dit is gewoon een artefact dat nog achterblijft door de manier van scannen.