A fuzzy decision making system for building damage map creation using high resolution satellite imagery

Recent studies have shown high resolution satellite imagery to be a powerful data source for post-earthquake damage assessment of buildings. Manual interpretation of these images, while being a reliable method for finding damaged buildings, is a subjective and time-consuming endeavor, rendering it unviable at times of emergency. The present research, proposes a new state-of-the-art method for automatic damage assessment of buildings using high resolution satellite imagery. In this method, at the first step a set of preprocessing algorithms are performed on the images. Then, extracting a candidate building from both preand post-event images, the intact roof part after an earthquake is found. Afterwards, by considering the shape and other structural properties of this roof part with its pre-event condition in a fuzzy inference system, the rate of damage for each candidate building is estimated. The results obtained from evaluation of this algorithm using QuickBird images of the December 2003 Bam, Iran, earthquake prove the ability of this method for post-earthquake damage assessment of buildings.


Introduction
Natural disasters such as earthquakes and floods have the consequence to affect many people, not only through the destruction they cause, but also through homelessness, injury and even death (Havidán, 2006).Earthquakes may be known as the most prevalent natural hazard and many with different magnitudes, which have caused various losses, have been recorded throughout history (USGS, visited 2011).Based on the recorded earthquakes on USGS website, more than 800 000 fatalities have been reported during the last decade alone.Even though earthquakes are not predictable by current technology and therefore no short-term preparedness is possible, any rescue activities that are performed quickly after an earthquake can decrease the number of fatalities.Rapid, accurate and comprehensive knowledge about the damaged area can, therefore, be very helpful during the response phase of disaster management.One of the most crucial pieces of information that can be used in the response phase after a natural disaster is a building damage map that shows the extent of damage for every individual building or, on a larger scale, for every district in an urban area.
There are many data sources such as satellite and aerial images, ground observation and LiDAR that can provide useful information for damage map generation (Li et al., 2008;Rezaeian, 2010).Among all available sources, images are more comprehensive and rapid-access for providing information about the damaged area.During the last decade, much of the research has focused on using this data source for post-earthquake damage assessment, which has led to image-based damage assessment trending amongst the hottest topics in photogrammetry and remote sensing (Chini et al., 2009;Matsuoka and Yamazaki, 2004;Rezaeian, 2010;Thomas, 2010;Turker and Sumer, 2008).Consequently, different methods and techniques have been reported by researchers, which can be classified based on various criteria.
The first criterion to categorize image-based damage assessment methods is based on the type of input data, which can be either airborne or space borne.In comparison with airborne data, rapid access and continuous coverage of satellite images have allowed most researchers to apply these images for damage assessment (Brunner et al., 2010;Vu and Ban, 2010;Yamazaki and Matsuoka, 2007); while other researchers have used aerial images for damage assessment (Li et al., 2008;Rezaeian, 2010;Thomas, 2010;Turker and Sumer, 2008).
One can also categorize damage assessment techniques based on the type of interpretation, which can be done visually or automatically.In visual interpretation, a human operator conducts images interpretation.Therefore, proficiency in working with the image (air/space-borne) increases the reliability and trueness of the generated damage map (Rezaeian, 2010).Although visual interpretation is a reliable tool, it is subjective and time-consuming (Ogava and Yamazaki, 2000), rendering it less useful at times of emergency.To eliminate these drawbacks, automatic interpretation of damage is introduced, in which interpretation and analyses based on pre-knowledge information are done by a computer.In automatic interpretation, damaged buildings are recognized by using different clues which are automatically extracted from the image(s).For example, information from adjacent pixels, such as edges and texture, were used by Turker and San (2003).
Alternatively, input images for damage assessment may be acquired by active sensors.One of the main advantages of radar images is that they can be used regardless of sunlight and weather conditions, and in cases of poor weather conditions, they may be the only available data (Matsuoka and Yamazaki, 2004;Rezaeian, 2010).However, in comparison with optical images, they are not easy to visually interpret.In other words, the methods that work based on radar images cannot precisely evaluate the rate of damage for each individual building (Chini et al., 2009;Dong et al., 2011;Matsuoka and Yamazaki, 2004;Yamazaki and Matsuoka, 2007).In Chini et al. (2009), SAR and QuickBird images of the city of Bam, Iran, were separately applied for damage assessment and the results were compared.In that study, the datasets included images acquired before and after the December 2003 earthquake.And, in contrast to SAR images, QuickBird images provided more accurate results for a single building.Furthermore, SAR images, high spatial resolution optical images and vector map were applied together for mapping earthquake damage at the block scale (Stramondo et al., 2006).
Even though a pre-event image by itself does not convey any information about the damage area, through a comparison with the post-event image, the interpreter (human or computer) is able to make a better decision about a building's condition.In the majority of studies, both pre-and postevent images have been used in damage assessment (Rejaie and Shinozuka, 2004).In other studies, however, only a postevent image has been applied (Dell'Acqua and Polli, 2011;Kohiyama and Yamazaki, 2005).In these studies, damaged area is recognized using features which allow for easy differentiation between damaged and intact roofs.For example, Turker and Sumer (2008) applied watershed segmentation technique to separate the damaged area from intact building using aerial images.The method was evaluated using the aerial images of Golcuk, one of the urban areas strongly hit by the 1999 Izmit, Turkey, earthquake.
Ancillary data such as vector maps can provide beneficial information in image-based damage assessment.Using a vector map, the roof prints of the candidate buildings can be effectively identified in the image and consequently can improve the efficiency, accuracy and performance of the damage assessment process (Chesnel et al., 2008).Fortunately, vector map data of almost all cities worldwide are available.However, the different accuracy parameters of objects in a vector map, such as positional and geometrical, are very important.In addition, due to the usual changes of urban areas, the vector map must be a recently-updated version.In Dong et al. (2011) and Samadzadegan and Rastiveis (2008), a vector map has been used as auxiliary data along with SAR and QuickBird images, respectively, for damage assessment of the 2008 Wenchuan, China, and 2003 Bam, Iran, earthquakes.
Regardless of whether or not ancillary data is used, image(s) interpretation may be done at pixel-or object level (Gusella et al., 2005).In pixel-level approaches, each pixel is examined as an individual object and is labeled a separate damage state based on its characteristics.On the other hand, in object-based approaches, images are firstly segmented into meaningful regions, which are called image objects, and all further analyses are performed on these image objects.In such s case, the rate of damage is separately assigned for each individual image object.In Kouchi and Yamazaki (2005), both pixel-and object-level damage assessment of the 2003 Boumerdes, Algeria, earthquake were performed, achieving more promising results with object-level assessment.The same result was gained by Matsumoto et al. (2006) in damage assessment of the 2006 Central Java, Indonesia, earthquake.Although object-level approaches provide better results in comparison to pixel-level, setting the appropriate parameters for generating proper image objects is one of the main challenges of these methods.
Change detection using stereo images and therefore height information has been developed (e.g.Chaabouni-Chouayakh and Reinartz, 2011), but has not been applied for earthquake monitoring.Also, up to now the availability of stereo data shortly after an earthquake is very seldom given.
From the aforementioned studies, it can be concluded that the variety of input images has resulted in numerous techniques for automatic damage assessment of earthquakes.However, due to the wide range of uncertainty in recognizing and classifying damaged buildings, these techniques are up to now not as accurate as visual interpretation of images.In other words, they cannot handle this uncertainty as perfectly as an expert who uses a lot of knowledge during decision making.Despite the seeming reliability of visual interpretation, it is not very useful for fast generation of results as in times of emergency, and the necessity for an accurate automatic method is incontrovertible.In this research, therefore, we propose a new method based on a fuzzy inference system for automatic damage assessment using high resolution satellite imagery.In the following sections, the paper describes the details of the proposed method and presents results obtained through its implementation.

Proposed method
The proposed method in this paper is based on the flowchart shown in Fig. 1.Both pre-and post-event high resolution satellite imagery of a damaged area are required for generating a building damage map.In addition, a relevant vector map of the area, as ancillary data, is also needed to find the location of the buildings.As seen in Fig. 1, the initial step in implementing this method involves the pre-processing of the satellite images.Selecting a candidate building, building areas on both pre-and post-event images are extracted using geo-referencing information.After detecting the initial roof within pre-event building areas, the intact roofs after the earthquake are detected in post-event building areas.The roof detection step is done through a segmentation approach based on extracted textural information of the building areas.After that, with the aim of removing noisy pixels, the detected roof data are modified using morphological operators.Finally, by considering the shape and other structural properties of the modified roof data in a fuzzy inference system, the rate of damage for the candidate building is estimated.The algorithm is executed on each individual building in the damaged area to create a final damage map.Further details of the proposed method are described in the following sections.

Pre-processing
Due to temporal resolution, the pre-and post-event images from a similar region have different illumination conditions.Therefore, a set of pre-processing algorithms should be performed on the images.For this purpose, first, atmospheric and solar illumination effects are eliminated through atmospheric correction of the images.Then, orthorectification of the images are done to compensate the ground elevation.Pan-sharpened pre-and post-event images are also created by fusing the MS (multispectral red, green and blue channels) and Pan QuickBird images.Next, histogram equalization and histogram matching are needed to increase the spectral similarity of the images.Finally, the images are accurately geo-referenced with digital vector map using enough control points.

Building areas extraction
Locating a candidate building on the images is an important step in the proposed method, and is easily achieved through the use of geo-referencing information.Using the geo-referencing information, each point of the map can be located on the images and vice versa.Therefore, in the first step of the algorithm, roof corners of a candidate building are found on both pre-and post-event images and consequently the building areas are extracted.Increased accuracy of georeferencing information leads to more reliable results.Because accurate geo-referencing of the images is of paramount importance prior to running the algorithm, the registration errors of the input images should be performed in sub-pixel level of accuracy.The reader is referred to Richards and Jia (2006) and Schowengerdt (2007) for a comprehensive overview of existing methods for geo-referencing of satellite images.

Roof detection
In visually interpreting images for damage assessment, an expert usually evaluates the damage rate of a candidate building based on the shape and the structure of the remaining roof on the post-event image.Comparing a roof to its former condition on the pre-event image renders the interpretation more explicit and consequently results in more reliable decisions.The main idea of the proposed method presented here is to automatically estimate the damage rate for all individual buildings based on the extracted pre-and post-event roofs.Therefore, the next step in the implementation of this method is to extract the roof areas.Figures 2a and b depict a single building and a city block, respectively, on an aerial photo.As can be seen from the figures, several objects such as chimney and cooler may exist on a roof area.After an earthquake, based on the magnitude of earthquake and condition of a building, some parts of the roof might be destroyed.The question here is how to automatically distinguish the roof pixels in building areas?In the proposed method, the extraction of the roof parts is performed in four consecutive processes of texture analysis, segmentation, roof detection and modification.For this purpose, first, through different texture analysis techniques, useful features are extracted to make a better distinction between the roof and non-roof pixels.Then, using a segmentation technique, the building area is stratified into a number of segments.Next, the roof segments are recognized, and finally, for a better understanding of the extracted roofs, they are modified using morphological operators.These processes are described in greater detail in the forthcoming sections.

Texture analysis
In this method, a segmentation technique is applied to stratifying the building areas into a number of segments.However, distinguishing roof pixels from other pixels in a building area may not be accurately possible merely based on spectral information (red, green and blue channels) of the pixels.Therefore, other information is needed to be used along with spectral information for extracting the roof areas.In this case, textural information is known as a powerful tool in image analysis.A wide range of studies have applied textural information in image-based damage assessment (Rezaeian, 2010;Vu and Ban, 2010).
In the proposed method, various textures are implemented and superior texture features are selected manually based on visual interpretation of the texture images.In other words, a user, by observing a resulted texture feature and comparing to the original images, decides to select or omit the texture.A superior feature here is one that better distinguishes between roof and non-roof pixels.Once the feature selection is performed, the selected list can be applied to the entire area.To improve efficacy, one can perform an optimal feature selection algorithm for finding the best features, which requires an accurate ground truth.
In this study, 21 features in four categories of statistical, Haralick, Gabor and semi-variogram textural features were implemented (see Table 1).From this list, based on visual observation six features were selected as textural features to be used along with spectral features in the segmentation process (features are formatted in bold in the table).The details of implementing these features are described below.

First-order statistical features
First-order textural features such as mean and variance are statistics that are calculated from image values and do not consider pixel neighborhood relationships.In this research, mean gray-level of adjacent pixels in a region is considered as a first-order statistical descriptor to be applied in the segmentation process.This feature can be calculated using Eq.(1).

Haralick features
Haralick features, which are the well-known and widely used texture features in image analysis, were proposed by Haralick in Haralick et al. (1973).These descriptors consider the pixel neighborhood relationships and are known as the second-order statistical features.The basis of the Haralick features is a two-dimensional co-occurrence matrix.This matrix, P, is a n × n matrix, where n is the number of graylevels within an image.The matrix acts as an accumulator so that P[i, j ] counts the number of pixel pairs having the intensities i and j .Pixel pairs are defined by a distance and direction that can be represented by a displacement vector d = (dx, dy), where dx represents the number of pixels moved along the x-axis, and dy represents the number of pixels moved along the y-axis of an image slice.Many features can be derived from the co-occurrence matrix, such as entropy, homogeneity, sum mean, variance, correlation, maximum probability, etc.Among them, "cluster tendency" and "sum mean" can make a meaningful distinction between intact roof and damaged pixels.These two features, which can be computed using Eqs.( 2)-( 3), are therefore selected for clustering.
Cluster Tendency

Gabor features
For a given image I (x, y), its discrete Gabor wavelet transform is given by a convolution (Tuceryan and Jain, 1998): where s and t are the filter mask size variables, ψ * mn is the complex conjugate of ψ mn and m and n specify the scale and orientation of the wavelet respectively, with m = 0, 1, . . ., M −1, n = 0, 1, . . ., N −1.Here, the mean and standard deviation of the magnitude of the transformed coefficients are used in clustering.

Semi-variogram features
Semi-variograms are the basic tool for geo-statistics and have been used in a wide range of remote-sensing applications such as damage assessment and change detection (Olmo and Hernández, 2006;Sertel et al., 2007).Different texture features can be extracted from semivariograms, e.g.simplevariogram, radogram, etc. (Olmo and Hernández, 2006).In this case, simple-variogram, which better differentiates between the roof and object pixels in the pre-event image and between damage area and the intact roof on the post-event image, was selected as an applied feature in the segmentation process.This feature can be calculated using Eq. ( 5).
where γ k (h) is the value of variogram with different variogram range h, DN are the digital values of pixels x i and x i + h and N (h) is the number of couple points whose distance is h in an image region.All of the aforementioned features are powerful descriptors to be used along with spectral information.However, other texture features may be considered for damage assessment.For a more comprehensive review on texture analysis techniques, the reader is referred to Tuceryan and Jain (1998).

Segmentation
The aim of this step is to apply the above-mentioned features to distinguish the roof and non-roof pixels.To achieve this, an image segmentation concept can be applied as a promising tool.Image segmentation may be defined as the process of stratifying a digital image into multiple segments, in which these segments cover the entire image.Image segmentation has been, and still is, one of the challenging topics in computer vision, and several segmentation methods, such as thresholding, dplit-and-merge, region growing, etc., have been proposed in the literature (Cufi et al., 2001).The reader is also referred to Zhang (2006) for a good review on the developed techniques.
Image segmentation can also be performed using clustering, which is one of the important tools in machine learning and computer vision.Clustering can be defined as the grouping of objects that are similar to each other.During clustering, objects based on their properties are categorized into a few clusters, in which similar objects belong to the same cluster and dissimilar objects are assigned to different clusters.In this research, therefore, based on the extracted spectral and textural information of pixels in a building area, the objects are divided into a few groups.
Different applications of clustering such as image segmentation and information retrieval have resulted in numerous techniques for data clustering.These techniques can be classified into several categories.A good survey of these techniques is available in Jain et al. (1999).
For our intention, the best clustering algorithm is one which is powerful in handling the uncertainty of distinguishing roof and non-roof pixels.For this purpose, the Fuzzy C-Means (FCM) method, which is more favorable in comparison to the traditional methods such as k-means algorithm at avoiding local minima, is applied.This method is discussed at length in the proceeding paragraphs.
Fuzzy c-means (FCM), which is frequently used in computer vision, is a method of clustering that allows one object to belong to two or more clusters (Bezdek, 1981).This method is based on minimization of the objective function J m in Eq. ( 6).
where m is a real number greater than 1, x i is the i-th measured data, u ij is the degree of membership of x i in the cluster j , c j is the center of cluster j and ||x i −c j ||is the distance measure between the data x i and the cluster center c j .Fuzzy partitioning is performed through an iterative process to optimize the objective function J m .In each iteration, membership u ij and the cluster centers c j are calculated using Eqs.( 7)-( 8): The iteration will continue until the convergence condition of max ij {|u The aforementioned steps of FCM algorithm can concisely be composed as follows: U(k) using Eq. ( 8).
3. Update U matrix using Eq. ( 7).Based on the preceding paragraphs, by using FCM algorithm the pixels on the building area can be stratified into a number of clusters.Using textural features along with spectral features causes pixels with similar texture to belong to the same cluster.Based on the homogeneity of the roof part of the building area, it can be expected that the clustering algorithm results in grouping of all the roof pixels in the same cluster.And, with regard to the number of clusters, other pixels would belong to other clusters.
One of the crucial parameters in any clustering is the user decision on the number of clusters prior to the algorithm.Here, only two clusters of roof and non-roof pixels are considered in clustering the pre-event building area.After an earthquake, based on the magnitude of earthquake and condition of a building, some parts of the roof might be destroyed.Therefore, an additional cluster in clustering of the post-event building area is considered.In other words, based on the interpretation of pre-and post-event building areas, the numbers of clusters should be assigned two and three, respectively.Different alternatives were tested for numbers of clusters.In all the cases, because the similarity of the roof pixels they tended to be included in the same clusters and no changes were appearing on roof clusters; therefore, extra clusters were included by non-roof pixels.In other words, the method is not sensitive to the number of clusters.Therefore, two and three clusters can be considered for pre-and post-event building area, respectively.
The next step of the roof detection process is to recognize the pre-and the post-event roof clusters among the resulting clusters, which is described in the forthcoming section.

Roof recognition
In comparison with supervised classification, which results in the class label for any object, clustering only groups similar objects and gives no information about the cluster label of the objects.In other words, in the segmentation step, in which a clustering technique is applied, the pre-and the postevent roof clusters are not distinguished.This problem for pre-event building area can be solved based on the fact that the number of roof pixels is considerably greater than the other pixels (see Fig. 2).Therefore, the following algorithm can recognize the roof cluster between two resulting clusters of the pre-event building area: pre-event roof cluster == cluster #1; if size (cluster #2) >size (cluster #1) then pre-event roof cluster == cluster #2 end As some parts of the roof might be destroyed after an earthquake, the abovementioned fact cannot be assumed for post-event building area.However, it is clear that the textural and spectral properties of the post-event roof cluster should be similar to the pre-event roof.In other words, among the three resulting clusters of the post-event building area, the closest cluster to the pre-event roof cluster may be the intact roof cluster.In this case, the Euclidian distance between the cluster centers in feature space can be considered as the similarity measure.However, if the total area of a roof has been destroved, none of the clusters in the post-event segmentation should be labeled as roof cluster.Hence, a threshold for the minimum distance should be assigned, whereby in the case of exceeding the distance from this threshold, the intact roof cluster would be assigned as an empty cluster.To sum up, the intact roof cluster can be recognized using the following algorithm: J * == argMin {d(C j , C pre−event roof cluster )}, j = 1,2,3 intact roof cluster = cluster #J* if d J * >threshold then intact roof cluster = Ø end As can be seen from the above algorithms, by applying the pre-event image along with the post-event image the intact roof can be automatically recognized without any training dataset.This fact can decrease the time of processing, which is very important for disaster management.Also, the roofs of the area are not assumed as the same because each intact roof is recognized using its pre-event roof cluster, and no training dataset is needed.

Roof modification
The obtained pre-and post-event roof clusters may involve some meaningless pixels that should be eliminated.The pixels of a roof area can be shown as a binary image, in which the roof pixels are represented by a value of 1 and non-roof pixels by a value of 0. The binary image simplifies the interpretation as well as modification of the roof areas.As morphological operators are powerful tools to deal with binary images, the modification of the roof clusters is performed using these operators.In this case, Opening and Closing are applied to smooth the extracted roof and to fill unexpected holes, respectively.
The structural element plays an important role in using morphological operators.In an urban area, a remained small roof area is not a meaningful roof and cannot be considered.So the structural element should be somehow assigned to eliminate the small areas.Smaller structural elements cannot eliminate meaningless parts of roof areas.On the other hand, a bigger structural element may eliminate some meaningful roof areas.Therefore, the structural elements should be somehow selected so that more reliable results would be obtained.It should be noted that the spatial resolution of a satellite image has straight influence on the size of the structural element.For example, one may consider the area less than 2 m 2 is a meaningless area to be taken into account in damage assessment.In this case, for an image with 1 m spatial resolution, 2 × 2 pixels structural elements can be an appropriate size.

Damage assessment
One can determine the amount of damage for each candidate building using the extracted post-event roof compared to its initial condition.For instance, the damage degree can be measured as the relation between the area of the intact roof and the initial area.This method of damage estimation based on a single parameter such as area is a rather simplistic procedure, given that on some occasions the estimated damage degree might not be in conformity with reality.For example, if an extracted post-event roof involves a number of small and meaningless intact roof parts, the sum of all the areas would result in a high value for the roof area and as a result the building would get the wrong damage grade.
1 Figure 3.A simulated pre-and post-event building areas are depicted by 2 shade areas.Parameters of bounding boxes and minimum area ellipses c 3 evaluate the rate of damages.Where "a" and "b" are the major and minor a 4 "φ" is the angular difference between major axes of two ellipses, "B" is the 5 "A" is the roof area.Zero and one indices indicate pre-event and post-ev 6 respectively.7 Fig. 3. Simulated pre-and post-event building areas are depicted by light-and dark-shade areas.Parameters of bounding boxes and minimum area ellipses can be applied to evaluate the rate of damages.Where "a" and "b" are the major and minor axes of ellipse "E", "ϕ" is the angular difference between major axes of two ellipses, "B" is the bounding box and "A" is the roof area.Zero and one indices indicate pre-event and post-event parameters, respectively.
To avoid this problem, one may only take the biggest roof part into account during the estimation but this means not using some parts of information.In another case, a very long and narrow-shaped post-event roof may results in a high area and, consequently, low level damage degree of a building, while the remaining roof may be of low relevance.
Therefore, the damage degree should be estimated based on a comprehensive observation of the shape and structure of the post-event roof area, which can be performed using shape analysis techniques.However, damage assessment based on these properties is not generally deterministic but is characterized by some level of fuzziness or uncertainty.Therefore, deterministic analysis of these descriptors may not produce a trustworthy result and the vague and uncertainty of the problem should be considered.
Fuzzy theory (Zadeh, 1965), which resembles human reasoning in its use of approximate information to generate decisions, is known as a useful tool in dealing with these types of problems (Cox, 1999;Zimmermann, 1996).In this research, therefore, a fuzzy rule-based system is introduced in analyzing the shape descriptors.In the following paragraphs, after describing the structural descriptors, the designed fuzzy rulebased system is described and implemented.

Shape analysis
Shape analysis may be defined as the process of extracting structural descriptors which can comprehensively describe the geometry and shape of a specific area inside an image.One can get a basic idea about a shape based on a number of primitive parameters such as area and length to width ratio.A useful shape descriptor for damage assessment should Table 2. Applied shape features for evaluating the building damage degree.Where "a" and "b" are the major and minor axes of ellipse "E", "ϕ" is the angular difference between major axes of two ellipses, "B" is bounding box and "A" is the roof area.Zero and one indices indicate pre-event and post-event parameters, respectively.

Name
Description Eq. either provide comprehensive information about the structure of the intact roof area or be able to measure differences between pre-and post-event roof areas.Simulated pre-and post-event roof areas for a sample building are shown in Fig. 3.As can be seen, many shape features can be extracted from a roof area.A very basic parameter to show the value of damage is the ratio of the postevent roof area to the pre-event roof area ((A 1 1 + A 2 1 )/A 0 ).A higher value of this parameter indicates a smaller change of the roof, and consequently a lower degree of damage.

RoofsAreaRatio
The ratio of the bounding boxes' area (B 1 /B 0 ) can similarly measure the change of the roof.However, this feature is not always compatible with reality.For example, when the intact roof includes separated small parts, the area of the post-event bounding box would be high.Therefore, to show the compactness of the post-event roof area, another measure is needed along with bounding boxes ratio.For this purpose, bounding box filling of the post-event roof can be used.This feature can be calculated by the ratio of the post-event roof area to the area of the bounding box ((A 1 1 + A 2 1 )/B 1 ).High value of bounding box filling means the bounding box ratio can be applied as the same as roofs area ratio in damage assessment.
The minimum area enclosed ellipses around the roofs can also provide useful information for damage assessment.The more similar the two ellipses are, the lower level of damage.Two ellipses can be compared using their structural parameters, such as major and minor axes.In this research, major axes ratio (a 1 /a 0 ), minor axes ratio (b 1 /b 0 ) and direction stability (1 − ϕ/90) are applied for describing the value of damage (see Fig. 3).Similar to the bounding box ratio, which is applied along with bounding box filling feature, ellipses similarity measures should be applied along with a complementary feature.Here, the ellipse filling measure ((A 1 1 +A 2 1 )/E 1 ) can be applied as a good feature.A more compact and unit post-event roof area is more meaningful to be considered in damage assessment.All the applied descriptors used in this research are listed in Table 2.
All the above-mentioned descriptors are very helpful for damage assessment, and by using these descriptors, an expert can precisely estimate the damage degree.However, the main goal of all damage assessment techniques is to increase the level of automation in damage degree estimation of buildings.In other words, these parameters should be automatically analyzed by a computer for rapid damage assessment.Such an analysis is not totally deterministic and involves a level of uncertainty.Therefore, a powerful decision making system that is able to resemble an expert's thought process in handling the uncertainty is required.For this purpose, a fuzzy rule-based system, which is known as a proper tool for handling uncertainty in solving problems, is applied for analyzing the shape features.

Damage assessment using a fuzzy decision making system
Every decision making process is not generally deterministic but is usually characterized by some level of fuzziness or uncertainty.Yet traditional decision making systems do not provide a good mechanism for coping with uncertainty.Fuzzy set theory, which was triggered by these considerations, provides a conceptual framework for solving non-deterministic problems in an ambiguous environment.In this research, we use a fuzzy rule-based system for analyzing the extracted shape features during the damage assessment process.Here, firstly a general structure of the fuzzy rule-based system is described and then more detail of the designed fuzzy decision making system for damage assessment is presented.
A fuzzy rule base (or fuzzy system) used for decision making is generally comprised of three principal steps of fuzzification, inference and defuzzification, as shown in space into fuzzy subspaces, each specified by a fuzzy membership function.Fuzzy rules are then generated from each fuzzy subspace.The second step, inference, requires the calculation of the strength of each rule being triggered.The final step, defuzzification, aggregates all triggered rules and generates a non-fuzzy output.

Fuzzification
The purpose of fuzzification is to partition the feature space into fuzzy subspaces and generate rules for each fuzzy subspace.Note that all fuzzy subspaces normally overlap each other to some degree.To carry out the process of fuzzification, one must first define membership functions in order to calculate the membership grade for the input elements.Although the fuzzy membership function can take any form (as long as the function can map the inputs onto the range [0, 1]), four kinds of fuzzy membership functions, known as monotonic, triangular, trapezoidal, and bell shaped, are the most frequently used in fuzzy rule base experiments (Tso and Mather, 2009).The selection of membership functions and the width of each fuzzy subspace are certainly case dependent.
In this research, seven extracted shape features from shape analysis step as input variables, and the degree of damage as output variable, are considered for damage assessment.Depending on the sensitivity of the variable, an expert assigns a number of linguistic labels to each variable (input/output), which reflect an interactively carried-out examination of all possible values of the variables.In practice, this assignment is mostly a mixture of expert knowledge and examination of the desired input-output data.

Inference
Fuzzy sets and fuzzy operators are the "subjects" and "verbs" of fuzzy logic (Samadzadegan et al., 2005).In order to create a useful statement, complete sentences have to be formulated.Conditional statements, IF-THEN rules, are statements that make fuzzy logic useful.A single fuzzy IF-THEN rule can be formulated according to: A and B are linguistic labels defined by fuzzy sets on the range of all possible values of x and y, respectively.The IF part of the rule "x is A" is called antecedent, the THEN part of the rule "y is B" is called consequent.The antecedent is an interpretation that returns a single number between 0 and 1, whereas the consequent is an assignment that assigns the entire fuzzy set B to the output variable y.The antecedent may integrate several inputs using logical AND and OR operators.
Fuzzy reasoning with fuzzy IF-THEN rules enables linguistic statements to be treated mathematically.For example, for estimating the damage degree of a building, one of the IF-THEN fuzzy rules might be the following: IF Roof-sAreaRatio is VeryLarge AND BoundingBoxFilling is Filled AND EllipseFilling is Filled, THEN DamageDegree is Neg-ligibleDamaged.This example reveals an important aspect of fuzzy reasoning, which is that the rule base should include observations of the important descriptors.Moreover, it reflects the fact that people may formulate similar "fuzzy statements" to characterize how they perceive negligible damage degree.
Formulating the rules is more a question of the expertise of an operator than of a detailed technical modeling approach.Given the rules and inputs, the degree of membership to each of the fuzzy sets has to be determined.For the above example, the input variables are: RoofsAreaRatio, BoundingBox-Filling and EllipseFilling, and the output variable is Dam-ageDegree.
If the membership grades are equal to one (i.e. the rule condition is fully satisfied), the THEN clause in the rule should be fully adopted (i.e. with full strength).On the other hand, if the rule condition is only partially satisfied, the THEN clause should be partially weighted.Two weighting approaches, known as multiplication and minimization, are commonly used.
In addition, the rules being triggered can be numerous because the fuzzy membership functions normally overlap.Hence, a feature value falling within the overlap area will simultaneously trigger several rules.Since the result of rule aggregation is a membership function, a defuzzification process has to be implemented in order to obtain a deterministic value.

Defuzzification
Several kinds of defuzzification strategies, such as the center of gravity and mean of maximum, have been suggested in the literature (Zimmermann, 1996).The most applied defuzzification method is to calculate the center of gravity, which determines the center of the area under the aggregated output function.The center-of-gravity method for discrete data can be calculated from the following equation: where n is the number of elements of the sampled membership function, and µ(s) is the membership grade of measurement s.
Because these fundamentals of fuzzy logic are well described in textbooks, we do not want to go into the theoretical background at this point.For a comprehensive study of fuzzy logic, please refer to Zimmermann (1996) and Cox (1999).
Using these three steps (fuzzification, inference, and defuzzification), one eventually reaches a deterministic value for the damage degree.Finally, regarding this value, each building can be represented by a specific color on a map in order to generate a damage map.

Experiments and results
To assess the efficiency of the proposed damage assessment method, two high resolution satellite images and 1 : 2000 relevant vector map of the city of Bam, Iran, are used.The before and after 26 December 2003 Bam earthquake images were acquired on 30 September 2003 and 3 January 2004, respectively, by the QuickBird satellite.
In the pre-processing step, the images were atmospherically corrected using the FLAASH atmospheric correction module of the ENVI ® image processing software package and the images were orthorectified by means of the SRTM Digital Terrain Model to compensate for the ground elevation.Then, pan-sharpened images were created using wavelet fusion technique.Also, histogram equalization and histogram matching were performed to increase the spectral similarity of the images.Next, the images were precisely registered to the map using 15 well distributed control points.The residuals of 10 well distributed check points did not exceed 43 and 48 cm on pre-and post-event images, respectively.Moreover, the co-registration accuracy of the images was measured by comparing check points coordinates on both images, where 7 cm of RMS was observed.Finally, a 2500 × 1900 pixels area, including 1136 buildings, was selected as the test area.The QuickBird images after preprocessing step and the applied vector map of the test area are shown in Fig. 5. igure 6. Extraction of building areas of a candidate building (shown by blue line on the vector map) from images by using geo-referencing information.

Geo-referencing Information
Candidate Building

Pre-event HRSI
Vector map Post-event HRSI Fig. 6.Extraction of building areas of a candidate building (shown by blue line on the vector map) from images by using georeferencing information.
After pre-processing, the algorithm is executed on every building in the test area.Using available geo-referencing information, corner points of a candidate building are located on both images and the building areas are extracted.Figure 6 shows the process of building areas extraction for a sample candidate building.
In order to detect the roof areas based on the segmentation technique, textural information along with spectral information of pixels are applied.In this research, among the 21 implemented texture features, six features (bold in Table 1) were selected to be used along with spectral features in the segmentation step.The textural features for each pixel were extracted using its 3 × 3 neighborhood pixels in the panchromatic images; however, a greater size such as 5×5 or 7×7 may be applied.During Gabor features extraction, a filter bank consisting of Gabor filters with three scales and four rotations were considered.
Except for the semi-variogram feature, all the textural features assign a higher value for roof pixels and a lower value for non-roof pixels.However, we applied the inverse value of this feature to make the simple-variogram feature consistent with the other features.The extracted features are applied in the segmentation step to distinguish the roof and non-roof pixels.
Number of clusters is one of the most crucial parameters in clustering.Here, the smallest acceptable numbers of clusters, two and three clusters for pre-and post-event building areas, were considered.Because of the similarity of the roof pixels extra clusters involve non-roof pixels.In other words, the numbers of clusters do not considerably influence the results and only small changes appear in the final damage degree.Therefore, considering these clusters numbers, FCM algorithm was executed on both pre-and post-event building areas.The segmentation results of the sample building are depicted in Fig. 7.As can be seen, the roof segment on the pre-event building is absolutely greater than the other cluster and can be easily distinguished.By calculating the Euclidian distance of the cluster center of this segment and the postevent roof cluster, the post-event roof cluster can be distinguished.The maximum accepted difference for the distance is considered as 0.2, and upon exceeding this threshold the building is labeled as a totally damaged building.
Moreover, the modifications of the extracted roofs are performed using opening and closing morphological operators.One of the most important parameters in this step is the size of the structural element.In this paper, based on the spatial resolution of QuickBird images (0.61 m), 2 × 2, 3 × 3, 4 × 4 and 5×5 square type structural elements were considered, in which more promising results was observed by applying 3×3 square type (equal to 1.83 × 1.83 m 2 ).In Fig. 7, the modified roof of the sample building is depicted.As can be seen, after modification, noisy roof pixels and holes are successfully cleared from the roofs.
Once the roofs are detected, the shape analysis steps are executed to extract proper features that can help to estimate the damage degree.For this purpose, the bounding box and the minimum area enclosing the ellipse around the roofs are found.The bounding box may be found using minimum and maximum coordinates of roof pixels.In this research, the minimum area enclosing the ellipse was calculated based on the Khachiyan algorithm (Khachiyan and Todd, 1993).In calculating this ellipse, one may use only the edge pixels of a roof area to simplify the calculation because inside pixels of the roof area do not influence the parameters of the ellipse.
In this research, the seven shape features of RoofsArea-Ratio, BoundingBoxRatio, BoundingBoxFilling, MajorAxes-Ratio, MinorAxesRatio, DirectionStability and EllipseFilling are calculated through shape analysis of pre-and post-event .Extracted shape features through shape analysis of the roof area for the sample building.roof areas.Figure 8 shows the extracted shape features for the sample building.
Extracted shape features are applied as input variables of the fuzzy inference system for estimating the damage degree, which is the output of the fuzzy system.A number of linguistic labels are assigned to each variable (input/output), which reflect an interactively carried-out examination of all possible values of the variables.Linguistic labels of input and output variables of our damage assessment FIS are listed in Table 3.
In this step, a membership function can be defined for each linguistic variable.Accurate definition of the membership functions is of high importance in any fuzzy decision making system.The shape and the values of the membership functions should be accurately defined by an expert based on his/her experience in damage degree estimation.Here, trapezoidal-and triangular-shaped functions were applied.These membership functions, depicted in Fig. 9, are defined to the system based on user experience.In this figure, membership functions of input and output variables are also depicted.In order to import user knowledge in the fuzzy reasoning system, 85 rules are constructed.Some of the employed IF-THEN rules are shown in Table 4. Here, Mamdani FIS, one of the most commonly used fuzzy engines, is used for making decisions using fuzzy rules.In this table two examples are given for each damage degree.
By applying the designed fuzzy decision making system for the sample building, a damage degree of 44.20 % was estimated.To evaluate the proposed method, the resulted damage value may be compared with manually-estimated damage degree by an expert.For this purpose, here the meaningful parts of the post-event roof area were manually detected and extracted by an expert for calculating the damage degree.For the sample building, which is shown in Fig. 10b, a damage value of 40.05 % was obtained by an expert.As can be seen, this roof is close to the resulted post-event roof from the algorithm (see Fig. 10a and c).The difference between the two resulted damage degrees for this building is 4.15 %, or about 10 % of the damage degree, which shows the ability of this method in assessing the damage value of this sample building.
The process of the proposed damage assessment method for three sample buildings with IDs 7, 150 and 541 are also shown step-by-step in Fig. 11.As can be seen from the table, the algorithm has successfully assessed the damage degree for these buildings.
According to the degree of damage, buildings can be categorized into different levels of damage for a better representation of damage map.Here, five degrees of damage were considered: 0-20 % damaged (shown by dark green), 20-40 % damaged (shown by light green), 40-60 % damaged (shown by yellow), 60-80 % damaged (shown by light red), and 80-100 % damaged (shown by dark red).The final resulting damage map was generated by implementing the proposed algorithm on the test area, which is depicted in Fig. 12.In the resulted damage map of the test area, a majority of buildings are labeled by fourth grade damage, while a minority are labeled as first grade.
Unfortunately, no field observation dataset was available for evaluating the results.Therefore, accuracy assessment of the proposed damage assessment technique was carried out using 325 randomly selected buildings, which had their preand post-event roofs manually measured by an expert.Distribution of these buildings is shown in  In this research, two different accuracy investigations were performed to assess the quality of the proposed method.In the first investigation, the manually extracted post-event roof of checked buildings and their counterpart, resulting from the proposed algorithm, were compared.Minimum, maximum and average of differences were 0, 34 and 9 pixels, respectively.Moreover, the mode of the differences was 6, showing that for most of the buildings the extracted post-event roof area from the algorithm was almost equal to their visually extracted roof area.In Fig. 14, some of the checked buildings with different damage degrees are illustrated.As can be seen, the algorithm has successfully extracted different postevent roof areas.
The confusion matrix is used as an indication of the properties of a classification which contains the number of elements that have been correctly or incorrectly classified for each class (Rokach, 2010).It can be seen on its main diagonal the number of observations that have been correctly classified for each class; the off-diagonal elements show the number of observations that have been incorrectly classified.
Based on the values of the confusion matrix, one can calculate a set of parameters such as overall accuracy and kappa coefficient to describe the classification results.
Using visually extracted post-event roof areas, damage value for each checked building was calculated and assigned pre-defined damage grades.Considering the resulted damage grade as a reference data and comparing to the algorithm results, a confusion matrix as another accuracy investigation was obtained, with an overall accuracy of 90.46 % and a kappa coefficient of 86.68 % using our method.The resulting confusion matrix is shown in Table 5.
As can be seen from Table 5, the reference data includes only four buildings in the first degree of damage, with all of them being successfully labeled by the proposed method.However, one extra building from the second degree of damage has mistakenly been labeled as the first degree of damage.Also, the confusion matrix shows that the algorithm did not produce very strong results in differentiation between fourth and fifth classes, while this is also often difficult for manual interpretation.On the other hand, the average user or producer's accuracy in the table proves the ability of this damage assessment method.
Logically, post-event roof boundary should be inside the pre-event boundary.However, this may not always happen.As can be seen from Fig. 14, in some cases the red boundary, which is the post-event roof boundary, exceeds the blue boundary.This may happen for two reasons, namely registration error or collapsed roof.In this research, this error, which has a mild effect on the final result, is not considered.
In comparison with previous studies, high promising results were obtained from the implementation of the proposed method.For example, In Chesnel et al. (2008), damage assessment of the Bam and Boumerdes earthquakes using postevent QuickBird images were performed through SVM classification method.The test area of the Bam earthquake included 2168 buildings.In that study, average performances close to 75 % when four damage classes were discriminated, up to 90 % for an intact/damaged detection, were reported.In another study, (Rezaeian, 2010), three different classification algorithms were applied for damage assessment of the Bam earthquake using QuickBird image.included 890 buildings, where 79 % overall accuracy and 67 % kappa coefficient were reported.In these studies, they used a training dataset for training their algorithm.On the one hand, the collection of a training dataset is time consuming and reduces the level of automation.On the other hand, using the same dataset for the whole area is only applicable for cities in which roof buildings are the same.This may be one of the causes for not achieving better results.In the proposed method, however, this drawback was handled by using the pre-event image.In this case, each roof area on the post-event image is automatically detected based on its pre-event information and no training dataset is needed.Moreover, none of the previous works took the shape and structural properties of the intact roofs into account during damage value estimation.As our results show, considering shape information causes a more realistic damage degree.

Conclusions
In this research, a new automatic method for post-earthquake damage assessment using pre-and post-event high resolution satellite images has been presented.In the proposed algorithm, after pre-processing of the satellite images, building areas are extracted using a 1 : 2000 vector map.Then, intact roofs of candidate buildings are extracted through a clustering algorithm by applying textural and spectral information.Finally, analysis of structural information of the intact roof areas using a fuzzy inference system allows for estimation of the degree of damage.
By evaluating the proposed method using the available dataset of the city of Bam, Iran, 90.46 % overall accuracy and a kappa coefficient of 86.68 % were obtained.These results prove the capability and high ability of this method for building damage map creation using high resolution satellite images.Compared to the same studies, the results are promising.Considering shape and structural information of the intact roof during damage assessment of the buildings, which was presented here for the first time in damage assessment, provides a realistic damage assessment and consequently more accurate results.Also, applying a pre-event image to eliminate the necessity of training dataset collection, which is a time consuming process, increases the level of automation in damage assessment.Moreover, using a fuzzy rule-based decision making system in damage assessment to handle the uncertainty is another selling point of the proposed method.
The inability to recognize totally collapsed buildings in cases when the roof has completely fallen down but is not fully destroyed may be one of the main drawbacks of this algorithm.This can be handled by using post-event DSM of the area.The proposed algorithm needs an accurate and updated vector map, which is available for almost all cities around the world, to locate buildings on the images.However, this method for areas where the vector map is outdated may not obtain reliable results.Therefore, the future studies may go into the way of not using vector map, in which case buildings may be recognized only from a pre-event image and located on the post-event image using co-registration information.Nonetheless, comparing a building from a postevent image to its counterpart on the pre-event image is the main reason for the case-independency of the algorithm.In other words, by accurately pre-processing images and using appropriate features in clustering, the algorithm should be reliable for other high resolution satellite imagery.Accordingly, the algorithm could be tested using various datasets with the same or different high spatial resolution sensor in future studies.

Figure 1 .Fig. 1 .
Figure 1.Flowchart of the proposed method.Pre-and post-event high resolution satellite 3 images (HRSI) along with vector map are the input data and damage map is generated.4 5 6

Figure 2 .Fig. 2 .
Figure 2. Building area on a 1:4000 scale aerial photo.(a) a single building (b) a city block.Several objects such as chimney and cooler may exist on a roof area.These objects are appeared as spots on satellite images.
the Bounding Box Around the Post−event Roof the Area of the Bounding Box Around the Pr e−event Roof B 1 B 0 BoundingBoxFilling the Area of the Post−event Roof the Area of the Bounding Box Around the Post−event Roof the Post−event Roof the Area of the minimum Area Enclosing Ellipse of the Post−event Roof of the Enclosing Ellipse Around the Post−event Roof The Major Axis of the Enclosing Ellipse Around the Pre−event Roof a 0 a 1 MinorAxisRatio The Minor Axis of the Enclosing Ellipse Around the Post−event Roof The Minor Axis of the Enclosing Ellipse Around the Pre−event Roof

Figure 7 .Fig. 7 .
Figure 7. Roofs detection process results for a sample building.(a) Building area (b) 1 Segmentation (c) Roof area recognition (d) Modification 2 3 Fig. 7. Roofs detection process results for a sample building.(a) Building area (b) Segmentation (c) Roof area recognition (d) Modification.

Fig. 8 .
Fig. 8. Extracted shape features through shape analysis of the roof area for the sample building.

Fig. 13. 40 Figure 9 .Fig. 9 .
Figure 9. Membership functions of the input and the output linguistic variables.1 2 Fig. 9. Membership functions of the input and the output linguistic variables.

Fig. 10 .
Fig. 10.Final damage assessment results of the candidate building.(a) Result of the proposed method.(b) Manually extracted roof by an expert.(c) Superimposed roofs on the post-event image.

Fig. 11 .
Fig. 11.Output results with detail about damage assessment for three candidate buildings.

2.
Final damage map of the test area using the proposed method.Depending on the degree, each building is shown by a specific colour.Five degrees of damage are : 0-20% damaged (shown by dark green), 20-40% damaged (shown by light green), damaged (shown by yellow), 60-80% damaged (shown by light red) and 80-100% damaged (shown by dark red).

Fig. 12 .
Fig. 12. Final damage map of the test area using the proposed method.Depending on the damage degree, each building is shown by a specific color.Five degrees of damage are considered: 0-20 % damaged (shown by dark green), 20-4 % damaged (shown by light green), 40-60 % damaged (shown by yellow), 60-80 % damaged (shown by light red) and 80-100 % damaged (shown by dark red).

Table 1 .
Implemented texture features in texture analysis step.From this list, six selected features which make a better difference between roofs and non-roof pixels are bold.

Table 3 .
Linguistic variables and labels for the fuzzy-based damage assessment process.Small, Medium, Large, Very Large BoundingBoxRatio Very Small, Small, Medium, Large, Very Large BoundingBoxFilling Not Filled, Moderately Filled, Filled EllipseFilling Not Filled, Moderately Filled, Filled MajorAxesRatio Very Small, Small, Medium, Large, Very Large MinorAxesRatio Very Small, Small, Medium, Large, Very Large DirectionStability Not Parallel, Nearly Not Parallel, Nearly Parallel, Parallel Output DamageDegree Negligible Damage, Moderate Damage, Substantial Damage, Heavy Damage, Complete Damage

Table 4 .
Fuzzy rules for estimating damage degree of a building.IF RoofsAreaRatio is VeryLarge AND BoundingBoxFilling is Filled AND EllipseFilling is Filled, THEN DamageDegree is NegligibleDamage.-IF RoofsAreaRatio is VeryLarge AND MajorAxesRatio is VeryLarge AND MinorAxesRatio is VeryLarge and Direction-Stability is Parallel AND EllipseFilling is Filled, THEN DamageDegree is NegligibleDamage.-IF RoofsAreaRatio is Large AND BoundingBoxFilling is Filled AND EllipseFilling is Filled, THEN DamageDegree is ModerateDamage.-IF RoofsAreaRatio is VeryLarge AND BoundingBoxRatio is VeryLarge AND BoundingBoxFilling is ModeratelyFilled AND EllipseFilling is ModeratelyFilled, THEN DamageDegree is ModerateDamage.-IF RoofsAreaRatio is Medium AND BoundingBoxFilling is Filled AND EllipseFilling is Filled, THEN DamageDegree is SubstantialDamage.-IF RoofsAreaRatio is Large AND BoundingBoxRatio is Large AND BoundingBoxFilling is ModeratelyFilled AND Ma-jorAxesRatio is Large AND EllipseFilling is ModeratelyFilled, THEN DamageDegree is SubstantialDamage.
-IF RoofsAreaRatio is Small AND BoundingBoxFilling is Filled AND EllipseFilling is Filled, THEN DamageDegree is HeavyDamage.-IF RoofsAreaRatio is Medium AND BoundingBoxRatio is Medium AND BoundingBoxFilling is NotFilled AND Ma-jorAxesRatio is Medium AND MinorAxesRatio is Medium AND EllipseFilling is NotFilled, THEN DamageDegree is HeavyDamage.-IF RoofsAreaRatio is VerySmall AND BoundingBoxFilling is Filled AND EllipseFilling is Filled, THEN DamageDegree is CompleteDamage.-IF RoofsAreaRatio is Small AND BoundingBoxRatio is Small AND BoundingBoxFilling is NotFilled AND EllipseFilling is NotFilled, THEN DamageDegree is HeavyDamage.

Table 5 .
Confusion matrix obtained by considering the manually observed damage grade as reference data and comparing it to the algorithm results.In this algorithm, overall accuracy of 90.46 %, kappa coefficient of 86.68 %, average producer accuracy of 94.25 % and average user accuracy of 86.80 % were obtained.