University of Southern California Dissertations and Theses

3D urban modeling from city-scale aerial LiDAR data

(USC Thesis Other)

3D urban modeling from city-scale aerial LiDAR data

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content 3DURBANMODELINGFROMCITY-SCALEAERIALLIDARDATA by Qian-YiZhou ADissertationPresentedtothe FACULTYOFTHEUSCGRADUATESCHOOL UNIVERSITYOFSOUTHERNCALIFORNIA InPartialFulfillmentofthe RequirementsfortheDegree DOCTOROFPHILOSOPHY (COMPUTERSCIENCE) August2012 Copyright 2012 Qian-YiZhou Acknowledgments I would like to express my sincere thanks to my advisor Prof. Ulrich Neumann. His insightfulguidance,inspiringadvice,andscientificknowledgeisaninvaluabletreasure for me during my five years’ study in the University of Southern California. I had a great time at the Computer Graphics and Immersive Technologies lab, and I hope our cooperationwillcontinueinthefuture. I would like to thank Prof. Suya You, Prof. C.-C. Jay Kuo, Prof. Jernej Barbiˇ c, and Prof. Aiichiro Nakano for taking their precious time to serve on the committee of my qualifying exam and thesis defense. I am very grateful to the people I have worked with or shared knowledge with. Special thanks to Prof. Shi-Min Hu from Tsinghua University, Prof. Tao Ju from Washington University in St. Louis, Thommen Korah fromNokiaResearchCenter,FlorentLafargefromINRIA,andthecurrentandprevious CGITlabmembers: RongqiQiu,KelvinChung,LucianoNocera,TaehyunRhee,Char- alambos ’Charis’ Poullis, Jonathan Mooser, Seon-Min Rhee, Lu Wang, Quan Wang, SangYunLee,TanasaiSucontphunt,ZhenzhenGao,WeiGuan,GuanPang,YiLi,Jing Huang, Jiaping Zhao, and Yili Zhao. I acknowledge Airborne 1 Corp., Sanborn Corp., and Chevron Corp. for providing data, and USC Provost’s PhD fellowship, Lockheed MartinCorp.,andNokiaCorp. forfinancialsupport. ii IamgratefultothesupportfrommyparentsinChina. IthankallmyfriendsatUSC for making my life here a happy and precious experience: Jing Zhao, Fan Sun, Chang Huang,RuiWang,XiaoxunSun,andmanyothers. Finally, a very special thanks to my fiancee Na Chen, a woman who is excellent in all aspects, beautiful, smart, and generous. I am very grateful for her encouragement, herpatience,andallherloveandtenderness. iii TableofContents Acknowledgments ii ListofTables vii ListofFigures ix Abstract xvii Chapter1: Introduction 1 Chapter2: RelatedWork 7 Chapter3: GeneralUrbanModelingApproach 10 3.1 UrbanModelingPipeline . . . . . . . . . . . . . . . . . . . . . . . 10 3.2 VegetationDetection . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2.1 Features . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2.2 Classifier . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.2.3 ConfigurationandResults . . . . . . . . . . . . . . . . . . 18 3.3 Segmentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.3.1 RegionGrowingImplementation . . . . . . . . . . . . . . . 20 3.3.2 AgglomerativeClusteringAlgorithm . . . . . . . . . . . . . 20 3.4 BuildingModeling . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.4.1 PlanarRoofPatchExtraction . . . . . . . . . . . . . . . . . 22 3.4.2 BoundaryDetection . . . . . . . . . . . . . . . . . . . . . . 23 3.4.3 BuildingModelCreation . . . . . . . . . . . . . . . . . . . 25 3.4.4 Non-PlanarRoofExtension . . . . . . . . . . . . . . . . . . 29 3.5 TerrainModeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.6 ExperimentalResultsonSmallSampleDataSets . . . . . . . . . . 31 Chapter4: StreamingUrbanModelingforLarge-ScaleDataSets 35 4.1 ReviewofStreamingApproaches . . . . . . . . . . . . . . . . . . . 36 4.2 StreamingFramework . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2.1 PointStreamsandFinalizer . . . . . . . . . . . . . . . . . . 37 iv 4.2.2 StreamingOperatorsandStates . . . . . . . . . . . . . . . . 39 4.2.3 StatePropagation . . . . . . . . . . . . . . . . . . . . . . . 41 4.3 StreamingUrbanModeling . . . . . . . . . . . . . . . . . . . . . . 43 4.3.1 StreamingModelingPipeline . . . . . . . . . . . . . . . . . 43 4.3.2 StreamingClassification . . . . . . . . . . . . . . . . . . . 45 4.3.3 StreamingSegmentation . . . . . . . . . . . . . . . . . . . 46 4.3.4 BuildingModeling . . . . . . . . . . . . . . . . . . . . . . 49 4.3.5 TerrainModeling . . . . . . . . . . . . . . . . . . . . . . . 50 4.4 ExperimentalResultsonCity-ScaleDataSets . . . . . . . . . . . . 51 4.4.1 StreamingVersusTiling . . . . . . . . . . . . . . . . . . . 54 Chapter5: 2.5DDualContouring 57 5.1 BriefReviewofDualContouring . . . . . . . . . . . . . . . . . . . 59 5.2 2.5DDualContouringPipeline . . . . . . . . . . . . . . . . . . . . 60 5.3 2.5DScanConversion . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.3.1 SurfaceHermiteData . . . . . . . . . . . . . . . . . . . . . 62 5.3.2 BoundaryHermiteData . . . . . . . . . . . . . . . . . . . . 63 5.4 AdaptiveCreationofGeometry . . . . . . . . . . . . . . . . . . . . 64 5.4.1 2.5DQuadraticErrorFunction . . . . . . . . . . . . . . . . 65 5.4.2 QuadtreeSimplificationwithQEFs . . . . . . . . . . . . . . 67 5.5 PolygonGeneration . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.5.1 SharpFeaturePreservingTriangulation . . . . . . . . . . . 69 5.6 Topology-SafeSimplification . . . . . . . . . . . . . . . . . . . . . 71 5.7 PrincipalDirectionSnapping . . . . . . . . . . . . . . . . . . . . . 73 5.8 ExperimentalResultsof2.5DDualContouring . . . . . . . . . . . 74 Chapter6: 2.5DBuildingTopology 78 6.1 ReviewofTopologyControlinVolumetricModeling . . . . . . . . 81 6.2 2.5DBuildingTopology . . . . . . . . . . . . . . . . . . . . . . . 81 6.2.1 TopologicalFeatureDefinitions . . . . . . . . . . . . . . . 82 6.2.2 ConnectionsbetweenTopologicalFeatures . . . . . . . . . 84 6.3 ContouringwithTopologyControl . . . . . . . . . . . . . . . . . . 86 6.3.1 Hyper-PointsClustering . . . . . . . . . . . . . . . . . . . 86 6.3.2 HandlingDegenerateCases . . . . . . . . . . . . . . . . . . 89 6.3.3 AdaptiveContouring . . . . . . . . . . . . . . . . . . . . . 90 6.4 ExperimentalResults . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.5 Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Chapter7: 2.5DBuildingModelingbyDiscoveringGlobalRegularities 97 7.1 ReviewofShape-from-SymmetryApproaches . . . . . . . . . . . . 99 7.2 GlobalRegularities . . . . . . . . . . . . . . . . . . . . . . . . . . 101 v 7.2.1 Roof-RoofRegularities . . . . . . . . . . . . . . . . . . . . 102 7.2.1.1 Orientationregularities . . . . . . . . . . . . . . . . . 102 7.2.1.2 Placementregularities . . . . . . . . . . . . . . . . . 103 7.2.2 Roof-BoundaryRegularities . . . . . . . . . . . . . . . . . 104 7.2.3 Boundary-BoundaryRegularities . . . . . . . . . . . . . . . 104 7.3 ModelingwithGlobalRegularities . . . . . . . . . . . . . . . . . . 105 7.3.1 PlanarRoofPatchExtraction . . . . . . . . . . . . . . . . . 107 7.3.1.1 Orientationalignment . . . . . . . . . . . . . . . . . 107 7.3.1.2 Placementalignment . . . . . . . . . . . . . . . . . . 108 7.3.1.3 Coarse-to-fineiteration . . . . . . . . . . . . . . . . . 108 7.3.2 Boundarysegmentproduction . . . . . . . . . . . . . . . . 110 7.3.2.1 Boundarysegmentinitialization . . . . . . . . . . . . 110 7.3.2.2 Segmentheightalignment . . . . . . . . . . . . . . . 110 7.3.2.3 Segmentpositionalignment . . . . . . . . . . . . . . 111 7.3.3 ModelGeneration . . . . . . . . . . . . . . . . . . . . . . . 111 7.4 ExperimentalResults . . . . . . . . . . . . . . . . . . . . . . . . . 112 Chapter8: ModelingResidentialUrbanAreas 116 8.1 RelatedWork . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 8.1.1 TreeDetectioninLiDAR . . . . . . . . . . . . . . . . . . . 118 8.1.2 TreeModeling . . . . . . . . . . . . . . . . . . . . . . . . 119 8.2 PointCloudClassification . . . . . . . . . . . . . . . . . . . . . . . 119 8.3 ModelingofUrbanElements . . . . . . . . . . . . . . . . . . . . . 122 8.3.1 TreeModeling . . . . . . . . . . . . . . . . . . . . . . . . 122 8.3.2 BuildingModeling . . . . . . . . . . . . . . . . . . . . . . 123 8.3.3 GroundModeling . . . . . . . . . . . . . . . . . . . . . . . 124 8.4 ExperimentalResults . . . . . . . . . . . . . . . . . . . . . . . . . 124 Chapter9: ConclusionandFutureWork 128 9.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 9.2 FutureWork . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 Bibliography 132 vi ListofTables 4.1 Statepropagationalgorithm . . . . . . . . . . . . . . . . . . . . . . 42 4.2 Streamingoperatorsandstatesforclassification . . . . . . . . . . . 45 4.3 Streamingunion-findalgorithm( 2 ) . . . . . . . . . . . . . . . . . 48 4.4 Streamingoperatorsandstatesforsegmentation . . . . . . . . . . . 48 4.5 Threedatasetswithdifferentsampleratesaretestedusingthestream- ing building reconstruction system on a consumer-level PC. This tablereportstherunningtimeandmaximummemoryusageineach pipeline module, i.e., Finalizer, Classifier, Splitter, and other com- ponents. The experimental results show the ability of the streaming systemtohandleextremelylargedatasetsinanefficientmanner. . . 52 5.1 Quantitative evaluation of four modeling approaches over models showninFigure5.12 . . . . . . . . . . . . . . . . . . . . . . . . . 77 6.1 Quantitative comparison between modeling with topology control and2.5DdualcontouringusingtheexperimentshowninFigure6.8. Thelastcolumnreportstheratioofcellcollapsesrejectedbytopol- ogy test among all rejected collapses. The topology test becomes dominantin2.5Ddualcontouringwithlargeerrortolerance. . . . . 95 7.1 QuantitativeresultsforthecomparisoninFigure7.6 . . . . . . . . . 112 7.2 Statistics of the experiments in Figure 7.7. Column 3 - 12 show statisticsoftheintermediateresults. Sincethesizeofconstraintsets is considerably smaller than the number of primitives (in bold), the solutionspace(andthusthecomplexity)ofthemodelingproblemis reducedsignificantly. . . . . . . . . . . . . . . . . . . . . . . . . . 114 8.1 Statisticsoftheexperimentsonthreedifferentdatasets . . . . . . . 125 8.2 Executiontimeofeachstepinthemodelingsystem . . . . . . . . . 126 vii ListofFigures 1.1 3DreconstructionforthecityofAtlantafrom683MLiDARpoints, using the streaming urban modeling system presented in Chapter 4. Bottom: fourcloseupsoftheurbanmodel,showing(a)anareawith large flat structures; (b) a downtown area; and (c,d) two residential areas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2 2.5D nature of the building modeling problem. Left: aerial LiDAR has dense sampling on rooftops but very sparse sampling on build- ing walls. Right: desired building models are usually composed of detailedroofstructuresandverticalwalls. . . . . . . . . . . . . . . 5 3.1 Thegeneralurbanmodelingpipelineincludesfourmodules: classi- fication, segmentation, building modeling, and terrain modeling. In particular, the building modeling module takes an individual point patch as input; detects plane patterns; extracts boundaries of these planes; and finally create polygonal building models based on the boundaries. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Illustrationofnormalvariationmeasurements: (a)normalsofpoints around a roof ridge; (b) normal distribution on a Gauss sphere, red/green/bluearrowspointtotheeigenvectorsofC n p withlengthes equal to corresponding eigenvalues; (c) normals of points in the neighborhood of a tree point. Both λ n 1 and λ n 2 are large due to the irregularityofnormals. . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3 Feature values of a small piece of aerial LiDAR point cloud. The top-leftsub-figureisrenderedwithz-valuesasgreyscalesatpoints. Others are rendered with colors representing corresponding feature values of the points. Blue points are more likely to be roof/ground, whileredpointsaremorelikelytobetrees. . . . . . . . . . . . . . 16 viii 3.4 (a) Input LiDAR point cloud with color intensity representing the heightateachpoint. (b)Greenpointsaredetectedastrees. (c)Seg- mentationresult;withtreesandnoises(blackpoints)removedfrom the input, the segmentation algorithm splits building roofs from the groundsuccessfully. . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.5 Agglomerativesegmentationimplementedusingtheunion-findalgo- rithmpresentedin[7] . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.6 An illustration of the boundary detection algorithm. Circles repre- sent the LiDAR points projected on the x-y plane. The boundary of this patch is composed of the boundary points (red circles) con- nectedbyboundarylines(thinredlines). . . . . . . . . . . . . . . . 24 3.7 Boundary extraction for part of an industrial site: (b) without any morphologicaloperation;(c)withanopeningmorphologicalopera- tion,mostartifactsareremoved. . . . . . . . . . . . . . . . . . . . 25 3.8 Buildingmodelingusingprincipaldirectionsnapping: (a)inputLiDAR point cloud chopped from Denver city; (b) histogram of boundary line directions, with four principal directions detected and marked inthefigure;(c)boundariesofplanarroofpatchesaresnappedonto theprincipaldirectionsandsimplebuildingmodelsarecreated. . . . 27 3.9 Reconstruction of an industrial site using non-planar roof exten- sion: (a)inputLiDARdata;(b)detectedobjectpatchesarerendered with bright colors, while ground and noise are rendered with dark- greyandblack;(c)thereconstructedmodels,withdifferenttypesof objectsgeneratedindifferentways. . . . . . . . . . . . . . . . . . . 30 3.10 Automatic urban modeling of a 1km×1km piece from the city of Denver: (a) input LiDAR data; (b,c) modeling result viewed from different perspectives; (d,e,f) closeups of the modeling result with initialpointcloudoverlayed. . . . . . . . . . . . . . . . . . . . . . 32 3.11 Urban modeling of an area in the city of Oakland: (a) input LiDAR pointcloud,withintensityrepresentingpointheight;(b)ahistogram with seven principal directions detected and illustrated as colored arrows in (a); (c) reconstructed building models; both orthogonal cornersandnon-orthogonalcornersarereconstructedcorrectly. . . . 33 4.1 Finalization result of the Oakland data set with cell colors in (b) representing the finalization time. Spatial coherence is revealed by theregularityofcolordistribution. . . . . . . . . . . . . . . . . . . 38 ix 4.2 An example of the state propagation algorithm. The numbers and colors denote states. When the state of a cell is changed (marked with black frame), it notifies its 1-ring neighborhood (red frame) to check if any of them is ready for the next operator. The whole processisarecursiveprocedure. . . . . . . . . . . . . . . . . . . . 42 4.3 An illustration of the streaming building reconstruction pipeline. A pre-processing module (which is called Finalizer in [27]) inserts finalization tags (yellow ovals) and generates a point stream which flows over the Classifier and the Splitter sequentially. Both com- ponents introduce a state-propagation mechanism so that only data with active states (solid colored cells in Classifier and Splitter) are loaded in-core. The Splitter finally generates a point stream and a building stream; which are converted into a terrain model and var- ious building models using the Terrain Generator and the Modeler respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.4 Thestreamingagglomerativeclusteringalgorithmstoresset-treeroots (marked with red frame) in a global hash table. It first performs a union operation on each neighbor point pair illustrated as the dot- ted line in the left figure. Then the algorithm flattens the set-trees shown in the middle figure and push new roots into the hash table (rightfigure). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.5 WiththeprincipaldirectiongridontheOaklanddataset,sixprinci- paldirectionsaredetectedforthebluecell(left),whiletwoprincipal directionsaredetectedfortheorangecell(right). . . . . . . . . . . 49 4.6 Reconstructed urban model of Atlanta city. Closeups of different areasareshownintherightsub-figures. . . . . . . . . . . . . . . . 50 4.7 MemoryusageduringexperimentonAtlantadataset . . . . . . . . 51 4.8 UrbanmodelofDenverwiththreecloseupsshownin(a,b,c). Although principal directions in these areas are different; with the principal direction grid, correct principal directions are generated for each of them. (d)Memoryusageduringprocessing. . . . . . . . . . . . . . 53 4.9 Comparisonbetweenstreamingmethod(left)andtilingmethod(right). Even with padding along tile boundaries and tiles batch processed one-by-one using an in-core modeling program (the blue square shown in the top-right sub-figure), artifacts can be generated along tileboundariesasshowninthebottom. . . . . . . . . . . . . . . . . 54 x 4.10 (a) Reconstructed urban model of Oakland downtown area. (b) Finalization result reveals the spatial coherence between cells; col- ors represent the finalization time. (c,d,e) Three snapshots during the streaming classification algorithm; only a small portion of cells areatactivestates(brightsolidcoloredcells). . . . . . . . . . . . . 56 5.1 Various kinds of building models are created using 2.5D dual con- touring. Fromlefttoright: twostadiummodelswithdifferentkinds ofnon-planarroofs;atypicalflatbuildingfromresidentialarea;and amodernhighbuildingfromurbanarea. . . . . . . . . . . . . . . . 57 5.2 Manually created models in Google 3D warehouse [21], showing the2.5Dnatureofbuildingstructuremodels . . . . . . . . . . . . . 58 5.3 Robust building modeling pipeline: (a) the input point cloud; (b) a 2D grid with surface Hermite data (gold arrows) and boundary Hermitedata(redarrows)attached;(c)hyper-points(turquoiseballs connected by red lines) generated by minimizing QEFs; (d) mesh model reconstructed via 2.5D dual contouring; and (e) final model withboundariessnappedtoprincipaldirections. . . . . . . . . . . . 61 5.4 Generating (a) surface Hermite data samples on grid points: the sample is assigned to the highest roof layer which covers the grid point;(b,c)boundaryHermitedatasamplesongridedges: themax- imum margin line (thin black lines) divides the lower surface Her- mitedatasamplefromthehigherrooflayer. . . . . . . . . . . . . . 63 5.5 2.5Ddualcontouringwithout(left)andwith(right)adaptivegeom- etrysimplification . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.6 (a,b)Creatingsurfacepolygons(coloredhollowpolygons)andbound- arypolygons(coloredsemitransparentpolygons)aroundhyper-point A. Viewing from top, (c) surface polygons are generated at grid points, while (d) boundary polygons are produced for grid edges whichexhibitarooflayergap. . . . . . . . . . . . . . . . . . . . . 68 5.7 Triangulation without (left) and with (right) the sharp feature pre- serving algorithm. The colors of input points represent the squared distancesfromthemesh. . . . . . . . . . . . . . . . . . . . . . . . 70 5.8 Comparisonbetweentopology-unsafesimplification(left)andtopology- safe simplification (right). Undesired features can be created by mergingleafcellsinatopology-unsafemanner. . . . . . . . . . . . 71 xi 5.9 Anunsafesimplificationcasedeniedbythetopologysafetyteststep 3. Since the center grid point has different roof layer assignments in these leaf cells, two different layers in the top-left leaf cell (a) belongtothesamerooflayerinthecoarsecell(b). Unsafemerging mayerasewallfeaturessuchastheoneshownin(c). . . . . . . . . 72 5.10 Roof layer boundaries (thick colored lines) are regularized using principaldirectionsnappingalgorithm. . . . . . . . . . . . . . . . . 73 5.11 Building reconstruction for a 2KM-by-2.5KM urban area of Los Angeles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75 5.12 Building models created using different approaches (from left to right): 2.5D dual contouring, plane-based method proposed in Sec- tion 3.4.3, general mesh simplification over a rasterized DEM, and manualcreation. Pointcolorsdenotethesquareddistancesbetween pointsandgeneratedmodels. Colorbarsunderthemodelsshowthe ratioofpointsatdifferentsquareddistancelevel. . . . . . . . . . . 76 5.13 Models of similar quality are generated with the same point cloud embeddedintogridsofdifferentsizesordifferentorientations. . . . 76 6.1 Building models reconstructed targeting to obtain (c) more precise geometry and (d) more precise topology respectively. Compared with (a) the input LiDAR and (b) unsimplified building model, the missingofthechimneymakestheformeronevisuallylessconvinc- ingthanthelatterone. . . . . . . . . . . . . . . . . . . . . . . . . . 79 6.2 Comparisonbetween(a)2.5Ddualcontouringand(b)2.5Dbuilding modeling with topology control. While the uniqueness of hyper- pointinonecellpreventsaflexiblesimplificationindualcontouring, the new method detects and controls building topology beyond the rigidquadtreestructure. . . . . . . . . . . . . . . . . . . . . . . . . 80 6.3 Topologicalfeaturesinanunsimplified2.5Dbuildingmodel . . . . 82 6.4 2.5Dbuildingmodelsmaycontainpointfeaturesinvolvingonlyone wall feature (left). These points are produced around grid edges (e.g., AB) which detect inconsistent roof layer assignments in two adjacentcells(right). . . . . . . . . . . . . . . . . . . . . . . . . . 83 xii 6.5 2.5D building topology is represented by topological features, and the associations between them in form of Equation (6.1) and (6.2). Examples include typical building structures such as (a) individual building blocks, (b) blocks with top attachments, (c) blocks with side attachments, (d) stair-shaped structures, and (e) combinations ofthesepatterns. . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 6.6 Foldingpointscanbepartofa1-layerhyper-point,a2-layerhyper- point, or a multi-layer hyper-point (from left to right). They are detectedandmarkedaspointfeaturesbeforeadaptivesimplification starts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 6.7 Hyper-pointclusterforestviewedfromobliqueandorthogonalper- spectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 6.8 A stadium model created using (a,b) 2.5D dual contouring, (c,d) themodelingmethodpresentedinthischapter,(e)manualcreation, and (f) the plane-based approach in Section 3.4. The relation curve between error tolerance and the triangle number of reconstructed models is illustrated in (g). With larger geometry error tolerance given,thenewmethodcanalwaysproducesimplermodelswithless triangles; while the overstrict topology test in 2.5D dual contouring creates numerous trivial triangles along thin roof features shown in closeupsof(a,b). . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.9 Modelevolutionwitherrortolerancegrowingfrom1.0toinfinite 93 6.10 2,099buildingmodelsarecreatedforanurbanareainDenverusing (top) building modeling with topology control, (middle) 2.5D dual contouring, and (bottom) plane-based method. The new method produces as few triangles as the plane-based method while recov- eringandpreservingthetopologicalfeaturesineachbuildingstruc- ture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 7.1 Themethodpresentedinthischapterautomaticallydiscoversglobal regularities from a noisy 2.5D point cloud, and uses them to cre- ate a convincing 2.5D building model. Two orthogonal projections illustratea subsetoftheglobalregularitiesinthismodel(lengthsin meters). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 xiii 7.2 Modelingresultsgeneratedfromthesameinputpointcloudbyman- ual creation, the modeling method proposed in this chapter, 2.5D dual contouring with principal direction snapping, and a primitive- based method [34]. The new modeling method creates the most visually convincing result among all three automatic methods since its resulting building model conforms to the most global regulari- ties. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 7.3 A typical gable-shaped building roof containing a set of 2.5D ele- ments(e.g.,planeprimitive,boundarysegments,andridges) . . . . . 101 7.4 Pipeline of the modeling approach: a 2.5D point cloud (top left) is transformed to a building model (bottom left) through a series of steps. Global regularities are discovered and enforced in each alignmentstep. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 7.5 Coarse-to-fineplanarroofpatchextraction . . . . . . . . . . . . . . 109 7.6 A comparison of geometric fitting errors (i.e., squared distances from input points to the model surfaces) between models created byfourdifferentapproaches . . . . . . . . . . . . . . . . . . . . . 113 7.7 Experiments on several building scans. By discovering global reg- ularities and enforcing them on the planar roof patches and their boundary segments (second column), the modeling algorithm cre- ates visually convincing 2.5D building modelings (third column). Aerialimageryand2.5Ddualcontouringresultsareincludedasref- erence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 8.1 Given (a) a dense aerial LiDAR scan of a residential area (point intensitiesrepresentheights),theresidentialurbanmodelingsystem reconstructs (b) 3D geometry for buildings and trees respectively. (c)Aerialimageryisshownasareference. . . . . . . . . . . . . . . 116 8.2 Local geometry features become unreliable when dealing with res- idential areas with rich vegetation. In closeups of (A) a tree crown region and (B) a rooftop, points are rendered as spheres while a locally fitted plane is rendered in yellow. Middle: classification results from the general urban modeling system in Chapter 3, trees in green, buildings in purple, and ground in dark grey. Right: mod- elingartifactsarecreatedbecauseofclassificationerrors. . . . . . . 117 xiv 8.3 While building structures have a 2.5D characteristic, trees do not possesssuchproperty. Denselaserscansmaycapturesurfacepoints underthetreecrown(right).. . . . . . . . . . . . . . . . . . . . . . 118 8.4 Ademonstrationoftheclassificationalgorithm: (b)thelowestlayer fragments faithfully capture the skyward surfaces of 2.5D struc- tures;(c)buildingandgroundlayerfragmentsarerenderedinpurple andgreyrespectively;(d)treesandoutliersareinblackwhilebuild- ingroofpatchesarerenderedinbrightcolors. . . . . . . . . . . . . 122 8.5 The residential urban modeling system is tested against multiple data sets from different sources. The system robustly reconstructs urban reality despite the variation in data resolution, building pat- terns,andtreetypes. . . . . . . . . . . . . . . . . . . . . . . . . . . 125 8.6 Urban models reconstructed from 5.5M aerial LiDAR points for a residentialareainthecityofAtlanta . . . . . . . . . . . . . . . . . 127 xv Abstract 3D reconstruction from point clouds is a fundamental problem in both computer vision and computer graphics. As urban modeling is an important reconstruction problem that has various significant applications, this thesis investigates the complex problem of reconstructing 3D urban models from aerial LiDAR (Light Detection And Ranging) pointcloud. Inthefirstpartofthisthesis,anautomaticurbanmodelingsystemisproposedwhich consistsofthreemodules: (1)theclassificationmoduleclassifiesinputpointsintotrees and buildings; (2) the segmentation module splits building points into different roof patches;(3)themodelingmodulecreatesbuildingmodels,ground,andtreesfrompoint patches respectively. In order to support city-scale data sets, this pipeline is extended into an out-of-core streaming framework. By storing data as stream files on hard disks and using main memory as only a temporary storage for ongoing computation, an effi- cient out-of-core data management is achieved. City-scale urban models are success- fullycreatedfrombillionsofpointswithlimitedcomputingresource. The second part of this thesis explores the 2.5D nature of building structures. The 2.5D characteristic of building models is observed and formally defined as “building structures are always composed of complex roofs and vertical walls”. Based on this observation,a2.5Dgeometryrepresentationisdevelopedforthebuildingstructures,and usedtoextendaclassicvolumetricmodelingapproachintoa2.5Dmethod,named2.5D xvi dual contouring. This algorithm can generate building models with arbitrarily shaped roofs while keeping the verticality of the walls. The next research studies the topology of 2.5D building structures. 2.5D building topology is formally defined as a set of roof features, wall features, and point features; together with the associations between them. Based on this research, the topology restrictions in 2.5D dual contouring are relaxed. The resulting model contains much less triangles but similar visual quality. To furthercapturetheglobalregularitiesthatintrinsicallyexistinbuildingmodelsbecause ofhumandesignandconstruction,abroadvarietyofglobalregularitypatternsbetween 2.5D building elements are explored. An automatic algorithm is proposed to discover and enforce global regularities through a series of alignment steps, resulting in 2.5D building models with high quality in terms of both geometry and human judgement. Finally,the2.5Dcharacteristicofbuildingstructuresisadoptedtoaid3Dreconstruction of residential urban areas: a more powerful classification algorithm is developed which adopts an energy minimization scheme based on the 2.5D characteristic of building structures. This thesis demonstrates the effectiveness of all the algorithms on a range of urban areascansfromdifferentcities;withvaryingsizes,density,complexity,anddetails. xvii Chapter1 Introduction Three dimensional urban models are very useful in a variety of applications such as urban planning, virtual city tourism, surveillance, disaster simulation, and computer games. Most of these applications require city-scale urban models composed of simple andelegantbuildingmodels,treemodels,andground. TheadvanceofacquisitiontechniqueshasmadeaerialLiDAR(LightDetectionAnd Ranging) data a powerful 3D representation of urban areas. Equipped on aeroplanes, laserscannersareabletocapturethesurfacegeometryoflargecitiesinanaccurateand efficient manner. Billions of LiDAR points are collected for city-scale data sets. These datasetsusuallycontaineverydetailinanurbanareaincludingbothimportantfeatures (e.g., detailed rooftops) andtrivial details(e.g., residualsensor noise, undesired vegeta- tion,andvehicles). Inaddition,themagnitudeofLiDARdatalimitsitsapplication. Theobjectiveofthisresearchistoreconstruct3DurbanmodelsfromaerialLiDAR data. Specifically,themodelingapproachdesirefollowingproperties: • Automation: Urban models should be generated in an automatic manner. Few userinteractionsareintroducedonlywhennecessary. • Discrimination: The method should be able to remove undesired features while preservingimportantelements, e.g.,buildings,andterrain. • Efficiency: The method should create city-scale urban models from huge data setswithinreasonabletimeandspace. 1 • Seamlessness: Theurbanmodelshouldbegeneratedinaseamlessmanner,mean- ingthattheLiDARdatashouldbeprocessedasoneundividedpiece. • Universality: The method should be able to process LiDAR data sets scanned from different cities. Data sets with different density and scales can be handled bythesameframeworkwithonlyafewparametersconfigured. • High-qualityoutputs: Themethodshouldproducesimplepolygonalmodelsfit- tingtheinputpointcloudinaprecisemanner. Sincebuildingsarethekeycompo- nentinurbanareas,therearethreespecificrequirementsforthebuildingmodels. 1. The reconstructed building models should contain small amount of vertices andtriangles,toenableapplicationssuchasefficientrendering. 2. Building models should fit the observation (e.g., the raw scan data and the aerialimagery)inaprecisemanner. 3. Building models should be visually convincing in terms of human judge- ment,becausehumansaretheultimatejudgeformodelqualityinmostappli- cations. The first part of this thesis proposes a general 3D urban modeling system towards these requirements. This system contains four main modules, namely, classification, segmentation, building modeling, and terrain modeling. Given a LiDAR point cloud as input, irrelevant features such as trees and noises are first classified from roof and ground points, and removed in subsequent processing. A segmentation algorithm then takes the remaining point cloud, and extracts individual buildings from terrain. These point patches are sent into a building modeling module and a terrain modeling module respectivelytocreateelegantpolygonalmodels. 2 a (a) (a) (b) (b) (c) (c) (d) (d) d b c Input LiDAR data Urban model 683M points 1.12M triangles (building) 8.87M triangles (terrain) Figure 1.1: 3D reconstruction for the city of Atlanta from 683M LiDAR points, using the streaming urban modeling systempresentedinChapter 4. Bottom: fourcloseupsof the urban model, showing (a) an area with large flat structures; (b) a downtown area; and(c,d)tworesidentialareas. Inordertohandleextremelylargedatesetsinaseamlessmanner,thegeneralurban modelingpipelineisextendedtoanout-of-coreexecutionarchitecture. Theout-of-core architecture benefits from a streaming process which utilizes free hard-disk space as main storage while allocating main memory only as temporary space to store data for ongoingcomputation. Inotherwords,thestreamingprocesstakesdataasastream(usu- ally, a disk file) - each system component reads from an input stream, loads necessary data in-core, processes it; and once data is no longer needed for furtherprocessing, itis writtenintoanoutputstream. 3 This automatic system is tested on four different data sets to show the universality. Urban models are automatically generated from aerial LiDAR data. The largest data set that has been processed is Atlanta LiDAR data containing 683M points stored in a 17.7GB disk file. The modeling system generates 1.12M triangles to represent the buildings and 8.78M triangles for terrain, by 25 hours of unattended processing using lessthan1GBmemory. TheresultingurbanmodelisshowninFigure1.1. Asillustrated in the closeups, different types of building models are created successfully over the entireurbanarea. Basedonthisurbanmodelingsystem,thesecondpartofthisthesisstudiesthemore fundamental problem about the 2.5D nature of building structures. As shown in Fig- ure 1.2, the LiDAR sensor captures the details of roof surfaces, but collects few points on building walls connecting roof boundaries; in addition, users desire building models composedofcomplexroofsandverticalwallsconnectingdifferentrooflayers. Boththe input and the objective of the building modeling problem show a 2.5D characteristic, causedbythe2.5Dnatureofbuildingstructuresduetohumandesignandconstruction. This 2.5D characteristic of building structures inspires a geometry representation of polygonal building models, which forces the verticality of wall polygons. A robust approach is developed in this thesis to create 2.5D building models from aerial LiDAR pointclouds,named2.5Ddualcontouring. Thismethodisguaranteedtoproducecrack- freemodelscomposedofcomplexroofsandverticalwallsconnectingthem. Byextend- ingclassicdualcontouringintoa2.5Dmethod,buildinggeometryissimultaneousopti- mized over the three dimensional surfaces and the two dimensional boundaries of roof layers. Thus,2.5Ddualcontouringcangeneratebuildingmodelswitharbitrarilyshaped roofs while keeping the verticality of connecting walls. An adaptive grid is introduced to simplify model geometry in an accurate manner. Sharp features are detected and preservedbyanovelandefficientalgorithm. 4 Dense sampling on roofs Sparse sampling on walls Vertical walls Detailed roofs Figure 1.2: 2.5D nature of the building modeling problem. Left: aerial LiDAR has dense sampling on rooftops but very sparse sampling on building walls. Right: desired buildingmodelsareusuallycomposedofdetailedroofstructuresandverticalwalls. While the building geometry describes where the building structures appear in the three dimensional space, the building topology determines the existence of structural piecesandtheconnectionsbetweenthem. Myfurtherresearchrevealsthathumanvision tendtobeverysensitiveintopologychangeseventherelatedstructuralpiecesaresmall. Therefore,2.5Dbuildingtopologyisexploredanddefinedasasetofrooffeatures,wall features, and point features; together with the associations between them. Based on this definition, 2.5D dual contouring is extended into a 2.5D modeling method with topology control. Comparing with 2.5D dual contouring, the new modeling method put less restrictions on the adaptive simplification process. The reconstruction results preserve significant topology structures while the number of triangle is comparable to thatofmanuallycreatedmodelorprimitive-basedmodels. Besides the 2.5D characteristic of building structures, other building natures due to human design and construction are further explored. They appear in the form of global regularities that encodes the orientation and placement similarities between planar ele- mentsinbuildingstructures. Anautomaticmodelingalgorithmisdesignedbasedonthe 5 study of 2.5D global regularities. This algorithm simultaneously detects locally fitted plane primitives and global regularities: while global regularities are extracted by ana- lyzing the plane primitives, they adjust the planes in return and effectively correct local fitting errors. By aligning planar elements to global regularities, the modeling method significantly improves the reconstruction quality in terms of both geometry and human judgement. When the urban modeling is applied to residential areas, rich vegetation becomes an important urban feature, and thus should be carefully handled. Therefore. a robust classification algorithm is proposed in this thesis, which effectively classifies LiDAR points into trees, buildings, and ground. The classification algorithm adopts an energy minimization scheme based on the 2.5D characteristic of building structures: buildings are composed of opaque skyward roof surfaces and vertical walls, making the interior of building structures invisible to laser scans; in contrast, trees do not possess such characteristic and thus point samples can exist underneath tree crowns. Once the point cloudissuccessfullyclassified,thesystemreconstructsbuildingsandtreesrespectively, resultinginahybridmodelrepresentingthe3Durbanrealityofresidentialareas. Therestpartofthisthesisisorganizedasfollows: Chapter2summarizestherelated literature to the urban modeling problem. Chapter 3 details a general urban modeling pipeline that is extended into an out-of-core streaming architecture in Chapter 4. The research based on the 2.5D nature of building structures include: 2.5D geometry rep- resentation and 2.5D dual contouring in Chapter 5; 2.5D building topology research in Chapter 6; 2.5D building modeling with global regularities in Chapter 7; and residen- tial urban area reconstruction in Chapter 8. The conclusion and future work is finally presentedinChapter9. 6 Chapter2 RelatedWork This chapter briefly reviews the related work in the scope of 3D urban reconstruction from aerial LiDAR data. Previous research related to each specific topic is reviewed individually in Section 4.1 for streaming approaches, Section 5.1 for classic dual con- touring,Section6.1fortopologycontrolinvolumetricmodeling,Section7.1forshape- from-symmetrymethods,andSection8.1fortreedetectionandmodelingalgorithms. Pioneerbuildingmodelingmethods[1,11,22,49,50,64]startbyconvertingLiDAR point cloud into a DEM (Digital Elevation Model), and then apply image processing algorithms on these depth images to detect building footprints, fit parametric models and reconstruct polygons. All of them share a similar building reconstruction pipeline withthreemajorsteps: classification,segmentation,andbuildingmodeling. Mostexistingresearchworkisbuiltuponthispipelineandadvancesthereconstruc- tion quality by improving individual steps. Fr¨ uh et al. [18] create building models with facade by integrating aerial LiDAR and ground based LiDAR. Hu et al. [24] employ aerial image and ground images to achieve a high resolution solution. Wang et al. [62] concentrate on building footprint extraction problem. Verma et al. [60] propose a roof- topologygraphtofindcomplexroofpatternsfromaerialLiDARdata. Lafargeetal.[33] present a two-stages method which can find optimal configuration of parametric mod- els via a RJMCMC sampler. Matei et al. [41] specialize segmentation for densely built areas. Zebedin et al.[65]detectplanesandsurfacesofrevolution. PoullisandYou[48] create simple 3D models by simplifying boundaries of fitted planes. Toshev et al. [58] propose parse trees as a semantic representation of building structures. Lafarge and 7 Mallet [34] combine primitives and a general mesh representation to achieve hybrid reconstruction. Since building models are usually considered to be the most important component in urban areas, these research efforts have an emphasis on building reconstruction. In general,mostoftheseeffortsattackthebuildingmodelingproblemwithprimitivefitting approaches. Specifically, local plane fitting is a popular strategy in extracting simple roof primitives [41, 48, 60]. Strong urban priors are frequently introduced to restrict the plane primitives. Typical priors include roof topology [60] and Manhattan-world grammars [41, 48]. Other research work aims at extending the ability of representing complicatedshapes,byintroducingadditionalprimitiveshapesandoptimizingjunctions betweenfittedprimitives[33,34,65]. Another branch of the urban reconstruction problem lies in vegetation detection, which is believed to be a difficult task in the modeling systems. Several research efforts[4,5,34,39,52,58,60,62]addressthisproblem. Ingeneral,theseeffortsemploy both geometry and reflection properties at each LiDAR sample point, and use a general classifier to assign each point a class label indicating its category. Typical classifiers includesimplethresholding[60],SVMclassifiers[52],orAdaboostclassifiers[39,62]. Although automatic solutions have been provided for various LiDAR data sets, to thebestofourknowledge,noneofthemcanprocessanextremelylargeLiDARdataset in a seamless manner. Instead, many of them (e.g. [41, 48]) partition huge LiDAR data into tiles, process them one at a time, and merge the partial results together to generate the aggregate model of a large scale city area. Artifacts can occur at seams even with paddingbetweentiles. Another limitation of these methods is that they are usually designed and tested for one specific data set. Prior knowledge learnt from this input data set is integrated into 8 thealgorithm. Thustheyloseuniversalityandmaynotperformwellwhendealingwith otherdatasets. Finally,thisthesisisthefirsttodiscoverandformallydefinethe2.5Dcharacteristic of building structures. Theories and algorithms are developed based on the 2.5D nature of buildings. The urban modeling results in this thesis are competitive in terms of both geometryfittingqualityandhumanvisualjudgement. 9 Chapter3 GeneralUrbanModelingApproach This chapter presents a fully automatic system which creates simple and elegant urban models from aerial LiDAR data. The system takes only the aerial LiDAR point cloud as input and works directly on the raw data without rasterizing it into a DEM (Digital Elevation Model). The output contains two parts: a polygonal model representing the terrain,andasetofsimpleandelegantbuildingmodels. Inalocalarea,roofboundaries arealignedtogetherasmuchaspossibletoimprovethequalityofthebuildingmodels. Section 3.1 proposes a general urban modeling pipeline, followed by the details of eachmoduleinsubsequentsections. Section3.6showsexperimentsonacoupleofsmall data sets chopped from different cities. These data sets contain small amount of points, thuscanbeloadedin-coreandprocessedatonce. 3.1 UrbanModelingPipeline Thissectionproposesanautomaticurbanmodelingpipelinewhichsequentiallyexecutes fourindividualprocessesasillustratedinFigure3.1: 1. Classification: A classification module classifies vegetation and noise points frombuildingandgroundpoints. Thedetectedvegetationandnoisepoints(green and black points in the second sub-figure in the top row) are considered as irrele- vantpartsandthusareremovedinsubsequentprocessing. 10 Input LiDAR data Classified points Classification Segmented patches Segmentation Urban model Terrain modeling Terrain Building modeling A building patch Detected planes Boundary points Polygonal model Plane fitting Boundary detection Modeling Figure 3.1: The general urban modeling pipeline includes four modules: classification, segmentation,buildingmodeling,andterrainmodeling. Inparticular,thebuildingmod- eling module takes an individual point patch as input; detects plane patterns; extracts boundaries of these planes; and finally create polygonal building models based on the boundaries. 2. Segmentation: Using a distance-based region growing approach, the remaining building and ground points are segmented into individual patches (points with differentcolorsinthethirdsub-figureinthetoprow). 3. Buildingmodeling: Meshmodelsarecreatedforeachindividualbuildingpatch, followingthreesteps: (a) Plane fitting: Planepatternsaredetectedandfittedfromthepointcloud. (b) Boundary detection: A watertight boundary is extracted for each plane detectedinthepreviousstep. 11 (c) Modeling: Mesh models are created from the boundary points. The roof boundaries are simplified and neighboring boundary edges are snapped togethertoreducethecracksinthefinalresult. 4. Terrain modeling: Taking the ground point patch as input, the terrain modeling module generates a polygonal model with rasterization. Holes caused by occlu- sionarefilledbysolvingaLaplace’sequation. To further improve the reconstruction quality of building models, I adopt the assumption that the roof boundary edges usually fall into a set of principal directions. An efficient method is proposed to learn the principal direction set from original data, and align roof boundary segments along these directions automatically. This mecha- nismenablesarbitraryanglesbetweenneighboringroofboundarysegments. Incontrast, many previous methods are limited by the pre-assumptions made on these angles, e.g., 90 ◦ ,45 ◦ and 135 ◦ in[62]. Another extension to this pipeline is the non-planar roof extension. In an industrial sitewhichcontainsseveraloiltinsandtanks,thegeometryofthesemodelsisintheform of cylinders and cones. Thus, the user has the option to determine the shape of each objectviaminimalinteractions. RANSAC[14]isappliedtoestimatetheparametersfor thepre-definedprimitives. 3.2 VegetationDetection Vegetation, mainly in the form of trees, is an irrelevant part in most downtown urban areas, thus should be eliminated in the first step. This section presents a novel classifi- cationalgorithmforvegetationdetection. 12 3.2.1 Features GivenasamplepointintherawLiDARpointcloudp∈P,N p ={q|q∈P,∥p−q∥ < δ}denotesthesetofpointswithinasmallneighborhoodofpointp,andpisthecentroid ofallpointsinN p . Note that if p lies in an extremely sparse area, it loses its geometry significance. Therefore,p is considered as a noise point. Specifically,p is detected as a noise point iff |N p |<κ,whereκisauser-giventhreshold,empirically,equalsto10. For a sample pointp with enough points within its neighborhood N p , five features aredefinedbasedonthelocaldifferentialgeometryproperties: 1. Regularity: Itmeasuresifthepointdistributionaroundpisregularbycalculating F 1 =∥p−p∥. (3.1) Intuitively, sample distribution around a roof point is more likely to be regular. Thus the distance between a roof point and its neighbors’ centroid should be smallerthanthatofavegetationpoint. 2. Horizontality: SinceaerialLiDARdatacapturesrooffromatopview,thenormal ataroofpointismorelikelytobevertical. Thesecondfeatureisdefinedas F 2 = 1−|n p ·e z |, (3.2) wheree z = (0,0,1) is the vertical direction, andn p is obtained through covari- anceanalysis[47],i.e.,tosolvetheeigenvectorproblemforthecovariancematrix: C p = 1 |N p | ∑ q∈Np (q−p)(q−p) T . (3.3) 13 z=0 plane z=0 plane (a) (b) (c) Figure3.2: Illustrationofnormalvariationmeasurements: (a)normalsofpointsaround a roof ridge; (b) normal distribution on a Gauss sphere, red/green/blue arrows point to the eigenvectors of C n p with lengthes equal to corresponding eigenvalues; (c) normals of points in the neighborhood of a tree point. Both λ n 1 and λ n 2 are large due to the irregularityofnormals. Thethreeeigenvaluesaresortedinascendingorder: λ 0 ≤ λ 1 ≤ λ 2 . Theeigenvec- tor v 0 corresponding to the smallest eigenvalue is the estimated normal at point p. 3. Flatness: With covariance analysis results, the flatness at this point, also known as surface variation[47],couldbeestimatedas: F 3 = λ 0 λ 0 +λ 1 +λ 2 . (3.4) Similartothepreviousfeatures,asmallerflatnessvalueindicatesmorepossibility forapointtobearoofpoint. 4. Normaldistribution: Oncethenormalateachpointisestimated,anothercovari- ance analysis can be further applied over these normals. The solution of the 14 eigenvectorproblemforanormalcovariancematrixwithinalargerneighborhood N n p ={q|q∈P,∥p−q∥ <η}: C n p = 1 |N n p | ∑ q∈N n p n T q ·n q , (3.5) resultsineigenvaluesλ n 0 ≤ λ n 1 ≤ λ n 2 withcorrespondingeigenvectorsv n 0 ,v n 1 ,v n 2 . AsGarland[20]andPauly[47]pointout,λ n 1 measuresthemaximumvariationof these normals on the Gauss sphere, thus could be regarded as a kind of normal variation. Therefore,oneofthenormaldistributionfeaturesisdefinedas: F 4 = λ n 1 . (3.6) Apparently,aroofpointprefersasmallernormalvariationthanatreepoint. Moreover,accordingtoprincipalcomponentanalysis[30],thesmallesteigenvalue λ n 0 measures the minimum normal variance, with its eigenvector v n 0 being the direction on Gauss sphere along which normals spread out the least. Specifically, λ n 0 determines whether the normals spread in a narrow band area on the Gauss sphere, or scatters irregularly. This property is particularly useful for sharp roof features at planar roof patch intersections. For example, Figure 3.2 shows a roof ridgebroadlyexistinginurbanarea. Normalsofthesepointsformtwoclusterson theGausssphere. Henceλ n 1 (illustratedasthelengthofthegreenvectorinFigure 3.2(b))islargewhileλ n 0 (lengthofthebluevectorinFigure3.2(b))isverysmall. Bycontrast,normaldistributionaroundatreepointisfairlyirregularandexhibits largeλ n 0 andλ n 1 (showninFigure3.2(c)). Therefore,thelastfeatureisdefinedas: F 5 = λ n 0 . (3.7) 15 1.0 Tree-like 0.0 Roof-like Input LiDAR data F1 F2 F3 F4 F5 Figure 3.3: Feature values of a small piece of aerial LiDAR point cloud. The top-left sub-figure is rendered with z-values as grey scales at points. Others are rendered with colors representing corresponding feature values of the points. Blue points are more likelytoberoof/ground,whileredpointsaremorelikelytobetrees. Figure 3.3demonstrates thesefeature valuesofan exampleurban area. The top-left sub-figure shows the input LiDAR data colored with its height value. The rest sub- figuresshowthefivefeaturesusingpointcolorsrespectively. Bluecolordenotesasmall feature value (thus, roof-like); while red color represents a large feature value (tree- like). The five features successfully distinguish vegetation points and building/ground pointevenwhentheyhavesimilarheight. 3.2.2 Classifier Despite the removed noise points, two classes remains in this classification problem, i.e.,treesandbuilding/ground. Thus,alinearclassificationhyperplaneinthe5Dfeature 16 spaceisintroducedtosplitthetwoclassesofpointsamples. Inparticular,adiscrimina- tivefunctionisdefinedasthelinearcombinationofthesefivefeatures: K =ω 0 +ω 1 F 1 +ω 2 F 2 +ω 3 F 3 +ω 4 F 4 +ω 5 F 5 , (3.8) where ω 0;1;2;3;4;5 are undetermined parameters. Once these parameters are determined, the classification algorithm simply computesK for each point and determines its cate- gorywithsgn(K). Toautomaticallydeterminetheseparameters,supervisedmachinelearningmethods are introduced which learn these parameters from a small but typical urban area with manual labeling. In particular, any linear classifier can serve as the learning module, including Linear Discriminant Analysis [15] and Support Vector Machine [2, 8]. For efficiency and accuracy, an unbalanced soft margin Support Vector Machine algorithm isadoptedasproposedin[8]. TheSupportVectorMachinealgorithmfindsahyperplane that maximizes the margin between the two classes. Due to the mislabeled examples and noise samples in the classification problem, there may exist no hyperplane that can strictly separates the two classes, thus a soft margin extension is introduced to split the examples as cleanly as possible, while still maximizing the margin distance. The soft margin Support Vector Machine algorithm is detailed in [8] and a third-party library SVM Light [29]isemployedinthesystemimplementation. To further improve the classification results, a voting algorithm is introduced as a post-processing step. Intuitively, points of a same category usually occur together in the space. Thus, once all the points are labeled by the classifier, pointp’s final label is determined by voting fromp’s neighboring points. That is,p belongs to tree category only if the percentage of tree points inp’s neighborhoodN p is greater than a threshold ω. Thisthresholdisalsolearntfromthelabeleddataset. 17 (a) input LiDAR data (b) classification result (c) segmentation result Figure 3.4: (a) Input LiDAR point cloud with color intensity representing the height at eachpoint. (b)Greenpointsaredetectedastrees. (c)Segmentationresult;withtreesand noises(blackpoints)removedfromtheinput,thesegmentationalgorithmsplitsbuilding roofsfromthegroundsuccessfully. 3.2.3 ConfigurationandResults Inthisclassificationalgorithm,allthefeaturesaredesignedbasedonthelocaldifferen- tial geometry properties. Without introducing data-set-related variants such as absolute heightandintensity,theproposedmethodachievesignificantadaptivitytodifferentdata sets. In the experiments, all the seven parameters ω 0;1;2;3;4;5 and ω are learnt from a 100m× 100m labeled area from the city of Oakland. Once they are determined, they workforallthetestingdatasets. The only parameters that need to be configured for each data set are the sizes of the neighborhood, i.e., δ andη. Since the neighborhood should be large enough to support an effective covariance analysis, δ is set to the value which makes the average point number in N p between 12 ∼ 20. As for η, since normals are estimated from points in N p and thus are smoothed in some sense; the size of N n p should be larger than that of N p togainsufficientgeometricsignificance. Thus,η istwotimesthevalueofδ. 18 Theclassificationalgorithmachievesanaccuracyabove95%onalltestingdatasets. Figure 3.4 shows a 300m× 300m urban piece chopped from Oakland data set. Trees and non-trees are correctly classified even at areas where points from both categories havesimilarheights,asshowninthecloseupofFigure3.4(b). Theweightsω 1;2;3;4;5 arelearntwithaSupportVectorMachinealgorithm,anddeter- minedonceandforall. Intheexperiments,theyare: ω 1;2;3;4;5 ={2.5, 0.1, 1.7, 5.2, 20.4}. (3.9) First, since all the weights are positive, each feature has a positive contribution to the linear classifier. This observation concurs with human intuition. Second, this result suggeststhatthenormalvariationsappeartobethemostimportantfactorintheclassifi- cationproblem,especiallytheminimaleigenvalueofthenormalcovariancematrix,i.e., F 5 . 3.3 Segmentation Once trees and noises are eliminated from the LiDAR data, the next task is to identify ground points and segment different building patches. This is achieved by a distance- basedsegmentationalgorithmbasedonthefollowingcriteria: apairofpointswiththeir distance smaller than a user-given threshold α are assigned to the same segment. With non-overlapping point segments generated from building and ground points, the largest patch is considered to be the ground patch. The rest point patches are identified as different roof patches. Roof patches that have neighboring projections on the x-y plane are considered to belong to the same building structure, and thus are merged into one buildingpatch. 19 AnexampleofthesegmentationapproachisillustratedinFigure3.4(c). Treepoints and noise points are rendered with black color and removed from the input beforehand. The segmentation algorithm splits the remaining points into several patches, with the largest patch identified as ground and rendered with dark-grey color. The rest patches aregroupedtogethertoformindividualbuildingroofsandrenderedwithdifferentbright colors in the figure. Since trees and points on the vertical walls (usually detected as noise)areeliminatedpreviously,therealwaysexistsaheightgapbetweenbuildingroofs andground,whichleadstoacorrectresultofdistance-basedsegmentation. 3.3.1 RegionGrowingImplementation The distance-based segmentation algorithm can be easily realized by a region growing implementation. Startingfromanunlabeledseedpointp,theregiongrowingalgorithm searchesitsα-neighborhoodN p andassignsallthepointswithinN p tothesamecluster asp,denotedasC p . SimilarprocessisappliedoverallpointsinC p repeatedlyuntilthe clustercangrownomore. Thenthesegmentationalgorithmshiftstothenextunlabeled pointq and grows it into a cluster C q . This algorithm runs iteratively until every point hasaclusterlabel. 3.3.2 AgglomerativeClusteringAlgorithm To improve time efficiency, an agglomerative clustering algorithm is introduced using theaunion-findimplementation[7]. Theagglomerativeclusteringalgorithmexecutesin abottom-upmanner. Itstartswitheachpointasanindividualcluster(orsegment),then traverses all the point pairs with point distance smaller than α, and merges the clusters towhichthetwopointsbelong. The union-find implementation utilizes a disjoint-set forests data structure as illus- tratedinFigure3.5. Eachtreeintheforeststoresthepointsinasamecluster. Eachpoint 20 (a) Input points (b) An intermediate state (c) Set union operation (d) Find operation with flattening Points Roots of set trees A point applied with a find operation Points on the root-seeking path Figure 3.5: Agglomerative segmentation implemented using the union-find algorithm presentedin[7] pholdsareferencer(p)pointingtoitsparentnodeinthetree. Thesegmentmergingis convenientlyimplementedbytheunionoperation,whichsimplyassignsonerootasthe parentofanotherroot,asshowninFigure3.5(c). To accelerate the process of retrieving root point for each cluster, the find operation includesaflattenprocesswhichlinkseachpointontheroot-seekingpathdirectlytothe root,asshowninFigure3.5(d). Thisflattenprocessisespeciallyusefulinthestreaming segmentationtoavoiddanglingpointers,asdetailedinSection4.3.3. 3.4 BuildingModeling This section presents the building modeling algorithm that generates a simple and ele- gant mesh model from each individual building point patch. As shown in the bottom sub-figure of Figure 3.1, the building modeling approach executes three steps, namely, planar roof patch extraction, boundary detection, and model creation. The details of these steps will be presented in the following sections. A non-planar roof extension withuserinteractionisproposedinSection3.4.4. 21 3.4.1 PlanarRoofPatchExtraction Givenabuildingpointcloud,thefirsttaskinbuildingmodelingistoextractplanarroof patches from the input. As suggested by Verma et al. [60], this task can be achieved using a similar distance-based segmentation as that proposed in Section 3.3, with dis- tancedefinedbasedontheinnerproductofpointnormalsinsteadofEuclideandistance. That is, as point q is in point p’s neighborhood N p , they are neighboring only if they satisfy: 1−|n p ·n q | < β. Empirically, any β between 1− cos(5 ◦ ) and 1− cos(10 ◦ ) workswell. However, in practice, this method may lose robustness when dealing with large and smoothcurvedsurfaces,suchastheroofsofstadiums. Inthesecases,thewholecurved surface is detected as one planar patch, which conflicts with human intuition. Thus, an iterativemethodisadoptedtosolvetheplanefittingproblemrobustly. Thegeneralidea issimilartotheEMalgorithm[17]. Starting from a planel determined by the position and normal of an unlabeled seed pointp,theplanefittingalgorithmiterativelyexecutes: • E-step: ApointsetLisdeterminedasthesetofunlabeledpointswithEuclidean distancestol smallerthanagiventhresholdt,andnormaldifferenceslessthanβ; • M-step: Planel isupdatedusingleastsquaresfittingwithrespecttopointsetL. Aplanarroofpatchisdetectedwhenthealgorithmconvergesorhasiteratedenough times. This algorithm is applied iteratively until each point is assigned to a plane. The labeling of points around boundaries between these plane patches are further refined usingak-means-likealgorithmnamedthek-proxyclustering,asproposedin[6]. With small patches removed as insignificant features, the remaining planar roof patchesareconsideredasimportantelementsincreatingbuildingpolygons. 22 3.4.2 BoundaryDetection The next step in the building modeling algorithm is to mark boundary points of each planar roof patch. These boundary points will be later used in footprint generation and polygon production. Some previous work, e.g., [62], finds boundary points by mea- suring certain characteristics of roof points. These methods, although efficient, cannot guaranteeawatertightboundary,thuslimitthesubsequentprocessing. Otherworkuses 2D delaunay triangulation to find polygonal boundaries. E.g., Verma et al. [60] intro- duce a plane ball-pivoting algorithm to triangulate roof points and detect boundaries. These algorithms are able to generate perfect boundaries, but sacrifice time efficiency. Inaddition,robustimplementationsforsuchalgorithmsarehardtoachieve. This section proposes an algorithm which combines the advantages of these two kinds of algorithm. The proposed method is inspired by part of my early work about a 3Dcontouringalgorithmin[68]. Initially, for each plane roof patch, all the roof points are projected on the x-y plane and embedded into a uniform grid. The grid cells with roof points inside are marked as object cells,andillustratedasthegreycellsinFigure3.6. Withthisuniformgridsetup, boundarypointsandboundarylinesaregeneratedrespectively: • Boundary points: For each grid edge l separating an object cell c in and a back- ground cell c out (thick red edges in Figure 3.6), the nearest LiDAR point to l in c in ismarkedasaboundarypointp(l),shownastheredcirclesinFigure3.6. • Boundarylines: Foreachgridcornerc(redsquaresinFigure3.6)sharedbytwo boundary grid edgesl 1 andl 2 , a boundary line (p(l 1 ),p(l 2 )) is created, shown as thethinredlinesinFigure3.6. 23 Figure 3.6: An illustration of the boundary detection algorithm. Circles represent the LiDAR points projected on the x-y plane. The boundary of this patch is composed of theboundarypoints(redcircles)connectedbyboundarylines(thinredlines). Intuitively, the result of this algorithm is the dual polygon of the object cell bound- ary along grid edges. Since the latter is a watertight polygon, its dual polygon is also watertight. Thispropertyisespeciallyusefulinthemeshcreation. The boundary extraction algorithm is efficient and robust. It produces a topology- correct boundary, but cannot guarantee the geometry completeness. For example, in Figure 3.6, there is one point outside the extracted boundary polygon. Since this sit- uation happens rarely and has little effect on the subsequent processing, the resulting boundaryisagoodapproximation. Morphological operations: Another advantage of this boundary detection algorithm is that it supports morphological operations in an easy manner: after marking object and background cells, one can treat this grid as a monochrome image and apply any morphological operation [54] onto it. Figure 3.7 illustrates an example from part of an industrial site. Directly extracting boundaries from the object patches produces numer- ous artifacts, shown in (b). Hence, an opening morphological operation is applied to removetheinsignificantfeatures,shownin(c). 24 (a) (b) (c) Figure 3.7: Boundary extraction for part of an industrial site: (b) without any mor- phological operation; (c) with an opening morphological operation, most artifacts are removed. Configuration: The only parameter in the boundary extraction algorithm is the unit lengthoftheuniformgrid. Thisparametercontrolstheconnectivityofthepatch. Similar to the ball radius parameter in ball-pivoting algorithm used by [60], a large cell length creates a coarse boundary, while a small cell length keeps details along the boundary as well as noise. Empirically, a cell length that equals to the neighborhood size δ is adopted. 3.4.3 BuildingModelCreation So far the building point cloud has been split into several planar roof patches with their boundaries marked out. A Na¨ ıve modeling approach takes each boundary as a roof polygon, and connects it to the ground by adding vertical walls. Solid building blocks are created in this way, which are combined together to form 3D building models. A min-maxpolygonsimplificationalgorithm[32]isintroducedtosimplifytheboundaries andreducethecomplexityofeachbuildingblock. Inpractice,duetotheresidualsensor noise, boundaries of neighboring planar roof patches are not strictly aligned together, 25 thus the Na¨ ıve approach may generate overlaps and cracks between neighboring build- ingblocks. Mystudyobservesthefactthatmostboundarylinesegmentsinalocalareafallintoa coupleofdirections,knownasprincipaldirections. Oncedetected,theyareusedtoalign boundaries before polygon simplification. In addition, if two segments from different roofpatchesarealignedtothesamedirectionandthedistancebetweentheirprojections on x-y plane is small, they are very likely to be the common boundary segment shared by the two roof patches (when their heights are close), or be the two ends of a vertical wall(whenthereisaheightgap). Basedontheseobservations,afour-stepsalgorithmis designedtoimprovethemodelingquality. Principaldirectiondetection Intuitively, the principal directions are the boundary edge directions which most commonly appear in the local area. Thus, the detection algorithm estimates the tangent direction at each boundary point using a 2D covariance analysis, and then employs a histogram to statistically analyze these directions. Principal directions are detected as thepeaksofthisdirectionhistogram. Forrobustness,aGaussianfilterisappliedtosmooththehistogrambeforedetecting peaks. Specifically,thehistogramismodifiedbyconvolutionwithaGaussianfunction g(x) = 1 √ 2π·σ ·e − x 2 2 2 . (3.10) The local maxima with sufficient samples of this smoothed histogram are detected aspeaks, i.e.,principaldirections. Figure 3.8(b) shows the direction histogram of a Denver city piece shown in Fig- ure3.8(a). Fourprincipaldirectionsaredetected: 0 ◦ ,43 ◦ ,90 ◦ ,133 ◦ ,whichcanbesepa- rated into two orthogonal pairs, illustrated as the colored arrows in Figure 3.8(a). This 26 0 0.5 1 1.5 2 2.5 3 0 50 100 150 (a) Input LiDAR data (b) Histogram of boundary edge directions (c) Modeling with principal direction snapping 1 2 3 4 Figure 3.8: Building modeling using principal direction snapping: (a) input LiDAR point cloud chopped from Denver city; (b) histogram of boundary line directions, with four principal directions detected and marked in the figure; (c) boundaries of planar roof patches are snapped onto the principal directions and simple building models are created. phenomenon concurs with the observation made by the previous work that orthogonal corners are a common pattern in building footprints. In addition, 0 ◦ and 90 ◦ directions contain less samples than the other two. They represents the chopping boundary of the area. Principaldirectionsnapping During this procedure, the boundaries of planar roof patches are snapped to the principal directions as much as possible without exceeding a small error toleranceϵ. A greedyalgorithmisappliedovereachpatchboundaryloopB,detailedasfollows. Everyboundarypointandprincipaldirectionpair(b i ,d j )definesalinein2Dspace: L(b i ,d j ) = {p|(p−b i )×d j = 0}. For each such pair, starting from the seed point b i , the longest continuous point sequence is detected alongB, in which all the points 27 have a distance toL smaller than an error toleranceϵ. I.e., in this continuous boundary segmentS(b i ,d j ),allpointssatisfy: distance(p,L(b i ,d j )) = |(p−b i )×d j | |d j | <ε. (3.11) The snapping algorithm then iteratively does the followings: (1) detect the point- directionpair(b max i ,d max j )whichmaximizesthesegmentlength|S(b i ,d j )|;(2)project (snap)thepointsinS(b max i ,d max j )ontolineL(b max i ,d max j );and(3)removeintermedi- ate points in S(b max i ,d max j ) fromB. The last step prevents overlaps between different segments. An exception is the end points of segments, which can belong to two bound- arysegments. Theiterativeprocessstopswhenthesizeofthelongestboundarysegment fallsunderacertainvalue(empirically,10points). Neighborsegmentssnapping After aligning segments with principal directions, the modeling results can be fur- ther improved by snapping neighbor segments to form vertical walls and avoid cracks. TwoboundarysegmentsL(b 1 ,d)andL(b 2 ,d)areneighborswhentheyhavethesame directiondandtheirdistanceissmallerthanagiventhreshold, i.e., distance(L(b 1 ,d),L(b 2 ,d)) = |(b 1 −b 2 )×d| |d| <ϵ l . (3.12) Foreachsuchneighborboundarysegmentpairs. • If they are from the same planar roof patch, and point toward the same direction, these two segments are snapped onto the same line. This operation handles the casethatalinesegmentisbrokenintoseveralsegmentsduetonoiseorocclusion. • Iftheyarefromdifferentplanarroofpatches,andpointtowardoppositedirections (i.e., are of the opposite orientation), the snapping algorithm first checks if they 28 overlapwhensnappedtothesameline. Ifso,theymustbethecommonboundary segment shared by the two roof patches (when their heights are close), or be the two ends of a vertical wall (when there is a height gap). Thus, they are snapped togethertocreatecorrespondingfeatures. Polygonalmeshcreation After boundary edge segments are detected and snapped to the principal directions, a polygon is created for each planar roof patch. Along the boundary loop, if two con- sequent segments share a common boundary point, the polygon simply connects them together. Otherwise, themin-maxpolygonalsimplificationalgorithm[32]isadoptedto simplifytheintermediatepointsequencebetweenthetwosegments. As mentioned previously, when all these polygons are created, vertical walls are generated to connect two snapped polygons, or a roof polygon and ground. The final modelisthecompositionofallthebuildingblocks. Figure3.8demonstratesawholepipelinetocreateabuildingmodel. Startingfroma buildingroofpatch,planepatchesaredetectedandboundariesareextracted. Thenthese boundaries are snapped to the principal directions which are detected in Figure 3.8(b). In addition, neighboring segments are aligned together to form vertical walls. The final result shown in Figure 3.8(c) contains 72 triangles, which is simplified from a building patchwith9,778roofpoints. 3.4.4 Non-PlanarRoofExtension Onelimitationofthisautomaticbuildingmodelingalgorithmisthatitonlysupportsflat roofs. However, in some cases, there are non-planar objects in the data set which need specialtreatments. Forexample,Figure3.9(a)showsanindustrialsiteinwhichthemost interesting objects are the oil tins and tanks represented in cylinder and cone shapes. 29 (a) (b) (c) Figure 3.9: Reconstruction of an industrial site using non-planar roof extension: (a) input LiDAR data; (b) detected object patches are rendered with bright colors, while ground and noise are rendered with dark-grey and black; (c) the reconstructed models, withdifferenttypesofobjectsgeneratedindifferentways. These shapes increase the complexity of the object reconstruction problem. However, fortunately,thesenon-planarshapesusuallyfallintoacoupleofspecificpatterns. There- fore, once the pattern types are known, typical pattern recognition algorithms such as RANSAC[14]canbeappliedtoestimateparametersandcreate3Dmodels. In particular, user interaction is introduced to identify the pattern types. The initial LiDAR data is first segmented into different patches and the noise/tree/ground points are removed from the data set. Then the object patches are exhibited to the user as Fig- ure3.9(b). Theusercanmovehismouseontoanyobjectpatch,clicktoselectthispatch, andpressakeytospecifyitspatterntype,e.g.,cone,cylinder,orplanar-roofobject. The modeling program automatically extracts boundary loops for each object, as presented in the previous sections; and a RANSAC algorithm [14] based on the identified pattern 30 type is applied to estimate the parameters of the model. The reconstruction result is showninFigure3.9(c). Objectsofdifferentpatternsarecreatedsuccessfully. 3.5 TerrainModeling The last part of the general urban modeling pipeline is terrain modeling, which is achieved by rasterizing the detected terrain point patch into a DEM (Digital Elevation Model). Specifically, all the terrain points are embedded into a uniform grid spanning over the entire urban area. One mesh vertex is created at the center of each grid cell with its height being the average height of all terrain points in this grid cell. Due to the occlusionbyovergroundobjectssuchasbuildings,trees,andvehicles,theuniformgrid may have empty grid cells at the occluded areas. Thus, the terrain modeling module fills holes by solving a Laplace’s equation∇ 2 z = 0 taking the heights in all empty grid cells as unknowns and the heights in all non-empty cells as the boundary condition. In particular,foreachemptycellat(i,j), 4z i;j = z i−1;j +z i+1;j +z i;j−1 +z i;j+1 . (3.13) The unknown heights are calculated by solving this equation array, and the ground meshisaDEMconnectingallthesevertices. 3.6 ExperimentalResultsonSmallSampleDataSets The general urban modeling system is tested on three different data sets: Denver, Oak- land,andanindustrialsite. Theresolutionofthesedatasetsvariesfrom6samples/m 2 to 17 samples/m 2 In order to process all the data in-core, a 1km-by-1km piece is chopped from Denver as shown in Figure 3.10(a) and a 600m-by-600m piece is chopped from 31 d e f (d) (d) (e) (e) (f) (f) (a) (b) (c) Figure3.10: Automaticurbanmodelingofa1km×1kmpiecefromthecityofDenver: (a) input LiDAR data; (b,c) modeling result viewed from different perspectives; (d,e,f) closeupsofthemodelingresultwithinitialpointcloudoverlayed Oakland as shown in Figure 3.11(a). All the experiments are made on a consumer- level laptop (Intel Core2 1.83GHz CPU, 2G memory) with an external hard disk. The time cost of the entire system is proportional to the number of input points, with a ratio around8minutes/millionpoints. Figure 3.10 shows the reconstruction results of an urban area in the city of Denver. Theresultingbuildingmodelsarecomposedofsimpleandcleantriangularmeshes,thus are ideal to some applications such as virtual city tourism. In addition, building blocks are aligned tightly together and form building models of good shapes. As shown in Figure 3.10(d,e,f), these building models fit the raw LiDAR point cloud in an excellent manner. 32 0 0.5 1 1.5 2 2.5 3 0 100 200 300 400 1 2 3 4 5 6 7 (a) (b) (c) Figure 3.11: Urban modeling of an area in the city of Oakland: (a) input LiDAR point cloud, with intensity representing point height; (b) a histogram with seven principal directions detected and illustrated as colored arrows in (a); (c) reconstructed building models;bothorthogonalcornersandnon-orthogonalcornersarereconstructedcorrectly. To further demonstrate the ability of supporting multiple principal directions, the systemistestedonachallengingareafromthecityofOakland,asshowninFigure3.11. Seven principal directions are automatically detected from the original data sets, which are illustrated as the arrows in Figure 3.11(a). Since the general modeling system has no pre-assumptions on the corner angles, correct models are generated with corners of anglesbetweenanytwoprincipaldirections,whichmaybe90 ◦ ,45 ◦ ,oranyotherangle; showninFigure3.11(c). The final testing data set is an industrial site shown in Figure 3.9(a). This data set is used to demonstrate the non-planar roof extension to the urban modeling system. Besides plane-shaped roofs, two new types of geometry patterns are involved, namely, a standing cylinder with a cone on top of it (an oil tin) and a lying cylinder (a tank). 33 Objectsofallthreepatternsaresuccessfullydetectedandreconstructedbytheproposed approach,asshowninFigure3.9(c). 34 Chapter4 StreamingUrbanModelingfor Large-ScaleDataSets Likemanypreviouswork,theautomaticurbanmodelingpipelineproposedinChapter3 isanin-coremethod. Inotherwords,itneedstoloadalltheLiDARdataintothememory before processing. Thus, there is a conflict between the increasing size of data sets and thelimitationofcomputerhardware. Acommonwaytoalleviatethisproblemistopartitionthewholedatasetintosmall tiles and process them one at a time (e.g. [41, 48]). By merging building models gen- erated from different tiles, this approach can produce 3D models from large LiDAR data sets. However, it may introduce artifacts alongside the boundaries between tiles. Although these boundary effects can be moderated by introducing extra processing on boundaryregions,theadditionalprocessingfortilepartitioning,boundaryhandlingand modeling merging is inefficient, tedious and may introduce ambiguity (e.g. large build- ingsthatspantheintersectionsofmultipletiles). This chapter presents a streaming framework to handle extremely large data sets in a seamless manner, meaning that the proposed method needs no special treatment for tile-boundaries. Bystoringdataasstreamfilesonharddisksandusingmainmemoryas onlyatemporarystorageforongoingcomputation,themethodachievesefficientout-of- core data management. This gives the urban modeling system the ability to handle data setswithhundredsofmillionsofpointsinauniformmanner. Forexample,byadapting the automatic urban modeling pipeline into this streaming framework, the whole urban 35 modelofAtlantaisreconstructedfrom17.7GBLiDARdatawith683Mpointsbyunder 25 hours of unattended processing using less than 1GB memory. As a comparison, an in-coreprogramwouldneedmorethan100GBmemorytoprocessthisdatainonepass. Section 4.1 briefly reviews the state-of-the-art streaming approaches. Section 4.2 presentsastreamingframework. InSection4.3thegeneralurbanmodelingsystempro- posed in Chapter 3 is adapted into the streaming framework. Experiments on different city-scaledatasetsaregiveninSection4.4. 4.1 ReviewofStreamingApproaches To solve the conflict between extremely large data sets and computer hardware limi- tation, streaming methods are developed in geometry modeling and computer graphics areas. They have succeeded in a board variety of applications, such as mesh process- ing [26] and compression [25], tetrahedral mesh simplification [61], level sets meth- ods [45], point cloud processing [46], LiDAR data rasterization [28], dynamic process- ing[35],Poissonequationsolver[53],anddelaunaytriangulation[27]. Here I highlight [27], [28] and [46]. Isenburg et al. [27, 28] reveal the local spatial coherence ofaerial LiDAR data andpropose agrid-based indexingstructureand a spa- tialfinalizationmechanism,whichisthebasisoftheapproachproposedinthischapter. Pajarola [46] performs a sequence of operations on a data stream. To allow data blocks tobeindifferentstatesanddealwithtransitionsbetweenstates,hearrangesdatapoints into a FIFO queue by sorting the data along one dimension of the largest extent, which islessefficientandgeneral,comparedwiththestatepropagationmechanismforsolving thesameproblem,asstatedinSection4.2. 36 4.2 StreamingFramework Generally, the streaming framework is an out-of-core architecture. A streaming pro- gram acts as a frontier through the data stream - it reads from an input stream, loads necessarydatain-core;processesdata;andoncethedataisnolongerneededforfurther processing, it is written into the output stream. To adapt the general modeling pipeline into this out-of-core architecture is particularly difficult, because each urban modeling module in the general pipeline involves several intermediate steps. Each step usually reliesononepoint’sneighboringpointsbeingprocessedbysomeprecedingsteps. E.g., in classification, normal covariance analysis over a pointp requires position covariance analysis results on all p’s neighboring points. Thus, these intermediate steps cannot be batch-processed,andexecutionorderofthesestepsmustbecarefullydetermined. Myresearchachievesthisgoalbydefining streamingoperatorsand steamingstates formally in Section 4.2.2, and introducing a state propagation algorithm to determine the execution order in Section 4.2.3. Section 4.2.1 gives the definition of point stream andreviewsthe finalizationmechanismasfirstproposedbyIsenburg et al.[27]. 4.2.1 PointStreamsandFinalizer As observed by Isenburg et al. [27] and Pajarola [46], spatial coherence, which either appears in the original data set or is the result of a resorting algorithm, can greatly improve the memory efficiency in an out-of-core algorithm. To exploit such spatial coherence,pointstreamisdefinedasthebasicformofdatathatisprocessedinastream- ingframework: Definition4.1 A point stream is a FIFO queue composed of point records and finaliza- tion tags. 37 Late Early Finalization time (a) (b) Figure4.1: FinalizationresultoftheOaklanddatasetwithcellcolorsin(b)representing thefinalizationtime. Spatialcoherenceisrevealedbytheregularityofcolordistribution. A point stream is generated by inserting finalization tags into original input data, which is done by the pre-processing module called Finalizer. A finalization tagf A is a symboltoindicatethatallthepointrecordsinaspatialareaAhasappearedinthepoint stream before it. This is necessary because the original data usually is not spatially ordered strictly. When a streaming program meets f A , it gets the guarantee that the informationwithinAisavailableandfurtheractionscanbetaken. Inparticular,theinputdataispartitionedinto2 k ×2 k uniformrectanglegrid,sothat each cell is the basic unit of spatial area in streaming processing. Therefore, the task of the Finalizer istoinsertonefinalizationtagforeachsuchcellintotheinputdata. Note that the finalization tags only providea mechanism to reveal the spatial coher- encebutnottogenerateit. Forexample,intheworstcase,thelastpointofeachgridcell appears at the very end of the input data; the finalization tags will then all be inserted at the end of the stream and no memory efficiency can be produced. Fortunately, the point sequence in the aerial LiDAR data shows remarkable spatial coherence to make significantmemoryefficiency[27]. Also,thespatialcoherencecanbefurtherenhanced bya chunkingalgorithmwhichpartiallyresortthepointrecords[27]. 38 Figure4.1showsthefinalizationresultoftheOaklanddataset. Thecolorsillustrate thetimewhenagridcellisfinalizedinthepointstream. Thespatialcoherencebetween gridcellsisrevealedbytheregularityofthecolordistribution. 4.2.2 StreamingOperatorsandStates The basis of the point stream processing framework is the following observation: most ofthecomplicatedlocalalgorithmscanbedecomposedintoaseriesofstreamingoper- ators,whichisageneralizedformofthe“streamoperators”definedin[46]. Definition 4.2 A streaming operator k (p i ) is a local operator which requires all the points inp i ’s local neighborhood N k (p i ) to be at streaming state s k−1 or higher state; k (p i )takesthesepointsasinputandtransitpointp i to streaming states k . Hereaseriesofstreamingstatess 0 ,s 1 ,...,s m isdefinedasfollows: thefirststates 0 is always “Unread” and the last states m is always “Written and released”. Except for the last state, the information in a point record at a lower state is always a subset of the information at a higher state. In a complete streaming process, each point sequentially experiencesstatesfroms 0 tos m . Astatetransitionfroms k−1 tos k canonlybeinvoked byastreamingoperator k . The streaming classification module is demonstrated here as an example. All the points are in the initial state s 0 =“Unread”. The first streaming operator 1 (p i ) reads p i from the point stream into memory and changes its state from s 0 to s 1 =“Read”. The second operator 2 (p i ) collects the positions of points in p i ’s local neighbor- hood N 2 (p i ), uses these information to estimatep i ’s normal, then transitsp i into state s 2 =“With normals”. During this process, all the points in N 2 (p i ) must be at least at streaming state s 1 , i.e., read from the stream and loaded in-core. In a similar manner, the following operators are executed sequentially until point p i finally gets into state 39 s m−1 . The last stream operator m (p i ) then writes it to the output stream, releases it from memory and turns its state into s m =“Written and released”. In this whole process, point records which are in the first and last states are stored in disk files, and only a small fraction of points in intermediate states need to be loaded in-core. These intermediatestatesarecalled active states. Todeterminewhenanoperator k canbeinvoked, scope radiusisdefined: Definition4.3ThescoperadiusR( k )istheradiusofpointp i ’sneighborhoodN k (p i ) requiredby k (p i ). R( k ) reflects the size of the area affecting k . For convenience, let R( k ) = 0 whenthestreamingoperator k doesnotneedinformationfromthelocalneighborhood, e.g.,a“readfromstream”operator. Inothercases,itisdeterminedbythecorresponding streaming operator. E.g., the “normal estimation” operatorofthe classificationmodule requiresascoperadiusequaltotheneighborhoodsizeδdefinedinSection3.2. Theonly exception is the scope radius of the last operator m , which is forced to be the largest scoperadiusofalltheotherstreamingoperators, i.e., R( m ) = max{R( 1 ),R( 2 ),...,R( m−1 )}, (4.1) because m is the only information-subtracting operator – once performed, the point record is no longer available in the memory. By forcingR( m ) to be the largest scope radius, m (p i ) is applied only when all the points in N m (p i ) are at least at state s m−1 (written or waiting to be written), so that no point still requires information fromp i to completeastreamingoperator k ,k < m. The scope radius is also helpful in determining the size of the spatial unit (cell). The side length of a grid cell must be no less than R( m ), so that the impact of any 40 streaming operator applied on a cellc is restricted within its 1-ring neighborhood. This isparticularlyconvenientforthestatepropagationalgorithminthefollowingsection. 4.2.3 StatePropagation Statepropagationisanalgorithmwhichperformsthestreamingoperatorsinthecorrect order. State propagation uses cell as the basic unit to perform an operator, and operator k performed on cell c i;j is denoted by k (c i;j ). As discussed before, to determine whether k (c i;j ) can be invoked, one only needs to check if all points in c i;j ’s 1-ring neighborhood are at least in state s k−1 . In addition, if a cell reaches a new state s k−1 by operator k−1 , only its 1-ring neighbor cells may receivethe direct impact from this transition, e.g., a neighbor cell may now satisfy the state prerequisite for k . Based on thisobservation,thekeyideaofthestatepropagationalgorithmistonotifyallthe1-ring neighborswheneveracell’sstateischangedbycompletinganewoperator. The algorithm starts by reading a cell c next from the input point stream S in . A recursivefunctioncellAction()isthencalledtoperformsteamingoperatorsinanorderly manner. Taking a cell c and a streaming operator k as input, cellAction() first checks if the state prerequisite for operator k (c) is satisfied. If not, it aborts the operation; otherwise,itperformsoperator k (c),transitcellctostates k andnotifieseachcellc ∗ in c’s1-ringneighborhood by recursivelycalling cellAction()forc ∗ andoperator k+1 . In this way, the state change of the initial cellc next is propagated in the grid and operators will be performed once they are ready. The pseudo-code of the algorithm is shown in Table4.1. Figure 4.2 shows an example of a state propagation procedure in a 4-states prob- lem. The state of each cell is denoted by the number and the color of the cell. In the first step, a cell c 1;1 (marked by the black frame) is read from the point stream, and its 41 /*************************** Mainprogram ***************************/ Input: apointstreamS in ;asetofstreamingstates{s 0 ,...,s m };andasetofstreaming operators{ 1 ,..., m }. Output: apointstreamS out . WhileS in isnotemptydo: • Readthenextcellc next fromS in ; • Callfunction cellAction(c next , 1 ). Endofmainprogram. /*********************** cellAction()function ***********************/ Input: agridcellcatstates k−1 ;andastreamingoperator k . Foreachcellc ∗ inthe1-ringneighborhoodofcdo: • ifthestateofc ∗ islowerthans k−1 ,thenreturn “not ready”. /*Ifnotreturned,allcellsinthe1-ringneighborhoodpassthestatetest. */ Execute k (c). /*takeactionandchangethestateofctos k */ Foreachcellc ∗ inthe1-ringneighborhoodofcdo: • ifthestateofc ∗ iss k ,thencallfunction cellAction(c ∗ , k+1 ). EndofcellAction()function. Table4.1: Statepropagationalgorithm 0 0 0 1 0 1 1 0 1 1 2 0 1 2 2 1 2 2 3 1 1 2 3 3 0 0 0 1 0 1 1 1 0 1 2 0 1 2 2 1 2 2 3 1 1 2 3 3 0 0 0 1 0 1 1 1 0 1 2 2 0 1 2 1 2 2 3 1 1 2 3 3 0 0 0 1 0 1 1 1 0 1 2 2 0 1 2 3 1 2 2 3 1 1 2 3 3 0 1 1 2 2 3 Figure 4.2: An example of the state propagation algorithm. The numbers and colors denote states. When the state of a cell is changed (marked with black frame), it notifies its1-ringneighborhood(redframe)tocheckifanyofthemisreadyforthenextoperator. Thewholeprocessisarecursiveprocedure. 42 state is transited to s 1 via a call to cell-Action() function, which then leads a number of streaming operators performed on other cells and state updates. Cells that reach the final state will be written to the output file. The propagation terminates when no more streamingoperatorispossible(i.e.,noneofthemfulfilltheirprerequisites),andthestate propagationalgorithmwillreadanewcellintotheactivesetandrepeatthepropagation untilallinputdataisfullyprocessed. 4.3 StreamingUrbanModeling This section shows how to adapt the general urban modeling system into this stream- ing framework. Section 4.3.1 presents the modified pipeline. The subsequent sections thendetailmodulesincludingstreamingclassification,streamingsegmentation,building modeling,andterrainmodeling. 4.3.1 StreamingModelingPipeline AnoverviewofthestreamingpipelineisdemonstratedinFigure4.3. TheinputLiDAR data(usuallystoredinalistofdiskfiles)issequentiallyreadbyapre-processingmodule calledFinalizer,whichinsertsfinalizationtagsintothedata,andproducesapointstream asdefinedpreviously. Taking this point stream as input, two specific streaming modules are performed sequentially: the Classifier whichclassifiesvegetationpointsfrombuildingandground points, and the Splitter which segments single building patches from the building and ground points. Both are implemented following the formulation of streaming operators andstatesinthelastsection. IntheClassifierandSplitterblocksofFigure4.3,thesolid colored cells denote active (i.e., in-core) data set with different colors corresponding to different streaming states; the dark red region denotes processed and released data; and 43 Input LiDAR Finalizer Classifier Splitter Points Points & spatial fin-tags Classified points & spatial fin-tags Modeler Boundary extraction Terrain Generator Urban model Ground points & spatial fin-tags Plane fitting Modeling Building points of each patch Building models Terrain Figure 4.3: An illustration of the streaming building reconstruction pipeline. A pre- processing module (which is called Finalizer in [27]) inserts finalization tags (yellow ovals) and generates a point stream which flows over the Classifier and the Splitter sequentially. Both components introduce a state-propagation mechanism so that only data with active states (solid colored cells in Classifier and Splitter) are loaded in-core. TheSplitter finallygeneratesapointstreamandabuildingstream;whichareconverted into a terrain model and various building models using the Terrain Generator and the Modeler respectively. thedarkbluedenotestheinputwaitingtoberead. Theactivesetprogressesasafrontier throughtheinputstreamuntilalldataisprocessed. The Splitter outputs two streams: a point stream with geometry information of all ground points and a building stream which consists of individual building patches. The pointstreamisconvertedtoaterrainmodelusingaTerrainGenerator;andaModeleris responsibleforturningeachbuildingpatchintoapolygonalbuildingmodel. Thewhole urbanmodelisfinallycreatedbycombiningthemtogether. 44 Streamingoperators Streamingstates 1 : Readdatafrominputpointstreamandallocate memory. s 0 : Unread s 1 : Read 2 : Applycovarianceanalysisonpositionsto estimatenormals;andcomputeF 1;2;3 . s 2 : Withnormals 3 : Applycovarianceanalysisonnormals;calculate F 4;5 ;andapplySVMclassifier. s 3 : Classified 4 : Refineclassificationbymakingpointsinlocal neighborhoodvoteonresult. s 4 : Refined 5 : Writepointrecordstooutputpointstream,and releasethemfrommemory. s 5 : Writtenandreleased Table4.2: Streamingoperatorsandstatesforclassification 4.3.2 StreamingClassification ThefirststepoftheclassificationalgorithmproposedinSection3.2isaSupportVector Machine training procedure [2, 8] based on local geometric features. This off-line step isprocessedonceandforall,thusneedsnottobeadaptedintothestreamingframework. With the training results, the classification algorithm introduces a linear classifier with five types of features based on differential geometry properties: regularity F 1 , horizontalityF 2 ,flatnessF 3 ,andtwonormaldistributionmeasurementsF 4 andF 5 . To computethesefeaturesforapointp,covarianceanalysisisperformedtwiceonitslocal neighborhood. The trained linear classifier then computes the classification result ofp from these features. Finally, the classification results on p’s neighbor points will vote forthefinallabelofpinarefinementstep. Sinceeachpartofthisalgorithmrequiresonlytheinformationwithinalocalneigh- borhood of point p, it is easily decomposed into a series of streaming operators and states shown in Table 4.2, which are placed into the streaming framework and form up the Classifier module. 45 Hash table Hash table Hash table Set-Union Set-Union Flatten trees c c c c c c Figure 4.4: The streaming agglomerative clustering algorithm stores set-tree roots (marked with red frame) in a global hash table. It first performs a union operation on each neighbor point pair illustrated as the dotted line in the left figure. Then the algorithm flattens the set-trees shown in the middle figure and push new roots into the hashtable(rightfigure). 4.3.3 StreamingSegmentation The segmentation module aims to divide the input points into a ground patch and indi- vidual building patches. When adapting the general segmentation approach (Section 3.3) into the streaming framework, the main difficulty lies in the global characteristic of the segmentation problem. Specifically, patches (e.g., the ground patch) may span over the entire area, and thus one cannot throw away all the information when data is written into the output stream. Instead, a small-but-sufficient global indexing structure is required, which is never released from memory to make the segmentation algorithm executeinacorrectmanner. IchoosetoextendtheagglomerativeclusteringalgorithmproposedinSection3.3.2 into a streaming segmentation approach, since it is a bottom-up algorithm and most of the structural information is stored as parent pointers attached to each point. The remainingstructuralinformationwhichneedstobeprocessedgloballyisthetemporary 46 roots of each trees, which are stored in a hash table as the global indexing structure. HereIcalltheseroots“temporary”becausetheymaybemergedtothesametreeinlater processing. However, in order to avoid dangling pointers, once a point is pushed into thehashtable,itisneverreleasedfromthememory. Even with this global indexing structure, there is still a danger of visiting dangling pointers. During streaming process, some points may be outputted to disk and released from the active set. Thus, if a point’s root-seeking path contains some released points, invalid pointers will be visited when retrieving the cluster of this point. To solve this problem, find operation is introduced to flatten the root-seeking path as proposed in Section 3.3.2. The idea is straightforward: every time the agglomerative clustering algorithm in a cellc finishes, a find operator is applied on all points inc’s 1-ring neigh- borhood, to flatten the root-seeking path of these points; after this, the newly created roots are pushed into the hash table. Thus, every point that has been touched in this processisguaranteedtohaveaparentpointerpointingtoarootintheglobalhashtable, whichwillnotbereleased. Thepseudo-codeofthisalgorithmisshowninTable4.3. Figure4.4showsanexam- ple that illustrates the streaming agglomerative clustering process. The algorithm starts withacellcwhose1-ringneighborhoodareallavailableinmemory. Thepairsofpoints involved in c whose distances are small enough for a merging operation are connected indashedlines. Toprocessc,aunionoperationisperformedoneachsuchpairofpoints andtheirsegmentsaremerged(middlefigure). Thealgorithmthenperformsafindoper- ation over all the points touched in the first step to flatten all the trees which have been changed. Finally, the new roots generated during this process are added into the hash table. 47 Input: a cell c at state s 1 , with its 1-ring neighborhood at state s 1 or higher; root hash tableH;andthedistancethresholdα. /*Apply unionoperationovercellc*/ Foreachpointpair (p,q)wherep∈cand||p−q||<α,do: • Callfunction union(p,q). /*Flattensettreesfromalltouchedpoints*/ Foreachpointpinthe1-ringneighborhoodofcdo: • Flattenp’sroot-seekingpathbycalling find(p). /*Putnewrootsintothehashtable*/ Foreachpointpinthe1-ringneighborhoodofcdo: • ifpisrootandp̸∈H,thenpushpintoH. Endofunion-findalgorithm. Table4.3: Streamingunion-findalgorithm( 2 ) Streamingoperators Streamingstates 1 : Readdatafrominputpointstreamandallocate memory. s 0 : Unread s 1 : Read 2 : Applystreamingunion-findalgorithmdescribed inTable4.3. s 2 : Segmented 3 : Writepointrecordstooutputpointstream,and releasethemfrommemory. s 3 : Writtenandreleased Table4.4: Streamingoperatorsandstatesforsegmentation Except the global hash table, this segmentation method is completely local, thus it canbedefinedasstreamingoperators. Alistofthestreamingoperatorsandcorrespond- ingstatesaregiveninTable4.4. As a result, the output point stream is decomposed into segments of points. The largest segment is taken as the ground patch and sent to the Ground Generator still in theformofapointstream. Therestsegmentsaresenttothebuildingmodelingmodule. 48 0 0.5 1 1.5 2 2.5 3 0 500 1000 1500 2000 2500 0 0.5 1 1.5 2 2.5 3 0 500 1000 1500 2000 2500 Figure4.5: WiththeprincipaldirectiongridontheOaklanddataset,sixprincipaldirec- tions are detected for the blue cell (left), while two principal directions are detected for theorangecell(right). 4.3.4 BuildingModeling Thebuildingmodelingalgorithmnowtakesoverthesebuildingpatches. Sincethenum- ber of points contained in a single building patch is small, the patches are loaded into the memory and processed one by one. In the experiments, the largest buildingpatch is thelargestructure showninFigure4.6(a),containing3.2Mpoints, whichtakes332MB ofmemorytoprocess. This section extends the automatic building modeling algorithm in Section 3.4 for building model reconstruction. Given a building patch and a set of principal directions asinput,thealgorithmautomaticallyfitsplanestothepointsandsnapstheplanebound- arysegmentsontotheprincipaldirections. 49 a b c d (a) (a) (b) (b) (c) (c) (d) (d) Figure 4.6: Reconstructed urban model of Atlanta city. Closeups of different areas are shownintherightsub-figures. In Section 3.4.3, the principal directions are extracted over the entire input data, which is problematic in city-scale input. To allow the principal directions to reflect different boundary directions within local regions (such as the downtown area in the cityofAtlantashowninFigure4.6(b)andtheresidentialareashowninFigure4.6(d)),a principaldirectiongrid(Figure4.5)isintroduced. Foreachcellinthisgrid,ahistogram ofthetangentdirectionsofallboundarypointsiscomputedwithinalocalneighborhood andthepeaksafterGaussianfilteringarefoundtobeprincipaldirections. 4.3.5 TerrainModeling The objective of the terrain modeling algorithm is to rasterise the ground point stream into a digital elevation model. Taking the point stream as input, the algorithm builds up a square grid (whose user-selected unit length determines the precision of the ter- rain mesh); and counts the lowest ground point in each grid cell. These points are later acceptedasverticesoftherasterisedterrainmodel. Theemptycellscanbefilledbysolv- ing a Laplace’s equation as proposed in Section 3.5. However, for performance reason, alinearinterpolationisappliedtothegapsalongeachcolumnonthegrid. Trivialdiffer- encesexistbetweenresultsgeneratedbythesetwomethods;however,theimprovement onefficiencyisremarkableforcity-scaledatasets. 50 0 1 2 3 4 5 6 7 8 x 10 4 0 200 400 600 800 1000 Time (s) Memory usage (MB) Finalizer Classifier Splitter Other components Figure4.7: MemoryusageduringexperimentonAtlantadataset 4.4 ExperimentalResultsonCity-ScaleDataSets Thestreamingurbanmodelingsystemistestedonthreedifferentdatasets,namely,Oak- land, Denver, and Atlanta. The problem scale varies from 16M points to 683M points. My program generates polygonal urban models for each of them. All the experiments are done on a consumer-level desktop PC (Intel Core2 2.4GHz CPU, with 2G mem- ory and 100GB free hard disk space). The running time and maximum memory usage are reported in Table 4.5. Although the average performance is affected by the charac- teristics of data sets, the average processing speed is faster than 3 minutes per million points. Benefiting from the streaming framework, the memory footprint during the experi- ments is kept at a low level. The whole pipeline consumes no more than 1GB memory at any time to process the largest data sets (Atlanta) in one pass. This memory-saving mechanism can be explained by Figure 4.10, showing the Classifier module for Oak- landdataset. First,thefinalizationresultinFigure4.10(b)revealsthespatialcoherence between streaming grid cells. Second, three snapshots are taken during the streaming 51 Model Oakland Denver Atlanta InputLiDAR data Urbanarea 1.2km-by- 0.8km 4km-by-3km 5.5km-by- 7.1km Point number 16M 73M 683M Filesize 437MB 1.90GB 17.7GB Grid resolution 64×64×64 512×512× 512 512×512× 512 Time (hh:mm:ss) Finalizer 24 4:43 34:58 Classifier 11:10 53:31 9:18:09 Splitter 3:14 2:03:56 11:54:29 Others 3:09 15:23 2:33:24 Total 17:57 3:17:33 24:21:00 Maximum memoryusage (MB) Finalizer 108 16 209 Classifier 276 156 888 Splitter 103 72 390 Others 23 71 332 Maximum 276 156 888 Output triangles Building 62K 182K 1.12M Terrain 1.92M 10.7M 8.78M Table 4.5: Three data sets with different sample rates are tested using the streaming building reconstruction system on a consumer-level PC. This table reports the running time and maximum memory usage in each pipeline module, i.e., Finalizer, Classifier, Splitter,andothercomponents. Theexperimentalresultsshowtheabilityofthestream- ingsystemtohandleextremelylargedatasetsinanefficientmanner. classification processing (Figure 4.10(c,d,e)), showing that only cells at active states (brightsolidcoloredcells)arestoredinmemory. Theremainingcellsareeitherwaiting intheinputstreamorhavebeenwrittentotheoutputstreamandreleasedfrommemory. Withthespatialcoherenceguaranteed,theactivecellsarealwaysasmallfractionofthe data;thusonlysmallamountofmemoryisrequired. 52 0 2000 4000 6000 8000 10000 0 50 100 150 200 Time (s) Memory usage (MB) Finalizer Classifier Splitter Other components a b c (a) (a) (b) (b) (c) (c) (d) (d) Figure 4.8: Urban model of Denver with three closeups shown in (a,b,c). Although principaldirectionsintheseareasaredifferent;withtheprincipaldirectiongrid,correct principal directions are generated for each of them. (d) Memory usage during process- ing. Figure4.6andFigure4.8demonstratethereconstructionresultsforAtlantaandDen- ver respectively. The memory usage in the process of reconstructing the urban model of Atlanta city is plotted in Figure 4.7. The Classifier is the most memory consuming module because it has more active states, thus stores more cells in memory; the Split- ter is the most time consuming module because of the overhead for saving segmented patches into files. For the Denver data set, its urban model and closeups are shown in Figure4.8. Bothdatasetsexhibitvariationofprincipaldirectionsacrossthewholecity. Nevertheless,itisnicelyhandledbythegrid-basedprincipaldirectionestimation. 53 Streaming Tiling Boundary artifacts Clean result Figure 4.9: Comparison between streaming method (left) and tiling method (right). Even with padding along tile boundaries and tiles batch processed one-by-one using an in-core modeling program (the blue square shown in the top-right sub-figure), artifacts canbegeneratedalongtileboundariesasshowninthebottom. 4.4.1 StreamingVersusTiling This section makes a comparison between the streaming method and the traditional tiling method [41, 48]. In particular, the tiling method is tested on Atlanta data set following the partition setup described in [41], i.e., 0.6km× 0.6km tiles with 0.2km padding at each edge. The tiling grid is shown in the top-right sub-figure of Figure 4.9, with the blue square representing an in-processing tile with its padding. Each tile is processed using the general urban modeling pipeline proposed in Section 3. Following observationsaremade: 54 1. Largebuildingstructuresthatarenotfullycapturedbyonesingletilestillexhibit boundary artifacts. These building structures include the largest building shown in Figure 4.6(a). Boundary artifacts are produced by the tiling method as shown inthebottom-rightsub-figureofFigure4.9. 2. Tiling disables global segmentation. In some tiles, ground is segmented into sev- eral pieces and some are incorrectly detected as building structures. E.g., with ramps connecting highways and local roads not in the same tile, the highway is often detected as an individual segment and reconstructed as a building; while with the streaming approach highways are appropriately detected as part of the ground. 3. Performanceissue: withsimilarurbanmodelingtechnique,thetilingmethodruns for over 60 hours mainly because of the additional overhead for padding areas (with padding, the data size grows to roughly 3 times to the original size); and requiresuserinteraction(e.g.,theuserneedstospecifywhichreconstructedbuild- ing structure is desired when the same building appears partially in neighboring tiles); while the streaming approach runs fully automatically for 24 hours on the samecomputer. Moreover, the streaming approach is especially useful under the current trend of rapidly increasing data resolution. E.g., the Atlanta data set is with 17 samples/m 2 resolution,comparedwith1sample/m 2 datasetin[41]and9samples/m 2 datasetin[60]. 55 (a) Early Late (b) (c) (d) (e) Figure 4.10: (a) Reconstructed urban model of Oakland downtown area. (b) Finaliza- tion result reveals the spatial coherence between cells; colors represent the finalization time. (c,d,e)Threesnapshotsduringthestreamingclassificationalgorithm;onlyasmall portionofcellsareatactivestates(brightsolidcoloredcells). 56 Chapter5 2.5DDualContouring Figure 5.1: Various kinds of building models are created using 2.5D dual contouring. Fromlefttoright: twostadiummodelswithdifferentkindsofnon-planarroofs;atypical flatbuildingfromresidentialarea;andamodernhighbuildingfromurbanarea. Thischapter focuseson the complicated problemofreconstructingbuildingmodels from aerial LiDAR point cloud. The aerial LiDAR point clouds are 2.5D data, i.e., the LiDAR sensor captures the details of roof surfaces, but collects few points on building walls connecting roof boundaries. In addition, manually created building models (Fig- ure 5.2) also show a 2.5D characteristic. Nearly all of them consist of complex roofs (green faces) connected by vertical walls (white faces). Thus, a 2.5D modeling method withthefollowingpropertiesisdesired: • Accuracy: Themethodshouldproducesimplepolygonalmodelsfittingtheinput pointcloudsinaprecisemanner. • Robustness: Regardless of the diversity and complexity of building roof shapes, the method should always generate crack-free models, even with the existence of undesiredelementssuchasresidualsensornoiseandsmallrooffeatures. 57 Figure 5.2: Manually created models in Google 3D warehouse [21], showing the 2.5D natureofbuildingstructuremodels • 2.5Dcharacteristic: Themethodshouldcreate2.5Dpolygonalmodelscomposed ofdetailedroofsandverticalwallsconnectingrooflayers. In Section 3.4.3, an automatic buildingreconstruction algorithm is developed based on planar roofs. Its extension supports non-planar roofs with the help of user interac- tions. However,bothmethodsarelimitedbythepre-definedroofpatterns,i.e.,planesor geometrypatternsinaprimitivelibrary. Thus,theycanhardlyproducebuildingmodels witharbitrarilyshapedroofs. This limitation is also a common disadvantage in previous research work. E.g., [41,48,60]createplanarbuildingroofsand[33,64,65]relyonasmallsetofuser-given primitives. Thesemethodsworkwellforbuildingscomposedofpre-definedshapes,but lose accuracy and robustness when dealing with arbitrary roof shapes such as those showninFigure5.1. Anotherwaytoattackthebuildingmodelingproblemiswithtraditionaldata-driven approaches. Polygonal models are first generated directly from input data using ras- terization or delaunay triangulation, then simplified with general mesh simplification algorithms. Thelatterstepsignificantlyreducestrianglenumberwhilepreservingalow fitting error. However, since the general simplification algorithms are usually “blind” to the 2.5D nature of the problem, they can hardly produce models satisfying the 2.5D requirement. 58 This chapter proposes a novel, data-driven approach to solve this problem, named 2.5Ddualcontouring. Liketheclassicdualcontouring[31],anadaptivegridisadopted as the supporting data structure. Geometry is reconstructed in each grid cell by min- imizing the quadratic error functions known as QEFs. Model simplification is easily achievedbymerginggridcellsandcombiningQEFs. In order to represent the detailed roof surfaces, 2.5D dual contouring works in a 3D space. However, unlike the classic 3D dual contouring, it uses a 2D grid as the sup- porting data structure. The 2.5D dual contouring approach generates a hyper-point in each grid cell, which contains a set of 3D points having the same x-y coordinates, but different z values. They can be regarded as a set of points intersected by a vertical line andmultiplerooflayers. Hence,theconsistencybetweenboundaryfootprintsofdiffer- entrooflayersisguaranteed,andverticalwallsareproducedbyconnectingneighboring hyper-pointstogether. The following part of this chapter is organized as follows: Section 5.1 reviews the classic dual contouring. Section 5.2 presents the 2.5D dual contouring pipeline, fol- lowed by sections describing the details of each step. Finally, experimental results are showninSection5.8. 5.1 BriefReviewofDualContouring Like many volumetric methods [10, 40], dual contouring [31] has proved to be a robust way of generating crack-free 3D models: input points are first scan-converted into a regularizedgrid;thengeometryandtopologyarecreatedrespectively. 59 Inparticular,thedualcontouringmethodtakesa3DgridofHermitedata(intheform ofpoint-normalpairs)asinput,andcreatesexactlyonemeshvertexineachminimalgrid cellbyoptimizingaquadraticerrorfunction(QEF),definedas: E(x) = ∑ i (n i ·(x−p i )) 2 , (5.1) where (p i ,n i ) are the Hermite data samples in this cell. Since the error function is defined based on both positions and normals at intersection points of the surface with grid edges, the optimized vertices have a trend to lie on the sharp features such as roof ridges and valleys. Based on this observation, Fiocco et al. [13] use classic 3D dual contouring to create 3D building models from both aerial and ground-based LiDAR, preservingsharpbuildingfeatures. The geometry simplification in dual contouring is achieved in an adaptive manner. Ju et al. [31] introduce an octree for 3D geometry simplification. They merge QEFs associatedwithleafnodeswhilecollapsingtheoctreestructure. Thesub-treecollapsing iscontrolledbytheresidualofthemergedQEF,whichisrequiredtobelessthanagiven tolerance. Finally,dualcontouringcreatespolygonsduringatraversalovertheadaptive grid. A topology safety examination is introduced to avoid possible topology changes duringsimplification. Nevertheless, the classic dual contouring approach is designed with regular 2D or 3Dgrids. Itdoesnotsatisfythe2.5Drequirementforbuildingmodels. 5.2 2.5DDualContouringPipeline Given a building point cloud as input, the 2.5D dual contouring modeling process exe- cutesfourstepsasillustratedinFigure5.3: 60 (a) (b) (c) (d) (e) Figure 5.3: Robust building modeling pipeline: (a) the input point cloud; (b) a 2D grid with surface Hermite data (gold arrows) and boundary Hermite data (red arrows) attached; (c) hyper-points (turquoise balls connected by red lines) generated by min- imizing QEFs; (d) mesh model reconstructed via 2.5D dual contouring; and (e) final modelwithboundariessnappedtoprincipaldirections. 1. Scan conversion: The point cloud is embedded in a uniform 2D grid. Surface Hermite data samples (gold arrows) are generated at grid points and boundary Hermitedatasamples(redarrows)areestimatedongridedgesconnectingdiffer- ent roof layers (Figure 5.3(b)). This 2D grid is also regarded as the finest level of thesupportingquadtree. 2. Adaptivecreationofgeometry: Ineachquadtreecell,ahyper-pointisestimated byminimizinga2.5Dquadraticerrorfunction(2.5DQEF).Geometrysimplifica- tion is achieved in an adaptive manner by collapsing subtrees and adding QEFs associatedwithleafcells(Figure5.3(c)). 3. Polygon generation: A watertight mesh model is created by connecting hyper- points with surface polygons (turquoise triangles) and boundary polygons (pur- ple triangles), which form building roofs and vertical walls, respectively (Fig- ure5.3(d)). 61 4. Principaldirectionsnapping: Theroofboundariesarerefinedtofollowtheprin- cipaldirectionsdefinedinSection3.4.3andSection4.3.4(Figure5.3(e)). 5.3 2.5DScanConversion Thefirststepofthe2.5Ddualcontouringalgorithmconvertstheinputpointcloudintoa volumetric form, by sampling Hermite data (a point-normal pair) over a 2D supporting grid. Withelementsbeingconsideredastheirinfiniteextensionsalongtheverticaldirec- tion,this2Dgridhasa3Dvolumetricconnotation. E.g.,agridcellrepresentsaninfinite threedimensionalvolume,whileagridpointcorrespondstoaverticallinecontainingit. 5.3.1 SurfaceHermiteData Given a 2.5D point cloud as input, the scan conversion algorithm first segments it into multiple roof layers using a local distance-based region growing algorithm 1 , as shown in Figure 5.4(a). Ideally, each vertical line passing through a grid point intersects with one and only one roof layer. The intersection point is taken as a surface Hermite data sample,andestimatedbyaveragingtheheightsandnormalsofitsk-nearestinputpoints within the same roof layer, illustrated as points marked with blue or purple outlines (takingk = 4)inFigure5.4(a). The only difficulty in this process is to robustly detect the right roof layer crossing the vertical line. Intuitively, a roof layer L covers a grid point g iff each of g’s four neighboring cells contains at least one input point p belonging to L or a higher cluster L ′ . For example, in Figure 5.4(a), point A is covered by no roof layers, and thus is assigned as ground; point B is only covered by and assigned to the dark-grey layer. 1 Therooflayersarealwayssegmentedinalocalarea,asglobalsegmentationmayeraselocalfeatures such as those shown in Figure 5.9(c). Specifically, the segmentation for grid point g is applied to all the inputpointsing’sfourneighboringcells. 62 Points of a high roof layer Points of a low roof layer Surface samples assigned to the high roof layer Surface samples assigned to the low roof layer Surface samples assigned as ground A B C Boundary Hermite data samples E F G H Support vectors (a) (b) (c) Figure 5.4: Generating (a) surface Hermite data samples on grid points: the sample is assigned to the highest roof layer which covers the grid point; (b,c) boundary Hermite data samples on grid edges: the maximum margin line (thin black lines) divides the lowersurfaceHermitedatasamplefromthehigherrooflayer. Note point C is covered by both the dark-grey layer and the light-grey layer. In this case, the scan conversion algorithm takes the highest roof layer covering point C, i.e., thelight-greylayer. 5.3.2 BoundaryHermiteData While surface Hermite data captures the surface geometry of building roofs, the shapes ofroofboundariesarerepresentedbythe boundary Hermite data. ConsideringagridedgeeconnectingtwogridpointswithsurfaceHermitedatasam- ples{s 0 ,s 1 } on different roof layerss 0 ∈ L 0 ,s 1 ∈ L 1 , 2 the vertical wall connectingL 0 andL 1 shouldsplittheirprojectionimagesonthex-yplane. Inspiredbythe2DSupport VectorMachinealgorithm[2],amaximum-marginlinel isdetectedwhichseparatesL 0 2 To avoid ambiguity, roof layers are determined again by a local segmentation over {s 0 ,s 1 }∪ P, whereP istheinputpointsetwithine’stwoadjacentcells. 63 and L 1 on the x-y plane. The boundary sample is estimated by intersecting line l and edgee. In practice, with the existence of residual sensor noise, the projections of different roof layers may overlap on the x-y plane. Since aerial LiDAR data is collected from a top view, more saliency is given to the higher roof layer L 1 (assuming height(L 0 ) < height(L 1 )). Thus the scan conversion algorithm takes the maximum-margin line l whichseparates{s 0 }andL 1 whilemaximizingdistance(s 0 ,l),shownasthethinblack lines in Figure 5.4(b,c). Empirically, this method is more robust than other methods includingthatusingamaximum-soft-marginlinedividingL 0 andL 1 . 5.4 AdaptiveCreationofGeometry Given a quadtree cellc (not necessarily being a finest-level leaf cell), the set of surface Hermite data samples on the grid points in c is denoted as S, and the set of boundary Hermite data samples on atomic grid edges in c is denoted as B. The roof layers in c are then determined by segmenting S into k clusters S = S 1 ∪···∪S k . Intuitively, if anatomicgridedgeinchasnoboundarysampleattached,itconnectstwosurfacesam- ples of the same roof layer. Thus, an agglomerative clustering algorithm is adopted to repeatedlycombinesurfacesamplesetsconnectedbyedgeswithoutboundarysamples. Now the task is to generate k vertices for the k roof layers, denoted as a hyper- point χ = {x 1 ,...,x k }. To maintain the consistency of roof layer boundaries, these k vertices are required to have the same projection on the x-y plane, i.e., they should have the same x-y coordinates, but different z values. Thus χ can be expressed as a k + 2 dimensional vector χ = (x,y,z 1 ,...,z k ). Let x 0 = (x,y,0) for convenience in followingdiscussions. 64 5.4.1 2.5DQuadraticErrorFunction The hyper-point χ is optimized by minimizing a 2.5D quadratic error function (2.5D QEF)definedasthelinearcombinationof2Dboundaryquadraticerrorsand3Dsurface quadraticerrors: E(χ) = ∑ (p;n)∈B (ωn·(x 0 −p)) 2 + ∑ i=1;:::;k ∑ (p;n)∈S i (n·(x i −p)) 2 (5.2) where ω is a user-given weight balancing between boundary samples and surface sam- ples. Empirically,aweightbetween1∼ 4satisfiesmostoftheexperiments. Due to the horizontality of boundary sample normals, the third coordinates of p and x 0 do not affect the 2D error term. However, I choose to write all these variables uniformlyin3D,inordertoexpresstheenergyfunctioninamatrixproductform: E(χ) = (Aχ−b) T (Aχ−b) (5.3) where A is a matrix whose rows come from normals in B,S 1 ,...,S k , with those in B multipliedbyω. Thex-yvaluesofeachnormalareplacedinthefirsttwocolumns,while thezvaluesofnormalsinS i areplacedinthe(i+2)-thcolumn. Theremainingentries inAarepaddedwithzeros. bisavectorcomposedofcorrespondinginnerproductsn·p withthefirst|B|entriesmultipliedbyω. 65 ThenumericalstabilityduringQEFoptimizationcanbeimprovedbyemployingthe QRdecompositionasproposedin[31]. I.e., (A b) =Q           ^ A ^ b 0 r 0 0 ... ...           (5.4) whereQisanorthogonalmatrixandEquation5.3canberewrittenas: E(χ) = (Aχ−b) T QQ T (Aχ−b) = ( ^ Aχ− ^ b) T ( ^ Aχ− ^ b)+r 2 . (5.5) Thus,E(χ)isminimizedbysolving ^ Aχ− ^ b = 0. Tohandlethepossiblesingularity of ^ A, an SVD decomposition is applied, following solutions in previous methods [31, 37]: ^ A =UV T . (5.6) Smallsingularvaluesinwithamagnitudeoflessthan0.1istruncated,andthepseudo- inverse + isadoptedtocomputethehyper-pointχas: χ = χ+V + U T ( ^ b− ^ A χ) (5.7) where χisaguessedsolutionwhosefirsttwocoordinatescomefromthecentroidofB, andthe(i+2)-thcoordinateisthemeanheightofsamplesinS i . IfB isempty,thefirst twocoordinatesequaltothoseofthecentroidofS. 66 5.4.2 QuadtreeSimplificationwithQEFs Taking a quadtree with QEF matrices pre-computed for all the finest-level cells, the geometry is simplified by collapsing leaf cells into parent cells and combining QEFs in a bottom-up manner. A user-given tolerance δ controls the simplification level by denyingsub-treecollapsewhentheresidualisgreaterthanδ. Combining four regular 3D QEFs can be simply achieved by merging the rows of their upper triangular matrices to form a 16× 4 matrix [31]. 2.5D QEF matrices are combined in a similar way, yet with the consideration of association between matrix columnsandrooflayers: asrooflayersinleafcellsmergeintoonerooflayerinthepar- ent cell, corresponding matrix columns are placed in the same column of the combined matrix. Specifically, the roof layer is re-segmented in the parent cell before merging matrices. Assuming thei-th rooflayer inaleaf cellbelongsto thej-throoflayer inthe parentcell,the(i+2)-thcolumnoftheleafcellmatrixisputintothe(j+2)-thcolumn of the combined matrix. 0-columns are used to pad the leaf cell matrices where no roof layersbelongtocertainrooflayersintheparentcell. Once again, the merged matrix is brought to the upper triangular form via a QR decomposition. Due to the orthogonality of involved transformation matrices, it repre- sentsthe2.5DQEFintheparentcell. Figure 5.5 shows a building model without and with geometry simplification. The triangle number decreases from 820 to 252 with an insignificant increase of the fitting error. 67 Sum of squared distances = 33.53 Sum of squared distances = 53.22 Figure5.5: 2.5Ddualcontouringwithout(left)andwith(right)adaptivegeometrysim- plification A A A A Surface polygons Boundary polygons Relevant grid points Relevant grid edges (a) (b) (c) (d) Figure 5.6: (a,b) Creating surface polygons (colored hollow polygons) and boundary polygons (colored semitransparent polygons) around hyper-pointA. Viewing from top, (c) surface polygons are generated at grid points, while (d) boundary polygons are pro- ducedforgridedgeswhichexhibitarooflayergap. 5.5 PolygonGeneration Giventhesimplifiedquadtreewithhyper-pointsestimatedineachleafcell,thenexttask istocreatepolygonsconnectingthesehyper-pointsintoamesh. Inparticular,twokinds ofpolygonsaregeneratedtosatisfythe2.5Dcharacteristic. 1. Surface polygons: At each grid point p, a surface polygon is created by con- necting vertices in the hyper-points on the same roof layer asp in its neighboring cells. 68 2. Boundary polygons: At each minimal quadtree edge e, a boundary polygon is generatedconnectingtwohyper-pointsegmentsintheadjacentcells. Figure 5.6 shows an example of polygon generation around a hyper-point A. The surfacepolygonsandboundarypolygonsarehighlightedwithcoloredoutlinesandcol- ored semitransparent polygons respectively. To avoid cracks generated within a hyper- point,aboundarypolygonsequentiallypassesthroughtheverticesinahyper-pointseg- ment in height ascending or descending order. E.g., the dark-blue boundary polygon in Figure 5.6 goes through all the three vertices in hyper-point A, from the top vertex to thebottomvertex. Thispolygongenerationmethodisguaranteedtoproducecrack-freemodels,which can be derived from the fact that except for the border edges created around the entire grid, the other mesh edges are contained by an even number of polygons. Proof is straightforward: a non-vertical mesh edge is either contained by two surface polygons, or by one surface polygon and one boundary polygon. As for the vertical mesh edges created within a hyper-point, considering all the boundary polygons around this hyper- point (e.g., the colored semitransparent polygons shown in Figure 5.6(a,b)): they go up and down though this hyper-point and finally return to the start vertex, forming up a closededgeloop. Thus,eachverticalmeshedgeinthishyper-pointappearseventimes. 5.5.1 SharpFeaturePreservingTriangulation By minimizing QEFs, 2.5D dual contouring has the ability to produce vertices lying on sharpfeatures,whichareacommonpatterninbuildingroofs. However,apoortriangu- lationofsurfacepolygonscanspoilthisadvantage,asshowninFigure5.7left. Tosolve thisproblem,anefficientsharpfeaturedetectionalgorithmisproposedtopreservethese featuresoncedetected. 69 0.0 1.0 Sum of distance = 635.10 m 2 2 Sum of distance = 36.69 m 2 2 Figure 5.7: Triangulation without (left) and with (right) the sharp feature preserving algorithm. Thecolorsofinputpointsrepresentthesquareddistancesfromthemesh. In a grid cell c containing only one roof layer, covariance analysis is applied over thenormalsofallsurfacesamples, i.e.,togettheeigenvaluesofmatrix: C = 1 N ∑ i n i ·n T i . (5.8) SincechasnoboundaryHermitedatasamples,Equation5.4and5.6canbeusedto simplifythismatrixas: C = 1 N A T A = 1 N ^ A T ^ A = 1 N V T V T . (5.9) Thus, the diagonal of matrix 1 N T gives the eigenvalues of C, while the columns of V are corresponding eigenvectors. As Pauly [47] suggests, the smallest eigenvalue λ 0 and the middle eigenvalue λ 1 estimate the minimal and maximal curvatures, as the corresponding eigenvectors v 0 ,v 1 point to the curvature directions. Therefore, ridges and valleys are detected by finding vertices with smallλ 0 and fairly largeλ 1 . v 0 is then adopted as the feature direction. Since the involved matrices have all been computed in previoussteps,theadditionaloverheadofthisalgorithmistrivial. 70 Figure 5.8: Comparison between topology-unsafe simplification (left) and topology- safe simplification (right). Undesired features can be created by merging leaf cells in a topology-unsafemanner. Specifically, for each diagonal e of a surface quad, the triangulation algorithm cal- culates: ∑ p∈eand 0 (p)< λ 1 (p)·|v 0 (p)·e| (5.10) andchoosesthediagonale ∗ whichmaximizesthisvaluetosplitthequadintotwotrian- gles. Hereτ isausergiventhreshold. Empirically,τ = 0.01. 5.6 Topology-SafeSimplification So far the quadtree simplification is completely built on QEFs, and the topology of output models may change during this process. Undesired features can be generated as showninFigure5.8left. Tosolvethisproblem,anadditionaltopologytestisintroduced rightbeforesub-treecollapsehappens. Itrejectscollapseifthereisadangeroftopology change. Regarding multiple roof layers as multiple materials, topology test algorithm in [31] is extended with an additional test (step 3) which prevents different roof layers in one leaf cell (top-left cell in Figure 5.9(a)) from merging into a same roof layer in the coarse cell (Figure 5.9(b)). This situation may cause removal of small vertical wall features(e.g.,Figure5.9(c)). 71 Surface samples with different roof layer assignments Surface polygons Boundary polygons (a) (b) (c) Boundary samples exhibiting roof layer gaps Figure 5.9: An unsafe simplification case denied by the topology safety test step 3. Since the center grid point has different roof layer assignments in these leaf cells, two different layers in the top-left leaf cell (a) belong to the same roof layer in the coarse cell(b). Unsafemergingmayerasewallfeaturessuchastheoneshownin(c). 1. Testwhethereachleafcellcreatesamanifold;ifnot,stop. 2. Testwhetherthecoarsecellcreatesamanifold;ifnot,stop. 3. Test whether any two roof layers in a same leaf cell belong to two different roof layersinthecoarsecell;ifnot,stop. 4. Testwhetherthetopologyofthedualcontourispreservedusingfollowingcriteria; ifnot,stop;otherwise,collapse. (a) Test whether the roof layer on the middle point of each coarse edge agrees withtherooflayeronatleastoneofthetwoedgeendpoints. (b) Testwhethertherooflayeronthemiddlepointofthecoarsecellagreeswith therooflayeronatleastoneofthefourcellcorners. 72 Principal directions Figure 5.10: Roof layer boundaries (thick colored lines) are regularized using principal directionsnappingalgorithm. 5.7 PrincipalDirectionSnapping The 2.5D dual contouring algorithm is completely data-driven, i.e., no pre-assumptions about the roof shapes have been made. Thus this algorithm can handle complex roofs in a robust manner. On the other hand, in some cases, prior knowledge of the urban areaisgivenanditisadesiretohavebuildingmodelsconcurringwithsuchknowledge. This section shows a post-processing refinement to the modeling results using the prior knowledgeofprincipaldirectionsasdetectedinSection3.4.3andSection4.3.4. The idea is straightforward: once the boundaries of individual roof layers are extracted, they can be snapped to the principal directions as much as possible without exceeding a small error tolerance. In order to maintain the consistency between bound- aries of different layers, the boundaries are handled one by one in height-descending order. I.e., when a roof layer boundary has been processed, the x-y coordinates of the touched hyper-points are fixed, which are then considered as constraints during the subsequent processing of lower roof layers. Figure 5.10 shows clean and simple roof boundariesgeneratedbytheprincipaldirectionrefinement. 73 5.8 ExperimentalResultsof2.5DDualContouring Figure5.11showsanurbanareaofLosAngelesreconstructedfrom26MLiDARpoints with7samples/m 2 resolution. ThestreamingurbanmodelingsystempresentedinChap- ter 4 is first adopted to remove irrelevant parts such as noises, trees, vehicles and even ground. The 2.5D dual contouring algorithm is then tested on point clouds of individ- ual buildings. This algorithm successfully creates 1,879 building models consisting of 857Ktriangleswithin 6 minutes on aconsumer-levellaptop(IntelCore21.8GHzCPU with 2GB memory). 2.5D building models with complex roofs are robustly generated intheentirearea. To further demonstrate the ability of handling various kinds of building models, the 2.5D dual contouring method is tested on a set of buildings from the city of Atlanta, as illustrated in Figure 5.1. Figure 5.12 shows a comparison between this method and the state-of-the-arts. In particular, I compare the average squared distance from input pointsetstothegeneratedmodels,andtheratioofpointswithsquareddistancesgreater than1m 2 InFigure5.12,pointcolorsdenotethesquareddistances,andthecoloredbars show the percentage of points at different squared distance levels. As the quantitative resultsinTable5.1illustrate,the2.5Ddualcontouringmethod(firstcolumn)isthemost accurate algorithm to produce 2.5D models. Plane-based approaches such as the one in Section 3.4.3 (second column) are unable to handle non-flat roofs (a,d) and small roof features (b,e). Cracks often exist when fitting is unsuccessful (c,d). A general mesh simplification over the DEM (third column) is competitive in the sense of fitting qual- ity. However, it cannot produce 2.5D models composed of roofs and vertical walls. In addition,thefittingqualityonroofboundariesisunsatisfactory(f,g,h). Thelastcolumn demonstrates point clouds aligning with manually created models. Designed without knowledge from real-world data, they often lack of accuracy even after registration to theinputpoints. 74 Figure5.11: Buildingreconstructionfora2KM-by-2.5KMurbanareaofLosAngeles Figure 5.13 finally demonstrates the influence of grid configuration. As an adaptive approach, 2.5D dual contouring is insensitive to the grid size (top row). In addition, it has the ability to place vertices at optimal positions, thus grid orientation affects the resultsinsignificantly(bottomrow). 75 0.0 Squared Distance 1.0 a b c f d g e h f d g e,h Figure 5.12: Building models created using different approaches (from left to right): 2.5Ddualcontouring,plane-basedmethodproposedinSection3.4.3,generalmeshsim- plificationoverarasterizedDEM,andmanualcreation. Pointcolorsdenotethesquared distances between points and generated models. Color bars under the models show the ratioofpointsatdifferentsquareddistancelevel. l = 0.7 m l = 1.0 m l = 1.4 m Tri. # = 835 Err avg = 0.009 Tri. # = 688 Err avg = 0.010 Tri. # = 584 Err avg = 0.016 Tri. # = 688 Err avg = 0.010 Tri. # = 878 Err avg = 0.016 Tri. # = 937 Err avg = 0.024 θ = 0 o θ = 30 o θ = 45 o Default grid configuration: Grid size: l = 1.0 m Grid orientation: θ = 0 o Tri. # = 688 Err avg = 0.010 Figure 5.13: Models of similar quality are generated with the same point cloud embed- dedintogridsofdifferentsizesordifferentorientations. 76 ModelsinFigure5.12 2.5Ddual contouring Plane- based method DEMsim- plification Manual creation [21] Firstrow (4679points) Trianglenumber 214 76 198 78 Avg. distance 2 0.016 0.599 0.061 0.058 Outlierratio 0.06% 12.37% 0.53% 0.83% Secondrow (684907points) Trianglenumber 8009 6262 8000 1227 Avg. distance 2 0.037 0.465 0.035 7.780 Outlierratio 0.44% 7.93% 0.87% 70.38% Thirdrow (198551points) Trianglenumber 12688 1619 12999 1558 Avg. distance 2 0.203 1.610 0.264 16.220 Outlierratio 2.03% 21.15% 3.08% 68.28% Table 5.1: Quantitative evaluation of four modeling approaches over models shown in Figure5.12 77 Chapter6 2.5DBuildingTopology This chapter studies 2.5D building topology, and extends 2.5D dual contouring into a 2.5Dbuildingmodelingmethodwithtopologycontrol. Themajorcontributionofthischapterisbasedontheobservationthathumanvision tend to be more sensitive to building topology rather than building geometry. Intu- itively, building topology determines the existence of structural pieces and the connec- tions between them; while building geometry describes where these structural pieces appear in the three dimensional space. Humans tend to be more aware of changes in topology even if the related structural pieces are small. For example, Figure 6.1(c,d) demonstrate two building models created targeting to achieve more precise geometry and more precise topology respectively. Although the left model fits the input point cloud better under typical geometrical error measurements (e.g., average quadratic dis- tance),itisvisuallylessconvincingthantherightonebecausearoofpiece(thechimney) ismissing. In 2.5D dual contouring, the topology issue is alleviated by introducing a topology testpresentedinSection5.6. Theadaptivesimplificationprocessfirstcollapsesquadtree cells and optimizes an anchor point completely based on geometric errors without con- sidering building topology; then rewinds the collapse operation if the topology test revealsapossibletopologychange. Thisstrategyperformswellunderstronggeometric control (i.e., with a small geometry error tolerance). However, in cases where simpler models are desired thus looser geometric control is given, the number of topology test failures increases rapidly and they become the dominant factor in preventing quadtree 78 (a) Input LiDAR (c) Geometrically more precise (d) Topologically more precise (b) Unsimplified model Figure6.1: Buildingmodelsreconstructedtargetingtoobtain(c)moreprecisegeometry and(d)moreprecisetopologyrespectively. Comparedwith(a)theinputLiDARand(b) unsimplifiedbuildingmodel,themissingofthechimneymakestheformeronevisually lessconvincingthanthelatterone. collapse. Figure6.2(a)showssuchanexampleinwhichtopologytestfrequentlydetects possiblerooflayercracksanddeniesthecellcollapse;therefore,numerousinsignificant triangles are produced along the thin long roof features as shown in the closeup. The deepreasonbehindthisproblemisthattheoptimizationprocessiscompletelyunaware of building topology. It produces exact one hyper-point 1 per quadtree cell without dis- crimination. Hence,themostcomplicatedtopologicalstructurethatcanexistinonecell isaconjunctionhyper-pointwithstar-shapedroofboundaries,asshowninFigure6.2(a) bottom right. The adaptive simplification becomes problematic in producing building 1 A hyper-point is defined as a series of 3D points having the same x-y coordinates but different z values,seeChapter5. 79 (a) 2.5D dual contouring (b) 2.5D contouring with topology control Denied Accepted Invoke rewind Figure 6.2: Comparison between (a) 2.5D dual contouring and (b) 2.5D building mod- eling with topology control. While the uniqueness of hyper-point in one cell prevents a flexiblesimplificationindualcontouring,thenewmethoddetectsandcontrolsbuilding topologybeyondtherigidquadtreestructure. structureswithtopologythatismorecomplexthanaconjunctionhyper-point. Collapse rewindisinvokedfrequently. This chapter proposes an extension to the 2.5D dual contouring method to enable building topology control. The key idea is to maintain multiple hyper-points in one quadtree cell. Therefore complicated in-cell building topology is allowed. With this extension,theadaptivemodelcreationprocedurebecomeslessrestrictive,andthusgen- erates simpler building models in a flexible manner, e.g., Figure 6.2(b). In particular, without changing building topology, the modeling method can produce building mod- els with triangles as few as manually created models or primitive-based models; while it still provides a similar geometric optimization scheme as the data-driven modeling approaches. 80 Section 6.1 reviews topology control approaches in volumetric modeling. Sec- tion 6.2 formally defines 2.5D building topology which is the basis of this chapter. Section 6.3 then extends 2.5D dual contouring into a building modeling method with topologycontrol. ExperimentalresultsareshowninSection6.4. 6.1 ReviewofTopologyControlinVolumetricModeling In classic 2D and 3D volumetric modeling methods, the topology issue is first noticed byJuetal.[31]. Theyproposeatopologytestmechanismtorejectsimplificationopera- tionsyieldingpossibletopologychanges. Thismechanismislaterextendedtothe2.5D dualcontouringmethodasdiscussedinSection5.6. The drawback of creating one vertex per octree/quadtree cell is noticed by researchers and different approaches have been proposed to solve this problem, e.g., [51, 59, 66]. Generally, they all allow one grid cell to have more than one vertices, in ordertotrackcontourcomponentswhicharetopologicallymorecomplicatedthanadisk (or in 2D, a line segment) in each cell. Although these methods share some similarities with the approach presented in this chapter, two key differences are noted: first, this research aims at 2.5D building modeling involving hyper-points that cannot be handled in classic 2D or 3D manner; second, various building topology features are defined and processedwhicharemorecomplicatedthandisk-likefeaturesinclassic2Dor3Dspace. 6.2 2.5DBuildingTopology Considering a 2.5D dual contouring process without adaptive simplification: taking aerial LiDAR data as input, the contouring method builds up a uniform grid with Her- mite data attached; creates one hyper-point (a series of points that are consistent on the x-y plane) in each cell by optimizing a quadratic error function; generates surface 81 Different roof features Markers representing different vertical wall features Typical point features Figure6.3: Topologicalfeaturesinanunsimplified2.5Dbuildingmodel polygons and vertical boundary polygons; and produces an unsimplified 2.5D building modelsuchastheoneshowninFigure6.3left. 6.2.1 TopologicalFeatureDefinitions The first observation about 2.5D building topology is that a typical building structural piece(e.g.,achimney)isusuallycomposedofoneuniqueroofpatchanditssurrounding walls. Hence, Definition 6.1 A roof feature R is a connected component composed of non-vertical surfacepolygons. Rooffeaturesarethekeyindetermining2.5Dbuildingstructures. Figure6.3utilizes green mesh pieces rendered with different color intensities to represent multiple roof features. Inparticular,neighboringroofpatchesexhibitaheightgapalongtheircommon boundary, which is sealed up by vertical boundary polygons (grey vertical polygons). ThesepolygonsformwallfeaturesthataremarkedinFigure6.3bycurveswithvarious colors. 82 0 0 0 1 0 0 1 1 0 0 0 0 A B Folding points R1 = R2 W b1 b2 F F Figure 6.4: 2.5D building models may contain point features involving only one wall feature (left). These points are produced around grid edges (e.g., AB) which detect inconsistentrooflayerassignmentsintwoadjacentcells(right). Definition 6.2 Given two roof features R 1 and R 2 , the corresponding wall feature W is defined as the connected component composed of vertical boundary polygons adjacent to R 1 and R 2 simultaneously. W intersects with R 1 and R 2 via non-identical roofboundarypolylines. Here a slight modification is made to 2.5D dual contouring, that the triangulation is disabled. Instead, surface polygons are treated as quads connecting all four vertices aroundagridcorner,whileboundarypolygonsareproducedwithtwonon-verticaledges linking a pair of neighboring hyper-points. Therefore, a vertical boundary polygon is “adjacent” to a roof feature as long as they share a non-vertical edge. By tracking the consecutive non-vertical edges, one roof boundary polyline is created regarding each adjacent wall-roof feature pair (W,R i ),i = 1,2, denoted as b i = W ∩ R i ,i = 1,2. According to definition 6.2, W is a valid feature whenb 1 andb 2 are not identical, even ifR 1 andR 2 refertothesameroofpatch. Definition6.3 GivenawallfeatureW sharingtwoconsecutiveroofboundarypolylines b 1 andb 2 withR 1 andR 2 respectively, point features are defined as hyper-points which containb 1 andb 2 ’sendpointsifthereisany. 83 Typically, a hyper-point shared by two neighboring wall features is a point feature. They are rendered in Figure 6.3 as golden balls connected by vertical lines. Note that definition6.3supportsthesecommonhyper-pointssharedbyneighboringwallfeatures; but is not limited to them. In a special case shown in Figure 6.4, the corner points of grid edge AB have the same roof layer assignment in one adjacent cell (right cell) but different assignments in another (left cell). Thus, a vertical boundary polygon is produced to reflect this significant topology feature. Specifically,R 1 andR 2 denote the same roof patch, which is adjacent withW alongb 1 andb 2 representing the upper roof boundary polyline and the bottom roof boundary polyline respectively. Sinceb 1 andb 2 are non-identical, W is a valid wall feature. Point feature F acts as a “folding point” whichconnectsb 1 andb 2 togetherandfoldsuptheboundaryofW. 6.2.2 ConnectionsbetweenTopologicalFeatures By projecting the 2.5D building model onto the x-y plane, the building topology can be viewed with a 2D cell complex representation 2 . In particular, roof features, wall features, and point features are projected onto the 2D x-y plane as 2-cells (regions), 1- cells (polylines), and 0-cells (points) respectively. High dimensional cells are always boundedbyasetoflowdimensionalcells. Given a projection operator P(·) which projects a set of 2.5D objects onto the x-y planeandaboundaryextractionoperator∂(·),theconnectionsbetweenrooffeatureset R,wallfeaturesetW,andpointfeaturesetP arerevealedwithfollowingequations: P(∂R)⊆P(W), foranyR∈R, (6.1) ∂P(W)⊆P(P), foranyW ∈W. (6.2) 2 Cellcomplexesarethebasicconceptsinalgebraictopology,see[23]fordetaileddiscussion. 84 R0 W0 P(∂R0) = P({W0}) ∂P(W0) = Ø R0 R1 W0 W1 P(∂R0) = P({W0, W1}) P(∂R1) = P({W1}) ∂P(W0) = Ø ∂P(W1) = Ø R0 R1 W0 W1 W2 P0 P1 P(∂R0) = P({W0, W1}) P(∂R1) = P({W1, W2}) ∂P(W0) = P({P0, P1}) ∂P(W1) = P({P0, P1}) ∂P(W2) = P({P0, P1}) R0 P0 P1 W0 W1 P(∂R0) = P({W0, W1}) ∂P(W0) = Ø ∂P(W1) = P({P0, P1}) R0 R1 R2 W0 W1 W2 W5 W4 P(∂R0) = P({W0, W1, W2, W3, W5}) P(∂R1) = P({W5}) P(∂R2) = P({W2, W4}) P0 W3 P1 P2 P3 ∂P(W0) = P({P1, P3}) ∂P(W1) = P({P1, P2}) ∂P(W2) = P({P2, P3}) ∂P(W3) = P({P0, P1}) ∂P(W4) = P({P2, P3}) ∂P(W5) = Ø (a) (b) (c) (d) (e) Figure6.5: 2.5Dbuildingtopologyisrepresentedbytopologicalfeatures,andtheasso- ciations between them in form of Equation (6.1) and (6.2). Examples include typical building structures such as (a) individual building blocks, (b) blocks with top attach- ments, (c) blocks with side attachments, (d) stair-shaped structures, and (e) combina- tionsofthesepatterns. These equations can be straightforwardly derived from the definitions of roof, wall, and point features. On the other hand, once topological features and their associations expressed in form of Equation (6.1) and (6.2) are fixed, the 2.5D building topology is determinedaccordingly. Figure 6.5 demonstrates typical building structures including standing-alone build- ingblocks,verticallyattachedblocks,horizontallyattachedblocks,stairshapes,andthe combinationsofthesepatterns. Nevertheless,the2.5Dbuildingtopologyrepresentation describestheminadeterministicanddifferentiablemanner. 85 In addition, the difference between 2.5D topology representation and classic 2D topology representation is noted as follows. The latter one can be achieved by project- ing all building elements onto the x-y plane at first, and treating different roof layers as multiple region materials. This representation, however, is problematic in handling wall features connecting the gap within one roof layer (e.g., W 1 in Figure 6.5(d)). It eliminates such wall features together with the folding point (e.g., P 0 ). In contrast, the topologyrepresentationpresentedinthissectionfaithfullypreservesallsignificant2.5D topologyfeatureswhicharethebasisofthetopologycontrolmethod 3 . 6.3 ContouringwithTopologyControl So far the 2.5D topological features are formally defined and the associations between them are well explored. These mechanisms can be naturally expanded from a uniform grid toa quadtree. Thus, the quadtree-based simplificationof 2.5D dual contouring can beextendedwithatopologycontrolmethodtomaintainthe2.5Dbuildingtopology. 6.3.1 Hyper-PointsClustering The core of the simplification algorithm is to optimize the geometry of hyper-points based on a 2.5D quadratic error function (see Section 5.4). These hyper-points can be categorizedbythenumberoftheirlayers. 1. 1-layer points: A hyper-point containing one vertex is optimized targeting the disk-likegeometryinagridcell. Inmostcases,itisconnectedtoverticescreated in its neighboring cells only by surface polygons. Such a hyper-point is an inner 3 Thisproblemcanalsoberegardedastheresultofchangingtheorderinapplyingboundaryextraction operator and projection operator. As 2D topology representation projects all elements onto the x-y plane at first, it attempts to replaceP(∂R) in Equation (6.1) with ∂P(R). This attempt is problematic because ∂P(R)̸≡P(∂R)asP(·)canabsorbwallfeaturessuchasW 1 inFigure6.5(d). 86 vertex of a roof feature. Thus it can be safely merged into a neighboring vertex without changing the 2.5D building topology. The accepting vertex can be either another 1-layer point or one vertex in a hyper-point with more than one layers. Theonlyexceptionisa1-layerfoldingpoint,whichisconnectedtotwodifferent layersofaneighboringhyper-pointbyaboundarypolygon,asshowninFigure6.4 andFigure6.6left. Inthiscasethe1-layerpointisapointfeature,thusshouldnot bemergedintootherpoints. 2. 2-layer points: A hyper-point containing two layers is a typical roof boundary anchor point, which is optimized with surface geometry and boundary geometry simultaneously. Typically, it is an inner element of a wall feature, thus can be merged into another hyper-point that is connected to it by a vertical boundary polygon. The accepting hyper-point can be either another typical 2-layer point or a multi-layer point. Similar to 1-layer points, folding points can exist within 2-layer points. Figure 6.6 middle showssuch an example where the 2-layer point isapointfeature. Itcannotbemergedintypicalmanner. 3. Multi-layerpoints: A hyper-point with more than two layers can be regarded as a conjunction point of more than two regions, if the problem is viewed on the 2D projectionplanewhereroofpatchesaretreatedasregionswithdifferentmaterials. Ifindthatanymulti-layerpointisapointfeature. Proofisstraightforwardasthey cannot be the inner elements of wall features; thus always stand at boundaries of wall features’ 2D projections. Therefore, multi-layer points can only accept 1- layer points and 2-layer points joining it, but cannot be merged with other point features. 87 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 1 0 1 0 0 0 0 0 0 0 1 1 0 1 0 2 1 2 1 0 0 2 0 2 0 0 2 2 2 2 2 Different roof patches Wall polygons Folding points Figure 6.6: Folding points can be part of a 1-layer hyper-point, a 2-layer hyper-point, or a multi-layer hyper-point (from left to right). They are detected and marked as point featuresbeforeadaptivesimplificationstarts. Since folding points can exist in 1-layer points, 2-layer points, and multi-layer points,asdemonstratedinFigure6.6,theyaredetectedandmarkedaspointfeaturesdur- ing pre-processing. The detection algorithm is implemented by uncovering grid edges withinconsistentrooflayerassignmentsinadjacentcells. Apart from folding points, typical 1-layer points, typical 2-layer points, and multi- layer points are denoted as p 1 , p 2 , and p m respectively. A set of hyper-point clustering operationsareintroducedtomergecomponentsconnectedbysurfacepolygons,denoted as 1→1 S , 1→2 S ,and 1→m S ,withfollowingfunctions: 1→1 S : {p 1 1 ,p 1 2 ,...,p 1 n } ⇒ p 1∗ , (6.3) 1→2 S : {p 1 1 ,p 1 2 ,...,p 1 n },p 2 ⇒p 2∗ , (6.4) 1→m S : {p 1 1 ,p 1 2 ,...,p 1 n },p m ⇒p m∗ . (6.5) Each operation merges a connected component within a roof feature, and produces one hyper-point (e.g., p 1∗ , p 2∗ , or p m∗ ) per cluster. In particular, the geometric coordinates ofoutputhyper-pointsareobtainedbyoptimizinga2.5DQEFmatrixwhichisthecom- binationofQEFmatricesfrominputhyper-points. Thethirdcolumnfromp 1 i ’smatrices 88 are placed with corresponding matrix column from p 2 or p m . Details of QEF matrices combination can refer to Section 5.4.2. Similarly, the clustering operations for compo- nentsthatareconnectedbyverticalboundarypolygonsare: 2→2 B : {p 2 1 ,p 2 2 ,...,p 2 n } ⇒ p 2∗ , (6.6) 2→m B : {p 2 1 ,p 2 2 ,...,p 2 n },p m ⇒p m∗ . (6.7) Theseoperationsmergeconnectedcomponentswithinawallfeature. 6.3.2 HandlingDegenerateCases Although these clustering operations aim at simplifying continuous roof and boundary features,theyriskinmakingtopologicalfeaturesdegenerate. Forexample,withintense geometry simplification, the building structure in Figure 6.5(a) may degenerate into a singleverticallinewithitstopvertexcollapsedfromR 0 . To address this problem, degenerate tests are adopted after each clustering opera- tion. Given f(C) and e(C) as the number of faces and edges in a cell complex C, the followingrequirementsmustbesatisfied: f(R)≥ 1, foranyR∈R, (6.8) e(P(W))≥ 1, foranyW ∈W. (6.9) Initially, these two criteria are fulfilled due to Definition 6.1 and 6.2. In adaptive sim- plificationphase,foreachpossiblehyper-pointclusteringoperation,adegeneratetestis applied to all the modified roof features and wall features, to check if these two criteria arestillsatisfied. Ifanyofthemisviolated,theclusteringoperationisrewound. 89 Inpractice,Ifoundthisdegeneratetestinefficientbecauseitinvolvespolygonrecre- ation and topological feature detection after every clustering operation. Thus an equiv- alent criterion is proposed to replace the former two, which is much easier to be imple- mented: Degenerate test A hyper-point clustering operation passes the degenerate test if the followingcriterionstaystrue: e(∂P(R))≥ 3, foranyR∈R. (6.10) The equivalence between this degenerate test and the one based on Equation (6.8) and(6.9)isprovedintheAppendixofthischapter. Since the degenerate test is irrelevant to typical 1-layer points, clustering opera- tions based on surface components are always allowed (i.e., 1→1 S , 1→2 S , 1→m S ). For clusteringamong2-layerpointsandmulti-layerpoints,thesimplificationalgorithmpre- computes boundaries of roof feature projections, i.e., ∂P(R), and keeps tracking of the edge numbers e(∂P(R)) during the complete simplification process. Boundary com- ponent clustering operations (i.e., 2→2 B and 2→2 B ) decrease the corresponding edge numbers. Onceaboundaryedgenumberislessthan3,thelatestoperationisrewound. 6.3.3 AdaptiveContouring 2.5D dualcontouring is extended by allowingmultiplehyper-pointsin eachgridcellof the quadtreeQ. In addition to the adaptive structure of quadtree, a hyper-point forest is maintainedtoallowtopology-preservingclusteringoperationswhicharedetailedinpre- vioussections. AsillustratedinFigure6.7,treesintheforestareconnectedcomponents that can be simplified via a series of clustering operations, and each of them is finally 90 Figure6.7: Hyper-pointclusterforestviewedfromobliqueandorthogonalperspectives represented by its root (blue and gold hyper-points) whose coordinates are determined byoptimizinga2.5DQEFcombininggeometricinformationfromleafpoints. This hyper-point cluster forest is built in a bottom-up manner. In a quadtree cell c composed of four leaf cells c 0;0 ,c 0;1 ,c 1;0 ,c 1;1 , assume each leaf cell has a set of cluster rootsthatareavailableforfurtherclustering(i.e.,withoutexceedingthegeometryerror tolerance or violating degenerate test). The simplification algorithm first traverses all the corner points in c that are shared by two of the four leaf cells. At each grid corner, four vertices in adjacent cells are connected via a surface polygon. The simplification algorithm retrieves the roots of these vertices and detects possible clustering operations based on surface component, i.e., 1→1 S , 1→2 S , or 1→m S . Similar approaches can be applied to minimal grid edges which exhibit roof layer gaps. They imply boundary- neighborshipsleadingtopossibleoperations 2→2 B and 2→m B . With possible clustering operation detected, they are sequentially tested against the geometry error tolerance and the degenerate test. Since the test sequence may affect themodelingquality,eachclusteringoperationisassignedwithapriority. Inparticular, high priority is assigned to 1→1 S and 2→2 B ; medium priority is assigned to 1→2 S and 1→m S ;andlowpriorityisassignedto 2→m B . Thereasonbehindthispriorityassignment 91 is that 1-hyper points and 2-hyper points are expected to be first clustered together to formmeaningfulgeometricpatterns(e.g.,roofridgesandstraightverticalwalls),before they are merged into key features with higher dimensional topology. As for cluster- ing operations with same priority, the test sequence is determined by the addition to quadraticerrorsinascendingorder. Polygon generation of 2.5D dual contouring is adapted to the hyper-point cluster forest in a straightforward manner. Considering the unsimplified polygonal model cre- atedfromtheuniformgrid, eachpointinthemodelisreplacedbytherootofitscluster (or the corresponding portion of the root if that has more layers than the leaf point). Numerous polygons become degenerate and are removed automatically, e.g., a triangle whose three vertices belong to the same cluster and thus map to the same root point. A simplepolygonalmodelwithsmallamountoftrianglesisproducedwhichhasthesame 2.5Dbuildingtopologyastheunsimplifiedmodel. 6.4 ExperimentalResults Figure 6.8 shows a stadium model reconstructed using different approaches, namely, 2.5Ddualcontouring(a,b),themodelingmethodpresentedinthischapter(c,d),manual creation (e), and the plane-based method presented in Section 3.4. In particular, the geometryerrortoleranceδisvariedinordertomaketrade-offsbetweenmodelscaleand fitting quality. The relation curve between δ and the number of triangles produced by differentapproachesisillustratedinFigure6.8(g). Quantitativemeasurementsaregiven inTable6.1. Witherrortoleranceδ increasing,thenewmethodconstantlydecreasesthe triangle number of reconstructed models. Reasonable cost is paid in fitting quality as the trade-off. On the contrary, 2.5D dual contouring reaches the simplification barrier around 3,000 triangles. This barrier can be explained by the last column of Table 6.1, 92 0.25 1 4 16 64 256 0 0.5 1 1.5 2 2.5 3 x 10 4 Geometry error tolerance # of triangles Manual creation baseline Plane-based model baseline 2.5D dual contouring Our method (a) (b) (c) (d) (e) (f) (g) δ = 4.0 Tri. # = 5191 δ = 4.0 Tri. # = 3290 δ = 64.0 Tri. # = 3315 δ = 64.0 Tri. # = 766 Tri. # = 1227 Tri. # = 6262 Figure 6.8: A stadium model created using (a,b) 2.5D dual contouring, (c,d) the mod- eling method presented in this chapter, (e) manual creation, and (f) the plane-based approach in Section 3.4. The relation curve between error tolerance and the triangle number of reconstructed models is illustrated in (g). With larger geometry error tol- erance given, the new method can always produce simpler models with less triangles; whiletheoverstricttopologytestin2.5Ddualcontouringcreatesnumeroustrivialtrian- glesalongthinrooffeaturesshownincloseupsof(a,b). ... δ 1.0 2.0 4.0 8.0 ∞ Tri. # = 169 Tri. # = 113 Tri. # = 80 Tri. # = 63 Tri. # = 52 ... Figure6.9: Modelevolutionwitherrortolerancegrowingfrom1.0toinfinite showing the percentage of unsuccessful collapses caused by topology test among all unsuccessful collapses. Many of them happen in small cells that create trivial triangles asshowninFigure6.8(a,b)closeups. Figure6.10showsthebuildingreconstructionfora5km-by-7kmurbanareaofDen- ver, from 73M input aerial LiDAR points with 6 samples/sq.m. resolution. The urban modeling pipeline in Section 4 is adopted to extract individual building patches, and 93 test the modeling method with topology control, 2.5D dual contouring, and a plane- based method independently. A fairly large error tolerance is utilized for both the new methodand2.5Ddualcontouring. Thenewmodelingmethodsuccessfullyreconstructs 2,0992.5Dbuildingmodelswithin5minutesonaconsumerlevellaptop(Inteli-7CPU 1.60GHz with 6GB memory). 227,566 triangles are produced for the building models which are rendered in the top rows of Figure 6.10. The output triangle number is com- parabletoplane-basedresults(181,752trianglesrenderedinthebottomrow). However, unliketheplane-basedmethod,themodelingapproachpresentedinthischapterdetects andpreserves2.5Dbuildingtopology,thusavoidsproducingcracksandinconsistencies between building blocks. E.g., the roof of the large structure shown in the bottom left closeup intersects with small features on top of it; while the new method does not have such problem. The middle row of Figure 6.10 shows 2.5D dual contouring result, it producestwiceasmanytriangles(551,341triangles)astheothertwoapproaches. Since the scale of the result is inversely proportional to geometry error tolerance δ. It is beneficial to study the evolution of building model with respect to δ. In particu- lar, 2.5D building models are created for the same LiDAR point cloud using an expo- nentially increasing δ, shown in Figure 6.9. Although the model geometry constantly becomes simpler, the building topology is faithfully preserved. Even in the extreme simplification case where δ = ∞, a model with 32 vertices and 52 triangles is cre- ated,whichcontainsthesmallestamountofverticesandtrianglesthatcanrepresentthe buildingtopologyofthismodel. 6.5 Appendix Proofofequivalencebetween(6.8)+(6.9)and(6.10): First,allthreeequationsstandtrueintheinitialunsimplifiedmodel. 94 Geometry error tolerance Thenewmethod 2.5Ddualcontouring Triangle# Avg. distance 2 Triangle# Avg. distance 2 Topologytest failurerate δ = 0.25 24161 0.0157 26776 0.0156 3.02%(82outof 2713) δ = 1.00 8864 0.0202 11280 0.0182 11.83%(109out of921) δ = 4.00 3290 0.0849 5191 0.0316 39.10%(149out of381) δ = 16.00 1374 0.1259 3644 0.0988 76.30%(161out of211) δ = 64.00 766 0.4812 3315 0.2104 90.45%(161out of178) Table 6.1: Quantitative comparison between modeling with topology control and 2.5D dual contouring using the experiment shown in Figure 6.8. The last column reports the ratio of cell collapses rejected by topology test among all rejected collapses. The topologytestbecomesdominantin2.5Ddualcontouringwithlargeerrortolerance. Afteraclusteringoperation,ifbothequation(6.8)and(6.9)aretrue,foreachR∈R, thereisf(P(R)) = f(R)≥ 1. I.e.,P(R)containsatleastone2-cell. Thus,theboundary ofP(R)containsatleast3edges,e(∂P(R))≥ 3. Conversely, when Equation (6.10) is true, P(R) has at least one 2-cell. Therefore, f(R) = f(P(R)) ≥ 1, i.e., (6.8). As for a wall feature W with|∂P(W)| ≥ 2, it is bounded bytwopointfeatures, and (6.9)istruebydefinitionsofhyper-pointclustering operations. Now considering a wall feature W with |∂P(W)| ≤ 1 (e.g., W 0 in Fig- ure6.5(d)andW 1 inFigure6.5(b)),P(W)isacloseloopon2Dspace. ItdividesRinto partitionR in andR out ,whereR in ̸=∅. Therefore: P(W) = ∂ ∑ R∈R in P(R). (6.11) Sincetherightpartsumsatleastone2-cellbeforeboundaryextraction,theboundaryof thecellcomplexcontainsatleast3edges, i.e.,e(P(W))≥ 3. 95 Figure 6.10: 2,099 building models are created for an urban area in Denver using (top) building modeling with topology control, (middle) 2.5D dual contouring, and (bottom) plane-based method. The new method produces as few triangles as the plane-based methodwhilerecoveringandpreservingthetopologicalfeaturesineachbuildingstruc- ture. 96 Chapter7 2.5DBuildingModelingbyDiscovering GlobalRegularities Besides data-driven building modeling approaches such as those proposed in Chap- ter 5 and Chapter 6, another popular strategy in attacking the 2.5D building modeling problem is to introduce primitives (e.g., planes, spheres, cones) to represent building shapes. In particular, planes receive the most attention since they are the common- est structures in man-made objects, especially in buildings. Planar roof patches are locally fitted from raw points, and are later combined with vertical facades aligning with roof boundaries, to construct a compact mesh model while maintaining low geo- metric fitting error rate. The main difficulty of this strategy is that local plane fitting can become unstable when dealing with noisy or incomplete point clouds. Artifactsare inevitably created from unreliable plane primitives. To alleviate this problem, existing methods typically introduce strong urban priors to prune the fitted planes, such as roof topology[60],Manhattan-worldgrammars[41,48],andprincipaldirectionsintroduced in Section 3.4.3 and Section 4.3.4. While prior knowledge successfully increases the robustness of these methods, it tends to be overstrict and thus limits their applicability when dealing with moderately complex building structures such as the one shown in Figure7.1. This chapter proposes global regularities, a moderate yet informative structure to organize planar roof patches and roof boundary segments. As illustrated in Figure 7.1 bottom,Iexplorebothorientationandplacementregularitiesbetweenplanarroofpatch 97 Aerial imagery 2.5D point cloud Final model View 1 View 2 View 1 View 2 116 o 116 o 107 o 107 o 107 o 116 o 130 o 130 o 11.9 6.4 6.4 5.6 5.6 90 o 90 o 8.7 8.7 6.4 11.9 11.9 90 o 90 o 90 o 5.6 Figure7.1: Themethodpresentedinthischapterautomaticallydiscoversglobalregular- itiesfromanoisy2.5Dpointcloud,andusesthemtocreateaconvincing2.5Dbuilding model. Two orthogonal projections illustrate a subset of the global regularities in this model(lengthsinmeters). pairs (e.g., slope angle equality), between roof boundary segment pairs (e.g., segment height equality), and between a planar roof patch and its boundary segments (e.g., orthogonality between their orientations). These global regularity patterns reveal the inter-elementsimilaritiesandrelationsthatintrinsicallyexistinbuildingmodelsbecause of human design and construction. With these patterns, the complexity of the building modelingproblemcanbesignificantlyreducedforcomplicatedbuildingmodelssuchas theoneinFigure7.1. 98 This chapter presents an automatic algorithm which detects global regularities and utilizesthemtocalibrateplaneprimitives. Unlikethestrongpriorsintroducedbyprevi- ousmethods,globalregularitiesofferamoreflexiblerepresentationoftheglobalknowl- edge in 2.5D building models, and thus enable the modeling algorithm to handle more complicatedbuildingshapes. Another significant advantage of global regularities is that they characterize the intrinsicstructuresofbuildingmodels,towhichhumanvisionissensitive. Forinstance, Figure 7.2 right shows two models created from plane primitives. Without comparing modelgeometrywithinput points, humanvisionimmediatelyfindsthetop-rightmodel moreconvincingsinceitconformstomoreglobalregularities. In this chapter, Section 7.1 summarizes shape-from-symmetry approaches. Sec- tion 7.2 explores various global regularity patterns in 2.5D building structures. Sec- tion 7.3 presents an automatic algorithm to discover and enforce global regularities through a series of alignment steps. Finally, Section 7.4 shows 2.5D modeling results withhighqualityintermsofbothgeometryandhumanjudgement. 7.1 ReviewofShape-from-SymmetryApproaches Inbothcomputervisionandcomputergraphics,symmetryhasbeenidentifiedasreliable global knowledge in 3D reconstruction. For instance, Fisher [16] introduces domain knowledge of standard shapes and relationships into reverse engineering problems. Thrun and Wegbreit [57] detect symmetries and utilize them to extend partial 3D mod- els into occluded space. Gal et al. [19] adopt 1D wires to carry both local geometry information and global mutual relationships in man-made objects. Li et al. [36] extract relationshipgraphsamongprimitivestoencodeintra-partrelationsandusethemtofur- therimprovethereconstructionquality. 99 Manual creation Our result 2.5D dual contouring Primitive- based method Figure 7.2: Modeling results generated from the same input point cloud by manual creation,themodelingmethodproposedinthischapter,2.5Ddualcontouringwithprin- cipaldirectionsnapping,andaprimitive-basedmethod[34]. Thenewmodelingmethod creates the most visually convincing result among all three automatic methods since its resultingbuildingmodelconformstothemostglobalregularities. These methods are similar in spirit to the method presented in this chapter. While previous research work focuses on 3D man-made objects, this research is the first to exploreglobalregularitiesin2.5Dbuildingmodelscomposedofroofpatchesandverti- calwalls. 100 x y z θ(n1) φ(n1) Ridge Boundary segments r1,2 b1 b2 b3 b4 b5 b6 n1 p1 p2 n2 n1 n2 P1 P2 Figure 7.3: A typical gable-shaped building roof containing a set of 2.5D elements (e.g.,planeprimitive,boundarysegments,andridges) 7.2 GlobalRegularities In 2.5D building models, global regularities characterize the inter-element similarities and relations arising from human design and construction. They are particularly use- ful in correcting plane primitives and creating more visually convincing building mod- els. This section explores various global regularity patterns that commonly exist in 2.5D building models, and demonstrates them using a typical gable-shaped building roofshowninFigure7.3. Considering a 2.5D building model composed of plane primitives including planar roof patches and planar facade patches, it can be equivalently represented by a set of planarroofpatchestogetherwiththeirboundarysegments;becausegiventhe2.5Dcon- straintsthatroofsurfacesarealwaysboundedbyverticalfacades,planarfacadepatches andlinearroofboundarysegmentshavethesameprojectiononthex-yplane. Inpartic- ular,theplanarroofpatchsetisdenotedasP ={P i : (v−p i )·n i = 0}inwhicheach plane P i is determined by a normal-position pair (n i ,p i ). A boundary segment set for eachplanarroofpatchiscollectedbyintersectingP i withitssurroundingplanarfacade patches, denoted asB i (e.g., in Figure 7.3,B 1 = {b 1 ,b 2 ,b 3 }). The following sections 101 explore global regularities among the 2.5D element setP∪(∪ i B i ) from three aspects: roof-roofregularities,roof-boundaryregularities,andboundary-boundaryregularities. 7.2.1 Roof-RoofRegularities Thissectionfocusesontwoclassesofcommonlyencounteredregularitiesbetweenroof planepair(P i ,P j )asfollows. 7.2.1.1 Orientationregularities In3Dmodels,theorientationregularitiesareusuallydefinedastheorthogonalityorpar- allelismbetweenplanenormals(e.g.,[36]). Thisdefinition,however,cannotbedirectly applied to 2.5D building models for two reasons: first, roof plane normals rarely show orthogonality or parallelism; second, roof inclination and direction are of more interest inbuildingmodeling. In2.5Dmodels,orientationregularitiesarenotdeterminedbythe angle between plane normals, but by the projections of normals on either the x-y plane orthez-axis. Forinstance,althoughn 1 andn 2 inFigure7.3donotexhibitorthogonality orparallelism,theyshowstrongorientationregularitiessincetheirprojectionsonthex- yplaneareopposite. Therefore,Ichoosetowriteplanenormalsinsphericalcoordinates (θ(n),φ(n)): θ(n) = arccos(n z ), (7.1) φ(n) = arctan(n y ,n x ), (7.2) whereθ(n)∈ [0,π/2)andφ(n)∈ [0,2π)(Figure7.3right). Intuitively, θ(n) determines the inclination of the planar roof patch, and φ(n) indi- cates the direction of the slope. As humans are particularly interested in roof patches having the same inclination and roof patches exhibiting regularized slope directions 102 (eitherorthogonal or parallel), fourtypical roof-rooforientationregularitiesaredefined accordingly: • θ-equalitywhenθ(n i ) = θ(n j ), • φ-equalitywhenφ(n i ) =φ(n j ), • φ-oppositewhenφ(n i ) =φ(n j )±π, • φ-orthogonalitywhenφ(n i ) = φ(n j )± 2 , 3 2 . For example, plane pair (P 1 ,P 2 ) in Figure 7.3 exhibits the same inclination and the oppositeslopedirection. Usingtheaboveformulation,thesecharacteristicsaredenoted asθ-equalityandφ-oppositerespectively. 7.2.1.2 Placementregularities Placement of roof planes (i.e., roof positions) by themselves do not contain much reg- ularity information. However, the placement of intersections between roof plane pairs maycarrymeaningfulstructuralinformationaboutthebuilding. Inparticular,ridgesare definedtorevealtheregularitiesofroofplaneplacements. Ridge definition: for a neighboring plane pair (P i ,P j ) satisfying both θ-equality and φ-opposite,theintersectionofP i andP j isdefinedasa ridge,denotedasr i;j . Thedirectionofridger i;j isuniquelydeterminedas d(r i;j ) = (sin(φ(n i )),−cos(φ(n i )),0) T . (7.3) Since d(r i;j ) is parallel to the x-y plane, the placement of r i;j can be parameterized by a pair of real numbers (h(r i;j ),p(r i;j )), denoting the height of r i;j and the distance from origin tor i;j ’s projection on the x-y plane respectively. They can be calculated by 103 solving an equation system with plane equations regarding (n i ,p i ) and (n j ,p j ). Two typesofplacementregularitiesforridgepair(r i;j ,r k;l )aredefinedas: • Ridge-height-equalitywhenh(r i;j ) = h(r k;l ), • Ridge-position-equalitywhend(r i;j )∥d(r k;l )andp(r i;j ) =p(r k;l ). 7.2.2 Roof-BoundaryRegularities I observe that the majority of boundary segments are aligned either orthogonally (e.g., b 2 in Figure 7.3) or parallel (e.g., b 1 ,b 3 ) to the normals of their owner planes (e.g., P 1 ), when projected on the x-y plane. Therefore, this section focuses on roof-boundary regularities between planeP i and its boundary segmentsB i . The direction ofP i on the x-yplaneisdenotedbya2Dvectoro i = (cos(φ(n i )),sin(φ(n i ))) T ,andthedirectionof a boundary segmentb j ’s x-y projection is denoted aso(b j ) = P(d(b j )),b j ∈B i , where Pistheprojectionoperatorandd(b j )istheb i ’sdirectionin3Dspace. Thereare: • o-parallelismwheno i ∥o(b j ), • o-orthogonalitywheno i ⊥o(b j ), 7.2.3 Boundary-BoundaryRegularities As the orientation regularities among boundary segments can be implied from roof- boundary regularities and roof-roof orientation regularities, this section focuses on placement regularities between boundary segment pairs. In particular, two significant regularity patterns are noted: first, boundary segments that are parallel to the x-y plane may have similar heights (e.g., (b 2 ,b 5 ) in Figure 7.3); second, when projected onto the x-y plane, boundary segments with the same direction may align to the same line (e.g.,(b 1 ,b 6 )and(b 3 ,b 4 )). Thus,boundary-boundaryregularitiesaredefinedasfollows, 104 Coarse-to-fine iteration φ h h Orientation alignment Plane fitting Placement alignment Boundary segment initialization Segment height alignment Segment position alignment Pruning and model generation hs hs p p p p Figure 7.4: Pipeline of the modeling approach: a 2.5D point cloud (top left) is trans- formed to a building model (bottom left) through a series of steps. Global regularities arediscoveredandenforcedineachalignmentstep. where h(b i ) is the height of boundary segment b i , and p(b i ) is the distance from origin tob i ’sprojectiononthex-yplane: • Segment-height-equality whenh(b i ) = h(b j ), and bothd(b i ) andd(b j ) are par- alleltothex-yplane, • Segment-position-equalitywheno(b i )∥o(b j )andp(b i ) = p(b j ). 7.3 ModelingwithGlobalRegularities Given a noisy 2.5D point cloud as input, this section presents an automatic method to simultaneously detect locally fitted plane primitives and global regularities. In general, a discover-then-align strategy is adopted: once initial plane primitives are identified, the modeling algorithm discovers global regularities from them, and then immediately 105 refinestheseinitialprimitivesbyaligningthemtotheglobalregularities. Thisoptimiza- tion strategy is applied individually to each type of global regularities defined in Sec- tion 7.2. It effectively corrects the geometric errors raised by local fitting approaches, andthussignificantlyimprovesthemodelquality. AnoverviewofthisapproachisshowninFigure7.4. Thebuildingmodelingsystem containsthreemainmodulestocreatea2.5Dbuildingmodel(bottomleft)fromanoisy aerialscan(topleft): 1. Planarroofpatchextraction: AsshowninFigure7.4top,withplaneprimitives detected via local fitting, two discover-and-align steps are sequentially executed to detect the roof-roof regularities and refine the planar roof patch, namely, ori- entation alignment and placement alignment. Both planar roof patches and the roof-roofregularitiesareiterativelygeneratedinacoarse-to-finemanner. 2. Boundary segment production: The modeling algorithm immediately enforces the roof-boundary regularities by creating a rectangular bounding box for each planar roof patch, and identifying boundary segments from bounding box edges, shown as the black lines in Figure 7.4 bottom. These boundary segments are furtherrefinedbydiscoveringandenforcingboundary-boundaryregularities. 3. Model generation: Vertical facades are automatically generated from boundary segments to connect roof patches and the ground. Rectangular roof patches are prunedbyneighboringelements. A2.5Dbuildingmodelisproducedbycombin- ing both planar roof patches and vertical facades as shown in Figure 7.4 bottom left. 106 7.3.1 PlanarRoofPatchExtraction Given a set of points equipped with normals 1 , a popular plane detection algorithm is adopted for aerial LiDAR scans [34, 60] to find plane primitives: a region-growing procedure is applied to find spatially connected point clusters with similar normals; then plane primitives are locally fitted to individual point clusters. The detected plane primitive set is denoted asP ={P i }. Orientation alignment and placement alignment areappliedtoP ={P i }sequentially. 7.3.1.1 Orientationalignment Byexpressingplanenormalsinsphericalcoordinates,theorientationregularitiescanbe categorizedintotwoclasses: θ-equalityfindsplaneswithsimilarslopeangles,whileφ- equality,φ-opposite andφ-orthogonality showregularized roof patch directions. These orientation regularities can be discovered by detecting clusters of = {θ(n i )} and clusters of = {φ(n i ) mod (π/2)} respectively. Each angle cluster implies a set of corresponding orientation regularities while the center of each cluster predicts the best alignment. In particular, the complete-linkage clustering algorithm [12] is adopted to identifyclustersinand. ClustercentersetsC andC aretakenasconstraintsinthe subsequentalignmentstage,inwhichθ(n i )andφ(n i )aresnappedtothecorresponding clustercentersinC andC . 2 1 Normalscanbeeffectivelyestimatedviacovarianceanalysis[47]. 2 In singular cases whereθ(n i )≈ 0,φ(n i ) becomes unstable. Thus,θ(n i ) is snapped to 0 andφ(n i ) isassignedtothemostpopularφ∈C . 107 7.3.1.2 Placementalignment To effectively deal with placement regularities, ridges are detected from neighbor- ing plane pairs. Similar to orientation alignment, the placement alignment algo- rithm decouples the placement alignment into two independent sub-problems: align- ing ridge heights towards ridge-height-equality and aligning ridge positions towards ridge-position-equality. Theseplacementregularitiescanbediscoveredbyfindingclus- ters of ridge height set H = {h(r i;j )}; and clusters of ridge position set S(d) = {p(r i;j )|d(r i;j ) ∥ d}, regarding each ridge directiond. Cluster center sets are denoted asC H andC S respectively, and used as regularity constraints henceforth. In the align- mentstage,ridgeheighth(r i;j )andpositionp(r i;j )arebothalignedtotheirclustercen- ters, resulting in modifications on plane positionp i andp j . In order to avoid conflicts betweenridgeheightalignmentandridgepositionalignment,theformeronlyaffectsthe z values of position vectors, while the latter makes modifications to the x and y coordi- nates. Therefore,theonlyconflictsourceliesinplanesthathavemultipleridges,where theridgescompeteinmodifyingtheplane’sposition. Inthiscase,onlythelongestridge isallowedtomodifytheposition. Theeffectsfromotherridgesareignored. 7.3.1.3 Coarse-to-fineiteration The planar roof patch extraction executes in a coarse-to-fine manner, as demonstrated in Figure 7.5. In particular, the modeling system fits planes to the input points, makes orientation alignment and placement alignment to the plane primitives, discards points already associated with existing plane primitives, and then iterates through these steps withthreemodificationsuntilnomoreplaneprimitivescanbefoundbyplaningfitting: 1. Plane-fitting criterion is loosened to accept smaller plane patches. In particular, normal variance allowance α is increased and the minimum number of points 108 Iteration 1 Iteration 3 Iteration 2 Figure7.5: Coarse-to-fineplanarroofpatchextraction requiredtovalidateaplaneprimitiveK isreduced. Empirically,αremainspi/12 whileK isreducedbyhalfineachiteration,from200to25. 2. In plane-fitting, instead of randomly picking region-growing seed, the plane detection algorithm starts from points with normals that are close to normal n(θ,φ),∀(θ,φ)∈C ×C . These points have great potential to grow into plane primitivesconformingtoexistingorientationregularityconstraints. 3. In alignment steps, the new plane primitives first attempt to snap to existing con- straintsets{C ,C ,C H ,C S }. Clusteringisperformedonlyforprimitivesthatcan- notaligntotheexistingconstraintsgivenadistancethreshold. Newclustercenters arecombinedwithexistingcenterstoformregularityconstraintsforthenextiter- ation. These modifications are driven by two observations: first, large plane primitives are morereliableinproducingglobalregularityconstraints,andthustheiterativealgorithm begins with large primitives and requires small primitive to be snapped to large primi- tives if possible (modification 3); second, the regularity constraints detected from large primitivescangreatlyimprovetherobustnessofplanefittingforsmallpatches(modifi- cation2). 109 7.3.2 Boundarysegmentproduction Given a set of plane primitives as shown in Figure 7.4 top right, initial boundary seg- ments are created for each planar roof patch with help of roof-boundary regularities. Segment height alignment and segment position alignment are performed based on boundary-boundaryregularities. 7.3.2.1 Boundarysegmentinitialization As discussed in Section 7.2.2, most boundary segments in 2.5D building models con- form to the roof-boundary regularities, i.e., when projected on the x-y plane, boundary segments are either orthogonal or parallel to the normals of their owner planes. On the other hand, boundary segments represent the borders of roof patches, and thus bound the patch content, i.e., points associated with the roof plane. Considering a planar roof patch P i and the set of points associated with it denoted asV i ; a boundary extraction algorithm computes a rectangularbounding boxR i ofpointprojectionsP(V i )onthex- y plane, with R i ’s orientation followingo i = (cos(φ(n i )),sin(φ(n i ))) T . The segment setB i isinitializedbybackprojectingR i ’sedgesontoP i . Given that the plane normals are already aligned to a small set of orientation con- straintsC ,thex-ydirectionsofboundarysegmentsfallintoafew2Ddirections,which can be regarded as building-scale principal directions (as discussed in Section 3.4.3). In the special case where|C | = 1, the boundary segments are in two orthogonal x-y directions,forming rectilinear contoursforbuildingmodels[41]. 7.3.2.2 Segmentheightalignment The modeling system now applies segment height alignment to horizontal boundary segments. Segment-height-equality is discovered and enforced by a clustering method similartopreviousalignmentstepswithtwoadditionalrules: 110 1. Heights of boundary segments from the same plane patch cannot belong to the same cluster, since snapping them together risks making the roof patch degener- ate. 2. Foreachplanepair(P i ,P j )thatcreatesaridger i;j ,theboundarysegmentsoppo- site tor i;j (e.g.,b 2 andb 5 in Figure 7.3) are tested with arelaxedcriteria, because eachridgeindicatesahighprobabilityofareflection-symmetry. Aboundarysegmentismarkedas“fixed”ifitsheightissnappedtoaclusterwithmore than one element. The position of each fixed segment is determined accordingly by substituting the modified height value into the plane equation. These positions act as placementconstraintsC p inthenextstep. 7.3.2.3 Segmentpositionalignment Segment position alignment is applied to the remaining boundary segments including bothnon-horizontalsegmentsandhorizontalsegmentsthatarenotmarkedas“fixed”in the previous step. Positions of these segments first attempt to snap to elements inC p , andonlyjointhepositionclusteringprocedurewhenthesnappingattemptfails. 7.3.3 ModelGeneration With planar roof patches and their boundary segments generated and aligned to the global regularities, a 2.5D building model is reconstructed by combining roof patches and vertical facades, as shown in Figure 7.4 bottom left. The facades are produced from boundary segments connecting roof patches to the ground, while rectangular roof patchesareprunedbyneighboringelementsincludingbothroofpatchesandfacades. 111 Metric Thismethod 2.5Ddual contouring Primitive-based method[34] Manual creation Triangle# 130 214 89 78 Avg. dis 2 0.012 0.016 0.018 0.058 Outlierratio 0.9‰ 10.0‰ 11.1‰ 17.3‰ Table7.1: QuantitativeresultsforthecomparisoninFigure7.6 7.4 ExperimentalResults This section first compares the modeling method proposed in this chapter with two existing approaches (i.e., 2.5D dual contouring and a primitive-based method [34]). Figure 7.2 shows a qualitative comparison between these methods. The modeling result with global regularities has the most similar visual appearance to manual cre- ation, because it conforms to the most global regularities that characterize the intrinsic structure of building models. In contrast, 2.5D dual contouring only considers bound- ary direction similarities by introducing the principal direction snapping algorithm, while [34] detects primitives but does not deal with the relations between them. In addition, quantitative comparison between the three methods are made using metrics shown in Figure 7.6 and Table 7.1. While the modeling result with global regularities shows comparable triangle number and average squared distance compared with pre- vious approaches, it significantly reduces the outlier ratio (i.e., the ratio of points with squared distances greater than 0.25m 2 ), because the modeling method can robustly fits small plane primitives through iterations with global regularities detected and updated progressively. The modeling method with global regularities is further tested on several LiDAR scans of buildings in the city of Atlanta, as illustrated in Figure 7.7. The input contains aerial LiDAR point cloud with 17 samples/m 2 resolution (first column). Planar roof patches and their boundary segments are detected and aligned with global regularities 112 Squared distance > 0.25 0.0 Manual creation Our result 2.5D dual contouring Primitive-based method Figure 7.6: A comparison of geometric fitting errors (i.e., squared distances from input pointstothemodelsurfaces)betweenmodelscreatedbyfourdifferentapproaches learnt from them (second column). 2.5D building models are reconstructed from these primitives by pruning the roof patches (see examples in the last three rows) and creat- ing vertical facades from boundary segments (third column). Given the aerial imagery as reference (fourth column), the new results are more “realistic” than the 2.5D dual contouring results (last column). Table 7.2 shows the statistics of the experiments in Figure7.7. Computationtime ismeasured onalaptopwith Inteli-7CPU1.60GHzand 6GBmemory. 113 Model Point # Plane # |C | |C | Ridge # |C H | |C P | Hor. segment # Seg. h-cluster # Non-hor. segment # Seg. p-cluster # Tri. # time (s) 1strow 3349 6 1 1 3 3 2 6 1 12 5 48 7.1 2ndrow 5603 6 2 1 3 3 2 6 2 12 6 48 6.8 3rdrow 4549 4 1 2 2 1 2 4 1 8 4 28 3.6 4throw 6143 13 3 1 5 3 3 18 7 24 13 110 16.6 5throw 6920 14 3 1 7 3 3 16 2 28 10 112 39.6 Table 7.2: Statistics of the experiments in Figure 7.7. Column 3 - 12 show statistics of the intermediate results. Since the size of constraint sets is considerably smaller than the number of primitives (in bold), the solution space (and thus the complexity) ofthemodelingproblemisreducedsignificantly. 114 Input point cloud Plane primitives aligned to global regularities Our modeling result Aerial imagery 2.5D dual contouring result for comparison Figure 7.7: Experiments on several building scans. By discovering global regularities and enforcing them on the planar roof patches and their boundary segments (second column), the modeling algorithm creates visually convincing 2.5D building modelings (third column). Aerial imagery and 2.5D dual contouring results are included as refer- ence. 115 Chapter8 ModelingResidentialUrbanAreas (a) Input aerial point cloud (c) Aerial imagery as a reference (b) Our modeling result Figure 8.1: Given (a) a dense aerial LiDAR scan of a residential area (point intensities represent heights), the residential urban modeling system reconstructs (b) 3D geometry forbuildingsandtreesrespectively. (c)Aerialimageryisshownasareference. Twonewchallengesemergewhentheurbanmodelingproblemextendstoresidential areas. First,asshowninFigure8.1(a),vegetationisamajorcomponentofurbanreality in residential areas. An urban modeling method for residential areas should detect and reconstruct both buildings and trees, e.g., as in Figure 8.1(b). The second challenge lies in the classification method: dense LiDAR scans capture the detailed geometry of tree crowns, which may have similar height and local geometry features as rooftops of residential buildings. Figure 8.2 shows such an example where part of the tree crown shows similar or even better planarity than part of the rooftop (see closeups illustrating local points as spheres together with the optimal plane fitted to them). Classification algorithms based on local geometry features may fail and produce significant modeling errors. E.g.,Figure8.2right. 116 A B Figure 8.2: Local geometry features become unreliable when dealing with residential areas with rich vegetation. In closeups of (A) a tree crown region and (B) a rooftop, pointsarerenderedassphereswhilealocallyfittedplaneisrenderedinyellow. Middle: classification results from the general urban modeling system in Chapter 3, trees in green,buildingsinpurple,andgroundindarkgrey. Right: modelingartifactsarecreated becauseofclassificationerrors. Toaddressthesetwochallenges,thischapterpresentsarobustclassificationmethod to classify input points into trees, buildings, and ground. Building models and trees are createdfromthesepointsusing2.5Ddualcontouringandanovelleaf-basedtreemodel- ing approach, respectively. The heart of the classification method is a simple, intuitive, but extremely effective measurement. In particular, I observe that residential buildings usually show a strong 2.5D characteristic, i.e., they are composed of skywards roofs and vertical walls; both are opaque and thus prevent the laser beams from penetrating the building structure. Therefore, there is no point sample inside the building structure. The rooftops (or ground) become the lowest visible surface at a certain x-y position, as illustratedinFigure8.3left. Incontrast,trees,composedofbranchesandleaves,donot have this 2.5D structure. With multiple passes of scanning from different angles, the point cloud captures not only the top surface of the tree crown, but also surfaces inside andunderneaththecrown,asshowninFigure8.3right. In this chapter, Section 8.1 reviews tree detection and tree modeling approaches respectively. Section8.2proposesaneffectivealgorithmtoclassifytrees,buildingroofs, and ground. In Section 8.3, a hybrid model containing both 2.5D building models and 117 Second pass First pass First pass Second pass Point samples Building Tree Ground Ground Figure 8.3: While building structures have a 2.5D characteristic, trees do not possess suchproperty. Denselaserscansmaycapturesurfacepointsunderthetreecrown(right). leaf-based tree models is generated in an automatic and robust manner. Experimental resultsareshowninSection8.4. 8.1 RelatedWork 8.1.1 TreeDetectioninLiDAR Inurbanmodelingsystems,treesareoftenrecognizedasoutliersandthusareclassified and removed in the first step. Most of the classification algorithms rely on point-wise featuresincludingheight[34,39,52,58]anditsvariation[5,39,52],intensity[39,52], and local geometry information such as planarity [34, 60], scatter [34, 58], and other local geometry features. Heuristics or machine learning algorithms are introduced as classifiers based on the defined feature set. To further identify individual building roof patches,segmentationiseitherintroducedinapost-classificationstep,orcombinedwith classificationintheformofenergyminimizationsuchas[34]. Nevertheless, the method proposed in this chapter is the first to introduce the 2.5D characteristic of building structures into the classification problem. A simple, efficient 118 and effective classification algorithm is presented, which gains great accuracy in resi- dentialareaswithrichvegetation. 8.1.2 TreeModeling Tree modeling is a missing part in most of the aerial LiDAR based urban modeling approaches. To the best of my knowledge, Lafarge and Mallet [34] is the only research work which addresses the tree modeling problem by matching simple ellipsoidal tem- plate to tree clusters. This method, however, is problematic when dealing with compli- catedtreestructuresinresidentialareas, e.g.,Figure8.1(a). Computergraphicsandremotesensingcommunitieshavemadegreateffortsinmod- eling trees from ground LiDAR and imagery, such as [9, 38, 44, 55, 56, 63]. A general tree model is broadly adopted in these literatures, composed of skeletal branches and leaves attached to them. Inspired by these efforts, this chapter proposes leaf-based tree modelingfromaerialLiDARscans. 8.2 PointCloudClassification Given an aerial LiDAR point cloud of a residential area as input, the objective of clas- sification is to classify points into three categories: trees, buildings, and ground. As mentionedpreviouslyandillustratedinFigure8.3,the2.5Dcharacteristicisthekeydif- ference between trees and buildings (or ground). In order to formulate this concept, the pointcloudisfirstembeddedintoauniform2DgridG. Ineachgridcellc,thepointset P(c)issegmentedintomultiplelayerfragmentsL(c),usinglocaldistance-basedregion growing. Ideally, a layer fragment l building ∈ L(c) lying on a 2.5D object (rooftop or ground) must have the lowest height among all layer fragments in L(c), because the 119 rooftop (or ground) is always the lowest visible surface to laser beams at a certain x- y position, as analyzed previously. On the other hand, a tree layer fragment l tree can exhibit any height. However, as there is usually a ground or rooftop surface underneath tree samples, l tree is not expected to be the lowest layer fragment in L(c). From an energy minimization perspective, these features can be described by a data energy term E d (x l )foreachl∈L(c)as: E d (x l ) =            α ifx l = buildingor ground,andl isnotthelowestinL(c) β ifx l = tree,andl isthelowestlayerfragmentinL(c) 0 otherwise (8.1) wherex l isthelabeloflayerfragmentl. Tofurtherdiscriminatebuildingandgroundintheenergyminimizationframework, elevation of layer fragment e(l) is introduced and defined as the height difference between l and the ground elevation at c’s center [34, 58] 1 . Another data energy term E g (x l )isdefinedaccordingly: E g (x l ) =            γ·max(1− e(l) ,0) ifx l =building γ·min( e(l) ,1) ifx l =ground 0 ifx l =tree (8.2) whereσ isthenormalizationfactor. Empirically,σ = 6m,assuggestedin[34]. 1 Thegroundelevationmapcanbeeasilyestimatedbyassigninga20m-by-20mcoarsegrid,estimating groundheightwiththelowestpointineachcell,andapplyinglinearinterpolationacrosstheentirecoarse grid. 120 WithasmoothenergyE s (x l 1 ,x l 2 )definedoverallneighboringlayerfragmentpairs (i.e., layer fragments belonging to neighboring cells and satisfying certain distance cri- teria), a Markov Random Field can be built, which leads to an energy minimization problemoverthelabelingxoftheentirelayerfragmentsetL: E(x) = ∑ l∈L (E d (x l )+E g (x l ))+λ ∑ (l 1 ;l 2 )∈N E s (x l 1 ,x l 2 ) (8.3) whereN is the set of neighboring layer fragment pairs, and smooth energyE s (x l 1 ,x l 2 ) isdefinedascharacteristicfunction1 x l 1 ̸=x l 2 . Withtheenergyminimizationproblembeingsolvedusingthewell-knowngraph-cut method[3],pointlabelsaredeterminedasthelabelofthecorrespondinglayerfragment. To further construct roof patches from building points, a region growing algorithm is applied based on certain distance criteria. While large building patches are adopted as rooftops,smallpatchesareconsideredasoutliersandremovedhenceforth. Figure8.4demonstratestheentireprocessofpointcloudclassification. Inputpoints are first discretized into layer fragments. For illustration purpose, Figure 8.4(b) shows onlythelowestlayerfragmentineachcell. Theyfaithfullycapturetheskywardsurfaces of 2.5D structures. By solving an energy minimization problem, building and ground layerfragmentsaredetectedandshowninFigure8.4(c). Becauseofthesmoothenergy term, eaves are segmented as part of the building roof, and small clusters of tree layer fragments with low heights are correctly detected, as illustrated in the closeups respec- tively. Finally,theclassificationresultisappliedtothepointcloudandaregiongrowing algorithmsuccessfullygroupsroofpatchesfrombuildingpoints, i.e.,Figure8.4(d). 121 Layer fragments representation (a) Input points (b) Lowest layer fragments in cells (c) Building and ground layer fragments (tree layer fragments not shown for clarity) (d) Points with classification results Figure 8.4: A demonstration of the classification algorithm: (b) the lowest layer frag- mentsfaithfullycapturetheskywardsurfacesof2.5Dstructures;(c)buildingandground layer fragments are rendered in purple and grey respectively; (d) trees and outliers are inblackwhilebuildingroofpatchesarerenderedinbrightcolors. 8.3 ModelingofUrbanElements Basedonthesuccessfulclassificationofinputpoints,differentmodelingapproachesare introducedfortrees,buildings,andgroundrespectively. 8.3.1 TreeModeling Modern tree modeling approaches adopt a general tree structure composed of skele- tal branches and leaves attached to them. Tree reconstruction usually begins with a 122 branch generation algorithm followed by a leaf modeling approach. However, unlike ground-based laser scans and imagery, aerial LiDAR date captures very few samples on branches, making branch generation a difficult task. Therefore, I choose to directly model tree leaves by fitting surface shapes around tree points having sufficient neigh- bors. In particular, for each tree point p with sufficient neighbors, Principal Component Analysis is applied to its neighboring point set N(p) to fit an ellipsoid. Eigenvectors v 0 ,v 1 ,v 2 and eigenvalues λ 0 ,λ 1 ,λ 2 of the covariance matrix represent the axes direc- tionsandlengthsoftheellipsoidrespectively. Theinscribedoctahedronoftheellipsoid isthenadoptedtorepresentthelocalleafshapearoundp. Specifically,anoctahedronis createdwithsixverticeslocatedat{v p ±sλ 0 v 0 ,v p ±sλ 1 v 1 ,v p ±sλ 2 v 2 },wherev p is thelocationofpandsisauser-givensizeparameter. AuniformsamplingoverthetreepointsetP tree canbeappliedtofurtherreducethe scaleofthereconstructedmodels. 8.3.2 BuildingModeling 2.5D dual contouring detailed in Chapter 5 is adopted to create building models from rooftoppatchesthroughthreesteps: (1)sampling2.5DHermitedataoverauniform2D grid,(2)estimatingahyper-pointineachgridcell,and(3)generatingpolygons. The only challenge in applying 2.5D dual contouring to residential area data lies in rooftop holes caused by occlusion. To solve this problem, a hole-filling step is added rightafter2.5DHermitedataissampledfrominputpoints. Inparticular,thehole-filling step scans the entire 2D grid to detect rooftop holes, and solves a Laplace’s equation ▽ 2 z = 0 to fill these holes, where z represents the heights of surface Hermite samples 123 at grid corners 2 . Existing surface Hermite samples serve as the boundary condition of theLaplace’sequation. 8.3.3 GroundModeling Ground models can be easily created by rasterizing ground points into a DSM (digital surfacemodel). Holesarefilledvialinearinterpolation. 8.4 ExperimentalResults The residential urban modeling system is tested on various data sets. For each data set, thefollowingparameterconfigurationisadoptedwithrespecttothedataresolution. The neighborhood radius r is set to 3 √ d given d as the point density in m −2 , to ensure suffi- cient samples in a point’s neighborhood. Octahedron sizes is chosen by the user in the interval[ 1 r , 3 r ]. Energyfunctionparameters,i.e.,{α,β,γ,λ}aresetto{1.0,2.0,0.5,4.0} empirically. Thisparameterconfigurationworkswellforallthedatasetsthathavebeen tested. Figure 8.6 shows the urban reconstruction results for a 520m-by-460m residential area in the city of Atlanta. The input contains 5.5M aerial LiDAR points with 22.9/m 2 resolution. The modeling system reconstructs 56K triangles for building models, and 53K octahedrons as tree leaves, in less than two minutes on a consumer-level laptop. As illustrated in the closeups of Figure 8.6, the classification algorithm successfully classifies points into trees, ground, and individual building patches (second column). A hybrid urban model is generated by combining 2.5D polygonal building models and 2 SurfaceHermitedataissampledpergridcorner,byintersectingaverticallineandtherooftopsurface. SeeSection5.3.1fordetails. 124 Atlanta area #2 Denver area Figure 8.5: The residential urban modeling system is tested against multiple data sets fromdifferentsources. Thesystemrobustlyreconstructsurbanrealitydespitethevaria- tionindataresolution,buildingpatterns,andtreetypes. leaf-based tree models (third column). Aerial imagery is given in the last column as a reference. The residential urban modeling system is then tested on another two data sets, with dataresolutionrangingfrom6/m 2 to19/m 2 . Visuallyappealingurbanmodelsarerecon- structedrespectively,despitethevariationinpointdensity,buildingmodelpatterns,and tree types, as shown in Figure 8.5. To quantitatively evaluate the modeling results, the falsepositives(unexpectedresults)andfalsenegatives(missingresults)ofbuildingsare counted by comparing modeling results with aerial imagery as a trusted externaljudge- ment. In all three experiments, no false negative is found, i.e., all building structures are successfully detected and reconstructed by the modeling system. In addition, false Dataset Point# Resolution Building # Building Tri. # Building errorrate Octahedron # Atlanta#1 5.5M 22.9/m 2 418 55,568 1.1% 52,924 Atlanta#2 4.0M 18.8/m 2 323 61,492 0.7% 29,151 Denver 1.0M 6.3/m 2 290 42,942 0.6% 17,054 Table8.1: Statisticsoftheexperimentsonthreedifferentdatasets 125 Dataset Classification Normal estimation Building modeling Tree modeling Totaltime Atlanta#1 9s 45s 54s <1s 109s Atlanta#2 6s 29s 32s <1s 68s Denver 3s 8s 11s <1s 23s Table8.2: Executiontimeofeachstepinthemodelingsystem positives exist in the form of small building-like trees and incorrectly classified ground components. The error rate is then calculated as the ratio of the number of triangles in false positives to the total number of triangles in all building models. Table 8.1 con- tainsstatisticsofthethreeexperiments,inwhicherrorratesaregenerallylow. Table8.2 shows the computation time, measured on a laptop with Intel i-7 CPU 1.60GHz and 6GBmemory. 126 (b) (b) (a) (a) (c) (c) Input point cloud Classification results Urban reconstruction Aerial Imagery Atlanta area #1 Figure8.6: Urbanmodelsreconstructedfrom5.5MaerialLiDARpointsforaresidential areainthecityofAtlanta 127 Chapter9 ConclusionandFutureWork 9.1 Conclusion Thisresearchstudiesthecomplicatedproblemofreconstructing3Durbanmodelsfrom aerial LiDAR point clouds. An automatic urban modeling system is proposed in Chap- ter 3, which divides this problem into four sub-problems, namely, classification, seg- mentation,buildingmodeling,andterrainmodeling. Automaticalgorithmsaregivento solve these problems individually. Two extensions based on prior knowledge are intro- ducedtoenhancethesystem,includingaprincipaldirectionsnappingmechanismanda non-planarshapemodelingmodule. In order to give the urban modeling system the ability to seamlessly handle huge data sets with billions of LiDAR points, Chapter 4 proposes a streaming framework which utilizes an out-of-core processing architecture. By decomposing each pipeline module into streaming operators and streaming states, the general urban modeling sys- temisadaptedintothestreamingframework. Experimentsaredoneonafewlarge-scale data sets, which no previous approach is able to process with such a small amount of resource. FromChapter5toChapter8,Istudythe2.5Dnatureofbuildingstructures. Theories andalgorithmsareproposedbasedonthe2.5Dcharacteristicofbuildingmodels. InChapter5,a2.5Dgeometryrepresentationisgivenforpolygonalbuildingmodels. A robust data-driven method is proposed to automatically create building models from aerial LiDAR point clouds. The reconstruction results are 2.5D models composed of 128 complex building roofs connected by vertical walls. By extending dual contouring into a2.5Dmethod,the2.5Ddualcontouringalgorithmoptimizesthesurfacegeometryand the boundaries of roof layers simultaneously. Sharp features are detected and faithfully preservedduringtriangulation. Chapter 6 defines 2.5D building topology as a combination of topological features and the associations between them. Convenient tools are given to change model geom- etry without modifying the topology. In addition, 2.5D dual contouring is extended with topology control strategies, to achieve a more flexible adaptive structure for sim- plification. The outputs have the same representability as models created by 2.5D dual contouring,butcontainfewerverticesandtriangles. Chapter 7, defines global regularities in 2.5D building models to characterize the intrinsic structure of building models. A primitive-based algorithm is proposed to dis- cover and enforce global regularities through a series of alignment steps. This building modeling method automatically integrates global regularities and local geometry, and thus creates 2.5D building models with high quality in terms of both geometry and visualjudgement. FinallyinChapter8,theurbanmodelingproblemisextendedfromdowntownareas to residential urban areas with rich vegetation. I observe the key difference between buildings and trees in terms of the 2.5D characteristic: while buildings are composed of opaque skyward rooftops and vertical walls, trees allow point samples underneath the crown. This feature enables a powerful classification algorithm based on an energy minimization scheme. By combing classification, building modeling and tree model- ing together, the residential urban modeling system automatically reconstructs a hybrid model composed of buildings and trees from the aerial LiDAR scan. The experiments demonstratetheeffectivenessandefficiencyofthissystem. 129 9.2 FutureWork Possiblefutureworklieswithinthreedirections. ModelingfromdenseLiDARdata With the fast advance of acquisition techniques, extremely dense LiDAR point clouds become available. E.g., a scan for the city of Vancouver in the year 2011 has approximately70samples/m 2 resolution,comparingwith6∼ 17samples/m 2 resolution in the experiments shown in this thesis. These extremely dense data sets bring two new challenges. First, besides typical urban features (buildings, trees, and ground), small urban features such as vehicles, poles, and signal lights are also accessible in the dense LiDAR data sets. It is a desire to reconstruct these features towards a better urban real- ity. Second, detailed geometry of building facades are partially captured by the dense LiDARscans. Buildingmodelingmethodcanbenefitfromthisadditionalinformation. Buildingmodelingwithbothrooftopsandfacades This thesis adopts the 2.5D characteristic of building structures, and describes the building facades as plain walls orthogonal to the x-y plane. However, in some applica- tion, building models are required to have both rooftops and detailed facade structures. Whilethisthesisfocusesonthebuildingroofreconstruction,manyresearcheffortshave attempted to model building facades from imagery (e.g., [42]) or ground based LiDAR (e.g. [43, 67]). Future research may seek a conjunction between roof modeling and facademodeling. Thehybridapproachshouldconnectthesetwotypesofbuildingstruc- tureswhilemaintainingtheconsistencyatroofboundaries. Inaddition,the2.5Dnature ofbuildingmodelsmayserveasasupportingstructureforfacademodelingapproaches. The global regularities in 2.5D building models can be further extended to handle the facadestructures. 130 IntegratingaerialLiDARwithotherdatasources Aerial LiDAR data has great advantages in the urban reconstruction problem, as demonstratedinthisthesis. However,itslimitationisalsonoticed: aspointsamplesare collected from nadir perspective, aerial LiDAR data hardly captures surfaces beneath opaque objects (e.g., bridges, highways, and vehicles). Future research may alleviate thisproblembyintroducingotherdatasourcestotheurbanmodelingsystem,including aerialimagery,groundbasedLiDAR,andpanoramaimages(streetview). 131 Bibliography [1] A. Alharthy and J. Bethel. Heuristic filtering and 3d feature extraction from lidar data. In ISPRS Commission III, Symposium,2002. [2] B. E. Boser, I. Guyon, and V. Vapnik. A training algorithm for optimal margin classifiers. Computational Learing Theory,pages144–152,1992. [3] Y. Boykov, O. Veksler, and R. Zabih. Fast approximate energy minimization via graphcuts. IEEE PAMI,2001. [4] M. Carlberg, P. Gao, G. Chen, and A. Zakhor. Classifying urban landscape in aeriallidarusing3dshapeanalysis. IEEE ICIP,2009. [5] G. Chen and A. Zakhor. 2d tree detection in large urban landscapes using aerial lidardata. IEEE ICIP,2009. [6] D.Cohen-Steiner,P.Alliez,andM.Desbrun. Variationalshapeapproximation. In ACM SIGGRAPH,2004. [7] T. H. Cormen, C. E. Leiserson, R. L. Rivest, and C. Stein. Introduction to Algo- rithms. MITPressandMcGraw-Hill,secondedition,2001. [8] C.CortesandV.Vapnik. Support-vectornetworks. Machine Learning,1995. [9] J.-F. Cˆ ot´ e, J.-L. Widlowski, R. A. Fournier, and M. M. Verstraete. The structural andradiativeconsistencyofthree-dimensionaltreereconstructionsfromterrestrial lidar. Remote Sensing of Environment,2009. [10] B.CurlessandM.Levoy. Avolumetricmethodforbuildingcomplexmodelsfrom rangeimages. In ACM SIGGRAPH,1996. [11] A. Elaksher and J. Bethel. Reconstructing 3d buildings from lidar data. In ISPRS Commission III, Symposium,2002. [12] B.S.Everitt,S.Landau,andM.Leese. ClusterAnalysis. Wiley,4thedition,2009. [13] M. Fiocco, G. Bostr¨ om, J. G. M. Gonc ¸alves, and V. Sequeira. Multisensor fusion forvolumetricreconstructionoflargeoutdoorareas. 3DIM,2005. [14] M.A.FischlerandR.C.Bolles. Randomsampleconsensus: aparadigmformodel fitting with applications to image analysis and automated cartography. Communi- cations of the ACM,1981. 132 [15] R.A.Fisher.Thestatisticalutilizationofmultiplemeasurements.AnnalsofEugen- ics,pages376–386,1938. [16] R.B.Fisher. Applyingknowledgetoreverseengineeringproblems. CAD,2004. [17] D.A.ForsythandJ.Ponce. ComputerVision: AModernApproach. PrenticeHall, 1999. [18] C.Fr¨ uh,S.Jain,andA.Zakhor.Dataprocessingalgorithmsforgeneratingtextured 3dbuildingfacademeshesfromlaserscansandcameraimages. IJCV,2005. [19] R. Gal, O. Sorkine, N. J. Mitra, and D. Cohen-Or. iwires: An analyze-and-edit approachtoshapemanipulation. In ACM SIGGRAPH,2009. [20] M. Garland. Quadric-based polygonal surface simplification. PhD thesis, Com- puter Science Department, Carnegie Mellon University, CMU-CS-99-105,1999. [21] Google. Google3dwarehouse. http://sketchup.google.com/3dwarehouse/. [22] T. L. Haithcoat, W. Song, and J. D. Hipple. Building footprint extraction and 3-d reconstructionfromlidardata. InIEEE/ISPRSJointWorkshoponRemoteSensing and Data Fusion over Urban Areas,2001. [23] A.Hatcher. Algebraic Topology. CambridgeUniversityPress,firstedition,2001. [24] J.Hu,S.You,andU.Neumann. Integratinglidar,aerialimageandgroundimages forcompleteurbanbuildingmodeling. In 3DPVT,2006. [25] M. Isenburg and S. Gumhold. Out-of-core compression for gigantic polygon meshes. In ACM SIGGRAPH,2003. [26] M.IsenburgandP.Lindstrom. Streamingmeshes. In IEEE Visualization,2005. [27] M. Isenburg, Y. Liu, J. Shewchuk, and J. Snoeyink. Streaming computation of delaunaytriangulations. In ACM SIGGRAPH,2006. [28] M. Isenburg, Y. Liu, J. Shewchuk, J. Snoeyink, and T. Thirion. Generating raster dem from mass points via tin streaming. In Proceedings of Geographic Informa- tion Science,2006. [29] T.Joachims. Svmlight. http://svmlight.joachims.org/,2004. [30] I.T.Jolliffe. Principal Component Analysis. Springer,secondedition,2003. [31] T.Ju, F.Losasso, S. Schaefer, andJ. Warren. Dualcontouring onhermite data. In ACM SIGGRAPH,2002. [32] Y.Kurozumiand W.Davis. Polygonal approximationbythe minimaxmethod. In Computer Graphics and Image Processing,1982. [33] F.Lafarge,X.Descombes,J.Zerubia,andM.Pierrot-Deseilligny. Buildingrecon- structionfromasingledem. In CVPR,2008. [34] F. Lafarge and C. Mallet. Building large urban environments from unstructured pointdata. In ICCV,2011. 133 [35] J.-F. Lalonde, N. Vandapel, and M. Hebert. Data structure for efficient dynamic processingin3-d. IJRR,2007. [36] Y. Li, X. Wu, Y. Chrysanthou, A. Sharf, D. Cohen-Or, and N. J. Mitra. Glob- fit: Consistently fitting primitives by discovering global relations. In ACM SIG- GRAPH,2011. [37] P.Lindstrom. Out-of-coresimplificationoflargepolygonalmodels. InACMSIG- GRAPH,2000. [38] Y. Livny, S. Pirk, Z. Cheng, F. Yan, O. Deussen, D. Cohen-Or, and B. Chen. Texture-lobesfortreemodelling. In ACM SIGGRAPH,2011. [39] S.K.Lodha,D.M.Fitzpatrick,andD.P.Helmbold. Aeriallidardataclassification usingadaboost. In 3DIM,2007. [40] W. Lorensen and H. Cline. Marching cubes: A high resolution 3d surface con- structionalgorithm. In ACM SIGGRAPH,1987. [41] B.Matei,H.Sawhney,S.Samarasekera,J.Kim,andR.Kumar. Buildingsegmen- tationfordenselybuilturbanregionsusingaeriallidardata. In CVPR,2008. [42] P.M¨ uller,G.Zeng,P.Wonka,andL.VanGool. Image-basedproceduralmodeling offacades. In ACM SIGGRAPH,2007. [43] L. Nan, A. Sharf, HaoZhang, D. Cohen-Or, and B. Chen. Smartboxes for interac- tiveurbanreconstruction. In ACM SIGGRAPH,2010. [44] B.Neubert,T.Franken,andO.Deussen. Approximateimage-basedtree-modeling usingparticleflows. In ACM SIGGRAPH,2007. [45] M. Nielsen, O. Nilsson, A. Soderstrom, and K. Museth. Out-of-core and com- pressedlevelsetmethods. In ACM SIGGRAPH,2007. [46] R.Pajarola. Stream-processingpoints. In IEEE Visualization,2005. [47] M.Pauly. Pointprimitivesforinteractivemodelingandprocessingof3dgeometry. PhD thesis, ETH Zurich,2003. [48] C.PoullisandS.You. Automaticreconstructionofcitiesfromremotesensordata. In CVPR,2009. [49] G.Priestnall,J.Jaafar,andA.Duncan. Extractingurbanfeaturesfromlidardigital surfacemodels. Computers, Environment and Urban Systems,2000. [50] F.Rottensteiner. Automaticgenerationofhigh-qualitybuildingmodelsfromlidar data. IEEE Computer Graphics and Applications,2003. [51] S. Schaefer, T. Ju, and J. Warren. Manifold dual contouring. IEEE Trans. Vis. Comput. Graph.,2007. [52] J. Secord and A. Zakhor. Tree detection in urban regions using aerial lidar and imagedata. IEEE Geoscience and Remote Sensing Letters,2007. 134 [53] X. Shi, H. Bao, and K. Zhou. Out-of-core multigrid solver for streaming meshes. In ACM SIGGRAPH Asia,2009. [54] P. Soille. Morphological Image Analysis: Principles and Applications. Springer- VerlagTelos,2002. [55] P. Tan, T. Fang, J. Xiao, P. Zhao, and L. Quan. Single image tree modeling. In ACM SIGGRAPH Asia,2008. [56] P. Tan, G. Zeng, J. Wang, S. B. Kang, and L. Quan. Image-based tree modeling. In ACM SIGGRAPH,2007. [57] S.ThrunandB.Wegbreit. Shapefromsymmetry. In ICCV,2005. [58] A. Toshev, P. Mordohai, and B. Taskar. Detecting and parsing architecture at city scalefromrangedata. In CVPR,2010. [59] G.Varadhan,S.Krishnan,Y.Kim,andD.Manocha. Feature-sensitivesubdivision andiso-surfacereconstruction. In IEEE Visualization,2003. [60] V. Verma, R. Kumar, and S. Hsu. 3d building detection and modeling from aerial lidardata. In CVPR,2006. [61] H.Vo,S.Callahan,P.Lindstrom,V.Pascucci,andC.Silva. Streamingsimplifica- tionoftetrahedralmeshes. IEEE Trans. Vis. Comput. Graph.,2007. [62] O. Wang, S. K. Lodha, and D. P. Helmbold. A bayesian approach to building footprintextractionfromaeriallidardata. In 3DPVT,2006. [63] H. Xu, N. Gossett, and B. Chen. Knowledge and heuristic-based modeling of laser-scannedtrees. ACM Trans. Graph.,2007. [64] S. You, J. Hu, U. Neumann, and P. Fox. Urban site modeling from lidar. In Proceedings, Part III,ICCSA,2003. [65] L.Zebedin,J.Bauer,K.Karner,andH.Bischof. Fusionoffeature-andarea-based informationforurbanbuildingsmodelingfromaerialimagery. In ECCV,2008. [66] N. Zhang, W. Hong, and A. Kaufman. Dual contouring with topology preserving simplificationusingenhancedcellrepresentation. In IEEE Visualization,2004. [67] Q. Zheng, A. Sharf, G. Wan, Y. Li, N. J. Mitra, B. Chen, and D. Cohen-Or. Non- localscanconsolidationfor3durbanscene. In ACM SIGGRAPH,2010. [68] Q.-Y.Zhou,T.Ju,andS.-M.Hu. Topologyrepairofsolidmodelsusingskeletons. IEEE Trans. Vis. Comput. Graph.,2007. 135

Abstract (if available)

Abstract 3D reconstruction from point clouds is a fundamental problem in both computer vision and computer graphics. As urban modeling is an important reconstruction problem that has various significant applications, this thesis investigates the complex problem of reconstructing 3D urban models from aerial LiDAR (Light Detection And Ranging) point cloud. ❧ In the first part of this thesis, an automatic urban modeling system is proposed which consists of three modules: (1) the classification module classifies input points into trees and buildings

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

City-scale aerial LiDAR point cloud visualization

PDF

Object detection and recognition from 3D point clouds

PDF

Rapid creation of photorealistic large-scale urban city models

PDF

Line segment matching and its applications in 3D urban modeling

PDF

Point-based representations for 3D perception and reconstruction

PDF

Accurate 3D model acquisition from imagery data

PDF

3D face surface and texture synthesis from 2D landmarks of a single face sketch

PDF

3D deep learning for perception and modeling

PDF

Efficient SLAM for scanning LiDAR sensors using combined plane and point features

PDF

Interactive rapid part-based 3d modeling from a single image and its applications

PDF

3D object detection in industrial site point clouds

PDF

Deep representations for shapes, structures and motion

PDF

Hybrid methods for robust image matching and its application in augmented reality

PDF

Automatic image matching for mobile multimedia applications

PDF

Towards building a live 3D digital twin of the world

PDF

Face recognition and 3D face modeling from images in the wild

PDF

Data-driven 3D hair digitization

PDF

Plant substructuring and real-time simulation using model reduction

PDF

Integrating complementary information for photorealistic representation of large-scale environments

PDF

CUDA deformers for model reduction

Asset Metadata

Creator Zhou, Qian-Yi (author)

Core Title 3D urban modeling from city-scale aerial LiDAR data

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Computer Science

Publication Date 07/26/2012

Defense Date 04/03/2012

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag 2.5D,global regularity,LiDAR,OAI-PMH Harvest,streaming,tree detection,urban modeling

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Neumann, Ulrich (committee chair), Barbic, Jerney (committee member), Kuo, C.-C. Jay (committee member), You, Suya (committee member)

Creator Email Qianyi.Zhou@gmail.com,qianyizh@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-70401

Unique identifier UC11289014

Identifier usctheses-c3-70401 (legacy record id)

Legacy Identifier etd-ZhouQianYi-1030.pdf

Dmrecord 70401

Document Type Dissertation

Rights Zhou, Qian-Yi

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA