2.1.

Extract the Pixel-Level Spectral Value ${\mathsf{F}}^{\mathsf{1}}$

Because buildings usually show a bright appearance in an image and their shadows are dark, a spectral value ${\mathit{F}}^1$ is used to describe this feature. Namely, for a given image $\mathit{Y}=({\mathit{Y}}^1,{\mathit{Y}}^2,\dots ,{\mathit{Y}}^P)$, each spectral channel ${\mathit{Y}}^t$ ($1\le t\le P$) is defined on an $M\times N$ rectangular lattice $\mathit{S}$, i.e., $\mathit{S}=\{s|s=(i,j),1\le i\le M,1\le j\le N\}$ and ${\mathit{Y}}^t={(y_s^t)}_{M\times N}$. Then, spectral value ${\mathit{F}}^1={(f_s^1)}_{M\times N}$ is defined as $f_s^1=\prod _{t=1}^P{y}_s^t$, which can describe the spectral value of each pixel $s$ on different channels.

2.2.

Extract the Region-Level Spectral Variance ${\mathsf{F}}^{\mathsf{2}}$

Different urban objects have various appearances, so their spectral variance should be relatively large. Hence, we design a region-level spectral variance ${\mathit{F}}^2$ to capture this feature. First, the initial objects are obtained using a mean shift method,²⁷ which constructs a probability density to reflect the underlying distribution of points in some feature space and to map each point to the mode of the density which is closest to it. Then, the given image $\mathit{Y}$ is divided into an over-segmented region set $\mathit{R}$, i.e., $\mathit{R}=\{R_1,R_2,\dots ,R_k\}$. Each $R_i$ of $\mathit{R}$ denotes an over-segmented region ($i=\mathrm{1,2},\dots ,k$), $R_i\cap R_j=\varnothing$ ($i\ne j$), and $k$ is the number of these regions. With the region set $\mathit{R}$, we can further define the neighborhood system $N=\{N_i|i=\mathrm{1,2},\dots ,k\}$ to describe the spatial context of regions. Here, each $N_i$ denotes the set of regions neighboring $R_i$. Let $M(R_i)$ be the mean value of pixels in $R_i$, and the local spectral variance between region $R_i$ and its adjacent regions can be calculated as follows:

In Eq. (1), every region has the same impact on $V(R_i)$. Intuitively, it may be preferable to determine the impacts in Eq. (1) using an adaptive way. Hence, the equation for $V(R_i)$ is revised as

In Eq. (2), $M^{*}(R_j,R_i)$ is defined as follows:

Based on the $V(R_i)$, ${\mathit{F}}^2={(f_s^2)}_{M\times N}$ is defined as $f_s^2=V{[R(s)]}^{1/2}$ to reflect the spectral variance among regions. Here, $R(s)$ is the region to which pixel $s$ belongs.

2.3.

Extract the Granularity Information ${\mathsf{F}}^3$

Urban areas have more different types of objects and more complicated appearances than nonurban areas. Therefore, in the over-segmented region set $\mathit{R}=\{R_1,R_2,\dots ,R_k\}$, objects of urban areas usually have smaller region sizes than objects of nonurban areas. In other words, the granularity of urban areas is finer than that of nonurban areas. Hence, we employ the region size and the spatial relationship among regions to define ${\mathit{F}}^3={(f_s^3)}_{M\times N}$, i.e.,

An example to illustrate these features is shown in Fig. 1, where Figs. 1(b), 1(d), and 1(e) are features ${\mathit{F}}^1$, ${\mathit{F}}^2$ and ${\mathit{F}}^3$ extracted from Fig. 1(a). From this example, one can see that the buildings in Fig. 1(b) are bright, which denotes a high ${\mathit{F}}^1$ value, and their shadows are of the low ${\mathit{F}}^1$ value. Similarly, the spectral variance ${\mathit{F}}^2$ of urban areas is larger than that of others areas, and urban areas have a small granularity ${\mathit{F}}^3$ value. Based on these features, we design $E^1={(e_s^1)}_{M\times N}$, $E^2={(e_s^2)}_{M\times N}$, $E^3={(e_s^3)}_{M\times N}$, and $E^4={(e_s^4)}_{M\times N}$ to describe the buildings’ spectral values, dark shadows’ spectral values, regional spectral variance, and granularity information, respectively. They are

Fig. 1

(a) Original aerial image. (b) Pixel-level spectral value. (c) Initial over-segmented region set $\mathit{R}$. (d) Region-level spectral variance. (e) Granularity information. (f) Histogram of ${\mathit{F}}^1$ and $\gamma$ (g) Histogram of ${\mathit{F}}^2$ and $\lambda$. (h) Histogram of ${\mathit{F}}^3$ and $\pi$. (i) Seed points extracted based on (b), (d), and (e).

$\gamma$,$\lambda$, and $\pi$ are the key parameters for the selection of seed points. The parameter $\gamma$ is used to make $E^1$ capture the spectral feature of buildings. Since buildings usually take a high spectral value, they are often expressed as the tail of the histogram of ${\mathit{F}}^1$. Hence, $\gamma$ is set to a high value to get the tail of the histogram of ${\mathit{F}}^1$, such as Fig. 1(f). Correspondingly, $E^2$ uses $1-\gamma$ to obtain the first peak of the histogram of ${\mathit{F}}^1$, which describes the dark shadows with a low ${\mathit{F}}^1$ value. For the same reason, $E^3$ and $E^4$ are set with a high $\lambda$ value and low $\pi$ value to catch the tail of the histogram of ${\mathit{F}}^2$ and the first peak of the histogram of ${\mathit{F}}^3$, respectively. These can extract buildings’ spectral variance and granularity features. An illustration of setting $\gamma$, $\lambda$, and $\pi$ is shown in Figs. 1(f)–1(h).

Then, by sequentially combining $E^1$, $E^2$, $E^3$, and $E^4$, seed points can be obtained. Namely, we first use $D^1={(d_s^1)}_{M\times N}$ to get pixels belonging to buildings and adjoining the shadows, or pixels belonging to shadows and adjoining the buildings. This is defined as

At last, seed points will be selected as the set $D=\{s|d_s^3=1,s\in S\}$.

For $r$ and $l$, these seed points are used to determine whether a local window $w(s,r)$ simultaneously contains pixels from $E^1$, $E^2$, $E^3$, and $E^4$ and whether pixels of each kind are not less than $l$. Because a building is spatially adjacent to its shadow, they can be effectively detected together using a relative small patch of the given image. Hence, by setting $r$ to 2 for $D^1$, $D^2$, and $D^3$, we use the local window $w(s,r=2)$ as the small patch to select seed points in the following. At the same time, if there are buildings and their shadows in the small patch, there will be at least one pixel labeled 1 in the patch for each $e_s^i$, $i=1$, 2, 3, 4. Therefore, $l$ is set to 1. It means that only a pixel which simultaneously possesses or neighbors $E^1$, $E^2$, $E^3$, and $E^4$ within the small local window $w(s,2)$ can be chosen as the seed point. An example is shown in Fig. 1(i). Note that one pixel would show different sizes of the Earth’s surface in remote sensing images with various spatial resolutions, which may affect the setting of parameter $r$. Namely, $r$ can be set to 1 for the low-spatial resolution remote sensing images and be set larger than 2 for extreme high-spatial resolution remote sensing images.

3.1.

MRF Model

Let $X=\{X_{R_i}|R_i\in \mathit{R}\}$ be the label random field defined on the over-segmented region set $R$. We use 1 to flag urban areas and 0 to flag nonurban areas, and each random variable $X_{R_i}$ takes a value of 1 or 0 to represent the label of region $R_i$ it belongs to. If $\mathit{x}=\{x_{R_i}|R_i\in \mathit{R}\}$ denotes the realization of $\mathit{X}$, the optimal realization $\widehat{\mathit{x}}$ can be obtained by maximizing the posterior probability, i.e.,

The energy form of Eq. (4) is

In Eq. (5), the likelihood function $P(\mathit{Y}|\mathit{X})$ is used to describe image features. In this paper, we assume that all $Y_{R_i}$ of $\mathit{Y}$ are independent given labels. That is

The distribution of random field $P(\mathit{X})$ is assumed to be of the Markovianity property, i.e.,

Therefore, Eq. (5) can be rewritten as

Due to the complexity caused by interactions among labels, it is difficult to find the solution of the MRF model. Hence, the local optimal solution $\widehat{\mathit{x}}=({\widehat{x}}_{R_i})$ can be obtained as follows:

3.2.

MRF-Based Region Growing

In this section, an MRF-based region-growing criterion is introduced to find the optimal realization $\widehat{\mathit{x}}$. To minimize the total energy of the MRF model, the proposed method will iteratively merge adjacent regions that could decrease the total energy. Namely, for neighboring regions $R_i$ and $R_t$, the total changed energy $E(R_i,R_t)$ is first calculated these two regions are merged. Based on Eq. (7), $E(R_i,R_t)$ equals the sum of the changed likelihood energy $E_f(R_i,R_t)$ and the changed label energy $E_l(R_i,R_t)$, i.e.,

Here,

Algorithm 1

The proposed criterion is different from traditional region-growing methods, as it does not begin from seed points but from nonseed points. We only consider the nonurban regions labeled 0 and their region sizes are less than the threshold. For each selected region $R_i$, the energy values are calculated between $R_i$ and its neighbor regions, respectively. Then, $R_i$ is merged with the one neighbor region that has the minimum energy value. Hence, $R_i$ merged with an urban region will lead to a larger urban region; in contrast, $R_i$ merged with a nonurban region will result a new nonurban region. Therefore, the rule of our approach is a competition rule of region growing for both urban and nonurban regions.

Urban areas can be extracted using the region-growing criterion. Namely, urban areas are first initialized using the label field $x=\{x_{R_i}|R_i\in R\}$ based on seed points $D$, i.e., set $x_{R_i}=1$ if $R_i\cap D\ne \varnothing$; or else, set $x_{R_i}=0$. Then, by increasing the thresholds, the growing criterion gradually updates the urban areas. Note that different sun angles may affect the shadow length and direction, but it does not change the spatial topological relationship between buildings and their shadows. Hence, the proposed method is robust for effective detection of varying urban areas contained in different remote sensing images.

3.3.

Parameter Setting

There are two parameters in the MRF-based region-growing criterion, i.e., $\beta$ and $T$. The potential parameter $\beta$ is used to balance the influence between $E_f(R_i,R_t)$ and $E_l(R_i,R_t)$. A high $\beta$ value emphasizes $E_l(R_i,R_t)$ and leads to results with large homogeneous objects. On the contrary, a low $\beta$ value emphasizes $E_f(R_i,R_t)$ and is suitable for getting results with many details. Hence, $\beta$ should select different values for various applications. However, as the relationship between urban and nonurban areas is quite stable, $\beta$ is fixed and is empirically set as 0.05 for simplifying the parameter setting.

The threshold $T$ is used to control the process of region growing. By gradually increasing $T$, small regions labeled nonurban are merged into larger urban regions or nonurban regions, then urban areas are extracted. In practice, we used $T=25$ as the initial threshold and doubled the threshold each time. The final termination threshold was determined by the change of the spectral variance. The assumptions supporting this threshold selection are that urban areas consist of various subobjects and their spectral variance should be large; if the nonurban areas are wrongly recognized as urban areas, an abrupt change of the spectral variance should be observed. Here, we use $\mathrm{CR}(i,i+1)$ to show the change rate of spectral variances, i.e.,

Fig. 2

Example of parameter $T$: (a) urban area with $T=25$; (b) urban area with $T=50$; (c) urban area with $T=100$; (d) urban area with $T=200$; (e) urban area with $T=400$; (f) urban area with $T=800$; (g) urban area with $T=1600$; (h) $\mathrm{Std}\_T(i)$; and (i) $\text{CR}(i,i+1)$.

4.1.

Experiments of Aerial Images

In this experiment, three aerial images, as shown in Fig. 3, are used to test our method and other urban extraction methods. These aerial images were acquired in 2009 and are located in Taizhou City, China. The three images have the same size of $500\times 500$, and the spatial resolution is 0.4 m. The test images contain plane agriculture fields and small villages, where urban objects show various spectral appearances and some nonurban objects are similar to seed points in terms of spectral characteristics. This makes urban detection challenging. Moreover, the following competitive methods are also considered for comparison:

Fig. 3

Experiments of aerial images: (a1)–(c1) aerial images; (a2)–(c2) traditional region growing; (a3)–(c3) Markov random field (MRF); (a4)–(c4) object-based MRF (OMRF); (a5)–(c5) support vector machine (SVM); (a6)–(c6) object-based SVM; and (a7)–(c7) MRF-based region growing.

For the sake of fairness, we chose the same seed points to train the urban areas for the traditional region-growing method and the two SVM methods and deliberately selected samples to train the nonurban areas for these SVM methods as well. We also tuned the parameters of these methods to get their optimal performances. For the traditional region-growing method, we chose the threshold parameter following the instructions in the literature.²⁸ For the two-class SVM, we set the radial basis function as the kernel type, the gamma in kernel function as 0.33, and the penalty function as 100, respectively. For the object-based SVM, we use 0.1% as the ratio of training samples based on the literature.²² Therefore, the comparison can demonstrate the difference between our model and other state-of-the-art methods.

Experimental results of aerial images are shown in Fig. 3. Here, the caption of Fig. 3 consists of two parts, where the first part using the alphabetical order denotes different test images and the second part using the number order denotes different detection methods. Detected urban objects are represented as yellow masks over the test images. From the comparative test, one can see that the proposed method exhibits a remarkable improvement for urban detection. Namely, the traditional region-growing method, as shown in Figs. 3(a2)–3(c2), still has huge misclassifications which belong to different object categories and have similarity spectral appearances. The main reason is that the traditional region-growing method only uses the spectral features which do not consider the spatial constraint. By employing the spatial context information, the classical MRF model has less misclassification of nonurban areas. However, this pixel-level generate model can just recognize the parts of the urban areas with similar appearances, since it cannot model the complex and macropatterns by incorporating the long-range interactions. It also wrongly labels some urban objects as nonurban, such as the roofs of buildings and vegetation. The OMRF model utilizes the regions to describe the macrospatial constraints and improves the classical MRF model, e.g., Figs. 3(a4) and 3(c4). However, the OMRF model usually leaves the characteristic of urban areas out of consideration, which may lead to some undesirable results such as Fig. 3(b4). The SVM method trains data to obtain urban areas. Although it can effectively recognize buildings, urban vegetation objects are sometimes classified as nonurban areas because of the lack of spatial information. The object-based SVM improves the pixel-based SVM and gets results that are more consistent by considering the object semantic information with regional features. Nevertheless, it still cannot sufficiently use spatial information whose results have some misclassifications. Compared with these methods, our MRF-based region-growing method first considers the urban characteristics when we select seed points, then employs the MRF defined on the region level to capture regional spatial constraints, and finally proposes a corresponding region-growing criterion that utilizes these features to detect urban areas. Hence, our method demonstrates a better performance than the other methods.

Experimental results are quantitatively evaluated by the overall accuracy (OA) and kappa coefficient $\kappa$. OA and $\kappa$ are the two indicators that measure the degree of similarity between two images.³¹ If $P_{ij}$ is the proportion of subjects that were assigned to the $i$’th class by the first image and the $j$’th class by the second image and denotes $P_{i•}=\sum _{j=1}^k{P}_{ij}$ and $P_{•j}=\sum _{i=1}^k{P}_{ij}$, then

The OA and $\kappa$ of aerial images are given in Table 1.

Table 1

Comparison of results.

From these quantitative indexes, we know that MRF-based region growing can enhance both the OA and kappa for each experimental image. This also shows that our method extracts a better scope of urban areas than do the other methods. In particular, when the topographic features are complex, the enhancement of indices is obvious. For clarity, the quantitative indices of Table 1 are illustrated in Fig. 4.

Fig. 4

Kappa and overall accuracy (OA) of experiments of aerial images: (a) kappa and (b) OA.

4.2.

Experiments of SPOT5 Images

The effectiveness of the proposed method is further tested in this section. Two SPOT 5 remote sensing images, as shown in Fig. 5, are employed for the next experiment. These test images are located on the Pingshuo area of China. Both sizes are $438\times 438$. These test images mainly consist of three object types, i.e., urban areas, cultivated land, and woodland. Among them, urban green space and woodland and urban building and cultivated land have similar spectral appearances, respectively. This phenomenon increases the difficulty of urban detection.

Fig. 5

Experiments of SPOT5 images: (a and d) Original SPOT5 image; (b and e) ground truth (red); and (c and f) results of MRF-based region growing (yellow).

Experiments of SOPT5 images are illustrated in Fig. 5. Compared with the ground truth, the MRF-based region-growing method performs well and the results are close to the ground truth. This demonstrates that our model can effectively extract urban areas from different datasets.