- Email: [email protected]
Segmentation and object-based classification for the extraction of the building class from LIDAR DEMs
ARTICLE IN PRESS
Computers & Geosciences 33 (2007) 1076–1087 www.elsevier.com/locate/cageo
Segmentation and object-based classification for the extraction of the building class from LIDAR DEMs$ George Miliaresisa,, Nikolaos Kokkasb a
Department of Geology, University of Patras, Rion 26504, Greece Department of Geomatic Engineering, University College London, Gower Street, London WC1E 6BT, UK
b
Received 10 August 2006; received in revised form 22 October 2006; accepted 28 November 2006
Abstract A new method is presented for the extraction of a class for buildings from light detection and ranging (LIDAR) digital elevation models (DEMs) on the basis of geomorphometric segmentation principles. First, seed cells and region growing criteria are specified. Then an object partition framework is defined on the basis of region growing segmentation. Size filtering is applied to objects and connected components labelling identifies background and foreground objects that are parametrically represented on the basis of elevation and slope. K-means classification reveals a set of clusters. The interpretation of the spatial distribution of clusters assisted by the interpretation of cluster centroids, allows for the identification of the building class, as well as building sub-classes with different geomorphometric characteristics. r 2007 Elsevier Ltd. All rights reserved. Keywords: Digital elevation model; LIDAR; Region growing; Urban segmentation; Terrain modelling; Terrain classification
1. Introduction High-resolution (spacing is less than 2 m) digital elevation models (DEMs) derived by airborne laser altimetry (LIDAR: light detection and ranging) are available for selected areas.1,2 The development of airborne laser scanning goes back to the 1970s (Irish $ Code available from server at http://www.iamg.org/ CGEditor/index.htm. Corresponding author. Tel.: +30 2610 996296; fax: +30 2610 991 900. E-mail addresses: [email protected] (G. Miliaresis), [email protected] (N. Kokkas). 1 Bare Earth Digital Elevation Model Generated from LIDAR Data. University of IDAHO, http://inside.uidaho.edu/ default.htm. 2 Maryland Department of Natural Resources, http://dnrweb. dnr.state.md.us/gis/data/lidar/.
0098-3004/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2006.11.012
and Lillycrop, 1999). By emitting a laser pulse and precisely measuring the return time to the source, the range can be calculated using the value for the speed of light.3 Elevations can be derived rapidly and at high resolution from LIDAR in comparison to manual reconstruction from photogrammetric techniques (time consuming, not a cost-effective solution). The LIDAR DEMs are extremely valuable in the field of earth sciences since elevation data from different flights might be compared to determine the patterns and magnitudes of coastal change4 3
LIDAR Tutorial, NASA, http://wwwghcc.msfc.nasa.gov/ sparcle/sparcle_tutorial.html. 4 Remote Sensing for Coastal Management, NOAA, http:// ekman.csc.noaa.gov/TCM/.
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
(erosion, over wash, etc.) and the loss of buildings and infrastructure.5 LIDAR DEMs are used in tectonic studies for fault recognition (Harding and Berghoff, 2000) and topographic change mapping.6 Additionally, airborne laser altimetry data were also used in an attempt to examine and model the surface morphology of landslides (Glenn et al., 2006). The building class extraction from high-resolution urban DEMs (spacing o2 m) is of primary importance in many applications, including urban planning, telecommunication network planning and vehicle navigation which are of increasing importance in urban areas (Kokkas, 2005). Nevertheless, the building class mapping is of great importance in the earth science field, for earthquake damage assessment. For example the difference DEM with 5 m spacing derived from pre- and post-earthquake aerial images revealed successfully the collapsed buildings, caused by the 1999 Izmit earthquake in Turkey (Turker and Cetinkaya, 2005). The critical element of detecting earthquake-induced heights changes through the DEM differencing method, was deciding where to place the boundaries between change and no-change pixels (Turker and Cetinkaya, 2005). Thus, if a building class was detected from the pre-earthquake DEM then the elevation differences caused by city ruins over the nonbuilding class (roads, city parks) could be omitted. Methods for building detection and reconstruction from laser altimeter data are applied to the point cloud (irregularly spaced points) and either directly derive the surface parameters in a parameter space by clustering the point cloud or segment a point cloud based on criteria like proximity of points or similarity of locally estimated surfaces (Vosselman et al., 2004; Heuel et al., 2000). Kokkas and Dowman (2006) introduced a method that employs a semi-automated technique for generating the building hypothesis by fusing LIDAR data with stereo matched points extracted from the air– photograph stereo model. The roof reconstruction is achieved by implementing a least squares-plane fitting algorithm on the LIDAR point cloud and subsequently neighbouring planes are merged using Boolean operations for the generation of solid features. 5
Hurricane and Extreme Storm Impact Studies: Coastal and Nearshore Mapping with Scanning Airborne Laser, US Geological Survey. 6 Active Tectonics Research Group, Arizona State University, http://lidar.asu.edu/data.html.
1077
On the other hand, quantitative techniques have been developed in order to automate the interpretation of terrain features from DEMs and various geomorphometric parameters were developed in an attempt to characterize the landscape (Miliaresis and Kokkas, 2004). Geomorphometric segmentation extracted mountains developed on different base levels (Miliaresis, 2001a) and fluvial landforms developed along mountain fronts (Miliaresis, 2001b) from moderate resolution DEMs. The technique is summarized as a region growing segmentation algorithm that uses seeds (e.g. the ridge cells) and growing criteria (usually based on slope and elevation). Whereas geomorphic segmentation defines both abstractions of landforms and an object partition framework of the landscape, the parametric representation of landscape objects on the basis of their spatial three-dimensional arrangement (Miliaresis and Argialas, 2002) allows the recognition of landscape objects) that have similar and contrasting ranges of characteristics, leading to a terrain classification scheme (Miliaresis and Illiopoulou, 2004) according to geomorphologic principles and understanding (Miliaresis et al., 2005). The aim of the current research effort is to design a new method for extraction of the building class from LIDAR DEMs on the basis of the geomorphometric segmentation principles. Thus, an object partition framework of the LIDAR DEM will be defined. Finally, objects representation on the basis of geomorphometric parameters combined by object classification is expected to allow the extraction and mapping of the building class. 2. Methodology First, the study area and its hypsometric characteristics are introduced. Then seeds cells and region growing criteria are defined. Finally, region growing segmentation is performed and thus, an object partition framework is defined. Connected components labelling identifies background and foreground objects that are parametrically represented on the basis of elevation and slope and classified. 2.1. Study area and data The LIDAR data used in this paper data were acquired by NPA Group using their in-house build LIDAR sensor during the leaf-on season. The
ARTICLE IN PRESS 1078
G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
Fig. 1. LIDAR DEM of study area and elevation frequency histogram. Elevation is in range 12– 215 m, the lighter a point is the greatest its elevation.
LIDAR point cloud represents the Bloomsbury area in central London bounded by co-ordinates (lower left: 529,174, 181,500, upper right: 530,199, 182,525) in the British National Grid. The point cloud has been post-processed and delivered with a 1 m 1 m (Fig. 1) sampling resolution covering an area of approximately 1 km2 (1024 rows and 1024 columns consist the derived regular grid DEM). The elevation frequency histogram (Fig. 1) indicates that the majority of DEM cells (98.94%) present elevation in the range 20–60 m. 2.2. Seed cells Median filtering was applied for seed pixels to be defined. The seed points were computed in a limited size kernel window. More specifically, the kernel size is 5 5 grid cells; note that each grid cell is 1 m apart. Due to the very small kernel size, it is impossible for the seed points to be affected from the regional topographic surface. Additionally, the regional topographic surface in the study area is a gently sloping plane. Even if the regional topographic surface was more complicated (either 2 or 3 order polynomial surface), the elevation differences computed in a 25 m2 area would be residual elevation anomalies and not regional elevation anomalies especially in an urban (city) environment. Median filtering removes very high and very low elevation values within the selected kernel size while
at the same time preserves edges (Mather, 1987). Then the Difference DEM (LIDAR DEMMedian filtered DEM) was derived. The Difference DEM is depicted in Fig. 2a. The Difference DEM revealed the cells of the LIDAR DEM where severe median filtering occurred. Note that elevation differences less than 1 m are depicted white while differences greater than 9 m are depicted black). The frequency of the elevation differences (Fig. 2b) indicates that the majority of DEM cells (97.1%) present difference in the range [4, 11]. The distribution is positively skewed (greater frequencies occur for positive elevation differences with respect to the corresponding negative elevation differences). Positive differences were interpreted to be associated to building edges. As it can be seen in the Figs. 2c and d, the greatest the elevation difference, the more the detected cells stands along the borders of the interpreted city buildings. Under a trial and error procedure though experimentation and visual interpretation of the results, was found that positive elevation differences in the range 5–9 m are more associated to the building edges (Fig. 2d). On the contrary differences greater than 9 m were interpreted to be associated to the internal building/roof structure (chimneys, statues, etc.). Thus, as seed cells were selected those presented elevation difference within the range of 5–9 m.
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
1079
Fig. 2. (a) The difference DEM (LIDAR DEM- Median filtered DEM). (b) Frequency of elevation differences. (c) Enlarged portion of difference DEM (differences less than 1 m are depicted white while differences greater than 9 m are depicted black). (d) Enlarged portion of difference DEM (differences less than 5 m are depicted white while differences greater than 9 m are depicted black).
2.3. Region growing segmentation For region growing segmentation, region growing criteria should be defined (Miliaresis and Argialas, 1999, 2000). In previous research efforts, aimed to landform segmentation (Miliaresis, 2001a b, 2006) the criteria were defined on the basis of slope and elevation. If the partial derivatives of elevation (H) along the east (x) and the north (y) direction are known then then slope and aspect (slope pointing direction) are computed from Burrough, 1987). sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 dH dH Slope ¼ , (1) þ dx dy Aspect ¼ arc tan
dH=dy . dH=dx
(2)
In order to estimate the partial derivates of elevation over a 3 3 kernel, Evans (1980) fitted a 6-parameter quadratic equation while Zevenbergen and Thorne (1987) used a 9-parameter quadratic equation. Evans (1980) method results to a polynomial surface which will not necessarily pass through the 9 grid cells included in the 3 by 3 kernel, while the opposite is true for Zevenbergen and Thorne (1987). A review and a comparison of these methods (Skidmore, 1989) proved that Evans (1980) method is more precise since Zevenbergen
and Thorne (1987) method is affected more from random noise and errors evident in DEMs. The slope frequency histogram (Fig. 3b) is completely different than the slope frequency distributions, we are used to observe in geomorphologic applications (the slope frequency is decreasing gradually for increasing slope values) where digital terrain models (DTMs) are used (Miliaresis and Kokkas, 2004; Miliaresis et al., 2005). DTM depicts ground elevation (the height of artificial structures like buildings is ignored) while in a DEM the elevation of the highest object (height of the trees and buildings) above the ground is included (Maune, 2001). So in DTMs (with spacing 30–1000 m), the ground elevation prevails and thus the slope frequency distribution is the outcome of geomorphologic (erosion, deposition) and geologic (tectonism) processes (Miliaresis, 2006). On the contrary in urban LIDAR DEMs, the building structure prevails and due to the high DEM resolution (p1 m), the slope frequency distribution we observe is the outcome of the interaction of the building surface (abrupt changes of elevation, high and abrupt changes of slope). The study of both the slope image (Fig. 3a) and the slope frequency histogram (Fig. 3b) combined by a trial and error procedure for selecting slope cut-off values, through experimentation, assisted by the photo-interpretation of the results, indicated
ARTICLE IN PRESS 1080
G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
Fig. 3. Slope image (the higher the slope value, the darkest a cell is). Frequency of slope in the interval 1– 901, (classes for slope o11 were excluded due to its great occurrence).
that the region growing criteria are defined for slope values in the interval [401, 901]. Then, an iterative region growing segmentation algorithm (Miliaresis, 2001a) was implemented. The seeds (Fig. 4) formed the initial set of building cells. If a cell of the DEM with slope in the interval [401, 901] was an 8-connected neighbour to the current set of cells, it was flagged as a new building cell and the current set of building cells was updated. Note that the 8-connected neighbours are the 8 pixels surrounding the central pixel in a 3 3 kernel. The segmentation stopped if no more cells were added during the current iteration. Finally, 448,865 cells were identified after 29 iterations (Fig. 4a). A connected component-labeling algorithm (Pitas, 1993) was applied and regions formed by 8connected background cells (black pixels) were identified and labelled. Many objects (the most of them very small in size) are evident in the segmented image. The very small-sized objects correspond to residual objects evident either on the roofs or in the non-building class. One could ask ‘‘what is the minimum size for an object to be clasified as a residual one?’’ an answer could be that building blocks in the study area should be far greater in size
than the ‘‘residual objects’’! In this particular case study (LIDAR DEM portion of the central London), we assume that the building blocks should be far greater in size than 120 m2 (120 cells). Small in size foreground (white) and background (black) objects were eliminated by applying object size filtering. More specifically, foreground objects with size less than 120 cells were merged to the background class (Fig. 4b). Six thousand six hundred and thirty-three cells (1.48% of the foreground cells class) cells were merged to the background and so the foreground cell class is formed by 441,828 cells. Then, 114,009 (18.8%) of the background class cells (black cells) were merged to the foreground (Fig. 4c) for an object size threshold of 120 cells (thus the foreground cell class is formed by 555,837 cells). Finally, the remaining foreground objects (Fig. 5a) and background objects (Fig. 5b) were labelled by the connected components algorithm (76 foreground and 291 background objects were identified) and merged to a single image (Fig. 5c). Region growing segmentation, connects components labelling, size filtering and the subsequent object parametric representation were performed with GeoLogic Shell software (Miliaresis, 2001b)
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
1081
Fig. 4. Segmentation and object size filtering. (a) Cells segmented (white pixels) by region growing segmentation algorithm (left) and 448,461 (42.77%) cells belong to foreground objects class. (b) Foreground objects (formed by white pixels) filtering with size less than 120 cells (middle image). (c) Background objects (formed by black cells) filtering with size less than 120 cells (right image).
Fig. 5. Foreground and background objects after size filtering. (a) Seventy-six foreground objects in left image (the darker an object is the greater the ID value is), (b) 291 background objects in middle image (the darker an object is the greater the ID value is), (c) merged image that contains background and foreground objects (right image). IDs 1–76 correspond to foreground objects while IDs 77–367 corresponds to background objects (the lighter an object is the greater the ID value is).
that is freely available for download on the internet from the WEB site of the International Association for Mathematical Geology.7 7 Computers and Geosciences, software code section, http:// www.iamg.org/CGEditor/cg2001.htm. Geologic Shell software that is part of Miliaresis (2001b) paper.
2.4. Parametric representation The objects in Fig. 5c are interpreted to belong either to the building class or to the non-building class (roads, vegetated areas, components of a building, etc.). Additionally, objects that are adjacent or the one encloses the other could both belong
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
1082
to the building class. It is assumed that if each object is parametrically represented by a set of geomorphometric attributes then unsupervised classification could define clusters of objects that belong either to the building or the non-building class. Four attributes were used for the parametric representation of objects (Fig. 5c):
H, mean elevation that was computed as the average elevation of the cells that belong to an object’s region. R, roughness, the standard deviation of elevation of the cells forming the object (a measure of the vertical variability of object’s height). S, object’s mean slope. SR, standard deviation of slope (a measure of the variability of object’s slope).
2.5. Classification The crosscorrelation between attributes is given in Table 1. It is observed that a partial correlation exists between S and SR, as well as R and SR. In other words, if the slope is high that for the most of the objects standard deviation of slope and elevation should be high, but this is not always true. Cluster analysis is a multivariate procedure, which is commonly used for regional classification. It is based on some measurement of distance among objects (for example Euclidean distance), which is calculated in a c-dimensional space, where c represents the number of attributes used in the clustering process (Mather, 1987). The centroid method was employed that requires a priori definition of the number of clusters. The main problem in cluster analysis has to do with the different value range of the different attributes (for example elevation and slope) included in the distance calculation. In order to solve this problem, the data were standardized (transformed to a new dataset in which the corresponding attribute mean value is
Table 1 Attribute correlation matrix H H R S SR
R
S
SR
1 0.56 0.63
1 0.72
1
1 0.28 0.12 0.03
zero and the corresponding attribute std. dev. is 1) before the exploratory (K-means) clustering algorithm was applied (Miliaresis and Argialas, 2002). Thus, the data values were transformed by subtracting the corresponding attribute mean value and dividing the outcome by the corresponding attribute standard deviation (Clark and Hosking, 1986). Ten clusters of objects were considered through experimentation and a trial and error procedure. Each object was assigned to a cluster on the basis of its Euclidean distance from the cluster centroid. Their cluster gravity centres are given in Table 2 while the graphical representation of the gravity centres was computed on the basis of the standardized co-ordinates (Fig. 6). The numbers of objects that belong to each cluster are also given in Table 2. The similarity among cluster centroids could be assessed either by the graphical representations of centroids (Fig. 6) or by divergence (Table 3). The cluster compactness (Table 3) was computed on the basis of compactness ¼ 100
mean st: dev: , Mean
(3)
where mean and std. dev. are the mean distance and the corresponding standard deviation of the objects from the gravity centre of the cluster. In order to evaluate the existence of statistically significant differences in the distances among the cluster centres the statistical technique of one-way analysis of variance (ANOVA) was used. ANOVA calculates the between-cluster mean square and the within-cluster mean square (Miliaresis and Illiopoulou, 2004). The ratio of these two mean squares is Table 2 Final cluster (standardized) centroids and the number of objects in each cluster (sixth row) Cluster ID
1 2 3 4 5 6 7 8 9 10
Standardized centroid co-ordinates H
R
S
SR
Number of objects
3.2498 1.1694 0.1284 4.2958 0.324 0.4666 0.6108 4.9048 6.3596 0.9844
8.5073 0.3474 0.3539 9.4762 0.932 2.5654 0.8661 0.1302 8.7856 0.3556
0.4618 0.5844 0.37 1.9155 1.638 2.2769 2.3381 0.3923 0.1835 0.5615
2.2419 0.3778 0.5094 2.3725 1.9751 1.3411 0.0694 0.025 2.2435 0.4224
1 81 119 1 56 2 19 2 1 85
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
the usual ANOVA F statistic F¼
variation between groups sum of squares=ðk 1Þ , variation within groups sum of squares=ðN kÞ
(4) where k is the number of groups and N is the total sample size. The logic of the test argues that for the class (cluster) means to be significantly different, interclass (between cluster) variance must be of much greater magnitude than intraclass (within clusters) variance (Clark and Hosking, 1986; Shaw and Wheeler, 1985). The results of the ANOVA procedure for the each centroid are shown in Table 4. Although the significance levels for the F statistic are below 0.001, the values of F in the case of cluster analysis should be mostly interpreted in terms of comparing
1083
differentiation of the cluster means for individual variables. It is concluded that clusters are grouped to two categories (a) clusters 2, 3 and 10, (b) clusters 5 and 7 on the basis of similarity/divergence (Table 3, Fig. 6). In the first category, the clusters 2, 3, and 10 are differentiated on the basis of mean elevation while the second category is formed by objects with almost equal mean elevation. Clusters 7 and 10 are less compact (Table 3) meaning that they could be divided possibly to sub-clusters. The spatial distribution of clusters is depicted in Fig. 7 in an attempt to assist the interpretation.
2.6. Discussion of the results Fig. 8 presents the nature of the LIDAR elevation data, which is quite noisy while at the same time residual elevation objects are observed over the roofs. These objects correspond either to arificial objects induced by noise effects and errors evident in the LIDAR DEM or to small elevated features that are difficult to recognise and assigned to a certain feature type, due to both the DEM resolution and the local elevation errors that are present.
Table 4 Analysis of variance
Fig. 6. Centroids (standardized) for clusters 2, 3, 5, 7 and 10 that include majority of objects (Table 2).
Variable Cluster MSa DFa Error MSa DFa F
Prob.
H R S SR
0 0 0 0
a
36.14 39.03 36.79 34.13
9 9 9 9
0.114 0.041 0.097 0.164
357 357 357 357
316.83 943.36 376.36 206.94
MS, mean square; DF, degrees of freedom.
Table 3 Distances between cluster centroids (divergence) and cluster compactness for clusters 2, 3, 5, 7, and 10 Distance between cluster centroids
Cluster
2
3
2 3 5 7 10
0 1.32 3.58 3.24 2.15
0 3.47 3.12 0.88
Statistics for the objects distance from the corresponding cluster centroid 5
0 2.05 3.74
7
0 3.56
10
Cluster
Mean
St. dev
Compactness %
0
2 3 5 7 10
0.570 0.444 0.654 0.634 0.614
0.263 0.217 0.341 0.381 0.383
53.8 51.2 47.8 39.9 37.6
ARTICLE IN PRESS 1084
G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
Fig. 7. Spatial distribution of 10 clusters with labels 1 to 10 is given in top left image (the greater the cluster ID the darkest an object is). In subsequent images, clusters 2, 3, 5, 7 and 10 (each one is depicted white) are presented in an attempt to assist the interpretation.
Fig. 8. LIDAR DEM visualisation. (a) Direction of line of sight (NW-SE) that is determined by position of observer (NW) and position of target (SE). Line of sight is superimposed over LIDAR DEM. (b) LIDAR DEM visualisation along line of sight.
Thus, by the size filtering implemented (Fig. 4), the building blocks were generalised with respect to roof appearance.
The visual interpretation (Fig. 7) indicates that the non-building class is formed by cluster 2 while the urban vegetation class corresponds to cluster 7.
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
The interpretation of cluster centroids (Fig. 6, Table 2) indicate that the non-building class is formed by the less elevated, most flat, less variable with respect to slope and elevation objects. The urban vegetation class is formed by objects with the highest slope and with significant slope variability (standardized centroid is 40). On the other hand, the mean standardized elevation is below 0 while roughness is high (1). Thus, the vegetation class is formed by objects with low mean elevation (higher than the elevation of the non-building class), while slope is maximized and the standard deviation of slope and elevation is high. The last remarks are functions of tree canopy properties and thus different vegetation types might present different parametric representation. The building class is formed mainly by the objects of cluster 5 (Fig. 7). For the objects of the building
1085
class slope variability is maximized (Fig. 6) while slope is among the maximum ones (Table 2). The interpretation given is that it happens due to the building edges that are in contact with the nonbuilding class (Fig. 7). Mean elevation (H) and roughness (R) are higher than the corresponding values of both building and vegetation classes. The clusters 10 and 3 correspond mainly to internal portions of buildings. That is why the roughness, slope and slope variability is among the minimum while mean elevation is maximized (Fig. 6). The outliers (clusters 1, 4, 6, 8 and 9) are interpreted on the basis of cluster centroids (Table 2). They certainly belong to the building class due to the extremely high mean elevation values presented and they are interpreted to correspond to elevated objects evident mainly in building roofs (chimneys, etc.).
Fig. 9. Visualisation of segmented building and vegetation classes. (a) Segmented building classes and vegetation classes near a city park as wells as direction of line of sight (S-N). (b) Visualisation of segmented classes (in city park) that were superimposed over LIDAR DEM. (c) Visualisation of LIDAR DEM in city park. (d) A portion of Fig. 5c that contains 367 segmented objects as wells as direction of line of sight (SE-NW). (e) Visualisation of the segmented objects that were superimposed over the LIDAR DEM.
ARTICLE IN PRESS 1086
G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
In Fig. 9a selected segmented building and vegetation objects are presented. The visualisation of these classes (Fig. 9b) as well as the visualisation of LIDAR DEM (Fig. 9c) indicate the two portions to which the most of the buildings were segmented to (the internal portion that is depict with the white tone and the building border that is depicted with a shade of gray) as well as the segmentation of the city park vegetation. These findings are also supported by the visualisation of the Fig. 5c of the segmented objects (before the objects merging according to their parametric representation to 10 clusters/building classes) in Figs. 9d and c. 3. Conclusion Geomorphometric region growing segmentation combined by median filtering for identifying seed cells, connected components labelling, size filtering and object labelling, object parametric representation on the basis of slope and elevation attributes and classification was proved capable of delineating the building class within the study area. The interpretation of cluster centroids allowed the identification building sub-classes with different geomorphometric characteristics that are associated with different building portions and present different hazard assessment. Difference DEM derived from pre- and postearthquake LIDAR DEMs can be used in order to reveal the earthquake collapsed buildings, if the building class is detected from the pre-earthquake DEM according to the method described. Nevertheless, the proposed method for building class detection has a major disadvantage. It requires a certain level of user interaction for some crucial parameters, which are very difficult to define automatically in varied situations and further research is required for full automation. Acknowledgements The authors are grateful for, and this paper was greatly benefited from, the thorough and evaluations of the two anonymous reviewers. The authors would like to thank NPA Group for providing the LIDAR data and express their gratitude for the financial support offered by the Greek Scholarship Foundation. References Burrough, P.A., 1987. Principles of Geographical Information Systems for Land Resources Assessment. Oxford University Press, Oxford, 194pp.
Clark, W., Hosking, P., 1986. Statistical Methods for Geographers. Wiley, New York, 518pp. Evans, I., 1980. An integrated system for terrain analysis and slope mapping. Zeitschrift fuer Geomorphologie N.F. Suppl. -Bd. 36, 274–290. Glenn, N., Streutker, D., Chadwick, J., Thackray, D.G., Dorsch, S., 2006. Analysis of LiDAR-derived topographic information for characterizing and differentiating landslide morphology and activity. Geomorphology 73, 131–148. Harding, D.J., Berghoff, G.S., 2000. Fault scarp detection beneath dense vegetation cover: airborne Lidar mapping of the Seattle fault zone, Bainbridge Island, Washington State. In: Proceedings of the American Society of Photogrammetry and Remote Sensing Annual Conference, Washington, DC, May, 2000. Heuel, S., Forstner, W., Lang, F., 2000. Topological and geometrical reasoning in 3D grouping for reconstructing polyhedral surfaces. International Archives of Photogrammetry and Remote Sensing 33 (B3/1), 397–404. Irish, J.L., Lillycrop, J.W., 1999. Scanning laser mapping of the coastal zone: the SHOALS system. ISPRS Journal of Photogrammetry and Remote Sensing 54, 123–129. Kokkas, N., 2005. City modeling and building reconstruction with Socet Set v.5.2. BAE Systems, Customer presentation in 2005 GXP Regional User Conference, Cambridge, UK, September 19–21. Kokkas, N., Dowman, I., 2006. Fusion of airborne optical and LIDAR data for automated building reconstruction. In: Proceedings of the American Society for Photogrammetry and Remote Sensing, Reno Nevada, May 1–7. Mather, P., 1987. Computer Processing of Remotely-sensed Images. Wiley, New York, 353pp. Maune, D. (Ed.), 2001. The DEM users manual. American Society for Photogrammetry and Remote Sensing, Bethesda, 549pp. Miliaresis, G., 2001a. Geomorphometric mapping of Zagros ranges at regional scale. Computers and Geosciences 27, 775–786. Miliaresis, G., 2001b. Extraction of bajadas from digital elevation models and satellite imagery. Computers and Geosciences 27, 1157–1167. Miliaresis, G., 2006. Geomorphometric mapping of Asia minor from globe digital elevation model. Geografiska Annaler 88A, 209–221. Miliaresis, G., Argialas, D., 1999. Segmentation of physiographic features from the global digital elevation model/GTOPO30. Computers and Geosciences 25, 715–728. Miliaresis, G., Argialas, D., 2000. Extraction and delineation of alluvial fans from DEMs and Landsat TM Images. Photogrammetric Engineering and Remote Sensing 66, 1093–1101. Miliaresis, G., Argialas, D., 2002. Quantitative representation of mountain objects extracted from the GTOPO30 DEM. International Journal of Remote Sensing 23, 949–964. Miliaresis, G., Illiopoulou, P., 2004. Clustering of Zagros ranges from the Globe DEM representation. International Journal of Applied Earth Observation and GeoInformation 5, 17–28. Miliaresis, G., Kokkas, N., 2004. Segmentation and terrain modeling of extra-terrestrial chasmata. Journal of Spatial Sciences 49, 89–99.
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087 Miliaresis, G., Sabatakakis, N., Koukis, G., 2005. Terrain pattern recognition and spatial decision making for regional slope stability studies. Natural Resources Research 14, 91–100. Pitas, I., 1993. Digital Image Processing Algorithms. PrenticeHall, London, 362pp. Shaw, G., Wheeler, D., 1985. Statistical Techniques in Geographical Analysis. Wiley, Chichester, 364pp. Skidmore, A., 1989. A comparison of techniques for calculating gradient and aspect from a gridded digital elevation model. International Journal of Geographical Information Systems 3, 323–334.
1087
Turker, M., Cetinkaya, B., 2005. Automatic detection of earthquake-damaged buildings using DEMs created from pre- and port-earthquake stereo aerial photographs. International Journal of Remote Sensing 26, 823–832. Vosselman, G., Gorte, B.G.H., Sithole, G., Rabbani, T., 2004. Recognizing structure in laser scanner point clouds. International Archives of Photogrammetry Remote Sensing and Spatial Information Sciences 36 (8/W2), 33–38. Zevenbergen, L., Thorne, C., 1987. Quantitative analysis of land surface topography. Earth Surface Processes and Landforms 12, 47–56.
Computers & Geosciences 33 (2007) 1076–1087 www.elsevier.com/locate/cageo
Segmentation and object-based classification for the extraction of the building class from LIDAR DEMs$ George Miliaresisa,, Nikolaos Kokkasb a
Department of Geology, University of Patras, Rion 26504, Greece Department of Geomatic Engineering, University College London, Gower Street, London WC1E 6BT, UK
b
Received 10 August 2006; received in revised form 22 October 2006; accepted 28 November 2006
Abstract A new method is presented for the extraction of a class for buildings from light detection and ranging (LIDAR) digital elevation models (DEMs) on the basis of geomorphometric segmentation principles. First, seed cells and region growing criteria are specified. Then an object partition framework is defined on the basis of region growing segmentation. Size filtering is applied to objects and connected components labelling identifies background and foreground objects that are parametrically represented on the basis of elevation and slope. K-means classification reveals a set of clusters. The interpretation of the spatial distribution of clusters assisted by the interpretation of cluster centroids, allows for the identification of the building class, as well as building sub-classes with different geomorphometric characteristics. r 2007 Elsevier Ltd. All rights reserved. Keywords: Digital elevation model; LIDAR; Region growing; Urban segmentation; Terrain modelling; Terrain classification
1. Introduction High-resolution (spacing is less than 2 m) digital elevation models (DEMs) derived by airborne laser altimetry (LIDAR: light detection and ranging) are available for selected areas.1,2 The development of airborne laser scanning goes back to the 1970s (Irish $ Code available from server at http://www.iamg.org/ CGEditor/index.htm. Corresponding author. Tel.: +30 2610 996296; fax: +30 2610 991 900. E-mail addresses: [email protected] (G. Miliaresis), [email protected] (N. Kokkas). 1 Bare Earth Digital Elevation Model Generated from LIDAR Data. University of IDAHO, http://inside.uidaho.edu/ default.htm. 2 Maryland Department of Natural Resources, http://dnrweb. dnr.state.md.us/gis/data/lidar/.
0098-3004/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2006.11.012
and Lillycrop, 1999). By emitting a laser pulse and precisely measuring the return time to the source, the range can be calculated using the value for the speed of light.3 Elevations can be derived rapidly and at high resolution from LIDAR in comparison to manual reconstruction from photogrammetric techniques (time consuming, not a cost-effective solution). The LIDAR DEMs are extremely valuable in the field of earth sciences since elevation data from different flights might be compared to determine the patterns and magnitudes of coastal change4 3
LIDAR Tutorial, NASA, http://wwwghcc.msfc.nasa.gov/ sparcle/sparcle_tutorial.html. 4 Remote Sensing for Coastal Management, NOAA, http:// ekman.csc.noaa.gov/TCM/.
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
(erosion, over wash, etc.) and the loss of buildings and infrastructure.5 LIDAR DEMs are used in tectonic studies for fault recognition (Harding and Berghoff, 2000) and topographic change mapping.6 Additionally, airborne laser altimetry data were also used in an attempt to examine and model the surface morphology of landslides (Glenn et al., 2006). The building class extraction from high-resolution urban DEMs (spacing o2 m) is of primary importance in many applications, including urban planning, telecommunication network planning and vehicle navigation which are of increasing importance in urban areas (Kokkas, 2005). Nevertheless, the building class mapping is of great importance in the earth science field, for earthquake damage assessment. For example the difference DEM with 5 m spacing derived from pre- and post-earthquake aerial images revealed successfully the collapsed buildings, caused by the 1999 Izmit earthquake in Turkey (Turker and Cetinkaya, 2005). The critical element of detecting earthquake-induced heights changes through the DEM differencing method, was deciding where to place the boundaries between change and no-change pixels (Turker and Cetinkaya, 2005). Thus, if a building class was detected from the pre-earthquake DEM then the elevation differences caused by city ruins over the nonbuilding class (roads, city parks) could be omitted. Methods for building detection and reconstruction from laser altimeter data are applied to the point cloud (irregularly spaced points) and either directly derive the surface parameters in a parameter space by clustering the point cloud or segment a point cloud based on criteria like proximity of points or similarity of locally estimated surfaces (Vosselman et al., 2004; Heuel et al., 2000). Kokkas and Dowman (2006) introduced a method that employs a semi-automated technique for generating the building hypothesis by fusing LIDAR data with stereo matched points extracted from the air– photograph stereo model. The roof reconstruction is achieved by implementing a least squares-plane fitting algorithm on the LIDAR point cloud and subsequently neighbouring planes are merged using Boolean operations for the generation of solid features. 5
Hurricane and Extreme Storm Impact Studies: Coastal and Nearshore Mapping with Scanning Airborne Laser, US Geological Survey. 6 Active Tectonics Research Group, Arizona State University, http://lidar.asu.edu/data.html.
1077
On the other hand, quantitative techniques have been developed in order to automate the interpretation of terrain features from DEMs and various geomorphometric parameters were developed in an attempt to characterize the landscape (Miliaresis and Kokkas, 2004). Geomorphometric segmentation extracted mountains developed on different base levels (Miliaresis, 2001a) and fluvial landforms developed along mountain fronts (Miliaresis, 2001b) from moderate resolution DEMs. The technique is summarized as a region growing segmentation algorithm that uses seeds (e.g. the ridge cells) and growing criteria (usually based on slope and elevation). Whereas geomorphic segmentation defines both abstractions of landforms and an object partition framework of the landscape, the parametric representation of landscape objects on the basis of their spatial three-dimensional arrangement (Miliaresis and Argialas, 2002) allows the recognition of landscape objects) that have similar and contrasting ranges of characteristics, leading to a terrain classification scheme (Miliaresis and Illiopoulou, 2004) according to geomorphologic principles and understanding (Miliaresis et al., 2005). The aim of the current research effort is to design a new method for extraction of the building class from LIDAR DEMs on the basis of the geomorphometric segmentation principles. Thus, an object partition framework of the LIDAR DEM will be defined. Finally, objects representation on the basis of geomorphometric parameters combined by object classification is expected to allow the extraction and mapping of the building class. 2. Methodology First, the study area and its hypsometric characteristics are introduced. Then seeds cells and region growing criteria are defined. Finally, region growing segmentation is performed and thus, an object partition framework is defined. Connected components labelling identifies background and foreground objects that are parametrically represented on the basis of elevation and slope and classified. 2.1. Study area and data The LIDAR data used in this paper data were acquired by NPA Group using their in-house build LIDAR sensor during the leaf-on season. The
ARTICLE IN PRESS 1078
G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
Fig. 1. LIDAR DEM of study area and elevation frequency histogram. Elevation is in range 12– 215 m, the lighter a point is the greatest its elevation.
LIDAR point cloud represents the Bloomsbury area in central London bounded by co-ordinates (lower left: 529,174, 181,500, upper right: 530,199, 182,525) in the British National Grid. The point cloud has been post-processed and delivered with a 1 m 1 m (Fig. 1) sampling resolution covering an area of approximately 1 km2 (1024 rows and 1024 columns consist the derived regular grid DEM). The elevation frequency histogram (Fig. 1) indicates that the majority of DEM cells (98.94%) present elevation in the range 20–60 m. 2.2. Seed cells Median filtering was applied for seed pixels to be defined. The seed points were computed in a limited size kernel window. More specifically, the kernel size is 5 5 grid cells; note that each grid cell is 1 m apart. Due to the very small kernel size, it is impossible for the seed points to be affected from the regional topographic surface. Additionally, the regional topographic surface in the study area is a gently sloping plane. Even if the regional topographic surface was more complicated (either 2 or 3 order polynomial surface), the elevation differences computed in a 25 m2 area would be residual elevation anomalies and not regional elevation anomalies especially in an urban (city) environment. Median filtering removes very high and very low elevation values within the selected kernel size while
at the same time preserves edges (Mather, 1987). Then the Difference DEM (LIDAR DEMMedian filtered DEM) was derived. The Difference DEM is depicted in Fig. 2a. The Difference DEM revealed the cells of the LIDAR DEM where severe median filtering occurred. Note that elevation differences less than 1 m are depicted white while differences greater than 9 m are depicted black). The frequency of the elevation differences (Fig. 2b) indicates that the majority of DEM cells (97.1%) present difference in the range [4, 11]. The distribution is positively skewed (greater frequencies occur for positive elevation differences with respect to the corresponding negative elevation differences). Positive differences were interpreted to be associated to building edges. As it can be seen in the Figs. 2c and d, the greatest the elevation difference, the more the detected cells stands along the borders of the interpreted city buildings. Under a trial and error procedure though experimentation and visual interpretation of the results, was found that positive elevation differences in the range 5–9 m are more associated to the building edges (Fig. 2d). On the contrary differences greater than 9 m were interpreted to be associated to the internal building/roof structure (chimneys, statues, etc.). Thus, as seed cells were selected those presented elevation difference within the range of 5–9 m.
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
1079
Fig. 2. (a) The difference DEM (LIDAR DEM- Median filtered DEM). (b) Frequency of elevation differences. (c) Enlarged portion of difference DEM (differences less than 1 m are depicted white while differences greater than 9 m are depicted black). (d) Enlarged portion of difference DEM (differences less than 5 m are depicted white while differences greater than 9 m are depicted black).
2.3. Region growing segmentation For region growing segmentation, region growing criteria should be defined (Miliaresis and Argialas, 1999, 2000). In previous research efforts, aimed to landform segmentation (Miliaresis, 2001a b, 2006) the criteria were defined on the basis of slope and elevation. If the partial derivatives of elevation (H) along the east (x) and the north (y) direction are known then then slope and aspect (slope pointing direction) are computed from Burrough, 1987). sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 dH dH Slope ¼ , (1) þ dx dy Aspect ¼ arc tan
dH=dy . dH=dx
(2)
In order to estimate the partial derivates of elevation over a 3 3 kernel, Evans (1980) fitted a 6-parameter quadratic equation while Zevenbergen and Thorne (1987) used a 9-parameter quadratic equation. Evans (1980) method results to a polynomial surface which will not necessarily pass through the 9 grid cells included in the 3 by 3 kernel, while the opposite is true for Zevenbergen and Thorne (1987). A review and a comparison of these methods (Skidmore, 1989) proved that Evans (1980) method is more precise since Zevenbergen
and Thorne (1987) method is affected more from random noise and errors evident in DEMs. The slope frequency histogram (Fig. 3b) is completely different than the slope frequency distributions, we are used to observe in geomorphologic applications (the slope frequency is decreasing gradually for increasing slope values) where digital terrain models (DTMs) are used (Miliaresis and Kokkas, 2004; Miliaresis et al., 2005). DTM depicts ground elevation (the height of artificial structures like buildings is ignored) while in a DEM the elevation of the highest object (height of the trees and buildings) above the ground is included (Maune, 2001). So in DTMs (with spacing 30–1000 m), the ground elevation prevails and thus the slope frequency distribution is the outcome of geomorphologic (erosion, deposition) and geologic (tectonism) processes (Miliaresis, 2006). On the contrary in urban LIDAR DEMs, the building structure prevails and due to the high DEM resolution (p1 m), the slope frequency distribution we observe is the outcome of the interaction of the building surface (abrupt changes of elevation, high and abrupt changes of slope). The study of both the slope image (Fig. 3a) and the slope frequency histogram (Fig. 3b) combined by a trial and error procedure for selecting slope cut-off values, through experimentation, assisted by the photo-interpretation of the results, indicated
ARTICLE IN PRESS 1080
G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
Fig. 3. Slope image (the higher the slope value, the darkest a cell is). Frequency of slope in the interval 1– 901, (classes for slope o11 were excluded due to its great occurrence).
that the region growing criteria are defined for slope values in the interval [401, 901]. Then, an iterative region growing segmentation algorithm (Miliaresis, 2001a) was implemented. The seeds (Fig. 4) formed the initial set of building cells. If a cell of the DEM with slope in the interval [401, 901] was an 8-connected neighbour to the current set of cells, it was flagged as a new building cell and the current set of building cells was updated. Note that the 8-connected neighbours are the 8 pixels surrounding the central pixel in a 3 3 kernel. The segmentation stopped if no more cells were added during the current iteration. Finally, 448,865 cells were identified after 29 iterations (Fig. 4a). A connected component-labeling algorithm (Pitas, 1993) was applied and regions formed by 8connected background cells (black pixels) were identified and labelled. Many objects (the most of them very small in size) are evident in the segmented image. The very small-sized objects correspond to residual objects evident either on the roofs or in the non-building class. One could ask ‘‘what is the minimum size for an object to be clasified as a residual one?’’ an answer could be that building blocks in the study area should be far greater in size
than the ‘‘residual objects’’! In this particular case study (LIDAR DEM portion of the central London), we assume that the building blocks should be far greater in size than 120 m2 (120 cells). Small in size foreground (white) and background (black) objects were eliminated by applying object size filtering. More specifically, foreground objects with size less than 120 cells were merged to the background class (Fig. 4b). Six thousand six hundred and thirty-three cells (1.48% of the foreground cells class) cells were merged to the background and so the foreground cell class is formed by 441,828 cells. Then, 114,009 (18.8%) of the background class cells (black cells) were merged to the foreground (Fig. 4c) for an object size threshold of 120 cells (thus the foreground cell class is formed by 555,837 cells). Finally, the remaining foreground objects (Fig. 5a) and background objects (Fig. 5b) were labelled by the connected components algorithm (76 foreground and 291 background objects were identified) and merged to a single image (Fig. 5c). Region growing segmentation, connects components labelling, size filtering and the subsequent object parametric representation were performed with GeoLogic Shell software (Miliaresis, 2001b)
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
1081
Fig. 4. Segmentation and object size filtering. (a) Cells segmented (white pixels) by region growing segmentation algorithm (left) and 448,461 (42.77%) cells belong to foreground objects class. (b) Foreground objects (formed by white pixels) filtering with size less than 120 cells (middle image). (c) Background objects (formed by black cells) filtering with size less than 120 cells (right image).
Fig. 5. Foreground and background objects after size filtering. (a) Seventy-six foreground objects in left image (the darker an object is the greater the ID value is), (b) 291 background objects in middle image (the darker an object is the greater the ID value is), (c) merged image that contains background and foreground objects (right image). IDs 1–76 correspond to foreground objects while IDs 77–367 corresponds to background objects (the lighter an object is the greater the ID value is).
that is freely available for download on the internet from the WEB site of the International Association for Mathematical Geology.7 7 Computers and Geosciences, software code section, http:// www.iamg.org/CGEditor/cg2001.htm. Geologic Shell software that is part of Miliaresis (2001b) paper.
2.4. Parametric representation The objects in Fig. 5c are interpreted to belong either to the building class or to the non-building class (roads, vegetated areas, components of a building, etc.). Additionally, objects that are adjacent or the one encloses the other could both belong
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
1082
to the building class. It is assumed that if each object is parametrically represented by a set of geomorphometric attributes then unsupervised classification could define clusters of objects that belong either to the building or the non-building class. Four attributes were used for the parametric representation of objects (Fig. 5c):
H, mean elevation that was computed as the average elevation of the cells that belong to an object’s region. R, roughness, the standard deviation of elevation of the cells forming the object (a measure of the vertical variability of object’s height). S, object’s mean slope. SR, standard deviation of slope (a measure of the variability of object’s slope).
2.5. Classification The crosscorrelation between attributes is given in Table 1. It is observed that a partial correlation exists between S and SR, as well as R and SR. In other words, if the slope is high that for the most of the objects standard deviation of slope and elevation should be high, but this is not always true. Cluster analysis is a multivariate procedure, which is commonly used for regional classification. It is based on some measurement of distance among objects (for example Euclidean distance), which is calculated in a c-dimensional space, where c represents the number of attributes used in the clustering process (Mather, 1987). The centroid method was employed that requires a priori definition of the number of clusters. The main problem in cluster analysis has to do with the different value range of the different attributes (for example elevation and slope) included in the distance calculation. In order to solve this problem, the data were standardized (transformed to a new dataset in which the corresponding attribute mean value is
Table 1 Attribute correlation matrix H H R S SR
R
S
SR
1 0.56 0.63
1 0.72
1
1 0.28 0.12 0.03
zero and the corresponding attribute std. dev. is 1) before the exploratory (K-means) clustering algorithm was applied (Miliaresis and Argialas, 2002). Thus, the data values were transformed by subtracting the corresponding attribute mean value and dividing the outcome by the corresponding attribute standard deviation (Clark and Hosking, 1986). Ten clusters of objects were considered through experimentation and a trial and error procedure. Each object was assigned to a cluster on the basis of its Euclidean distance from the cluster centroid. Their cluster gravity centres are given in Table 2 while the graphical representation of the gravity centres was computed on the basis of the standardized co-ordinates (Fig. 6). The numbers of objects that belong to each cluster are also given in Table 2. The similarity among cluster centroids could be assessed either by the graphical representations of centroids (Fig. 6) or by divergence (Table 3). The cluster compactness (Table 3) was computed on the basis of compactness ¼ 100
mean st: dev: , Mean
(3)
where mean and std. dev. are the mean distance and the corresponding standard deviation of the objects from the gravity centre of the cluster. In order to evaluate the existence of statistically significant differences in the distances among the cluster centres the statistical technique of one-way analysis of variance (ANOVA) was used. ANOVA calculates the between-cluster mean square and the within-cluster mean square (Miliaresis and Illiopoulou, 2004). The ratio of these two mean squares is Table 2 Final cluster (standardized) centroids and the number of objects in each cluster (sixth row) Cluster ID
1 2 3 4 5 6 7 8 9 10
Standardized centroid co-ordinates H
R
S
SR
Number of objects
3.2498 1.1694 0.1284 4.2958 0.324 0.4666 0.6108 4.9048 6.3596 0.9844
8.5073 0.3474 0.3539 9.4762 0.932 2.5654 0.8661 0.1302 8.7856 0.3556
0.4618 0.5844 0.37 1.9155 1.638 2.2769 2.3381 0.3923 0.1835 0.5615
2.2419 0.3778 0.5094 2.3725 1.9751 1.3411 0.0694 0.025 2.2435 0.4224
1 81 119 1 56 2 19 2 1 85
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
the usual ANOVA F statistic F¼
variation between groups sum of squares=ðk 1Þ , variation within groups sum of squares=ðN kÞ
(4) where k is the number of groups and N is the total sample size. The logic of the test argues that for the class (cluster) means to be significantly different, interclass (between cluster) variance must be of much greater magnitude than intraclass (within clusters) variance (Clark and Hosking, 1986; Shaw and Wheeler, 1985). The results of the ANOVA procedure for the each centroid are shown in Table 4. Although the significance levels for the F statistic are below 0.001, the values of F in the case of cluster analysis should be mostly interpreted in terms of comparing
1083
differentiation of the cluster means for individual variables. It is concluded that clusters are grouped to two categories (a) clusters 2, 3 and 10, (b) clusters 5 and 7 on the basis of similarity/divergence (Table 3, Fig. 6). In the first category, the clusters 2, 3, and 10 are differentiated on the basis of mean elevation while the second category is formed by objects with almost equal mean elevation. Clusters 7 and 10 are less compact (Table 3) meaning that they could be divided possibly to sub-clusters. The spatial distribution of clusters is depicted in Fig. 7 in an attempt to assist the interpretation.
2.6. Discussion of the results Fig. 8 presents the nature of the LIDAR elevation data, which is quite noisy while at the same time residual elevation objects are observed over the roofs. These objects correspond either to arificial objects induced by noise effects and errors evident in the LIDAR DEM or to small elevated features that are difficult to recognise and assigned to a certain feature type, due to both the DEM resolution and the local elevation errors that are present.
Table 4 Analysis of variance
Fig. 6. Centroids (standardized) for clusters 2, 3, 5, 7 and 10 that include majority of objects (Table 2).
Variable Cluster MSa DFa Error MSa DFa F
Prob.
H R S SR
0 0 0 0
a
36.14 39.03 36.79 34.13
9 9 9 9
0.114 0.041 0.097 0.164
357 357 357 357
316.83 943.36 376.36 206.94
MS, mean square; DF, degrees of freedom.
Table 3 Distances between cluster centroids (divergence) and cluster compactness for clusters 2, 3, 5, 7, and 10 Distance between cluster centroids
Cluster
2
3
2 3 5 7 10
0 1.32 3.58 3.24 2.15
0 3.47 3.12 0.88
Statistics for the objects distance from the corresponding cluster centroid 5
0 2.05 3.74
7
0 3.56
10
Cluster
Mean
St. dev
Compactness %
0
2 3 5 7 10
0.570 0.444 0.654 0.634 0.614
0.263 0.217 0.341 0.381 0.383
53.8 51.2 47.8 39.9 37.6
ARTICLE IN PRESS 1084
G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
Fig. 7. Spatial distribution of 10 clusters with labels 1 to 10 is given in top left image (the greater the cluster ID the darkest an object is). In subsequent images, clusters 2, 3, 5, 7 and 10 (each one is depicted white) are presented in an attempt to assist the interpretation.
Fig. 8. LIDAR DEM visualisation. (a) Direction of line of sight (NW-SE) that is determined by position of observer (NW) and position of target (SE). Line of sight is superimposed over LIDAR DEM. (b) LIDAR DEM visualisation along line of sight.
Thus, by the size filtering implemented (Fig. 4), the building blocks were generalised with respect to roof appearance.
The visual interpretation (Fig. 7) indicates that the non-building class is formed by cluster 2 while the urban vegetation class corresponds to cluster 7.
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
The interpretation of cluster centroids (Fig. 6, Table 2) indicate that the non-building class is formed by the less elevated, most flat, less variable with respect to slope and elevation objects. The urban vegetation class is formed by objects with the highest slope and with significant slope variability (standardized centroid is 40). On the other hand, the mean standardized elevation is below 0 while roughness is high (1). Thus, the vegetation class is formed by objects with low mean elevation (higher than the elevation of the non-building class), while slope is maximized and the standard deviation of slope and elevation is high. The last remarks are functions of tree canopy properties and thus different vegetation types might present different parametric representation. The building class is formed mainly by the objects of cluster 5 (Fig. 7). For the objects of the building
1085
class slope variability is maximized (Fig. 6) while slope is among the maximum ones (Table 2). The interpretation given is that it happens due to the building edges that are in contact with the nonbuilding class (Fig. 7). Mean elevation (H) and roughness (R) are higher than the corresponding values of both building and vegetation classes. The clusters 10 and 3 correspond mainly to internal portions of buildings. That is why the roughness, slope and slope variability is among the minimum while mean elevation is maximized (Fig. 6). The outliers (clusters 1, 4, 6, 8 and 9) are interpreted on the basis of cluster centroids (Table 2). They certainly belong to the building class due to the extremely high mean elevation values presented and they are interpreted to correspond to elevated objects evident mainly in building roofs (chimneys, etc.).
Fig. 9. Visualisation of segmented building and vegetation classes. (a) Segmented building classes and vegetation classes near a city park as wells as direction of line of sight (S-N). (b) Visualisation of segmented classes (in city park) that were superimposed over LIDAR DEM. (c) Visualisation of LIDAR DEM in city park. (d) A portion of Fig. 5c that contains 367 segmented objects as wells as direction of line of sight (SE-NW). (e) Visualisation of the segmented objects that were superimposed over the LIDAR DEM.
ARTICLE IN PRESS 1086
G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087
In Fig. 9a selected segmented building and vegetation objects are presented. The visualisation of these classes (Fig. 9b) as well as the visualisation of LIDAR DEM (Fig. 9c) indicate the two portions to which the most of the buildings were segmented to (the internal portion that is depict with the white tone and the building border that is depicted with a shade of gray) as well as the segmentation of the city park vegetation. These findings are also supported by the visualisation of the Fig. 5c of the segmented objects (before the objects merging according to their parametric representation to 10 clusters/building classes) in Figs. 9d and c. 3. Conclusion Geomorphometric region growing segmentation combined by median filtering for identifying seed cells, connected components labelling, size filtering and object labelling, object parametric representation on the basis of slope and elevation attributes and classification was proved capable of delineating the building class within the study area. The interpretation of cluster centroids allowed the identification building sub-classes with different geomorphometric characteristics that are associated with different building portions and present different hazard assessment. Difference DEM derived from pre- and postearthquake LIDAR DEMs can be used in order to reveal the earthquake collapsed buildings, if the building class is detected from the pre-earthquake DEM according to the method described. Nevertheless, the proposed method for building class detection has a major disadvantage. It requires a certain level of user interaction for some crucial parameters, which are very difficult to define automatically in varied situations and further research is required for full automation. Acknowledgements The authors are grateful for, and this paper was greatly benefited from, the thorough and evaluations of the two anonymous reviewers. The authors would like to thank NPA Group for providing the LIDAR data and express their gratitude for the financial support offered by the Greek Scholarship Foundation. References Burrough, P.A., 1987. Principles of Geographical Information Systems for Land Resources Assessment. Oxford University Press, Oxford, 194pp.
Clark, W., Hosking, P., 1986. Statistical Methods for Geographers. Wiley, New York, 518pp. Evans, I., 1980. An integrated system for terrain analysis and slope mapping. Zeitschrift fuer Geomorphologie N.F. Suppl. -Bd. 36, 274–290. Glenn, N., Streutker, D., Chadwick, J., Thackray, D.G., Dorsch, S., 2006. Analysis of LiDAR-derived topographic information for characterizing and differentiating landslide morphology and activity. Geomorphology 73, 131–148. Harding, D.J., Berghoff, G.S., 2000. Fault scarp detection beneath dense vegetation cover: airborne Lidar mapping of the Seattle fault zone, Bainbridge Island, Washington State. In: Proceedings of the American Society of Photogrammetry and Remote Sensing Annual Conference, Washington, DC, May, 2000. Heuel, S., Forstner, W., Lang, F., 2000. Topological and geometrical reasoning in 3D grouping for reconstructing polyhedral surfaces. International Archives of Photogrammetry and Remote Sensing 33 (B3/1), 397–404. Irish, J.L., Lillycrop, J.W., 1999. Scanning laser mapping of the coastal zone: the SHOALS system. ISPRS Journal of Photogrammetry and Remote Sensing 54, 123–129. Kokkas, N., 2005. City modeling and building reconstruction with Socet Set v.5.2. BAE Systems, Customer presentation in 2005 GXP Regional User Conference, Cambridge, UK, September 19–21. Kokkas, N., Dowman, I., 2006. Fusion of airborne optical and LIDAR data for automated building reconstruction. In: Proceedings of the American Society for Photogrammetry and Remote Sensing, Reno Nevada, May 1–7. Mather, P., 1987. Computer Processing of Remotely-sensed Images. Wiley, New York, 353pp. Maune, D. (Ed.), 2001. The DEM users manual. American Society for Photogrammetry and Remote Sensing, Bethesda, 549pp. Miliaresis, G., 2001a. Geomorphometric mapping of Zagros ranges at regional scale. Computers and Geosciences 27, 775–786. Miliaresis, G., 2001b. Extraction of bajadas from digital elevation models and satellite imagery. Computers and Geosciences 27, 1157–1167. Miliaresis, G., 2006. Geomorphometric mapping of Asia minor from globe digital elevation model. Geografiska Annaler 88A, 209–221. Miliaresis, G., Argialas, D., 1999. Segmentation of physiographic features from the global digital elevation model/GTOPO30. Computers and Geosciences 25, 715–728. Miliaresis, G., Argialas, D., 2000. Extraction and delineation of alluvial fans from DEMs and Landsat TM Images. Photogrammetric Engineering and Remote Sensing 66, 1093–1101. Miliaresis, G., Argialas, D., 2002. Quantitative representation of mountain objects extracted from the GTOPO30 DEM. International Journal of Remote Sensing 23, 949–964. Miliaresis, G., Illiopoulou, P., 2004. Clustering of Zagros ranges from the Globe DEM representation. International Journal of Applied Earth Observation and GeoInformation 5, 17–28. Miliaresis, G., Kokkas, N., 2004. Segmentation and terrain modeling of extra-terrestrial chasmata. Journal of Spatial Sciences 49, 89–99.
ARTICLE IN PRESS G. Miliaresis, N. Kokkas / Computers & Geosciences 33 (2007) 1076–1087 Miliaresis, G., Sabatakakis, N., Koukis, G., 2005. Terrain pattern recognition and spatial decision making for regional slope stability studies. Natural Resources Research 14, 91–100. Pitas, I., 1993. Digital Image Processing Algorithms. PrenticeHall, London, 362pp. Shaw, G., Wheeler, D., 1985. Statistical Techniques in Geographical Analysis. Wiley, Chichester, 364pp. Skidmore, A., 1989. A comparison of techniques for calculating gradient and aspect from a gridded digital elevation model. International Journal of Geographical Information Systems 3, 323–334.
1087
Turker, M., Cetinkaya, B., 2005. Automatic detection of earthquake-damaged buildings using DEMs created from pre- and port-earthquake stereo aerial photographs. International Journal of Remote Sensing 26, 823–832. Vosselman, G., Gorte, B.G.H., Sithole, G., Rabbani, T., 2004. Recognizing structure in laser scanner point clouds. International Archives of Photogrammetry Remote Sensing and Spatial Information Sciences 36 (8/W2), 33–38. Zevenbergen, L., Thorne, C., 1987. Quantitative analysis of land surface topography. Earth Surface Processes and Landforms 12, 47–56.