Csaba Benedek

Xavier Descombes

Fellow, IEEE Josiane Zerubia

Building Development Monitoring in Multitemporal Remotely Sensed Image Pairs with Stochastic Birth-Death Dynamics

Keywords: Building extraction, change detection, Marked Point Processes, Multiple Birth and Death Dynamics !

In this paper we introduce a new probabilistic method which integrates building extraction with change detection in remotely sensed image pairs. A global optimization process attempts to find the optimal configuration of buildings, considering the observed data, prior knowledge, and interactions between the neighboring building parts. We present methodological contributions in three key issues:

(1) We implement a novel object-change modeling approach based on Multitemporal Marked Point Processes, which simultaneously exploits low level change information between the time layers and object level building description to recognize and separate changed and unaltered buildings. (2) To answering the challenges of data heterogeneity in aerial and satellite image repositories, we construct a flexible hierarchical framework which can create various building appearance models from different elementary feature based modules.

(3) To simultaneously ensure the convergence, optimality and computation complexity constraints raised by the increased data quantity, we adopt the quick Multiple Birth and Death optimization technique for change detection purposes, and propose a novel non-uniform stochastic object birth process, which generates relevant objects with higher probability based on low-level image features.

INTRODUCTION

F OLLOWING the evolution of built-up regions is a key issue of aerial and satellite image analysis. Although the topic has been extensively studied since the 80's, it has had to continuously face the challenges of the quickly evolving quality and quantity of remotely sensed data, the richness of different building appearances, the data-heterogeneity in the available image repositories and the various requirements of new application areas.

Input Data

Numerous methods in the bibliography address building extraction at a single time instance [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF], [START_REF] Lafarge | Structural approach for building reconstruction from a single DSM[END_REF], [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF]. It is common to use multiview inputs [START_REF] Noronha | Detection and modeling of buildings from multiple aerial images[END_REF], [START_REF] Bignone | Automatic extraction of generic house roofs from high resolution aerial imagery[END_REF] to exploit 3-D information in building modeling. Detection in densely populated areas can be efficient by working on stereo-or lidar-based Digital Elevation/Surface Models (DEM/DSM), where the silhouettes of the building footprints can be separated from the ground planes by the estimated height data [START_REF] Lafarge | Structural approach for building reconstruction from a single DSM[END_REF], [START_REF] Ortner | A marked point process of rectangles and segments for automatic analysis of digital elevation models[END_REF], [START_REF] Champion | 2D building change detection from high resolution aerial images and correlation digital surface models[END_REF], [START_REF] Karantzalos | Large-scale building reconstruction through information fusion and 3D priors[END_REF]. Other benefits are provided by multiple sensor inputs such as fusion of aerial images with color infrared (CIR) [START_REF] Rottensteiner | Building detection by fusion of airborne laser scanner data and multi-spectral images: Performance evaluation and sensitivity analysis[END_REF], or laser data [START_REF] Jaw | Building roof reconstruction by fusing laser range data and aerial images[END_REF]. However, several image repositories from city suburbs and smaller settlements lack stereo or special sensor information. We address this case in the paper: building identification becomes here a challenging monocular object recognition task based on purely optical data [START_REF] Shufelt | Performance evaluation and analysis of monocular building extraction from aerial imagery[END_REF].

Multitemporal Information

Remote sensing image databases often contain multitemporal image samples from the same geographical areas. Exploiting the temporal information, change recognition and classification are nowadays important aspects of urban scene analysis. Several recent building change detection approaches [START_REF] Champion | 2D building change detection from high resolution aerial images and correlation digital surface models[END_REF], [START_REF] Rottensteiner | Automated updating of building data bases from digital surface models and multi-spectral images: Potential and limitations[END_REF] assume that a topographic building database is already available for the earlier time layer, thus the process can be decomposed into old model verification and new building exploration phases. On the other hand, when dealing with image repositories without any meta data, the task requires automatic building detection for each image.

In this paper, we solely use as input a registered pair of (2-D) images taken at several years time difference. Applying conventional stereo-matching algorithms in this case may face several difficulties. First, the scene content, the viewpoints and the image qualities of the two views may be significantly different, which can corrupt feature matching algorithms (e.g. corner point tracking) needed for 3-D structure extraction. Moreover, several databases contain images which are created by mosaicking separately taken aerial photos, and the components undergo different geometric corrections.

In the proposed approach we apply a 2-D building and change detection technique [START_REF] Tanathong | Object oriented change detection of buildings after the Indian ocean tsunami disaster[END_REF] which is less influenced by the above effects than stereo based approaches. Moreover, at a higher processing level, the proposed footprint extraction step may also contribute to 3-D building reconstruction [START_REF] Lafarge | Structural approach for building reconstruction from a single DSM[END_REF], [START_REF] Ortner | A marked point process of rectangles and segments for automatic analysis of digital elevation models[END_REF].

Change detection methods frequently rely on the assumption that changes occur very rarely, thus they can be identified through outlier detection using global image statistics [START_REF] Bruzzone | An adaptive semiparametric and context-based approach to unsupervised change detection in multitemporal remote-sensing images[END_REF]. However, in dynamically improving (sub-)urban areas this hypothesis is often invalid, and there is a need for solutions which are insensitive to the quantity of differences.

An object oriented change detection technique is introduced in [START_REF] Tanathong | Object oriented change detection of buildings after the Indian ocean tsunami disaster[END_REF] and applied to the extraction of damaged buildings after natural disasters. This method follows the Post Detection Comparison (PDC) approach, as independent building detection processes are applied for the two images, followed by object level comparison. However, the object detection phase can be corrupted by image noise, irregular structures or occlusions by vegetation [START_REF] Tanathong | Object oriented change detection of buildings after the Indian ocean tsunami disaster[END_REF] which may present missing or only partially extracted buildings to the object matching module. Moreover, this comparison may be affected by further intensity artifacts caused by shadows or altered illumination conditions.

Following another approach, several low level change detection methods have been proposed for remote sensing [START_REF] Benedek | Change detection in optical aerial images by a multi-layer conditional mixed Markov model[END_REF], [START_REF] Fournier | A contrast equalization procedure for change detection algorithms: applications to remotely sensed images of urban areas[END_REF], which work without using any explicit object models. They extract image regions which have been altered in an irregular way based on an appropriately selected set of features, such as color difference, texture or block correlation. Although these techniques are usually considered as preprocessing filters, there have not been many attempts given to justify how they can support the object level investigations. We take a step forward in this paper, and exploit interaction between object extraction and local textural image-similarity information in a unified probabilistic model. It will be shown that we can obtain additional evidences for the presence of new, modified or demolished buildings through detecting changes in relevant low level feature domains. As for unchanged buildings, the images of the two time instances provide multiple views of the same objects, which may increase the detection accuracy compared to relying on a single time layer.

Object and Configuration Models

Another important issue is related to modeling the building entities. The conventional bottom-up techniques [START_REF] Müller | Robust building detection in aerial images[END_REF] construct the objects from primitives, like roof blobs, edge parts or corners. Although these methods can be fast, they may fail if the primitives cannot be reliably detected. To increase robustness, it is common to follow the Hypothesis Generation-Acceptance (HGA) scheme [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF], [START_REF] Song | Building roof detection from a single high-resolution satellite image in dense urban area[END_REF]. Here the accuracy of object proposition is not crucial, as false candidates can be eliminated in the verification step. However, objects missed by the generation process cannot be recovered later, which may result in several false negatives. On the other hand, generating too many object hypotheses (e.g. applying exhaustive search) slows down the detection process significantly. Finally, conventional HGA techniques search for separate objects instead of global object configurations, disregarding population-level features such as overlapping, relative alignment, color similarity or spatial distance of the neighboring objects [START_REF] Lafarge | Structural approach for building reconstruction from a single DSM[END_REF].

To overcome the above drawbacks, recent inverse methods [START_REF] Descombes | Marked point processes in image analysis[END_REF] assign a fitness value to each possible object configuration, and an optimization process attempts to find the configuration with the highest confidence. This way, flexible object appearance models can be adopted, and it is also straightforward to incorporate prior shape information and object interactions. Marked Point Processes (MPP) [START_REF] Descombes | Marked point processes in image analysis[END_REF] are good candidates for addressing these challenges, since they can efficiently model the geometry of objects and deal with an unknown number of entities [START_REF] Ortner | A marked point process of rectangles and segments for automatic analysis of digital elevation models[END_REF], [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF], [START_REF] Lafarge | Geometric feature extraction by a multi-marked point process[END_REF]. However, this inverse approach needs to perform a computationally expensive search in a high dimensional population space, where local maxima of the fitness function can mislead the optimization. Due to the large databases, the optimization issue plays a particular role in remote sensing applications. In previous techniques [START_REF] Ortner | A marked point process of rectangles and segments for automatic analysis of digital elevation models[END_REF], [START_REF] Lacoste | Point processes for unsupervised line network extraction in remote sensing[END_REF], [START_REF] Lafarge | Geometric feature extraction by a multi-marked point process[END_REF] the optimization has been performed using a Reversible Jump Markov Chain Monte Carlo (RJMCMC) scheme, with implementations where each iteration perturbs one or a couple of objects, and the rejection rate, especially for the birth move, induces a heavy computation time. Besides, one should be very careful when decreasing the temperature, because at low temperature, it is difficult to add objects to the population.

Taking a different approach, we adopt here the Multiple Birth and Death Dynamic technique (MBD) [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF] for the change detection purposes. Unlike following a discrete jumpdiffusion scheme like in RJMCMC, the MBD optimization method defines a continuous time stochastic evolution of the object population, which aims to converge to the optimal configuration. The evolution under consideration is a birthand-death equilibrium dynamics on the configuration space, embedded into a Simulated Annealing (SA) process, where the temperature of the system tends to zero in time. The final step is the discretization of this non-stationary dynamics: the resulting discrete process is a non-homogeneous Markov chain with transition probabilities depending on the temperature, energy function and discretization step. In practice, the MBD algorithm evolves the population of buildings by alternating purely stochastic object generation (birth) and removal (death) steps in a SA framework. In contrast to the above RJMCMC implementations, each birth step of MBD consists of adding several random objects to the current configuration, which is allowed due to the discretization trick. Using MBD, there is no rejection during the birth step, therefore high energetic objects can still be added independently of the temperature parameter. Thus the final result is much less sensitive to the tuning of the SA temperature decreasing process, which can be achieved faster. Due to these properties, in selected remote sensing tasks (bird and tree detection) [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF] the optimization with MBD proved to be around ten times faster than RJMCMC with similar quality results. In addition, MBD has already been applied in different application areas, such as cell counting [START_REF] Hatsuda | Automatic cell identification using a multiple marked point process[END_REF] and video surveillance [START_REF] Utasi | A 3-D marked point process model for multiview people detection[END_REF].

Another key point is the probabilistic approach for object proposal. In several previous MPP applications [START_REF] Ortner | A marked point process of rectangles and segments for automatic analysis of digital elevation models[END_REF], the generation of object candidates followed prior (e.g. Poisson) distributions. On the contrary, we apply a data driven birth process to accelerate the convergence of MBD, which proposes relevant objects with higher probability based on various image features. In addition, we calculate not only a probability map for the object centers, but also estimate the expected object appearances through low-level descriptors. This approach uses a similar idea to the Data Driven MCMC scheme of image segmentation [START_REF] Tu | Image segmentation by Data-Driven Markov Chain Monte Carlo[END_REF]. However, while in [START_REF] Tu | Image segmentation by Data-Driven Markov Chain Monte Carlo[END_REF] the importance proposal probabilities of the moves are used by a jump-diffusion process, we should embed the data driven exploration steps into the MBD framework.

In this paper, we propose a novel multitemporal MPP (mMPP) model and an efficient bi-layer MBD (bMBD) optimization algorithm for the 2-D building change detection problem in remotely sensed image pairs. The present approach has been partially introduced in [START_REF] Benedek | Building extraction and change detection in multitemporal remotely sensed images with multiple birth and death dynamics[END_REF], [START_REF]Building detection in a single remotely sensed image with a point process of rectangles[END_REF], and further demonstrating figures and experimental results are provided in [START_REF]Building extraction and change detection in multitemporal aerial and satellite images in a joint stochastic approach[END_REF]. Due to its modularity, the method could be easily adapted to different object level change detection applications, for instance tree or road detection. On the other hand, we attempt to focus on the task specific issues as well. We present a broad feature library, which can be appropriate for the detection of a large set of buildings, expecting various image properties. For this reason, in the following section we provide an overview on the state-of-the art methods for monocular building extraction.

Related Works in Monocular Building Detection

A SIFT key point based method has been presented in [START_REF] Sirmac | Urban-area and building detection using SIFT keypoints and graph theory[END_REF] for urban area extraction and building detection. This technique assumes that the building structures in a given image can be efficiently characterized by a couple of template buildings (here two templates: a bright and a dark one) which are used for training. However, images containing a high variety of buildings may need a huge template library, where the overlap between the building and background domains in the descriptor space may be hard to control. A recent model based on Gabor filters (Gabor) [START_REF]A probabilistic framework to detect buildings in aerial and satellite images[END_REF] represents building positions in the image as joint probability density functions of four different local feature vectors, and performs data and decision fusion in a probabilistic framework to detect building locations. It is important to note that in [START_REF] Sirmac | Urban-area and building detection using SIFT keypoints and graph theory[END_REF] and [START_REF]A probabilistic framework to detect buildings in aerial and satellite images[END_REF] the goal is building localization, but the roof outlines are not extracted, which makes it difficult to apply the method for change detection.

A stochastic MRF framework is introduced in [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF] for detecting building rooftops from single images, which combines 2-D and 3-D information. This approach is based on hierarchical grouping of extracted edge segments to form continuous lines, junctions and finally closed curve hypotheses. However, several restrictions are applied for buildings: it is assumed that they have uniform height, they are composed of planar surfaces with parallel sides, and each building casts its shadow on a locally flat surface. Similarly to [START_REF] Song | Building roof detection from a single high-resolution satellite image in dense urban area[END_REF], [START_REF] Saeedi | Automatic building detection in aerial and satellite images[END_REF] the method needs a reasonable edge map, because missing large side parts and false edges inside and around the buildings may corrupt the edge grouping process. Edge based building detection is also applied in [START_REF] Porway | A hierarchical and contextual model for aerial image parsing[END_REF] as a part of a complete scene interpretation process. This approach deals with different object categories, and implements multi-level interactions within a scene and between different object types as well.

Combining roof color, shadow and edge information has been suggested in [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF] in a two-step process which we refer to later as the Edge Verification (EV) approach. In EV, color and shadow are used first for coarse built-in candidate area estimation; thereafter, fitting the building rectangles and verification of the proposals are based purely on the Canny edge map of the obtained candidate regions. As a drawback, this sequential approach is sensitive to the failure of each individual feature. A corrupted edge image causes unreliable corner detection and edge mask stretching, meanwhile, without shadow and color information, the building search area should be extended for the whole image, increasing the processing time and the appearance of false edge patterns.

Segment-Merge (SM) techniques follow an approach different from edge based methods, as they consider building detection as a region level or image segmentation problem [START_REF] Müller | Robust building detection in aerial images[END_REF], [START_REF] Khoshelham | A split-and-merge technique for automated reconstruction of roof planes[END_REF], [START_REF] Song | A region-based approach to building detection in densely build-up high resolution satellite image[END_REF]. In [START_REF] Song | A region-based approach to building detection in densely build-up high resolution satellite image[END_REF] the authors assume that buildings are homogenous areas w.r.t. either color or texture, which can be used for training-based background subtraction. Hence, elementary constraints for shape and size are used to group the candidate regions into building objects. This method can fail, if the background and building areas cannot be efficiently separated with the chosen color or texture descriptors, thus several building and background parts are merged in the same regions of the oversegmented map. On the other hand, for homogenous buildings (see BEIJING, Fig. 4) or salient roof colors (see BUDAPEST red roofs, Fig. 26, top) region features are often more robust than weak or ragged edge maps.

Beside probabilistic models [START_REF] Lafarge | Structural approach for building reconstruction from a single DSM[END_REF], [START_REF] Ortner | A marked point process of rectangles and segments for automatic analysis of digital elevation models[END_REF], variational techniques [START_REF] Karantzalos | Recognition-driven two-dimensional competing priors toward automatic and accurate building detection[END_REF], [START_REF] Peng | An improved snake model for building detection from urban aerial images[END_REF] have been proposed recently for building extraction through energy minimization. Similarly to our method, the Recognition-Driven Variational (RDV) framework of [START_REF] Karantzalos | Recognition-driven two-dimensional competing priors toward automatic and accurate building detection[END_REF] is based on data and prior term decomposition. However, they focus principally on the prior shape modeling issue and use a simplified image-dependent model part, which assumes that the building and background regions can be roughly separated through considering them as locally homogenous intensity classes. In cases where this data term cannot detect probable building regions, the algorithm naturally fails.

From another point of view, the prior models of [START_REF] Lafarge | Structural approach for building reconstruction from a single DSM[END_REF], [START_REF] Karantzalos | Recognition-driven two-dimensional competing priors toward automatic and accurate building detection[END_REF] contain libraries of complete object shapes, while other approaches [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF], [START_REF] Ortner | A marked point process of rectangles and segments for automatic analysis of digital elevation models[END_REF] construct the objects from elementary building blocks (rectangles or line segments), and the higher level shape information is encoded by interaction constraints of the nearby components. While global description of RDV [START_REF] Karantzalos | Recognition-driven two-dimensional competing priors toward automatic and accurate building detection[END_REF] can be efficient if all objects of the scene can be characterized by a restricted number of prototype shapes, the algorithm fails to detect the boundaries accurately, if a given building cannot be sufficiently represented by any shape from the database, using any possible planar projection. On the other hand, the constructive approach -which we follow in the current paper -is preferable if the prior models of the buildings are partially unknown or largely diverse.

As for image data modeling, the above overviewed methods are based on image-or scene-specific hypotheses, such as unique roof colors [START_REF] Sirmac | Urban-area and building detection using SIFT keypoints and graph theory[END_REF], shadows [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF], [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF], strong edges [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF], [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF], [START_REF] Song | Building roof detection from a single high-resolution satellite image in dense urban area[END_REF], [START_REF] Saeedi | Automatic building detection in aerial and satellite images[END_REF], [START_REF] Porway | A hierarchical and contextual model for aerial image parsing[END_REF], homogeneous roofs [START_REF] Müller | Robust building detection in aerial images[END_REF], [START_REF] Khoshelham | A split-and-merge technique for automated reconstruction of roof planes[END_REF], [START_REF] Song | A region-based approach to building detection in densely build-up high resolution satellite image[END_REF], or a limited number of 2-D [START_REF] Sirmac | Urban-area and building detection using SIFT keypoints and graph theory[END_REF] or 3-D [START_REF] Lafarge | Structural approach for building reconstruction from a single DSM[END_REF] building templates. The obvious limitations of these techniques come from the nature of the varying image data, and the lack of adaptivity to different circumstances. To develop more generic models, besides the extraction of the descriptors, feature integration

PROBLEM DEFINITION

The input of the proposed method consists of two co-registered aerial or satellite images which were taken from the same area with several months or years of time difference. Thus a single photo is available at each time instance, and we cannot exploit additional meta-information such as maps or topographic building databases. We expect the presence of registration or parallax errors, but we assume that they only cause distortions of a few pixels. We consider each building to be constructed from one or many rectangular building segments, which we aim to extract by the model described in the following. As output we provide the size, position and orientation parameters of the detected building segments, and give information which objects are new, demolished, modified/rebuilt or unchanged.

Let us denote by S the common S W ×S H pixel lattice of the input images and by s ∈ S a single pixel. Let u be a building segment candidate assigned to the input image pair, which is jointly characterized by geometric and temporal attributes. We consider the center of each building, c = [c x , c y ] as a point in [0, S W ]×[0, S H ] ⊂ R 2 , which can be projected to S by simple discretization: c → [⌊c x ⌋, ⌊c y ⌋]. Let the rectangle R u ⊂ S be the set of pixels corresponding to u. Apart from the center, R u is described by the e L , e l side lengths, and θ ∈ [-90 • , +90 • ] orientation parameters as shown in Fig. 1.

For purposes of dealing with multiple time layers, we assign to each u an index flag, ξ(u) ∈ {1, 2, * }, where ' * ' indicates an unchanged object (i.e. present in both images), while '1' and '2' correspond to building segments which appear only in the first or second image respectively. We will denote the set of all the possible object records u=(c x , c y , e L , e l , θ, ξ) by H. The output of the proposed model is a configuration of building segments, ω ∈ H n , where n, the number of objects is also unknown.

The method exploits rough preliminary knowledge about the object sizes, which will be introduced in two steps for easier interpretation. In the first part of the discussion, we assume that the side length parameters of the building segments in the scene have the same order of magnitude and can be constrained by e L (u) ∈ [e min L , e max L ] and e l (u) ∈ [e min l , e max l ]. Later in Section 5 we present a multi-scale extension of the process, which enables us to handle of image inputs which contain buildings with significantly different sizes.

FEATURE SELECTION

In this section, we introduce different image features for building and change recognition. Since the proposed model obtains the optimal object configuration through stochastic birth-death iterations, two essential questions should be answered based on the image data. First, how can we efficiently generate relevant objects during the birth process? Secondly, how can it be ensured that the adequate objects survive the death step? To keep focus on both challenges, we utilize low level and object level features in parallel.

Low level features are extracted around each pixel as typical color, texture and local similarity between the time layers. They are principally used in the birth step, to estimate where the buildings might be located, what they might look like, and where changes should be expected. As a consequence, objects are generated with higher probability in the estimated built-up regions, considering the estimated appearance models.

On the other hand, object level features evaluate a building hypothesis for each proposed oriented rectangle. The choice of preserving or killing an object in the death step strongly depends on object descriptors, thus their accuracy is crucial.

Low level features for building detection

We begin the discussion with low level features extracted from individual images. For the purposes of built-in area estimation, at each pixel s we calculate a pair of birth probabilities, P

b (s) and P

(2) b (s), which give the likelihood of s being an object center in image 1, and 2, respectively. The nomination refers to the fact that in the birth step the frequency of proposing an object at s will be proportional to the local birth probabilities. On the other hand, we also assign expected orientation µ (i) θ (s), and side length values µ

(i) L (s) resp. µ (i)
l (s) to the image pixels, which help in estimating the θ, e L and e l parameters of objects centered at s based on various descriptors from the i th image (i ∈ {1, 2}). Since the calculation of birth, orientation and mean side length maps are the same for both time layers, we simplify the notation by ignoring the image index in the following part of this section. Later on, we will denote the time stamp again by a superscript index in parentheses wherever necessary.

Local Gradient Orientation Density

The first feature exploits the fact that regions of buildings should contain edges in perpendicular directions. This property can be robustly characterized by local Gradient Orientation Density Functions (GODF) [START_REF] Kumar | Detection in natural images using a causal multiscale random field[END_REF]. Let ∇g s be the intensity gradient vector at pixel s with magnitude ||∇g s || and angle ϑ ∇ s . Let W l (s) be the rectangular l × l sized window around s, where l is chosen as W l (s) can cover an average building from the training set narrowly. For each s we calculate the weighted ϑ ∇ s density of W l (s):

λ s (ϑ) = 1 N s r∈W l (s) 1 h • ||∇g r || • k ϑ -ϑ ∇ r h (1)
where N s = r∈W l (s) ||∇g r ||, and k(.) is a kernel function with a bandwidth parameter h. We use uniform kernels for quick calculation. If W l (s) covers a building, the λ s (ϑ) function has two peaks, located at a distance of 90 • from each other in the ϑ-domain (see Fig. 2). This property can be measured by correlating λ s (ϑ) with an appropriately matched bi-modal density function:

α(s, m) = λ s (ϑ)η 2 (ϑ, m, d λ ) dϑ (2)
where η 2 (.) is a mixture of two Gaussians with mean values m, resp. m + 90 • , and deviation d λ for both components (d λ is a parameter of the process set by training). Offset m s and value α s of the maximal correlation can be obtained as:

m s = argmax m∈[-90 • ,0] α(s, m) α s = α s, m s
Pixels with high α s are more likely to be centers of buildings, which can be encoded in a gradient-based birth map P gr b (s) = α s / r∈S α r . For the sample image in Fig. 3, a thresholded P gr b map is shown in Fig. 3(b). Furthermore, let us observe that offsets m s and m s + 90 • estimate the dominant gradient directions in the W l (s) region. Thus, for a building with center s, we expect its θ parameter around a mean orientation value µ θ (s), defined as:

µ θ (s) = m s if λ s (m s ) > λ s (m s + 90 • ) m s + 90 • otherwise (3)
For this reason if the birth step proposes an object u at pixel s, its orientation is set as θ(u) = µ θ (s) + η θ , where η θ is a random value, generated for each object independently according to a zero-mean Gaussian distribution with a small deviation parameter σ θ .

Roof color filtering and shadow evidence

Several types of roofs can be identified by their typical colors [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF]. Let us assume that based on a roof color hypothesis, we extract an indicator mask ̺ co (s) ∈ {0, 1} (e.g. by thresholding a chrominance channel), where ̺ co (s) = 1 marks that s has roof color. Many roof pixels are expected around building centers, thus for each s we calculate the accumulated ̺ co -filling factor in its neighborhood:

Γ s = r∈W l (s) ̺ co (r).
The color birth map value is obtained as A supplementary evidence for the presence of buildings can be obtained through their cast shadows [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF], [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF]. In several types of remote sensing scenes, a binary shadow mask ̺ sh (s) can be derived by filtering pixels from the darkblue color domain [START_REF] Tsai | A comparative study on shadow compensation of color aerial images in invariant color models[END_REF]. The relative alignment of shadows to the buildings is determined by the global Sun direction, which can be set with minor user interaction or calculated automatically [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF]. Consequently, we can identify the building candidate areas as image regions lying next to the shadow blobs opposing the Sun direction (see Fig. 3(d) and later Fig. 10). As for the shadow based birth map, we use a constant birth rate P sh b (s) = p sh 0 within the obtained candidate regions and a significantly smaller constant on the outside. It is also important to note that for building detection only the cast shadows (i.e. shadows on the ground) are relevant, while self shadows (i.e. weakly or not illuminated building parts) should be ignored. However, as pointed out in [START_REF] Dare | Shadow analysis in high-resolution satellite imagery of urban areas[END_REF], in most cases cast and self shadows have different intensity values, since the shadowed object parts are mostly illuminated by secondary light sources such as reflections from surrounding buildings.

P co b (s) = Γ s / r∈S Γ r . Note

Roof homogeneity

As illustrated in Fig. 3, the P gr b (s) and P sh b (s) birth maps usually give a quite coarse estimation of the built-up regions, which is hardly appropriate for building separation and size estimation. Although we may obtain notably accurate footprints through roof color filtering (Fig. 3(c)), it can only be used for a limited subset of the images and objects. On the other hand, in high resolution images provided by satellites such as Ikonos and Quickbird, a significant part of the roof tops can be identified as homogeneous blobs in the coarsely detected building candidate regions. In this section we investigate how roof homogeneity can be exploited for building region detection and refinement.

The feature extraction algorithm consists of the following steps (illustration for the BEIJING image is shown in Fig. 4):

• Candidate Region Filtering: for a given input image (Fig. 4(a)) obtain the coarse preliminary building candidate (PBC) regions based on the gradient and/or shadow features, as explained in Sec.

P ho b (s) = k R ||s -c(R min s )|| h R (4)
with a k R (.) kernel function, and h R bandwidth parameter. Besides marking the candidate regions of the building centers, the {R i |i = 1 . . . t} set provides local estimations for the side length parameters: µ L (s) = e L (R min s ) and µ l (s) = e l (R min s ). Of course, we can only assume this information to be reliable in pixel positions with high homogeneity birth factors. Thus, for an object u proposed at s, we set the side length values with a probability proportional to P ho b (s) as:

e L (u) = µ L (s) + η L (s), e l (u) = µ l (s) + η l (s)
where η L (s) and η l (s) are independent zero mean Gaussian random variables. Note that side length estimates can be similarly extracted from the color feature map. This preliminary calculation is particularly significant if the object sizes show a large variety, since sampling the side length parameters of the proposed objects according to a prior distribution with a wide support can slow down the speed of the iterative birth-death process critically.

Low level change feature

Up to this point, we have used various descriptors to estimate the location and appearance of the buildings in the individual images. However, the gradient orientation statistics also offers a tool for low level region comparison, which can be directly involved in the scenario model. Let us consider the λ s (.) orientation density introduced in Sec. 3.1.1. Matching the λ

s (.) and λ

s (.) functions from the two time layers can be interpreted as low level similarity checking of the areas around s in the two images, based on "building-focused" textural features (see Fig 5), which are independent of illumination and coloring effects and robust regarding parallax and registration errors. For measuring the dissimilarities, we use the Bhattacharyya distance:

b(s) = -log λ (1) s (ϑ) • λ (2) s (ϑ)dϑ (5)
Choosing an appropriate b 0 threshold [START_REF] Bruzzone | An adaptive semiparametric and context-based approach to unsupervised change detection in multitemporal remote-sensing images[END_REF], the binarized pixel level change mask is obtained as: As shown in Fig. 6, the above comparison separates efficiently the image regions which contain the changed and unchanged buildings, respectively. Knowing that l 2 is the area of window W l (s), the probability of change around pixel s is derived as:

̺ ch (s) = 1 if b(s) > b 0 0 if b(s) ≤ b 0 (6)
P ch (s) = r∈W l (s) ̺ ch (r)/l 2 (7)
Considering the change feature, we can exploit an additional information source for scene interpretation, which is independent of the object recognizer.

Integration of the different birth maps

Since the main goal of the combined birth map in each image is to keep focus on all building candidate areas, we derive it with the maximum operator from the birth maps of the features. For example, when gradient, color and shadow are simultaneously used, we obtain the final field as P b (s) = max P gr b (s), P co b (s), P sh b (s) ∀s ∈ S. For input, without shadow or color information, we can ignore the corresponding feature in a straightforward way, or exchange the P co b (s) component to the homogeneity birth value, P ho b (s). In the birth step of the bMBD process, the birth maps of both time layers, P

Object-Level Features

Besides efficient object generation, the second key point of the applied birth-death dynamics based approach is to validate the proposed building segment candidates. In this section, we construct a ϕ (i) (u) : H → [-1, 1] energy function, which calculates a negative building log-likelihood value of object u in the i th image (hereafter we ignore the i superscript). By definition, a rectangle with ϕ(u) < 0 is called attractive object, and we aim to construct the ϕ(u) function so that attractive objects correspond exclusively to the true buildings.

The process consists of three parts: feature extraction, energy calculation and feature integration. First, we define different f (u) : H → R features which evaluate a building hypothesis for u in the image, so that 'high' f (u) values correspond to efficient building candidates. In the second step, we construct energy subterms for each feature f , by attempting to satisfy ϕ f (u) < 0 for real objects and ϕ f (u) > 0 for false candidates. For this purpose, we project the feature domain to [-1, 1] with a monotonously decreasing function shown in Fig. 8:

ϕ f (u) = Q f (u), d f 0 , D f where Q(x, d 0 , D) = 1 -x d0 , if x < d 0 exp -x-d0 D -1, if x ≥ d 0 (8)
Observe that the Q function has two parameters: d 0 and D. While D f performs data-normalization, d f 0 is the object acceptance threshold concerning feature f : u is attractive according to the ϕ f (u) term iff f (u) > d f 0 . Finally, we must consider, that the decision based on a single feature f can lead to a weak classification, since the buildings and the background may overlap in the f -domain. Therefore, in the third step (Sec. 3.4.2), the joint energy term ϕ(u) must be appropriately constructed from the different ϕ f (u) feature modules.

Feature Models

We begin with gradient analysis. Below the edges of a relevant rectangle candidate R u , we expect the magnitudes of the local gradient vectors (∇g s ) to be high and the orientations to be close to the normal vector (n s ) of the closest rectangle side (Fig. 7). The f gr (u) feature is calculated as:

f gr (u) = 1 # ∂R u s∈ ∂Ru ∇g s • n s (9)
where '•' denotes scalar product, ∂R u is the dilated edge mask of rectangle R u , and # ∂R u is the number of pixels in ∂R u .

The dilation of the R u mask outline is necessary to tolerate slightly imperfect edge alignment and minor registration errors between the images. The data-energy term is calculated as:

ϕ gr (u) = Q(f gr (u), d gr , D gr ).
The calculation of the roof color feature is shown in Fig. 9. We expect the image points to have dominantly roof colors inside the building footprint R u , while the T u object-neighborhood (see Fig. 9) should contain a majority of background pixels. Hence we calculate the internal f co in (u) and external f co ex (u) filling factors, respectively, as: Here #X denotes the area of X in pixels and ̺ co (s) is the color mask value by s. We prescribe that u should be attractive according to the color term if it is attractive both regarding the internal and external subterms. Thus the color energy term is obtained as:

f co in (u) = 1 #R u s∈Ru ̺ co (s); f co ex (u) = 1 #T u s∈Tu 1-̺ co (s)
ϕ co (u) = max [Q(f co in (u), d co in , D co in ), Q(f co ex (u), d co ex , D co ex )]
We continue with the description of the shadow term. This step is based on the binary shadow mask ̺ sh (s), extracted in Sec. 3.1.2. Using the shadow direction vector v sh (opposite of the Sun direction vector) we identify the two sides, AB and BC, of the rectangle R u which are supposed to border on cast shadows, where A, B and C denote the corresponding vertices as shown in Fig. 10. (Note that if v sh is parallel to one of the rectangle sides, we have only one shadow-object edge). Then, we check the presence of shadows in parallelograms (A, A + ε sh , B + ε sh , B) and (B, B + ε sh , C + ε sh , C). Here ε sh is a scalar so that ||ε sh • v sh || approximates the shadow width of the shortest buildings in the scene. The union of the two parallelograms forms the T sh u shadow candidate region as shown in Fig. 10. Thereafter, similarly to the color feature, expect low shadow presence f sh ex (u) in the R u internal and a high one f sh ex (u) in the T sh u external region:

f sh in (u) = 1 #R u s∈Ru 1-̺ sh (s) ; f sh ex (u) = 1 #T sh u s∈T sh u ̺ sh (s)
As for the energy term:

ϕ sh (u) = max Q(f sh in (u), d sh in , D sh in ), Q(f sh ex (u), d sh ex , D sh ex )
Note that this approach does not require accurate building height information, since we do not penalize it, if shadow blobs of long buildings exceed the T sh u regions. The roof homogeneity feature can also be exploited at object level. Fig. 11 shows an example of how to describe two-sided roofs. After extracting the symmetry axis of the object candidate u, we can characterize the peakiness of the dark (d) and bright (b) side histograms by calculating their kurtosis f ho d (u), and f ho b (u), respectively. Denoting by g s the gray value of pixel s, and by R d u and R b u the dark and bright regions of R u object rectangle, we get:

f ho d (u) = R d u g 4 s R d u g 2 s 2 ; f ho b (u) = R b u g 4 s R b u g 2 s 2 (10)
If the roof parts are homogeneous, the f ho d (u) and f ho b (u) kurtosis values should be high. However, as shown in Fig. 12, the homogeneity feature may have false maxima for incomplete roofs, since parts of a homogeneous roof are homogeneous as well. Therefore we characterize roof completeness in the following way. We derive the F u floodfill mask of u, which contains the pixels reached by floodfill propagations from the internal points of R u . If the homogeneous roof is complete, F u must have low intersection with the NH u , resp. NV u , 'horizontal', and 'vertical', neighborhood regions of R u (see Fig. 12). Finally, the ϕ ho (u) energy term can be constructed from the kurtosis and completeness descriptors in a manner to the previous attributes.

Feature integration

Usually, the individual features are in themselves inappropriate for modeling complex scenes, which is illustrated in Fig. 13. For the Ground Truth (GT) buildings (see Fig. 13(b)), we can follow here the effectiveness of the gradient (Fig. 13(c),(d)), shadow (Fig. 13(e),(f)) and color (Fig. 13(g),(h))) descriptors, respectively. The gradient and shadow maps are considerably noisy in the right upper image part, however, the roofs can be detected here fairly by extracting image regions with high a* color component values in CIE L*a*b* color space representation. Conversely, 'non-red' buildings in the bottom-left regions can be efficiently detected by edge and shadow features.

To answer the challenges of such object or data heterogeneity problems, the proposed framework enables flexible feature integration depending on the available image inputs. From the feature primitive terms introduced in Sec. 3.4, first we construct building prototypes. For each prototype we can prescribe the fulfillment of one or many feature constraints whose ϕ f -subterms are connected with the max operator in the joint energy term of the prototype (logical AND in the negative log-likehood domain).

Additionally, several building prototypes can be detected simultaneously in a given image pair, if the prototype-energies are joined with the min (logical OR) operator. Thus the final object energy term is derived by a logical function, which expresses some prior knowledge about the image and the scene, and it is chosen on a case-by-case basis. For example, in the BUDAPEST pair we use two prototypes: the first one prescribes the edge and shadow constraints, the second one the roof color, thus the joint energy is calculated as:

ϕ(u) = min max {ϕ gr (u), ϕ sh (u)}, ϕ co (u) . (11)
Similarly, for the BEIJING images (see Fig. 26, bottom) we use gradient (ϕ gr ) & shadow (ϕ sh ) and homogeneity (ϕ ho ) & shadow (ϕ sh ) prototypes.

CONFIGURATION MODEL AND OPTIMIZATION

In this section we transform the building change detection task into an energy minimization problem. Following our definitions from Sec. 2, the u building segment candidates (i.e. objects) live in a bounded parameter space H. Since we aim to extract building populations from the images, we need to propose a configuration space Ω, which is able to deal with an unknown number of objects:

Ω = ∞ n=0 Ω n , Ω n = {u 1 , . . . , u n } ⊂ H n (12)
Hereafter we will use the notation ω ∈ Ω for an arbitrary object configuration, thus ω = ∅, or ω = {u 1 , . . . , u n } for an n ∈ {1, 2, . . .} and u i ∈ H : ∀i ∈ {1, 2, . . . , n}.

Configuration Energy

The Marked Point Process framework enables to characterize whole populations instead of individual objects, through exploiting information from entity interactions. Following the classical Markovian approach, each object may only affect its neighbours directly. This property limits the number of interactions in the population and results in a compact description of the global scene, which can be analyzed efficiently. To realize the Markov-property, one should define first a ∼ neighborhood relation between the objects in H. In our model, we say that u ∼ v if their rectangles R u and R v intersect.

Let us denote by D the union of all image features derived from the input data. For characterizing a given ω object population considering D, we introduce a non-homogenous data-dependent Gibbs distribution on the configuration space:

P D (ω) = 1 Z • exp -Φ D (ω) (13)
with a Z normalizing constant: Z = ω∈Ω exp -Φ D (ω) , and Φ D (ω) configuration energy:

Φ D (ω) = u∈ω A D (u) + γ • u,v∈ω u∼v I(u, v) (14)
Here A D (u) ∈ [-1, 1] and I (u, v) ∈ [0, 1] are the data dependent unary and the prior interaction potentials, respectively, and γ > 0 is a weighting factor between the two energy terms. Thus the Maximum Likelihood (ML) configuration estimate according to P D (ω) can be calculated as ω ML = argmin ω∈Ω Φ D (ω) . Unary potentials characterize a given building segment candidate u = {c x , c y , e L , e l , θ, ξ} as a function of the local image data in both images, but independently of other objects of the population. This term encapsulates the building energies ϕ (1) (u) and ϕ (2) (u) extracted from the 1 st , resp. 2 nd , image (Sec. 3.4) and the low level similarity information between the two time layers which is described by the ̺ ch (.) change mask (Sec. 3.2).

We remind the reader that our approach marks each building segment u with an image index flag from the set {1, 2, * }, depending on that u appears in one [ξ(u) ∈ {1, 2}] or both [ξ(u) = * ] of the input images. In this way, the classification of the building segment u is straightforward: u is unchanged iff ξ(u) = * ; new iff ξ(u) = 2 and ∄v ∈ ω : {ξ(v) = 1, u and v overlap}; and demolished iff ξ(u) = 1 and ∄v ∈ ω : {ξ(v) = 2, u and v overlap}. Modified buildings are considered as two objects u 1 and u 2 , so that ξ(u 1 ) = 1, ξ(u 2 ) = 2.

The following soft constraints are considered by the potential terms in the various cases:

• unchanged building u: we expect low object energies in both images, and penalize textural differences (i.e. pixels with ̺ ch (s) = 1) under its footprint R u . • demolished or modified building in the first image: we expect low ϕ (1) (u), and ϕ (2) (u) is indifferent. We penalize high similarity under the footprint. • new or modified building in the second image: we expect low ϕ (2) (u), and ϕ (1) (u) is indifferent. We penalize high similarity under the footprint.

A D (u) = I [ξ(u)∈{1, * }] • ϕ (1) (u) + I [ξ(u)∈{2, * }] • ϕ (2) (u)+ + I [ξ(u)= * ] • 1 #R u s∈Ru ̺ ch (s)+ + I [ξ(u)∈{1,2}] • 1 #R u s∈Ru 1 -̺ ch (s) (15)
On the other hand, interaction potentials realize prior geometrical constraints: they penalize intersection between different object rectangles sharing the time layer (see Fig. 14):

I(u, v) = I [ξ(u)≃ξ(v)] • #(R u ∩ R v ) #(R u ∪ R v ) (16)
where

ξ(u) ≃ ξ(v) relation holds iff ξ(u) = ξ(v), or ξ(u) = * , or ξ(v) = * .
Since ∀u, v : I(u, v) ≥ 0, the optimal population should exclusively consist of objects with negative data terms (i.e. attractive objects): if A D (u) > 0, removing u from the configuration results in a lower Φ D (ω) global energy [START_REF] Bruzzone | An adaptive semiparametric and context-based approach to unsupervised change detection in multitemporal remote-sensing images[END_REF]. Note also that according to eq. ( 14), the interaction term plays a crucial role by penalizing multiple attractive objects in the same or strongly overlapping positions. Note that in the introduced probabilistic model, it is also possible to involve additional prior knowledge about the layout of settlements, by adding further prior terms to the global energy function Φ D (ω). For example, in a town, buildings are usually aligned, hence, we can use the geometric interaction terms of [START_REF] Ortner | A marked point process of rectangles and segments for automatic analysis of digital elevation models[END_REF], where the alignment constraint favors small angle difference and low distance between appropriately matched corners of neighboring segments, while the paving constraint favors parallel rectangles that are located side by side inducing clean arrangements of buildings.

Bi-layer Multiple Birth and Death Optimization

By fixing the A D (u) and I(u, v) potential terms, the Φ D (ω) configuration energy is completely defined, and the optimal ω ML building population can be obtained by minimizing eq. ( 14). For this purpose, we have developed the bi-layer Multiple Birth and Death (bMBD) algorithm, the main steps can be followed in Fig. 15. The bMBD method extends the conventional MBD technique by handling two time layers, thus it encapsulates change and object information simultaneously. Pairs of consecutive birth and death processes are iterated until convergence is obtained in the global configuration. In the birth step, multiple object candidates are generated randomly according to the birth maps P (i) b (s), and as a further novelty, also considering the change probabilities P ch (s) with the expected parameter maps µ

(i) θ (s), µ (i) L (s) and µ (i) l (s) i ∈ {1, 2}.
The death process attempts to eliminate the inappropriate objects based on the global configuration energy.

MULTI-SCALE GENERALIZATION

In Sec. 2, we have used a single size hypothesis for all buildings in the image. However, in scenes where the sizes

Bi-layer Multiple Birth and Death (bMBD) algorithm

1) Initialization: calculate the

P (i) b (s), P ch (s), µ (i) θ (s), µ (i) L (s) and µ (i) l (s) (i ∈ {1, 2}
) birth maps, and start with an empty population ω = ∅.

2) Main program: initialize the inverse temperature parameter β = β0 and the discretization step δ = δ0 and alternate birth and death steps:

• Birth step: for each pixel s ∈ S, if there is no object with center s in the current configuration ω, pick up ξ ∈ {1, 2, * } randomly, let be For each object u taken in this order, compute ∆Φω(u) = ΦD(ω/{u}) -ΦD(ω), derive the death rate as follows:

P b = P ch (s) • P (ξ) b (s) if ξ ∈ {1, 2} 1 -P ch (s) • max {P
dω(u) = δaω(u) 1 + δaω(u) , with aω(u) = e -β•∆Φω (u)
and remove u from ω with probability dω(u). Note that according to eq. ( 14), ∆Φω(u) depends only on u and its neighbours in ω, thus dω(u) can be calculated locally without computing the global configuration energies ΦD(ω/{u}) and ΦD(ω). • Convergence test: if the process has not converged, increase the inverse temperature β and decrease the discretization step δ by a geometric scheme and go back to the birth step. Convergence is obtained when all the objects added during the birth step, and only these ones, have been killed during the death step.

SZADA Aerial F ÖMI ‡ √ √ × √ × C. D'AZUR Satellite Google Earth × √ √ √ × BODENSEE Satellite Google Earth × √ √ √ × NORMANDY Satellite Google Earth × √ √ √ × MANCHESTER Satellite Google Earth × √ √ √
× †indicate if multiple time layers are available for change detection ‡Hungarian Inst. of Geodesy, Cartography and Remote Sensing of the footprints are notably diverse (see Fig. 17 and18), this simplification may prove to be inefficient. Considering multiple scales concerns the low level feature extraction part and the Birth step of the bMBD process particularly (see Fig. 15), since the object level features (Sec. 3.4) calculated in the Death step are normalized either with the object area or with the perimeter. In this section, we demonstrate a generalization of the method through providing a multi-scale extension of the Gradient Orientation Density Function (GODF), λ s (.), and introducing further modifications which are necessary in the proposed bMBD algorithm. We also note that most of the other low level features can be handled in a similar manner.

Let us assume that we have J different building size hypothesizes: ∀u : u ∈ ∪ J j=1 Υ j where u ∈ Υ j iff e L (u) ∈ [e min L,j , e max L,j ] and e l (u) ∈ [e min l,j , e max l,j ]. We remind the reader that the λ s (.) feature has been calculated over a rectangular image region W l (s) around pixel s, where the l window side length has been set according to the estimated average object size. In the multi-scale extension, we calculate the local GODF for J different window sizes (l = l 1 , . . . , l J ) corresponding to the J size hypotheses. Fig. 16 demonstrates this process with J = 3, l 1 = 14, l 2 = 24 and l 3 = 40 parameter settings. At each scale j = 1 . . . J, we calculate the P gr b (s, j) birth probabilities and µ θ (s, j) mean orientation estimates separately, and get the final gradientbirth-map as P gr b (s) = max j P gr b (s, j). A minor modification should also be inserted into the birth map summarization step: we store at each pixel s the Fig. 17. Detection results in case of significantly different building scales, using the features from Fig. 16 dominant feature term, χ(s) = argmax χ∈{gr,co,sh,ho} P χ b (s). This indicator is utilized by the modified Birth process. If we decide to generate a new object u at pixel s based on P b (s), we select its scale j(u) randomly, so that the probability of choosing scale j is P

χ(s) b (s, j)/ J i=1 P χ(s) b
(s, i). Then we set the orientation θ(u) following the η(., µ θ (s, j(u)), σ θ ) distribution, and the side lengths e L (u) and e l (u) according to uniform distributions around the expected length values at scale j(u).

PARAMETER SETTINGS

We can divide the parameters of the proposed mMPP method into three groups corresponding to the prior model, data model and the bMBD optimization.

The prior model parameters, such as the number of the examined scales (J in Sec. 5), l (or l j ) window sizes for GODF calculation (see Sec. 3.1.1) and maximal/minimal rectangle side lengths at the difference scales, depend on image Edge Verification method [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF] • Building Candidate Region (BCR) extraction with a morphological approach -shadow and roof color information exploited solely in preprocessing • Canny edge detection inside the BCRs • Roof corner estimation by detection of perpendicular edgejunctions in the BCRs • Rectangle fitting for the BCR-edge map around each corner candidate • Hypothesis acceptance/rejection Segment-Merge method [START_REF] Müller | Robust building detection in aerial images[END_REF] • Building segment estimation by seeded region growing • Region merging and shadow evidence verification • Filtering based on geometric and photometric features • Polygon approximation of the building blocks Fig. 19. Main steps of the Edge Verification [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF] and Segment-Merge [START_REF] Müller | Robust building detection in aerial images[END_REF] methods used for comparison resolution and expected object dimensions. They are set based on sample objects. We used a constant γ = 2 weight between the data term and the overlapping coefficient in (eq. [START_REF] Bruzzone | An adaptive semiparametric and context-based approach to unsupervised change detection in multitemporal remote-sensing images[END_REF].

The parameters of the data model are estimated based on training image regions containing Ground Truth building segments {u gt 1 , u gt 2 , . . . , u gt n }. Consider an arbitrary f (u) feature from the feature library (e.g. f gr (u) gradient descriptor). We remind the reader that each f (u) of our model is a noisy quality measure and the corresponding energy term is obtained as

ϕ f (u) = Q(f (u), d f 0 , D f ) (see Sec. 3.4).
Here we set the normalizing constant as D f = max j f (u gt j )min j f (u gt j ). Exploiting that the Q transfer function is monotonously decreasing with a sole root

f (u) = d f 0 , object u is attractive in image i (i.e. ϕ (i) f (u) < 0) iff f (u) > d f 0 .
Consequently, increasing d f 0 may decrease the false alarm rate and increase the missing alarms corresponding to the selected feature. Since in the proposed model we can simultaneously utilize several building prototypes, our strategy for setting d f 0 is to minimize the false alarms for each prototype, and eliminate the missing buildings using further feature tuples.

Finally, regarding the relaxation parameters, we followed the guidelines provided in [START_REF] Descombes | Object extraction using a stochastic birth-and-death dynamics in continuum[END_REF], and used δ 0 = 20000, β 0 = 50, and geometric cooling factors 1/0.96. , Gabor [START_REF]A probabilistic framework to detect buildings in aerial and satellite images[END_REF], MRF [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF], Edge Verification (EV) [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF], Segment-Merge (SM) [START_REF] Müller | Robust building detection in aerial images[END_REF] methods, and to the Ground Truth. Circles denote completely missing or false objects. SIFT and Gabor only extract building centers.

EXPERIMENTS

The goal of this section is to validate the three key developments of the paper and compare them to the state of the art: (i) the proposed multiple feature based building appearance model, (ii) the joint object-change modeling framework and (iii) the non-homogeneous object birth process based on low level features.

We have evaluated our method using eight significantly different data sets whose main properties are summarized Table 1. Four image collections contain multitemporal aerial or satellite photos from the monitored regions, which enables testing both the building extraction and the change detection abilities of the proposed mMPP model. The remaining four data sets contain standalone satellite images acquired from Google Earth, which are only exploited in the evaluation of the building appearance model (Sec. 7.1). To guarantee the heterogeneity of the test sets, we have chosen completely different geographical regions as listed in Table 1. We collected samples from densely populated suburban areas, and built a manually annotated database for validation. For parameter settings, we have chosen in each data set 2-8 buildings (≈ 5%) as training data, while the remaining Ground Truth labels have only been used to validate the detection results. Qualitative results are shown in Fig. 17 We perform quantitative evaluation both at object and pixel levels. On one hand, we measure how many buildings are recognized or classified incorrectly in the different test sets, by counting the missing and falsely detected objects (MO and

Building Segment Description

Although the proposed model handles multiple time layers simultaneously, the building description module introduced in Sec. 3.4 works on single image inputs (birth maps and object energies are calculated in the two images independently). In this subsection, we evaluate solely the object recognition part, therefore we use temporarily a simplified data term A D (u) = ϕ (2) (u), i.e. we detect buildings only in the second image independently of the first one.

In this section, we present numerical and qualitative comparison results versus single-view building detection techniques from the state-of-the-art, which were briefly introduced in Sec. 1.4. On one hand, we have evaluated three recent methods in collaboration with their authors: MRF [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF], SIFT [START_REF] Sirmac | Urban-area and building detection using SIFT keypoints and graph theory[END_REF], and Gabor [START_REF]A probabilistic framework to detect buildings in aerial and satellite images[END_REF]. Source codes of SIFT and Gabor have been provided us for the experiments, while the authors of MRF have tested their model with our image sets. In addition, we have also implemented two methodologically orthogonal methods: the Edge Verification (EV) technique [START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF] and the Segment-Merge (SM) model [START_REF] Müller | Robust building detection in aerial images[END_REF], their main steps are listed in Fig 19. Since EV and SM use similar image features (gradient, shadow, color, homogeneity) to our framework, by considering them in the comparison, we can focus purely on validating the model structures instead of special inputdependent descriptors. Sample output images of the reference methods can be found in Fig. 20 and21.

Quantitative evaluation on the database was performed with SIFT, Gabor, EV, SM and the proposed MPP models, results are shown in Table 2. Since SIFT and Gabor extract the building centers instead of estimating the outline, they are only involved in the object level comparison. Numerical results confirm that the proposed model surpasses all references with 10-26% at object level and with 5-18% at pixel level. Table 3 lists the computational time requirements of the test images with the different methods. From this viewpoint, the Gabor technique is dominantly the most efficient, since its Matlab version outperforms the other C ++ coded methods. On the other hand we can observe that the proposed MPP model is competitive with most reference techniques regarding the average running time as well. Note that as Table 3 demonstrates, the computational complexity for the different images depends in parallel on various factors, such as image size, dominance of the color map, and diversity of the building side length values.

Apart from the low complexity, a significant advantage of the Gabor model [START_REF]A probabilistic framework to detect buildings in aerial and satellite images[END_REF] is that it can deal with various images by changing only a single scale parameter. Conversely, our method uses 1-2 free parameters for each feature, thus it is more dependent on the training set. However, while the Gabor algorithm is successful in detecting compact shaped buildings, it faces difficulties with long segments (see BEI-JING). Some additional problems appear considering dark buildings (BODENSEE), where the gradient directions point away from the building center [START_REF]A probabilistic framework to detect buildings in aerial and satellite images[END_REF]. Both the Gabor and SIFT methods have better performance on panchromatic satellite images, while in aerial photos false positives appear due to many redundant local features extracted in the background. In general, buildings in rural regions are efficiently detected with these models, but in densely populated areas, the false alarm rate increases. On the contrary, inserting various building hypotheses into our MPP framework is straightforward, thus it can efficiently deal both with high contrast satellite images and aerial photos where color is a more dominant feature.

Since our datasets do not contain any metadata about intrinsic and extrinsic camera parameters, several benefits of [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF] cannot be exploited. Most of the observed artifacts of [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF] (see Fig. 20(d)) derive from the sensitivity of the edge detection algorithm, which is alleviated by using a strong assumption: the detected buildings should have a uniform height. Resolution of the images affects strongly the quality of the results. The authors confirmed that the minimal resolution, at which the method is able to operate is only fulfilled in the HR versions of the Budapest and Côte d'Azur images. They reported fundamental problems with the remaining images, partially due to small size, missing color information and low quality of the edges. Buildings with irregular shapes are not detected because they do not fulfill a second assumption, that rooftops are polygons with pairwise parallel sides.

We continue the discussion with the EV and SM reference methods. Both of them follow the deterministic hypothesis generation-acceptance scheme, where buildings ignored by the As in Sec. 1.4, the sequential EV technique is sensitive to the quality of the individual feature maps used in the consecutive algorithmic steps. Similarly to [START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF], the edge mask exploited in object proposition and verification can be strongly corrupted or misleading in cases of low contrast, textured background or large noise (Fig. 20(e)). In addition, examples in Fig. 21(c) and (d) show that curved or scalloped roof edges do not fit appropriately to the straight sides of the model rectangles and result in poor detection on the left part of the image. For the same reasons, we must expect similar artifacts by using other edge-critical methods [START_REF] Song | Building roof detection from a single high-resolution satellite image in dense urban area[END_REF], [START_REF] Saeedi | Automatic building detection in aerial and satellite images[END_REF].

As for the region based approach, Fig. 20(f) and Fig. 21(e) -(f) illustrate that the SM technique fails in detecting roofs which are inhomogeneous both in color and texture. On the other hand, building-like homogeneous blobs may result in false positive objects, while low contrast buildings can be merged with the background and missed during the segmentation process. These problems can also appear using similar methods [START_REF] Khoshelham | A split-and-merge technique for automated reconstruction of roof planes[END_REF], [START_REF] Song | A region-based approach to building detection in densely build-up high resolution satellite image[END_REF].

In summary, the tests confirm that the proposed model surpasses the reference techniques, particularly due to its two key properties: the stochastic object generation process and the parallel utilization of multiple features in the building description module.

Finally, we have tested the sensitivity of the proposed model against the parameters of various feature extraction steps. Fig. 22 shows the pixel level F-scores of detection on the BUDAPEST image, where we perturbated the chrominance threshold (τ cr ) of roof color filtering, the shadow darkness threshold (τ sh ) and the gradient acceptance threshold (d gr ) with maximum ±30% around the optimal value. Results show that the performance varies around 10% in these parameter domains, most significant is the dependence on τ cr .

Joint Object-Change Model

After testing the introduced building detector module in single images, we continue with the validation of the proposed joint object-change classification framework. The mMPP model evaluates a given building segment candidate by simultaneously considering its bi-temporal ϕ (1) (u) and ϕ (2) (u) object energies and the low level change information under the footprint. This approach is compared to the conventional Post Detection Comparison (PDC) [START_REF] Tanathong | Object oriented change detection of buildings after the Indian ocean tsunami disaster[END_REF] technique, where the buildings are separately extracted from the two image layers, and the change information is a posteriori estimated through comparing the location, geometry and spectral characteristics of the detected objects. In the latter case the object-change decision is sequential, thus less information can be exploited by the individual object extraction and change classification steps, respectively. Table 4 confirms, that the PDC method causes more false change alarms than mMPP.

To understand the reasons for the differences between PDC and mMPP, a few illustrative examples are shown in Fig. 25. First, the layer-by-layer detector has missed two object candidates: one in the top of the (a) image (edges are partially hidden by the trees), and one in the bottom-left corner of the (b) image (low contrast). These errors result automatically in false changes by using the PDC approach. However, the joint mMPP model produces appropriate detection results (images (c) and (d)), exploiting that both imperceptible buildings are in certainly unchanged image parts according to the low level ̺ ch -change feature, meanwhile the given objects have been correctly and confidently detected in the other images. On the other hand, false objects appearing in the background regions in PDC have been eliminated by the mMPP model, exploiting that the corresponding local similarity is high again, but the 'twin' object cannot be found at the other time instance.

Feature Based Birth Process

Although the ω ML configuration estimate does not depend on the birth maps, the exploration strategy in the population space affects the speed of optimization notably. In the bMBD algorithm (see Fig. 15) the most significant part of the computational time corresponds to calculating the A D (u) data term for each generated object, since the interaction potential I(u, v) only needs to calculate the intersection area of rectangles which can be solved efficiently in an analytical way. For this reason, the complexity of the approach can be characterized by the number of object birth steps.

The key objective of the proposed non-uniform Feature Based Birth (FBB) procedure is to generate relevant objects with higher probability, so that we need to deal with less inefficient building segment candidates, and high quality configurations can be reached more quickly. We should note here that the P l (s) birth maps are calculated only once before starting the iterative algorithm, and using dynamic programming techniques its computational need is negligible [START_REF]Building extraction and change detection in multitemporal aerial and satellite images in a joint stochastic approach[END_REF] considering the cost of the whole Simulated Annealing (SA) process.

For evaluation, we compared the convergence speed of the bMBD optimization using the proposed FBB and the conventional Uniform Birth (UB) processes. In the UB case, the P (i) b (s) and P ch (s) maps follow a uniform distribution and the side length/orientation parameters are also set as uniform random values. In Fig. 23, the object-errors are shown as a function of the birth steps: the FBB approach reaches the final error rate with 3 times less birth calls than the UB. The difference is even more significant at pixel level. As Fig. 24 shows, with the UB process the pixel level accuracy rates converge much slower than the object errors; to reach the 75% DA rate, we need to generate 400, 000 objects with the UB map, and only 24, 000 building candidates with the proposed FBB map. This observation means that the appearing object silhouettes in the uniform approach are usually notably inaccurate in the beginning, and considerable time is needed to reach the optimum.

CONCLUSION

We have proposed a multitemporal Marked Point Process (mMPP) framework for building extraction and change monitoring in remotely sensed image pairs taken with significant time differences. The method incorporates object recognition and low level change information in a joint probabilistic approach. A global optimization process attempts to find the optimal configuration of buildings, considering the observed

Fig. 1 .

 1 Fig. 1. Definition of the rectangle parameters

Fig. 2 .

 2 Fig. 2. Kernel density estimation of the local gradient orientation histogram around two selected pixels: a building center s and an empty site r.

Fig. 3 .

 3 Fig. 3. Building candidate regions obtained by the low level (b) gradient (c) color and (d) shadow descriptors

  Fig. 4. Preliminary building estimation based on roof homogeneity. Missing and false alarms -denoted by ( * ) in image (f) -are eliminated later in the process.

  3.1.1 and 3.1.2 (Fig. 4(b)-(c)). • Intensity based segmentation: we (over)segment the PBC regions of the input image into homogenous components, and ignore the blobs smaller than 20% of the expected mean building area. This step results in the homogenous building candidate (HBC) region map (Fig. 4(d)). • Orientation based clustering: we re-cluster the HBC map based on the µ θ (s) dominant local gradient orientation values obtained in the regions of interest, and call the result GHBC image as shown in Fig. 4(e). Each uniform component of GHBC is considered in the following as a building segment candidate. • Candidate parameter estimation: we estimate the center and the bounding box (Fig. 4(f)) parameters for each building segment candidate through morphological box fitting techniques. Let us denote the candidate rectangles (Fig. 4(f)) obtained in the previous filtering process by R i , i = 1 . . . t, and let c(R i ) be the center of R i . Then, for each pixel, we determine the closest rectangle R min s = arg min i ||sc(R i )|| and calculate the homogeneity birth value as:

Fig. 5 .

 5 Fig. 5. Comparing the λ(.) functions in the two image layers regarding two selected pixels. s corresponds to an unchanged point and r to a built-up change.

regions of image 2 Fig. 6 .

 26 Fig. 6. Low level change detection: (a) and (b) input images, (c) Bhattacharyya change mask ̺ ch

( 1 )

 1 b (s) and P(2) b (s), and the change map P ch (s) are utilized in parallel. We propose an unchanged object at s with a probability proportional to (1 -P ch (s)) • max i∈{1,2} P (i) b (s), while at the same location, the likelihood of generating a changed building segment is P ch (s) • P (i) b (s) for image i.

Fig. 8 .

 8 Fig. 7. Utility of the gradient feature

Fig. 9 .Fig. 10 .Fig. 11 .

 91011 Fig. 9. Utility of the color roof feature

  Fig. 13. Illustration of the feature maps in the BUDAPEST 2008 image. Gradient and shadow features are relevant in the left-bottom regions, while the color descriptor is efficient in the top-right image parts. In image (d), the gradient feature is shown under the GT object borders and the background color is equal to the average gradient value.

Fig. 14 .

 14 Fig. 14. Calculation of the I(u, v) interaction potentials: intersections of rectangles are denoted by striped areas

b

  (s)} if ξ = * and execute the following birth process with probability δ P b : -generate a new object u with center s and image index ξ -set the eL(u) e l (u) side length parameters as follows:* with a probability P h b (s)/ P b : set the parameters according to η ., µ(ξ)L (s), σL resp. η ., µ (ξ) l (s), σ l Gaussian distributions as explained in Sec. 3.1.3. (Notes: (i) If the homogeneity feature is ignored, P h b (s) is considered as a constant zero map. (ii) If ξ = * , we choose between µ (1) L/l and µ (2)L/l randomly.) * otherwise: set the parameters randomly between prescribed maximal and minimal side lengths, following a uniform distribution -set the orientation θ(u) following the η(., µ (ξ) θ (s), σ θ ) Gaussian distribution as shown in Sec. 3.1 -add u to the current configuration ω • Death step: Consider the configuration of objects ω = {u1, . . . , un} and sort it from the highest to the lowest value of AD(u).

Fig. 15 .

 15 Fig. 15. Pseudo code of the bi-layer Multiple Birth and Death (bMBD) algorithm

Fig. 16 .

 16 Fig. 16. Multi-scale investigations: demonstrating the dependency of gradient orientation histograms on the l window size. Maximal α feature can be obtained with l = 14 for the small building (see on the left) and with l = 40 for the large building (see on the right)

Fig. 18 .

 18 Fig. 18. Detection results in a densely built-in part of the ABIDJAN image set

Fig. 20 .

 20 Fig. 20. Evaluation of the single view building model. Comparing the proposed MPP model to the SIFT[START_REF] Sirmac | Urban-area and building detection using SIFT keypoints and graph theory[END_REF], Gabor[START_REF]A probabilistic framework to detect buildings in aerial and satellite images[END_REF], MRF[START_REF] Katartzis | A stochastic framework for the identification of building rooftops using a single remote sensing image[END_REF], Edge Verification (EV)[START_REF] Sirmac | Building detection from aerial imagery using invariant color features and shadow information[END_REF], Segment-Merge (SM)[START_REF] Müller | Robust building detection in aerial images[END_REF] methods, and to the Ground Truth. Circles denote completely missing or false objects. SIFT and Gabor only extract building centers.

  , 18, 20, 21, 25 and 26.

Fig. 21 .

 21 Fig. 21. Limitations of the EV and SM methods: compared to the proposed model (im. (b)), weak edge map (c) results in weak EV matching (d); while textured buildings on the right do not appear as homogenous blobs in the floodfill map (e), and are ignored by SM detection (f)

3 Fig. 23 .

 323 Fig. 23. Evolution of the detection performance over the iteration steps during optimization on the BUDAPEST image pair. Object level error is given as a function of the applied birth steps.

b

  (s), P ch (s), µ

Fig. 24 .

 24 Fig. 24. Evolution of the detection performance over the iteration steps during the optimization on the BUDAPEST image pair. Detection Accuracy DA (i.e. pixel level Fscore) is given as a function of the applied birth steps.

Fig. 25 .

 25 Fig. 25. Results on ABIDJAN images (source: DGA c France). Left: image from 1996, right: image from 1997. Top: Post Detection Comparison (PDC) (errors are highlighted by circles), Bottom: proposed joint mMPP model

Fig. 26 .

 26 Fig. 26. Results on BUDAPEST (top, image part -provider: Andr ás G ör ög) and BEIJING (bottom, provider: Liama Laboratory CAS, China) image pairs, marking the unchanged (solid rectangles) and changed (dashed) objects

TABLE 1

 1 Main properties of the test data sets, and applicable features from Sec. 3 ( √ =Yes, ×=No).

	Data Set	Type	Source	Ch. det † co	Usable features sh gr

TABLE 2

 2 Numerical object level and pixel level comparison of the SIFT, Gabor, EV, SM and the proposed methods (MPP) on each test data set (best results in each row are typeset by bold.) is ignored from the summarization due to weak performance with most of the methods

	Data Set		SIFT [29]	Object level performance Gabor [30] EV [3] SM [17]	Prop. MPP	EV [3]	Pixel level performance SM [17] Prop. MPP
	Name	#obj. MO	FO	MO	FO	MO	FO	MO	FO	MO	FO	Pr	Rc	Pr	Rc	Pr	Rc
	BUDAPEST	41	20	10	8	17	11	5	9	1	2	4	0.73	0.46	0.84 0.61	0.82	0.71
	ABIDJAN	21	8	5	0	1	2	0	2	1	1	0	0.91	0.73	0.84 0.79	0.83	0.74
	BEIJING	17	7	2	9	8	2	3	4	2	1	0	0.59	0.26	0.71 0.72	0.93	0.71
	SZADA	57	17	26	17	23	10	18	11	5	4	1		0.62	0.79 0.71	0.93	0.75
	C ÔTE D'AZUR	123	55	9	12	24	14	20	20	25	5	4	0.73	0.51	0.75 0.61	0.83	0.69
	BODENSEE	80	34	9	32	8	11	13	18	15	7	6	0.56	0.30	0.59 0.41	0.73	0.51
	NORMANDY	152	69	14	24	14	18	32	30	58	18	1	0.60	0.32	0.62 0.55	0.78	0.60
	MANCHESTER	171	NA NA	53	85	46	17	53	42	19	6	0.64	0.38	0.60 0.56	0.86	0.63
	Overall F-score *		0.663	0.799	0.842	0.798	0.944		0.537		0.668		0.743

*

  On the contrary, our proposed model uses a stochastic birth process, where the P b (s) maps control the frequency of local object propositions only. This approach is more robust, since even in the image regions with erroneously low P b (s) values, the objects are not completely absent; they merely appear later during the bMBD iterations or the silhouettes may be slightly inaccurate.

	F-score	0.75 0.8							Settings: =30 t sh d =22 gr	F-score	0.75 0.8							Settings: =16 t cr d =22 gr
	level	0.65 0.7								level	0.65 0.7
	Pixel	0.6	10	12	14	16	18	20	22	Pixel	0.6	21	24	27	30	33	36	39
	Crominance thresh. of roof color ( ) t cr	Shadow darkness	thresh.	( ) t sh
						F-score	0.7 0.75 0.8						Settings: =16 t cr =30 t sh
						level	0.65
						Pixel	0.6 Gradient accept. threshold ( ) 16 18 20 22 24 26 28 d gr
	Fig. 22. Pixel level performance (F-scoore) in case of
	different parameter settings
	hypothesis generator modules appear automatically as missing
	objects.

TABLE 3

 3 Computational time of the different w.r.t. image sizes (in kPixels)

	Data Set	Size (kPix)	SIFT	Computational time (seconds) Gabor EV SM	MPP
	BUDAPEST †	280	197.3	14.5	120.4	18.4	25.2
	ABIDJAN	148	110.1	7.1	12.0	6.3	10.7
	BEIJING	515	391.1	37.2	155.5	12.9	52.2
	SZADA	1472	200.9	49.8	30.3	89.1	31.5
	C. D'AZUR	723	416.0	57.7	324.3	47.2	68.5
	BODENSEE	536	90.0	27.9	30.4	35.1	66.2
	NORMANDY	1116	236.7	67.2	109.6	72.3	46.3
	MANCHESTER	1073	NA	137.1	132.1	54.0	65.9
	Average Implementation language: 733	234.6 Matlab	37.4 Matlab	111.8 C ++ C ++ 40.2	42.9 C ++
	†Test of [1] with full resolution (1068kPix) needed 45 minutes

ACKNOWLEDGEMENT

The authors acknowledge the test data provided by András Görög (BUDAPEST images), the French Defense Agency (ABIDJAN) and Véronique Prinet from LIAMA Laboratory of CAS Beijing (BEIJING). We are grateful to the authors of the reference methods and their colleagues for their help in the evaluation, especially to Dr. Beril Sirmac ¸ek and Cosmin Mihai. We thank our colleagues at MTA SZTAKI for linguistic corrections and checking the notations of the manuscript.

Her current research interests are in image processing using probabilistic models and variational methods. She also works on parameter estimation and optimization techniques.