Cedrigue Boris

Djiongo Kenfack

Olivier Monga
email: olivier.monga@ird.fr

Serge Moto Mpong

René Ndoundam

Quaternion-Based Texture Analysis of Multiband Satellite Images: Application to the Estimation of Aboveground Biomass in the East Region of Cameroon

Keywords: Discrete quaternion Fourier transform, Quaternion-based texture indices, Multiband satellite images, Aboveground biomass, Cameroon

Within the last decade, several approaches using quaternion numbers to handle and model multiband images in a holistic manner were introduced. The quaternion Fourier transform can be efficiently used to model texture in multidimensional data such as color images. For practical application, multispectral satellite data appear as a primary source for measuring past trends and monitoring changes in forest carbon stocks. In this work, we propose a texture-color descriptor based on the quaternion Fourier transform to extract relevant information from multiband satellite images. We propose a new multiband image texture model extraction, called FOTO++, in order to address biomass estimation issues. The first stage consists in removing noise from the multispectral data while preserving the edges of canopies. Afterward, color texture descriptors are extracted thanks to a discrete form of the quaternion Fourier transform, and finally the support vector regression method is used to deduce biomass estimation from texture indices. Our texture features are modeled using a vector composed with the radial spectrum coming from the amplitude of the quaternion Fourier transform. We conduct several experiments in order to study the sensitivity of our model to acquisition parameters. We also assess its performance both on synthetic images and on real multispectral images of Cameroonian forest. The results show that our model is more robust to acquisition parameters than the classical Fourier Texture Ordination model (FOTO). Our scheme is also more accurate for aboveground biomass estimation. We stress that a similar methodology could be implemented using quaternion wavelets. These results highlight the potential of the quaternion-based approach to study multispectral satellite images.

Introduction

The growth of emissions of carbon dioxide constitutes a prime contributor to global warming, which is nowadays a key ecological issue. Natural factors, particularly those linked to the forests, can play an important role to reduce these emissions. The main reason is that forests constitute one of the biggest carbon pools, which explains that mapping and monitoring biomass in tropical regions have attracted a great deal of attention within past decades [START_REF] Goetz | Mapping and monitoring carbon stocks, with satellite observations: a comparison of methods[END_REF][START_REF] Verhegghen | Mapping Congo basin vegetation types from 300 m and 1 km multi-sensor time series for carbon stocks and forest areas estimation[END_REF]. Therefore, estimating the amount of forest biomass becomes nowadays a major environmental challenge. Indeed, it is still an open issue for the tropical and the subtropical regions. Forest biomass represents living mass such as trees, roots, shrubs, and the dead mass of fine and coarse litter associated with soil [START_REF] Dengsheng | The potential and challenge of remote sensing-based biomass estimation[END_REF]). Forest biomass is often split in two parts: the aboveground and the belowground. As most research works coping with forest biomass, the present study deals with aboveground biomass estimation. Traditional methods (such as destructive sampling and allometric equations) based on forest survey data constitute the most accurate way to estimate the biomass of a specific area, but they remain time-consuming and labor intensive [START_REF] Ketterings | Reducing uncertainty in the use of allometric biomass equations for predicting above-ground tree biomass in mixed secondary forests[END_REF][START_REF] Chave | Error propagation and scaling for tropical forest biomass estimates[END_REF][START_REF] Chave | Improved allometric models to estimate the above ground biomass of tropical trees[END_REF][START_REF] Basuki | Allometric equations for estimating the above-ground biomass in tropical lowland Dipterocarp forests[END_REF][START_REF] Saatchi | Benchmark map of forest carbon stocks in tropical regions across three continents[END_REF][START_REF] Feldpausch | Tree height integrated into pantropical forest biomass estimates[END_REF][START_REF] Vieilledent | A universal approach to estimate biomass and carbon stock in tropical forests using generic allometric models[END_REF][START_REF] Fayolle | Tree allometry in Central Africa: testing the validity of pantropical multi-species allometric equations for estimating biomass and carbon stocks[END_REF][START_REF] Hunter | Tree height and tropical forest biomass estimation[END_REF][START_REF] Kearsley | Conventional tree height-diameter relationships significantly overestimate aboveground carbon stocks in the Central Congo Basin[END_REF][START_REF] Mitchard | Uncertainty in the spatial distribution of tropical forest biomass: a comparison of pan-tropical maps[END_REF][START_REF] Ekoungoulou | Tree allometry in tropical forest of Congo for carbon stocks estimation in aboveground biomass[END_REF][START_REF] Ekoungoulou | Evaluating the carbon stock in above-and below-ground biomass in a moist Central African forest[END_REF][START_REF] Picard | Reducing the error in biomass estimates strongly depends on model selection[END_REF][START_REF] Gaia | Aboveground biomass estimation in an African tropical forest with lidar and hyperspectral data[END_REF][START_REF] Re | Tree allometric equations in mixed forest plantations for the restoration of seasonal semi deciduous forest[END_REF]. However, they fail in yielding accurate spatial distributions of biomass for large extents. Within this context, remote sensing represents an opportunity to fix this task and to monitor large forests attributes with reasonable accuracy levels. The data provided by remote sensing instruments appear therefore as a primary source for measuring past trends and monitoring forest biomass changes. Two main approaches are available to extract relevant information from remote sensing data for aboveground biomass estimation: the reflectance-based approach and the image texture-based approach. The main idea of the reflectance-based approach consists in establishing a relationship between reflectance values of satellite images and forest structure parameters. Hence, [START_REF] Roy | Biomass estimation using satellite remote sensing data an investigation on possible approaches for natural forest[END_REF] proposed a linear regression model aiming to establish a linear relationship between biomass and vegetation indices, such as the wetness or the Normalized Difference Vegetation Index (NDVI), derived from satellite images. [START_REF] Thenkabail | Hyperion, IKONOS, ALI, and ETM+ sensors in the study of African rainforests[END_REF] assessed the ability to predict biomass and other structural stand parameters for the Cameroonian rainforest from vegetation indices derived from hyperspectral and hyperspatial satellite images. They reached to the conclusion that vegetation indices from hyperspatial satellite images, such as IKONOS or Landsat ETM+, underestimate the biomass. Hence, these vegetation indices provided by reflection-based approaches are not well-adapted to derive accurate biomass estimations in rainforests, where higher biomass values are encountered. The biomass estimation model proposed by [START_REF] Gong | Estimation of forest aboveground biomass in California using canopy height and leaf area index estimated from satellite data[END_REF] constitutes an alternative to address this issue. It is a simple linear model of prediction that combines the Leaf Area Index (LAI), extracted from high spatial resolution images of the Landsat Thematic Mapper sensor, with canopy maximum height information provided by the Geoscience Laser Altimeter System sensor. This simple parametric model provided good results in the case of the Canadian forest, but is found to be not scale-invariant. One major limitation encountered in the reflectance-based analysis of satellite image data is due to the fact that a reflectance is always saturated when a forest has a high value of biomass. Texture-based approaches constitute an interesting alternative to deal with this issue. These methods rely on the assumption that estimating biomass from satellite data can be tackled by image texture classification. Most texture-based approaches consist in extracting canopy properties of trees from one spectral band. A spectral band represents the electromagnetic energy reflected or emitted by the vegetation for a limited range of wavelengths. One of the most popular methods to address aboveground forest biomass estimation is the FOurier-based Texture Ordination (FOTO) [START_REF] Couteron | Quantifying change in patterned semi-arid vegetation by Fourier analysis of digitized aerial photographs[END_REF][START_REF] Couteron | Predicting tropical forest stand structure parameters from Fourier transform of very high-resolution remotely sensed canopy images[END_REF]. The basic principle of FOTO is to link canopy granularity to image texture features derived from Fourier spectrum analysis. It allows ordering canopy images according to a monochromatic (panchromatic satellite data) image coarseness-finesses texture gradient in a way congruent with the visual appraisal [START_REF] Couteron | Predicting tropical forest stand structure parameters from Fourier transform of very high-resolution remotely sensed canopy images[END_REF][START_REF] Barbier | The variation of apparent crown size and canopy heterogeneity across lowland Amazonian forest[END_REF][START_REF] Barbier | Linking canopy images to forest structural parameters: potential of a modeling framework[END_REF][START_REF] Ploton | Assessing aboveground tropical forest biomass using Google Earth canopy images[END_REF][START_REF] Proisy | Predicting and mapping mangrove biomass from canopy grain analysis using Fourier-based textural ordination of IKONOS images[END_REF]. Monochromatic data, as those processed by FOTO, are usually the panchromatic products yielded by satellite sensors. The panchromatic band is a single band recorded over a large spectral band. Nevertheless, satellite images often have more than one spectral band. Thus, using the total multiband information can improve the assessment of biomass estimation and other structural forest parameters, such as quadratic mean tree diameter or tree density. Despite of numerous works dealing with multiband image processing, only few methods based on multispectral satellite image analysis have been proposed to address biomass issues.

In recent years, the first approaches using quaternion numbers to handle and model multispectral images in a holistic manner as a compact entity were introduced [START_REF] Sangwine | Fourier transforms of colour images using quaternions or hypercomplex, numbers[END_REF]; [START_REF] Shi | Quaternion colour texture[END_REF]. Thanks to the use of the quaternion Fourier transform, up to four bands of multispectral satellite images can be efficiently and easily processed.

In this work, we are dealing with multispectral satellite data acquired using passive sensors. Our method, called FOTO++, starts by removing noise from multispectral data while preserving edges. Then, color texture indices are extracted using the discrete form of the quaternion Fourier transform, and finally the support vector regression method (SVR) is used to derive the biomass estimation.

The paper is organized as follows. Section 2 presents the texture-based image analysis approach to estimate biomass called FOTO, which has been introduced by [START_REF] Couteron | Quantifying change in patterned semi-arid vegetation by Fourier analysis of digitized aerial photographs[END_REF]. Section 3 deals with quaternion algebra tools which are useful for multiband image processing. We also present a brief review of color texture analysis. The fourth section is devoted to the FOTO++ model description. The fifth section investigates the validation of FOTO++ on synthetic and real data. Section 6 gives future perspectives.

Texture-Based Analysis Image for Forest Biomass Estimation: The FOTO Method

FOurier-based Texture Ordination (FOTO) is an aboveground biomass estimation method performed on data obtained from monochromatic satellite canopy image texture analysis. The basic principle is to link canopy grain to image texture features derived from Fourier spectrum analysis. This section deals with the presentation of the method and its limitations.

Description of the Method

The FOTO method proposed by [START_REF] Couteron | Quantifying change in patterned semi-arid vegetation by Fourier analysis of digitized aerial photographs[END_REF] orders canopy images according to a coarseness-fineness texture gradient, in a way congruent with the visual appraisal [START_REF] Couteron | Predicting tropical forest stand structure parameters from Fourier transform of very high-resolution remotely sensed canopy images[END_REF][START_REF] Barbier | The variation of apparent crown size and canopy heterogeneity across lowland Amazonian forest[END_REF]. It shows promising results for the characterization and measurement of canopies and for biomass estimation in tropical rainforests [START_REF] Proisy | Predicting and mapping mangrove biomass from canopy grain analysis using Fourier-based textural ordination of IKONOS images[END_REF]. The FOTO method is summarized in Fig. 1 and detailed below.

• Image splitting The first step consists in splitting the initial image into subimages representing a sample unit measured in the field. For instance, it can be one hectare as in this work. This process depends on the image size and the splitting size. • Feature extraction Afterward, for each sub-image, 2D Fourier transform periodograms are calculated to characterize textural image properties. It takes into account the spatial frequency information contained in each sub-image and provides the so-called R-spectrum. The R-spectrum is taken from the polar representation of the sub-image's Fourier coefficients by averaging the coefficients having the same radius value in the polar space. For each sub-image, the R-spectrum represents the breakdown of variance accounted for by successive intervals of spatial frequencies [START_REF] Barbier | The variation of apparent crown size and canopy heterogeneity across lowland Amazonian forest[END_REF]. The R-spectra values are standardized and ordinated using Principal Component Analysis (PCA) in order to iden-Fig. 1 The FOTO method algorithm tify the main gradients of the canopy textural variation. The PCA score of the most prominent axes is used as texture index for feature extractions.

Limitation of the Method

The FOTO method has shown promising results in the case of mangrove forests [START_REF] Proisy | Predicting and mapping mangrove biomass from canopy grain analysis using Fourier-based textural ordination of IKONOS images[END_REF], and for linking canopy images to forest structural parameters [START_REF] Barbier | Linking canopy images to forest structural parameters: potential of a modeling framework[END_REF]). Nevertheless, the FOTO method can handle only monochromatic (panchromatic data) images for texture analysis, although multispectral information (e.g. visible colors combined with infrared) can be of great importance. Indeed, multispectral data contain more precise spectral reflectance signatures than panchromatic data. These spectral signatures can be used to detect fundamental vegetation properties such as phenology, canopy structure, vegetation stand structure, and physiological characteristics that link physical properties to ecological theory [START_REF] Ustin | Remote sensing of plant functional types[END_REF]. In addition, the feature extraction process and the classification scheme (regression model) of the FOTO method can be improved by using more robust and powerful non-linear statistical tools [START_REF] Tapamo | Linear versus non-linear methods: a comparative study for forest above ground biomass estimation from texture analysis of satellite image[END_REF]. These limitations can be reduced using advanced techniques of color image processing. Indeed, when the spectral bands are pooled together, they provide one multispectral image that needs to be addressed using multiband image processing methods. In the last decade, some approaches using quaternion numbers to handle and model multiband images were introduced [START_REF] Sangwine | Fourier transforms of colour images using quaternions or hypercomplex, numbers[END_REF][START_REF] Shi | Quaternion colour texture[END_REF]. These numbers allow processing up to four multispectral images in a holistic manner as a compact entity. Furthermore, the choice of an appropriate color space can improve significantly the quality of the image processing methods.

Color Images Processing and Quaternion Algebra

This section presents useful concepts on color image processing and quaternion algebra. We also present a brief state of the art of color image modeling.

Color Spaces

A color is the way that the human vision system (HVS) perceives the light in the visible region of the electromagnetic radiation. The HVS uses colors to sense environment, recognize objects and convey information. These colors are perceived through a set of detectors of which each is sensitive to a specific spectral wavelength of the electromagnetic radiation [START_REF] Fernandez-Maloigne | Digital color imaging[END_REF]. For color sensation, at least two detectors with different sensitivities are needed. The HVS embeds three types of color detectors, named cell cones, that are sensitive to the red (R), the green (G) and the blue (B) wavelengths [START_REF] Noor | Understanding color model: a review[END_REF][START_REF] Shuy | Fundamentals of color[END_REF]. Hence, these three colors constitute the three primary color channels used by the HVS for color representation.

The set of all the colors obtained from the primary ones (red, green and blue) forms the RGB color space. The main drawback of the RGB color space is the high correlation between its components. Other color representations can be derived from the RGB color space using either linear or non-linear transformations. One of them is the HSI (hue-saturation-intensity) color space which is more intuitive to human vision. The HSI system separates specific color information in an image according to its intensities. There exist some variants of this color space such as the HSB (hue-saturation-brightness), the HSL (hue-saturation-lightness) and the HSV (huesaturation-value) color spaces. In these color spaces, the hue and the saturation components represent the color information, while the intensity component describes either the amount of light (the brightness and the value component) or the amount of white (the lightness component) of the image. In the HSI system, the hue component denotes the dominant wavelength in the spectral distribution of light. It represents the basic color such as yellow, blue or black. The saturation component corresponds to a measure of the purity of the color. Also, the Commission Internationale de l'Éclairage (CIE) proposed the Lab color space, which meets the psychological perception of colors by a human observer. Color is described in this space by using three components: L, a and b. The component L denotes the lightness or the intensity of the color while components a and b describe the color information. The main feature of the Lab color space is that it is a perceptually equal space, in which a Euclidian distance between two colors is strongly correlated with the human's perceptual difference of these colors [START_REF] Cheng | Color images segmentation: advances and prospects[END_REF]. It means that the difference between two color coordinates in the Lab color space has the same meaning as in the way a human perceives it. "Appendices 2a and b" respectively present how to compute the Lab and HSV components from the RGB color space. Figure 2 presents the HSV and the Lab color spaces as described in [START_REF] Noor | Understanding color model: a review[END_REF].

Like the HVS, a sensory satellite is commonly used to sense and study environment phenomena. It embeds many detectors, each one being sensitive to one spectral wavelength of the electromagnetic spectrum. These detectors, also called sensors, are able to detect and record the signal intensity over a large range of the electromagnetic spectrum rather than the small visible part that HVS can handle. Satellites produce generally two types of images: panchromatic images and multispectral images. Panchromatic images consist in only one band recorded over a large range of spectral bands with a higher spatial resolution as compared to multispectral images, and is usually viewed as a gray scale image or a monochromatic image. The multispectral images are constituted of several bands of data having a small spectral range with a relatively low spatial resolution. Depending on the chosen spectral bands, multispectral images can be viewed as true color images, pseudo-natural color images or color composite images. True color images are obtained by putting multispectral satellite images in the RGB color space with the red, the green and the blue bands assigned to red, green and blue channels respectively. In the case where the near infrared (NIR) band is used jointly with the red and the green bands and are assigned to the red, the green and the blue channels respectively, we have a pseudo-natural color image. Here, near infrared refers to the near infrared wavelength band. In the case of IKONOS satellite images, the near infrared band is acquired by a sensor within the spectral range of 780-920 nm. In all other cases, images obtained are said to be color composite images. Depending on the specific spectral characteristics of the features under study, spectral bands can be selected and rendered as color images to enhance these features. Figure 3 shows canopy images of a forest in the South-west region of Cameroon, acquired using the IKONOS sensors, represented in RGB and HSV color spaces. The RGB color image provided here is a pseudo-natural color representation of IKONOS images. We get the near infrared spectral band coded in the red channel, the red spectral band coded in the green channel and the green spectral band coded in the blue channel.

In this paper, as we deal with vegetation and know that their reflectance values are greater on the NIR, red and green wavelengths than in other wavelengths, we retained these bands from satellite products and processed the resulting images as pseudo-color images. Furthermore, we choose among the existing color spaces the RGB, HSV and Lab color spaces to represent color. Indeed, none of the proposed spaces dominates the others, for various kinds of color images, so selecting the best one remains a main difficulty in color image processing.

Color and Texture Features Extraction

Image contents are usually described using color and texture features. Texture is an image feature that represents the local spatial organization of spatially varying spectral values that extend to a region of larger spatial scale [START_REF] Chen | The handbook of pattern recognition and computer vision[END_REF]. Thus, the perception of texture is a function of spatial and radiometric scales. Texture analysis is a major step in the texture classification task, which aims to extract relevant attributes for a specific application from color images [START_REF] Fernandez-Maloigne | Digital color imaging[END_REF]. Within the last decades, color has become a necessity to address pattern recognition and computer vision problems. Because acquisition and processing hardware for color images has also become more accessible to deal with the computational complexity of the high-dimensional color space, color features extraction becomes more effective. The simultaneous use of color and texture information as color-texture features can improve significantly the performance of image analysis [START_REF] Luccheseyz | Color image segmentation: a state-of-the-art survey[END_REF][START_REF] Fernandez-Maloigne | Digital color imaging[END_REF]. Classically, three main approaches are used to extract color and texture features: the statistical approach, the stochastic model approach and the spatio-frequential approach. The statistical-based approach defines texture features regarding to its color variation, and gives the relationship between a pixel and its neighbors. The stochastic model-based approaches aim to describe spatial dependencies from a probabilistic point of view by means of a Markov random field. The spatio-frequential approach consists in characterizing the texture in terms of number of transitions per unit area (spatial domain), or in terms of transitions per unit wavelength (frequency domain), or both. It is possible to transform an image into a frequency domain by means of a Fourier transform and manipulate this frequency domain representation for information extraction. The spatio-frequential approaches aim to characterize the color and texture of multiband images in the frequency domain because it is possible to express all the information present in the image using just a small number of coefficients. Within a spatial domain, a color image is represented as a matrix of three-dimensional vectors (R, G, B). However, in the frequency domain a monochromatic image is represented by its spatial frequencies having each a magnitude and a phase. The classical Fourier transform is not well adapted to multispectral images of dimension greater than 2 (color images, multiband images, hyperspectral images), because complex numbers have only two components. The quaternion numbers can be used to overcome this drawback.

The Quaternion Algebra

Quaternion algebra, generally denoted by ℍ , was introduced by Hamilton in 1843 as an extension of complex algebra (Hamilton 1866). Quaternion algebra is a noncommutative algebra (Hamilton 1866; [START_REF] Meister | A concise quaternion geometry of rotations[END_REF][START_REF] Ell | Quaternion Fourier transforms for signal and image processing[END_REF] in which a quaternion q ∈ ℍ is a four-dimensional hypercomplex number having one real part and three imaginary parts. Its Cartesian form is given by where q 0 , q 1 , q 2 , q 3 are real numbers. The three imaginary units i, j, k are square roots of -1 and obey to the following rules:

The Cartesian form of a quaternion represents its expression in a specific 4D basis of the algebra ℍ , namely in the basis {1, i, j, k} . We present below some prop- erties of quaternion numbers that are useful to understand our work.

1. (Scalar and vector parts) Any quaternion q = q 0 + q 1 i + q 2 j + q 3 k can be decomposed into a scalar part S(q) and a vector part V(q) as:

where S(q) = q 0 and V(q) = q 1 i + q 2 j + q 3 k. 2. (Pure quaternion) If the scalar part of a quaternion q is equal to zero, that is S(q) = 0, then q is said to be a pure quaternion. 3. (Conjugate) The conjugate of quaternion q = q 0 + q 1 i + q 2 j + q 3 k, denoted by q , is given by Eq. 3.4: 4. Let p = p 0 + p 1 i + p 2 j + p 3 k and q = q 0 + q 1 i + q 2 j + q 3 k be two quaternion numbers, then:

(3.1) q = q 0 + q 1 i + q 2 j + q 3 k,

(3.2) ⎧ ⎪ ⎨ ⎪ ⎩ i 2 = j 2 = k 2 = -1, ij = -ji = k, jk = -kj = i, ki = -ik = j. (3.3) q = S(q) + V(q), (3.4) q = q 0 -q 1 i -q 2 j -q 3 k.
(a) (Sum) The sum p + q of p and q is: (b) (Product) The product of p and q , denoted as p•q or pq, is:

where × denotes the cross-product:

(Modulus)

The modulus |q| of a quaternion q = q 0 + q 1 i + q 2 j + q 3 k is given by Eq. 3.8:

6. (Norm) The norm ‖q‖ of quaternion q = q 0 + q 1 i + q 2 j + q 3 k is given by expres- sion 3.9:

7. (Unit quaternion) If ‖q‖ = 1 , then q is called a unit quaternion. 8. (Inverse) Any non-zero quaternion has an inverse given by Eq. 3.10:

If q is a unit quaternion, then its inverse is equal to q -1 = q 9. (Orthogonality) Two quaternions p and q are said to be orthogonal, denoted by p ⊥ q, if the product pq is a pure quaternion. 10. (Polar form) In addition to the Cartesian form, any quaternion q ≠ 0 can also be written in the polar form using Euler's formula, as shown in Eq. 3.11:

In Eq. 3.11, |q| is the modulus, μ q a unit vector called axis, and ϕ q is the phase of the quaternion q, and are determined by expression 3.12. From this expression it follows that these only exist for those q for which both S(q) ≠ 0 and V(q) ≠ 0. In "Appendix 1" we present Euler's formula for those q ≠ 0 for which either S(q) = 0 or V(q) = 0.

(3.5) p + q = (p 0 + q 0 ) + (p 1 + q 1 )i + ( p 2 + q 2 ) j + (p 3 + q 3 )k.

(3.6)

p ⋅ q = p 0 q 0 - ( p 1 q 1 + p 2 q 2 + p 3 q 3 ) + p 0 V(q) + q 0 V(p) + V(p) × V(q), (3.7) V(p) × V(q) = ( p 2 q 3 -p 3 q 2 ) i + ( p 3 q 1 -p 1 q 3 ) j + ( p 1 q 2 -p 2 q 1 ) k. (3.8) |q| = √ q 2 0 + q 2 1 + q 2 2 + q 2 3 . (3.9) ‖q‖ = √ q 2 0 + q 2 1 + q 2 2 + q 2 3 .
(3.10)

q -1 = q ‖ ‖ q 2‖ ‖ . (3.11) q = |q| ( cos q + q sin q ) = |q|e q q .
11. (Symplectic form) Another representation of quaternions is obtained using the symplectic form. Consider the basis {1, , , } of ℍ , where μ and η are two orthogonal unit pure quaternions. Then, any quaternion q = a + bi + cj + dk can be written as:

where

The symplectic form is useful for practical calculations of the quaternion Fourier transform using the classical Fourier transform code [START_REF] Sangwine | Fourier transforms of colour images using quaternions or hypercomplex, numbers[END_REF].

Texture-Based Analysis of Color Images for Biomass Estimation: The FOTO++ Method

The FOTO++ method consists in an extension of the FOTO method for handling multispectral data. We propose to represent multiband images by means of quaternion matrices. Indeed, in recent years several approaches using quaternion numbers to process color images in a holistic manner have been introduced. [START_REF] Pei | A novel block truncation coding of color images by using quaternion-moment preserving principle[END_REF] have presented one of the very first works using quaternion numbers to model color images. Powerful tools for multiband images processing have followed: Quaternion Fourier Transform [START_REF] Sangwine | Fourier transforms of colour images using quaternions or hypercomplex, numbers[END_REF][START_REF] Sangwine | Hypercomplex auto-and cross-correlation of color images[END_REF][START_REF] Pei | Efficient implementation of quaternion Fourier transform, convolution, and correlation by 2D complex FFT[END_REF][START_REF] Ell | Hypercomplex Fourier transforms of color images[END_REF], Auto-correlation [START_REF] Sangwine | Hypercomplex auto-and cross-correlation of color images[END_REF][START_REF] Moxey | Hypercomplex correlation techniques for vector images[END_REF], Cross-correlation [START_REF] Sangwine | Hypercomplex auto-and cross-correlation of color images[END_REF], Convolution [START_REF] Pei | Efficient implementation of quaternion Fourier transform, convolution, and correlation by 2D complex FFT[END_REF], Color Texture analysis [START_REF] Shi | Quaternion colour texture[END_REF]. The aim of FOTO++ is

(3.12) ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ |q| = √ q 2 0 + q 2 1 + q 2 2 + q 2 3 , q = q 1 i+q 2 j+q 3 k √ q 2 1 +q 2 2 +q 2 3 , q = arctan ( √ q 2 1 +q 2 2 +q 2 3 q 0 ) ∈ ] - 2 , 2 [ . (3.13) q = ( a ′ + b ′ ) + ( c ′ + d ′ ) , ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ a ′ = S(q), b ′ = -1 2 (V(q) + V(q)), c ′ = -1 2 (V(q) + V(q)), d ′ = -1 2 (V(q) + V(q)).
to model texture in multiband images thanks to Quaternion Fourier Transform. We present within the following subsections the principal outlines of FOTO++.

Image Preprocessing: Nagao-Matsuyama and Median Combined Filter

The Nagao-Matsuyama filter, proposed by [START_REF] Nagao | Edge preserving smoothing[END_REF], is a simple non-linear image filtering method designed to remove image noise. First, this preprocessing technique defines local directional neighborhoods of a pixel for various orientations. Afterward, the pixel's gray level is replaced by the average of the gray levels of the pixels of the neighborhood having the least standard deviation. We use a variant of this method by introducing a median filtering touch. The median filter is a non-linear filter. Its principle consists in calculating the median value within a specified neighborhood in an image and puts this value in the output image. This is accomplished by collecting and sorting the pixel values in the neighborhood of a pixel. The median value lies in the middle of the sorted list. The median filter has many properties such as edge preserving and impulse noise removing. Indeed, we replace the pixel's gray level by the median value of the neighborhood instead of the mean value. Figure 4 presents the neighborhoods used by the Nagao-Matsuyama filter. The main effect of this filter is the reduction of the smoothing of edges. It yields relatively homogeneous regions while preserving edges.

Our combination of median and Nagao filters gives a new filtering scheme with good edge-filtered response while removing noises (see Fig. 5). This non-linear filtering scheme is particularly adapted to satellite images of forests where linear filtering performs generally not well. In this paper we apply on each image band successively a Nagao "median" filter with a mask of size 5 × 5 pixels. This process (median filtering + Nagao filter) is repeated until convergence is reached, as in the classical Nagao filter. Convergence here is defined as idempotency, which means that the result obtained by the process of filtering does not change anymore if the filtering process is applied once more. When applied on a high spatial resolution multiband image, it yields a canopy image with homogeneous regions and well-delimited edges. This operation improves the feature extraction stage of the FOTO++ method. Figure 5 presents a result of this filtering stage on a real multispectral image. The result is an image where features are now well distinguishable.

Multiband Texture-Based Quaternion Features Using Discrete Quaternion

Fourier Transform [START_REF] Pei | A novel block truncation coding of color images by using quaternion-moment preserving principle[END_REF] were among the first to introduce quaternion numbers to model color images in a holistic manner. They proposed to encode the three color components into the vector part of a quaternion number, leaving the scalar part equal to zero. Thus, a color image, disregarding the color space used, is described by a matrix of pure quaternions. This scheme has enabled the development of powerful tools for multiband image processing. Indeed, a color image can then be processed as a single and unique entity using the algebraic properties of the quaternion set ℍ.

Table 1 Quaternion Fourier transforms

The quaternion numbers , 1 and μ 2 are pure unit quaternion square roots of -1, f is a quaternionvalued function, (x, y) represents the coordinates in the spatial domain, and (u, v) the coordinates in the frequential domain. We recall that the quaternionic product is not commutative

Left Right Sandwich Single-axis ∫ ℝ 2 e -(ux+ y) f (x, y)dxdy ∫ ℝ 2 f (x, y)e -(ux+ y) dxdy ∫ ℝ 2 e -ux f (x, y)e -y dxdy Dual-axis ∫ ℝ 2 e -( 1 ux+ 2 y) f (x, y)dxdy ∫ ℝ 2 f (x, y)e -( 1 ux+ 2 y) dxdy Factored ∫ ℝ 2 e -1 ux e -2 y f (x, y)dxdy ∫ ℝ 2 f (x, y)e -1 ux e -2 y dxdy ∫ ℝ 2 e -1 ux f (x, y)e -2 y dxdy
The quaternion Fourier transform treats the signals as vector fields, and generalizes the conventional Fourier transform. Following the classical convention, we present (x, y) as spatial variables, and (u, v) as frequencies. Then, the image of spatial variables under a Fourier transform are interpreted as frequencies. A classical way to extend the two-dimensional Fourier transform to quaternion images consists in replacing the complex number i with a pure quaternion square root of -1, denoted by μ. Due to the non-commutative property of the quaternion multiplication, and to the expression of the quaternion square root of -1, there are eight different types of quaternion Fourier transforms (QFT). Table 1 presents each QFT type. [START_REF] Ell | Quaternion Fourier transforms for signal and image processing[END_REF] proposed an interesting discussion about these quaternion Fourier transforms, and a way to derive one formulation from the others. In the present work, we choose the left single-axis form of the quaternion Fourier transform and its discrete form as proposed by [START_REF] Sangwine | Fourier transforms of colour images using quaternions or hypercomplex, numbers[END_REF] for color image processing. The left singleaxis form of the quaternion Fourier transform was selected because it allows the QFT to be decomposed into two complex Fourier transforms. Under this particular form, the QFT can be performed using existing efficient implementations of the complex fast Fourier transform, as shown in [START_REF] Ell | Quaternion Fourier transforms for signal and image processing[END_REF] and in [START_REF] Sangwine | Quaternion toolbox for Matlab ® , Ver.2 with support for octonions[END_REF], rather than a custom-written code.

In the sequel, we express the equations using the RGB color space, but our scheme can be applied directly for other color spaces (HSV, Lab). In the RGB color space, a color image I is represented as a three-dimensional vector matrix. Therefore, image I is a (M × N) × 3-matrix where M is the number of lines, N refers to the number of columns and 3 is the number of bands. Thus, to each pixel with coordinates (m, n) (with (0 ≤ m ≤ M -1) and (0 ≤ n ≤ N -1) ) of I is attached a three-dimensional color vector (R(m, n), G(m, n), B(m, n)). Let f be a quaternionvalued function of the color image I . For any pixel of I with coordinates (m, n) , its corresponding quaternion representation f(m, n) is given by Eq. 4.1:

where i, j, k satisfy Eq. 3.2, and where R(m, n), G(m, n) and B(m, n) are the color components in the RGB color space.

The quaternion Fourier transform of f is:

with = i+j+k √ 3 .
The quaternion Fourier inverse transform is:

(4.1) f (m, n) = R(m, n)i + G(m, n)j + B(m, n)k, (4.2) F(u, v) = 1 √ MN M-1 ∑ m=0 N-1 ∑ n=0 e -2 ( mu M + nv N ) f (m, n), (4.3) f (m, n) = M-1 ∑ u=0 N-1 ∑ v=0 e 2 ( mu M + nv N ) F(u, v).
Equation 4.2 can be expressed in the symplectic form (see Eq. 3.13) using three orthogonal unit quaternion numbers 1 , 2 and 3 (see Eq. 4.4 for the expression of these orthogonal unit quaternions). The symplectic form allows the use of existing discrete complex Fourier transform codes to implement the discrete quaternion Fourier transform (DQFT).

The sympletic form of f (m, n), using the three orthogonal axes given by Eq. 4.4, is:

f (m, n) = c 1 (m, n) + c 2 (m, n) 2 , with c 1 (m, n) = a ′ + b ′ 1 and c 2 (m, n) = c ′ + d ′ 1 .
For the calculation of a′, b′, c′, d′ one can refer to Eq. 3.13. By using the symplectic form of f(m, n) in Eq. 4.2 it becomes: where Algorithm 1, described below, presents the discrete quaternion Fourier transform (DQFT) obtained using the classical Fourier transform from the symplectic form of the DQFT as proposed by [START_REF] Sangwine | Quaternion toolbox for Matlab ® , Ver.2 with support for octonions[END_REF]. The quaternion Fourier transform of multiband images is a matrix of full quaternion values q, each of which satisfies S(q) ≠ 0 and V(q) ≠ 0 (see Eq. 3.3). The arrangement of coefficients within the quaternion spectrum follows the same layout as the one within the complex Fourier transform spectrum [START_REF] Ell | Hypercomplex Fourier transforms of color images[END_REF]. The quaternion Fourier spectrum, as in the complex Fourier transform, can be divided into four quadrants, and in each quadrant the coordinates correspond to exactly one of the four combinations of positive and negative horizontal and vertical spatial frequencies. The quaternion spectrum can be represented in the quaternion polar form using Euler's formula (see Eqs. 3.11,3.12). This representation yields three components: modulus, phase and axis. The modulus and the phase component are both scalar values.

(4.4)

⎧ ⎪ ⎨ ⎪ ⎩ 1 = (i + j + k)∕ √ 3, 2 = (j -k)∕ √ 2, 3 = 1 2 . (4.5) F(u, v) = 1 √ MN M-1 ∑ m=0 N-1 ∑ n=0 e -1 2 ( mu M + nv N ) ( c 1 (m, n) + c 2 (m, n) 2 ) = 1 √ MN ( M-1 ∑ m=0 N-1 ∑ n=0 e -1 2 ( mu M + nv N ) c 1 (m, n) + M-1 ∑ m=0 N-1 ∑ n=0 e -1 2 ( mu M + nv N ) c 2 (m, n) 2 ) = F 1 [u, v] + F 2 [u, v] 2 , (4.6) F i [u, v] = 1 √ MN M-1 ∑ m=0 N-1 ∑ n=0 e -1 2 ( mu m + nv N ) c i (m, n), i = 1, 2.
Their distributions can be visualized as a gray level image (see Fig. 6). The axis can be displayed as a color image (see Fig. 6).

To estimate the texture features, a periodogram analysis is performed using the quaternion Fourier transform of function f . Periodograms are calculated from the moduli of the Fourier transform of f. Periodograms analysis enables to identify the dominant cyclical behavior that might explain the variation pattern (color and texture) within the observed data. In this work we apply the zero-padding process to function f before taking the quaternion Fourier transform to deal with image border issues. The zero-padding process aims to obtain a square quaternion matrix from f by extending the matrix with lines and/or columns containing only zeros (Hamey [START_REF] Leonord | A functional approach to border handling in image processing[END_REF]. The polar representation of the DQFT spectrum of f allows to define the R-spectrum, phase spectrum, and axis spectrum. As in Eq. 3.11, it can be expressed as:

where R = |F DQFT [R, μ, θ]| is the DQFT modulus, μ the axis, θ the phase and (u, v) are coordinates in the frequency domain.

Let S be the number of lines of the image that we obtain thanks to the application of the zero-padding process to f . We characterize the scale texture pattern of f by means of the R-spectrum. We compute the mean values of the periodogram on concentric circles of radii r (with 0 ≤ r ≤ S 2 ) , as described in [START_REF] Couteron | Quantifying change in patterned semi-arid vegetation by Fourier analysis of digitized aerial photographs[END_REF]. Indeed, we characterize the color image texture pattern by performing a periodogram analysis of the DQFT spectrum as delineated by [START_REF] Mugglestone | A practical guide to the spectral analysis of spatial point processes[END_REF] for the classical Fourier transform. Periodograms corresponds to the number of Fourier spatial frequencies repetitions over a given distance, (generally the image's size). If this distance is equal to 100 m then it is expressed in cycles∕hectometer(cycles/hm) . Periodograms obtained from the DQFT modulus hold in a holistic manner all three data bands. R-spectra values are given by Eq. 4.8. b is the length between two consecutives concentric circles, and ∈ [0, 2 ] . In this paper we choose b = 1 when binning the radius in polar coordinates, and in each interval bin we sum the periodograms for all orientations on concentric rings. Indeed, the R-spectrum is a vector such that each component corresponds to the average value of the amplitudes of frequencies lying between the concentric circles of radius r and r + b with r ∈ {0, 1, 2, …} . The value of one component of the R-spectrum is given by Eq. 4.8.

(4.7) F DQFT [R, , ](u, v) = R(u, v)e (u,v) (u,v) , (4.8) k R (r) = 1 n r r+b ∑ r ′ =r 2 ∑ =0 | | | F DQFT [ r ′ , , ] | | | ,
Algorithm 1: computation of the DQFT from a 2D quaternion function [START_REF] Sangwine | Fourier transforms of colour images using quaternions or hypercomplex, numbers[END_REF])

Input: image f , 1 = (i + j + k)∕ √ 3, 2 = (j -k)∕ √ 2. Output: DQFT of f.

Begin

Step 1: Decomposition of image f(m, n) in symplectic form using equation 3.13:

f (m, n) = f 1 (m, n) + f 2 (m, n) 2 ,
Step 2: Expand symplectic components:

f 1 (m, n) = h 1,1 (m, n) + h 1,2 (m, n) 1 , f 2 (m, n) = h 2,1 (m, n) + h 2,2 (m, n) 1 ,
Step 3: Construction of equivalent complex images:

f ′ 1 (m, n) = h 1,1 (m, n) + h 1,2 (m, n)i, f ′ 2 (m, n) = h 2,1 (m, n) + h 2,2 (m, n)i,
Step 4: Perform 2D Fast Fourier transform of f 1 ′ (m, n) and f 2 ′ (m, n):

F ′ 1 [u, v] = H 1,1 (m, n) + H 1,2 (m, n)i, F ′ 2 [u, v] = H 2,1 (m, n) + H 2,2 (m, n)i,
Step 5: Construction of simplex and perplex parts of Fourier transform:

F 1 (u, v) = H 1,1 (m, n) + H 1,2 (m, n) 1 , F 2 (u, v) = H 2,1 (m, n) + H 2,2 (m, n) 1 ,
Step 6: Construction of symplectic form of the DQFT:

F[u, v] = F 1 (u, v) + F 2 [u, v] 2 ,
Step 7: Reconstruction of the DQFT in Cartesian form:

F[u, v] = G 0 (u, v) + G 1 (u, v)i + G 2 (u, v)j + G 3 (u, v)k. End

Support Vector Machine: Non-linear Learning Tool to Handle Color-Texture Indices

Nowadays, machine learning is becoming a major field of artificial intelligence. This domain aims to permit computers to acquire knowledge from the real world. Its final goal is to establish a learning computational theory in order to build systems that can perform in various application areas. The support vector machine proposed by [START_REF] Vapnik | Support vector method for function approximation regression estimation and signal processing[END_REF] is one of such methods widely used in artificial intelligence. It exists of two types of support vector machines: the support vector machine for classification (SVC) and the support vector machine for regression (SVR). SVC aims at classifying new unknown variables into one or more categories. SVR attempts to predict a real-valued output given a set of input variable data, as in our specific case. The SVR method consists in minimizing a quadratic objective function within a convex domain. The objective function is given by the combination of a loss function with regularization terms [START_REF] Smola | A tutorial on support vector regression[END_REF]. The loss function corresponds to the conventional least squares error criterion as used in traditional regression procedures. Thus, the SVR method attempts to minimize the generalized error bound to achieve a generalized performance rather than minimizing the observed training error. We begin this subsection with an overview of the linear SVR theory. We end it by taking into account the non-linearity properties of the phenomenon under study.

Principle of SVR

Let us consider the training sample set given by Eq. 4.9.

where

( x i , y i ) ∈ × ℝ, i = 1 … n;
is the input-space ( can for instance be equal to ℝ d , with d ∈ ℕ ), and y i , i = 1, … , n are the components of the predicted variable. We want to determine a function f from to ℝ which enables to predict the output values with a given precision ∈ ℝ with respect to training sample S and also as flat as possible. Then, for each input sample x i , we allow a deviation less than between the corresponding output sample y i and the predicted value f(x i ) as stated in Eq. 4.10.

Imposing the flatness of the function allows to minimize the complexity of the model and enhances its performances regarding generalization. In this subsection we describe how to get a linear function for function f and in the next subsection we will show how to cope with non-linearity using the kernel concept.

In the linear function case, the expression of f is given by Eq. 4.11.

where ⟨., .⟩ denotes the dot product in , x ∈ is the input vector, ∈ represents the vector of parameters to be determined, b ∈ ℝ is an unknown bias term to be found and y is the predicted variable. Hence, both ω and b need to be determined in Eq. 4.11 in order to get the expression of f. To obtain a maximal flatness of function f, we need to seek the minimal value of the parameter norm ‖ ‖ . Hence, we find f by resolving the convex optimization problem given by Eq. 4.12.

(4.9)

S = {( x 1 , y 1 ) , ( x 2 , y 2 ) , … , ( x n y n )} , (4.10) | | | y i -f ( x i ) | | | ≤ , i = 1, … , n. (4.11) f (x) = y = ⟨ , x⟩ + b, (4.12) { minimize 1 2 ‖ ‖ 2 ‖ ‖ subject to | | | y i -f ( x i ) | | | ≤ , i = 1, … , n.
In the SVR method, a minimal prediction error value is obtained through the constraints of the optimization problem. It supposes that the convex optimization problem is feasible, that is, there exist a function f that approximate the entire training sample with a precision , but sometimes this is not always the case in practice. Gen- erally, in the presence of important noise or outliers within the data, one can allow some errors using the concept of soft margin. The soft margin concept introduces two positives slack variables η i -and η i + in Eq. 4.12 in order to cope with infeasible constraints of the convex optimization problem. The error function in this case is called an -insensitive function, denoted by g and given in Eq. 4.13. It can be graphically interpreted as a tube of insensitivity with a radius as shown by Fig. 7. Hence, if the predicted value is within the tube then the error value is zero, while if the predicted value is outside the tube then the error is equal the magnitude of the difference between the predicted value and the radius of the tube. The slack variable η i + represents the distance between a point ( x i , y i ) and the bor- der above the tube and η i -the distance the distance between a point ( x i , y i ) and the border below the tube. It can be noted that for the same point, these values do not exist simultaneously. We rewrite Eq. 4.12 using the slack variables and obtain Eq. 4.14:

(4.13) g ((y -f (x)) = { 0 if |y -f (x)| ≤ , y -f (x) -if |y -f (x)| > .
The constant C > 0 defines the trade-off between the flatness of f and the amount up to which deviations larger than are tolerated. The Eq. 4.14 is easily solved when expressed in its dual formulation. The key idea is to construct a Lagrange function from both the objective function and the corresponding constraints. This is achieved by introducing a dual set of variables. An objective function is an equation to be optimized given certain constraints. [START_REF] Vapnik | Support vector method for function approximation regression estimation and signal processing[END_REF] proceeds as follows:

where the Lagrange multipliers + i , - i , + i and - i , i = 1, … , n, satisfy the non-nega- tivity constraints. Differentiating Eq. 4.12 with respect to b, , - i and η i + and setting the derivatives to zero for optimality yields: From Eq. (4.21), we can see that ω can be completely described as a linear combination of the training inputs x i , i = 1, … , n . Hence, the complexity of the SVR is independent of the dimensionality of the input space and depends only on the numbers of support vectors. That is one of the reasons why we choose -SVR for the prediction step of FOTO++.

(4.14) ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ minimize 1 2 ‖ ‖ 2 ‖ ‖ + C n ∑ i=1 ( - i + + i ) , subject to ⎧ ⎪ ⎨ ⎪ ⎩ y i -⟨ , x i ⟩ -b ≤ + + i , i = 1, … , n, ⟨ , x i ⟩ + b -y i ≤ + - i , i = 1, … , n, - i ≥ 0, i = 1, … , n, + i ≥ 0, i = 1, … , n. (4.15) L = 1 2 ‖ ‖ ‖ 2 ‖ ‖ ‖ + C n ∑ i=1 ( - i + + i ) - n ∑ i=1 - i ( + - i + y i -⟨ , x i ⟩ -b ) - n ∑ i=1 + i ( + + i -y i + ⟨ , x i ⟩ + b ) - n ∑ i=1 ( - i - i + + i + i ) , (4.16) L b = n ∑ i=1 ( - i -+ i ) = 0, (4.17) L = - n ∑ i=1 ( - i -+ i )x i = 0, (4.18) L - i = C - ( - i + - i ) = 0
The expression of b can be obtained by exploiting the so-called Karush-Kuhn-Tucker (KKT) conditions, which state that for the optimal solution the product between dual variables and constraints has to vanish [START_REF] Karush | Minima of functions of several variables with inequalities as side constraints[END_REF][START_REF] Kuhn | Nonlinear programming[END_REF]. In the support vector case we get:

There are three useful conclusions that can be made from the KKT conditions. First of all, only samples (x i , y i ) with corresponding α

i -= C or α i + = C lie outside (4.19) L + i = C - ( + i + + i ) = 0. (4.20) ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ maximize n ∑ i=1 y i ( + i -- i ) -1 2 n ∑ i,j=1 ( + i -- i ) ( + j -- j ) ⟨ x i , x j ⟩ - n ∑ i=1 ( + i + - i ) , subject to ⎧ ⎪ ⎨ ⎪ ⎩ n ∑ i=1 ( + i -- i ) = 0, + i , - i ∈ [0, C], i = 1, … , n. (4.21) = n ∑ i=1 ( + i -- i )x i . (4.22) y ′ = f ( x ′ ) = b + n ∑ i=1 ( + i -- i )⟨ x i , x ′ ⟩ . (4.23) ⎧ ⎪ ⎨ ⎪ ⎩ - i ( + - i + y i -⟨ , x i ⟩ -b ) = 0, + i ( + + i -y i + ⟨ , x i ⟩ + b ) = 0, - i - i = ( C -- i ) - i = 0, + i + i = ( C -+ i ) + i = 0, i = 1, … , n.
the -insensitive tube around f . Secondly, there can never exist a set of dual vari- ables - i , + i which are simultaneously non-zero, as this would require non-zero slacks in both directions. Finally, for α i -∈ (0, C) and + i ∈ (0, C) , we have - i = 0 and + i = 0 . Hence, b can be computed as given by Eq. 4.24. Its value is unique for any α i -, α i + such that 0 < - i and + i < C (Smola and Schölkopf 1998). The value of b is not unique only in the case where the optimality criteria are not hold that is when - i or + i ly on the boundary values 0 or C [START_REF] Keerthi | Improvements to Platt's SMO algorithm for SVM classifier design[END_REF]).

From Eq. (4.21) and the previous conclusions, it appears that ω does not need all the x i , i = 1, … , n in its expressions. The set of x i that come with non-zero coef- ficients is called the set of Support Vectors.

Non-linear SVR for FOTO++

The non-linearity of training patterns x i , i = 1, … n is taken into account in a non- linear support vector machine regression using a kernel [START_REF] Smola | A tutorial on support vector regression[END_REF]. The main idea consists in making a preprocessing of the inputs with a map Φ ∶ → into some feature space , and then to apply the standard linear sup- port vector regression procedures as described above. This approach appears to be directly not tractable from a computational point of view [START_REF] Smola | A tutorial on support vector regression[END_REF][START_REF] Basak | Support vector regression[END_REF]). One solution is obtained by finding a kernel function k (

x i , x ′ ) , i = 1, … , n that can replace the dot product in Eq. (4.22). We then get:

with

+ i ≥ 0, - i ≥ 0, i = 1, … , n, Φ ∶ → , k ( x i , x ′ ) = ⟨ ( x i ) , ( x ′ )⟩ , i = 1,

… , n, and

We also need to evaluate the value of b by replacing the dot product used in Eq. 4.24 by the kernel function. Our selection of the SVR model as a prediction tool is based on [START_REF] Tapamo | Linear versus non-linear methods: a comparative study for forest above ground biomass estimation from texture analysis of satellite image[END_REF], highlighting that it can provide promising results in biomass estimation and that it was the best one among other non-linear learning tools, such as random forest or k-nearest-neighbor. In this paper we used a non-linear SVR regression model with the popular radial kernel to estimate biomass values from color-texture features. The radial kernel expressions is given by Eq. (4.27).

(4.24)

{ b = y i -⟨ , x i ⟩ -for - i ∈ (0, C) or b = y i -⟨ , x i ⟩ + for + i ∈ (0, C). (4.25) = n ∑ i=1 ( + i -- i )Φ ( x i ) , (4.26) y ′ = b + n ∑ i=1 ( + i -- i ) k ( x i , x ′ ) .
where σ is a free parameter defining the kernel's behavior, and ‖⋅‖ is the Euclidian norm.

The FOTO++ Model

FOTO++ is a quaternion-based texture analysis model to address biomass estimation issues. FOTO++ characterizes multiband image textures directly in multispectral space. Our multispectral space includes: near infrared (NIR), red (R) and green (G). Indeed, we consider the texture of spatial distribution of pure quaternions without separating gray level textures (NIR, G, R + G + B, etc.) and color information. Therefore, our model merges in an elegant mathematical way classical texture properties and color information. We use the non-linear learning method (SVR) on multispectral space texture descriptors to infer biomass values. We notice that our methodology could be easily adapted to tackle other issues related to the analysis of multiband texture images. Algorithm 2 presents in details the texture image features computation, and Fig. 8 describes the corresponding flowchart. The biomass estimation scheme is given Algorithm 3, while Fig. 9 shows its flowchart.

(4.27) k ( x ′ , x * ) = e -‖ x ′ -x * 2 ‖ 2 2 ,

Output: Vector of color texture indices

Begin

Step 1: Filtering image I with combined Nagao and median filter,

Step 2: Converting filtered image into color space (RGB, HSV or Lab). Let C 1 C 2 C 3 be the resulting color image,

Step 3: Construct the equivalent quaternion image:

Q(x, y) = C 1 (x, y)i + C 2 (x, y)j + C 3 (x, y)k,
Step 4: Perform 2D Quaternion Fast Fourier Transform of Q(x, y) using Algorithm 1,

Step 5: Compute the modulus of the 2D Quaternion Fast Fourier Transform obtained in Step 4 and construct the R-spectra vector using equation 4.8,

Step 6: Save R-spectra vector as color texture indices.

End

Algorithm 3: Biomass estimation model

Input: A set M of multispectral satellite images, M = {M 1 , M 2 , …, M n }.
A set B of biomass values sampled on the fields,

B = { B 1 , B 2 , … , B n } .
There exists one biomass value for each image in M.

Output: non-linear model BiomassModel for biomass estimation.

Begin

Step 1: For each multispectral image M i i = 1, … , n, compute its color texture indices using Algorithm 2,

Step 2 Construct a table where each row corresponds to the color texture indices of one image,

Step 3: Compute the non-linear learning model BiomassModel using SVR with the R-spectra table values as predictive variables and B as predictors.

End

Fig. 10 Allostand output. The figure displays some Allostand output generated from field data. We obtain a simple 3D forest model built using allometric equation and physical properties of tree stands

Experimental Results

This section presents experimental results for biomass estimation provided with FOTO++ on real data. We have implemented FOTO++ model on Matlab ® version 2014a. We also use the Quaternion toolboxes release developed by Sangwine and Le Bihan [START_REF] Sangwine | Quaternion toolbox for Matlab ® , Ver.2 with support for octonions[END_REF].

We calibrate FOTO++ on synthetic data and validate it on high-resolution multispectral satellite images. A sensitivity analysis and a comparison with the classical FOTO method show the FOTO++ accuracy.

Data

Synthetic Multispectral Images

We produce one dataset of synthetic data using the Allostand and the DART frameworks. Allostand framework was developed by [START_REF] Barbier | Linking canopy images to forest structural parameters: potential of a modeling framework[END_REF], as a Matlab script to provide simulated forest stands from tree dendrometric data. Allostand is a forest structure model that yields simplified 3D forest simulations from a class of trunk Diameter at the Breast Height (DBH). Figure 10 presents an example of 3D forest images output from Allostand. DART software was developed in CES-BIO which aims at simulating the interactions between scene components and electromagnetic signals of varying wavelengths, in different sun-scene-sensor Dart frameworks, we have produced 1920 different multispectral images with 100 pixels per row and per column. Each image was simulated at high spatial resolution with three spectral bands that are well known to be strongly correlated with the forest characteristics: the near-infrared (NIR) band, the red band (R) and the green band (G). In order to include variability in our dataset, we have chosen different maximum DBH values. Among these dominant trees, we set among all trees on the field unit, a different proportion value of trees in the maximum DBH class. Maximum DBH varied from 50 to 100 cm by steps of 10 cm. The proportion of trees in the maximum DBH values was set to 0.33, 0.5, 1 and 2, as used in [START_REF] Barbier | Linking canopy images to forest structural parameters: potential of a modeling framework[END_REF] to study the bidirectional texture function of simulated forest canopy images. It represents the proportion of trees in the maximum DBH class among all the trees on the field unit. ( 1 ha in our experiments).We also produced synthetic images using different acquisition parameters. For the acquisition parameters, the view azimuth (VA) was set between 0° and 180° by steps of 20°, and the view zenith (VZ) between 5°and 40° by steps of 5°. Figure 11 shows some multiband images produced by DART, where we observe texture variability between images due to the variation of density maximum DBH values and acquisition parameter values.

Forest Surveys and Real Multispectral Satellite Data

Field stand structure measures of 12 plots located in Cameroonian forest were used to validate the FOTO++ method. Plots are located in the Cameroonian side of TRIDOM: three plots in the Kumu subdivision (district in Cameroon) within the South region, three plots in the Ngoila subdivision within the East region, and six plots in the Mindourou (Alpicam industry). The TRIDOM area is the tropical forest which lies in the Congo basin forest crossing Cameroon, Congo and Gabon. TRI-DOM is delimited by the Dja, Odzala, and Minkébé protected areas representing 141, 000 km 2 area. Figure 12 shows the geographical localization of these plots. On each plot, the DBH of each tree greater or equal to 10 cm was measured, and the Chave's allometric equation for wet forests was used to compute the aboveground biomass [START_REF] Chave | Improved allometric models to estimate the above ground biomass of tropical trees[END_REF]. Table 2 presents the biomass values of each plot.

We use three high spatial resolution IKONOS images covering these plots. The spectral bands processed by FOTO++ were Near InfraRed (NIR), Red (R) and Blue (B), and we tested FOTO on the panchromatic band as it was used in [START_REF] Couteron | Predicting tropical forest stand structure parameters from Fourier transform of very high-resolution remotely sensed canopy images[END_REF] and [START_REF] Ploton | Assessing aboveground tropical forest biomass using Google Earth canopy images[END_REF]. In most papers, FOTO appears to be processed with this band, mainly because of its high spatial resolution.

Practically, we needed a pre-processing in order to use the IKONOS image data. Below follows the description of the pre-processing step used to convert IKONOS digital numbers to surface reflectance values more suitable for analysis. We use the conversion method proposed by [START_REF] Thenkabail | Hyperion, IKONOS, ALI, and ETM+ sensors in the study of African rainforests[END_REF]. • Step 1 Firstly, the 11-bit IKONOS digital numbers provided by satellite sensors were converted to radiance values measured in watt per steradian per square meter (W sr -1 m -2 ) using Eq. 5.1 [START_REF] Dial | IKONOS satellite imagery, and products[END_REF]). • Step 2 Secondly, using Eq. 5.2, radiance values were converted to albedo [START_REF] Dial | IKONOS satellite imagery, and products[END_REF], a dimensionless number, and • Step 3 Finally, these values were corrected, from atmospheric effects to provide surface reflectance values.

where the i, k, 1 ≤ i ≤ MN , k = 1, 2, 3 , represent the band indices; M and N are the numbers of lines and columns respectively; L i,k and DN i,k are respectively the (5.1) in-band radiance at sensor aperture and the image product digital value of the ith pixel in the kth band; CalCoef k is the in-band radiance calibration coefficient for the band indicated by k. We have that k = 1 corresponds to band 1 (Green), k = 2 that corresponds to band 2 (Red) and k = 3 that corresponds to band 3 (NIR). The values of CalCoef k were 727 for band 1, 949 for band 2 and 843 for band 3 [START_REF] Dial | IKONOS satellite imagery, and products[END_REF].

L i,k = DN i,k × CalCoef -1 k Table 2 Aboveground biomass (AGB)
with k = 1, 2, 3 the band indices; k the at-satellite albedo of band k ; L k the radiance of band k; d the earth to sun distance in astronomic unit at the acquisition date; Sol λk is the solar flux of band k, and θ s is the solar zenith angle.

Figures 13, 14 and 15 present satellite images of each plot in the Kumu, Ngoila and Alpicam field sites respectively, after pre-processing.

(5.2)

k = L k d 2 × ( Sol k cos ( S )) -1 ,

Model-Performance Statistics

In order to precisely evaluate the contribution of the FOTO++ model in aboveground biomass estimation, we calculated and compared the performances of the classical FOTO model, FOTO with SVR used as prediction tool model [START_REF] Tapamo | Linear versus non-linear methods: a comparative study for forest above ground biomass estimation from texture analysis of satellite image[END_REF]) and the FOTO++ model. The same parameters of SVR was used both for the synthetic multiband images and the satellite multiband images. We use = 0.1 , C = 1 and the radial function as the kernel.

We investigated two model-performance statistics commonly found in literature. The first one is the root mean square error (RMSE) whose formula is given by Eq. 5.3: wherein model-derived biomass estimates (BiomassPred i , i = 1, 2, …, n) are compared with the corresponding biomass field values (BiomassField, i = 1, 2, … , n) that are judged to be reliable.

RMSE is frequently used to measure the differences between values predicted by a model and the values actually observed. When evaluating model performances, the smallest RMSE correspond to the best models. Thus, if the ratio between the RMSE of two models is less than 1, then the model with the RMSE in the numerator is the best. Otherwise, the model with the RMSE in the denominator is the best. Subtracting this ratio from 1 yields the amount (in percentage) of improvement gained by the best model. Greater errors have a greater influence on the RMSE criterion when comparing models. We also used jointly another model performance index in order to get a more complete performance evaluation.

(5.3) RMSE = √ √ √ √ 1 N N ∑ i=1 ( BiomassPred i -BiomassField i ) 2 ,
The second performance criterion we used in this work is the dimensionless refined index model d r proposed by [START_REF] Willmott | A refined index of model performance[END_REF]. It lies in the range [-1, 1]. In order to compute this index, we calculated the mean absolute error (MAE) as given by Eq. 5.4 and the mean absolute deviation (MAD) given by Eq. 5.5. Willmott's index formula is given by Eq. 5.6. The interpretation of this index as stated by Willmott et al. (2015) depends on the sign of the index. When the index value d r lies in the positive domain, it describes the proportion to which a set of model predictions is, on average, error free. In the other case, when the index value lies in the negative domain, it describes the proportions of the real data error that is underrepresented by the average-error magnitude of the model.

(5.4) that are judged to be reliable, and BiomassField represents the mean of biomass field values. The MAD value of our dataset is equal to 93.8033 kg/ha.

MAE = 1 N N ∑ i=1 | | BiomassPred i -BiomassField i | | , (5.5) MAD = 1 N N ∑ i=1 | | | BiomassField i -BiomassField | | | ,

Sensitivity of Acquisition Parameters Analysis

We have conducted a sensitivity analysis of the FOTO++ method according to acquisition parameters and canopy grains in order to study the influence of instrumental bias. We have computed the texture indices of each image of the simulated dataset by means of the FOTO++ model. The corresponding texture indices of the classical FOTO model were computed on the NIR spectral band of each image. We have compared R-spectra computed by FOTO++ and FOTO within their prediction steps. It appears that color-texture indices gave more information on images. Texture features provided by FOTO++ and FOTO are shown in Fig. 16 for several acquisition parameters. We have observed that FOTO++ indices hold more variability (both color and texture) than FOTO indices. It seemed that using simultaneously color and texture information to describe image contents can be a valuable tool. This fact is confirmed when evaluating the RMSE produced by these models to estimate aboveground biomass. This allowed us to state that quaternionic R-spectra used as texture indices can be a valuable tool for texture analysis of multiband satellite data.

We evaluated the accuracy of FOTO++ regarding biomass estimation versus the FOTO improved method proposed by [START_REF] Tapamo | Linear versus non-linear methods: a comparative study for forest above ground biomass estimation from texture analysis of satellite image[END_REF], called here FOTO + SVR. We trained our model on 1344 synthetic multiband images and assessed its performances on the remaining 576 synthetic multiband images. We report in Table 3 only the best accurate model based on the RMSE with respect to parameter acquisitions.

We observed that the FOTO++ method performed better in 56.67% of the cases, and when it is not the case, its RMSE values are very close to those of FOTO + SVR, as shown on Figs. 17 and 18. These figures present the RMSE yielded for view azimuth angles respectively of 15° and 40°, and show the accuracy variations with respect to acquisition parameters. The fact that FOTO++ is sometimes equivalent to FOTO + SVR in terms of estimation of biomass was due to the fact that in some generated images, nearly all the structural and the color information were concentrated only in the NIR band. Colorless information in this experimental dataset did not allow us to point out the most suitable color space to encode color information. We will see in the next experiments on colorful real datasets that color can contribute significantly to the accuracy of biomass estimation methods.

(5.6) We compare the accuracy of each method (FOTO, FOTO + SVR [START_REF] Tapamo | Linear versus non-linear methods: a comparative study for forest above ground biomass estimation from texture analysis of satellite image[END_REF]) and FOTO++). We report below the most accurate method. In parenthesis we have the RMSE value of the most precise method. In bold we have the RMSE of the FOTO++ method indexed by the color spaces used. SVR means here the FOTO + SVR method. HSV, Lab and RGB mean the FOTO++ method processed respectively with the HSV, the

d r = { 1 -MAE 2MAD if MAE < 2MAD, 2MAD MAE -1 if MAE > 2MAD.

Biomass Model Comparison: FOTO versus FOTO++

In this subsection, we assess the performance of our model FOTO++ and also evaluate the contribution of each main steps of the model. This evaluation is conducted on 'real' multiband satellite images. We trained the model on 4 satellite images (1 for the Kumu dataset, 1 for the Ngoila dataset and 2 for the Alpicam dataset) and investigated its performances on the remaining 8 satellite images. We present below the results obtained and their discussion.

Contribution of Support Vector Machine for Regression in Biomass Estimate Accuracy

In this section we compare FOTO with FOTO + SVR as a prediction tool in order to assess the statement that using non-linear machine learning methods can improve the performance of FOTO, as stated in [START_REF] Tapamo | Linear versus non-linear methods: a comparative study for forest above ground biomass estimation from texture analysis of satellite image[END_REF]. (2014), and confirmed that machine learning tools can sensibly improve FOTO biomass estimates. We needed also to evaluate which part of this improvement is due to the color-textures features extraction process using quaternions tools. The next subsection presents the evaluation of this contribution.

Contribution of Color-Texture in Biomass Estimate Accuracy

This subsection aims to evaluate the contribution of image indices derived from multiband images processed in a holistic way using the quaternion algebra tool. We used color-texture indices extracted from multiband satellite data by means of the discrete quaternion Fourier transform (Algorithm 1). We calculated and compared RMSE values and Willmott's model performance index of FOTO-SVR and FOTO++ without taking into account the filtering step. These experiments were conducted on a real satellite dataset. Table 5 summarizes these experiments.

The ratio of RMSE values for the best FOTO++ without filtering and the one for FOTO + SVR is 0.7598. Thus, FOTO++ without filtering improved by 24.02% the accuracy of FOTO + SVR. Willmott's performance index also confirmed the improvement of the precision using FOTO++ while getting 0.8170 instead of 0.785 for FOTO + SVR. We can then conclude, based on these results, that color-texture indices extracted from multiband satellite data significantly improve biomass estimates. As the FOTO++ complete method used an image filter in its first step, an assessment of its contribution to the final improvement of the method is required. The next subsection presents this evaluation.

Contribution of Filtering in Biomass Estimate Accuracy

In this last subsection, we evaluated how a filtering step using the Nagao-median filter contributes to biomass estimate accuracy in the FOTO++ model. We also assessed this evaluation according to color representations. Table 6 summarizes the results yielded by this scheme. We observed, based on RMSE, that a filtering step can be used to improve biomass estimate accuracy in both color representations. For our data, RMSE values of FOTO++ using the initial filtering step vary from 43.366 kg for the Lab color space, to 42,331 kg for the HSV color space. We calculated the amount of improvement obtained while filtering data in each color space and we found 6.36% for the RGB color space, 5.27% for the HSV color space and 5.42% for the Lab color space. Contrariwise, Willmott's index showed that the Lab color space yielded a smallest performance model. It seems that a filtering step enhanced both color and texture information for our dataset. Our experiments showed that a filtering step does not change significantly the Euclidian distance of colors on both non-filtered and filtered data, but changes seemed to occur in perception of colors. According to this index, the HSV color space representation yielded the most adapted model for biomass estimates. We observed here that filtering data contributes effectively to the improvement the final performance of the model. We can conclude, based on these experiments, that color-texture indices extracted with the quaternion Fourier transform combined with a machine learning tools and with the use of filtering techniques can significantly improve the accuracy of biomass estimates compared with the classical FOTO model.

Conclusion

Satellite image data appear as a primary source for measuring past trends and monitoring changes in forest carbon stocks. Thus, the processing of these data represents a fundamental challenge. Image processing techniques applied on remote sensing data are valuable tools to fix this issue. The problem of estimating biomass from satellite data can be viewed as image texture analysis and classification. One such approach is the FOurier-based Texture Ordination (FOTO) method proposed by [START_REF] Couteron | Quantifying change in patterned semi-arid vegetation by Fourier analysis of digitized aerial photographs[END_REF] for aboveground forest biomass estimation. It ordered canopy images according to coarseness-finesses texture gradient and used a simple regression model to estimate aboveground biomass values. The limitation of this method lies on its inefficiency to process multiband satellite images. We note that classical image processing algorithms are not able to simultaneously deal with multiband satellite images. Instead, when dealing with these data, previous methods first separate it into several channels and then apply traditional image processing algorithms to these channels. In recent years, some approaches using quaternion numbers to handle and model multiband images in a holistic manner were introduced. By introducing quaternions, quaternion Fourier transform allows multidimensional data such as multiband images to be processed efficiently and easily. In this paper we propose the FOTO++ model that takes into account in a compact elegant mathematical way all information hold in multiband satellite images. We used this information to address aboveground biomass issues in tropical regions. FOTO++ starts by removing noises in multiband images while preserving the edges of objects. After this preprocessing stage, multiband texture indices were extracted using the discrete form of the quaternion Fourier transform, and finally the support vector regression method was used to derive biomass estimations from these multiband texture indices. We used quaternion-based tools to extract useful texture features from multiband satellite images. Texture features are represented by a vector composed by the radial spectrum obtained from the amplitude of the quaternion Fourier transform. The sensitivity analysis of the FOTO++ model yielded that it is often more efficient with regards to the variation of acquisition parameters than the FOTO model. A comparative study between the two models based on model-performance statistics showed that FOTO++ is more accurate than FOTO + SVR. Color representation, the filtering step as well as machine learning seemed to be valuable tools to significantly improve biomass estimates. We observed indeed that, over all, FOTO++ provided 54.49% of improvement compared with FOTO. Also, 82.11% of biomass values yielded was on average free of error (according to Willmott's index), and biomass values obtained using by FOTO + SVR was on average 59.5% free of error. We investigated in detail the contribution of each step of FOTO++ to this final improvement and found that when using SVR, FOTO accuracy is improved by 36.8% with on average 78.5% of values free of error. When applying SVR on quaternion-based color-texture, an improvement of 24.02% is observed with on average 82.70% of biomass value free of error. The filtering step of FOTO++ contributed to 5.27% of the performance and yielded on average 83.11% of biomass values free of error.

Hence, quaternion-based color-textured and filtering are valuable tools that can be used to improve biomass value estimates.

We also observed that whatever is the color space used, FOTO++ has more often a better performance than both FOTO and FOTO + SVR used as prediction tool. These results highlight the potential of the quaternion-based approach to classify multiband texture. We validated this approach on forest biomass estimation, both on simulated and real image data. Our methodological scheme can be further improved while handling other information derived from the DQFT, such as phase spectrum and the axis spectrum. A similar methodology could be investigated with quaternion wavelets and is currently being tackled to address multiband satellite texture classification issues. Also, other hypercomplex numbers such as octonions can be used to handle up to eight spectral bands, but from our knowledge the theoretical part of this algebra is not well advanced to derive such methods. On the other hand, Clifford algebra can be used to process any number of spectral bands. Clifford algebra represents an opportunity to handle hyperspectral images with similar methodology as presented in this paper.

Fig. 2

 2 Fig. 2 Color spaces. The left figure corresponds to the HSV color space and the right figure presents the Lab color space

Fig. 3

 3 Fig. 3 IKONOS images of canopies in the Korup forest (South-west region of Cameroon). The spatial resolution is 4 m. Each image has a dimension of 200 × 200 pixels which corresponds to 16 ha in the field. a RGB color representation of Korup. The color image provided here is a pseudo-color image with the near infrared band coded in the red channel, the red band coded in the green channel and the green band coded in the blue channel. b Image in HSV color space viewed as a RGB color image. c The hue component. d The saturation component, and e the value component. (Color figure online)

Fig. 4

 4 Fig. 4 Nagao neighborhoods. The black square corresponds to the pixel of interest (central pixel) of the neighborhood. The shaded squares represent the neighborhoods of the central pixel

Fig. 6

 6 Fig. 6 DQFT of canopy forest image. a Original image, b modulus part, c axis part and d phase part of DQFT

Fig. 7

 7 Fig.7Regression with -insensitive tube. The solid line represents the predicted function f(x), and the dotted lines the boundary of the tube of radius around the predicted function The y-axis corresponds to the dependent variable y, and the x-axis to the independent variable x, for the case where = ℝ . Any point lying inside this tube has a loss value of zero, while all other points, either above or under the tube boundary, are required to have a loss value of η + or -respectively. All points lying exactly on the boundaries are called the support vector

Fig. 8

 8 Fig. 8 FOTO++: color-texture indices extraction

Fig. 11

 11 Fig. 11 Multispectral images simulated using Allostand and Dart frameworks. Each row corresponds to the density factor of the largest DBH (here DBH max = 50 cm), and each column corresponds to various acquisition parameters (VA View Azimuth, VZ View Zenith) of the same simulated stand

Fig. 12

 12 Fig. 12 Sample of field where dendrometric values were collected. They are localized in the South and East region of Cameroon

  values of each field site Chave's allometric equation for tropical wet forest[START_REF] Chave | Improved allometric models to estimate the above ground biomass of tropical trees[END_REF]) was used to estimate the biomass from data collected on the

Fig. 13

 13 Fig. 13 IKONOS satellite images of field samples in the Kumu subdivision

Fig. 15

 15 Fig. 15 IKONOS satellite images of field samples in the Alpicam forest industry

Fig. 16

 16 Fig. 16Comparison of R-spectra using different acquisition parameters: x-axis represents R-spectra bins and y-axis represents the amplitude of the R-spectra. RGB's graph, HSV's graph and Lab's graph are sketched using FOTO++ method where images are represented in RGB, HSV and Lab color spaces, respectively. Periodograms corresponds to the number of Fourier spatial frequencies repetitions over a given distance, (generally the image's size). In our experiment, this distance is equal to 100 m and then it is expressed in cycles/hectometer (cycles/hm)

Fig. 17

 17 Fig. 17 Comparison of RMSE using FOTO + SVR (green line) and FOTO++ (blue dots) on different azimuthal view angles for view zenith angle equal to 15°: x-axis represents azimuthal angles and the y-axis represents the RMSE. (Color figure online)

•

  Biomass estimation model Texture indices are used as predictors in a simple regression model where biomass is a dependent variable. The biomass values of the field are used during this step to build an estimation model.

  ,

	Substituting Eqs. (4.16), (4.17), (4.18), (4.19) into Eq. (4.15) yields the dual optimi-zation problem given by Eq. (4.20) (Xu et al. 2014):
	Equation (4.17) can be rewritten as follows:

By substituting Eq. (4.21) in Eq. (4.11), we get for a new input x′ a new prediction y′ as follows:

Table 3

 3 Root of mean square error (RMSE)

Table 4

 4 Non-linear statistical learning tools and biomass estimate accuracyWe evaluated RMSE and Willmott's index d r for FOTO and FOTO + SVR. We tested these models on a real satellite dataset

		FOTO	FOTO + SVR
	RMSE (kg)	92.998	58.826
	MAE (kg)	75.851	40.175
	d r	0.595	0.785

  Table 4 presents the RMSE and Willmott's model-performance indices. We observed that FOTO yields a RMSE of 92,998 kg, and that when SVR is used as prediction tool, its RMSE decreases down to 58,826 kg. The ratio RMSE FOTO+SVR /RMSE FOTO is then equal to 0.632. Hence, FOTO + SVR model is better than the FOTO model, and the improvement is about 36.8%. Willmott's index of FOTO + SVR, getting a value of 78.5%, also showed that biomass values obtained by FOTO + SVR are on average 78.5% free of error. Biomass values given by FOTO are only on average 59.5% free of error. Hence SVR applied on FOTO texture indices yielded more accurate biomass estimates. These results were in accordance with those found in Tapamo et al.

Table 5

 5

	Contribution of color-texture indices to biomass estimate accuracy	FOTO + SVR	FOTO++ without filtering step RGB HSV Lab
	RMSE (kg)	58.826	45.522	44.689	45.963
	MAE (kg)	40.175	34.642	34.319	35.276
	d r	0.785	0.8153	0.8170	0.8119
	We evaluate the RMSE and the Willmott's index d r of the FOTO + SVR and FOTO++ models. We tested these models on a real satellite dataset
	Table 6 Evaluation of filtering step contribution to biomass estimate accuracy	FOTO++ applied on non-filtered data RGB HSV Lab	FOTO++ applied on filtered data RGB HSV Lab
	RMSE (kg) 45.522 44.689 45.852 42.624 42.331 43.366
	MAE (kg) 34.642 34.319 35.276 31.756 31.675 33.997
	d r	0.8153 0.8170 0.8119 0.8307 0.8311 0.8187
	We have evaluated the RMSE and the Willmott's d r indices of FOTO++ using non-filtered and filtered image data. We have con-ducted these experiments on a real satellite image dataset

Acknowledgements

The authors would like to thank Pierre Ploton for his important contribution to data collection and for fruitful discussions on tropical forest structures. Many thanks go also to Nicolas Barbier, Nicolas Picard and Bonaventure Sonke for their expertise about tropical forest trees, dendrometric properties, forest dynamic processes and radiometric properties of trees, and to the anonymous reviewers for their careful reading of our manuscript and their many insightful comments and suggestions. Finally, they thank all the members of the SAM's team of UMMISCO-Cameroon. The BIOFORAC Project and the PDI (Programme Doctoral International) of IRD and Paris 6 University have funded this research.

Appendix 1

We present here a complete scheme to compute the polar form of any quaternion numbers.

In Sect. 3.3, when dealing with quaternion algebra, we have shown how to compute the polar representation of a quaternion only in the case where it is non-equal to zero. We present here how to compute the polar form when either the scalar or the modulus of its vector part is equal to zero. Let with |q| ≠ 0 . The polar representation of the quaternion q can be expressed as:

where α is a pure quaternion with |V( )| ≠ 0 . It follows that:

ii. If cos(|V(q)|) < 0 and sin (|V(q)|) > 0 with arctan ( |V(q)|

S(q)

) ∈

( 2 + 2k , (2k + 1)

iii. if cos (|V(q)|) < 0 and sin(|V(q)|) < 0

Appendix 2a

We describe here below an algorithm to convert a RGB-color image to Lab color space.

if S(q) = q 0 ≠ 0

Appendix 2b

We present here an algorithm to convert a RGB-color image to HSV color space.