Mahyar Shahsavari
email: mahyar.shahsavari@lifl.fr

Faisal Nadeem

Arash Ostadzadeh

Philippe Devienne

Pierre Boulet

M Faisal Nadeem

S Arash Ostadzadeh

Unconventional digital computing approach: memristive nanodevice platform

Keywords: material implication, nanoscale crossbar, memristor nanodevice, digital computing

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1 Introduction Memristor has recently drawn wide attention of scientists and researchers due to non-volatility, better alignment, and excellent scalability properties [START_REF] Helen | Compact model of memristors and its application in computing systems[END_REF]. Memristor remembers its last state after the last power plugging and has simple physical structure, high-density integration, and low-power consumption. These features make memristor an attractive candidate for building the next generation non-volatile memories [START_REF] Ho | Nonvolatile memristor memory: device characteristics and design implications[END_REF]. From highperformance computing point of view, memristor has potential capability to conquer the memory bottleneck issue, by utilizing computational unit next to the memory [START_REF] Ventra | [END_REF].

The logic computing applications of memristor have been investigated by several researchers [START_REF] Raja | Digital logic implementation in memristor-based crossbars -a tutorial[END_REF][START_REF] Borghetti | [END_REF][START_REF] Lehtonen | Implication logic synthesis methods for memristors[END_REF][START_REF] Rose | Proceedings of the IEEE[END_REF][START_REF] Vourkas | [END_REF][START_REF] Kvatinsky | Memristor-based imply logic design procedure[END_REF][START_REF] Lehtonen | [END_REF][START_REF] Lehtonen | Stateful implication logic with memristors[END_REF][START_REF] Borghetti | [END_REF][13]. For instance, Borghetti et al. [START_REF] Borghetti | [END_REF] used material implication (IMP) logic operation to carry out logic computation by using memristors. In IMP approach, A IMP B (A→B) operation means 'if A, then B' and can be read as A implies B. IMP together with a FALSE gate are able to form a functionally complete set (any boolean function can be expressed). By applying memristor as a digital switch, a high memristance (memristor resistance) is considered as logic '0' and a low memristance is considered as logic '1'. An-other approach to make gates by a two-terminal device as a switch is the programmable crossbar architecture [START_REF] Vourkas | [END_REF]13]. The crossbar nanowire array architecture can be used to compute logic functions by using memristor as a switch between two nanowires [START_REF] Raja | Digital logic implementation in memristor-based crossbars -a tutorial[END_REF].

The main contributions of this paper are:

-To demonstrate that memristor has appropriate characteristics to be applied as a switch in both IMP and crossbar array methods, we proposed our electrical model of the memristor behavior. Modeling is performed by using Verilog-A Hardware Description Language (HDL). -Material implication (IMP) logic function and crossbar array architecture are investigated as the two novel and most promising methods for creating logical functions by memristors.

-The crossbar array architecture approach is developed in the term of using only memristors as switches between nanowires.

The reminder of the paper is organized as follows. The electrical model is presented in Section 2. Then, we apply our Verilog-A model in the HSPICE circuit simulator to observe memristor behavior in the circuits. In Section 3, the stateful logic operations via IMP is described. In Section 4, the crossbar array architecture is presented by using memristor as a junction switch. Evaluations of the two approaches are investigated in Section 5. Section 6 concludes the paper.

2 Electrical model When an electric field is applied to the terminals of memristor, the shifting in boundary between its doped and undoped regions leads to variable total resistance of the device. In Figure 1.a, the electrical behavior of memristor can be modeled as follows [14]:

v(t) = R mem i(t) (1)
R mem = R ON w(t) D + R OF F (1 - w(t) D ) (2)
where w(t) is the width of the doped region, D is the overall thickness of the TiO 2 bi-layer, R ON is the resistance when the active region is completely doped (w = D) and R OF F is the resistance, when the TiO 2 bi-layer is mostly undopped (w→ 0).

dw(t) dt = µ v R ON D i(t) (3)
which yields the following formula for w(t):

w(t) = µ v R ON D q(t) (4)
Where µ v is the average dopant mobility. By inserting Equation (4) into Equation ( 2) and then into Equation [START_REF] Helen | Compact model of memristors and its application in computing systems[END_REF] we obtain the memristance of device, which for R ON R OF F simplifies to:

M (q) = R OF F (1 - µ v R ON D 2 q(t)) (5)
Equation [START_REF] Borghetti | [END_REF] shows the dopant drift mobility µ v and semiconductor film thicknesses D values, which are two factors with crucial contributions to the memristance magnitude. Subsequently, we can write Kirchoff's voltage law for memristor given by:

v(t) = M (q)i(t) (6)
By using Verilog-A HDL, we simulate the behavior of memristor, based on its behavioral equations. To investigate the characteristics of memristor in electrical circuits, the Verilog-A model of memristor behavior must be applied as a circuit element in the HSPICE netlist. In the HSPICE circuit, we apply a sinusoidal source to observe the memristor reaction in a simple circuit consisting of memristor and sinusoidal source. plot of memristor terminals that we measured in our simulation. This i -v plot, which is the most significant feature of memristor [15], is namely called "pinched hysteresis loop". The i -v characteristic demonstrates that memristor can [U+201C]remember[U+201D] the last electric charge flowing across it by changing its memristance. Therefore, we can use memristor as a latch to save data and also as a switch for computing. Moreover, in Figure 1.b, it is depicted that the pinched hysteresis loop is shrunk by increasing frequency. In fact, when the frequency increases toward infinity, memristor behaves similar to a linear resistor.

3 Stateful implication logic One of the basic potential applications of memristors is to utilize them in building blocks of logic gates. Therefore, by applying a digital pulse voltage to the memristor terminal, we have a switch with ON or OFF state. Unlike conventional CMOS, in memristor-based gates, data will be stored as a resistance rather than a voltage. In this case, the latches are nonvolatile. Thus, R ON displays logic '1' which means closed switch and R OF F displays '0' for presenting open switch. In contrast to the three-terminal CMOS-based transistor as a switch, in a two-terminal switch, there is no terminal to control ON or OFF states of the switch. Consequently, instead of conventional boolean logic, we should find other substitutes to create a gate and perform computing.

IMP (Figure 2.a) is a way to use one memristor to control the other one. IMP is recognized as a promising method for making gates by memristors [START_REF] Kvatinsky | Memristor-based imply logic design procedure[END_REF][START_REF] Lehtonen | [END_REF][START_REF] Lehtonen | Stateful implication logic with memristors[END_REF][START_REF] Borghetti | [END_REF]. In IMP structure, memristors have different roles in different stages of the computing process: input, output, computational logic element, and a latch depending upon which write, read, computing and storing processes are taking place, respectively.

To figure out how IMP operates, imagine A as a question and B as the answer to that question. If the question is wrong, any answer (wrong or correct) makes a true output (logic 1), as depicted in Figure 2.b. The only case for the false (logic 0) output would be a wrong answer to a correct question. So the implication logic is equivalent to function (¬A)∨B. Figure 2.a shows the basic circuit of memristors A and B to perform implication logic, which are formed by the vertical nanowire crossing over the horizontal nanowire connected to a load resistor R G . After the operation of material implications, the result is stored as the state of switch B (B ) while A is left unchanged. To switch between logic 1 to logic 0 (and vice versa), we need a tri-state voltage driver with a high impedance output state when it is undriven. V set is a negative voltage which should be applied to its corresponding tri-state driver. V set can switch memristor to conductive state with low resistance R ON . Similarly, the positive voltage V clear is required to change the memristive switch state to lowconductance (high-resistance) state R OF F . It is important to mention that the magnitude of V set and V clear must be larger than device threshold voltage for switching ON and OFF. In order to remain in a specific state (line 3 in truth table 2.b), the V cond is applied as a negative voltage with a magnitude smaller than V set . Consequently, tri-state drive -since is not in high-impedance state-is pulsed by one of the V set , V clear or V cond . By applying V cond and V set to A and B simultaneously, the memristive IMP operates properly. Although the conditional voltage (V cond ) is not necessary, except for the case AB='10' (third line in the truth table Figure 2.b), it is possible to either apply V cond or use high-impedance (HZ) for all other cases. If V set is applied to B alone, it would be unconditionally 'logic 1', nevertheless applying V cond by itself to A does not change its state. On the other hand, if both voltages are applied together, the present state of switch A determines the next state of switch B. If A='0', it means memristor A is in high resistance state (R OF F ). Therefore, there is a small voltage drop across R G . In this case, B will be set and A is left unchanged. Alternatively, if present state of A='1', switch A is in low resistance state and V cond drops across R G , so both A and B remain unchanged. It should be noted that the R G value must be chosen such that R ON <R G <R OF F , where R ON and R OF F are resistance states of 'ON' and 'OFF' switches, respectively. In boolean logic algebra, a set of logic operations is called functionally complete if all other logic functions can be constructed from combining the members of this set. The single-element sets NAND (AND, NOT) and NOR (OR, NOT) are functionally complete operations. In addition, it has been demonstrated that all boolean expressions can then be written in one of the standard normal forms using only a (IMP) operation and a false (Inverse) operation [START_REF] Lehtonen | Stateful implication logic with memristors[END_REF][START_REF] Borghetti | [END_REF]. In fact, to show that the memristive implication logic is a functionally complete operation; the easiest way is to synthesize the NAND function with it.

A circuit with three memristors A, B, and C is illustrated schematically in Figure 3.a. Assume two implication operations being performed subsequently, first A IMP C and then B IMP C. Hence A and B are inputs and C is output. The final output result is represented with variable C in Figure 3.c. Firstly, V clear should be applied to switch C to create the false operation (C = (C IMP 0)). By applying V cond and V set pulses to A and C respectively, the second step would be performed. Finally, C = (B IMP C) is yielded by applying V cond to V B and V set to V C (see Figure 3.b and Figure 3.c). In other words, the resulting state of switch C can be written as:

C = B IMP (A IMP C) = ¬B ∨ ((¬A) ∨ C) if C = 0 initially then C = ¬B ∨ ((¬A) ∨ 0) = (¬B) ∨ (¬A) = ¬(A ∧ B)
which is a NAND operation. Similarly, we can produce equal IMP structure of other logic operations. In Table 1 different boolean logic gates are listed. Obviously, we can demonstrate all logic relations in Table 1, not only by applying boolean logic rules but also by checking the truth tables of both sides of relations.

A A IMP 0 3 A AND B {A IMP (B IMP 0)} IMP 0 4 A NAND B A IMP (B IMP 0) 3 A OR B (A IMP 0) IMP B 3 A NOR B {(A IMP 0) IMP B} IMP 0 4 A XOR B (A IMP B) IMP {(B IMP A) IMP 0} 3
4 Crossbar architecture Programmable crossbar architectures have been proposed as a promising approach for future computing architectures because of their simplicity of fabrication and high density, which support defect tolerance [13,16]. In such architectures, assume that each junction within the crossbar can be utterly configured to activate an electronic device, such as a resistor, diode, transistor, or recently memristor. In fact, attractive features of memristor, such as simple physical structure, non-volatility, high-density, and unlimited endurance make this nano-device one of the best choices to play a switch role in the crossbar junctions. The memristive-based crossbar opens new windows to explore advanced computer architecture, different from the classical Von Neumann architecture [16]. On the other hand, in crossbar architecture, memory and logic operators are not separated. The memory can perform logic implementations on the same devices which store data. This is because during the operation process, control signals determine which elements act as logic gates and which ones act as memory cells.

Memristive switches in crossbar architectures

Each junction in a crossbar could be connected or disconnected by replacing a memristor as a switch in the junction point between two vertical and horizontal nanowires as depicted in Figure 4.a. Such a switch is in high resistance (R OF F ) (Figure 4.c) for open state or low resistance (R ON ) (Figure 4.d) for close state, similar to the states of memristor in Section 3. The memristive switch retains its state as long as the voltage drop across it, is not more than the required threshold voltage to change the memristor state. In Figure 4.a the input voltage (v in ) either could be connected to the output voltage of an external CMOS circuitry or be connected to the output of another latch from the previous stage. As we have already mentioned, data in our architecture is saved as the resistance of the memristive latch. However, for data transferring, an input voltage would be necessary. This input voltage is driven by the state of the input impedance (R in in Figure 5), which is a memristive switch. This voltage can be disconnected which means the floating state (for instance, a very high impedance R OF F ). The floating state input happens when the input memristor (R in in Figure 5.a) is OFF. Subsequently, the input voltage can also be a fixed negative voltage (-V f ) with less magnitude than the threshold, while the input memristor (R in in Figure 5.b) is ON.

The memristive switches can be configured as an inverting or non-inverting mode as it is depicted in Figure 5.c and Figure 5.d, respectively. If the stored data in the receiving switch is the logical complement of the input data, the configuration is in the inverting mode. If the stored data is the same as the input data the configuration is in the non-inverting mode. Different control signals make different configurations. Inverting configuration requires three steps to perform the appropriate latch operation: 1) We preset the switch unconditionally open (Figure 5.a) by applying a positive voltage (more than the positive threshold) to the control of the vertical nanowire and also by forcing input to the high impedance mode. A diode is also required to provide a low-impedance path to the ground to protect the junction.

2) If the input voltage is logic '1' then R in =R ON and -V f drops across R G , so the voltage across memristor is not enough to close it and the switch remains OFF. On the other hand, when input is logic '0', then R in =R OF F and memristor switch turns ON. Thus, the junction has held the inverted state of the input (see Figure 5.b).

3) The input signal must be in high impedance (disconnected). The state of the switch is read out onto the horizontal nanowire. This is accomplished by driving the vertical nanowire with a small voltage whose magnitude is less than the threshold control voltage (can not change the switch state). We call this voltage v R (Read voltage).

For non-inverting configuration, as depicted in Figure 5.d, the pull-down resistor R G is not required when writing the state of the driving switch to the receiving one. The steps for appropriate latch operation are the same as the 'inverting configuration' except grounding the vertical nanowire of the driving switch rather applying a negative voltage in conditionally close step. In this case, by applying negative voltage to the vertical control line, the receiving switch will be close if the driving switch is close (low-impedance path to the ground) and the receiving switch remains open if the driving switch is open. Thus, the state of the first switch is replicated in the second switch.

Configurable crossbar crray logic gates

In this section, we demonstrate how to make a gate by using the crossbar approach. We start using crossbar to make a 3-input AND gates. Figure 6.a shows a crossbar array which consisting of 8 memristors. We need three inputs A, B, and C, each of which is assumed to be impedance encoded. Basically, various voltages are applied to vertical nanowires to control the memristive switches. The inputs are connected to the horizontal wires as represented in Figure 6.a. The crossbar array functions as an AND gate by applying the following sequences.

1.All junctions are unconditionally opened by applying a positive voltage higher than the threshold voltage to the v in , v and , v out control vertical lines. 2.By driving v in with a negative voltage, the input data (A, B, C) are latched in the input switches. The inputs are in the inverting configuration thus '0' and '1' inputs, close and open the switches, respectively (Figure 6.b). 3.In this step, by driving the v in and K input to the ground, as well as v and to the negative voltage (Fig-ure 6.c) the input data are latched to the wired-AND junctions. 4.In this step, the vertical control lines of the input and output (v in , v out ) are activated by a negative voltage to capture the result to the output memristor (S8). The voltage in point X in Figure 6.d can determine output switch state. If there is at least one closed input (also wired-AND) in the route from v in to the X point, the negative voltage efficacy causes the output switch to stay open. The reason is that the potential difference across the output switch is not enough to close it. It is noteworthy to mention that the voltage at X point is yielded from the voltages dividing between the resistance of switches in the route and S7 that connected to the ground in the last horizontal nanowire line. 5.By driving a positive voltage to v and , all junctions in this column will be opened (Figure 6.e). 6.The AND of the three inputs is stored in the output switch and is ready to be read out.

If the output switch is at inverse mode (inverting configuration), the crossbar array becomes a 3-input NAND gate. Moreover, the crossbar array can also become a 3-input NOR gate, if R G 's are removed while input data are applied to the circuit (Non-inverting configuration). By inverting the last output with the recent situation of the inputs, the crossbar architecture becomes a 3-input OR gate. Consequently, all logic gates are created except XOR and XNOR, which require a little different structures. First, we create an exclusive-NOR gate, subsequently, by inverting the output, XOR is obtained. To implement XNOR, two minterms should be OR'ed together. The crossbar array operation for creating XNOR is similar to AND crossbar array in Figure 6. However, three additional steps are required. To simplify the crossbar array, we apply two variables A and B (Figure 7). In the following, we discuss these sequences step by step.

1.All junctions are unconditionally opened by applying a positive voltage to the v in , two v AN D 's, and v out control lines. 2.By driving v in with a negative voltage, the input data (A, B, Ā, B) are latched in the input switches. The inputs are in the inverting configuration, therefore, '0' and '1' inputs close and open the switches, respectively. 3.In this step, by driving the v in and K input to the ground, and v AN D1 to the negative voltage, the input data are latched to the wired-AND for the first minterm. 4.In this step, the vertical control lines of input and output (v in , v out ) are activated by a negative voltage for first minterm to capture the result ( Ā B) to the output. 5.By applying the positive voltage to v AN D1 (the wired-AND), junctions become open for the first minterm. 6.For the second minterm implementation, step 3 should be repeated by applying a negative voltage to v AN D2 . 7.In this step, the control lines of input and output (v in , v out ) are activated by a negative voltage for the second minterm to capture the result (AB) to the output. It is worth to note that performing OR of the two minterms is sequential rather than concurrent. It means that if the output switch is set into the close state by the first minterm, the value of the second minterm has no effect. Otherwise, changing the state of the output switch is dependent on the second minterm. 8.By applying the positive voltage to the v AN D2 , the wired-AND junctions become open for the second minterm. 9.The XNOR of two inputs A and B is stored in the output switch and it is ready to be read out.

It is also possible to make other boolean logic functions by applying NAND combinations. However, it requires a larger number of memristive switches. Furthermore, it nullifies the most important capability of memristors, i.e. configurability.

Evaluation

In this section, we present an evaluation of both the IMP and crossbar array approaches. Both approaches have significant advantages over the CMOSbased logic gates, such as having computing and storing on the same physical unit, non-volatility, and small scalability because of the nature of memristive switches. By applying the IMP method, the number of memristors to make a gate can be saved. The less number of memristive switches causes less power consumption and cost. Table 2 provides the number of memristors required to make various logic gates in both cases. We note that the number of memristors to make gates in the IMP approach in Table 1 is different from Table 2. When logic operations listed in Table 1 operate on two variables A and B, the original logic would be lost during the implication process. For instance, in case of NAND, only B is changed and in case of XOR both A and B will change. Therefore, to store data we require additional number of memristive switches. For a more detailed explanation see [START_REF] Bickerstaff | Memristor-based arithmetic[END_REF]. Despite the overhead introduced in the crossbar approach in terms of number of memristors, it can dynamically adapt its logic function by controlling voltages. Therefore, based on the run-time requirements, the crossbar array approach is adaptable and reconfigurable.

6 Conclusions In this paper, the electrical model of the memristor is represented to show how physical properties of this device can be utilized for digital circuit applications. Memristor is a two-terminal device, therefore, we need new approaches to apply it to operate as a switch. IMP logic with memristor is studied as the first method to create digital gates. We demonstrated that by using IMP and the inverse function (NOT), it is possible to produce all digital boolean functions. In the IMP approach, the number of memristors to make gates are fewer than the crossbar array approach. We purposed the crossbar array architecture as the second novel promising technique for nanocomputing paradigm. The crossbar array method with memristive switches has been investigated comprehensively. Reconfigurability is the most significant advantage of using the crossbar array architecture. It is interesting to note that in both computing techniques, the storing and computing units are physically the same. This property has the potential to overcome the memory bottleneck. For future work, one can organize a Programmable Logic Device (PLA) platform in such a way that the crossbar array operates as a programmable AND array beside IMP as a fixed connection for the OR array.

Figure 1 .Figure 1

 11 Figure 1 Memristor schematic and behavior: a) the memristor structure, the difference in applied voltage changes doped and undoped regions, b) current versus voltage diagram, which demonstrates hysteresis characteristic of memristor, in the simulation we apply the sinusoidal input wave with an amplitude of 1.5v, different frequencies, RON = 100Ω, ROF F = 15kΩ, D = 10nm, µv = 10 -10 cm 2 s -1 V -1 .

Figure 2

 2 Figure 2 Memristor-based IMP: a) circuit schematic, b) IMP truth table.

Figure 3

 3 Figure 3 NAND configuration with IMP: a) circuit schematic, b) required voltages for controlling the process, c) sequential truth table to obtain NAND.

Figure 4

 4 Figure 4Different states of a switch in a crossbar array.

Figure 5

 5 Figure 5 Different states and configurations of the memristive switches in a crossbar array.

Figure 6

 6 Figure 6 The crossbar array architecture for AND function with memristor switches: a) by applying the positive voltages all switches become open (clear), b) capture data in the input latches, c) transfer data to the wired-AND switches, d) if all inputs are '1' (open) then the X spot is not negative so the output switch has enough voltage across it to be ON, e) open all wired-AND switches to be prepared to read output, f) read output.

Figure 7

 7 Figure 7 Crossbar array architecture for exclusive-NOR function.

Table 1

 1 Different logic operations made by IMP operations and

	the number of required memristors
	Logic Operation Equal IMP functions	#Device
	NOT

Table 2

 2 The number of memristive switches to make logic gates for the imp and crossbar array approaches.

	Logic Operation	#IMP #Crossbar Array
	S ← NOT A	2	3
	S ← A AND B	4	6
	S ← A NAND B	3	6
	S ← A OR B	6	6
	S ← A NOR B	6	6
	S ← A XOR B	7	11

Copyright line will be provided by the publisher


Victor Picheny

Rodolphe Le Riche

Revisiting the paper helicopter project using an adaptive surrogate-based approach

The paper helicopter project has been given to graduate students of the French engineering school Ecole des Mines de Saint Etienne, as a practical work of two lectures: design of experiments and global optimization. The project consists of designing a very simplified paper helicopter for maximum flying time. Its main difficulty lies in the very high level of uncertainty, as flying times for the same design can vary very substantially due to experimental conditions or manufacturing imprecisions. The students were asked to combine design of experiments and surrogate modeling techniques (in order to deal with noise) with global optimization strategies. Despite its apparent simplicity, the project was found quite challenging and the strategies experimented by students were eventually close to what would require a real-case problem.

Problem description

The paper helicopter problem has been initially proposed by [START_REF] Box | Teaching engineers experimental design with a paper helicopter[END_REF] to introduce design of experiments and regression concepts. Since then, it has been re-used and modified many times (see for instance [START_REF] Antony | Teaching the taguchi method to industrial engineers[END_REF][START_REF] Annis | Rethinking the paper helicopter: Combining statistical and engineering knowledge[END_REF][START_REF] Viana | Efficient global optimization with experimental data: revisiting the paper helicopter design[END_REF], as it offers many pedagogical and practical qualities.

We present here a variant of the problem, as it has been given to graduate students of the French engineering school Ecole des Mines de Saint Etienne, as a practical work of two lectures: design of experiments and global optimization.

The helicopter is composed of two superimposed rectangles (a slightly simpler version than the original of [START_REF] Box | Teaching engineers experimental design with a paper helicopter[END_REF]; the upper rectangle is cut in two and folded to constitute the rotor and the lower rectangle constitutes the tail. The rotor wings make a right angle with the tail. For stability, two paper clips are attached at the bottom of the tail (Figure 1). The objective is to design the helicopter for maximum flying time.

The helicopter must be cut from an A4 sheet without gluing; hence, we have the following constraints:

t h + w h ≤ 29.7 cm 2w w ≤ 21 cm
The helicopter should hold two large paper clips (3 × 1cm) placed vertically next to each other, which gives a minimal size for the tail:

t w ≥ 2 × 1 cm t h ≤ 3 cm
To ease in particular the design of experiments process, we decided to simplify those constraints to interval constraints on the design variables:

3.0 ≤ w h ≤ 13.0 cm 1.0 ≤ w w ≤ 7.0 cm 3.5 ≤ t h ≤ 10.0 cm 2.2 ≤ t w ≤ 6.0 cm
The flying time T is measured by throwing the helicopters from a fixed height, as large as possible to minimize timing error (typically, from two storey high), preferably in a wind free environment (a stairwell). This time is a random variable, due mainly to three factors: the experimental conditions (wind, throwing gesture, ...), the measuring error and the manufacturing uncertainties (in particular, the folding angle between the tail and the wings). Hence, we choose as an objective the maximization of the expected flying time T , which is measured by the average time of five throws.

The optimization problem is then expressed as follow: max The overall objective of the project is the global optimization of the paper helicopter based on design of experiments and metamodeling techniques. As the class is divided in groups, each group will experiment a different strategy, all the results being eventually compared. However, all the strategies use the following general scheme:

1. construction of an initial design of experiments (DOE) (50 designs);

2. construction of a metamodel based on the initial set of observations;

3. improvement of the DOE and the metamodel (20 designs);

4. optimization of the helicopter based on the metamodel (5-10 designs).

A "cheaper" variant (timewise) is to consider only three repetitions, 40 initial experiments and 10 improving experiments. Two types of metamodels are proposed: kriging and polynomial response surface (PRS). The DOE improvement and the final optimization strategy depends on the metamodel.

Steps 1, 2 and 4 are standard in surrogate-based optimization, while step 3 is more original. It consists of adding a set of observations that are chosen based on the initial information, in particular in order to improve the quality of prediction of the metamodel. At that stage, the design space can also be reduced if some regions are clearly identified as non-optimal. it would be possible to suppress this step and transfer the experimental effort (20 observations) to either step 1 or step 4; however, this step is motivated by several reasons:

• it is more efficient than step 1 only, as part of the observations are model-oriented and within a smaller region;

• step 4 requires to perform experiments one-by-one, which is very time consuming while they are done all at once in step 3;

• it is also much more robust than step 4 when the number of observations is too small to fit an accurate model.

Construction of the initial design of experiments and metamodel

The students are asked to implement space-filling methods: centroidal Voronoi tesselations (CVT) and latin hyper-cubes designs (LHD) (optimized with respect to the maximin criterion), and generate DOEs of 50 points (in four dimensions). The initial DOE for the project is the best generated according to several criteria (discrepancy, maximin, etc.). An R script that creates a LateX file from a matrix is used to generate the drawings of the helicopters, which are then printed, cut and folded. The initial set of measurements is performed (250 throws), and the 50 observations are taken as the empirical means over the five repeated throws ( Ti , 1 ≤ i ≤ 50). Note that keeping the 250 measurements is also a possibility. To account for the heteroskedasticity of the data, it is also possible to compute the estimated variance of Ti , which is the empirical variance of the five measurements divided by five. This information is not critical when classical linear regression is used, but appears to be very useful for the kriging model.

The initial metamodels are fit using the R functions lm() (library stats) for PRS and km() (library DiceKriging ) for kriging. On each case, the model is adjusted as best as possible to the data, by choosing the polynomial terms based on ANOVA and cross-validation for PRS, and by choosing the trend and covariance functions for kriging.

Improvement of the DOE and the metamodel

For this problem, 50 observations is quite insufficient to build an accurate metamodel. However, it may be sufficient to narrow some of the bounds and remove some of the non-optimal regions, which can be seen as a pre-optimization step.

Criteria The initial DOEs are improved by adding a set of 20 observations using model-oriented criteria, that is, the new observations are added to improve the quality of the initial model. For PRS, the D-criterion (determinant of the Fisher information matrix (see [START_REF] John | D-optimality for regression designs: a review[END_REF]) is used. For kriging, the IMSE criterion is used (see [START_REF] Sacks | Design and analysis of computer experiments[END_REF]). On both cases, the criterion is of the form:

C(x p+1 , . . . , x p+n ) = φ(x 1 , . . . , x p+n ), (1)
where x = (t w , t h , w w , w h ) denotes a design, x 1 , . . . , x n is the current DOE, x p+1 , . . . , x p+n is a set of candidate additional observations, and φ is a scalar function defined by the model. Calculation details are given in 4.

Improving the DOE amounts to solve the following optimization problem: max C(x p+1 , . . . , x p+n )

s.t. x p+i ∈ D ∀1 ≤ i ≤ p. (2)
Optimization of the criteria Solving this problem is challenging, as it can be of high dimension (4 × p), and the objective is strongly multimodal (it is for instance invariant by a permutation of any x p+i , x p+j ). Therefore, it constitues an excellent application for global optimization algorithms.

The students have first implemented different strategies: Nelder-Mead, CMA-ES, DIRECT, etc. Then, they are tested on two analytical functions (a quadratic function, and Michalewicz's function) before applying it to the DOE problem.

Global optimization of the helicopter

Polynomial response surface The final step consists of performing an experiment where the response surface is maximum. The metamodel can be updated with this observation, and the procedure repeated a couple of times if the maximum changes significantly. Finding the maximum of the surface requires solving a global optimization problem (in dimension 4). This is done using the algorithm previously implemented.

Note that here, the strategy assumes that the metamodel is sufficiently accurate to provide a good estimate of the maximizer of the objective function.

Kriging Contrarily to PRS, kriging is well-suited for sequential design, and the final step consists of several iterations of the Efficient Global Optimization (EGO) algorithm (see [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]). Similarly to the PRS case, this strategy amouts to find the maximizer of a function defined by the kriging mean and variance called the Expected Improvement, defined as:

EI n (x) = (m n (x) -y max )Φ m n (x) -y max s n (x) -s n (x)φ m n (x) -y max s n (x) , (3)
where Φ and φ denote the Gaussian cumulative distribution function and probability density, respectively, m n and s n denote the kriging mean and variance and y max is the maximum of the current set of observations. As the observations are noisy, the original EGO strategy should not be applied, so a heuristic modification (called the reinterpolation procedure) is used instead (see [START_REF] Forrester | Design and Analysis of" Noisy" Computer Experiments[END_REF]).

On both cases, as doing experiments one at a time is particularly time-consuming, only a few designs are tested.

Discussion

Scientific challenges and potential advanced topics In our setting, a relatively large proportion of a space-filling design (say, 20 to 30%) experience unstable flights, which result in a lot shorter recorded times. Those unstability are not systematic, and often happen two or three times out of five.

A first consequence is that noise is largely heteroskedastic, which requires a statistical treament. It even seems reasonable to infer a noise model, either continuous (noise varying smoothly with some design variables) or discrete (by partitioning the space into stable/unstable regions).

Another interesting point is that optimal designs seem to be at the limit of stability (helicopter as big as possible but stable). An alternative approach would be to define a constrained optimization problem. However, the stochasticity of the constraint evaluation (stability) makes the problem particularly challenging.

Finally, it seems that a critical factor, hidden in our approach, is the angle between the wing and the body. The students were told to ply the paper to an approximate 90 degree angle, but this is quite difficult to control. As shown in [START_REF] Annis | Rethinking the paper helicopter: Combining statistical and engineering knowledge[END_REF], this angle has a strong influence on the helicopter velocity, and its uncertainty may explain a large proportion of the noise in the observations. Controlling this angle (by changing the helicopter design), or coupling the statistical approach with a physical model [START_REF] Annis | Rethinking the paper helicopter: Combining statistical and engineering knowledge[END_REF] might provide interesting solutions.

Pedagogical aspects Despite its apparent simplicity, this was found quite rich and challenging regarding both statistical and optimization aspects, which makes it an excellent tool for teaching.

First, each observation has an actual cost (timely speaking), hence it really motivates the students to double-check their codes, designs of experiments, etc., before performing any actual experiment. Besides, contrarily to computer experiments for instance, the study cannot be redone, so the students have to deal with the data they have acquired, which sometimes lead to interesting compromises.

Second, this problem share many characteristics with current challenging computer experiments and optimization problems: it can be seen as an expensive black-box, stochastic, multi-physics. Hence, while the primary purpose of the problem was to teach engineers statistics [START_REF] Box | Teaching engineers experimental design with a paper helicopter[END_REF][START_REF] Antony | Teaching the taguchi method to industrial engineers[END_REF], we feel that it can also serve as an introduction to advanced research topics.

Appendix

D criterion

Let the polynomial reponse surface be of the form:

y(x) = p k=1 f k (x) + ε (4)
For a design of experiments x 1 , . . . , x n , the corresponding Fisher information matrix is given by:

M = F T F, (5)
with

F ij = f j (x i ). The D-criterion is then D = det(F T F).
In our context, part of the DOE is fixed (the observations already performed), and we would like to optimize the criterion with respect to a set of p new observations. The form of M allows us to avoid computing the expensive determinant of the new (n + p) × (n + p) matrix for every set of candidate new observations. Indeed, let F n be the part of F corresponding to the current DOE and F p the part corresponding to the candidate points. We have:

F = F n F p (6)

IMSE criterion

The IMSE criterion is defined as the integral of the kriging prediction variance over the design space:

IM SE = D s 2 n (x)dx, (9)
where s 2 n () denotes the prediction variance based on n observations. In our context, p observations are to be added to the design, so the criterion to optimize is:

IM SE(x p+1 , . . . , x p+n ) = D s 2 n+p (x)dx. ( 10
)
The kriging prediction variance doe not depend on the value of the observations, hence s n+p can be computed by updating the current kriging model with dummy observation values. Note that kriging update equations can be use for speed-up, but they are not considered here. The integral over D must be done numerically, using Gauss quadrature or Quasi Monte-Carlo methods for instance. As the integral calculation is embedded in an optimization loop, it is important to always use the same integration points, which creates some bias but avoid having to treat numerical noise in the objective function.

  det(F T F) = det F T n F n + F T p F p (7) = det F T n F n det I p + F p F T n


Ahlem Othmani
email: ahlem.othmani@u-bourgogne.fr

Alexandre Piboule
email: alexandre.piboule@onf.fr

Lew F C Lew

Yan Voon

HYBRID SEGMENTATION OF DEPTH IMAGES USING A WATERSHED AND REGION MERGING BASED METHOD FOR TREE SPECIES RECOGNITION

Keywords: Forest inventory, terrestrial laser scanning, single tree species recognition, 3D geometric texture classification

Tree species recognition from Terrestrial Light Detection and Ranging (T-LiDAR) scanner data is essential for estimating forest inventory attributes in a mixed planting. In this paper, we propose a new method for individual tree species recognition based on the analysis of the 3D geometric texture of tree barks. Our method transforms the 3D point cloud of a 30 cm segment of the tree trunk into a depth image on which a hybrid segmentation method using watershed and region merging techniques is applied in order to reveal bark shape characteristics. Finally, shape and intensity features are calculated on the segmented depth image and used to classify five different tree species using a Random Forest (RF) classifier. Our method has been tested using two datasets acquired in two different French forests with different terrain characteristics. The accuracy and precision rates obtained for both datasets are over 89%.

INTRODUCTION

The development of software tools for the automatic measurement of forest inventory attributes has become a major research field due to the increasing use of Light Detection and Ranging (LiDAR) technology in the forestry domain [START_REF] Dassot | The Use of Terrestrial LiDAR Technology in Forest Science: Application fields, Benefits and Challenges[END_REF]. Recently, the National Office of Forests ("Office National des Forêts -ONF") has started to develop a software tool named "Computree" for the automatic measurement of forest inventory attributes using T-LiDAR Scanner (TLS) data [START_REF] Othmani | Towards Automated and Operational Forest Inventories with T-LiDAR[END_REF]. Among the important attributes that "Computree" should be able to determine is the species of each tree in the ground plot. The identification of the species in a mixed planting is essential for estimating forest inventory attributes according to the species and not in a global manner. But, to the best of our knowledge not much has been done regarding single tree species identification based on TLS data. In this paper we present our work on single tree species identification using ground plot TLS data. The common criteria used to identify the species of a single tree are the shape of the leaves, the general shape of the crown and the variations in geometry across the surface of the bark known as the geometric texture of the bark. Since forest inventory data are mostly acquired during winter in order to reduce occlusions due to leaves and also in order to do the measurement outside the growing period of the trees, it cannot be envisaged to recognize the species based on the shape of the leaves. Finally, the bark is probably the most discriminating feature of the species even if it is subject to changes during the tree's life because of age, injuries and modified growth pattern due to environmental disturbances. Figure 1 shows the bark texture of the five most dominant species in European forest. They are the hornbeam, the oak, the spruce, the beech and the pine.

Figure 1: Example of 3D mesh of the five tree species

One can notice that each of the five species has a distinguishable bark feature: the beech has a relatively smooth surface, the spruce is less smooth compared to the beech and it has circular scars, the pine and the oak are rough with vertical strips but the growth pattern is different, and the hornbeam is smooth with an undulating texture. We thus propose a method that analyzes the geometric variation patterns or geometric texture of the 3D surface of the bark in order to determine the species. The first step of our method, as depicted in section 2, consists thus in generating a 2D depth map representation of the 3D geometric texture from the 3D point cloud of a 30 cm segment of the trunk called a "patch". Then, a hybrid depth image segmentation method is used to extract texture features for tree species classification using the Random Forest technique. Section 3 describes the experimental setting and discusses about the results obtained on two datasets, D1 and D2, consisting of respectively 16 and 30 patches per species.

METHODOLOGY

Our method consists of four steps as represented in Figure 2. Each of the steps, except the first one which is a denoising step using the RapidForm™ software (http://www.rapidform.com) for removing ghost points, will be described in the subsequent sections.

, , is determined using the efficient Aligned Axis Bounding Box (AABB) tree structure [START_REF] Alliez | AABB Tree[END_REF].

The smoothed mesh is computed using Taubin's λ/μ smoothing algorithm [START_REF] Taubin | Geometric Signal Processing on Polygonal Meshes[END_REF]. It consists in basically performing the Laplacian smoothing two consecutive times with different scaling factors noted λ and μ. A first step with λ > 0 (shrinking step) and a second step with a negative scaling factor μ<-λ<0 (unshrinking step). Laplacian smoothing consists in iteratively moving each of the vertices of the mesh to a new position that corresponds to the weighted average position of the neighboring vertices. The new position v' i of a vertex i is given by:

(3)
where v i is the current position, λ is a scalar that controls the diffusion speed and ∆ the Laplacian operator, a weighted sum of the difference between the current vertex v i and its neighbors v j , given by the following equation:

, ( 4
)
where i* is the set of all the neighbors of the vertex v i .

Taubin's λ/μ smoothing algorithm is run with equal weights w i,j for each of the neighbors such that ∑ , 1 and with λ and μ respectively equal to 0.6307 and -0.6732 (values suggested by Taubin). It is run iteratively until the mesh is sufficiently smoothed. To measure the smoothness of the mesh the maximum curvature is computed for all the points of the mesh and the median of the maximum curvature values considered. An experimental study that we have conducted showed that the mesh is sufficiently smoothed when the difference in the median of the maximum curvature values between two successive runs of Taubin's algorithm is less than or equal to 1%. We thus use this criterion as the smoothing stopping criterion.

The geometric details, a 3D map or a dataset (x, y, z, d) in the space, is next transformed into a depth image, a dataset in the (X, Y, I) space that can be modeled as:

, , , (5)
where d i is the details for a point v i of coordinates , , given by equation (1) and , , ,

. is the Isomap dimensionality reduction algorithm proposed by Tenenbaum et al [START_REF] Tenenbaum | A Global Geometric Framework for Nonlinear Dimensionality Reduction[END_REF] followed by a quantization step. The latter step is necessary since Isomap transforms nonlinear manifolds in to Euclidean manifolds in whereas we would like to transform the 3D coordinates of a non linear manifold to the 2D coordinates (X, Y) ∈ of the pixels of the depth image. Isomap is particularly suitable in our case since it minimizes the geodesic distortion between points on the 3D surface which results in a minimization of shape distortion when transforming from to In equation ( 5) the mean of the distances d i is computed since there may be several 3D points , , that map to the same pixel at coordinates (X, Y) in the depth image. Also, it is possible that no 3D point , , maps to a pixel of coordinates (X, Y). In that case a bicubic interpolation is done in order to estimate a value for the pixel. Examples of the depth images in the gray level scale [0, 255] for each of the five species to identify are shown in Figure 4.

Hybrid edge-region segmentation

The proposed method is a 2-step approac uses watershed to produce a primary segm input depth images, while the second weakness boundary based region-mergin obtain the final segmented map [START_REF] Haris | Hybrid Image Segmentation Using Watersheds and Fast Region Merging[END_REF]. The flow proposed method is shown in Figure 5. where F is the depth image and ○ is th opening operation with an octagonal struc of size 30 pixels from the origin to the side It is a contrast detector suitable for enhanci in the image. A post processing stage removing small objects from the image and by an opening operation followed by a c resulting segmented tree trunk depth imag Figure 6.

Region merging based on we

Despite all the processing in the overcome the over-segmentation pr neighboring regions that could b meaningful region. In our case, t segmented adjacent regions R 1 homogeneous, but, they nevertheles region. An edge dissimilarity mea weakness of the boundaries will b criteria. The Region Adjacency Gra in Figure 7, will be used for the me partition is defined as an undirected R = {1, 2, .., K} is the set of nodes a edges. Each region is represented adjacent regions R i and R j are rep i, j ϵ R and connected by the edge measure between the two regions the graph edge. The region-merging algorithm segmentation. The corresponding e image is extracted. Two adjacent re E i,j lower than a fixed threshold (0.5

,
where n is the number of pixels in adjacent regions, is the number black one (grayscale < 0.8), is considered as white pixels (graysc = 0.25. The final segmentation is

Feature extraction and classi

We assign each region a set quantitative attributes that describe directionality. We used the regionpr to calculate the properties of r perimeter, orientation, and so on descriptors, we used other shape circularity ratio, rectangularity, co solidity as described in [START_REF] Mingqiang | A Survey of Shape Feature Extraction Techniques[END_REF].

eak boundaries initial segmentation to roblem, there still remain be merged into a more there can be two overand R 2 that are not ss belong to the same real asure that evaluates the be used as the merging aph (RGA), as illustrated erging Severe injury occurs during the life of the tree which can damage the tree bark and change its discriminating criteria. Also, the size of the scars depends strongly on the age of the tree and the distance to the scanner which can result in a dispersion of the features. To evaluate the intraclass variability, we consider the mean, the standard deviation and the median of each feature of all the regions. Classification is done using the Random Forest (RF) classifier proposed by Breiman [START_REF] Breiman | Random Forests[END_REF] in order to classify the textures and thus the tree species.

EXPERIMENTAL RESULTS AND CONCLUSION

We used two different datasets D1 and D2 to validate our approach experimentally. The test site of D1 is a mixed state forest of Montiers-sur-Saulx, France. The second test site of D2 is a mixed forest stands of grove and coppice-undergrove forests in Lorraine, France. The T-LiDAR data was collected with a Faro photon 120™ scanner. A segmentation of the 3D plots as described in [START_REF] Othmani | Towards Automated and Operational Forest Inventories with T-LiDAR[END_REF] is done and 3D patches are extracted. The two datasets contain five different species of trees (Figure 1): 16 patches per species in D1 and 30 patches per species in D2. Two tests are done and in each test we train a separate RF classifier built with 1000 decisions trees. For both D1 and D2, we split the dataset into 10 equal size subsets and performed a 10-fold crossvalidation. The confusion matrices are reported in Tables 1 and2 For both datasets, each with different terrains and architecture characteristics, we obtain good accuracy and precision rates. The achieved overall accuracy and precision in tree species recognition are respectively 95% and 89% for D1 and 97% and 92% for D2. We can differentiate two families from the five tested species (species that have straps or cracks: hornbeam, oak and pine) and the family of random shapes (smooth species as beech and spruce). We note from the confusion matrices that our algorithm mainly mixes intra-family species but also inter-family species. Future work involves finding better features to discriminate inter and intra family species and the validation of the performance of our method when changing the distance to the scanner.

Figure 2 :

 2 Figure 2: Flow diagram of the method 2.1 Depth image computation After the denoising step the 3D point cloud of the tree trunk segment is transformed into a 2D depth image that represents the 3D geometric details of the bark. 3D geometric details are defined as the local variations of the original mesh M o with respect to a smoothed version of the same mesh M s as represented in Figure 3.

Figure 3 :

 3 Figure 3: Illustration of 3D geometric detailsThus, the 3D geometric details is a dataset (x, y, z, d) in the space where (x, y, z) are the coordinates of the 3D points or vertices of the smoothed mesh and d, known as the details, the Euclidean distance between a point of the smoothed mesh Ms of coordinates (x, y, z), denoted by , , , and its nearest neighbor, denoted by , , , in the original mesh Mo. ,

Figure 4 :

 4 Figure 4: Examples of depth im2.2 Hybrid edge-region segmentationThe proposed method is a 2-step approac uses watershed to produce a primary segm input depth images, while the second weakness boundary based region-mergin obtain the final segmented map[START_REF] Haris | Hybrid Image Segmentation Using Watersheds and Fast Region Merging[END_REF]. The flow proposed method is shown in Figure5.

Figure 5 :

 5 Figure 5: Flow diagram of the segmenta 2.2.1 Initial segmentation using Watershed Watershed transformation [6, 7] is a powerf segmentation. It considers the gray level altitude in the relief which is exactly th images. This is why we have chosen this me The watershed transform of the gradie the image very often leads to an over-segm image. It gives poor results in the case of t of tree barks even when foreground object markers are used because of noises in th Consequently, we use a function based o criterion (contrast between facets due to b the watershed. It is a morphological trans "Top-Hat" transformation which is defined b

Figure 6

 6 Figure 6: Watershed segmented

Figure 7 : 1 )

 71 Figure 7: 1) Five-partition o corresponding RAG, and 3) m

Figure 8 :

 8 Figure 8: Segmented images a2.3 Feature extraction and classiWe assign each region a set quantitative attributes that describe directionality. We used the regionpr to calculate the properties of r perimeter, orientation, and so on descriptors, we used other shape circularity ratio, rectangularity, co solidity as described in[START_REF] Mingqiang | A Survey of Shape Feature Extraction Techniques[END_REF].

  . The RAG of a Kd graph, G = (R, E) where and E ⊂ R×R is the set of d by a graph node. Two presented by two nodes e (i, j). The dissimilarity represents the weight of of an image, 2) the merging of two nodes.starts with the initial edge map in the original egions with an edge cost 5 in our case) are merged.(8)the edge linking the two r of pixels considered as s the number of pixels cale > 0.8), = 0.75 and shown in Figure8. after region merging ification of values representing es the geometry and the rops() function of Matlab regions including area, n. In addition to those e features: compactness, oncavity, and convexity,

Table 1 : Confusion matrix for D1 cross-validation (1: hornbeam, 2: oak, 3: spruce, 4: beech, 5: pine)

 1 .

		1	2	3	4	5
	1 14.428	0.571	0	0	1
	2 0.428	12.571	1.714	0	1.285
	3 0.285	0	15.142	0.428	0.142
	4	0	0	0.714	15.142	0.142
	5 0.714	0.142	0.857	0	14.285
		1	2	3	4	5
	1 27.166	0.166	0	1.833	0.833
	2	0	29	0	0	1
	3 0.166	1.166	27.333	1.333	0
	4	1	0	1.333	27.333	0.333
	5	0	0	0.166	1.5	28.333

Table 2 : Confusion matrix for D2 cross-validation

 2

ACKNOWLEDGEMENT

This work is financially supported by the "Conseil Regional de Bourgogne" under contract N° 2010-9201AAO048S06469 and N° 2010-9201CPERO007S06470, and the "Office National des Fôrets".


Ahlem Othmani
email: ahlem.othmani@u-bourgogne.fr

Alexandre Piboule
email: alexandre.piboule@onf.fr

Oscar Dalmau
email: dalmau@cimat.mx

Nicolas Lomenie
email: nicolas.lomenie@mi.parisdescartes.fr

Lew F C Lew

Yan Voon

Tree species classification based on 3D bark texture analysis

Keywords: Tree species classification, 3D pattern recognition, 3D bark texture analysis, forest inventory

Terrestrial Laser Scanning (TLS) technique is today widely used in ground plots to acquire 3D point clouds from which forest inventory attributes are calculated. In the case of mixed plantings where the 3D point clouds contain data from several different tree species, it is important to be able to automatically recognize the tree species in order to analyze the data of each of the species separately. Although automatic tree species recognition from TLS data is an important problem, it has received very little attention from the scientific community. In this paper we propose a method for classifying five different tree species using TLS data. Our method is based on the analysis of the 3D geometric texture of the bark in order to compute roughness measures and shape characteristics that are fed as input to a Random Forest classifier to classify the tree species. The method has been evaluated on a test set composed of 265 samples (53 samples of each of the 5 species) and the results obtained are very encouraging.

Introduction

TLS is today a well established technique for the acquisition of precise and reliable 3D point clouds from which forest inventory attributes can be calculated at the single tree level [START_REF] Dassot | The Use of Terrestrial LiDAR Technology in Forest Science: Application fields, Benefits and Challenges[END_REF] [START_REF] Othmani | Towards Automated and Operational Forest Inventories with T-LiDAR[END_REF]. Numerous work on the calculation of the Diameter at Breast Height (DBH), the height of the tree, the volume of wood and so on can be found in the literature. However, to the best of our knowledge not much has been done in the field of tree species recognition at the single tree level although it is a very important issue if one would like to analyze the 3D data of each of the species of a mixed planting separately. One can only find some work concerning tree species recognition using a combination of TLS data and hyperspectral [START_REF] Puttonen | Tree species classification from fused active hyperspectral reflectance and LIDAR measurements[END_REF] or panoramic images [START_REF] Haala | Combination of terrestrial Laser Scanning with high resolution panoramic Images for Investigations in Forest Applications and tree species recognition[END_REF] [START_REF] Reulke | Tree species recognition with fuzzy texture parameters[END_REF]. The aim of our work is to recognize tree species based on TLS data only for two major reasons. The first one is to save acquisition time and/or to avoid the use of multiple acquisition systems for capturing different types of datasets. The second reason is to avoid the need for data coregistration and/or fusion techniques in order to simplify data processing.

Using TLS data only, we can only analyze 3D geometrical shape features such as the shape of the leaves, the general shape of the crown and the variations in geometry across the surface of the bark known as the 3D geometric texture of the bark in order to recognize the tree species. In our case, we have to exclude leaf shape analysis because our forest inventory data are mostly acquired during winter when the trees are leafless. The reason is to reduce occlusions due to leaves for more accurate wood volume calculation and also to do the measurement outside the growing period of the trees. The general shape of the canopy is a good characteristic feature of the species of a tree for isolated trees. However, in a forest planting, the management type and the density have a big impact on the canopy. They render it polymorphic so that it is difficult to use the shape of the canopy as a discriminating criterion for tree species recognition. Finally, the bark is probably the most discriminating feature of the species even if it is subject to changes during the tree's life because of age, injuries and modified growth pattern due to environmental disturbances.

Fig. 1 illustrates the 3D geometric texture characteristics of the bark of the five most important tree species that we have to recognize. One can notice that each of the five tree species has a distinguishable 3D geometric bark texture characteristic. The beech has a relatively smooth surface, the spruce is less smooth compared to the beech and it has circular scars, the pine and the oak are rough with vertical strips but the growth pattern is different, and the hornbeam is smooth with an undulating texture. We thus propose a method that analyzes the 3D geometric texture of the bark in order to classify and recognize the tree species. For the analysis a 30 cm long segment of the tree trunk at about 1.3 m from the ground (breast height) called a patch is considered. 30 cm is a good trade-off between a small patch for rapid processing times and a long enough segment that contains sufficient texture patterns for recognition.

Our method consists of several steps. Firstly, a 3D deviation map is computed from the 3D point cloud of the 30cm long segment of the tree trunk at breast height. The first step is described in details in section 2. Secondly, the 3D deviation map is transformed into a 2D deviation map or height map that is next segmented in order to reveal characteristic shape features of the tree species. The second step is presented in section 3. Finally, classification features are computed from both the 2D deviation map and the segmented 2D deviation map and fed as input to a Random Forest classifier for tree species classification. The classification features and the selection of the most pertinent features are presented in section 4. In section 5, we describe the test set used to evaluate our method and discuss about the classification results obtained before concluding in section 6.

3D deviation map of the tree bark

We define the 3D geometric texture of a 3D surface as the local variations of the original meshed surface denoted by M o with respect to a smoothed version of the same meshed surface denoted by M s as illustrated in Fig. 2 for some tree trunk segments. To extract the 3D geometric texture of a patch, we first apply a denoising filter, implemented in the RapidForm software, to the 3D point cloud in order to remove ghost points. Next, the 3D point cloud is meshed and smoothed. Finally, the deviation between the original mesh and the smoothed mesh is computed. This yields a 3D deviation map representation of the geometric texture that can be modeled by a dataset DM 3 = {(x, y, z, d) : x, y, z, d ∈ R} where (x, y, z) are the coordinates of the 3D points or vertices of the smoothed mesh and d, the Euclidean distance between a point of the smoothed mesh M s of coordinates (x, y, z), denoted by v(x, y, z), and its nearest neighbor in the original mesh M o , denoted by ṽ(x , y , z ). d is computed according to equation 1.

d(v, ṽ) = (x -x ) 2 + (y -y ) 2 + (z -z ) 2
(1) ṽ(x , y , z ) is determined using the efficient Aligned Axis Bounding Box (AABB) tree structure [START_REF] Alliez | AABB Tree[END_REF]. It is equal to

ṽ = argmin p i ∈Mo p i -v (2)
The smoothed mesh is computed using Taubin's λ/µ smoothing algorithm [START_REF] Taubin | Geometric Signal Processing on Polygonal Meshes[END_REF]. It consists in basically performing the Laplacian smoothing two consecutive times with different scaling factors denoted by λ and µ. A first step with λ > 0 (shrinking step) and a second step with a negative scaling factor µ < -λ < 0 (unshrinking step). Laplacian smoothing consists in iteratively moving each of the vertices of the mesh to a new position that corresponds to the weighted average position of the neighboring vertices. The new position v' i of a vertex i is given by Eq. 3:

v i = v i + λ ∆v i (3)
where v i is the current position, λ a scalar that controls the diffusion speed and ∆v i the Laplacian operator given by Eq. 4. It is a weighted sum of the difference between the current vertex v i and all its neighbors v j .

∆v i = j:vj ∈i * w i,j (v j -v i ) ( 4
)
where i * is the set of all the neighbors of the vertex v i . Taubin's λ/µ smoothing algorithm is run with equal weights w i,j for each of the neighbors such that j:vj ∈i * w i,j = 1 and with λ and µ equal to 0.6307 and -0.6732, respectively (values suggested by Taubin). It is run iteratively until the mesh is sufficiently smoothed.

The smoothness of a 3D surface is usually quantified by the minimum, maximum, mean and Gaussian curvatures of each of the points of the mesh. Consequently, we have studied the evolution of the median of these curvature values as a function of the number of iterations of Taubin's algorithm for several samples of the five species to classify. Our study showed that there is no significant difference in the curvature against number of iterations curve for the four types of curvature values. We have thus decided to consider only the median of the maximum curvatures curve, as shown in Fig. 3, to determine the smoothing stopping criterion.

We have chosen to stop the smoothing process when the slope of the tangent to the curve is less than or equal to -0.01. For this value of the slope we have noticed that the mesh is sufficiently smoothed while preserving the main structures of the trunk. One can find in the literature several non linear techniques that are more appropriate for our 3D deviation map. Among all these techniques we have chosen the Maximum Variance Unfolding (MVU) dimensionality reduction algorithm proposed by Weinberger et al. [START_REF] Weinberger | Nonlinear dimensionality reduction by semidefinite programming and kernel matrix factorization[END_REF] because it is fast and it gives good results. An example of the 2D deviation map obtained using MVU is illustrated in Fig. 4. The 2D deviation map is a set of points in a 3D space. It is different from the point cloud in the sens that the third dimension represents a distance from a plane surface defined by the other two dimensions. It is like a height map or relief map defined by a set of points. The idea now is to cluster the points that are above a certain distance or "height" from the plane in order to define clusters or regions whose shape will allow us to classify the different species. The result of the clustering is a segmented 2D deviation map. To acheive this, a first thresholding step is done in order to keep the most salient features, e.g. the points with the highest deviation values. The threshold value is empirically defined as the median value of the deviation values. In this way only all the points in the 2D deviation map that have a height value greater than the median value are kept.

Secondly, DBSCAN a density based algorithm for discovering clusters proposed by Ester et al. [START_REF] Ester | A density-based algorithm for discovering clusters in large spatial databases with noise[END_REF] is used to cluster the points in order to build regions. DBSCAN is chosen because it is a fast and efficient algorithm even for large spatial set of points which is our case. DBSCAN is based on the notion of density of points in an -neighborhood. It has two required parameters: the neighborhood size ( ) and the minimum number of points in the -neighborhood (minP ts). For our application = 0.6 and minP ts = 8.

An example of the segmented 2D deviation map for each of the five species is illustrated in Fig. 5.

Classification features

From the 2D deviation map and the segmented 2D deviation map, we can compute about 128 features for classification. These features are, for example :

-The number of clusters per segmented 2D deviation map.

-The roughness features described in section 4.1.

-The principal component analysis features described in section 4.2.

-The shape and intensity features presented in section 4.3. Each cluster is defined by a set of shape and intensity features. We calculate the median, the mean, the standard deviation, the minimum and the maximum of the intensity and the shape features of all the clusters and use them as classification features.

Roughness features

The 2D deviation map represents the geometric details of the surface of the bark from which roughness measures can be computed. Many types of measures can be found in the literature. The most common ones are statistical values such as the root mean square value, the arithmetical mean and the standard deviation respectively given by Eqs. ( 5)- [START_REF] Taubin | Geometric Signal Processing on Polygonal Meshes[END_REF]. 6)

s q = 1 n n i=1 d(v i , ṽi ) 2 (5)
s a = 1 n n i=1 |d(v i , ṽi )| (
s D = 1 n n i=1 d(v i , ṽi ) 2 - 1 n n i=1 d(v i , ṽi ) 2 (7)

Principal component analysis features

Principal Component Analysis (PCA) is a technique for dimension reduction and feature extraction. It uses linear transformations to map data from a high dimensional space to a low dimensional space. The new variables in the low dimensional space are called principal components. Some of the features that can be calculated for each of the clusters of points of the segmented deviation map using PCA are listed below.

-The percentage of the total variance explained by each principal component.

-The maximum and the median distances between the observations and the center of the data set as well as the ratio between the maximum and the median distances. -The longest and shortest diameter (length of the major and the minor axis).

-The aspect ratio defined as the major axis length divided by the minor axis length. -The orientation or the direction of the first and the second principal component.

The minimum, maximum, mean, median and standard deviation values of these PCA features for all the clusters of the segmented 2D deviation map are used as features for classification.

Shape and intensity features from segmented 2D deviation map

To characterize the clusters of points in the segmented deviation map, shape and intensity features are calculated for each cluster. The intensity features are the maximum and the median intensity of the points of the cluster. To calculate the shape features, the alpha shape (α-shape) algorithm proposed by Edelsbrunner et al. [START_REF] Edelsbrunner | On the Shape of a Set of Points in the Plane[END_REF] is used to compute the ∞-shape that corresponds to the convex hulls of the clusters of points. The 0.6-shape that best represents the real shape has also been computed. From the ∞-shape and the 0.6-shape, shape features such as perimeter, area, solidity, roundness, compactness, RFactor, shape, convexity and concavity are computed and used as classification features. The mathematical expressions of some of these features are given by Eqs. ( 8)-(12). To discriminate between the strips of the pine that are close and narrow and those of the oak that are touching each other, the ratio between the area of the convex hull and the area of the 0.6-shape is used as feature.

Solidity = Area ConvexArea (8)
Roundness = 4 × Area π × M ajorAxisLength 2 (9)
Compactness = 4 π × Area M ajorAxisLength (10)
RF actor = P erimeterOf ConvexHull M ajorAxisLength × π (11)
Shape = P erimeter 2 Area (12)

Experiments and results

The 3D TLS data used to evaluate our method were captured using either a FARO Photon 120 or a FARO Focus3D scanner with a resolution of about 6mm at 10m. All trees at a distance of about 6m and with a Diameter at Breast Height (DBH) of about 30cm are manually located in the scans. Next, non occluded tree trunk segments of about 30cm long and at a height of about 1.30m from the ground are extracted using a software that we have developed, the "Computree" software [START_REF] Othmani | Towards Automated and Operational Forest Inventories with T-LiDAR[END_REF], in order to constitute the evaluation database. We used two different datasets D1 and D2 to validate our approach experimentally. The test site of D1 is a state mixed forest in Montiers-sur-Saulx, France. The second test site of D2 is a mixed forest stands of grove and coppiceunder-grove forests in Lorraine, France. Both datasets contain the five different species of trees represented in Fig. 1. D1 and D2 are composed of 20 patches per species and 33 patches per species, respectively. Classification is done using the R Language implementation of the Random Forest (RF) classifier proposed by Breiman [START_REF] Breiman | Random Forests[END_REF]. The RF classifier is built with recommended values for the number of decision trees (1000) and the number of features used to split the node in the decision tree growing process denoted by M try (M try = √ D where D is the feature vector size). We have tested 128 features and used RF to select the 30 most pertinent ones to evaluate our method. So, the value of D is 30.

Three tests are done using datasets D1, D2 and a combination of D1 and D2, respectively. For each test, stratified 10-fold cross-validation is performed 1)-( 3) and the accuracy rates for the three tests are summarized in Table 4. For both datasets, each with a different terrain and architecture characteristics, we obtained good classification results. One can note from the confusion table that the accuracy ranges from 83% to 100%. The worst accuracy of 83% is for the pine that is misclassified as any one of the other four species in the three tests. Pine is mainly confused with hornbeam: the two species have strips. We can differentiate two families of tree species from the five tested species: species that have straps or cracks such as the hornbeam, the oak and the pine, and smooth surface species such as the beech and the spruce. We note that our algorithm mixes mainly intra-family species (spruce and beech) but also inter-family species. Future work will focus on finding more pertinent features to discriminate inter and intra family species.

We can note that training and testing with dataset D3 (test 3) gives overall better accuracy rate than the other tests. This may be due to the fact that D3 contains more samples. Nevertheless, the result is relatively good enough and we can consider that the classifier performs reasonably well since the datasets D1 and D2 are composed of patches extracted from 3D data acquired from two different forest sites.

Conclusions

We have proposed a method for classifying five different tree species using TLS data only. The method is based on the analysis of the 3D geometric texture of the bark in order to identify the species class to which pertains the tree under analysis. In our method, the 3D geometric texture is transformed into a 2D deviation map on which roughness measures and shape features are computed and used as features for classification using the Random Forest classifier. Results obtained on a dataset composed of 265 samples with equal number of samples of each species are quite good (83% to 100%). In future work we plan to do more tests with other datasets in order to verify this observation. Moreover, we would also like to study the influence of the distance to scanner and of the DBH on the results. Indeed, in the current dataset the samples are all extracted from trees located at about 6m from the scanner and with a DBH of about 30cm.

Fig. 1 .

 1 Fig. 1. 3D point clouds of the five tree species

Fig. 2 .

 2 Fig. 2. 3D geometric texture model

Fig. 3 .

 3 Fig. 3. The median of the maximum curvatures against the number of iterations for each of the five tree species

Fig. 4 .

 4 Fig. 4. Example of 2D deviation map for each of the five species

Fig. 5 .

 5 Fig. 5. An example of segmented 2D deviation map for each of the five species

Fig. 6 .

 6 Fig. 6. (a) 3D point cloud captured by FARO Photon 120 scanner (b) Computree Software screenshot

Table 1 .

 1 Confusion matrix for D1 cross-validation (1: hornbeam, 2: oak, 3: spruce, 4:

	beech, 5: pine)
		1	2	3	4	5	Accuracy
	1	20	0	0	0	0	100 %
	2	0	20	0	0	0	100 %
	3	0	0.8	19.2	0	0	96 %
	4	1.8	0	0.2	18	0	90 %
	5	2.2	0	1	0.2	16.6	83 %

Table 2 .

 2 Confusion

		matrix for D2 cross-validation (1: hornbeam, 2: oak, 3: spruce, 4:
	beech, 5: pine)
		1	2	3	4	5	Accuracy
	1	33	0	0	0	0	100 %
	2	0.6	31.8	0	0	0.6	96.3 %
	3	0	0	31.6	0.6	0.8	95.7 %
	4	1.8	0	0.6	32.4	0	98.1 %
	5	0.8	1	0	0	31.2	94.5 %

Table 3 .

 3 Confusion matrix for D1 and D2 cross-validation (1: hornbeam, 2: oak, 3:

	spruce, 4: beech, 5: pine)
		1	2	3	4	5	Accuracy
	1	53	0	0	0	0	100 %
	2	0	52.6	0	0	0.4	99.3 %
	3	0	0	51.6	0.4	1	97.4 %
	4	0.4	0	0.4	52.2	0	98.5 %
	5	1.4	0.4	0.4	0.4	50.4	95.1 %

Table 4 .

 4 Accuracy rates for the three tests

		Min	Average	Max	σ
	Test 1	83%	93.8%	100%	6.5%
	Test 2	94.5%	96.92%	100%	1.9%
	Test 3	95.1%	98.06%	100%	1.7%

Acknowlegments

This work is financially supported by the "Conseil Régional de Bourgogne" under contract Nº 2010 9201AAO048S06469 and Nº 2010 9201CPERO007S06470, and the 'Office National des Fôrets".