Building Reconstruction of Digital Height Models with the Markov Chain Monte Carlo Method

(1)

Master of Science Thesis in Electrical Engineering

Department of Electrical Engineering, Linköping University, 2018

Building Reconstruction of

Digital Height Models with

the Markov Chain Monte

Carlo Method

(2)

Carlo Method Mats Nilsson LiTH-ISY-EX--18/5130--SE Supervisor: Karl Holmquist

isy, Linköpings universitet

Bertil Grelsson

Saab Dynamics

Examiner: Maria Magnusson

isy, Linköpings universitet

Computer Vision Laboratory Department of Electrical Engineering

(3)

Abstract

Data about the earth is increasing in value and demand from customers, but it is difficult to produce accurately and cheap. This thesis examines if it is possible to take low resolution and distorted3ddata and increase the accuracy of

build-ing geometry by performbuild-ing buildbuild-ing reconstruction. Buildbuild-ing reconstruction is performed with a Markov chain Monte Carlo method where building primitives are placed iteratively until a good fit is found. The digital height model and pixel classification used is produced by Vricon. The method is able to correctly place primitive models, but often overestimate their dimensions by about 15%.

(4)

(5)

Acknowledgments

I would like to thank my supervisors at Saab Dynamics and Linköping University, Bertil and Karl for invaluable discussions as well as their advice and knowledge on the subject. I would also like to thank my examiner Maria for the help and feedback on the thesis.

Linköping, June 2018 Mats Nilsson

(6)

(7)

Notation

Sets and mathematical expressions Notation Interpretation

N+ The set of positive integers, {1, 2, 3, . . . }.

A◦ B Hadamard product of matrices A and B of same size.

Abbreviations

Abbreviation Interpretation dsm Digital Surface Model dtm Digital Terrain Model

dhm Digital Height Model

lod Level of Detail

mc Markov chain

mcmc Markov chain Monte Carlo

sa Simulated Annealing

mbr Minimum Bounding Rectangle

lir Largest Inscribed Rectangle

(10)

(11)

1

Introduction

Accurate data of the earth is becoming ever more valuable for customers, but it is often difficult to produce. Building reconstruction is the process of taking low resolution and distorted height data of the earth and reconstruct the real world buildings as close as possible.

1.1 Motivation

Digital surface models (dsm) is one way to digitize the height of objects on the earth’s surface. In figures 1.1b–1.1c a dsm of a part of Linköping can be seen and in figure 1.1a an orthophoto, which is similar to a satellite photo, of the same area is shown. A more thorough description of the dataset can be found in section 2.1 and section 2.2. These3dmaps have a large potential in very diverse military and

consumer applications. The military can do intelligence gathering, analysis and mission planning in remote and hostile areas. Consumer applications include urban and telecommunications planning, pollution propagation simulations and aerial navigation. [8] [11]

The dsm records the height of points on the earth’s surface. The measurement can be done with several methods, light detection and ranging (lidar), structure-from-motion or stereo photography all being able to provide the depth informa-tion. Stereo photography or lidar from an aeroplane on low altitude have very high resolution but have the disadvantage that it takes a lot of time to fly over large areas and that the areas of interest can be dangerous to fly in. Satellites have lower resolution but can provide data of larger areas faster than by air planes or drones as well as not having to fly on low altitude over high risk areas. [19] The company providing the dsm for this thesis is Vricon, a company owned by

(12)

(a)Orthophoto. (b) dsm. (c) dsmside view.

Figure 1.1:Orthophoto and digital surface model of the city of Linköping.

the Swedish defence contractor Saab and DigitalGlobe. DigitalGlobe owns the largest private satellite fleet that is gathering the photos and Vricon owns the image processing algorithms creating dsms from the stereo photos. [1] [2] Creating dsms from satellite imagery results in lower resolution with more noise and distortions compared to aerial lidar or stereo photography. Certain applica-tions require high precision between the3dmodel and the real world which dsms

created from satellites currently lack. For example it is currently not possible to determine building dimensions exactly since the buildings appear smudged, see figure 1.1c.

Building reconstruction utilise the data provided by the dsm and prior informa-tion about how most buildings look. For example most buildings are built with straight walls that is perpendicular to the ground. The dataset contains buildings in Linköping which usually have either a flat, gable or hipped roof, see figure 3.2. The combined data and prior information can be used to create an automatic building reconstruction algorithm which tries to model how the building look in the real world.

1.2 Purpose

The purpose of this thesis is to measure the geometric accuracy improvement of an automatic building reconstruction algorithm on the dsm by utilizing the dsm, orthophoto and pixel terrain classification provided by Vricon.

The objective is to achieve building models of lod2, see section 2.2.3 and

fig-ure 2.3c. The output model should have correct footprint and the correct roof type. After processing an area, the resulting models should be saved to allow for reuse of the refined buildings. Human interaction with the system should be kept to a minimum. The goal is complete automation of the building reconstruction process.

(13)

1.3 Problem Formulation 3

1.3 Problem Formulation

The motivation and purpose of the thesis give rise to two principal questions to answer:

• What measurement accuracy improvement is achieved after the building reconstruction process?

• Is the proposed algorithm computationally effective enough to be used for practical applications?

1.4 Limitations

The algorithm is only required to reconstruct buildings up to lod2, see figure 2.3,

and handle three roof types. The three building roof types is flat, gable and hipped, see figure 3.2. Buildings containing more than one primitive should correctly place sufficiently many primitives, but the primitives do not need to be merged or optimised globally. Terraced buildings is not an algorithm require-ment.

(14)

(15)

2

Theory

This chapter presents background information and theory behind the chosen method. Section 2.1 and section 2.2 present a brief description of the data in this thesis. Following is an introduction to Markov chains (mc) and Markov chain Monte Carlo methods (mcmc) is given in section 2.3 and section 2.4. The sec-tions also include a motivation to why the proposed method should be able to solve the problem. The chapter concludes with a description of simulated anneal-ing in section 2.5. This explains how the transition between different states in the optimisation procedure is decided.

2.1 Orthophoto

An orthophoto is a satellite or aerial photo that lacks perspective view, therefore looking just like a map but with the advantage of visual features like a photo, an example of an orthophoto is seen in figure 1.1a. All parts of an orthophoto looks like it has been taken from above that location, making all projection rays parallel to each other and orthogonal to projection plane, as seen in figure 2.1b. Orthophotos containing buildings will only show the roof of the buildings in the photo and not the walls, since they are normally perpendicular to the datum plane. [3]

An orthophoto is impossible to achieve from only one photo. One photo will have only one point where it is captured with no perspective view, perpendicular from the sensor, see figure 2.1a. An orthophoto is instead synthetically created from several photos and the dsm stitched together. This makes it possible to correlate orthorectified satellite photos with dsms making processing much more effective. [3]

(16)

(a)Normal photo. (b)Orthophoto.

Figure 2.1:Different perspectives in photos.

2.2 Digital Models

A model is a simplified representation of something else. Here follows how the earth and its surface are modelled and how building models are categorised.

2.2.1 Earth Models

To model the Earth one could assume a perfect sphere with some radius r. This can suffice for many applications, but some require higher accuracy since the Earth is not a perfect sphere. The dsm used in this thesis is related to the WGS84 ellipsoid which is defined by the U.S. Department of Defence. [14]

2.2.2 Surface Models

With photos taken from different viewpoints a depth image can be created. This is why predators in the animal kingdom often have two forward-facing eyes, to accurately estimate distances to preys. If many photos are stitched together a depth surface can be created as in figure 2.2a, this is called a digital surface model, dsm.

Depending on the application different digital models can be used. In literature there is no clear consensus on the naming of different types of digital models, the notation used in this thesis is from Li et al. [12]. Sometimes it is more interesting just to look at the the terrain of the Earth without trees and buildings. This model is called digital terrain model, dtm, see figure 2.2b. It is important to note that the terrain model is an approximation with no definite answer since it is impossible to know the true value of the ground level underneath buildings and trees which cover the ground. The height difference between the dsm and the dtm is the digital height model, dhm, see figure 2.2c. The dhm should in theory only contain the values of everything above the defined ground level in

(17)

2.3 Markov Chain 7

the dtm. The dhm is the better choice of surface model instead of the dsm to use in this applicaiton since the height is normalised to the ground level. This makes the height of the primitive equal to the height of the building in the real world, which makes processing easier and more precise. Another reason for choosing the dhm is that the contribution of pixel values outside the actual building and near the ground should be negligible. The difference in height between the top of the building and the ground should therefore be as high as possible, which is the result when using the dhm.

Height [m] 90 100 110

(a)Digital surface model.

Height [m] 88 90 92

(b)Digital terrain model.

Height [m]

0 10 20

(c)Digital height model.

Figure 2.2:Differences between three digital models: dsm, dtm and dhm.

2.2.3 Building Models and their Level of Detail

Reconstructed buildings can be categorized according to their level of detail (lod), see figure 2.3 for examples of different lods ranging from just the building foot-print (lod0) to also including details on the roof (lod3). Generally speaking the

output lod of the building reconstruction is dependent on the resolution and quality of the input data. A more thorough treatment of the building categorisa-tion subject can be found in Biljecki et al. [4] and Kolbe et al. [10].

2.3 Markov Chain

A stationary discrete Markov chain (mc) is a stochastic process, Xt, describing a

sequence of events in discrete time where the probability of a certain next event only depends on the current state. This description will only treat stationary Markov chains, which means that the probabilities of the next state is not de-pendent on the time, only on the current state, i.e. the Markov chain is time-homogeneous. This memoryless property is called the Markov property and can be applied to other stochastic processes as well. A formal mathematical

(18)

descrip-(a) lod0 (b) lod1 (c) lod2 (d) lod3 Figure 2.3:Building level of detail specification. tion of the Markov property is therefore defined as

Pr (Xt+1= xt+1| Xt= xt, ..., X0 = x0) = Pr (Xt+1= xt+1| Xt= xt) . (2.1)

An example would be the stochastic process X which has two states i and j. The probability of moving from state i to the next state j is Pr (j | i) = pi,j. All

tran-sitions between states can be summarised in a transition matrix T . The element p in row i and column j of the transition matrix is the probability of moving from state i to j. The sum of every row must therefore be equal to one, since the probability of jumping to the next state must happen with probability one. An example is the stochastic process, X, which has two states, 1, 2 and a transition matrix T , see figure 2.4. [5] [17]

1 2 T = p1,1 p1,2 p_2,1 p_2,2 ! p_1,2 p_2,1 p1,1 p2,2

Figure 2.4:Example of a Markov chain with two states.

Another example would be to model the weather as a Markov chain where the weather for the next day is only dependent on the current weather conditions. For simplicity the weather in this example has only three possible states, sunny, rainy and cloudy, see figure 2.5 for a visualization. The number adjacent to the directed edges is the probability for that state transition to happen.

2.4 Markov Chain Monte Carlo

Markov chains can be used to simulate complex probability distributions with the mcmcmethod, usually where the distributions have high dimensionality. mcmc is a way of simulating these complex probability distributions. The probability distributions are denoted π = (π1, π2, . . . , πd) where πk is one parameter of the

(19)

2.4 Markov Chain Monte Carlo 9 Sunny Cloudy Rainy 0.3 0.3 0.1 0.2 0.5 0.4 0.5 0.4 0.3

Figure 2.5:Example of the weather modelled as an Markov chain.

chain where the distribution in the limit is π. A sample from π is then the output from the Markov chain if run long enough so that it (almost) converges. One popular algorithm for simulating and then sampling from these distributions is the Metropolis-Hastings algorithm, see algorithm 1. From the current state, xt

a proposed state x0 _{is generated. The candidate, x}0_{, is usually sampled from a} symmetric distribution q, which makes the fraction q (x | x0_{) /q (x}0_{| x) = 1. An} example of such a q (x0_{| x) is the gaussian distribution centered at x. [5] [17]}

Algorithm 1Metropolis-Hastings 1: Initialize x0and t ← 0 2: whilet < N do

3: Generate next state candidate x0_{from q (x}0_{| x}_t₎ 4: Acceptance ratio α = min1,_π(xπ(x0)

t)

q(x_{| x}0)

q(x0_{| x)}

5: Sample from the uniform distribution u ∼ U(0, 1) 6: ifu < α then

7: Accept candidate xt+1← x0

8: else

9: Reject candidate and keep current state xt+1← xt

10: Set t ← t + 1

In the scope of this thesis, an example would be to initialise a building at a proba-ble location with the parameters Θ. Then generate some candidate building with a slight change in parameters, Θ0_{, and see how well the model fits the data. If} Θ0 fits better, accept it as current Θ, otherwise accept it with some probability dependent on how well it fits the data.

(20)

2.5 Simulated Annealing

Simulated annealing (sa) is a method of approximating the global optimum of a function where an exhaustive search is infeasible. Metropolis et al. introduced the ideas in a 1953 paper [13] and Kirkpatrick et al. extended it in 1983 [9]. The main idea behind the technique comes from controlled and slow cooling of metals to reduce the defects and therefore increasing the strength of the metal. The probability, P , of a system being in a state, i, is a function of the energy of that state, εi, Boltzmann constant kB, and the temperature, T , of the system with

a total of M states as Pi = exp− εi k_BT P_M j=1exp − εj kBT  . (2.2)

A state with high energy is less probable than a state with low energy. By gradu-ally lowering the total energy of the system, the system will naturgradu-ally converge to a state requiring low energy, with few defects.

This method can be applied to an optimisation problem with many local extrema using the mcmc method described previously. The prerequisite is that there ex-ists a cost function, E(M, D), dependent on the model, M, and data, D, which is in analogy with the energy of that state. The cost function used in this thesis is defined in section 3.5. The probability of going from state i − 1 to i is defined as

Pi = min 1, exp −∆Ei T , (2.3) where ∆Ei = Ei − min {E1, . . . , Ei−1} , (2.4) and T being the temperature of the system.

The starting temperature, Ts, and final temperature, Tf, is calculated from

T = −1

ln(p), (2.5)

where p is the probability of choosing a worse candidate in either the starting state or the final state. The constant annealing factor, τ, is then calculated as

τ = Tf Ts

!1

N

, (2.6)

where N is the total number of iterations. The temperature of the next state is the annealing factor times the current temperature as

Ti+1= τ · Ti. (2.7)

If the temperature T is high, the probability of accepting the proposed state and making the jump is greater compared to when the temperature is low. If instead the proposed state has a lower cost E than seen before, therefore being a more

(21)

2.5 Simulated Annealing 11

likely candidate for optimum, the jump to that state is made with probability 1. Otherwise depending on E the proposed state will be accepted with some probability P according to equation (2.3). By starting with a high temperature and gradually lowering it, the function is exploring many possible states early. Then it can settle on the best state where the cost function is the lowest. [15]

(22)

(23)

3

Method

This chapter will describe the method used to solve the problem posed earlier.

3.1 Method structure

The method is based on the works of Huang et al. [7] and Partovi et al. [16]. The algorithm structure can be broken down as follows:

1. Preprocessing of terrain classification data. 2. Get probable initialisation from dhm. 3. Iterate primitive placement with mcmc.

3.2 Data Preprocessing

The data provided by Vricon is the following: • Digital Surface Model (dsm)

• Digital Terrain Model (dtm) • Digital Height Model (dhm) • Orthophoto

• Terrain classification image

The data covers a residential part of the city of Linköping, Sweden, and has a reso-lution of approximately 4 pixels/m2_{. The terrain classification image is an image} of the same size and resolution as the dsm and the orthophoto. Each pixel in the

(24)

classification image contains information of what the pixel has been classified as, for example “building” and “road”, see figure 3.1a.

Due to some noise in the classification image for buildings, figure 3.1b, it was pro-cessed by first filling holes, using the MATLAB™ command imfill. To remove small objects, the morphological operation opening, imopen, with a structuring element being a disk of radius 8, strel(’disk’,8), was performed. Buildings with an area less than 125 m2_{(500 pixels) were removed with bwareaopen. The} reason was that buildings smaller than 125 m2 _{were deemed too hard to} recon-struct at current resolution and noise level. The resulting image can be seen in figure 3.1c. Each connected component should be a separate building, which is then processed individualy. To allow the model to extend slightly outside the classification, the data is padded with zeros and is henceforth called D.

(a)Classification. (b)Buildings. (c)Building noise removed.

Figure 3.1:Removing noise from building labels in a classification image.

3.3 Primitive Library

A primitive library is constructed similar to Huang et al. [7] and Partovi et al. [16] with rectangular footprints. The building primitives are parametrized with

Θ= {x, y, α, w, l, z1, z2, h1, h2, h3, h4}, (3.1) and are related to the primitives as shown in figure 3.3. Depending on the build-ing complexity, a different amount of parameters is needed. The flat roof in fig-ure 3.2a needs only 6 parameters and the gable and hip roof in figfig-ures 3.2b–3.2c uses a maximum of 11 parameters.

For flat roofs z2 = 0, and the size of hips h1,...,4 do not matter so they are set to zero. The method assumes that a z2> 0 indicate a gable or hipped building. This is not correct in all cases as there exist a lot more different roof types. The method also assumes that buildings are symmetric. A gable and hip roof are assumed to have the ridge line in the middle, and the size of the smaller hips on each size is equal. This is described mathematically as if z2> 0 and w < l, then h3= h4= w/2

(25)

3.4 Initialisation 15

and h1 = h2. Otherwise, if z2 > 0 and l < w, then h1 = h2= l/2 and h3 = h4, see figure 3.3a.

(a)Flat. (b)Gable. (c)Hip.

Figure 3.2:Roof primitive library.

y x (x, y) w l h1 h2 h3 h4 α

(a)Side parameters.

z1 z2

(b)Height parameters.

Figure 3.3:Building parameters.

3.4 Initialisation

The algorithm requires initialisation parameters before the optimisation proce-dure. The initialisation should preferably be along the ridge of the building and as close to the real length, width and height dimensions as possible. The building ridge line is assumed to be the local maxima of the building classification blob. This is located with

h(r, c) = r+η X i=r−η c+η X j=c−η        1, if dhm(i, j) − ∆z ≤ dhm(r, c), 0, otherwise. (3.2)

Consequently h(r, c) is the number of pixels that the height specified in dhm(r, c), figure 3.4a, is higher than in some neighbouring area, (2η + 1)2_{, η} _{∈ N}+_, sub-tracted some ∆z. The local maxima function, h, produces an image as in fig-ure 3.4b. To achieve the image h, η = 10, ∆z = 0.8 and the threshold value was 400. This image is masked with the building classification image C (figure 3.1c)

(26)

to remove everything but the buildings, see figure 3.4c. It is then thresholded at a value dependent on the size of η. After thresholding the masked local maxima im-age, the skeleton function, bwmorph, is applied to retrieve a line. See figure 3.5a with the orthophoto as background to ease understanding. The Hough transform, houghpeaks, houghlines, is applied to each building’s local maxima to get the start and end points of the ridge line. For a description of functions used, see the MATLAB™ documentation or for example Serra 1983 [18].

The angle of the building is initialised with the angle α of the Hough line. The building length and width are calculated from the classification. The classifica-tion is rectified by rotating an angle −α around the center of the Hough line. The median number of non-zero pixels in each row and column gives the initialisa-tion of the length and width. When initializing a hipped roof, the size of the hip is set equal to the distance from the building edge to the ridge, h1= h2= h3= h4.

(a) dhm. (b)Local maxima image,h. (c) Local maxima masked with just buildings,h◦ C.

Figure 3.4:Obtaining local maxima from dhm.

3.5 Cost Function

For a given set of parameters Θ, a model M can be generated. The model is given in a coordinate system with the same size as the data, D, see section 3.3. Since the data is a bit larger than the actual building contained in it, and the model has the same size as D, the model M is not just the building, but also some zero padding around it. The cost function mentioned in section 2.5, E(M, D), is defined as the sum of the absolute difference in height for every pixel between the dhm and the model as E(M, D) =X r X c |M(r, c) − D(r, c)|, (r, c) ∈ Ω, (3.3) where Ω is the total area, represented by the number of data and model rows and columns, respectively. The area Ω will always be larger than the model since padding around the classification is always applied. Therefore if the dhm has a non-zero value z where the model is zero, the contribution for that pixel to the

(27)

3.6 Parameter Proposition 17

(a) After thresholding and using the skeleton function.

(b)Hough lines.

Figure 3.5:After thresholding the Hough lines can be calculated.

cost function will be exactly z. The goal is to minimize this cost function, which is equivalent with minimizing the difference between the model and data and the “system energy” as defined in section 2.5.

3.6 Parameter Proposition

The initialisation parameters are used to place a primitive. Using the cost func-tion, section 3.5, and Boltzmann probability, section 2.5, parameters are pro-posed and either accepted or rejected. The building parameters (see figure 3.3) are proposed with a truncated normal distribution, makedist, truncate. The normal distribution for the position parameters (x, y) are truncated so that it is en-sured that (x, y) are within the area D, mentioned in section 3.2, with a marginal of 2 pixels. The primitive cannot extend outside the edge of D since the data is already padded with zeros and is assumed to be big enough. The size of the normal distribution for the length, l, and width, w, is therefore dependent on the center position, (x, y), of the primitive and is calculated as

l ∈ [1, min {y, ymax− y}] , (3.4)

w∈ [1, min {x, xmax− x}] . (3.5)

The height to the roof, z1, is limited by reasonable values for the area tested, 1 to 30 metres,

z1∈ [1, 30] . (3.6)

(28)

and the algorithm tend to overestimate the height of the roof. The bound was chosen as z2∈ 0, min5,2 3z1 . (3.7)

3.7 Markov Chain Monte Carlo

The mcmc algorithm is run for N = 3000 iterations for each of the building roof types flat, gable and hipped. The temperature was calculated with starting probability ps = 0.9 and final probability pf = 0.2. To check if the placement

is okay, the average deviation is used. Average deviation takes into account only the pixels inside the model boundary and measure the distance to the dhm and divides by the area, ω, with side lengths l and w, see figure 3.3a, as

∆zavg. dev.(M, D) = P r P c|M(r, c) − D(r, c)| l · w , (r, c)∈ ω. (3.8)

If the average deviation of the currently tested roof type is ∆zavg. dev.> 3 m then the algorithm is applied again since it is assumed the algorithm did not con-verge correctly. This sometimes happens due to the randomness of the parameter proposition.

The question whether a roof is flat or not is determined by the cost function and average deviation. First a flat roof is initialised and optimised, then a gable roof is initialised and optimised. An increase in parameter space should mean a sig-nificant improvement in the fit of model, otherwise a simpler model is probably the best choice.

The two primitives are denoted Mflatand Mgable, the dhm is denoted D, and if E(Mgable, D)

E(Mflat, D) ·

∆zavg. dev.(Mgable, D)

∆zavg. dev.(Mflat, D) < 0.87, (3.9) the gable roof is significantly better than the flat roof. The corresponding formula is then used to decide between hipped roof and gable roof.

Figure 3.6 and figure A.3 show the evolution of the mcmc algorithm for three iterations, i ∈ {1, 1000, 3000}. The cost function decreasing as the temperature gets lowered is visualized in figure 3.7 and figure A.4.

3.8 Composite Buildings

After placing a primitive, the algorithm has to check whether there needs to be more primitives to accurately represent the entire building. This is done by check-ing the total cost reduction and intersection of model and data. The cost of the data with no model, E(D), is just the sum of heights for every pixel in the image. Since the classification in most cases is an overestimation of the building area there needs to be a stricter limit. The stricter intersection, I, is defined as the

(29)

3.8 Composite Buildings 19

(a)i = 1 (b)i = 1000 (c)i = 3000

(d)i = 1 (e)i = 1000 (f)i = 3000 Figure 3.6:Building 2 for iterationi∈ {1, 1000, 3000}

intersection of the model and data above the mean of the data as ¯

d = mean(D > 0), (3.10)

I(M, D) = M∩ (D > ¯d)

D > ¯d , (3.11)

see figure 3.8. A building is assumed to not require any more primitives if the cost function with primitive and data is 40% of the cost with just the data and covers 90% of the intersection limit. This is described by

E(M, D)

E(D) < 0.4 (3.12)

and

I(M, D) > 0.9, (3.13)

where M contains every primitive.

If this is not the case, the building is a composite building which requires more than one primitive. The classification and dhm is masked with the footprint of the first primitive and the process is repeated with detecting a new ridge line and iterative primitive placement.

(30)

100 ₁₀1 ₁₀2 ₁₀3 0 0.2 0.4 0.6 0.8 1 Iteration C ost of primitiv e [%]

Figure 3.7:Building 2 for iterationi∈ {1, . . . , 3000}

(a)Unprocessed classification. (b)Classification above mean height.

Figure 3.8:Difference between unprocessed and the processed classification.

3.9 Evaluation

The smudged and noisy edges of the dhm and classification create a problem defining a good measurement of the unprocessed building. To measure the edges of a rectangular building the boundaries must first be found. A proposed method is measuring the area of the minimum bounding rectangle (mbr) and largest in-scribed rectangle (lir) of the classification, see figure 3.9 and figure 3.10. The length and width of the rectangles are then checked against the dimensions of the proposed model from the mcmc algorithm and measured ground truth. By normalising every data point so the ground truth is equal to one the results can be compared. This evaluation method only works for buildings that can be viewed as approximately rectangular.

Ground truth data of building heights is not included in maps and measuring by hand is difficult and very time-consuming on large building, therefore some qualitative measure of the results is needed. When z1and z2heights are available it is easy to compare, otherwise a qualitative analysis must be performed. The

(31)

3.9 Evaluation 21

mbr lir Building

Figure 3.9: Example of how minimum bounding rectangle and largest in-scribed rectangle is defined from a building contour.

(a)Orthophoto. (b) Classification, mbrand lir.

(c)Orthophoto. (d) Classification, mbrand lir.

Figure 3.10: Orthophoto and an image with building classification (green), minimum bounding rectangle (light blue) and largest inscribed rectangle (yellow) overlaid for two buildings.

building roof is qualitatively evaluated by projecting all points of the dhm and primitive inside the classification to a plane. The normal of the plane is aligned with the ridge line of the proposed primitive, see figure 3.11. The roof width and type are also compared with ground truth data.

(32)

0 5 10 15 20 0 5 10 15 Width [m] Height [m] dhm Model

(33)

4

Results

The following chapter presents the results of the building reconstruction algo-rithm. The description contains quantitative measurements of some buildings and their measured ground truth for footprint, roof and angle. Some qualitative results of the building position are presented as well.

4.1 Ground Truth Data Collection

The building footprint dimensions, length and width, come from two distinct data sources. The bulk of the data comes from a measurement tool provided for free by the Linköping municipality through the web service “kartan”: http:// kartan.linkoping.se/spatialmap. Some buildings were measured in the real world to identify any scaling problems with the online tool. Buildings were measured by the provided online tool two times for each side and the average was rounded to nearest half metre.

Each building has been designated a unique number. In figure A.1 and figure A.2 the number and their digitally measured dimensions can be viewed. The dimen-sions can also be viewed in table A.1. Digital measurement were checked by measuring some buildings in the real world, the results can be seen in table A.2. The difference is at or below 0.3 m, with three outliers marked in bold and one missing value marked with “–”. The reasons for the outliers are probably mea-surement error or disagreement on where the building starts and ends. Building 3 have a significant length difference of almost 7 m. The building has garages on each side which may cause the large disagreement between the digital map and real world measurements. The available heights of the buildings are shown in table 4.1, where z1 and z2 are defined as in figure 3.3b. The roof height, z2,

(34)

are only available for gable roofs and not hipped roofs due to how the data was acquired.

4.2 Building Dimensions

The algorithm processed 28 buildings in total, 12 larger buildings, see figure A.1, and 16 smaller buildings, see figure A.2. They were chosen so that the unpro-cessed classification did not connect with nearby other buildings. The calculated length and width was normalized so that the ground truth is equal to one. For the unprocessed classification a lower bound and upper bound for the foot-print were calculated with lir and mbr respectively. Each of the 28 buildings have a length and a width, which makes a total of 2 · 28 = 56 data points for each of the three measurements, lower and upper bound as well as the mcmc algo-rithm. The data is visualised in figure 4.1 with a histogram and fitted normal distribution.

As can be seen in figure 4.1 the lower bound has a mean µlir= 1.01 whereas the mcmcmodel has a mean µmcmc = 1.15. This means that the proposed mcmc algorithm puts the dimensions about 15% too big. Both lir and mcmc have stan-dard deviations which is about the same, σlir = 0.15 and σmcmc = 0.14. The upper bound mbr performs significantly worse with a mean µmbr = 1.52, i.e. the dimensions are around 50% larger than ground truth and have higher uncer-tainty with a standard deviation of σmbr= 0.43.

4.3 Building Angle

Building angles were also collected from the same set of buildings used in sec-tion 4.2. Angles were measured digitally from the same map mensec-tioned previ-ously as well as from the mcmc algorithm, lir and mbr. The digitally measured angles was rounded to nearest integer in degrees, since high accuracy is impossi-ble to acheive by hand. The lir method also has integer degree resolution while the mcmc and mbr have higher resolution, so the last two was rounded to near-est tenth of an degree. The building angle set has been defined as α ∈ [0◦_{, 45}◦_{] so} that it is defined uniquely for these building primitives. The resulting values can be viewed in table A.3.

Figure 4.2 shows a histogram of the difference between the assumed true angle and predicted angle value for the three different methods. Largest inscribed rect-angle method has a normal distribution with a mean µlir= 0.07◦and a standard deviation of around σlir = 4.1◦. The proposed mcmc algorithm has a mean µmcmc = −0.18◦, which is slightly worse than for lir. However, mcmc has a significantly smaller standard deviation σmcmc = 2.1◦ compared to lir with σlir = 4.1◦. The mbr method performs much worse than the other two. It has a mean far from zero with µmbr= 3◦and a standard deviation of σmbr= 12◦, three times the lir method.

(35)

4.4 Building and Roof Height 25

4.4 Building and Roof Height

The available height data is more scarce than the footprint data and can be seen in table 4.1. There is one hipped roof, building 1, and 5 gable roofs, buildings 2–6. The 6 buildings were processed by the mcmc algorithm. Figure 4.3 shows the same phenomena as in figure 4.1, the mcmc algorithm usually overestimates the width of the building. The model roof height, z2, is also estimated to be larger than in the real world. There seems to be a large accuracy deviation between total estimated height and true height.

Table 4.1: Data measurements of building heights, wherez1andz2 are de-fined as in figure 3.3b andz = z1+ z2.

Real World [m] mcmc[m] Difference [m]

Building z1 z2 z z1 z2 z z1 z2 z 1 9.1 – – 9.2 5.0 14.2 0.1 – – 2 9.7 3.0 12.7 10.1 4.4 14.5 0.4 1.4 1.8 3 7.4 1.9 9.3 6.0 3.0 9.0 -1.4 1.1 -0.3 4 7.5 4.4 11.9 7.4 4.7 12.1 -0.1 0.3 0.2 5 7.8 4.0 11.8 8.3 5.0 13.3 0.5 1.0 1.5 6 7.7 3.8 11.5 7.2 4.8 12.0 -0.5 1.0 0.5

4.5 Building Classification

Building classification tries to determine what roof type the specified building has, flat, gable or hipped. Results from the processing show that the algorithm correctly classified 23 of 28 buildings, about 82%, as seen in table A.4. The large buildings have an accuracy of 10/12 ≈ 83% and the small buildings have an accuracy of 13/16 ≈ 81%. No clear distinction in accuracy between small and large building can be observed from the test data set.

4.6 Qualitative Measure of Building Position

A qualitative measurement of how well a certain model fits the building footprint can be viewed in figures 4.4–4.5. The models have correct roof type, gable in the first case and hipped in the second.

Figure 4.4 shows a larger gable building. The position of the building primitive seems correct in the dhm but slightly off in the orthophoto. By looking at the orthophoto in figure 4.4d and figure 4.5d, the model appears wider than ground truth. The reason is probably the smudging effect in the dhm.

In figure 4.5 a smaller hipped building is shown. According to the height con-tours in figure 4.4c the hip appears to be about the correct size. Position of the

(36)

primitive looks centered on the building though the length and width appear slighly too large. This supports the hypothesis from figure 4.1 where the nor-mal distribution indicated that the mcmc algorithm overestimated the building dimensions.

4.7 Composite Buildings

Building reconstruction of composite buildings do not have a quantitative evalu-ation. Instead only a qualitative evaluation is performed. The results of a com-posite building reconstruction is found in figure 4.6. All primitives have approx-imately the correct dimensions. The angle of the side primitives appear slightly off the correct one, whereas the main one looks right. Only the main primitive has the correct hip roof type. Since no global optimization is performed, the ridge lines do not connect and the sides of the building are discontinued. The general position and number of primitives is correct according to the orthophoto.

4.8 Terraced Buildings

The thesis do not include a requirement to successfully reconstruct terraced build-ings but to evaluate the feasibility of the proposed method, it should include a test of terraced buildings located in the Linköping suburbs, see figure 4.7. The algorithm correctly identifies the ridge lines of the six gable primitives in the classification area. Instead of propose multiple primitives, the reconstruction algorithm proposed one building primitive covering all buildings. In reality there should be 6 gable primitives and 5 flat primitives to correctly reconstruct the building.

(37)

4.8 Terraced Buildings 27 0 0.5 1 1.5 2 2.5 3 0 5 10 15 20

Normalized length and width

Frequency

(a)Lower bound,µlir= 1.01, σlir= 0.15.

0 0.5 1 1.5 2 2.5 3 0 5 10 15 20

Frequency (b)Upper bound,µmbr= 1.52, σmbr= 0.43. 0 0.5 1 1.5 2 2.5 3 0 5 10 15 20

Frequency

(c) mcmcmodel,µmcmc= 1.15, σmcmc= 0.14.

Figure 4.1:Histograms and fitted normal distributions of calculated length and width normalized to ground truth equal to one.

(38)

−15 −10 −5 0 5 10 15 0 2 4 6 8 10

Angle difference compared to ground truth

Frequency

(a)Angle from lir,µlir= 0.07◦,σlir= 4.12◦.

−40 −30 −20 −10 0 10 20 30 40

0 5 10 15

Frequency (b)Angle from mbr,µmbr= 3.17◦,σmbr= 12.05◦. −15 −10 −5 0 5 10 15 0 2 4 6 8 10

Frequency

(c)Angle from mcmc model,µmcmc= −0.18◦,σmcmc= 2.11◦.

Figure 4.2:Histograms and fitted normal distributions of the difference be-tween true building angle and predicted value for different methods.

(39)

4.8 Terraced Buildings 29 0 5 10 15 20 25 0 5 10 15 20 Width [m] Height [m] (a)Building 1. 0 5 10 15 20 25 0 5 10 15 Width [m] Height [m] (b)Building 2. 0 5 10 15 20 25 0 5 10 15 Width [m] Height [m] (c)Building 3. 0 5 10 15 20 25 0 5 10 15 Width [m] Height [m] (d)Building 4. 0 5 10 15 20 0 5 10 15 Width [m] Height [m] (e)Building 5. 0 5 10 15 20 0 5 10 15 Width [m] Height [m] (f)Building 6.

Figure 4.3: dhm(blue), primitive (orange) and ground truth (green) heights viewed along the model ridge line.

(40)

(a)Orthophoto and classification. (b) dhm.

(c)Data and primitive height contour. (d)Orthophoto and primitive height contour.

(e)Data surface. (f)Primitive surface.

(41)

4.8 Terraced Buildings 31

(c)Data and primitive height contour. (d) Orthophoto and primitive height contour.

(42)

(e)Data surface. (f)Primitive surfaces.

(43)

4.8 Terraced Buildings 33

0 10 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 Length [m] Height [m]

(g)Primitive side view.

(44)

(45)

5

Discussion

This chapter discusses and highlights important parts from previous chapters.

5.1 Results

The proposed method manages to successfully reconstruct buildings consisting of one or more primitives where the primitives are connected on the short side in accordance to the limitations set up in section 1.4.

The algorithm tends to favour higher and wider roofs than is the ground truth. This is probably because the buildings in the dhm are smudged, so a bigger roof probably minimises the cost function better. A possible solution to this would be to try to add another cost function, for example the average height deviation or the horizontal distance to the building walls.

Neither of the bound measures, mbr and lir, take into account the position, (x, y), of the building. The dimensions given by lir are very close to correct dimen-sions, but lir does not give enough information about the building to be fully reconstructed. Nevertheless the dimensions and angle results from lir are quite astounding as it is a fairly fast and simple method. It performs on par with the proposed mcmc algorithm on both dimension and angle estimation.

The mcmc algorithm has trouble separating gable and hipped roofs due to small height differences between them. Terraced buildings are currently not possible to reconstruct correctly since the algorithm favours one large primitive instead of several smaller primitives. Since the primitive angle can be identified with high accuracy from the Hough lines, a possible solution to this is to not allow the building primitive to move very far from the initialisation angle. This could

(46)

maybe force the algorithm to optimise with more primitives.

5.2 Method

The cost function currently used makes it difficult to separate large buildings, length ≥ 30 m, with gable roof or hipped roof with high confidence. Since the size of the hip is small compared to the rest of the building, the cost difference between the different roof types is often within the margin of error. This usu-ally does not happen with smaller houses, where the cost difference of the hip is significantly different between a gable roof and a hipped roof.

At the current state the orthophoto is utilized only in the visualisation. If one could use an edge detector to detect the ridge line and footprint and incorporate this information in the optimisation procedure it is possible that it would increase the algorithm precision and speed.

Some kind of height analysis should be part of the preprocessing to decide if the building is a flat or a more complex structure. In the dataset available, there were no true flat roofs to test on, either they were not flat enough in the real world, or the data was too noisy for the algorithm to make a good prediction. If the algorithm could decide beforehand if the building is flat or not, this could significantly cut the algorithm runtime.

A better design of the parameter proposals should also be implemented. Trun-cated normal distributions is not the ideal choice since it slightly breaks the pro-posal symmetry described in section 2.4. A log-normal or beta distribution could have been a better choice for some parameters since it is strictly positive and for the beta distribution bounded as well.

To increase the success of the algorithm in areas with more diverse and exotic buildings a more complex parametrisation could be required. If another height is added with more hips a lot more primitives is obtained. This algorithm could also be used with ellipse shapes to parametrise rounded buildings and roofs. During collection of ground truth data some problems defining the height of a building arose. Often, the ground underneath a building is sloped. The difference in height between two sides of the building could be larger than a metre. This makes the comparison from the reconstructed roof and ground truth somewhat unreliable. The dtm should take care of this, but some noise is unavoidable. To improve the optimisation procedure it would be beneficial to implement the Reversible Jump Markov chain Monte Carlo (rjmcmc) introduced by Green in 1995 [6]. This makes it possible to jump between the different parameter spaces of different roof types and the number of primitives effectively, which is also how Huang et al. [7] implemented their algorithm.

For buildings containing more than one primitive, a global optimisation should be part of the algorithm. If one of the primitives is a little off the correct angle

(47)

5.2 Method 37

it will propagate to the other primitives which often tries to line up with the al-ready placed primitive. With a global optimisation step, the primitives should align correctly. After the optimisation step a merging algorithm should be de-signed where discontinuities of ridge lines and walls is handled. Huang et al. [7] proposed a merging method which successfully merged primitives to avoid discontinuities.

(48)

(49)

6

Conclusions

This chapter presents conclusions regarding the performance and viability of the proposed method as well as thoughts on future work to improve the method.

6.1 Principal Questions

In the problem formulation two questions was asked:

• What measurement accuracy improvement is achieved after building recon-struction?

• Is the proposed algorithm computationally effective enough to be used for practical applications?

The algorithm correctly places building primitives at roughly the correct location, often the dimensions are a bit larger than the ground truth. The results show that the lower bound normal distribution have a mean µlir≈ 1 and an upper bound at around µmbr ≈ 1.5. The mcmc algorithm places the primitive at correct position with dimensions slightly larger, µmcmc ≈ 1.15, than the true value. Since the mcmccontains all information to fully reconstruct a building wheras lir does not, it is still the better option for building reconstruction.

Currently the execution time for the mcmc algorithm is around one minute per primive with a MATLAB™ implementation on a Xeon processor. With an L- or U-shaped building the runtime increases linearly with the number of primitives required to fully reconstruct the building.

(50)

6.2 Future Work

It would be interesting to change the cost function to incorporate the average deviation and maybe punish large deviations harder. A global optimisation and merging step is required for buildings consisting of more than one primitive to reconstruct with satisfaction. Since the lir method performed so well future work could utilize it in the initialisation step in the mcmc algorithm.

(51)

(52)

(53)

A

Building Data

1

66

10

2

48

12

3

10.5

51.5 _49.5

4

12.5

5

40.5

12

6

60.5

12 7 59

12

8

70

13

9

54

13 10 37

11

45

11

12

28

11

13

26

13

14

24

14.5

15

24

13

16

24

13

17

23

13

18

23

13

19

23

13

20

23

13

21

23 ₁₃

22

23 ₁₃

23

24 ₁₃

24

22 _10.5

25

22 _10.5

26

22 _10.5

27

22 _10.5

28

22 _10.5

Figure A.1:Larger building numbers and their dimensions in metres.

(54)

1

66

10

2

48

12

3

10.5

51.5 _49.5

4

12.5

5

40.5

12

6

60.5

12 7 59

12

8

70

13

9

54

13 10 37

11

45

11

12

28

11

13

26

13

14

24

14.5

15

24

13

16

24

13

17

23

13

18

23

13

19

23

13

20

23

13

21

23 ₁₃

22

23 ₁₃

23

24 ₁₃

24

22 _10.5

25

22 _10.5

26

22

10.5

27

22 _10.5

28

22 _10.5

(55)

45

(a)i = 1 (b)i = 1000 (c)i = 3000

(d)i = 1 (e)i = 1000 (f)i = 3000 Figure A.3:Building 7 for iterationi∈ {1, 1000, 3000}.

100 ₁₀1 ₁₀2 ₁₀3 0 0.2 0.4 0.6 0.8 1 Iteration C ost of primitiv e [%]

(56)

Table A.1:Building dimensions from the digital map, mcmc, mbr and lir.

Digital Map [m] mcmc[m] mbr[m] lir[m]

Building Length Width Length Width Length Width Length Width

1 66.0 10.0 69.1 17.2 69.8 16.6 67.5 10.5 2 48.0 12.0 49.2 13.6 52.7 15.7 49.0 10.5 3 51.5 10.5 46.6 13.2 52.9 30.7 50.0 11.0 4 49.5 12.5 52.5 16.2 57.3 19.1 51.0 12.0 5 40.5 12.0 41.3 13.9 43.0 15.4 40.0 11.0 6 60.5 12.0 59.9 13.4 59.5 13.9 56.5 8.0 7 59.0 12.0 59.9 16.5 64.4 20.4 61.0 13.5 8 70.0 13.0 72.2 16.3 79.0 23.4 69.5 11.0 9 54.0 13.0 56.2 15.5 65.0 33.3 57.0 12.5 10 37.0 11.0 46.9 15.1 56.6 22.0 50.0 10.0 11 45.0 11.0 39.2 12.9 44.1 17.1 39.0 9.5 12 28.0 11.0 28.9 14.5 30.9 16.5 26.0 12.5 13 26.0 13.0 28.2 16.0 41.1 18.8 35.0 14.0 14 24.0 14.5 22.6 15.4 28.3 20.7 24.0 15.0 15 24.0 13.0 24.6 15.1 31.7 27.2 26.0 12.5 16 24.0 13.0 25.3 15.1 33.9 23.1 19.5 16.0 17 23.0 13.0 25.3 15.0 35.5 24.9 28.0 13.0 18 23.0 13.0 25.5 15.5 40.0 33.1 26.0 15.5 19 23.0 13.0 26.7 16.0 34.5 22.0 19.5 16.5 20 23.0 13.0 25.9 16.4 36.0 21.9 26.5 16.0 21 23.0 13.0 23.9 15.2 28.0 27.8 20.5 16.0 22 23.0 13.0 25.8 14.9 28.6 22.6 23.0 15.0 23 24.0 13.0 26.6 15.4 31.2 29.9 23.5 13.5 24 22.0 10.5 24.3 13.9 22.1 22.1 20.0 10.0 25 22.0 10.5 23.0 15.0 25.5 15.2 23.0 9.0 26 22.0 10.5 24.0 14.7 39.2 22.3 23.5 9.0 27 22.0 10.5 23.6 15.2 25.2 14.5 20.0 11.5 28 22.0 10.5 23.7 12.1 28.0 12.5 26.0 5.0

(57)

47

Table A.2:Building dimensions from real world measurements, digital map measurements and differences between them.

Real World [m] Digital Map [m] Difference [m] Building Length Width Length Width Length Width

1 66.1 10.0 66.0 10.0 -0.1 0.0 2 47.8 12.1 48.0 12.0 0.2 -0.1 3 44.6 10.5 51.5 10.5 6.9 0.0 4 49.2 10.2 49.5 12.5 0.3 2.3 5 40.2 – 40.5 12.0 0.3 – 6 60.3 12.9 60.5 12.0 0.2 -0.9 14 24.1 14.5 24.0 14.5 -0.1 0.0 15 24.0 13.2 24.0 13.0 0.0 -0.2 22 23.2 13.3 23.0 13.0 -0.2 -0.3

(58)

Table A.3:Building angles from digital map. Angle, α, [◦_]

Building Digital Map mcmc lir mbr

1 12 13.4 12 12.0 2 12 12.0 12 12.2 3 12 11.8 12 36.3 4 12 12.1 12 12.7 5 15 16.5 17 15.1 6 12 10.0 9 11.1 7 12 12.1 12 14.4 8 15 15.4 16 13.8 9 15 14.2 17 17.5 10 11 12.1 10 11.9 11 12 10.9 14 8.1 12 12 14.8 16 7.1 13 18 16.8 22 16.7 14 8 3.2 6 4.1 15 9 6.3 4 3.0 16 9 8.4 15 13.3 17 8 10.6 8 21.5 18 8 5.6 8 20 19 8 8.0 10 0 20 8 10.5 16 12.1 21 45 42.4 38 17.9 22 45 44.5 35 17.0 23 45 48.5 47 16.0 24 31 32.3 27 11.3 25 31 33.5 31 26.6 26 31 35.2 33 14.3 27 31 32.6 35 35.5 28 39 37.2 30 35.8

(59)

49

Table A.4:True building roof types and mcmc predictions Roof Type

Building True Predicted

1 Hipped Gable 2 Gable Gable 3 Gable Gable 4 Gable Gable 5 Gable Gable 6 Gable Gable 7 Gable Gable 8 Gable Gable 9 Gable Gable 10 Gable Gable 11 Gable Gable 12 Gable Hipped 13 Gable Gable 14 Hipped Hipped 15 Hipped Hipped 16 Hipped Hipped 17 Hipped Hipped 18 Hipped Hipped 19 Hipped Hipped 20 Hipped Hipped 21 Hipped Flat 22 Hipped Hipped 23 Hipped Hipped 24 Gable Gable 25 Gable Gable 26 Gable Gable 27 Gable Hipped 28 Gable Flat

(60)

(61)

Bibliography

[1] Digital Globe – Vricon 3D and Elevation. URL https://www. digitalglobe.com/products/vricon-3d-elevation. Accessed: 2018-04-15. Cited on page 2.

[2] Vricon Company Overview. URL https://www.vricon.com/company/ overview/. Accessed: 2018-04-15. Cited on page 2.

[3] Fahmi Amhar, Josef Jansa, Christine Ries, et al. The generation of true or-thophotos using a 3d building model in conjunction with a conventional DTM. International Archives of Photogrammetry and Remote Sensing, 32: 16–22, 1998. Cited on page 5.

[4] Filip Biljecki, Hugo Ledoux, and Jantien Stoter. An improved lod specifica-tion for 3d building models. Computers, Environment and Urban Systems, 59:25–37, 2016. Cited on page 7.

[5] Robert P Dobrow. Introduction to stochastic processes with R. John Wiley & Sons, 2016. Cited on pages 8 and 9.

[6] Peter J Green. Reversible jump markov chain monte carlo computation and bayesian model determination. Biometrika, 82(4):711–732, 1995. Cited on page 36.

[7] Hai Huang, Claus Brenner, and Monika Sester. A generative statistical ap-proach to automatic 3d building roof reconstruction from laser scanning data. ISPRS Journal of photogrammetry and remote sensing, 79:29–43, 2013. Cited on pages 13, 14, 36, and 37.

[8] Johan Källström. Building and tree parameterization in partially occluded 2.5D DSM data. Master’s thesis, Linköping University, Department of Elec-trical Engineering, 2015. URL http://urn.kb.se/resolve?urn=urn: nbn:se:liu:diva-125130. ISRN: LiTH-ISY-EX–15/4917–SE. Cited on page 1.

[9] Scott Kirkpatrick, C Daniel Gelatt, and Mario P Vecchi. Optimization by simulated annealing. science, 220(4598):671–680, 1983. Cited on page 10.

(62)

[10] Thomas H Kolbe, Gerhard Gröger, and Lutz Plümer. Citygml: Interopera-ble access to 3d city models. In Geo-information for disaster management, pages 883–899. Springer, 2005. Cited on page 7.

[11] Florent Lafarge, Xavier Descombes, Josiane Zerubia, and Marc Pierrot-Deseilligny. Structural approach for building reconstruction from a single dsm. IEEE Transactions on pattern analysis and machine intelligence, 32(1): 135–147, 2010. Cited on page 1.

[12] Zhilin Li, Qing Zhu, and Chris Gold. Digital terrain modeling: principles and methodology. CRC press, 2004. Cited on page 6.

[13] Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Au-gusta H Teller, and Edward Teller. Equation of state calculations by fast computing machines. The journal of chemical physics, 21(6):1087–1092, 1953. Cited on page 10.

[14] WM Mularie. Department of Defense world geodetic system 1984, its defi-nition and relationships with local geodetic systems. National Geospatial-Intelligence Agency, Tech. Rep, 152, 2000. Cited on page 6.

[15] Alan R Parkinson, R Balling, and John D Hedengren. Optimization methods for engineering design. Brigham Young University, 5, 2013. Cited on page 11.

[16] Tahmineh Partovi, H Huang, Thomas Krauß, H Mayer, and Peter Reinartz. Statistical building roof reconstruction from worldview-2 stereo imagery. The International Archives of Photogrammetry, Remote Sensing and Spa-tial Information Sciences, 40(3):161, 2015. Cited on pages 13 and 14. [17] Christian P Robert and George Casella. Monte carlo statistical methods.

Springer-Verlag, New York, 3, 2004. Cited on pages 8 and 9.

[18] Jean Serra. Image analysis and mathematical morphology. Academic Press, Inc., 1983. Cited on page 16.

[19] Richard Szeliski. Computer vision: algorithms and applications. Springer Science & Business Media, 2010. Cited on page 1.

Building Reconstruction of Digital Height Models with the Markov Chain Monte Carlo Method

Master of Science Thesis in Electrical Engineering

Department of Electrical Engineering, Linköping University, 2018