INTERACTIVE DIGITAL TERRAIN MODEL ANALYSIS IN ATTRIBUTE SPACE

: The use of high-resolution digital terrain model derived from airborne LiDAR system becomes more and more prevalent. Effective multi-scale structure characterization is of crucial importance for various domains such as geosciences, archaeology and Earth observation. This paper deals with structure detection in large datasets with little or no prior knowledge. In a recent work, we have demonstrated the relevance of hierarchical representations to enhance the description of digital elevation models (Guiotte et al., 2019). In this paper, we proceed further and use the pattern spectrum, a multi-scale tool originating from mathematical morphology, further enhanced by hierarchical representations. The pattern spectra allow to globally and efﬁciently compute the distribution of size and shapes of the objects contained in a digital elevation model. The tree-based pattern spectra used in this paper allowed us to analyse and extract features of interest. We report experiments in a natural environment with two use cases, related to gold panning and dikes respectively. The process is fast enough to allow interactive analysis.


INTRODUCTION
Data analysis is usually achieved through a representation of the data in an appropriate feature space.In the context of digital images, the histogram of graylevels has been used for decades as a simple and efficient probability density function to characterize image contents.While today feature learning achieved with deep neural networks has shown great success including for remote sensing data, it still requires the availability of massive amounts of data, either with or without annotations (supervised and unsupervised learning, respectively).In the context of interactive data analysis, efficient methods that do not impose such requirements remain appealing.
Among existing methods for representing an image content in a predefined feature space, the pattern spectra (Maragos, 1989) have established as an advanced solution to describe the image through a probability distribution function of some attributes measured on the image parts.More precisely, and conversely to the histogram of graylevels, the analysis is not conducted at the pixel level but rather on all components or regions present in the image.The underlying scale-space is efficiently built using either level sets or multiscale segmentation, through inclusion or partitioning trees respectively (Bosilj et al., 2018).As far as the attributes are concerned, pattern spectra give access to a wide range of properties beyond the distribution of graylevels, e.g.related to the size or shape of images components.Let us note that the attributes can be combined to provide a multidimensional attribute space in which a further analysis can be conducted.For instance, area, non-compactness and Shannon entropy have been successfully used for aerial image retrieval (Bosilj et al., 2016).Interactive filtering of satellite images has been explored in (Ouzounis, Soille, 2011, Gueguen, Ouzounis, 2012).
However, to the best of our knowledge, pattern spectra have * Corresponding author never been used in the context of digital terrain model (DTM) data analysis yet.

FROM DTM TO ATTRIBUTE SPACE
In this paper, we will use an attribute space to characterize structures present in a DTM.
Direct visual interpretation of DTM is feasible for medium and large structures, but it can be difficult especially in mountainous areas where a high range of vertical information is present (cf. Figure 1a).While hillshades appear as an alternative way for visual interpretation that is especially efficient for microreliefs, it remains difficult to perceive the scale and the depth of objects through hillshades (Figure 1b).We claim that pattern spectra can be used as a relevant alternative, since they provide an automatic tool to deal with the high range of vertical values while preserving micro-reliefs.

Morphological Hierarchy
An interesting mean to represent the hierarchy in an image is to use morphological operators.Among them, min-trees, maxtrees and pattern spectra (Salembier et al., 1998) provide reliable models and descriptors, and will be used in this paper.Min-and max-trees form a hierarchical decomposition of an image X from a domain E ⊂ R 2 with E → Z or R based on level sets of flat zones and are briefly described below.
Max-tree is composed of nodes, a set of flat zones N k h linked together in relation to their level h.The flat zones are the peak components P k h (X) valued with the level h.Peak components are determined according to thresholds T h (X) over the points x of the image X such that:  k is an index over the set T h .The root node of the tree includes the whole image X with the lowest level h of X, while the leaves contain local maxima.All flat zones are nested for decreasing values of h and this results in the so-called max-tree.Min-tree is constructed conversely using lower level sets.Because of duality, the min-tree hierarchy can also be constructed as a max-tree of the inverted image (−X).Tree of shapes introduced by (Monasse, Guichard, 2000), describe the image in a self-dual way, similar to the merging of min-and max-tree.The tree characterizes both local maxima and local minima.
A full image can be reconstructed directly from its min-, maxtrees and tree of shapes.Efficient techniques exist to construct such trees at a low computational cost.In this paper, the elevations of the DTM are viewed as graylevels to construct min-, max-trees and tree of shapes (note that other trees can be used, including fine-grained partition trees such as the α-tree).
These morphological hierarchies capture the structures inside the DTM.To do so, we rely in this paper on the characterization functions known as pattern spectra.

Pattern Spectra
For each node N k h of a tree, many criteria (a.k.a.attributes) related to the properties of the peak components (e.g.shape or size, see next section for a discussion of available attributes) can be computed.A 1D pattern spectrum can be viewed as the probability density function related to the probability that a component with a given attribute is present in the image.In a similar way, a 2D pattern spectrum can be viewed, for a whole image, as the joint probability density function of peak components with two attributes.A direct link between a pair of attribute values and the corresponding areas in the image can be established and depending of the chosen attributes, the different bins of pattern spectra can highlight very meaningful areas.
The pattern spectra can be seen as a histogram describing the distribution of properties (e.g.sizes or shapes) present in the tree.The size and shape attributes are first split into ranges of size or shape classes.We used two rules to create such classes: Linear partitioning that consists in a straightforward range partitioning between the minimum and the maximum of the attribute with a fixed step.We use this for attributes with a normal distribution (e.g.compactness).
Geometric partitioning where a logarithmic range partitioning appears suitable to describe wide distributions while keeping the ability to characterize small values (e.g.area, height, volume).
In this paper, 2D pattern spectra will be used with specific attributes to highlight key information in DTM.
The construction of the 2D pattern spectra S is as follow: 1. Choose the size i × j of the spectrum.
2. Choose two attributes (e.g.area A and compactness C).
3. Define the classes of the two attributes according to the size i and j and the previous partitioning rules.
Then for each node N k h of the tree: 1. Compute the classes ci and cj of the two attributes.

A DTM PERSPECTIVE ON ATTRIBUTES
From the chosen tree, we recursively compute attributes of the nodes.We consider here different attributes: the surface area of the node (i.e. the pixel count in the node).perimeter: P (N k h ) the perimeter of the node (i.e. the pixel count of the node contour).compactness: C(N k h ) degree of compactness of the shape of the node, defined as Compactness ranges from 1 for compact shapes (e.g.circles) to 0 for non-compact shapes.height: H(N k h ) difference between the elevation of the parent of the node N kp hp and the elevation of the deepest node in the subtree rooted in the node.volume: V (N k h ) of the node is the area times the elevation difference with its parent (so A(N k h ) × δ h ), plus the sum of the volumes of all children of the node.mean altitude: M (N k h ) of the pixels contained in the node.altitude dynamics: D(N k h ) of the node, the difference between the altitude of the deepest minima of his children and the altitude of his closest ancestor that has one of his children with deeper minima.
One major advantage of using hierarchical representations on DTM is the direct link between the altitude of the node in the tree and the altitude of the objects in the DTM.As a result, in addition to the conventional area attribute that can be formulated in square meters, the height and volume attributes can be translated with physical meaning (i.e.respectively in meters and cubic meters).
Depending on the application, we choose one or several attributes to characterize the structures in the DTM.

EXPERIMENTS
We carried out experiments on a large dataset to assess the efficiency of our method.

Dataset
The study area covers 255 km 2 of tropical forest.The data was acquired by an aerial multi-echo LiDAR system at 100 m above ground level.The point cloud is first processed into a DTM of 1 m 2 pixel resolution.To remove trees and vegetation, only ground points are kept.The raster is processed through triangulated interpolation.The resulting dataset is a large DTM with 32 bits floating-accuracy elevations.For the sake of readability, we used for this paper a sample of this dataset with a size of 5000 × 5000 pixels (Fig. 1a).The ground truth of the areas of interest has been provided by geologists and geophysicists who are interested in automatically detecting natural lineaments or man-made structures.In the first case, it will be a question of locating structures rich in raw materials, in the second, the zones of gold panning which are a threat to the environment.
Two classes are presented as examples: gold panning zones, shown in Fig. 3a, and dikes.

Characterization of gold panning sites
In order to illustrate the potential of our method, we deal here with the characterization of gold panning areas.
The first step is to compute the hierarchical representations of the DTM.We used the Higra library (Perret et al., 2019) to process several hierarchical representations and corresponding attributes.We chose several component trees: max-tree, min-tree and tree of shapes.For each node of a tree, we compute the attributes listed in Sec. 3. We then obtain a distribution of attribute values to be described.The range of the attributes were divided into classes.
We processed a spectrum by choosing two attributes following the procedure detailed in Sec.2.2. Figure 2 shows the spectrum of area and compactness with the max-tree.
The next step is to find the nodes corresponding to the area of interest.We used the ground truth (Fig. 3a) as a pixel activation map.To select the nodes of the tree corresponding to the gold panning areas, we traversed the tree in depth, from leaves to root.A node is selected if all the pixels of the activation map included in the node are active.If the node contains a disabled pixel, the node and its ancestors are not selected.These selected nodes were used to process a new spectrum of the areas of interest.
To assess the relevance of a pair of attributes, we then searched for separability between areas of interest and background in the spectrum.We defined a metric based on the intersection of the spectrum of selected nodes and the spectrum of background nodes.The intersection was normalized with the weights of the selected nodes such as: • intersection = 1 if the selected nodes are fully merged in the background nodes spectrum.• intersection = 0 if the selected nodes do not share a single bin with the background spectrum.
This intersection metric is very severe for spectra since an object well segmented is likely to have numerous child nodes tree count Max-tree 10 Min-tree 23 Tree shapes 6 Table 1.Count of significant intersections (below 0.95) for each tree type.
mixed with small noise-like objects inherent to hierarchical representations.However, this intersection metric is good enough for ranking the spectra among themselves.
We ran series of spectrum intersections by combining the attributes as well as the trees (i.e.max-tree, min-tree and tree of shapes).We compared the number of meaningful spectra per tree by counting the number of spectra with intersection below 0.95 (Table 1).The min-tree achieves good results in this context (Fig. 5).The gold panning sites have pits in the ground well characterized by the min-tree (i.e.structures that are lower than their neighborhood).Usually max-tree performs well on DTM and appeared more locally consistent in the spectra (Fig. 4).
Tree of shapes seemed ill-suited for this experiment (Fig. 6).However, tree of shapes performances can be explained by our lack of rule to create dual classes with a normal distribution (e.g.height attributes ranges from −40 to +40 in the tree of shapes).Their dual representation calls for further investigations.
To get the most meaningful spectra, we selected the attribute combinations with the lowest intersection scores (Table 2).The best combination was achieved by the mean altitude and altitude dynamics spectrum (Figure 4a).The best attribute combination had in common the mean altitude and the altitude dynamics.Both are strongly correlated with nodes altitude.In this sample zone, the gold panning areas shared the same elevations.To limit biases, we selected a subset of attributes unrelated to absolute altitude.These attributes are mostly related to shape (e.g.area, perimeter, compactness, height and volume).The best intersection scores of these attributes are visible in Table 4.Among them, we found that the best combination was the area and height spectrum (Figure 4b An interesting spectrum combines the compactness and volume attributes (Figure 4c).This last spectrum is an underlying combination of area and height (via volume) together with compactness attributes.
We evaluated 3D pattern spectra following the same methodology.To build the 3D spectrum, we chose 3 attributes and create the subsequent classes.We ran the experiments with the previously selected attributes.The intersection metrics are available in Table 4.We can notice an overall better separability.The best attribute combination was area, height and compactness (Figure 8).In this figure, one can very well perceive a 3D cluster that characterizes well the gold panning areas.
Let us note that the 2D and 3D pattern spectra presented here can be easily extended to higher dimensions, but they will then require more advanced visualisation techniques.

Dikes Extraction
In this section, we discuss the feasibility of structure extraction from DTM without any prior knowledge.Our goal was to extract structures of interest such as dikes.
We used the max-tree according to the observations made in the previous section.In addition, our aim was to characterize above-ground structures.
We plot the spectrum in Fig. 7 (left).We can then interactively select bins in the spectrum.For this example, we made a rectangular selection in green at the bottom center of the spectrum to select structures with height between 35 m and 60 m and compactness around 0.10 (i.e.rather non-compact shapes).We pruned the tree to keep only the nodes from the selected classes.• Retrieve the footprint of the structures.It is useful for visual inspection by displaying the results as an overlay over the DTM or hillshade visualizations (on right side of Figure 7).
• Reconstruct the altitudes of the pruned tree with direct filtering (Salembier et al., 1998).
• Reconstruct the selected structures from the tree with the subtractive filtering.Subtractive filtering introduced by (Urbach, Wilkinson, 2002) allows to sum the altitude differences of the nodes with their direct parents.The filtering result, when used with DTM, is the height of structures selected in the tree.
For dike detection, a false positive structure detection is visible at the bottom left of the image (Fig. 7 right).Since compactness is a ratio involving the perimeter, it is expected but undesired to get such results with non-convex shapes.In future works, we will investigate the use of elongation attributes such as the first moment invariant of Hu (Hu, 1962) for this use case.

CONCLUSION
This study focuses on DTM structure characterization with the pattern spectra.The pattern spectra offer many benefits when used with digital elevation model (DEM).On the one hand, the pattern spectra offer a global and multi-scale description of the objects contained in the DEM.On the other hand, the direct link between DEM values and level definition used in the pattern spectra gives an intrinsic meaning of the attributes commonly used on mathematical morphology.Furthermore, the use of underlying hierarchical representations to compute the pattern spectra and reconstruct characterized objects enables interactive applications.
Future works will be carried out on the use of new attributes suitable with DEM and user interface for 3D spectrum exploration.
Figure 1.Visualizations of a DTM below vegetation in a mountainous area generated from multi-echo LiDAR system.

Figure 2 .
Figure 2. Pattern spectrum (area and compactness) of the DTM with the max-tree.
Figure 3. Visualizations of the ground truth used for the experiments.The ground truth is displayed over the hillshades of Fig. 1b.Yellow surfaces are gold panning areas and the purple square is the closeup (b).

Figure 4 .
Figure 4. Activation maps of pattern spectra over the gold panning areas with the max-tree.The bin color represents the percentage of nodes defined as gold panning in the global spectrum.Values range from 0% (purple) to 100% (yellow).The more consistent the spectrum, the better the characterization.

Figure 5 .
Figure 5. Activation maps of pattern spectra over the gold panning areas with the min-tree.The bin color represents the percentage of nodes defined as gold panning in the global spectrum.Values range from 0% (purple) to 100% (yellow).

Figure 6 .
Figure 6.Activation maps of pattern spectra over the gold panning areas with the tree of shapes.The bin color represents the percentage of nodes defined as gold panning in the global spectrum.Values range from 0% (purple) to 100% (yellow).

Figure 7 .
Figure 7. Left: Interactive selection on pattern spectrum and compactness) of the DTM with the max-tree.Right: Structures of the DTM corresponding to the selection are displayed with a green overlay on top of the hillshade visualization.

Figure 8 .
Figure 8. Activation map of 3D pattern spectrum (area, height and compactness).The bin color represents the percentage of nodes defined as gold panning in the global spectrum.Values range from 0% (purple) to 100% (yellow).

Table 3 .
). 2D spectra intersection of gold panning and background with selected attributes (lower is better).

Table 4 .
3D spectra intersection of gold panning and background with selected attributes (lower is better).
The results can be displayed in different ways: