DERIVATION OF BUILDING STRUCTURES FROM NOISY DIGITAL SURFACE MODELS

: In this work we present a novel approach for segementation of a noisy DSM to building structures and other non-building structures – normally trees – and the modeling of them. Mostly Digital Surface Models (DSMs) from only a few aerial images or only from one pair of satellite images tend to be very noisy and lack good quality especially in shadow areas. Since actual methods for deriving roofs rely on a valid height information by joining areas of same slope to a roof-plane these fail regularly with such noisy DSMs. In our presented approach we use a slope map of the DSM only to detect flat regions. Since those regions on top of roofs are mostly good illuminated we can derive the ridges of roofs and flat roofs and also ground areas. All narrow, flat, elevated areas are ridges and may occur on roofs or on trees. After connecting ridges in ridge-directions there remain two types of ridges: long, straight ridges of roofs and mixed short ridges in many directions for the trees. Fitting symmetric planes through the roof-ridge-lines gives finally the roof-planes reducing the effects of noise on shadowed parts of the roof. Taking the other tree-ridges as seeds for a watershed transformation will give the trees. Finally the proposed method is applied to a noisy DSM and the results will be discussed.


INTRODUCTION
Creating Digital Twins (DTs) of urban areas only from Digital Surface Models (DSMs) is mostly very challenging due to (a) much noise in DSMs created only from very few stereo images like e.g.only two views for satellite imagery and (b) the intermixture of different elevated objects like buildings and trees.Due to (a) it is normally very complicated to distinguish between building and tree objects.Existing methods are mainly based on the detection of areas of homogeneous slopes for the extraction of roofs from precise DSMs like LiDAR-DSMs (Orthuber and Avbelj, 2015).But in practice in DSMs mentioned above the surface of a slanted roof on the shadowed side consists only of noise without any detectable homogeneous slope.Also LiDAR pointclouds gives the great advantage in delivering first-and last-pulse information, which allows a easy discrimination of solid objects like buildings and vegetation like trees.
Previous works rely mostly on high quality DSMs like provided with aerial laser scanning as shown in (Brenner, 2003), (Rottensteiner and Briese, 2002), (Haala and Kada, 2010) or (Perera et al., 2012).Further analysis also including other urban objects like trees are discussed in (Haala and Brenner, 1999) or for discrimination of trees and buildings using also the RGBimages as in (Krauß, 2019).The newest methods incorporate deep learning methods and are mostly based on fusing height information from DSMs with image data like in (Partovi et al., 2017).
Deeper analysis of the available noisy airborne DSMs reveals interesting details in derived filtered slope-maps: Due to the noise mainly originating from matching errors in dark, shadowed regions, bright areas like flat roofs or the ridges of gabled roofs are very stable areas and such the derived slopes in these areas are also much more reliable than those from shadowed areas.So the here presented method is based on the extraction of flat areas from noisy DSMs.

METHOD
The presented method consists of four main parts which are described in detail in the following sub-sections:  Due to the noise in the DSM there are also many negative outliers below the ground surface.Hence approaches deriving a DTM from the "lowest neighbours" fail in this case dramatically.So we developed a method incorporating an outlier filtering to remove too extreme values in the noisy DSM.For the outlier filtering used in this work we calculate a mean-(µ) and standard-deviation-image (σ) of the noisy DSM using a filter radius of 2.0 m (with the given ground sampling distance (GSD) of 10 cm of the images this is a filter size of 41 px).Afterwards  all DSM values lower than µ − σ or higher than µ + 2σ are removed and interpolated from neighbours as shown in fig. 2. The resulting filtered DSM is changed too heavy for further processing but perfect to derive the DTM.
For deriving the DTM from the outlier-filtered DSM we use the improved classical morphological approach as described in (Krauß et al., 2011).The combination of a low percentile filter (5 %) followed by a heigh percentile filter (95 %) resembles a more robust morphological opening.For both a filter radius of 100 m is used.From this DTM and the DSM a so called normalized digital elevation model or nDEM can be derived as nDEM = DSM − DT M as shown in fig. 4.

Classification of flat areas
The second step in our processing chain extracts and classifies flat areas in the DSM as shown in fig. 3. First a slope/aspect map is calculated from the pre-filtered DSM.To pre-filter the DSM a median filter with a window-size of 5 is applied to level out the noise, followed by a Gaussian filter with σ = 2.0.The results show low slope areas on flat roofs and on ridges of slanted roof buildings.Extracting all areas with slopes below 10 • gives us the candidates of flat areas and the others ("slanted").Using the DSM and the previously derived DTM we can calculate a so called normalized DEM (nDEM) which contains the heights of objects above ground.Using this nDEM allows us to separate high flat areas from such on the ground ("flat low") using a height threshold of h = 2.5 m.The high flat (fig.5) areas are further subdiveded to "ridges" if the detected flat area is elongated and narrow (red in fig.5, right) and "flat roofs" otherwise (white in fig.5, right).The decision which flat area is low or high is based on the nDEM and a height threshold of e.g.2.5 m.So all flat areas below 2.5 m will become "low flat", the others "flat high" or "ridges".
To distinguish between narrow and not narrow high areas a binary opening with a structuring element of a small size like e.g.50 cm is applied.Ridges will vanish after the binary opening, larger areas will still remain.
The flat roofs can directly be modeled at this step.For this we apply iteratively a binary dilation and add all dilated pixels to the flat roof if their height is within a small threshold (e.g. 25 cm) around the mean height of the original flat high area.
The iteration stops if no more pixel could be added.
The result of this first step are flat areas of the three types "flat low", "flat high" and "ridges".

Modeling of slanted roofs
Afterwards we can derive the 3D vector models from the classified "ridges" as shown in fig.6.It starts with the transformation of a bitmap-ridges-map to a vectorized skeleton of all ridges as shown in fig.7, center.
Connecting ridge areas with neighbours of same orientation and height and lying in the correct direction allows for joining of discontinous ridge areas.The result is a graph containing all connected ridge-lines of one object (see fig. 7, right).Since we have a noisy DSM many of these ridge-lines will be only small sections continuing after a missing part, so the connecting step is mandatory.
Here a first quality check occurs by checking the consistency of the height values along these extracted 3D ridge lines.Getting the roof ridges splits the graphs to real roof-ridges (a few long, straight or rectangular lines in the graph) and tree-ridges (many small ridges in random directions).Afterwards the tree-ridges get removed from the roof-ridges.
For the remaining roof-ridges profiles perpendicular to valid ridge lines are extracted and merged along each straight skeleton-vector to build up the slanted roof from planes on both sides of the ridge (see fig. 8, left).Therefore height profiles parallel to the ridge line in a distance xi are extracted until the median of the height profile is lower than the DTM plus the height-threshold or the step from the previous to the actual profile is positive (raising again) or negative and steeper than 1 m.These extracted height values for the planes are fitted through the roof-ridge-line to give the slant angle of the slanted roof.
For this we minimize the errors of a linear function yi = axi + b with yi as the measured heights from the DSM and xi the perpendicular distances to the ridge line.a gives the slope, b the offset fixed to the mean height of the ridge-line.So we will have to minimize the error ε in eq. 1 for a general fixed point (x f ix , y f ix ) giving b = y f ix − ax f ix (in our special case the ridge line is defined as x f ix = 0 and y f ix as the mean height of the ridge line).
The minimum of the error ε will be given by using a value of a which can be calculated from the derivative with respect to a of eq. 1 as shown in eq. 3.
So we have finally fitted two planes -left and right -through the ridge line with the common length of the ridgeline l, a individual width w l,r and slope a l,r .Now these planes undergo a second quality assessment.First the planes should have a width of more than 1 pixel and they should not be wider than 1.2 times the length of the ridge.Second the planes should be about the same width (one plane not being smaller than 0.2 the width of the other plane).Third the values for the slopes a of the planes left and right of the ridge should not differ more than 10 • and both should be larger than 10 • going downward (not flat and descending).
If all these quality checks are passed the planes are reduced to the minimum width of both (w = min(w l , wr)) and the slopes are set to the mean of the both slopes (a = (a l + ar)/2).Finally the four corners of each plane are derived as the two ridgepoints at mean ridge-height H and the two eaves points at the height of the eaves derived from the plane width and slope as h l,r = H + aw.The result for the region shown in fig.7 is given in fig.8.

Modeling of trees
To model the trees first a inverted watershed-transformation (Beucher, 1982, Haithcoat et al., 2009) is applied to the Gaussian filtered (σ = 5.0) noisy DSM providing a segmentation by the local maxima as shown in fig.6, upper part and in fig. 10, left.
For each of the watershed segments a high-mask (e.g. with a height of 4 m above the DTM) is calculated.If the class "not flat" dominates this high-mask it is taken as a tree.This is shown in fig. 10, center which adds a class "tree-top" (in green) to the flat classification shown in fig. 5.This new class -the green areas -represent the areas higher than 4 m above ground in each watershed segment.The center of gravity of each of these areas gives the trunk position.Using the area A of the high mask gives the radius of the tree as r = A/π.The height of the tree is given by the nDEM height at the tree position.
The derived vectorized trees (as circles), buildings (as rectangles), flat roofs (in grey) and flat ground areas (in sand color) are finally shown in fig. 10, right.

Data
The method is applied to a noisy DSM from a flight campaign over Braunschweig which was flewn on 2020-07-31 from 9:30 to 11:30.As test area an area of 180 × 220 m 2 around the Tostmannplatz is selected as shown in fig.5, left.The corresponding result is shown in fig. 10, right.
The test area was selected such it contains smaller and larger normal slanted roof buildings, flat roof buildings, trees and also two challenging buildings: the church "Dankeskirche" just below the center of the area and the school "Astrid-Lindgren-Schule" in the northern part of the area.

DISCUSSION
As shown in fig. 10 the result of the presented method gives alread very good results.Most of the buildings are modeled correctly, the tree areas are also mapped correctly.But some problems and possible improvements to the method can be detected.Fig. 11 shows the area around the "Dankeskirche".Due to a missing clearly detectable ridge line on top of the church-tower (see terrestrial image in fig.2, right) the detection as roof-ridge fails and the remaining slanted (tower roof) area will be classified as tree-tops.North and south of the church small hedges and bushes get classified as flat or small slanted roofs.The large stairs landing in front of the church (east) raise nearly 2 m above ground and is such mis-classified as "flat roof" (cf. the terrestrial image in fig.11, right).
The latter two errors originate from the fixed height threshold of 2.5 m above the DTM to distinguish between high and low flat areas and in turn for excluding the low flat areas from modeling.Using a lower value will remove also most of the garages from the modeling.Due to the noisy DSM also the DTM has some virtual "sinks" (see fig. 4, left: left of the church) which result in too high values for the heights of some objects in the nDEM.Fig. 12 shows a detail containing smaller and larger buildings with dormers (and a garage, bushes and hedges).These dormers are mostly missed completely (top, on the small buildings) or mis-classified as tree-tops (bottom, left on the large building).Fig. 7 shows the skeletonization results of the small building on the top right.The classification of the ridges of the dormers is there still existing.But later-on they are dropped due a missing special handling of ridges contained in roof segments.Such an extension to the algorithm may be added in future work.The correct calculation of vector intersections of roofs may be a main extension to the method in future work.
Fig. 14 shows a detail of the complex roof of the "Astrid-Lindgren-Schule" located in the top of the test area.As can be seen both a T-shaped join as well as a zig-zag-path of roof ridges can not be modeled correctly.The T-shape is replaced by the longer, dominating roof while the zig-zag-shape get's a strange mixture of a hipped end overlaying a slanted roof part.How such rather complex roof shapes will be modeled is also left to future research until now.
Summarizing the discussion of the results the following points are worth remarking: • Hipped roof parts are often missed or wrongly detected • Problems arise with complicated roofs like at the school building • The method cannot detect a roof if no flat ridges can be derived (tower of the church) • The vectorization of flat roofs is missing in the method until now • The handling of dormers has to be included -ridges inside of detected roof segments • Some small areas are wrongly detected as flat or slanted roofs (e.g. the hedges around the church or the stairslanding of the church) • A correct joining of roof-parts is still missing (e.g. the building in the lower right edge)

CONCLUSIONS AND OUTLOOK
In this paper we presented a method to derive building-and tree-objects based only on a noisy digital surface model (DSM).
The process derives from a noisy DSM a digital terrain model (DTM), flat roof areas, slanted roof areas, trees and low flat areas.The first results of this work look already very promising but there is still some further research needed.
Future work should include the vectorization of the flat roof areas and the proper detection of hipped roofs and dormers.Afterwards also a correct intersection of all detected roof planes should be implemented.Also the mis-classification of small areas to flat or slanted roofs should be solved.
A major gain of the quality of the results may be achieved by also including true ortho images into the classification.Using such ortho imagery as shown in fig.15 may give additional hints for ridge locations or the discrimination of actually wrongly classified small hedges or other objects.But due to the noisy DSM also some of the noise may be introduced to the ortho images -especially if off-nadir imagery is used (cf. the gable ends of the church in fig.15).
For segmentation of a noisy DSM to building structures and other non-building structures -normally trees -we start with a preprocessing which gives us beneath the noisy DSM also a smooth DTM -a digital terrain model -representing the ground of the area as shown in fig.1.

Figure 1 .
Figure 1.Preprocessing of the DSM and derivation of a DTM (images and elevation models are depicted as rounded squares, processing steps as rectangles)

Figure 4 .
Figure 4. Same detail region as in fig.2, left: derived DTM, right: nDEM containing only elevated objects

Figure 5 .
Figure 5. Classification of the DSM to flat areas, left: DSM, right: derived flat classes (blue: not flat, sand: flat low, white: flat high, red: narrow flat and high, ridges)

Figure 6 .
Figure 6.Modeling of slanted roofs and trees (vector objects are shown as yellow squares)

Figure 7 .
Figure 7. Skeletonization of ridges, left: class image with ridges in red, center: extracted skeleton, right: joined graph

Figure 9 .
Figure 9. Section (approx.50 × 60 m) of DSM (left) and resulting classification (right) located in the left top edge of the test area As shown in fig. 9 flat areas remain still to be vectorized.Most of the hipped roof ends (top right of fig.9) are missed until now.Complex trees get splitted to many smaller trees (left top) but some trees overlaying the flat roofs are detected correctly (lower right).

Figure 10 .
Figure10.Results of the presented method based on the test area, left: DSM segmented using the watershed-transformation, center: classification results after tree-segmentation (a class "tree-top" (green) is added for all high and non-flat objects per segment), right:results (buildings and trees as vector-objects, remaining classes: gray: flat roofs, sand: low flat areas, blue: not flat areas)

Figure 13 .
Figure 13.Section (50 × 50 m) of DSM (left) and resulting classification (right) at the southern edge of the test area