Persistent Homology for Land Cover Change Detection

Simplicial Complexes

Definition 1: A set of $\text{k}+1$ points, ${\text{x}}_0,{\text{x}}_1,\dots ,{\text{x}}_{\text{k}}$, is affinely independent if the $\text{k}$ vectors, ${\text{x}}_1-{\text{x}}_0,{\text{x}}_2-{\text{x}}_0,\dots ,{\text{x}}_{\text{k}}-{\text{x}}_0$, are linearly independent.

Definition 2: A $\text{k}$-simplex $\sigma$ is the convex hull of $\text{k}+1$ affinely independent points. The value of $\text{k}$ is the dimension of the $\text{k}$-simplex $\sigma$ and the points ${\text{x}}_0$ to ${\text{x}}_{\text{k}}$ are called vertices.

Simplices of dimension $0,1,2,3$ correspond to vertices, edges, triangles, tetrahedral, respectively.

Definition 3: A face of $\sigma$ is a simplex spanned by a subset of the vertices of $\sigma$.

Definition 4: A simplicial complex is a finite collection of simplices, $\text{K}$, such that:

a.

every simplex $\sigma \in \text{K}$, every face of $\sigma \in \text{K}$;

b.

for every two simplices $\sigma$, $\varphi \in \text{K}$, the intersection, $\sigma \cap \varphi$, is either empty or a face of both simplices.

Open in new tab

Figure 3. The entire set of modules involved in the proposed approach.

The dimension of $\text{K}$ is the largest dimension of any simplex that is an element of $\text{K}$. The Euler characteristic of $\text{K}$ is the alternating sum of simplex numbers:

$$\text{χ}(\text{K})={\text{s}}_0-{\text{s}}_1+{\text{s}}_2-\dots \pm {\text{s}}_{\text{k}},$$(1)

where $\text{k}$ denotes the dimension of $\text{K}$ and ${\text{s}}_{\text{i}}$ corresponds to the number of $\text{i}$-simplices, with $0\le \text{i}\le \text{k}$.

Vietoris-Rips Simplicial Complex

The geometrical constructors known as VR simplicial complex are defined as follows.

Definition 6: Letting $\text{R}\ge 0$ be a real number, and $\text{S}$ a finite set of points in ${ℝ}^2$, the VR complex of $\text{S}$ and $\text{R}$, denoted as Vietoris-Rips(R), consists of all abstract simplices in $2\text{S}$ whose vertices are at most a distance $2\text{R}$ from one another.

In other words, any two vertices at a distance at most $2\text{R}$ from each other are connected by an edge, and a triangle (or higher-dimensional simplex) is addended to the complex if all its edges are in the complex (Figure 4).

Open in new tab

Figure 4. The Vietoris-Rips complex associated to six points equidistant on the unit circle.

Persistent Homology and Betti Numbers Computation

The topological characterization of simplicial complexes is performed through the concept of persistent homology whereby Betti numbers are computed. The theory of persistence is inspired by Morse theory and associated with spectral sequences. It brings a quantitative constituent to the Betti numbers that is steady under disturbances. It aims to disclose the different scales under which a cloud of points can be perceived in one and only one coherent formalism.

Let C be a simplicial complex such as (a) isolated point, (b) line segment, (c) triangle, and (d) tetrahedron obtained using VR constructor. Larger simplicial complexes are made up of these basic ones using different values of d (Figure 5). Therefore, several simplicial complexes are exhibited in order to approximate the true shape of the cloud of points (data set). The ultimate goal of persistent homology consists of determining the one that is the most representative of the true shape formed by the cloud of points. It corresponds to the one that produces the qualitative noiseless features.

Open in new tab

Figure 5. 2D simplicial complexes constructed using five points with different values of $\text{d}$. Their building blocks are: (a) points, (b) line segments, and (c) triangles. Betti numbers for each simplicial complex are depicted. This drawing has been generated via MATLAB JAVAPLEX (Homology) ToolBox.

Birth and Death

The concept of birth and death is addressed in this section. The starting point is a sequence of homomorphisms connecting the homology groups of the complexes in the filtration. For $\text{i}\le \text{j}$, ${\text{K}}_{\text{i}}$ represents a subcomplex of ${\text{K}}_{\text{j}}$, which can be viewed as an injective map ${\text{f}}^{\text{i},\text{j}}\text{:}{\text{K}}_{\text{i}}\_\to {\text{K}}_{\text{j}}$. It maps a simplex in ${\text{K}}_{\text{i}}$ to the same simplex in ${\text{K}}_{\text{j}}$. This is why it is called the inclusion map. It carries over to an inclusion map on the $\text{p}$-cycles, ${\text{f}}^{\text{i,j}}{}_{\text{p}}\text{:}{\text{Z}}_{\text{p}}({\text{K}}_{\text{i}})\_\to {\text{Z}}_{\text{p}}({\text{K}}_{\text{j}})$. This in turn creates a map on homology: ${\mathrm{\varphi }}^{\text{i},\text{j}}{}_{\text{p}}:{\text{H}}_{\text{p}}({\text{K}}_{\text{i}})\to {\text{H}}_{\text{p}}({\text{K}}_{\text{j}})$, which does not often represent an inclusion map. More precisely, if $\mathrm{\delta }$ is a class in ${\text{H}}_{\text{p}}({\text{K}}_{\text{i}})$, and $\text{z}\in {\text{Z}}_{\text{p}}\left({\text{K}}_{\text{i}}\right)$ is a representative cycle, let ${\varphi }^{\text{i},\text{j}}{}_{\text{p}}(\mathrm{\delta })$ be the class in ${\text{H}}_{\text{p}}({\text{K}}_{\text{j}})$ that contains ${\text{f}}^{\text{i,j}}{}_{\text{p}}(\text{z})$. It accepts as input a $\text{p}$-cycle in ${\text{K}}_{\text{i}}$ and sends it to ${\text{K}}_{\text{j}}$ For example, if $\text{δ}$ encloses a hole in ${\text{K}}_{\text{i}}$ that fills up at the time ${\text{K}}_{\text{j}}$ is reached, then ${\varphi }^{\text{i},\text{j}}{}_{\text{p}}$ maps $\text{δ}$ to $0\in {\text{H}}_{\text{p}}\left({\text{K}}_{\text{j}}\right)$. The image of ${\varphi }^{\text{i},\text{j}}$ p is called a persistent homology group. It includes all $\text{p}$-dimensional homology classes that already have representatives in ${\text{K}}_{\text{i}}$. A class $\text{δ}\in {\text{H}}_{\text{p}}\left({\text{K}}_{\text{i}}\right)$ is born at ${\text{K}}_{\text{i}}$ if $\text{δ}$ is not present in the image of ${\varphi }^{\text{i}-\text{1},\text{i}}{}_{\text{p}}$, and a class born at ${\text{K}}_{\text{i}}$ dies entering ${\text{K}}_{\text{j+1}}$ if ${\varphi }^{\text{i},\text{j}}{}_{\text{p}}(\mathrm{\delta })$ is absent in the image of ${\varphi }^{\text{i}-\text{1},\text{j}}{}_{\text{p}}$ but ${\varphi }^{\text{i},\text{j}+\text{1}}{}_{\text{p}}(\text{δ})$ is present in the image of ${\varphi }^{\text{i}-\text{1},\text{j}+\text{1}}{}_{\text{p}}$.

Betti Numbers

A simplicial complex exhibits a group structure through the addition of $\text{k}$-simplices. The resulting free group is called the chain group, ${\text{C}}_{\text{k}}$. The notion of Betti numbers emanates from the quotient group ${\text{H}}_{\text{k}}={\text{Z}}_{\text{k}}/{\text{B}}_{\text{k}}$, in which ${\text{Z}}_{\text{k}}$ denotes the group structure of all $\text{k}$-chains that have empty boundary; and ${\text{B}}_{\text{k}}$, subgroup of ${\text{C}}_{\text{k}}$, is also a subgroup of ${\text{Z}}_{\text{k}}$, called the boundary group. The number of distinct equivalent classes of ${\text{H}}_{\text{k}}$ corresponds to the $\text{k}$-th Betti number denoted by ${\beta }_{\text{k}}$. Indeed, the $\text{k}$-th Betti number tallies the number of $\text{k}$-dimensional holes in the simplicial complex. In the case where $\text{k}=0$, ${\beta }_0$ counts the number of path-connected components of the simplicial complex. In the case of subsets of the Euclidean space ${ℝ}^3$, ${\beta }_{\text{1}}$ counts the number of independent tunnels (holes) and ${\beta }_{\text{2}}$ counts the number of cavities, or enclosed voids. For example, in the case of the solid torus, the Betti numbers are ${\beta }_0=\text{1}$, ${\beta }_{\text{1}}=\text{1}$, and ${\beta }_{\text{2}}=0$. Therefore, one can state that a torus is topologically represented by the series $1,1,0$. However, the surface of a torus has ${\beta }_0=\text{1}$ (one connected component), ${\beta }_{\text{1}}=\text{2}$ (the one in the center and the one in the middle), and ${\beta }_{\text{2}}=\text{1}$ (the inside of the donut). The idea of describing local structures via Betti numbers is much more efficient than the one using Euler characteristic (a connectivity measure that is a zero-dimensional Minkowski functional).

Barcodes

Betti numbers are computed for each simplicial complex (by varying $\text{d}$). From each stage to the next, pairing up the births and deaths, a set of intervals or bars called the barcode of the filtration is obtained. Each bar represents a class in one of the homology groups and thus has a definite dimension. Figure 6 depicts an example of a cloud of five points representing a “house with 1 hole.” The persistent homology (barcodes) discloses for $2\le \text{d}\le 2.8$, one connected component $({\text{β}}_0=1)$, one hole $({\text{β}}_1=1)$, and no “voids” or “cavities” $({\text{β}}_2=0)$, yielding a series of Betti numbers $1,1,0$.

Open in new tab

Figure 6. Chain of simplicial complexes with their barcodes. The persistent homology discloses one connected component $({\text{β}}_{0=1})$, 1 hole $({\text{β}}_{1=1})$.

Traditional HMM

The HMM approach views an observation sequence $\mathcal{O}=o_1{o}_2{o}_3\dots o_T$ as a series of symbols produced by a sequence of hidden states with a certain emission probability distribution. Each symbol is emitted from one single hidden state using this probability distribution, and the hidden state distribution is a Markov chain. The following definition exhibits the statistical parameters of a tradition HMM.

$\text{B}=[{\text{b}}_{\text{j}}(\text{t})]$ is the state conditional probability matrix of the emitted symbol ${\text{o}}_{\text{t}}$ given a hidden state ${\text{q}}_{\text{j}}$, with $({\sum }_{\text{t}}{\text{b}}_{\text{j}}(\text{t})=1)\forall \text{j}$, $\text{T}$ is the number of observations ${\text{o}}_{\text{i}}$, each one of length ${\text{r}}_{\text{i}}$.

As an example, Figure 8 depicts a simple HMM graph.

Open in new tab

Figure 8. A simple HMM Graph with two hidden states and three symbols.

Three problems are assigned to an HMM: (a) probability evaluation, (b) statistical decoding, and (c) parameter estimation (or training). In summary, the evaluation problem in the HMM consists of evaluating the probability for the model $\text{λ}=[\text{π},\text{A},\text{B}]$ that most likely produces the sequence $\mathcal{O}$. This probability can be expressed as

Using state conditional independence of the observation sequence $\mathcal{O}$, one obtains

The evaluation problem is based on the observation state conditional independence; however, there are numerous applications in which the conditional independence condition does not hold.

Topological HMM: Mathematical Framework

Because the proposed approach extracts qualitative features (using homology theory) from local shapes described by simplicial complexes, it is called a topological hidden Markov model (THMM). To provide a complete introduction of the THMM formalism, this section discusses the THMM context, background on the VR simplicial complexes constructor, and search and characterization of local structures.

THMM Context

In order to capture shapes, the THMM approach views an observation sequence $\mathcal{O}=O_1{O}_2\dots O_T$ of length $\text{T}$ as a series of constituents $O_i$ made of strings of symbols $o_{ij}\in \sum$ interrelated in some way. Each observation sequence $\mathcal{O}$ is not only one sequence in which all constituents are conditionally independent but also a sequence that is divided into a series of $\text{T}$ strings $O_i$$(O_i=o_{i_1}{o}_{i_2}\dots o_{i{r}_i})$, $(1\le \text{i}\le \text{T})$, ($r_i$ is the length of $O_i$). Each $O_i$ is produced from a set of local structures ${\text{S}}_{\text{i}}$. The observed symbols $o_{ij}$ are emitted from hidden states ${\text{q}}_{\text{ij}}$.

Parameter Estimation

The parameters $\text{π}$, $\text{A}$, and $\text{B}$ are computed in the same way as in traditional HMM using a Baum-Welch parameter estimation (or expectation maximization), however the matrices $\text{C}=[{\text{c}}_{\text{ij}}]$ and $\text{D}=[{\text{d}}_{\text{ij}}]$ that unfold the local structures are estimated as follows:

where $\left|{\text{cluster S}}_{\text{i}}\right|$ designates the cardinality of cluster ${\text{S}}_{\text{i}}$ and $<{\text{S}}_{\text{i}},+>$ represents the number of times the cluster ${\text{S}}_{\text{i}}$ occurs in the training set.

Remote Sensing Feature Generation

Before any classification task such as change detection can take place, one has to extract discriminative and less costly features. It is well known that the classifier power generalization depends in part on the dimension of the feature space. Feature generation (aka extraction) in change detection represents a challenging and tedious task due to diverse regional characteristics such as (a) area scattering: urban regions usually account for small areas and are scattered in different locations; (b) large regions: vegetation and bare soil region represent large and homogeneous surfaces; and (c) thematic similarity: spectral features between some of the thematic classes such as “less dense urban” and “bare soil” are very much alike especially when the spatial resolution is limited. In traditional approaches, feature extraction in change detection is usually performed using software such as ENVI and eCongnition. These programs invoke spatial features (e.g., texture, area, orientation, size) as well as spectral attributes (e.g., histogram depicting occurrence frequencies of gray values, pixel radiometric values in different bands). Because this work is driven by the extraction of local structures that account for image spatial changes between two points in time, it is deemed sufficient to rely on radiometric values at different bandwidths to represent a pixel in an image.

Building VR Simplicial Complexes and Local Structures from a Pixel

The amount of variation of pixels’ radiometric values is an indicator of the type of change in the image. For example, in the case of a disaster, the image exhibits a drastic radiometric change, whereas in the case of urban growth, this radiometric variation is rather smooth. Each image pixel contributes to the formation of a simplicial complex. An observation sequence that represents a series of radiometric vectors in three bands is associated with a pixel. This sequence represents a cloud of points from which a simplicial complex is built. A two-dimensional discrete wavelet transform (DWT) is applied to the entire image of dimension $(\text{r,c})$ for each band out of three. Therefore, this decomposition produces an average image and a sequence of detailed images in each scale. The inverse function of the wavelet transform is applied to these images in order to restore the original dimension $(\text{r,c})$. The goal of this decomposition is to associate each pixel of the original image to its corresponding pixel in a one-to-one mapping.

Let $\text{D}\text{=}\left[{\text{A}}_{\text{j}},\left\{{\text{D}}_{\text{j}}^{\text{k}}\right\},\text{j}=1,\dots ,\text{J and k}=1,\dots ,\text{K}\right]$ be the decomposition in which ${\text{A}}_{\text{j}}$ refers to the approximation image at the $\text{j}$-th scale and ${\text{D}}_{\text{j}}^{\text{k}}$ is the detail image at scale $\text{j}$ and orientation $\text{k}$. For the work described herein, a wavelet decomposition with four levels $\text{(scale}\text{J}=4)$ and three orientations $\text{(K}=3)$ for each level, was undertaken. For each scale, the $\text{i}$-th pixel in the original image is associated with a vector containing its three values in the detail images. This procedure is repeated for each of the four scales and then for each band. Hence, DWT decomposition produces a set of 12 points in the details vector space. Figure 10 shows an example of a simplicial complex (built from the 12 points generated from a single pixel using DWT) in the 3D frequency bands. Betti numbers’ qualitative features are extracted from this simplicial complex and then assigned to its corresponding local structure cluster (a class of Betti numbers). This information is fed into the THMM framework for the purpose of classification.

Open in new tab

Figure 10. An image pixel is mapped to a cloud of 12 3-D points of the transformed frequency band vector space to build a pixel simplicial complex.

THMM in Change Detection

Change detection is a classification problem with only two classes—Changed and Unchanged—that applies to a data sequence extracted from the pixel of the remote satellite image at two different points in time, ${\text{T}}_1$ and ${\text{T}}_2$. The THMM formalism fits seamlessly with the change detection application at hand. Because a pixel is an input, a VR simplicial complex is extracted from this pixel.

Open in new tab

Figure 11. THMM associated with the same pixel of an image at two different points in time, ${\text{T}}_1$ and ${\text{T}}_2$.

he entire sequence is $\mathcal{O}={\text{O}}_1{\text{O}}_2,$ whereby each ${\text{O}}_{\text{i}}={\text{o}}_{\text{i}1}{\text{o}}_{\text{i}2}\dots {\text{o}}_{\text{i}12}$ designates the 12 3-D points associated with a pixel of the image at time ${\text{T}}_{\text{i}}(\text{i}=1,2)$. A local latent structure $\text{Si}$ mapped to each ${\text{O}}_{\text{i}}(\text{i}=1,2)$ corresponds to the Betti number cluster computed from the simplicial complex generated by the VR constructor. In conclusion, the following instantiations are made:

•

The classes ${\text{ω}}_1$ and ${\text{ω}}_2$ correspond to Changed and Unchanged categories, respectively.

The sequence ${\text{O}}_{\text{i}}$ corresponds to the 12 3-D points generated for a pixel using DWT.

The symbol $\text{Si}$ represents the Betti number class (a cluster, aka local structure) assigned to the sequence ${\text{O}}_{\text{i}}.$A symbol ${\text{o}}_{\text{ij}}$ of a sequence ${\text{O}}_{\text{i}}$ represents a 3D point among these 12 points.

•

A hidden state q_i is a frequency band $\text{Bi}$ selected during the DWT procedure. A set of frequency bands has been selected from which the 12 points have been generated.

Finally, maximum likelihood is applied in the computation of the best ${\text{ω}}^{*}$ among ${\text{ω}}_1$ and ${\text{ω}}_2$ according to the following criterion:

$$\underset{\text{ω}}{\max }\text{Prob}\left(\mathcal{O}\text{|ω}\right)=\underset{\text{ω}}{\max }\left[\sum _{\text{S}}\text{Prob}\left(\mathcal{O},\text{S|ω}\right)\right].$$(10)

It is worth noting that maximum likelihood-based technique is the one most commonly used for land cover classification especially for medium resolution imagery. However, classification based on maximum likelihood alone has many limitations in separations of remote sensing data such as salt-and-pepper effects and challenges caused by mixed pixels, especially for high-resolution images. To overcome this drawback, maximum likelihood is used only in the computation of the best class ${\text{ω}}^{*}$ among ${\text{ω}}_1$ and ${\text{ω}}_2$ that optimizes the criterion defined in Equation 10. It is worth noting that during the supervised training, the sequences $\mathcal{O}={\text{O}}_1{\text{O}}_2$ are tagged changed or unchanged, depending on whether a pixel has changed or not in the image viewed at time ${\text{T}}_2$. The ratio of the number of changed pixels to the total number of pixels in the image represents the magnitude of change. Figure 11 illustrates the instantiation process of the THMM for change detection.

Data Collection

To analyze the efficiency of the proposed model, two change detection data sets were selected, namely the Change Detection dataset and Sztaki AirChange Benchmark. The Change Detection data set comprises 1,000 pairs of images with $(800\times 600)=480,000$ pixels each. An image pair corresponds to one reference image and one test (or moving) image. Each pair is associated with a binary ground truth image with the same resolution. The images were constructed by the serious game Virtual Battle Station 2 with Bohemia Interactive Simulations. The data set consists of 100 different scenes containing several thematic (vegetation, water, buildings) and moderate ground reliefs. Each scene was viewed from five different angles, where camera positions were varied at steps of 10 degrees according circle of radius 100 meters at approximately 250 meters high, and with a fixed tilt of about -70 degrees. All images were acquired with a resolution of about 50 cm per pixel. To obtain more information about this data set, please refer to the “Change Detection dataset”. Figure 12 depicts a sample of image pairs at two points in time ${\text{T}}_1$ and ${\text{T}}_2$ and their ground binary truth images, highlighting the changed area with white dots.

Open in new tab

Figure 12. Samples of images at time ${\text{T}}_1$ (targets: column 1), same images at time ${\text{T}}_2$ (moving: column 2), and binary ground truth corresponding images (column 3). White dots correspond to the changed areas.

The Sztaki AirChange Benchmark is the second database used. It comprises 13 aerial image pairs of size $952\times 640$ and with a resolution of 1.5 meters and their corresponding ground of truth. Images constituting each pair were acquired at different times over several years, allowing the generation of relevant changes, such as new built-up regions, building operations, planting of large groups of trees, and so on. This database also contains challenging problems in which reversible changes are considered as unchanged, where variations of the vegetation greenery or even shadows are labeled as unchanged in the ground truth. Figure 13 illustrates some samples of image pairs and their corresponding masks (ground truth images).

Open in new tab

Figure 13. Samples of images at time ${\text{T}}_1$ (first column), same images at time ${\text{T}}_2$ (second column), and their corresponding ground truth images (third column).

Training Phase

To measure the generalization power of the THMM classifier, a three-fold cross-validation estimation technique is performed in the entire set of 1,000 pairs of pixels. This set is divided into three subsets of 333 pairs of images each, and selected two subsets (666 pairs) for training and one subset for testing (333 pairs). All feature vectors are extracted using the MATLAB platform (MATLAB, 2015). The visible observation sequences $(\mathcal{O}=O_1{O}_2)$ are extracted from these 666 pairs of images. This procedure is performed using pixels from both classes, Changed and Unchanged, separately. Five Gaussian mixtures $(\text{M}=5)$ are used for the estimation of the emission probability values $\text{Prob}(\text{o}{}_{\text{t}}|{\text{q}}_{\text{t}})$. Therefore, two optimal models ${}_1^{*}=[{\pi }_1^{*},A_1^{*},B_1^{*},C_C^{*},D_1^{*}]$ and ${}_2^{*}=[{\pi }_2^{*},A_2^{*},B_2^{*},C_2^{*},D_2^{*}]$ are obtained during each training phase. Testing is conducted following each training phase. This procedure is repeated three times with three different subsets for testing. The predicted classification for each pixel is compared to the corresponding pixel in the ground binary truth image. The results in each training/testing procedure are finally averaged to obtain an overall classification performance.

Testing: Decision and Classification

It is well known that a classification problem taking into account a reject option (whereby the classifier has the possibility to hold back its affectation) incurs a cost lower than misclassification (Duda, Hart, & Stork, 2001). Indeed, there are a significant number of cases in which radiometric values of a group of pixels have changed slightly because of a modification in the weather conditions (e.g., presence of clouds). These image regions are often classified as Changed, but no change has actually occurred. The geometric characterization of local structures assigned to these regions can only partially mitigate this false positive effect. A reject threshold can be determined via the risk minimization given by loss matrix. However, the loss matrix element $\text{λ}$ representing the loss incurred for choosing the reject option is computed using a selected validation set. Furthermore, when the THMM classifier makes a decision for a test observation sequence, it should take into account both the potential gain and potential loss. An accepted low-risk pixel increases profit, while a rejected high-risk pixel decreases loss. The situation is much more critical and far from symmetric in the area of remote sensing change detection. Denote ${\text{α}}_{\text{i}}(\text{i}=\text{1},\text{2})$ the action of assigning a single pixel (viewed via $\mathcal{O}$)to class ${\text{ω}}_{\text{i}}$ (${\text{ω}}_1$ = Changed and ${\text{ω}}_2$ = Unchanged, respectively) and ${\text{α}}_{\text{3}}$ the action assigned to the “reject” (or doubt) decision. Likewise, designate ${\lambda }_{ik}$ as the loss incurred for taking action ${\alpha }_i$ when the input actually belongs to class ${\omega }_k$. For an optimal classifier performance, statistical decision theory advocates the action ${\alpha }_i$ with a minimum conditional risk:

A correct classification is assigned a zero loss; however, the classification errors are not equally costly in the remote sensing change detection problem. In fact, wrongly classifying the changed pixels is more critical than misclassifying unchanged pixels. To discriminate wrong decision losses, the parameters ${\text{λ}}_{12}$ and ${\text{λ}}_{21}$ are set to 1 and 0.5, respectively. The loss matrix is then expressed as

Therefore, the conditional risk of reject is

To precisely evaluate the system performance, the receiver operating characteristic (ROC) curve is invoked (Li & Guo, 2014). The results from the folds are then averaged for each performance measure using k-fold cross-validation estimation technique to produce a single estimation.

Experimental Results and Interpretation

The THMM approach has been compared with the traditional HMM. In the HMM framework, a hidden state is assigned a discrete value representing one of the thematic classes (such as “water,” “grass,” etc.). The image was first segmented into thematic classes and each pixel belongs to one thematic class. The THMM has been implemented using the VR approach using the cross-validation criterion that computes the optimal $\text{R}*$ with the maximum accuracy (Figure 14). It is worth noting that the classification performance depends on the quality of features and the classifier’s parameters. Different connectivity parameters have been tested during the training phase (parameter tuning phase). The $\text{R}$-value $(\text{d}=2\text{R})$ that has produced the highest accuracy (95%) has been selected $({\text{R}}^{*}=0.9)$.

Open in new tab

Figure 14. Cross- validation graph showing the optimal ${\text{R}}^{\text{*}}$ value equal to 0.9 that produced the highest accuracy equal to 95% as well as the persistent interval [Birth..Death] = [0.6..0.9] in red.

This selection of a unique $\text{d}$ value represents one way to proceed because a filtration (for continuous $\text{d}$-values) of a finite simplicial complex in which persistent homology exhibits qualitative features (barcodes) can also be performed. The choice of selecting a unique $\text{d}$ value is justified by the risk of extracting data inherent noise or artifacts. Furthermore, an experimental classification reject threshold value is set to $\text{λ}=0.02$. Figure 15 illustrates an example of change detection result obtained on a test image applied to a change detection data set. Figure 16 depicts a sample of qualitative results over three aerial data sets, namely the Szada, Archive, and Tiszadob. The obtained binary images correspond perfectly to the change mask of the ground truth. Furthermore, these results show that the proposed approach produces homogeneous and smoothing effect. One of the main strengths of the proposed method based on VR simplicial complex representation is its ability to distinguish between real changes and false changes caused by shadows, or seasonal weather effects.

Open in new tab

Figure 15. Samples of images at time $\text{T}1$ (targets: column 1), same images at time $\text{T}2$ (moving: column 2), and binary ground truth corresponding images (column 3). White dots correspond to the changed areas.

Open in new tab

Figure 16. Qualitative results over three data sets. Column (a), image at ${\text{T}}_1$; column (b), image at ${\text{T}}_2$; column (c), ground truth; column (d), the proposed model results.

For comparison purposes, three classifiers have been implemented: (a) an HMM based only on radiometric values (b) and an HMM that takes into account textural information assigned to a pixel. In order to conduct a comparison between these three classifiers, different values of the classification reject threshold $\text{λ}$ that allows different loss matrices to be exhibited are considered. The ROC curves of these classifiers are computed and depicted in Figure 17.

As a second experiment task, the THMM has been compared with other competitive techniques in the field of machine learning (see Dai & Khorram, 1999; Volpi et al., 2013). These techniques were brought inside a MATLAB platform. To ensure reproducibility of the results, this comparative work is conducted on the same database, Change Detection dataset. In other words, some results from the state-of-the-art benchmarked machine learning techniques are provided to derive a global assessment and to widely position the contribution of the THMM formalism. To promote reproducibility and perform a fair model comparison between different change detection approaches, the same set of features prior to classification are extracted.

Open in new tab

Figure 17. ROC curves for different values of the reject threshold $\text{λ}$ of the HMM $(\text{AUC}1=0.8548)$, HMM (considering textural features) $(\text{AUC}2=0.8748)$, and THMM $(\text{AUC}3=0.9083)$.

Tables 1 and 2 indicate that the THMM outperforms some state-of-the-arts machine learning methodologies when applied to LCCD on both databases (Change Detection dataset and Sztaki AirChange Benchmark). For the purpose of comparison, experiments were conducted on both publicly available data sets, allowing for a fair evaluation. In Figure 17, the area under curve assigned to the THMM is larger than the one assigned to HMMs. The results shown in Tables 1 and 2 clearly underscore the advantage of the THMM over traditional change detection classifiers. The THMM outperforms the other methods because the proposed formalism exploits geometric information, which is highly sensitive to real changes (expressed by geometric representation of pixel sequence). However, these geometric clues cannot always distinguish between real and false changes. According to the geometric formalism, many reversible changes are also labeled as changes. The topological aspect was considered by introducing Betti numbers in extracting features to separate false changes from detected pixels.