BROCELIANDE: a comparative study of attribute profiles and featureprofiles from different attributes

: Morphological attribute proﬁles (APs) are among the most effective spatial-spectral methods to perform multilevel image description based on hierarchical tree-based representation. They have been widely applied to the processing and characterization of remote sensing images, in particular to tackle classiﬁcation task, in the literature. Recently, a novel extension of APs called FPs has been proposed by replacing pixel gray-levels with some statistical and geometrical features when forming the output proﬁles. FPs have been proved to be more efﬁcient than the standard APs when generated from both inclusion and partition trees. The motivation of this article is to conduct a comparative study of APs and FPs using different attributes including some novel ones that have not been used in the literature. We also present our developed library called Broceliande, which proposes efﬁcient implementation of APs and FPs to perform remote sensing image classiﬁcation, with various choices of tree structures as well as attributes. We perform our experiments on two high resolution optical image data sets and provide comparative results of APs and FPs, showing and conﬁrming their effectiveness to describe and classify remote sensing images.


INTRODUCTION
Earth observation has become essential for activities connected with the environment, whether to understand it or to preserve it. One of the most efficient manners remains the use of remote sensing imagery. However, due to the increasing amount of data, we have to find both fastest processing and relevant results. Among a great number of spatial-spectral approaches in the literature, morphological attribute profiles (APs) (Dalla Mura et al., 2010, Pham et al., 2018b have been widely used to process and describe remote sensing images in the literature. The reason relies on their powerful multilevel modeling of the image content and their efficient implementation via tree representation. In fact, the core idea of APs is to explore hierarchical information from a sequence of filtered versions of an image. These images are obtained by the application of different filter rules (called attributes) characterizing the size and shape of objects present in the image. Recently, a novel extension of APs called feature profiles (FPs) (Pham et al., 2017a) has been proposed by replacing pixel gray-levels with the attributes themselves when forming the output profiles. FPs have been proved to be more efficient than the standard APs and the authors showed that APs are indeed a part of the general FPs, in case that the gray-level is used as attribute to form the output profiles. A recent paper has investigated the performance of APs and FPs generated from different kinds of tree structures (min-tree, max-tree, tree of shapes and partition trees such as α-tree or ω-tree) (Pham et al., 2018a). Another one has studied the filtering rules (Das et al., 2020) during AP construction. Nevertheless, tree formation and pruning are not the only factors that affect the performance of APs and FPs. The choice of attributes is even more important in most cases since they decide how the tree should be pruned as well as which output profiles would be formed.
In this article, our motivation is to conduct a comparative study of APs and FPs using different attributes to tackle the classification of optical remote sensing images. Our analysis could help future AP and FP users to have good selection of attributes in order to to perform their multilevel characterization and processing of remote sensing images. In the remainder of this paper, Section 2 briefly summarizes some backgrounds of hierarchical representation using APs and FPs; then presents our developed libraries including TRISKELE and Broceliande as well as the efficient implementation of some promising attributes proposed by our libraries. In section 3, we present our experimental study to perform and compare the performance of APs and FPs using different types of attributes on two remote sensing data sets including one panchromatic image acquired by the IKONOS satellite at 1-m resolution, and one large-scale multispectral image acquired by the Pleiades sensor at 2-m resolution. Section 4 finally concludes the paper and discusses some future works.

Hierarchical representation with APs and FPs
Component trees are structured multilevel representations of the components/regions of an image. They are usually divided into two categories that are both explored in this research: inclusion trees and partition trees. Among the inclusion trees, we use the well-established max-tree and min-tree, and the newly proposed median-tree, which aims to approximate a tree-of-shapes through a linear time complexity algorithm. Hierarchical representation is particularly suited for handling remote sensing images, as attested by the vast research on this field. For instance, component trees are the basis for computing the morphological attribute profiles (APs), which have been successfully applied to the classification of high resolution remote sensing images (Dalla Mura et al., 2010, Pham et al., 2018b. For a reminder, APs are multilevel image description tools obtained by successively applying a set of morphological attribute filters (AFs). Given a grayscale image X : E → Z, E ⊆ Z 2 , the standard generation of APs on X is achieved by applying a sequence of AFs {φ k } K k=1 based on a tree model (which could be a min-tree, max-tree, tree of shapes, median-tree or other partition trees). The AP descriptor of each pixel p in the definition domain of X is written as: AP (p) = X(p), φ1(X) (p), φ2(X) (p) . . . , φK (X) (p) (1) where φ k (X) is the filtered image obtained by applying the attribute filtering φ with regard to the threshold k.
Pruning a tree to take the node value rather than the pixels it refers to provides an AP. When we perform the same reconstruction, not with level depending on the structure of the tree, but with another attribute value we form a feature profile (FP) (Pham et al., 2017a). Specifically, for each pixel p, AP of p obtained by an arbitrary AF φ k is the gray value X ′ (p), where X ′ = φ k (X) is the image reconstructed from the filtered tree (cf. Eq (1)). Now, let Γp(X) be the connected component (CC) of X containing p and let f be a feature or an attribute, i.e. a function admitting a CC and outputting a real value, to be extracted. The FP of p will be f [Γp(X ′ )]. More formally, the generation of FPs are defined as follows: Unlike in AP technique where only one profile is produced from a pruned tree, several features can be simultaneously extracted and stacked to form the final FP. For more information about the extraction of APs and FPs, readers are invited to read their related papers (Dalla Mura et al., 2010, Pham et al., 2017a.

Developed tools
We have developed tools to enable the hierarchical and multilevel image representation and processing: TRISKELE 1 (Merciol, al., 2017) for tree building and pruning, Broceliande 2 (Merciol, al., 2018) using TRISKELE and Shark Random Forest (Igel et al., 2008) for image classification. The aforementioned component trees are efficiently implemented in TRISKELE and can be used to build real-time user interfaces. All attributes managed by Broceliande can be used both for tree pruning and for image reconstruction from a pruned tree.
Since the first utilization of APs, the standard attributes exploited in the filtering step have been area, moment of inertia and standard deviation (Dalla Mura et al., 2010, Pham et al., 2017b, Pham et al., 2018b. Those attributes are suitable for the processing of large-scale high resolution (HR) images since they can be computed in one pass over the nodes of a component tree. Keeping the software up-to-date requires adding the attributes according to user requests. Therefore, some attributes including perimeter, compactness, and rectangularity have been recently added to our library. All of these attributes are computed in linear time on a component tree.

Efficient implementation
The added attributes are introduced with the TRISKELE constraints. To maintain an efficient implementation, all algorithms integrated into the software have linear or quasi-linear time complexity. This subsection explains how we can achieve this goal. It starts with a brief presentation of the hidden data structure with compact and linear information that help us to design attribute algorithms. First, with the moment of inertia attribute, we illustrate the bottom-up approach for each attribute process in order to limit access by only one read and at most two writes per element. Then, with the perimeter attribute, we explain how to take advantage of the hierarchical structure to reduce time by taking into account neighborhood pixels. Lastly, with the standard deviation attribute, we present another way to obtain faster results close to the original definition.

Tree and attributes implementation
First of all, we use the fundamental properties of hierarchical data representation. Consider a gray scale image (left side of Figure 1). This image is a matrix of pixels. We can build different types of trees (min-tree, max-tree, . . . ) to manage the image. We choose to build a max-tree (right side of figure 1).
This tree is a hierarchical representation. All pixels have parent. These parents are called nodes. All nodes (except the root) also have a parent node. This defines an inclusion tree. Because we choose a max-tree, the nodes are organized from light gray (low level of the tree) to dark gray (high level). All nodes can have properties. We represent them with colored bullets in Figure  TRISKELE uses a compact representation of this data. The tree is defined only with an array of parents (see Figure 2). All pixels and nodes are associated with a unique index. The indexes under "pixel cardinality" refer to the pixels from the first row to the last row of the image. Indexes with higher "cardinality of pixels", refer to tree nodes. Thanks to this building process, all nodes are sorted according to their level. So the light-gray nodes are on the right and the darker-gray ones are on the left. Attributes are only defined for nodes. The attributes can be simple scalar (surface, perimeter, gray value) or multivalued (centroid, bounding box), and they are also stored in another array. The position of an attribute in an array gives the position of the associated node. Here are our efficient calculations of new attributes.

MoI computation
The moment of inertia is a morphological attribute. It depends on the centroid distance among the nodes and the child nodes.
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) We first compute barycenter with Then we computed the moment of inertia The issue is that a strict application of this formula involves the processing of nodes from the leaves to the root and, then, it goes back each time to retrieve children values. The computation time impact is significant. The key is therefore to consider each item only once and store its contribution on the fly. The consequence is only one read and only one write per item (pixels and nodes). The complexity is to be O(n).
Because the moment of inertia is used with a random forest, only trends are needed. So we can save time by forgetting the square root step. In that case, the Algorithm 1 actually gives M oI 2 .
We assume that the average position of each node (centroid) is mainly stored in G. All attribute algorithms in TRISKELE are designed with the same approach (bottom up).

Perimeter computation
The perimeter attribute is a good example to see how to take into account the array data structure to implement the tree structure. The perimeter is the set of pixels that define the border of a node. It depends on the connectivity. Figure 3 illustrates that the higher connectivity (C8) increases the number of pixels involved in perimeter. Note that, even if two areas share the same border, the perimeter depends on their convex (or concave) sides. The perimeter of an embedded node is always shorter than the enclosing node. A simple approach to finding the perimeter could be to create a parent map. In this case, we replace the gray level of each pixel with the index value of its parent. We must consider in this map all the pairs of neighbors (depending on the connectivity) to be sure not to miss any pixel border. On the other hand, we have to be sure not to consider a pixel twice to play this role. This simple approach requires an additional map. Its size is the same as the number of pixels in the original image. With huge images (Giga pixels), an index of parent needs at least 32 bits. The complexity of such algorithm is also inappropriate for handling huge images.
Our approach starts with morphological observations. Consider a simple image (a) with only 3 nodes in Figure 4. We assume it is a gray-scale image, but colored here to highlight the nodes. We build an inclusion tree (e.g. max-tree). The orange pixels have the max value, the green ones have the min values and the blue ones in the middle. We also consider some no-data pixels. The hierarchical representation is depicted in (b) in the same figure. Figure 5 represents the pixel borders with C4 connectivity. Pay attention on the fact when a pixel becomes a perimeter pixel in a node, and when it loses this property in a parent node. We now develop these different cases, represented by 6 numbered pixels from the figure: 1 is orange border with blue node but no longer in parent nodes.
2 is orange border with blue node and blue border with green node, but no longer in parent nodes.
3 is border everywhere with no-data.
4 is blue border with green node but no longer in parent nodes. Each pixel joins the tree at a base node. At this moment, it could be surrounded by the same value in flat area and never be considered as a border. Or it may be close to another value. If The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition) Algorithm 2 increasePerimeter (f rom, to) for c ← f rom to to do perimeter(c) ← perimeter(c) + 1 c ← parent(c) end for the value belongs to an undernode (near the leaves), it should simply wait to be joined by this neighbor and never be a border. If the value belongs to a uppernode (close to the root), it will play the role of border until it joins this node. So, the general approach is for each pixel, we consider the higher ancestor of all its neighbors. The algorithm 3 is quasi-linear. Time per pixel depends on efficient way to find ancestor.

Algorithm 3 Compute perimeter
for all np ∈ npSet do max ← ancestor(max, np) end for increaseP erimeter(pp, max) end if end for end if end for The properties of the array 'parents', described in the previous section, help us to implement an efficient method to find ancestors. We just have to compare the index as in Algorithm 4.
The Table 1 gives execution time obtained with the remote sensing image described in the next section. Except for the perimeter attribute, the 13 megapixels of the image are processed in around 10 seconds (depending of attributes). So the throughput is around 1 megapixel/ second. This table and the Figure 6 are produced with the command in our library: channelGenerator src.tif dst.tif -b 3 -f ֒→ < FeatureProfileName > -t Med -a ֒→ Area --thresholds 10 ,100 ,5000 --֒→ auto --time 2.3.4 Alternative SD To keep this efficiency we also propose alternative production of the standard deviation attribute. For each node, the standard deviation depends on its mean gray value (X).
The computation of the mean gray values takes time (see Table  1). If we consider the weight of a node, it also depends on the gray value of all pixels included in this node. The weight is the maximum value (resp. minimum value) of all pixels in a min-tree (resp. max-tree). If we consider median-tree with well balanced values, the weight (distance to the median) could be close to the average of gray values the one area. The attribute becomes: where W is the gray weight of the node accord the type of tree. As you see with the visual representation in Figure 6, the SDW results are close to the SD results. We have not yet compared deeply these two approaches, but we suggest it as future works. That is the reason why we implemented this alternative algorithm, witch saves time and memory of mean gray value evaluation.

AP Area
Rectangularity SD

Data sets
The first data set is a panchromatic image of size 628 × 700 pixels acquired by the IKONOS Earth imaging satellite with 1m resolution in Reykjavik, Iceland. This data consists of six The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B3-2020, 2020 XXIV ISPRS Congress (2020 edition)  Fig. 7. The second data set is a multispectral image composed of 107, 531, 923 pixels acquired by the Pleiades satellite in Podgorica, Montenegro. It is composed of four spectral bands (Blue, Green, Red and Near-Infrared) at 2-m spatial resolution. This image has been partially annotated  and includes a single class of interest: the small wood features (SWF). To perform classification on this data, a subset of background pixels, i.e. pixels that do not belong to the SWF class, was selected automatically. In total, the ground truth is composed of 417, 285 SWF sample pixels and 17, 211 background pixels. Then, 16, 000 ground truth pixels were randomly selected for training a random forest classifier and the remaining pixels were used for testing. Figure 8 presents this image in false color composition: (Near infrered, Red and Green) as RGB.

Evaluation method
Recently, area attribute has been used for the classification of small woody features (Merciol et al., 2019). In this study, we will use a similar data, the high resolution Podgorica image, plus the Reykjavik dataset to compare the different attributes with the FPs. The evaluation on the Reykjavik data is based on the common approach used in the literature: a FP is computed on the data, the feature profiles of the training pixels are used to train a random forest and, then, the classification is performed on the entire ground truth. Regarding the classification of the small woody features, here is our evaluation process: • We start by checking the consistency of the classification samples. The set of reference values is separated: (A) for training and (B) for prediction. We eliminate negative predictions in (B). They are due to an imprecision of the contours of the reference labeling. This provides us a cleaned subset called (C).
• We compute the accuracy for the SWF class of each AP and FP pair. The set (C) is randomly separated: (D) for the training and (E) for prediction. The quality of prediction corresponds to the accuracy between (E) predicted and the reference.
• We use the same input pixels set (E) for all tests. Only the attributes for the cut and for output FPs are different from one test to another.
Thanks to the partnership (IRISA, SIRS-CLS), we drive our study with the real industrial images. The result will be used to improve actual production. Specific hardware requirement is not necessary since the accuracy depends only on the choice of attributes. We provide computation time only to validate the real-time property of our approach.

Experimental setup
In this section, we explore APs and FPs for the classification of the remote sensing data described in section 3.1. From each image, we compute profiles using the following set of attributes: A = {area, compactness (Li et al., 2013), moment of inertia (of a region) (Li et al., 2013), rectangularity (Rosin, 2003), square of the standard deviation of pixel gray levels (Kumar, Gupta, 2012), perimeter}.
For most pairs (A1, A2) of attributes in the set A, we computed a FP using A1 as a filtering rule and A2 for projection. Our goals are (a) to compare the overall accuracy of FPs and APs in the classification of HR images and (b) to investigate the adequacy of each pair of attributes in this context.
The threshold values adopted in the filtering stage were selected manually. For the area filtering, we considered the same set of threshold values used in the computation of FPs in (Pham et al., 2017a) for the Pavia data, and the threshold values used in (Merciol et al., 2019) for the Podgorica data: λa,Rey = {25, 100, 500, 1000, 5000, 10000, 20000, 50000, 100000, 150000} λ a,P od = {1000, 2500, 5000, 7500} For the perimeter attribute, we adopted the same threshold values used for area. Then, for the scale invariant attributes (moment of inertia, standard deviation and compactness), we considered the following set of thresholds commonly used in the literature: The random forests trained on the Reykjavik and Podgorica data are composed of 100 and 64 trees, respectively. The experiments were performed on a machine with Intel Xeon Gold 6136 CPU, 3.00GHz and 263 GB of memory.

Results and discussion
Tables 2 and 3 present the classification scores and execution times for the Reykjavik and the Podgorica data, respectively.
In Table 2, we present the overall accuracy (OA), the average accuracy (AA) over all classes, and the kappa coefficient (κ) of the classification of the Reyjavik data using random forests. We also present the learning (train) and prediction times using this approach. Our baseline is the classification based solely on the panchromatic values (Pan) of this data. As mentioned previously, we evaluate the FPs based on different combinations of attributes in the filtering and projection steps. We can see that all tested APs and FPs outperform the baseline approach. In most cases, except for the FPs based on the 'moment of inertia' filtering, using other attributes for projection (instead of gray values) provides significant improvements. This outcome has already been discussed in (Pham et al., 2018a) with respect to the attributes area, moment of inertia and standard deviation.
Here, we show that this is also the case for rectangularity, perimeter and compactness. For example, the FP computed using the compactness attribute both for filtering and for projection outperforms the compactness AP by 7.7%, 13.62% and 10.85% in terms of OA, AA and κ, respectively.
Regarding our evaluation on the Podgorica data, our two baselines are the classification based only on the spectral data (composed of four spectral bands), and on the spectral information concatenated to the NDVI channel and to the Sobel gradient of each band. For each method, we show in Table 3 the precision and recall with respect to the SWF class. While that including the NDVI and Sobel bands provides little improvement in recall, using APs and FPs improves the recall by up to 4.66%. Similar to the Reykjavik data, the highest scores are achieved using the FPs obtained with area and perimeter filtering. On the other hand, FPs with 'standard deviation' filtering does not give as good results on this data when compared to Reykjavik. The later may be due to the choice of threshold values, which is left for future research.

CONCLUSION AND FUTURE WORKS
We have conducted a comparative study of attribute profiles and features profiles using different attributes (from which some new ones have been proposed) applied to remote sensing image classification task. Results obtained on both 1-m panchromatic Reykjavik image (IKONOS sensor) and 2-m multispectral Podgorica image (Pleiades sensor) have confirmed the effectiveness of both APs and in particular FPs in terms of classification performance as well as computation time. As a conclusion, FPs are more flexible than APs with their several output options. The choice of attribute for tree pruning plays a crucial role in such an approach. Area attribute has been proved to be the best one in our work as well as in several literature studies, while the newly proposed perimeter and compactness could provide other alternative options. All experiments have been conducted using our two libraries, TRISKELE and Broceliande, which are publicly available for the research community.
As mentioned in Section 2.3.4, our libraries include both SDW and SD as options to compute the standard deviation. Since the computation of SDW is much more efficient than SD, one of our future works could be to conduct a deep comparative study on their qualitative performance to determine if SDW could permanently replace SD in future use of this attribute. Another work is to integrate the APs and FPs without thresholds (threshold-free) into the libraries as well as make them available to other types of remote sensing data such as time series, Radar and Lidar data.