*2.3. Methods*

2.3.1. Absorption Wavelength Mapping

To map absorption features, namely absorption wavelengths and depths, new tools, which are described here for the very first time, were programmed to allow automatic detection of multiple absorption feature parameters. The tools—e.g., called QUANTools—have been created using IDL programming language (ENVI/IDL: version 5.0 and higher, [62]).

The tools can process both spectral libraries and hyperspectral image data. They consist of basic modules (GUIs) (Figure 2) allowing users to:


The only decisions made by an operator/expert are to define the spectral region (Figure 2, step 1) and the number of desired absorption features mapped in each spectral region (Figure 2, step 4). All the other steps are done automatically when processing the datasets. A detailed explanation of the processing steps available in the toolbox is given in Section 2.3.2, below.

#### 2.3.2. The Toolbox (QuanTools) Description

QuanTools can be used for spectral absorption band mapping using spectral libraries or hyperspectral image data (Figure 2). A user first defines the spectral range to be analysed within the visible (VIS), near-infrared (NIR), shortwave infrared (SWIR) or longwave infrared (LWIR) regions (Figure 2, step 1). Different spectral ranges can be defined and analysed after one another. Next comes continuum removal (CR, [63]), by which the continuum—that is, a convex hull of straight-line segments, is fitted over a reflectance spectrum and subsequently removed by division or rationing (Figure 2, step 2). The next step is noise detection and removal using local minima of the spectral curve defined as:

$$R = \left\{ \begin{array}{l} 1 \ if \rho\_{\mathcal{CR}}(\lambda - 1) \ge \rho\_{\mathcal{CR}}(\lambda) \\ 0 \text{ if } \rho\_{\mathcal{CR}}(\lambda - 1) < \rho\_{\mathcal{CR}}(\lambda) \end{array} \right\} \tag{1}$$

$$\inf\left[R(\lambda - 1) = 1\right] \land \left[R(\lambda) = 0\right] \Rightarrow \text{local minimum},\tag{2}$$

where *ρCR*(*λ*) and *ρCR*(*<sup>λ</sup>* − 1) are the values of the spectral curve after Continuum Removal and *R* is an auxiliary variable. When the local minima are found, it has to be decided whether the value represents noise or not. The decision is based on two sets of values (a defined number of surrounding bands, e.g., called a spectral window)—derived from bands placed before and after the local minimum:

$$A(\lambda) = \begin{Bmatrix} \rho\_{\mathcal{CR}}(\lambda - nn); (\lambda - nn) \\ \vdots \\ \rho\_{\mathcal{CR}}(\lambda - 2); (\lambda - 2) \\ \rho\_{\mathcal{CR}}(\lambda - 1); (\lambda - 1) \end{Bmatrix} \tag{3}$$

$$\rho\_c(\lambda) = \left\{ \begin{array}{c} \rho\_{\mathcal{CR}}(\lambda+1); (\lambda+1) \\ \rho\_{\mathcal{CR}}(\lambda+2); (\lambda+2) \\ \vdots \\ \rho\_{\mathcal{CR}}(\lambda+nn); (\lambda+nn) \end{array} \right\},\tag{4}$$

where *A*(*λ*) represents the set of surrounding bands before the local minimum, *B*(*λ*) represents the set after the local minimum, and *nn* is the number of surrounding bands. Then the minimum values of these sets are detected as:

$$A\_{\min}(\lambda) = \min \rho\_{\mathbb{C}\mathcal{R}} \in A(\lambda) \tag{5}$$

$$B\_{\min}(\lambda) = \min \rho\_{\text{CR}} \in B(\lambda). \tag{6}$$

The minimum of the spectral window (*Amin*(*λ*), *Bmin*(*λ*)) is compared to the local minimum *ρCR*(*λ*). Using this comparison, the noise is detected as follows:

$$\left[B\_{\min}(\lambda) < \rho\_{\mathcal{CR}}(\lambda)\right] \land \left[A\_{\min}(\lambda) > \rho\_{\mathcal{CR}}(\lambda)\right] :\tag{7}$$

$$if\rho\_{CR} \in A(\lambda) < \rho\_{CR\max} \in B(\lambda) \Rightarrow \text{print}\left\{\begin{array}{l} \text{Err}(\lambda - 1) \\ \text{Err}(\lambda - 2) \\ \vdots \\ \text{Err}(\lambda - n) \end{array} \right\}; n \leq mn \\\\ if\left\{ \begin{array}{l} (\lambda + 1) \\ (\lambda + 2) \\ \vdots \\ (\lambda + m) \end{array} \right\} \in B(\lambda) < \rho\_{CR\max}(\lambda) \in B(\lambda) \Rightarrow \text{print}\left\{ \begin{array}{l} \text{Err}(\lambda + 1) \\ \text{Err}(\lambda + 2) \\ \vdots \\ \text{Err}(\lambda + n) \end{array} \right\}; n \leq mn \\\\ \left[A\_{\text{min}}(\lambda) < \rho\_{CR\max} \in A(\lambda) \Rightarrow \rho\_{CR\max}(\lambda)\right] : \\\\ if\rho\_{CR} \in B(\lambda) < \rho\_{CR\max} \in A(\lambda) \Rightarrow \text{print}\left\{ \begin{array}{l} \text{Err}(\lambda + 1) \\ \text{Err}(\lambda + 2) \\ \vdots \\ \text{Err}(\lambda + n) \end{array} \right\}; n \leq mn \\\\ if\left\{ \begin{array}{l} (\lambda - 1) \\ (\lambda - 2) \\ \vdots \\ (\lambda - m) \end{array} \right\} \in A(\lambda) < \rho\_{CR\max}(\lambda) \in A(\lambda) \Rightarrow \text{print}\left\{ \begin{array}{l} \text{Err}(\lambda - 1) \\ \text{Err}(\lambda - 2) \\ \vdots \\ \text{Err}(\lambda - n) \end{array} \right\}; n \leq mn,\end{aligned}$$

where *Err*(*λ*) represents noise ("bad" bands).

These error values are recalculated using the values of the surrounding bands (their number is defined by a spectral window size). For sample spectral data, a graphical interphase can be used to set up a spectral window size and check the corrected spectral curves (Figure 3, step 3). This allows the user to compare the results of different settings and finally employ a correction that is optimally tuned for the data under analysis.

After a noise-cleaned image is retained, the operator decides how many absorption features are to be derived (Figure 2, step 4) and the analysis is done automatically (Figure 2, step 5) using the following processing:

$$
\mathbb{C}Rdepth(\lambda) = 1 - \rho\_{\mathbb{C}R}(\lambda). \tag{9}
$$

Then the trend of the spectrum is analysed in a similar way to noise detection. First the saddle points (the local absorption maxima) are detected using these conditions:

$$R\_{i\bar{j}} = \left\{ \begin{array}{l} 1 \text{ if } \mathsf{CRdepth}(\lambda)\_{i\bar{j}} \ge \mathsf{CRdepth}(\lambda - 1)\_{i\bar{j}} \\ 0 \text{ if } \mathsf{CRdepth}(\lambda)\_{i\bar{j}} < \mathsf{CRdepth}(\lambda - 1)\_{i\bar{j}} \end{array} \right\} \tag{10}$$

*Remote Sens.* **2017**, *9*, 1006

$$\operatorname{if}\left[\mathcal{R}(\lambda - 1)\_{i\bar{j}} = 1\right] \land \left[\mathcal{R}(\lambda)\_{i\bar{j}} = 0\right] \Rightarrow \operatorname{print}\left[\operatorname{Loc\\_max}\_{i\bar{j}}\right]$$

where *R* is an auxiliary variable that records increase (1) or decrease (0). During this step, local absorption maximum depths and wavelengths are registered:

$$\operatorname{Loc}\_{\operatorname{max}ij} = \left\{ \begin{array}{c} \operatorname{Loc}\_{\operatorname{max}1} \operatorname{depth}(\lambda)\_{ij} \operatorname{Loc}\_{\operatorname{max}1}(\lambda)\_{ij} \\ \operatorname{Loc}\_{\operatorname{max}2} \operatorname{depth}(\lambda)\_{ij} \operatorname{Loc}\_{\operatorname{max}2}(\lambda)\_{ij} \\ \vdots \\ \operatorname{Loc}\_{\operatorname{max}max} \operatorname{depth}(\lambda)\_{ij} \operatorname{Loc}\_{\operatorname{max}nn}(\lambda)\_{ij} \end{array} \right\}. \tag{11}$$

A set number of desired absorption features is detected as the most pronounced local absorption maxima (a set number of the most pronounced absorptions according to the absorption depths). If the desired number is larger than the number of local maxima detected, the final array is completed by zero values:

$$\operatorname{Loc}\_{\operatorname{max}ij} = \left\{ \begin{array}{c} \operatorname{Loc}\_{\operatorname{max}1} \operatorname{depth}(\lambda)\_{ij} \operatorname{Loc}\_{\operatorname{max}1}(\lambda)\_{ij} \\ \operatorname{Loc}\_{\operatorname{max}2} \operatorname{depth}(\lambda)\_{ij} \operatorname{Loc}\_{\operatorname{max}2}(\lambda)\_{ij} \\ \vdots \\ \operatorname{Loc}\_{\operatorname{max}Tol\_{\operatorname{am}}} \operatorname{depth}(\lambda)\_{ij} \operatorname{Loc}\_{\operatorname{max}n}(\lambda)\_{ij} \\ \operatorname{0,0} \\ \vdots \\ \vdots \\ \operatorname{0}\_{\operatorname{min}} \operatorname{0}\_{\operatorname{min}} \end{array} \right\}. \tag{12}$$

Finally, the local absorption maximum features are assigned to image matrices. For *z* absorption features, two matrices with *z* bands are created; one has the absorption wavelengths assigned, while the second one has corresponding absorption depths assigned, respectively.

**Figure 2.** QUANTools: a simplified processing scheme (CR: continuum removal, λ : wavelength loc\_mac: local absorption feature maximum). The detected absorption feature parameters are sorted in ascending order by wavelength (λ1, λ2 − <sup>λ</sup>*z*).

**Figure 3.** Bad band detection (spectral range 0.450–1.200 μm) in QuanTools where the bad bands detected are shown by the green square points and the corrected curve as a dashed line (set spectral window: three neighbouring bands).

#### 2.3.3. Specific Setting of the QUANTools

The main aim was to test if the major absorption features present though the VIS/NIR/SWIR and LWIR ranges can be (i) derived and integrated into one raster dataset; or (ii) successfully used for a final mineral mapping. The spectral ranges from which the absorption feature parameters were derived were defined as follows: VIS/NIR (HyMap): 0.450–1.200 μm, SWIR (HyMap): 2.100–2.400 μm and LWIR (AHS): 8.500–12.500 μm. Noise detection and correction was employed to the VIS/NIR reflectance (0.450–1.200 μm) and the spectral window was set to 3 (three spectral bands before and after the local minimum were analysed). The number of absorption features to be mapped in each spectral region was set to 2. As a result, two new raster datasets were derived for each spectral range, one having assigned the two major absorption wavelengths and the second raster dataset with corresponding absorption depths.

2.3.4. Integration of Absorption Feature Information Detected In VIS/NIR/SWIR and LWIR Data and Further Classification

After the wavelength mapping was employed to different spectral ranges, it was necessary to find a way to further integrate the absorption feature mapping results (Figure 4), more specifically the absorption wavelength and depth matrices. The MNF transformation [64] was employed to compress the data variability of these image matrices (six raster absorption wavelength matrices and six raster absorption depth matrices).

The intention was to test what the advantage will be of adding LWIR data (AHS) to the further mineral classification compared to using only the VIS/NIR/SWIR data (HyMap). Therefore, the MNF transformation was employed for two different scenarios:


For each scenario, the first three MNF bands were visually analysed further to identify training areas (ROIs) representing different material/surfaces. The ROIs were defined as representative pixels of different colour clusters/regions. These were easy to identify when MNF1/MNF2/MNF3 images were displayed as RGB (Figure 5). The ROIs were then used for further supervised classification. In this case the simple non-parametric supervised (parallelepiped) classification was employed [65–67], which uses a simple decision rule to classify multispectral data. The parallelepiped classifier uses the threshold of each class signature to determine if a given pixel falls within the class or not. The thresholds specify the dimensions (in standard deviation units) of each side of a parallelepiped surrounding the mean of the class in the feature space. The best results (e.g., class separability, lowest number of unclassified pixels) were achieved when the upper and lower limits of each parallelepiped were set to ±1.5 standard deviations. In the last step, the classifications were overlaid over the original HyMap and AHS image data and the average spectrum from the original HyMap and AHS image data was computed for each class to (i) ensure that each class represents different surface material (is represented by a unique spectral signature) and (ii) be able to describe the mineral composition of each class when interpreting the spectral property.

**Figure 4.** Scheme of the processing allowing integration of the absorption feature information detected in VIS/NIR/SWIR and LWIR data and further classification.

**Figure 5.** Lítov: MNF1/MNF2/MNF3 images displayed as RGB: ( **A**) Scenario 1 (absorption feature parameters derived only from the HyMap data were used); (**B**) Scenario 2 (absorption feature parameters derived from both datasets, HyMap and AHS, were used), the MNF bands were further visually analysed to identify training areas (ROIs) representing different material/surfaces.
