Farah Mourad
email: farah.mourad@utt.fr.

Paul Honeine
email: paul.honeine@utt.fr.

Hichem Snoussi
email: hichem.snoussi@utt.fr.

Polar interval-based localization in mobile sensor networks

Keywords: mobile sensor networks, multi-coordinate systems, polar intervals, state estimation. I

This paper considers the problem of localization in uncontrolled mobility sensor networks. Based on connectivity measurements, the problem is solved using polar intervals. Computation is performed, in several polar coordinate systems, using both polar coordinates and interval analysis. Position estimates are thus partial rings enclosing the exact solution of the problem. Simulation results corroborate the efficiency of the proposed method compared to existing methods, especially to those handling single coordinate systems.

Of Arrival (TDOA) [START_REF] Eickhoff | 3D-accuracy improvements for TDoA based wireless local positioning systems[END_REF], Angle Of Arrival (AOA) [START_REF] Tai | Three-dimensional non-line-of-sight localisation in an indoor multipath environment[END_REF], or Received Signal Strength Indicator (RSSI) [START_REF] Zanca | Experimental comparison of RSSI-based localization algorithms for indoor wireless sensor networks[END_REF], [START_REF] Li | Research on indoor wireless localization system for radioactive sources based on zigbee[END_REF]. TOAbased methods measure the travel times of signals exchanged between the sensors, whereas TDOA and AOA methods measure the difference of arrival times or the angles at reception of exchanged signals. While TOA, TDOA and AOA-based methods need extra hardware such as timers and synchronous clocks, RSSI-based methods are simpler, inexpensive and lessenergy consuming using only powers of exchanged signals, but they are more challenging because of the reflection, the diffraction and the scattering of signals.

The alternative category to range-based localization consists of range-free algorithms. These techniques yield coarsegrained location estimates, leading to more than one possible solution to the problem. In this category, one could refer to connectivity-based techniques, where RSSI is compared to a strength threshold leading to bounds over distances [START_REF] Mourad | Anchor-based localization via interval analysis for mobile ad-hoc sensor networks[END_REF]. Other range-free techniques are based on hop counts between anchors and nodes. Having per-hop distance, obtained through anchors communication, one could estimate the distances separating nodes from anchors [START_REF] Niculescu | Dv based positioning in ad hoc networks[END_REF], [START_REF] Yi | Multi-hop-based monte carlo localization for mobile sensor networks[END_REF]. Distances information are then combined using different computation techniques, such as the Robust Optimization approach [START_REF] Bertsimas | Theory and applications of Robust Optimization[END_REF], the particle filter based on Monte-Carlo [START_REF] Hu | Localization for mobile sensor networks[END_REF], the variational filter [START_REF] Teng | Decentralized variational filtering for simultaneous sensor localization and target tracking in binary sensor networks[END_REF], interval analysis [START_REF] Mourad | Interval-based localization usinf RSSI comparision in MANETs[END_REF], Sigma-Point Kalman Smoothers as in [START_REF] Paul | RSSI-based indoor localization and tracking using Sigma-Point Kalman Smoothers[END_REF], etc.

This paper considers the problem of localization in uncontrolled mobility sensor networks. The proposed technique is a range-free anchor-based method. Using connectivity measurements, the proposed method also takes advantage of the mobility of the nodes to set the problem. The novelty of this method is to use polar intervals to solve the problem. Indeed, position estimation is performed using polar coordinates in the interval framework. Estimates are thus partial rings, called polar boxes, enclosing the exact solution of the problem. Based on interval analysis [START_REF] Jaulin | Applied interval analysis[END_REF], computation is handled in four different polar coordinate systems, leading to four position estimates at each time-step. Only one polar box is then selected at a considered step. It corresponds to the estimate having the smallest area, being the solution encloser presenting the least incertitude. The three remaining polar boxes are also kept in the memory to be used in the following time-step computation. Simulation results using Matlab corroborate the efficiency of the proposed method compared to other methods, especially to those handling computation in a single coordinate system.

The rest of the paper is organized as follows. Section II introduces the localization problem. Section III describes the proposed algorithm to solve the problem. Simulation results are given in Section IV whereas Section V concludes the paper.

II. PROBLEM STATEMENT

Assume that the network is composed of N a anchors and N u nodes, deployed in a two-dimensional square area, denoted by Ω and having u 0 as side lengths. Let a i,t , i ∈ {1, ..., N a }, and u j,t , j ∈ {1, ..., N u }, be the respective positions of these sensors at time t. The aim of the method is to estimate all u j,t positions using a i,t information, as well as previous estimates. In order to reduce the communication costs, one assumes that nodes exchange information only with anchors. For this reason and without loss of generality, only one node u t is considered in this paper and the index j is withdrawn. This paper handles computation using polar coordinates of sensors. In other words, the aim of the method is to estimate the polar coordinates of the node using the polar coordinates of some anchors. For this reason, one Polar Coordinate System (PCS), denoted by PCS , is first considered. Without loss of generality, the origin O of PCS is located at the low-left corner of the surveillance area Ω. An illustration is shown in Fig. 1. Here, the sensors polar coordinates are given by a i,t = ρ i,t , θ i,t and u t = (ρ t , θ t ), where ρ denotes the distance from the origin O to the sensor and θ denotes the angle measured anticlockwise from the horizontal x-axis to the line joining O to the sensor. According to this definition,

θ i,t ∈ [0, π 2 ], θ t ∈ [0, π 2 ], ρ i,t ∈ [0, ρ m ] and ρ t ∈ [0, ρ m ],
where ρ m = √ 2.u 0 is the length of the diagonal of the surveillance area. It is worth noting that an infinite number of PCSs could be considered to define the localization problem. In the following, the localization problem is first described in PCS , and afterwards, it is reformulated in other PCSs.

A. Problem statement in PCS

The proposed method takes advantage of the mobility of the node to estimate nodes positions. The localization problem is thus defined using a mobility model, in addition to observation constraints. Any available information about the mobility of the node could be used to set the mobility constraints. This paper considers a general model assuming that the maximal velocity of the node is known and denoted by v m . If ∆t is the localization period, then the maximal distance that could be traveled by the node between two consecutive time-steps is equal to D m = ∆t.v m . The mobility constraint is then given by the following,

d t ≤ D m , (1)
where d t is the value of the distance traveled by the node between time-steps t -1 and t. Let u t-1 = ρ t-1 , θ t-1 be the position of the node in PCS at time t -1. Then d t expressed in PCS is the length of the side u t-1 u t of the triangle having the origin O , u t and u t-1 as vertices. In this triangle, the lengths of O u t and O u t-1 are given by ρ t and ρ t-1 respectively, as shown in Fig. 1. The distance d t could then be computed using the generalized Pythagorean theorem as follows,

d 2 t = ρ t 2 + ρ t-1 2 -2ρ t ρ t-1 cos θ t -θ t-1 , (2)
where |θ t -θ t-1 | is the angle measured at the vertex O in the considered triangle. The mobility constraint is then rewritten as follows,

u t u 0 u 0 θ = 0 θ = π 2 u t-1 a i,t O Ω ρ i,t ρ t ρ t-1 θ i,t θ t θ t-1 d i,t d t
ρ t 2 + ρ t-1 2 -2ρ t ρ t-1 cos θ t -θ t-1 ≤ D 2 m . (3)
Graphically, this constraint yields a disk, called mobility disk, having the previous position as center and D m as radius.

In addition to the mobility model, the proposed method is a range-free method based on RSSI information. At each time step, every anchor broadcasts signals in the network with the same initial power. According to the Okumura-Hata model [START_REF] Medeisis | On the use of the universal okumura-hata propagation prediction model in rural areas[END_REF], the strengths of the signals decrease monotonically with the increase of their traveled distances as follows,

ξ i,t = ξ 0 d 0 d i,t α , (4)
where ξ i,t is the strength of the signal emitted by the anchor i and received by the mobile node at time t, ξ 0 is the strength measured at a reference distance d 0 from the anchor i, d i,t is the Euclidian distance between the anchor i and the node at time t and α is the path loss exponent. In practice, the RSSI of a signal could be modified due to the reflection, the diffraction or the scattering of the signal. Moreover the values of ξ 0 and α may vary from an anchor to the other. This may lead to inaccurate distances estimates. For this reason, the proposed method uses connectivity information, instead of using distances estimates. In other words, received strength values are only used to be compared to a threshold ξ r , corresponding to the sensing range r of the sensors. If ξ i,t ≥ ξ r , the anchor i is assumed to be within the sensing range of the node at time t. Otherwise, the anchor i is assumed to be too far and its information is not used. Connectivity measurements are then one-bit information generated as follows,

z i,t = 1 if ξ i,t ≥ ξ r 0 otherwise , i ∈ {1, ..., N a }. (5)
Let I t be the set of indices of all anchors having z i,t = 1 at time t. The anchors denoted in I t are assumed to be within the sensing range of the node, and thus they are located at distances from the node less than r,

∀i ∈ I t , d i,t ≤ r. (6)
One could obtain the expression of the distance d i,t between the anchor i and the considered node in PCS by considering the triangle having the origin O , a i,t and u t as vertices, as shown in Fig. 1. The lengths of the sides O u t and O a i,t are given by ρ t and ρ i,t respectively, whereas d i,t is the length of the side u t a i,t in the triangle. d i,t could then be computed using the generalized Pythagorean theorem as follows,

d 2 i,t = ρ t 2 + ρ i,t 2 -2ρ t ρ i,t cos θ t -θ i,t , (7)
with |θ t -θ i,t | being the angle measured at the vertex O in the considered triangle. The observation constraints are thus rewritten as follows,

ρ t 2 + ρ i,t 2 -2ρ t ρ i,t cos θ t -θ i,t ≤ r 2 , i ∈ I t . (8)
This leads to a set of disks, called observation disks, having the detected anchors as centers and the communication range r as radii.

The localization problem is then defined in the considered PCS by both the mobility and the observation models as follows,

ρ t 2 + ρ t-1 2 -2ρ t ρ t-1 cos θ t -θ t-1 ≤ D 2 m , ρ t 2 + ρ i,t 2 -2ρ t ρ i,t cos θ t -θ i,t ≤ r 2 , i ∈ I t , ρ t ∈ [0, ρ m ], θ t ∈ [0, π 2 ], (9)
where ρ i,t , θ i,t , i ∈ I t , r, D m and ρ m have known values, and ρ t-1 and θ t-1 are estimated at time t-1. Graphically, the problem at a given time t consists of overlapping the mobility disk, centered on the previous position and having D m as radius, with a set of observation disks, having r as radii and the detected anchors as centers. An example of such a problem with three detected anchors is shown in Fig. 2. The solution of the problem is given by the overlapping area of all disks, as shown in dark gray in the plot.

B. Problem statement in several PCSs

In the previous section, the localization problem is defined in a specific PCS, denoted by PCS , having its origin O at the low-left corner of Ω. In PCS , all distances to the origin ρ are included in [0, ρ m ] and all angles θ are included in [0, π 2 ]. One could define the localization problem in an infinite number of PCSs over the surveillance area Ω. However, only three PCSs other than PCS could be set with ρ ∈ [0, ρ m ] and θ ∈ [0, π 2 ]. These PCSs, denoted by PCS , PCS and PCS , have their respective origins O , O and O at the low-right, the up-right and the up-left corners of Ω respectively. Also, they are rotated anticlockwise respectively by angles of π 2 , π and 3π 2 with respect to PCS . An illustration of these PCSs is given in Fig. 3. According to this definition, the bounds over the distances to origins and the angles remain unchanged. Observation and mobility constraints are also valid, which leads to the following formulation of the localization problem in any of these PCSs,

u t u t-1 a 1,t a 2,t a 3,t O Ω ρ t θ t r r r Dm Fig. 2. An example of a localization problem in PCS at time t. PCS u t θ = 0 θ = π 2 a i,t O ρ i,t ρ t θ i,t θ t PCS u t θ = 0 θ = π 2 a i,t O ρ i,t ρ t θ i,t θ t PCS u t θ = 0 θ = π 2 a i,t O ρ i,t ρ t θ i,t θ t PCS u t θ = 0 θ = π 2 a i,t O Ω Ω Ω Ω ρ i,t ρ t θ i,t θ t
ρ ⋄ t 2 + ρ ⋄ t-1 2 -2ρ ⋄ t ρ ⋄ t-1 cos θ ⋄ t -θ ⋄ t-1 ≤ D 2 m , ρ ⋄ t 2 + ρ ⋄ i,t 2 -2ρ ⋄ t ρ ⋄ i,t cos θ ⋄ t -θ ⋄ i,t ≤ r 2 , i ∈ I t , ρ ⋄ t ∈ [0, ρ m ], θ ⋄ t ∈ [0, π 2 ], (10)
where the superindex ⋄ ∈ { , , , }, r and D m are known, ρ ⋄ i,t and θ ⋄ i,t are anchors coordinates expressed differently in the given PCSs and ρ ⋄ t-1 and θ ⋄ t-1 are estimated at time t -1, each in its corresponding PCS. Assume that anchors coordinates are given in the first PCS, denoted by PCS . One could deduce geometrically the expressions of these

ρ i,t = ρ i,t 2 + u 2 0 -2ρ i,t u 0 cos θ i,t , θ i,t = arctan u0-ρ i,t cos(θ i,t ) ρ i,t sin(θ i,t ) , ρ i,t = ρ i,t 2 + 2u 2 0 -2ρ i,t u 0 cos θ i,t + sin θ i,t , θ i,t = arctan u0-ρ i,t sin(θ i,t ) u0-ρ i,t cos(θ i,t ) , ρ i,t = ρ i,t 2 + u 2 0 -2ρ i,t u 0 sin θ i,t , θ i,t = arctan ρ i,t cos(θ i,t ) u0-ρ i,t sin(θ i,t )
.

(11) Once the localization problem is set in PCS , PCS , PCS and PCS , one could choose a single PCS among these four PCSs to solve the problem. Another way to handle the problem consists of performing computation in all the defined PCSs and then select the best resulting solution, as it is shown in the following section. A glossary of all terms is given in TABLE I.

III. POLAR-INTERVAL LOCALIZATION ALGORITHM

Solving the localization problem consists of estimating the coordinates of the considered node at each time step, given the constraints of [START_REF] Tharmarasa | Decentralized sensor selection for large-scale multisensor-multitarget tracking[END_REF]. One could choose a single PCS, for instance PCS , to compute the solution. In a different manner, this paper proposes to handle the problem in all defined PCSs and then choosing the best solution between the resulting ones. Computation is thus performed in each PCS of PCS , PCS , PCS and PCS , leading to four estimates of the position of the node at each time-step. In the following, the proposed method, based on interval analysis [START_REF] Jaulin | Applied interval analysis[END_REF], is first described, and then the localization algorithm is presented.

A. Description of the method

The solution of the problem is proposed using interval analysis [START_REF] Jaulin | Applied interval analysis[END_REF]. Instead of computing an exact position estimate at each time step, the method consists of performing an outer estimation of the solution. In other words, it aims at bounding the coordinates of the node, in the way to cover all possible solutions of the problem. Consider that computation is performed in a specific PCS, denoted by PCS ⋄ with ⋄ ∈ { , , , }. The solution is then given by a two-dimensional interval, denoted by [u

[u ⋄ t ] = [ρ ⋄ t ] × [θ ⋄ t ] = [ρ ⋄ t , ρ ⋄ t ] × [θ ⋄ t , θ ⋄ t ], (12)
where

ρ ⋄ t = inf ([ρ ⋄ t ]) and ρ ⋄ t = sup ([ρ ⋄ t ]
) denote respectively the lower and the higher endpoints of the real interval

[ρ ⋄ t ] and θ ⋄ t = inf ([θ ⋄ t ]) and θ ⋄ t = sup ([θ ⋄ t ]
) denote respectively the lower and the higher endpoints of the real interval

[θ ⋄ t ].
Having the localization problem of (10) at time t, solving the problem consists of finding the minimal polar box [u ⋄ t ] including all possible solutions. Starting with an initial box [u] 0 , the proposed method aims at minimizing the widths

ρ ⋄ t -ρ ⋄ t and θ ⋄ t -θ ⋄ t of [ρ ⋄ t ] and [θ ⋄ t ]
respectively according to the constraints of [START_REF] Tharmarasa | Decentralized sensor selection for large-scale multisensor-multitarget tracking[END_REF]. The initial polar box [u] 0 could be defined by [0,

ρ m ]×[0, π 2 ], since ρ ⋄ t ∈ [0, ρ m ] and θ ⋄ t ∈ [0, π 2 ]
as shown in Section II.

Graphically, the solution box is a partial ring, having the origin O ⋄ of the considered PCS as center and ρ ⋄ t and ρ ⋄ t as inner and outer radii respectively. Moreover, it is defined between the lines starting at the origin and having θ ⋄ t and θ ⋄ t as angles. The best solution box corresponding to the problem of Fig. 2 is given in thick black line in Fig. 4. Here the illustration is shown in PCS and the previous estimate is assumed to be exact. In fact, according to the proposed method, the previous estimate is a polar box denoted by [u t-1 ]. Propagating [u t-1 ] using the mobility disk leads to a larger domain, defined by the union of all the disks obtained by the propagation of each point of [u t-1 ] with the mobility model. An illustration of this domain, called mobility domain, is given in Fig. 5. Actually, instead of using the mobility disk as in Fig. 4 and since the previous solution is a box, one should overlap the mobility domain with the observation disks to define the solution polar box in the considered PCS. An illustration of the solution box of PCS is given in black in Fig. 6. It is clear that this box is larger than the one illustrated in Fig. 4. This way the resulting box is guaranteed to contain the exact solution of the problem since it takes previous and actual incertitude into consideration.

Fig. 7 illustrates the polar boxes obtained in PCS , PCS , PCS and PCS . Here the previous estimates are assumed to be exact for simplicity of illustration. In fact, for each PCS of superindex ⋄ ∈ { , , , }, the corresponding previous polar box [u ⋄ t-1 ] should be employed to compute the actual estimate

[u ⋄ t ].
Once the four estimates are obtained, areas of these estimates are computed as follows,

Ar ([u ⋄ t ]) = θ ⋄ t -θ ⋄ t ρ ⋄ t 2 -ρ ⋄ t 2 2 , (13)
where

[u ⋄ t ] = [ρ ⋄ t , ρ ⋄ t ] × [θ ⋄ t , θ ⋄ t ]
. The values of these areas represent the incertitude obtained while computing the minimal encloser of the solution in the corresponding PCS. For this reason, and since all computed boxes include the solution area of the problem, the polar box having the minimal area is selected to be the solution box at the considered time-step. The three remaining boxes are also kept in the memory to be used with the mobility model in the following time-step.

[u t ] u t-1 a 1,t a 2,t a 3,t O Ω ρ t ρ t θ t θ t Dm r

B. Proposed algorithm

Consider the problem of (10) defined in one PCS, denoted by PCS ⋄ , ⋄ ∈ { , , , }. In order to compute the solution box [u ⋄ t ] at a given time t, one should set all available constraints on ρ ⋄ t and θ ⋄ t according to [START_REF] Tharmarasa | Decentralized sensor selection for large-scale multisensor-multitarget tracking[END_REF]. The first general constraints to be set are given by the dimensions of the surveillance area,

0 ≤ ρ ⋄ t ≤ ρ m and 0 ≤ θ ⋄ t ≤ π 2 . ( 14
)
Moreover, since the cosine of an angle is always less than 1 and ρ ⋄ t and ρ ⋄ i,t are always positive for i ∈ I t , then

ρ ⋄ t -ρ ⋄ i,t 2 ≤ ρ ⋄ t 2 + ρ ⋄ i,t 2 -2ρ ⋄ t ρ ⋄ i,t cos θ ⋄ t -θ ⋄ i,t .
Hence, each observation constraint of [START_REF] Tharmarasa | Decentralized sensor selection for large-scale multisensor-multitarget tracking[END_REF], for i ∈ I t , leads to the following constraint,

ρ ⋄ t -ρ ⋄ i,t 2 ≤ r 2 ⇔ ρ ⋄ i,t -r ≤ ρ ⋄ t ≤ ρ ⋄ i,t + r. ( 15
)
[u t ] [u t-1 ] a 1,t a 2,t a 3,t O Ω ρ t ρ t θ t θ t Dm r
Fig. 6. Polar box obtained in PCS for the problem of Fig. 2 taking the previous polar box into consideration. In the same manner, the mobility constraint leads to the following,

Ω Ω Ω Ω PCS [u t ] O PCS [u t ] O PCS [u t ] O PCS [u t ] O t -1 t -1 t -1 t -1
ρ ⋄ t-1 -D m ≤ ρ ⋄ t ≤ ρ ⋄ t-1 + D m , (16)
where the exact value of ρ ⋄ t-1 is unknown, but according to the previous description,

ρ ⋄ t-1 ∈ [ρ ⋄ t-1 ] with [ρ ⋄ t-1 ] = [ρ ⋄ t-1 , ρ ⋄ t-1 ]. Then, ρ ⋄ t-1 -D m ≤ ρ ⋄ t ≤ ρ ⋄ t-1 + D m . ( 17
)
Combining all constraints over ρ ⋄ t leads then to the following,

max 0, ρ ⋄ t-1 -D m , max i∈It ρ ⋄ i,t -r ≤ ρ ⋄ t ≤ min ρ m , ρ ⋄ t-1 + D m , min i∈It ρ ⋄ i,t + r . (18)
On the other hand,

ρ ⋄ t 2 + ρ ⋄ i,t 2 -2ρ ⋄ t ρ ⋄ i,t cos θ ⋄ t -θ ⋄ i,t = ρ ⋄ t -ρ ⋄ i,t cos θ ⋄ t -θ ⋄ i,t 2 + ρ ⋄ i,t 2 sin 2 θ ⋄ t -θ ⋄ i,t . (19)
Then each observation constraint of (10) leads inevitably to the following,

ρ ⋄ i,t 2 sin 2 θ ⋄ t -θ ⋄ i,t ≤ r 2 ⇔ - r ρ ⋄ i,t ≤ sin θ ⋄ t -θ ⋄ i,t ≤ r ρ ⋄ i,t . ( 20
) If r ρ ⋄ i,t ≤ 1, for i ∈ I t , and since θ ⋄ t -θ ⋄ i,t ∈ [-π 2 , π 2 ]
where the sine function is monotonically increasing, the previous constraint leads to

θ ⋄ i,t -arcsin * r ρ ⋄ i,t ≤ θ ⋄ t ≤ θ ⋄ i,t + arcsin * r ρ ⋄ i,t , (21)
where arcsin * (x) = arcsin (min (1, x)). In the same manner, the mobility constraint leads to the following,

θ ⋄ t-1 -arcsin * D m ρ ⋄ t-1 ≤ θ ⋄ t ≤ θ ⋄ t-1 + arcsin * D m ρ ⋄ t-1 . ( 22
) Since ρ ⋄ t-1 ∈ [ρ ⋄ t-1 ] and θ ⋄ t-1 ∈ [θ ⋄ t-1 ] with [ρ ⋄ t-1 ] = [ρ ⋄ t-1 , ρ ⋄ t-1 ] and [θ ⋄ t-1 ] = [θ ⋄ t-1 , θ ⋄ t-1 ]
, and since the arcsine function is monotonically increasing on the considered domain, then

θ ⋄ t-1 -arcsin * D m ρ ⋄ t-1 ≤ θ ⋄ t ≤ θ ⋄ t-1 + arcsin * D m ρ ⋄ t-1 . (23
) Combining all constraints over θ ⋄ t leads then to the following,

max 0, θ ⋄ t-1 -arcsin * Dm ρ ⋄ t-1 , max i∈I t θ ⋄ i,t -arcsin * r ρ ⋄ i,t ≤ θ ⋄ t ≤ min π 2 , θ ⋄ t-1 + arcsin * Dm ρ ⋄ t-1 , min i∈I t θ ⋄ i,t + arcsin * r ρ ⋄ i,t .
(24) One could then define the first solution box by

[u ⋄ t ] (1) = [ρ ⋄ t (1) , ρ ⋄ t (1) ] × [θ ⋄ t (1) , θ ⋄ t (1) ], (25)
with

ρ ⋄ t (1) = max 0, ρ ⋄ t-1 -Dm, max i∈I t ρ ⋄ i,t -r , ρ ⋄ t (1) = min ρm, ρ ⋄ t-1 + Dm, min i∈I t ρ ⋄ i,t + r , θ ⋄ t (1) = max 0, θ ⋄ t-1 - * Dm ρ ⋄ t-1 , max i∈I t θ ⋄ i,t -arcsin * r ρ ⋄ i,t , θ ⋄ t (1) = min π 2 , θ ⋄ t-1 + arcsin * Dm ρ ⋄ t-1 , min i∈I t θ ⋄ i,t + arcsin * r ρ ⋄ i,t .
(26) Fig. 8 shows in thick black line the first solution box [u ⋄ t ] (1) , obtained for the localization problem of Fig. 2, with ⋄ = . It is obvious that this box is not minimal. It is indeed larger than the solution box that should be obtained, as illustrated in Fig. 6. For this reason, one could set more constraints on ρ ⋄ t and θ ⋄ t , with bounds being functions of θ ⋄ t and ρ ⋄ t respectively. Once [u ⋄ t ] (1) is computed, these constraints could be used to contract it at maximal leading to the final solution box. On the first hand, observation constraints of ( 10) and equation ( 19) lead to the following constraints with i ∈ I t ,

[u t ] (1) [u t-1 ] a 1 a 2 a 3 O Ω ρm ρ t-1 -Dm ρ 1,t -r ρ 2,t -r ρ 3,t -r ρ t-1 + Dm ρ 1,t + r ρ 2,t + r ρ 3,t + r θ t-1 -arcsin * Dm ρ t-1 θ 1,t -arcsin * r ρ 1,t θ
ρ ⋄ i (t) cos θ ⋄ t -θ ⋄ i,t -r 2 -ρ ⋄ i,t 2 sin 2 θ ⋄ t -θ ⋄ i,t ≤ ρ ⋄ t ≤ ρ ⋄ i,t cos θ ⋄ t -θ ⋄ i,t + r 2 -ρ ⋄ i,t 2 sin 2 θ ⋄ t -θ ⋄ i,t . (27)
In the same manner, the mobility constraint leads to the following,

ρ ⋄ t-1 cos θ ⋄ t -θ ⋄ t-1 -D 2 m -ρ ⋄ t-1 2 sin 2 θ ⋄ t -θ ⋄ t-1 ≤ ρ ⋄ t ≤ ρ ⋄ t-1 cos θ ⋄ t -θ ⋄ t-1 + D 2 m -ρ ⋄ t-1 2 sin 2 θ ⋄ t -θ ⋄ t-1 , (28)
with

ρ ⋄ t-1 ∈ [ρ ⋄ t-1 ] and θ ⋄ t-1 ∈ [θ ⋄ t-1
]. On the second hand, observation constraints of (10) could be reformulated as follows,

cos θ ⋄ t -θ ⋄ i,t ≥ ρ ⋄ t 2 + ρ ⋄ i,t 2 -r 2 2ρ ⋄ t ρ ⋄ i,t . (29)
Then,

θ ⋄ i,t -arccos ρ ⋄ t 2 +ρ ⋄ i,t 2 -r 2 2ρ ⋄ t ρ ⋄ i,t ≤ θ ⋄ t ≤ θ ⋄ i,t + arccos ρ ⋄ t 2 +ρ ⋄ i,t 2 -r 2 2ρ ⋄ t ρ ⋄ i,t . (30)
In the same manner, the mobility constraint leads to the following,

θ ⋄ t-1 -arccos   ρ ⋄ t 2 +ρ ⋄ t-1 2 -D 2 m 2ρ ⋄ t ρ ⋄ t-1   ≤ θ ⋄ t ≤ θ ⋄ t-1 + arccos   ρ ⋄ t 2 +ρ ⋄ t-1 2 -D 2 m 2ρ ⋄ t ρ ⋄ t-1   , (31)
with

ρ ⋄ t-1 ∈ [ρ ⋄ t-1 ] and θ ⋄ t-1 ∈ [θ ⋄ t-1 ].

One could use bounds of θ ⋄

t with ( 27) and (28) to compute new bounds of ρ ⋄ t and bounds of ρ ⋄ t with (30) and (31) to compute new bounds of θ ⋄ t , which might contract [u ⋄ t ] (1) . Practically, in order to contract the polar box [u ⋄ t ] (1) , the constraints of ( 27), (28), (30) and (31) are iterated in the interval framework using the forward-backward contractor [START_REF] Jaulin | Applied interval analysis[END_REF], [START_REF] Mourad | Anchor-based localization via interval analysis for mobile ad-hoc sensor networks[END_REF]. This contractor iterates all constraints, while using interval notations, without any prior order until no contraction is possible. The polar box obtained at time t using this contractor would be at best the box illustrated in Fig. 6, given the problem of Fig. 2. It is worth noting that the computation is performed separately in the four PCSs, in order to select at last the resulting box having the minimal area. The proposed algorithm at a considered time t is given in appendix. If an exact estimate is needed, this estimate would correspond to the barycenter of the selected polar box at each time-step.

IV. SIMULATIONS

This section illustrates the performances of the proposed method, denoted PIL4 for Polar Interval-based Localization using 4 PCSs. For this aim, a single mobile node moving over 100s in a 100m × 100m square area is considered. The maximal velocity of the node is set to its exact value, equal to 3.17m/s, leading to a maximal distance of 3.17m with a localization period of 1s. Anchors are assumed to be static either randomly or uniformly deployed over the surveillance area. The sensing area of the sensors is assumed to be circular with a 10m sensing range. It is worth noting that connectivity measurements at a given time step are generated by computing the distances between the mobile node and all anchors at this step and then setting to one the measurement corresponding to the anchor having a distance to the node less than 10m. In the following, all simulations are performed on an Intel(R) Core(TM) i5-2520M CPU (2.5GHz, 4.00GB RAM) using MATLAB 7.10.0.499.

A. Comparison of PIL4 to a one-PCS-based method

This paragraph illustrates the effectiveness of the proposed method, compared to another version of the method, denoted PIL1, using only one PCS which is PCS . Here anchors are first randomly deployed with a density set to 0.015 anchor per m 2 , leading to 150 anchors randomly deployed over the square area, as shown in Fig. 9. With a 10m circular sensing range, the node detects at average 5 anchors at each time step, with the number of detected anchors going from 1 to 11. Fig. 9 shows the polar boxes obtained using PIL4 in straight line and PIL1 in dashed line. The plot also shows the estimated exact positions with PIL4 given at the centers of the estimated boxes. Let the estimation errors be the distances separating the centers of the estimated boxes from the exact positions. Fig. 10 shows the ratios of the areas of the boxes obtained with PIL4 over those of PIL1 in the top plot and the ratios of estimation errors PIL4 over PIL1 in the bottom plot. As illustrated, the estimates of PIL4 are more accurate than those of PIL1. 150 anchors are then randomly deployed in 5 different manners. Both methods are then performed using these distributions and results are averaged over the 5 attempts in the following. TABLE II shows the average boxes areas, the average estimation errors and the average computation times per time-step obtained with PIL4 and PIL1. According to these results, using four PCSs leads to a decrease of around 22% of the boxes areas and 23% of the estimation errors. In other words, performing computation in more than one PCS leads to more accurate estimates with less incertitude, at the cost of the computation time that increases but remains negligible.

PIL4 is then compared to PIL1 with a different anchors configuration. Here anchors are assumed to be uniformly

B. Effectiveness of the mobility model

In this section, the proposed method is compared to another version of the method using only the observation constraints, denoted PIL4o. 150 anchors are first randomly deployed in the surveillance area in 5 different manners. TABLE IV shows the computation results obtained with both methods averaged over the 5 attempts. As expected, the use of the mobility model increases the accuracy of the estimation. 100 uniformly deployed anchors are then considered. TABLE V shows the results obtained with both methods. Here also using the mobility model reduces the incertitude of estimation. Indeed, with the mobility model, the mobility domain is overlapped with the intersection area of observation disks, which reduces the resulting box as shown for instance in Fig. 11. While PIL4o yields the smallest polar box covering the intersection of the observation disks, PIL4 yields a smaller box reduced by the mobility domain. It is worth noting that the mobility box illustrated in gray in the plot corresponds to the polar box covering the exact mobility domain.

C. Effectiveness of the contraction phase

In this section, the proposed method PIL4 is compared to its relaxed version, denoted PIL4r, yielding the first solution box [u ⋄ t ] (1) . In other words, PIL4r is a simplified version of PIL4 not performing the contraction phase that iterates the combined constraints in the forward-backward contractor. The anchors are first assumed to be randomly deployed with a density of 0.015 anchor per m 2 . The distribution of anchors is generated 5 times and all results are averaged over the 5 cases. TABLE VI shows the average boxes areas, the average estimation errors and the average computation costs per timestep obtained with PIL4 and PIL4r. The use of the contraction phase in PIL4 leads thus to a slightly higher accuracy, at the cost of a high increase of the computation time. 100 uniformly deployed anchors are then considered. The results of both methods are illustrated in TABLE VII. As expected, PIL4 yields at least as much accuracy as PIL4r, with higher computation times due to the iteration of more constraints. In the following, the relaxed version PIL4r of the method is considered, since it leads to almost the same accuracy as PIL4 with less computation time.

D. Impact of the anchors density

This section studies the impact of the anchors density on the estimation results of the proposed method PIL4r. The number of anchors is varied from 25 to 225 and the distribution of anchors is assumed to be uniform. With the increase of the anchors density in the network, the average number of detected anchors per time step varies from 0.78 to 7.1. In other words, the number of constraints to be considered in the algorithm is increased, leading to an increase of the estimation accuracy as well as an increase of the computation time. Fig. 12 illustrates the average number of detected anchors, the estimation error, the boxes areas and the computation time as a function of the number of anchors uniformly deployed in the network. As expected, the more we have anchors in the network, the more estimation is accurate and the more the computation time increases.

E. Comparison to another interval-based method

In this section, the proposed method is compared to an interval-based method using cartesian coordinates, as shown in [START_REF] Mourad | Anchor-based localization via interval analysis for mobile ad-hoc sensor networks[END_REF]. Using both observation and mobility models, this method, denoted by CIL, yields rectangular boxes estimates covering all possible solutions. 100 uniformly deployed anchors are first considered. Fig. 13 shows the boxes obtained using the PIL4r and the CIL methods, whereas TABLE VIII illustrates the computation results. It is obvious that with a uniformly deployed distribution of anchors, PIL4r yields much higher accuracy than CIL. 150 anchors are then randomly deployed

Fig. 1 .

 1 Fig. 1. Polar coordinates in PCS .

Fig. 3 .

 3 Fig. 3. Polar coordinates in four different PCSs.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Polar box obtained in PCS for the problem of Fig. 2 with an exact previous estimate.

Fig. 7 .

 7 Fig. 7. Solutions of the problem in the four predefined PCSs.

Fig. 9 .

 9 Fig. 9. Estimated boxes obtained using PIL4 and PIL1.

Fig. 10 .

 10 Fig. 10. Ratios of boxes areas PIL4/PIL1 in the top plot. Ratios of estimation errors PIL4/PIL1 in the bottom plot.

Fig. 11 .

 11 Fig. 11. Estimated boxes obtained with PIL4 and PIL4o at time-step 70.

Fig. 14 .

 14 Fig. 14. Estimated boxes with PIL4r and particles obtained using MCL.

TABLE I GLOSSARY

 I OF TERMS.Polar coordinates vector of the anchor i in PCS ⋄ at time t

	Terms	Definitions
	⋄	Designation of the considered PCS, ⋄ ∈ { , , , }
	ρ ⋄ t	1 st polar coordinate of the node in PCS ⋄ at time t
	θ ⋄ t	2 nd polar coordinate of the node in PCS ⋄ at time t
	u ⋄ t	Polar coordinates vector of the node in PCS ⋄ at time t
	ρ ⋄ i,t	1 st polar coordinate of the anchor i in PCS ⋄ at time t
	θ ⋄ i,t	2 nd polar coordinate of the anchor i in PCS ⋄ at time t
	a ⋄ i,t
	Dm	Maximal distance traveled by the node
	r	Communication range
	ρm	Length of the diagonal of the surveillance area
	u0	Length of a side of the surveillance area
	coordinates in the remaining PCSs as follows,

  ⋄ t ]. Also called polar box, [u ⋄ t ] is defined by the cartesian product of two real intervals [ρ ⋄ t ] and [θ ⋄ t ], defined over the polar coordinates ρ ⋄ t and θ ⋄ t respectively,

TABLE II COMPARISON

 II OF PIL4 TO PIL1 WITH RANDOMLY DEPLOYED ANCHORS.

	Methods	Boxes areas	Estimation errors	Computation times
	PIL4	64.44m 2	1.99m	0.0064s
	PIL1	82.95m 2	2.59m	0.0032s

TABLE III COMPARISON

 III OF PIL4 TO PIL1 WITH UNIFORMLY DEPLOYED ANCHORS.

	Methods	Boxes areas	Estimation errors	Computation times
	PIL4	59.77m 2	2.01m	0.0037s
	PIL1	71.02m 2	2.21m	0.0016s
	deployed with a density of 0.01 anchor per m 2 , leading to
	100 anchors regularly spaced over the surveillance area as
	shown in Fig. 13. With this configuration, the node detects
	at average 3.13 anchors per time-step with the number of
	detected anchors varying between 2 and 4 anchors. TABLE III
	shows the average boxes areas, the average estimation errors
	and the average computation times per time-step obtained with
	PIL4 and PIL1.

According to these results, PIL4 remains more accurate than PIL1 but more time consuming. It is worth noting that the computation times here are less than the ones obtained above since less anchors are detected leading to less observation constraints at each time-step.

TABLE IV COMPARISON

 IV OF PIL4 TO PIL4O WITH RANDOMLY DEPLOYED ANCHORS.

	Methods	Boxes areas	Estimation errors	Computation times
	PIL4	64.44m 2	1.99m	0.0064s
	PIL4o	85.99m 2	2.67m	0.0067s

TABLE V COMPARISON

 V OF PIL4 TO PIL4O WITH UNIFORMLY DEPLOYED ANCHORS.

	Methods	Boxes areas	Estimation errors	Computation times
	PIL4	59.77m 2	2.01m	0.0037s
	PIL4o	74.85m 2	2.26m	0.0038s

TABLE VI COMPARISON

 VI OF PIL4 TO PIL4R WITH RANDOMLY DEPLOYED ANCHORS.

	Methods	Boxes areas	Estimation errors	Computation times
	PIL4	64.44m 2	1.99m	0.0064s
	PIL4r	66.08m 2	2.11m	0.00088s

TABLE VII COMPARISON

 VII OF PIL4 TO PIL4R WITH UNIFORMLY DEPLOYED ANCHORS.

	Methods	Boxes areas	Estimation errors	Computation times
	PIL4	59.77m 2	2.01m	0.0037s
	PIL4r	60.17m 2	2.02m	0.00055s

TABLE X COMPARISON

 X OF PIL4R TO MCL WITH UNIFORMLY DEPLOYED ANCHORS.

	Methods	Estimation errors	Computation times
	PIL4r	2.02m	0.00055s
	MCL (N = 50)	2.77m	0.021s
	MCL (N = 20)	2.87m	0.019s

TABLE XI COMPARISON

 XI OF PIL4R TO MCL WITH RANDOMLY DEPLOYED ANCHORS.

	Methods	Estimation errors	Computation times
	PIL4r	2.11m	0.00088s
	MCL (N = 50)	2.65m	0.034s
	MCL (N = 20)	2.67m	0.031s

over the surveillance area, in 5 different manners. TABLE IX shows the computation results averaged over the 5 cases. As shown in the results table, PIL4r still provide more accurate estimation but here the difference between both results is minor. Indeed, the accuracy of both methods is tightly related to the shape of the solution area, which depends on the distribution of the anchors in the network.

F. Comparison to a Monte-Carlo-based method

This section compares the proposed method to a Monte-Carlo-based method [START_REF] Hu | Localization for mobile sensor networks[END_REF]. This method, called MCL, yields at each time step a fixed number N of positions, called particles, in the way to best cover the solution area. The estimated positions using MCL correspond to the centers of the computed particles. N is first set to 50. Fig. 14 shows the estimated particles using MCL as well as the estimated boxes using PIL4r with 100 uniformly deployed anchors. Let the estimation errors of MCL be the distances between the centers of the particles and the actual positions. TABLE X shows the average estimation errors and the average computation times per time-step obtained with PIL4r and MCL while using 50 and 20 particles. According to these results, the MCL method leads at some time-steps to particles covering precisely the solution area but it remains less accurate than PILr4 with more computation time even with less particles. It is also more consuming in terms of memory resources. One could reduce

V. CONCLUSION

This paper proposes an original technique for sensors localization in wireless sensor networks. The proposed approach is an anchor-based method using connectivity measurements, as well as the mobility of the nodes. Based on interval analysis, the solution is provided using polar intervals. An outer approximation of the solution is thus performed in four different polar coordinate systems. Instead of exact estimates, the proposed method yields partial rings including for sure all the possible solutions of the problem. Simulation results corroborate the efficiency of the proposed method compared to other methods based on Monte-Carlo or cartesian intervals. Future works will consider the problem under imperfect circumstances, where a robust extension of the proposed method would be proposed. Realistic simulations based on real sensors with real RSSI measurements would also be considered.


Member, IEEE Paul Honeine
email: paul.honeine@utt.fr

C Richard
email: cedric.richard@unice.fr.

Geometric unmixing of large hyperspectral images: a barycentric coordinate approach

In hyperspectral imaging, spectral unmixing is one of the most challenging and fundamental problems. It consists of breaking down the spectrum of a mixed pixel into a set of pure spectra, called endmembers, and their contributions, called abundances. Many endmember extraction techniques have been proposed in the literature, based on either a statistical or a geometrical formulation. However, most if not all these techniques for estimating abundances use a least squares solution. In this paper, we show that abundances can be estimated using a geometric formulation. To this end, we express abundances with the barycentric coordinates in the simplex defined by endmembers. We propose to write them in terms of a ratio of volumes or a ratio of distances, which are quantities that are often computed to identify endmembers. This property allows us to easily incorporate abundance estimation within conventional endmember extraction techniques, without incurring additional computational complexity. We use this key property with various endmember extraction techniques, such as N-Findr, vertex component analysis, simplex growing algorithm, and iterated constrained endmembers. The relevance of the method is illustrated with experimental results on real hyperspectral images.

I. INTRODUCTION

The emergence of hyperspectral imaging sensors in recent years has brought new opportunities and challenges in image analysis. Hyperspectral images are cubes of data, measuring spectral composition within a spatial view. As opposed to conventional three-channel color systems, the spectral information provides in-depth analysis of the composition of objects in the image scene. Since the early 90's, hyperspectral imaging has been adopted as an airborne technique for military and environmental remote sensing, and more recently to identify chemical species on the surface and the atmosphere of Mars [START_REF] Moussaoui | On the decomposition of mars hyperspectral data by ICA and Bayesian positive source separation[END_REF]. Hyperspectral imaging provides abundant information about ground composition, thanks to an improved resolution in the spatial dimensions, with a high sensitivity in the spectral range, up-to hundreds of contiguous bands covering visible and short-wavelength infrared spectral range.

Even with high spatial resolution, a single pixel may consist of different materials. Spectra in an hyperspectral image can be seen as a mixing of some spectral signatures of pure physical components over the ground [START_REF] Keshava | Spectral unmixing[END_REF]. Despite the potential of nonlinear unmixing techniques [START_REF] Guilfoyle | A quantitative and comparative analysis of linear and nonlinear spectral mixture models using radial basis function neural networks[END_REF], the linear mixing model has a suitable physical interpretation widely accepted within the remote sensing community. In this case, pure spectral signatures, called endmembers, are linearly combined with some abundance fractions. The key challenges for spectral unmixing are twofold: find the collection of endmembers; and estimate their abundances for each pixel in the hyperspectral image.

To identify the endmembers, early techniques were based on a pre-defined library of laboratory minerals, unrelated to the image under investigation [START_REF] Clark | USGS digital spectral library splib06[END_REF]. Because these spectra are rarely acquired under the same conditions, this brings several difficulties in atmospheric correction and variations in the sensor intrinsic parameters. Current approaches avoid such problems by inducing the endmembers from the scene, which falls into the class of blind source separation problems (see [START_REF] Lawton | Self-modeling curve resolution[END_REF] for an early study, and [START_REF] Chan | A convex analysis framework for blind separation of non-negative sources[END_REF] for a more recent one). Therefore, a natural tool to blindly unmixing is the independent component analysis (ICA). Concerning hyperspectral data, the assumption that sources are statistically independent is violated since, due to physical constraints in the data acquisition process, the sum of the abundance fractions is constant, implying dependence among abundances [START_REF] Nascimento | Does independent component analysis play a role in unmixing hyperspectral data ?[END_REF]. Such limitations have contributed to the development of new unmixing techniques, roughly divided into two classes, geometric approaches and statistical ones (see [START_REF] Parente | Survey of geometric and statistical unmixing algorithms for hyperspectral images[END_REF] for a recent survey). While the latter requires probabilistic foundations with some prior knowledge, geometric approach exploits the theory of convex set. In this paper, we give particular attention to the geometric formulation. Prior to the endmember extraction, most of these techniques are never used on the original spectra. To correct the redundancy in the spectral dimension as opposed to the low number of endmembers, a dimensionality reduction transformation is often applied, such as the principal component analysis (PCA) or the minimum noise fraction (MNF).

From the linear mixing model, mixed pixel spectra are linear mixtures of endmembers. As a consequence, they can be viewed as the vertices of a simplex enclosing all mixed spectra. Geometric techniques exploit this insight, by seeking the vertices of a simplex englobing all the image data. One of the most widely used automatic endmember extraction method is Winter's N-Findr algorithm [START_REF] Winter | N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF]. It starts with a random set of candidate endmembers (see [START_REF] Zortea | A quantitative and comparative analysis of different implementations of N-FINDR: A a quantitative and comparative analysis of different implementations of nfindr: A fast endmember extraction algorithm[END_REF], [START_REF] Wang | Construction of fast and robust N-FINDR algorithm[END_REF] for several implementations of N-Findr). Iteratively, one at a time, each pixel spectrum is considered as a candidate to replace each endmember. It is accepted as a new endmember if the volume of the simplex of endmembers increases. The N-Findr is a fast endmember extraction algorithm, successfully applied to many hyperspectral images. Sequential algorithms can provide further reduction in computational complexity, by growing the simplex gradually, vertex by vertex. The Simplex Growing Algorithm (SGA) is a sequential version of the N-Findr [START_REF] Chang | A new growing method for simplex-based endmember extraction algorithm[END_REF]. Starting with one vertex, it finds a simplex with maximum volume by sequentially adding a new (optimal) vertex at a time. The Vertex Component Analysis (VCA) algorithm is another sequential algorithm [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF]. It exploits that orthogonal subspace projection of a simplex is also a simplex, and thus vertices of the latter are vertices of the initial one. The VCA gradually selects vertices by performing orthogonal subspace projections, sequentially starting from a one-dimensional subspace. The Iterated Constrained Endmembers (ICE) algorithm fits an enclosing simplex to the data cloud while penalizing its volume [START_REF] Berman | ICE: a statistical approach to identifying endmembers in hyperspectral images[END_REF]. To this end, an alternating optimization scheme is used by minimizing, on the one hand a quadratic cost function for fitness, and on the other hand the distance between endmembers for penalization. A sparsity-promoting version of ICE is presented in [START_REF] Zare | Hyperspectral band selection and endmember detection using sparsity promoting priors[END_REF].

Two physical constraints are generally imposed on the linear mixing model, enforcing constraints on the fractional abundances for each pixel: (1) the sum-to-one (or equality) constraint, indicating that the endmembers contribution must add up to 100%, and (2) the non-negativity constraint, since negative contributions are physically unrealistic. The most straightforward approach to estimate abundances is by an unconstrained least-squares minimization, with the abundances of a given pixel determined by a matrix inversion [START_REF] Strang | Introduction to Linear Algebra[END_REF]. For each pixel, the computational complexity is cubic in the number of endmembers. The sum-to-one constraint can be imposed to the least-squares solution, using for instance Lagrangian multipliers. The non-negativity constraint is more difficult to address, since it does not have a closed-form solution. Many iterative approaches have been proposed to this end, including non-negative matrix factorization techniques [START_REF] Miron | Joint processing of the parallel and crossed polarized Raman spectra and uniqueness in blind nonnegative source separation[END_REF], [START_REF] Jia | Constrained nonnegative matrix factorization for hyperspectral unmixing[END_REF], as well as general iterative techniques [START_REF] Lawson | Solving Least Squares Problems[END_REF], [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF]. A fully constrained solution requires more advanced optimization methods [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF], [START_REF] Cardenas | Fully constrained linear spectral unmixing: Analytic solution using fuzzy sets[END_REF]. Solutions to these problems involve computational extensive algorithms, not appropriate for large number of hyperspectral pixels (for implementation issues, see for instance [START_REF] Robila | Considerations on parallelizing nonnegative matrix factorization for hyperspectral data unmixing[END_REF], [START_REF] Schmidt | Implementation strategies for hyperspectral unmixing using bayesian source separation[END_REF]. It is worth noting that, if the endmembers were properly identified, i.e., the vertices of a simplex englobing all the data, both sum-to-one and non-negativity constraints would be naturally fulfilled. For this reason, extraction techniques often apply directly the unconstrained least-squares method [START_REF] Craig | Minimum-volume transforms for remotely sensed data[END_REF], [START_REF] Winter | N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF], [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF].

In this paper, we give a geometric formulation to solve the abundance estimation problem. To this end, we redefine the fractional abundances as barycentric coordinates (or areal coordinates) in the reduced dimensional space. We derive two expressions of these coordinates, either as a volume-ratio of simplexes, or as a ratio of distances. It is amazing to observe that these volumes and distances are often computed by well-known endmember extraction methods in order to identify the endmembers. This is the case of N-Findr and SGA algorithms which compute volumes, while VCA and ICE estimate distances. By taking advantage of this situation, we show that one can easily incorporate abundance estimation within conventional endmember extraction techniques, without essentially incurring additional computational cost.

The remainder of this paper is organized as follows. In Section II, we present barycentric coordinates for the estimation of fractional abundances. In Section III, we demonstrate that our approach is natural to many endmember extraction techniques, and we outline its application to N-Findr, SGA, VCA, and ICE. The effectiveness of the method is illustrated in Section IV, by comparing it to the fully constrained least-squares solution. Finally, Section V summarizes this work, and describes ongoing efforts. But before, we give a mathematical formulation of the linear mixing/unmixing model.

THE LINEAR (UN-)MIXING MODEL

Given a hyperspectral image, let x ℓ be the spectrum of the ℓ-th pixel, each pixel being indexed using a single index. The linear mixing model takes the form

x ℓ = n i=1 α i x ω i , (1)
where {x ω 1 , x ω 2 , . . . , x ωn } denotes the collection of the n endmembers. In matrix form, the linear mixing model is given by

x ℓ = Xα, (2)
where

X = [x ω 1 x ω 2 • • • x ωn ], and α = [α 1 α 2 • • • α n ]
⊤ the vector of abundances to be determined. We suppose for now that the endmembers have been identified, using any off-the-shelf endmember extraction technique. The fractional abundance α i represents the contribution in x ℓ of the pure material defined by x ω i . In order to adopt a physical interpretation of this mixing model, the abundances must satisfy two constraints:

• the sum-to-one constraint (also called the equality constraint), with n i=1 α i = 1;

• the non-negativity constraint, with α 1 , α 2 , . . . , α n ≥ 0.

These constraints can be written in matrix form, as 1 ⊤ α = 1 and α ≥ 0, where 1 is the unit column-vector of n entries and the inequality is taken component-wise.

Ignoring these constraints, the optimal solution of the problem (2), with optimality in the least-squares sense, is given by the normal equation

α LS = (X ⊤ X) -1 X ⊤ x ℓ . (3)
The equality constraint is the simplest constraint to enforce, and can be incorporated by using Lagrangian multipliers, resulting into the optimal solution

α eqLS = α LS - 1 1 ⊤ (X ⊤ X) -1 1 (X ⊤ X) -1 1(1 ⊤ α LS -1).
(4) Enforcing the non-negativity constraint is not as easy to address as the equality constraint, since no closed-form solution exists. In the least-squares sense, the problem consists of minimizing

f (α) = 1 2 x ℓ -Xα 2 , while maintaining α ≥ 0. Let ∇f (α) = X ⊤ (Xα -x ℓ ) denotes the gradient of f (α).
Then, the Karush-Kuhn-Tucker optimality conditions [START_REF] Nocedal | Numerical Optimization[END_REF] for this problem are α ≥ 0, ∇f (α) ≥ 0, and ∇f (α) ⊤ α = 0. This can be solved using iterative techniques, such as an alternating strategy applied to both expressions α ⊤ X ⊤ (Xαx ℓ ) = 0, and α ≥ 0. (5)

The most commonly used algorithm to compute the solution is the Lawson and Hanson's algorithm [START_REF] Lawson | Solving Least Squares Problems[END_REF], which uses the active set scheme. Other iterative algorithms include multiplicative iterative strategies [START_REF] Lantéri | A general method to devise maximum-likelihood signal restoration multiplicative algorithms with non-negativity constraints[END_REF]. The problem of estimating the abundances using the linear model subject to both equality and non-negativity constraints, the so-called fully constrained solution, requires advanced iterative optimization methods [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF]. All these solutions are optimal, in the least-squares sense, and can be considered even when the number of bands (dimension) is greater than the number of endmembers.

In practice, the spectral dimension is reduced to n -1 in order to have the mixing model (1) satisfied. Such representations allow to extract the endmembers. We shall now give a direct geometric formalism to estimate the abundances in the same low-dimensional representation. To this end, we solve the equality-constrained optimization problem geometrically, and give a straighforward interpretation to the violation of the non-negativity constraint.

II. BARYCENTRIC COORDINATES FOR ABUNDANCE ESTIMATION

Suppose for now that the endmembers have been extracted, using for instance any of the above mentioned techniques, such as N-Findr or SGA algorithms. We seek to estimate the abundances, α 1 , α 2 , . . . , α n for each spectrum in the image. Given any x ℓ , we consider the linear mixing model with the equality constraint, namely

x ℓ = n i=1 α i x ω i , (6)
subject to

n i=1 α i = 1. (7)
Since both are equality expressions, we combine them into a single matrix expression, with the augmented linear system

1 1 • • • 1 x ω 1 x ω 2 • • • x ωn α = 1 x ℓ . (8)
It is worth noting that this is a linear system with as many equations as unknowns, i.e., n. Linear algebra provides many elegant methods to solve such systems [START_REF] Strang | Introduction to Linear Algebra[END_REF].

We propose to solve this problem using the Cramer's rule [16, page 259]. The consequences of this choice will be clear throughout this paper. Cramer's rule states that the solution of a linear system can be expressed in terms of the determinants of the system matrix, and of matrices obtained from it with one column substituted by the right-handside vector. Applied to the above augmented linear system [START_REF] Parente | Survey of geometric and statistical unmixing algorithms for hyperspectral images[END_REF], the solution can be written as

α 1 = det 1 1 • • • 1 x ℓ x ω 2 • • • x ωn det 1 1 • • • 1 x ω 1 x ω 2 • • • x ωn , Fig. 1.
The volume of a tetrahedron (three-dimensional simplex) is given by the third of the area of the base times the height from the base to the apex. Since both (blue and red) tetrahedrons share the same base (gray color), the ratio of their volumes are proportional to the ratio of their heights (dashed lines).

α 2 = det 1 1 • • • 1 x ω 1 x ℓ • • • x ωn det 1 1 • • • 1 x ω 1 x ω 2 • • • x ωn , . . . α n = det 1 1 • • • 1 x ω 1 x ω 2 • • • x ℓ det 1 1 • • • 1 x ω 1 x ω 2 • • • x ωn , ( 9
)
where det is the determinant operator. The amazing thing about Cramer's rule is that one only needs to compute determinants in order to solve a linear system, as opposed to conventional matrix inversion techniques. In fact, both methods are equivalent. Moreover, applied to the proposed augmented linear system, we get an intuitive geometric interpretation, natural in the context of simplexes with vertices defined by hyperspectral data.

A geometric formulation

Now, we recall the expression of the volume of a simplex. Let X = {x ω 1 , x ω 2 , . . . , x ωn } be the set of linearly independent vectors in a (n -1)dimensional space, thus defining a non-degenerated simplex in that space. This assumption is true in the case of endmembers. We define the oriented volume x ω 1

x ω 1 x ω 3 x ω 2 x ℓ x ω 1 x ω 3 x ω 2 x ℓ
x ω 3 x ω 2 x ℓ δ(x ℓ ) δ(x ω1 ) x ω 1 x ω 3 x ω 2 x ℓ δ(x ℓ ) δ(x ω1 )
Fig. 3. Illustration in 2D of α1 computation using the heights. In this case, α1 = δ(x ℓ )/δ(xω 1 ). In the left figure, x ℓ and xω 1 are on the same side of the line defined by xω 2 and xω 3 , and α1 is non-negative. In the right figure, α1 is negative since x ℓ and xω 1 lie on both sides.

of the simplex with vertices X, as follows

V X = 1 (n -1)! det 1 1 • • • 1 x ω 1 x ω 2 • • • x ωn .
(10) This expression gives a signed value, either positive or negative depending on the order of the simplex vertices, i.e., the sequence x ω 1 , x ω 2 , . . . , x ωn . The consequences of the use of the oriented (signed) volume, rather than the use of the conventional positive-valued expression |V X |, will be justified.

Let \ denotes the set difference operator, with X\{x ω i }∪{x ℓ } the set defined by the set X where x ω i has been removed and x ℓ added. From Cramer's rule [START_REF] Winter | N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF], each coefficient α i can be written as a ratio of two oriented volumes of simplexes, on the one hand the simplex with vertices X given by the endmembers, and on the other hand the same simplex with x ω i replaced by x ℓ . In a compact form, we have

α i = V X\{xω i }∪{x ℓ } V X , (11)
for all i = 1, 2, . . . , n. As shown below, these volumes are often computed by the endmember extraction techniques, as in the case of the N-Findr algorithm for instance. Therefore, one can estimate abundances with essentially no computational cost.

A distance-based formulation

All the abundances can be estimated using [START_REF] Wang | Construction of fast and robust N-FINDR algorithm[END_REF], or equivalently from [START_REF] Winter | N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[END_REF]. This illustrates the fact that, for each abundance, the two considered simplexes have n -1 identical vertices and only one

x ω 1 x ω 3 x ω 2 x ℓ
x ω 1

x ω 3

x ω 2 x ℓ Fig. 4. Illustration of endmember extraction incremental algorithms in a 2-dimensional Euclidean space. From two previously extracted endmembers, xω 1 and xω 2 , the third endmember is determined (left) using SGA by comparing simplex volumes, or (right) using VCA by comparing distances to the subspace spanned by previous endmembers. The abundance α3 is given by (left) the ratio of volumes or (right) the ratio of distances. different vertex, x ω i , which has been substituted in the numerator by x ℓ . Since these n -1 vertices provide a common base for both simplexes, we can further simplify the expression (11) (see Figure 1 for an illustration in three-dimensions). In fact, remember that the volume of a simplex of n vertices can be given in terms of the content (length, area, volume, ...) of a base of n -1 vertices and the distance of the remaining vertex to the subspace including the base. Namely, we have

V X = 1 n-1 δ(x ω i ) V X\{xω i }
, where δ(x ω i ) is the (signed) distance between the vertex x ω i and the subspace spanned by the other vertices in X. By using this notation, we get

α i = δ(x ℓ ) δ(x ω i ) , (12)
for all i = 1, 2, . . . , n. This formulation is more adapted to incremental endmember extraction techniques, as opposed to the volume ratio in [START_REF] Wang | Construction of fast and robust N-FINDR algorithm[END_REF]. This is mainly due to the fact that this expression does not require explicit relations between previously selected endmembers. Only distances to the subspace spanned by the latter are considered. This results in low computational complexity algorithms, as illustrated for instance with the VCA.

On the non-negativity constraint

As defined above, the oriented volume of a simplex can be either positive or negative, depending on the order of the sequence defined by its vertices. In order to study the sign of any coefficient α i , we need to compare the volumes in [START_REF] Wang | Construction of fast and robust N-FINDR algorithm[END_REF], or equivalently the corresponding determinants. The only difference resides in one column, the i-th column, either x ω i or x ℓ . Therefore, one only needs to compare the orientation of these simplexes. If both simplexes have the same orientation, the coefficient is positive; otherwise, when they have opposite orientations, it is negative.

To illustrate this point, we consider a twodimensional example in Figure 2. Consider the first coefficient, given by α 1 = V X\{xω 1 }∪{x ℓ } /V X . Therefore, we compare the simplexes X = {x ω 1 , x ω 2 , x ω 3 } and {x ℓ , x ω 2 , x ω 3 }. Figure 2 (left) illustrates the case where both simplexes have the same orientation, counterclockwise in the sequences x ω 1 x ω 2 x ω 3 and x ℓ x ω 2 x ω 3 . In this case, α 1 is nonnegative. Figure 2 (right) illustrates the case of opposite orientations, with a clockwise orientation in the sequence x ℓ x ω 2 x ω 3 as opposed to the counterclockwise in x ω 1 x ω 2 x ω 3 . This results in α 1 < 0. Applying the same argument for the other coefficients, we find that α 2 and α 3 are non-negative in both cases (left and right figures).

Consider now signed-distances in equation ( 12). As illustrated in Figure 3, x ω 2 x ω 3 divides IR 2 into two regions. Coefficient α 1 is positive if both x ℓ and x ω 1 lie on the same side (left figure). Otherwise, it is negative (right figure).

These results obtained for all the coefficients are in accordance with the definition of a convex hull, saying that x ℓ is inside the simplex if and only if all the α i 's are non-negative. Moreover, we can derive a simple inclusion test as follows: x ℓ is inside the simplex if and only if all the determinants in ( 9) have the same sign (either positive or negative); otherwise, if it is outside, there exists at least one determinant with opposite sign. This provides an elegant inclusion test, where only determinant evaluation is required, a step carried out by the endmember extraction stage as derived in the next section.

It is worth noting that the non-negativity constraints are not imposed in the proposed method. However, the violation of these constraints, i.e., the existence of at least one negative α i , means that x ℓ is outside the simplex. Such a result implies that both equality and non-negativity constraints cannot be satisfied by the resulting linear combination of endmembers. This illustrates the limitations of geometric approaches for endmember extraction (described hereafter) for properly identifying the endmembers 1 .

III. ABUNDANCE ESTIMATION USING STATE-OF-THE-ART ENDMEMBER EXTRACTION TECHNIQUES

Many endmember extraction techniques determine endmembers as vertices of a simplex englobing all the spectra. In this section, we show how one can easily take advantage of some of these techniques, in order to estimate the abundances using 1 Such a result is also related to the uniqueness of the unmixing solution (or the endmember identifiability condition as given in [27, Theorem 1]). In fact, let x ℓ 1 , x ℓ 2 , . . . , x ℓn be the spectra corresponding to some arbitrary pixels. Then we have from ( 8)

1 • • • 1 x ℓ 1 • • • x ℓn = 1 • • • 1 xω 1 • • • xω n [α ℓ 1 • • • α ℓn ],
where the equality also holds for the volumes as defined in [START_REF] Zortea | A quantitative and comparative analysis of different implementations of N-FINDR: A a quantitative and comparative analysis of different implementations of nfindr: A fast endmember extraction algorithm[END_REF], namely

|V {x ℓ 1 ,...,x ℓn } | = |V {xω 1 ,...,xω n } | | det[α ℓ 1 • • • α ℓn ]|.
The determinant in this expression can be upper-bounded using Hadamard's inequality, with

| det[α ℓ 1 • • • α ℓn ]| ≤ n j=1 α ℓ j ,
where the inequality becomes equality when the matrix is a permutation matrix. If the sum-to-one and the non-negativity constraints are satisfied, then this upper bound equals 1. Therefore, the volume of the simplex defined by any arbitrary set cannot be greater than the one defined by the endmembers.

barycentric coordinates. To keep the presentation as simple as possible, we restrict ourselves to classical, well-known, techniques.

A. N-Findr algorithm

The N-Findr algorithm seeks the simplex with the largest volume, in an iterative manner by visiting each pixel to inflate the simplex. In a pre-processing stage, a dimensionality reduction method, such as PCA or MNF, is applied to get data into a (n -1)-dimensional space, n being the number of endmembers. At initialization, an initial set of n endmember candidates is selected randomly 2 . Let X = {x ω 1 , x ω 2 , . . . , x ωn } be the set of corresponding spectra.

The following process is iterated for each pixel, x ℓ being its spectrum. One at a time, each endmember is replaced by the spectrum under investigation, and the oriented3 volume of the resulting simplex is evaluated. We get n new volumes, V X\{xω i }∪{x ℓ } , as well as the initial volume V X . Comparing these volumes, we get the following decision rule:

• if the tested volumes are less than the initial one, in absolute values, namely max i |V X\{xω i }∪{x ℓ } | < |V X |, then the initial set of endmembers remains unchanged; • otherwise, in order to inflate the volume, some entry of the endmember set is replaced by x ℓ . The outgoing spectrum x ω i is identified as follows: ω i = arg max |V X\{xω i }∪{x ℓ } |. It turns out that the volumes already computed for endmember extraction can be used to estimate the abundances using expression [START_REF] Wang | Construction of fast and robust N-FINDR algorithm[END_REF], namely by defining the contribution of x ω i in x ℓ by

α i = V X\{xω i }∪{x ℓ } V X .
As given in this expression, the contribution of any endmember requires only an arithmetic division for each pixel. Therefore, estimating all the abundances for all the pixels can be done with a simple arithmetic division, applied nN times.

B. Simplex growing algorithm

The simplex growing algorithm (SGA) improves the N-Findr in two main aspects: (1) at initialization, the starting set is not randomly selected, and (2) the simplex is growing one endmember at a time, until all endmembers are extracted.

At initialization, a single dimension is considered, using for instance PCA, where a simplex is a nontrivial line segment. The first endmember x ω 1 is determined by maximizing |V {xω 1 ,x 0 } |, where x 0 is randomly selected. Iteratively, the endmembers are extracted one at a time. At step i for 2 ≤ i ≤ n, the dimensionality being reduced to i -1, the i-th endmember is determined by the index ω i as follows

ω i = arg max ℓ |V X∪{x ℓ } |,
where X corresponds to the set of endmembers previously extracted. This is illustrated in Figure 4 (right).

In the same spirit as the N-Findr, the volumes computed by the SGA allow to estimate the abundances, as given by expression [START_REF] Wang | Construction of fast and robust N-FINDR algorithm[END_REF]. The difference resides in that, at a given step i, the computed volumes give the contributions of ω i in all x ℓ in the (i -1)-dimensional space, and not in the final (n -1)-dimensional space defined by the n-order mixing model.

C. Vertex component analysis

The vertex component analysis (VCA) exploits the fact that an affine transformation of a simplex is also a simplex, in particular by projection operation. The incremental scheme projects the data onto a direction orthogonal to the subspace spanned by previously selected endmembers. The extreme of this projection specifies the new endmember. The algorithm iterates until it reaches the desired number of endmembers.

Roughly speaking, the VCA considers the distance of candidate data to the subspace spanned by the endmembers already determined. The new endmember corresponds to the farest data, as illustrated in Figure 4 (left). It turns out that these distances are the same as those introduced in (12), and consequently allow us to estimate abundances, as follows

α 3 = δ(x ℓ ) δ(x ω 3 ) .
Once again, only one arithmetic division is required to estimate the contribution of any endmember for each pixel. Moreover, the VCA selection criteria, i.e., x ω 3 = arg max x ℓ δ(x ℓ ), gives a natural framework for both equality and non-negativity constraints. If there exists some x ℓ such that δ(x ω 3 ) < δ(x ℓ ), then from the above expression α 3 > 1, and consequently the constraints are violated.

D. Connections to the ICE method

The iterated constrained endmembers (ICE) method seeks the simplex englobing the data, by minimizing the residuals in a least-squares sense. Essentially, the algorithm minimizes, on the one hand, the mean squared error

1 N N ℓ=1 x ℓ - n i=1 α i x ω i 2 ,
and on the other, the sum of squared distances between all endmembers, namely n i,j=1 i =j

x ω i -x ω j 2 .
The former enlarges the simplex englobing the data, while the latter constrains its volume. ICE uses an alternating minimization scheme to extract endmembers (from a quadratic objective function) and estimate their contributions (solving a quadratic programming problem).

The ICE method considers (squared) distances between all couples of data, endmembers and nonendmembers. Using distances between vertices of a simplex, we can easily compute its volume thanks to the Cayley-Menger determinants [START_REF] Blumenthal | Theory and Applications of Distance Geometry[END_REF]. The square volume of a simplex in an (n-1)-dimensional space is given by

V 2 X = c n-1 det         0 1 1 1 • • • 1 1 0 δ 2 1,2 δ 2 1,3 • • • δ 2 1,n 1 δ 2 2,1 0 δ 2 2,3 • • • δ 2 2,n 1 δ 2 3,1 δ 2 3,2 0 • • • δ 2 3,n . . . . . . . . . . . . . . . . . . 1 δ 2 n,1 δ 2 n,2 δ 2 n,3 • • • 0        
, where δ i,j denotes the distance between two vertices and c n = (-1) n+1 /2 n (n!) 2 . While the Cayley-Menger expression computes the determinant of a (n + 1)-by-(n + 1) matrix,

x ω 1 x ω 2 x ωn x ω k x ℓ 1 1-α k n i=1 i =k α i x ω i α k 1 -α k δ(x ℓ ) δ(xω k )
Fig. 5. By applying the intercept theorem for the two dashed lines (blue color) in the (red) triangle, we get

α k = δ(x ℓ )/δ(xω k ).
its particular form results into lower computational cost. To show this, consider the case of n -1 = 2 dimensions. Then the Cayley-Menger determinant leads to the Heron's formula of triangles, with

V 2 X = 1 16 (δ 1,2 + δ 1,3 + δ 2,3 )(δ 1,2 + δ 1,3 -δ 2,3 )(δ 1,2 - δ 1,3 + δ 2,3 )(-δ 1,2 + δ 1,3 + δ 2,3
), where only four arithmetic multiplications are required.

E. Connection to subspace projection techniques

The Orthogonal Subspace Projection (OSP) method is a general technique for detection, discrimination and classification in hyperspectral imaging [START_REF] Chang | Orthogonal subspace projection (OSP) revisited: a comprehensive study and analysis[END_REF]. Applied for endmember extraction, the abundances are often given using a least-squares estimate [START_REF]Further results on relationship between spectral unmixing and subspace projection[END_REF]. It is described by writing (1) as

x ℓ = α k x ω k + n i=1 i =k α i x ω i ,
for any k = 1, 2, . . . , n. To extract x ω k , an operator (matrix) is designed to annihilate the second term in the right-hand-side. Let U = X\{x ω k }, then one can consider the projection operator

I -U (U ⊤ U ) -1 U ⊤ . ( 13
)
This operator projects the data into the null space of U , as shown in [START_REF] Luo | Null space spectral projection algorithm for endmember extraction in hyperspectral remote sensing image (in chinese)[END_REF]. In [START_REF] Luo | Null subspace analysis for spectral unmixing in hyperspectral remote sensing[END_REF], the abundances are estimated using a geometric point of view, in the same spirit as the barycentric approach proposed in this paper. The connection between both techniques is illustrated in Figure 5 where, by applying the intercept theorem, one gets relation [START_REF] Chang | A new growing method for simplex-based endmember extraction algorithm[END_REF], namely

α k = δ(x ℓ ) δ(x ω k ) .
It is worth noting that our method does not require any matrix inversion, as opposed to subspace projection techniques as illustrated in expression [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF].

IV. EXPERIMENTATIONS In this section, we study the relevance of the proposed abundance estimation scheme. To this end, we compare it to the results obtained with the least-squares solutions, with the unconstrained (3), the equality-constrained (4), the non-negativity constraint [START_REF] Lawton | Self-modeling curve resolution[END_REF], and the fully constrained [START_REF] Heinz | Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[END_REF] solutions. We emphasize on the fact that all these leastsquares solutions are computed from the spectral data, without any dimension reduction, as opposed to our approach which exploits the PCA (or MNF) representation, already computed for endmember extraction. As we work in a lower dimension, one should expect to have lower performance.

The study of the optimality of the dimensionality reduction technique and of the endmember extraction method are beyond the scope of this paper.

A. Real hyperspectral image

The studied hypersepectral image is the scene of the Cuprite mining district in western Nevada, USA, captured by the NASA's Airborne Visible/Infrared Imaging Spectrometer (AVIRIS). Over the past years, this site has been extensively studied by the remote sensing community [START_REF] Goetz | Mineralogical mapping in the cuprite mining district, nevada[END_REF], providing highaccuracy ground truth map [START_REF] Swayze | Ground-truthing aviris mineral mapping at cuprite, nevada[END_REF], [START_REF] Swayze | The hydrothermal and structural history of the cuprite mining district, southwestern nevada: An integrated geological and geophysical approach[END_REF]. The AVIRIS sensor collects 224 contiguous spectral bands covering wavelengths ranging from 0.38 to 2.5 µm, with a spectral resolution of approximately 0.01 µm. Atmospheric correction is applied (removing the bands 1-2, 104-113, 148-167, and 221-224), which yields a total of 188 bands. Previously studied in [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF], [START_REF] Miao | Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[END_REF] for endmember extraction, the subimage consists of 250 lines of 191 pixels each, as shown in Figure 6.

In the experiments, we used N-Findr to extract the endmembers, set here to n = 3 for illustration purpose (two-dimensional representation and RGB map presented below). This is not the optimal number of endmembers for this image, with an estimated virtual dimensionality equals to 12 as studied in [START_REF] Nascimento | Vertex component analysis: A fast algorithm to unmix hyperspectral data[END_REF], or 9 for almost the same scene in [START_REF] Miao | Endmember extraction from highly mixed data using minimum volume constrained nonnegative matrix factorization[END_REF]. The identified endmembers are the first three endmembers obtained from the VCA, namely Alunite, Kaolinite and Sphene, a result conform to the ground truth information. The spectra of these endmembers are given in Figure 7.

We applied barycentric coordinates, i.e., using [START_REF] Wang | Construction of fast and robust N-FINDR algorithm[END_REF], to estimate the abundances of these endmembers within the hyperspectral image. These results are illustrated on the map in Figure 8, using a threechannel color composite, with red, green and blue for the first, second and third extracted endmembers. With both the endmember extraction and the abundance estimation operated in the two-dimensional space given by the PCA, we represented data in Figure 10. With (30%) of the data outside the simplex defined by these endmembers, 14 586 data (out of 47 750) violated the non-negativity constraint, with at least one estimated negative α i . The spatial distribution of the pixels violating the non-negativity constraint indicate specific areas where the linear model does not hold (see for instance [START_REF] Moussaoui | On the decomposition of mars hyperspectral data by ICA and Bayesian positive source separation[END_REF], [START_REF] Duran | Spectral unmixing with negative and superunity abundances for subpixel anomaly detection[END_REF]). As illustrated in Figure 9 where the largest negative α i 's are shown (pixels with zero value correspond to pixels with all abundances non-negative), the violation of the linear model clearly exhibits a spatial structure, which could result from regions with undetected endmembers. These negative abundances were set to zero for the RGB map in the threechannel color Figure 8. This can be considered as a drawback of using a dimensionality reduction In order to measure the performance of the abundance estimation technique, we considered the spectral angle. For each data, it is defined between the initial spectrum, x ℓ , and the one computed with the linear mixture using the estimated abundances, x ℓ , with θ(x ℓ , x ℓ ) = cos -1

x ℓ , x ℓ x ℓ x ℓ .

The spectral angle is invariant with respect to the energy, making it a widely used measure of error for hyperspectral data [START_REF] Kruse | The spectral image processing system (SIPS): Interactive visualization and analysis of imaging spectrometer data[END_REF]. The performance of the proposed approach was measured using the spectral angle, and compared with the fully constrained least squares. For this purpose, we considered the spectral angle between the initial spectra and, on the one hand, the ones reconstructed using the fully constrained least squares and, on the other, the barycentric coordinates. Figure 11 shows the histograms of both distributions, demonstrating the small angular error for either techniques, with a mean angle of 0.0866 and 0.0848, respectively. The study of the spectral angle between spectra estimated by each technique showed a very small difference, as illustrated in Figure 12.

V. CONCLUSION AND ON-GOING WORK

In this paper, we presented a new abundance estimation method, based on pure geometry as opposed to conventional matrix inversion techniques. We showed that, once a dimensionality reduction technique such as PCA is applied, it is useless to perform a least squares estimation. Thanks to barycentric coordinates, the abundances are expressed as a ratio of volumes or as a ratio of distances. Since these quantities are often computed for endmembers identification, we showed that this method is inherent to well-known endmember extraction techniques, such as N-Findr, SGA, VCA and ICE. Experimental results showed the relevance of the proposed approach.

The goal of our ongoing work is two-fold. On the one hand, we are working on imposing the non-negativity constraint to barycentric coordinates. This can be applied using simple geometry, in the same spirit of this paper. On the other hand, we are interested in a parallel implementation, using graphics hardware (see for instance [START_REF] Setoain | Parallel morphological endmember extraction using commodity graphics hardware[END_REF]). It turns out that barycentric coordinates have been largely used in computer graphics applications, such as defining an inclusion test in rasterization or linear interpolation for shading. Therefore, a parallel implementation of the proposed method is natural to graphics hardware.

Fig. 2 .

 2 Fig. 2. Illustration of the simplex (here triangle) in a 2-dimensional Euclidean space, to study the the sign of coefficient α1. In the left figure, x ℓ and xω 1 are on the same side of the line defined by xω 2 and xω 3 , thus α1 is non-negative. In the right figure, α1 is negative since they lie on opposite sides.

Fig. 6 .

 6 Fig. 6. The AVIRIS Cuprite Nevada hyperspectral subimage, at wavelength 1.3 µm (band 100).

3 Fig. 7 .

 37 Fig. 7.The three extracted endmembers, representing mainly Alunite, Kaolinite and Sphene.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Three-channel color composite, with red •, green • and blue • for the first, second and third extracted endmembers. The location of each endmember, indicated by x, is coherent with ground truth measurements.

Fig. 10 .

 10 Fig. 10. The Cuprite hyperspectral data (•) in the two-dimensional space defined by the PCA, with endmembers (•) identifying a triangle enclosing 70% of the data.

Fig. 12 .

 12 Fig. 12. Histogram of the spectral angle of the estimated spectra, between the ones obtained from the barycentric coordinates and the fully constrained least-squares techniques.

Fig. 11 .

 11 Fig. 11. Histogram of the spectral angle between the initial spectra and the ones obtained using the abundances estimated with the fully constrained least-squares (left) and the proposed technique (right).

TABLE I NUMBER

 I OF HYPERSPECTRAL DATA FROM THE CUPRITE SCENE (OUT OF 47 750) WITH ESTIMATED ABUNDANCES VIOLATING THE NON-NEGATIVITY AND THE SUM-TO-ONE CONSTRAINTS.

			n
		∃ i : α i < 0	α i = 1
			i=1
	Unconstrained least-squares (3)	43 065 (90%)	47 746 (99%)
	Sum-to-one least-squares (4)	13 997 (29%)	0
	Non-negative least-squares [19]	0	47 746 (99%)
	Fully constrained least-squares [21]	0	0
	Barycentric coordinates [this paper]	14 586 (30%)	0
	technique, here PCA. However, applying the sum-
	to-one constrained least-squares (4) on the hyper-
	spectral data, without any dimensionality reduction,
	yielded up to 29% of the data violating the non-
	negativity constraint. Table I presents a quantitative
	comparison of the constraints violation for different
	abundance estimation techniques.

The initial set influences the resulting endmembers, and therefore the estimated abundances. More efficient selection techniques may also be applied to construct the initial set. See, e.g.,[START_REF] Plaza | An improved N-FINDR algorithm in implementation[END_REF]. We keep the exposition clear with the original N-Findr algorithm.

The original algorithm uses a positive-valued volume, here recasted into a signed one. While this turns out to be very useful, the algorithm remains essentially the same.