Simon Henrot

Charles Soussen

David Brie

Fast Positive Deconvolution of Hyperspectral Images

Keywords: Image restoration, hyperspectral images, regularized least squares, positivity, cross-spectral prior information

In this correspondence, we provide an efficient scheme for performing deconvolution of large hyperspectral images under a positivity constraint, while accounting for spatial and spectral smoothness of the data.

I. INTRODUCTION

In many imaging applications, the acquisition process induces various degradations on the collected image which can prevent accurate post-processing of the data of interest. For instance, the image of a biological sample by a microscope is blurred due to diffraction and corrupted by stochastic noise mainly caused by photon-to-electron conversion [START_REF] Sarder | Deconvolution methods for 3D fluorescence microscopy images[END_REF]. Restoration methods such as deconvolution are thus used whenever post processing requires it.

Dealing with hyperspectral images (that is, a stack of images representing the same scene captured at many wavelengths) adds the challenge of accounting for cross-spectral information within the restoration process [START_REF] Zhang | A hyperspectral image restoration technique[END_REF], [START_REF] Bongard | 3D deconvolution of hyper-spectral astronomical data[END_REF]: the intensity values of a given pixel at two neighboring wavelengths should be close. When the data size is large, advanced deconvolution methods become computationally prohibitive and linear methods such as Tikhonov or Wiener-filtering can be used to process the set of images [START_REF] Galatsanos | Digital restoration of multichannel images[END_REF], [START_REF] Galatsanos | Least squares restoration of multichannel images[END_REF]. Within this framework, the problem is formulated as the minimization of a compound criterion which possesses a closed-form solution. Up to the circular approximation of the observation matrix, this solution can be efficiently computed in the Fourier domain using 2D Fast Fourier Transforms (FFTs) [START_REF] Thiébaut | Introduction to image reconstruction and inverse problems[END_REF].

A drawback of these methods is their inability to restore the original object's high frequency components beyond the PSF bandwidth. This typically results in ringing artifacts around edges which can generate negative-valued pixels in the deconvolved image [START_REF] Shah | Deconvolution algorithms for fluorescence and electron microscopy[END_REF]. When the image is known to be positive (e.g., acquired by a photon counting process), a positivity constraint needs to be applied to yield a physically meaningful estimate. In [START_REF] Giovannelli | Positive deconvolution for superimposed extended source and point sources[END_REF], the authors use the augmented Lagrangian method to solve this constrained optimization problem for 2D images ; another fast 2D restoration method using variable splitting is given in [START_REF] Afonso | Fast image recovery using variable splitting and constrained optimization[END_REF]. Other popular algorithms are the Image Space Reconstruction Algorithm (ISRA) in the case of Gaussian noise and the Richardson-Lucy algorithm for Poisson noise: these methods can be applied to the hyperspectral positive restoration problem but do not lead to fast computations [START_REF] Benvenuto | Nonnegative leastsquares image deblurring: improved gradient projection approaches[END_REF].

To the best of our knowledge, no hyperspectral restoration algorithm accounts for both cross-spectral information and positivity in a fast manner. The main contribution of this paper is to propose an algorithm that incorporates these two priors while remaining computationally efficient.

The remainder of this paper is organized as follows. Section II introduces the observation model and notations used throughout the paper. The deconvolution problem is formulated in section III. In section IV, we detail the procedure we use to solve the problem. Section V presents experimental results on synthetic data and some insight into the influence of the hyperparameters. Finally, we conclude in section VI. This work was supported by the French Agence Nationale de la Recherche, with reference number ANR-09-BLAN-0336-04. The authors are with the Université de Lorraine and the CNRS at the CRAN (Research Center for Automatic Control), UMR 7039, Campus Sciences, BP 70239, 54506 Vandoeuvre cedex, France, e-mail: firstname.lastname@cran.uhp-nancy.fr.

II. OBSERVATION MODEL AND NOTATIONS

The hyperspectral image to be restored X is assumed to be a stack of N matrices {X , = 1, . . . , N } of size M × M (we use square matrices for simplicity and without loss of generality), where subscript refers to the wavelengths or channels. Using lexicographical order, the stack can be reshaped into a collection of M 2 × 1 vectors {x , = 1, . . . , N } or a single M 2 N × 1 vector x with x t = [x t 1 , . . . , x t N ] and we use similar notations for other images.

A. Within-channel observation model

Within a given channel , X generates data Y according to the two-dimensional (2D) observation model (see for instance [START_REF] Bertero | Image deblurring with Poisson data: from cells to galaxies[END_REF])

Y = H * 2D X + N ( 1
)
where H is the M × M convolution kernel (possibly comprising null entries) encoding the -th channel Point-Spread Function (PSF):

H =    H 11 . . . H 1M . . . . . . . . . H M 1 . . . H M M    , (2)
* 2D denotes the discrete 2D convolution product carried out in the image domain, yielding an output of size M × M and N is a noise term assumed to be independent and identically distributed (i.i.d) and Gaussian. Equation (1) rewrites as a linear system:

y = H x + n (3)
where H is a M 2 ×M 2 block-Toeplitz matrix with M ×M Toeplitz blocks, e.g. for M = 3, H reads

H =   H 2 H 1 0 H 3 H 2 H 1 0 H 3 H 2   with H i =   H 2i H 1i 0 H 3i H 2i H 1i 0 H 3i H 2i   .
Working assumption: We impose periodic boundary conditions on PSFs H 's. This approximation allows H to be rewritten as a blockcirculant matrix with circulant blocks [START_REF] Galatsanos | Least squares restoration of multichannel images[END_REF], a structure we will denote by circulant-block-circulant (CBC):

H =   H 2 H 1 H 3 H 3 H 2 H 1 H 1 H 3 H 2   with H i =   H 2i H 1i H 3i H 3i H 2i H 1i H 1i H 3i H 2i   .
A well-known property of CBC matrices is that they can be diagonalized by a two-dimensional discrete Fourier transform (2D-DFT): this property is used in section IV.

B. Hyperspectral observation model

The hyperspectral observation model reads y = Hx+n where the N M 2 ×N M 2 observation matrix H is block-diagonal since imaging in a given channel does not interfere with other channels = :

H =       H1 0 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 HN       . ( 4
)
We further assume that the noise process n for a given channel is independent of the other channels = , and that the noise variance is common to all channels. Using this assumption, the noise process n (accounting for all channels together) obeys an i.i.d. Gaussian distribution.

III. PROBLEM FORMULATION

A. Linear restoration problem

The deconvolution problem can be formulated as the estimation of x given data y and the system response H. In the framework of least square restoration [START_REF] Galatsanos | Least squares restoration of multichannel images[END_REF], it is usually written as the unconstrained minimization of a composite criterion

J (x) = y -Hx 2 + µ Dx 2 + ν Ex 2 (5)
where the first term measures the goodness of fit between the data y and the observation model for a candidate solution x and the second and third term regularize the solution, e.g. by penalizing large differences between spatially and/or spectrally adjacent pixels. Hyperparameters µ and ν are positive and balance the relative strength of the three terms of J (depending on the application and the level of blur, one of the parameters may be zero). D and E are respectively the spatial and spectral regularization matrices with M 2 N columns:

• For each channel , spatial regularization is performed by selecting a convolution kernel D representing some kind of difference operator, e.g. the Laplacian filter:

  0 -1 0 -1 4 -1 0 -1 0   . (6)
One can then construct the corresponding M 2 ×M 2 CBC block D similar to the procedure given in section II-A. The hyperspectral spatial regularization matrix D has a block-diagonal structure:

D =       D1 0 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 DN       . (7)
• Spectral regularization occurs across channels, so matrix E is not block-diagonal. A simple choice to penalize large amplitude differences at adjacent wavelengths is to use a first order derivative filter E01 in the wavelength domain, e.g. E0 = (1, -1), which yields the 1D convolution matrix:

E0 =         -1 1 0 . . . 0 0 -1 . . . . . . . . . . . . . . . . . . . . . 0 0 . . . . . . -1 1         (8)
of size (N -1) × N . The spectral regularization matrix reads

E = E0 ⊗ IM ( 9
)
where IM is the identity matrix of size M ×M , and ⊗ stands for the Kronecker product between matrices, consisting in replacing the 0's and 1's in (8) by 0 and IM , respectively. The method proposed hereafter is not specific to the regularization operators defined in (7 -9); it can be applied to any convolution operators D 's with periodic boundary conditions and to any E0 (without imposing specific boundary conditions). Although matrices D and E are very large, they are mainly introduced for formal reasons. As we will see in Section IV, they are not explicitly stored in practice because the terms Dx 2 and Ex 2 can be efficiently computed using two-dimensional Fourier transforms [START_REF] Henrot | Fast deconvolution of large hyperspectral images[END_REF].

B. Positive restoration problem

The positive deconvolution problem is formulated in the framework of constrained optimization

min x J (x) s.t. x ≥ 0 (10)
This is a classic constrained minimization problem which has no closed-form solution. Constrained optimization methods are usually classified as interior point, active set and exterior penalty methods [START_REF] Nocedal | Numerical Optimization[END_REF]. The quadratic penalty method belongs to the latter class and is often used in computations for its simplicity [13, chapter 17].

The rationale for choosing this technique is that it allows us to use a slightly modified implementation of the unconstrained minimization of J , based on fast computations in the Fourier domain, in order to yield the constrained estimate. This results in an algorithm with a relatively low complexity, which is especially useful when processing large hyperspectral images. The minimization procedure is detailed in section IV.

IV. NUMERICAL OPTIMIZATION

The quadratic penalty method simply consists in replacing the constraints by quadratic penalty terms in the objective function. Handling inequality constraints can be done by introducing M 2 N slack variables si and replacing the inequalities {xi ≥ 0 [START_REF] Benvenuto | Nonnegative leastsquares image deblurring: improved gradient projection approaches[END_REF] is replaced by the surrogate problem:

, i = 1 . . . M 2 N } by {xi -si = 0, si ≥ 0, i = 1 . . . M 2 N }. The original problem
min x,s K(x, s) s.t. s ≥ 0 (11)
where the new objective function reads

K(x, s) = y -Hx 2 + µ Dx 2 + ν Ex 2 + ξ x -s 2 . ( 12
)
When ξ tends to infinity, the entries of x identify with the corresponding entries of s: the constraint x ≥ 0 is asymptotically fulfilled. The surrogate problem is unconstrained w.r.t. x and involves the constraint s ≥ 0. However, the positivy constraint is much easier to handle than in the original problem since the objective function is quadratic and separable w.r.t. s. Minimization of K w.r.t. (x, s) is performed by alternating unconstrained minimization w.r.t. x ( §IV-A), constrained minimization w.r.t. s ( §IV-B) and increasing ξ ( §IV-C). The three steps are done once per iteration until some appropriate stopping criterion is met ( §IV-D).

A. Unconstrained minimization

At a given iteration, the first step consists in solving the unconstrained quadratic minimization problem:

min x K(x, s) (13)
We first rewrite the objective function in terms of its channel components:

K(X, S) = N =1 ( Y -H * 2D X 2 + µ D * 2D X 2 + ν {E0 * 1D X} 2 + ξ X -S 2 ) ( 14
)
where * 1D denotes the discrete 1D convolution product carried out in the spectral (wavelength) domain, i.e., there are as many 1D convolution products as the number of image pixels, and . is the matrix Frobenius norm. Using the 2D-DFT (denoted by the usual hat notation) and Parseval's Theorem, criterion [START_REF] Press | Numerical recipes: The art of scientific computing[END_REF] reads

K ( X, Ŝ) = N =1 ( Ŷ -Ĥ × X 2 + µ D × X 2 + ν {E0 * 1D X} 2 + ξ X -Ŝ 2 ) (15)
where × denotes entry-wise multiplication, . now refers to the matrix Frobenius norm in the spatial frequency domain and * 1D still refers to the 1D convolution in the wavelength domain. To obtain [START_REF] Chouzenoux | Efficiency of liensearch strategies in interior point methods for linearly constrained signal restoration[END_REF], we have exploited that for a given , {E0 * 1D X} is a linear combination of images {X , = 1 . . . N } weighted by the nonzero values of E0. Thus, the 2D-DFT of {E0 * 1D X} is the linear combination of images X 's with identical weights.

To ensure that all images X 's are real-valued, it is necessary to impose that X 's are Hermitian matrices. Thus, it can readily be shown that the problem reduces to the minimization of K over half the 2D spatial frequency plane, e.g. the right-half plane and the whole minimizer is deduced using Hermitian symmetry.

We remark that in [START_REF] Chouzenoux | Efficiency of liensearch strategies in interior point methods for linearly constrained signal restoration[END_REF], each squared norm reads as a sum over the spatial frequencies f of the related matrix elements (e.g., the first term is the sum over f of | Ŷ (f ) -Ĥ (f ) X (f )| 2 ). Because the convolution {E0 * 1D X} is done in the wavelength domain, the objective function K ( X, Ŝ) is separable with respect to each spatial frequency f . To be more specific, the problem may be recast as a collection of M 2 N/2 independent quadratic problems, each corresponding to a point f in the right space frequency half-plane, with N unknowns:

min x f { y f -∆H(f )x f 2 +µ ∆D(f )x f 2 +ν E0x f 2 +ξ x f -s f 2 }
(16) where complex vectors x f , y f and s f and complex diagonal matrices ∆H(f ) and ∆D(f ) index the following quantities

x f { X (f ), = 1 . . . N } y f { Ŷ (f ), = 1 . . . N } s f { Ŝ (f ), = 1 . . . N } ∆H(f ) diag{ Ĥ (f ), = 1 . . . N } ∆D(f ) diag{ D (f ), = 1 . . . N }
and E0 is the 1D convolution matrix induced by kernel E0. For each f in the right half-plane, the solution of ( 16) is then simply given by

x f = T f -1 (∆H(f ) y f + ξs f ) (17)
where stands for the complex conjugate and T f is the real N × N matrix given by

T f = ∆H(f ) ∆H(f ) + µ∆D(f ) ∆D(f ) + νE t 0 E0 + ξIN . ( 18
)
Choosing E0 = (1, -1) leads to (8) which implies that T f is a tridiagonal matrix.

B. Constrained minimization

The second step consists in solving the following separable quadratic minimization problem:

min s≥0 K(x, s) (19)
where the positivity constraint on s is easy to handle since the only term of K in which s is involved is the last term xs 2 . This term is separable with respect to si's and the solution reads

si = max (0, xi). (20)
C. Update of ξ Increasing ξ at each iteration ensures that the algorithm converges to the constrained minimizer of J . A broad range of schemes is available and a common choice is to use a linear rule ξ (k+1) = αξ (k) where α > 1 [13, p. 493].

D. Proposed algorithm

The algorithm procedure can be summed up as follows (τ is a stopping parameter):

1) for = 1 . . . N , compute the 2D-DFTs of data Y and PSFs H (H 's are usually theoretically known or experimentally measured); set ξ (1) = 1 and α > 1 and τ < 0; initialize s = 0; 2) repeat a) for each , compute S from s using 2D-FFTs; b) for each f in the right-half plane, compute { X (f ), = 1 . . . N } using ( 17) and ( 18); c) for each f in the left-half plane, deduce X (f )'s using Hermitian symmetry; d) for each , compute x from X using inverse 2D-FFTs; e) compute s from (20); f) update coefficient ξ (k+1) = αξ (k) . until the stopping criterion is met, e.g. the lowest-valued pixel of X is superior to some tolerance threshold τ .

E. Computational complexity

Using the FFT, the complexity of a 2D-DFT is O(M 2 log 2 (M )). The complexity of each iteration is dominated by FFTs since each N ×N tridiagonal linear inversion reduces to O(N ) using the Crout's LU decomposition algorithm [START_REF] Press | Numerical recipes: The art of scientific computing[END_REF]. Hence, each iteration is done in O(N M 2 log 2 (M )) operations.

Additionally, no expensive storage is required since the largest matrix involved in computation T f is sparse and of size N × N .

V. EXPERIMENTAL RESULTS

We illustrate the performance of the algorithm with experimental results on synthetic data. We use a linear mixing model to generate a multichannel cube of data representing simulated objects against a uniform background, captured at several wavelengths. In this model, the spectrum of one pixel (i, j) is given by

X (i, j) = a1(i, j)u1( ) + a2(i, j)u2( ) (21)
where endmember u1 is the segment [0, 1] sampled on N points and u2 is the vector obtained by flipping u1 in the left -right direction. By setting the corresponding abundances a1 and a2 to an arbitrary M × M spatial pattern, one may synthesize objects that respectively fade in and out when the channel index increases from 1 to N . Blurring is simulated by convolving each channel of the image by a M × M Airy disc computed by the 2D DFT of a disk whose radius linearly depends on the channel . The stack is further corrupted by an additive white Gaussian noise term so as to obtain a SNR of 15 dB. The resulting M ×M ×N images are not meant to be realistic but provide a simple way to evaluate the performance of the algorithm from a computational standpoint.

A. Comparison with other constrained optimization methods

In this section we investigate constrained optimization methods to solve problem [START_REF] Benvenuto | Nonnegative leastsquares image deblurring: improved gradient projection approaches[END_REF]. Nonnegativity constraints are most commonly handled by the log-barrier interior point method [START_REF] Chouzenoux | Efficiency of liensearch strategies in interior point methods for linearly constrained signal restoration[END_REF], [START_REF] Nocedal | Numerical Optimization[END_REF]. We use the MATLAB optimization toolbox to write two different implementations of this technique, using either matrices with the quadprog function and a Fourier domain formulation with fmincon. We also compare these results with the unconstrained estimate [START_REF] Henrot | Fast deconvolution of large hyperspectral images[END_REF]. A pair of hyperparameters is arbitrarily selected and results are given for a 2.4 Ghz Intel Core 2 Duo processor with a RAM of 4 gigabytes.

The unconstrained estimate performs worse in terms of MSE (defined as xtruexest 2 2 ) : see figure 1. We observe that all three constrained optimization methods do converge to the same estimate, yielding very similar objective values and identical reconstruction errors. However, our method strongly outperforms the others in terms of computation time and memory requirements, as shown in table I.

The main reason for the increased computation time is that the log-barrier objective function is not quadratic, hence its gradient is nonlinear and thus minimizing the objective w.r.t. to x must be done in an iterative fashion, e.g. using a quasi-Newton step. When the data size is large, closed-form minimizations at each iteration make the restoration problem tractable. Furthermore, we observe that computing convolution matrices quickly reaches memory limitations.

B. Sensitivity to the choice of hyperparameters

Running the algorithm involves choosing the hyperparameter pair (µ, ν). Because of the tradeoff between spatial and spectral regularization and their interplay with the positivity constraint, this is not a trivial problem and its analytical resolution is still an open problem. Furthermore, the choice strongly depends on the application, e.g. the restored image should minimize the reconstruction error as a preprocessing step to spectral unmixing, but applying more smoothing (hence selecting higher parameters) is often desirable for visual inspection. In this section we select the MSE as the relevant criterion. We attribute a [10 -4 , 10 4 ] range to each hyperparameter and discretize it on a 2D grid with logarithmic scales. For each grid point, we restore the data using the proposed method and we compute the MSE. The results corresponding to the least MSEs are presented in figure 1 for M = 25 and N = 3.

We observe that the best constrained estimate in the MSE sense is obtained for lower hyperparameters, i.e. requires less smoothing than the best unconstrained estimate: see for instance [START_REF] Bardsley | The stabilizing properties of nonnegativity constraints in least-squares image reconstruction[END_REF]. In this simulation, the estimation is not very sensitive to the variation of hyperparameters, the ratio between maximum and minimum MSEs on the grid being 1.3 for the constrained and 1.6 for the unconstrained estimate. What is more, the former has a lower reconstruction error than the latter and presents slightly less restoration artifacts: see figure 1. Notice how more ringing artifacts are located near the edges of the unconstrained estimate as opposed to the more "uniform" distribution of errors for the constrained estimate.

In practical applications we obviously do not have access to the true image and hyperparameters may best be selected by trial and error. However, this process remains tractable in practice since the proposed method remains fast even for large data and numerous hyperparameter candidate values.

C. 3D + λ Fluorescence microscopy simulation

In this section we demonstrate the performance of the algorithm on synthetic data resembling actual fluorescence microscopy images for visual inspection and to get an idea of the processing time of 4D data for a MATLAB implementation. We create 3D abundance maps of physical size 25µm × 25µm × 5µm to mimic the random spatial distribution of bacteria using spheres of 1µm radius, yielding 255 × 255 × 5 data cubes. Pseudo fluorescent spectra modeled as sums of Gaussian signals in the 300 -600nm range are fed into the linear mixing model to create the 255 × 255 × 5 × 30 data stack. We use the scalar diffraction-limited Stokseth model [START_REF] Pankajakshan | Blind Deconvolution for Confocal Laser Scanning Microscopy[END_REF] to generate a widefield microscopy point-spread function with a 63X/1.4 oil immersion objective and a CCD camera pixel size of 0.1µm. The image is blurred and an additive white Gaussian noise term is added to obtain a SNR of 15 dB: in widefield microscopy, the SNR is high enough so that the Poisson distributed photo-counting noise term can be accurately approximated by a Gaussian process [START_REF] Pankajakshan | Blind Deconvolution for Confocal Laser Scanning Microscopy[END_REF]. Results are presented in figure 2. Due to the use of 4D tensors, the running time has increased to 15 min for ten iterations of the restoration algorithm. Deconvolution achieves a 77% decrease of the MSE and the proposed estimate thus has a significantly greater resolution. By using information from neighboring wavelengths, our method can restore severely degraded image planes which cannot be tackled by simpler schemes such as PSF fitting.

VI. CONCLUSION

This paper presents a method to solve a hyperspectral image deconvolution problem where one needs to jointly account for spatial, spectral and non-negativity information. The problem is formulated as the minimization of a regularized least-squares criterion under a positivity constraint. While a closed-form solution to the unconstrained problem is known, we provide an intuitive implementation based on FFTs. We show how to easily incorporate a positivity constraint into the estimate. The main advantages of the method are its ease of implementation since it does not require computation of large matrices and its speed, mainly due to the use of 2D-FFTs and low-complexity tridiagonal inversions. Furthermore, the algorithm can make use of parallel processing on each pixel of the stack of spectra using Graphics Processing Units (GPUs) [START_REF] Plaza | Parallel hyperspectral image and signal processing[END_REF] with no additional implementation issues. The method can also be extended to hyperspectral images with three spatial dimensions in a straightforward manner. Future investigations might include accounting for Poisson noise in accordance with confocal microscopy observation models.

11 Fig. 1 .Fig. 2 .

 1112 Fig. 1. True, degraded and restored 25 × 25 × 3 toy images for a SNR of 15. Each column corresponds to a channel.

TABLE I COMPARISON

 I OF THE COMPUTATION TIMES OF THE DIFFERENT OPTIMIZATION METHODS

Convolution kernel E 0 is written with a subscript to remain coherent with other kernels D .

ACKNOWLEDGMENT

The authors thank Dr. S. Moussaoui and Prof. J. Idier from IRCCyN (UMR CNRS 7597, Nantes, France) for fruitful discussions regarding interior point algorithms.


Stéphane Pomportes

Anthony Busson

Joanna Tomasik
email: joanna.tomasik@supelec.fr

Véronique Vèque

Resource Allocation in Ad Hoc Networks with Two-Hop Interference Resolution

The multi-user medium access mechanism OFDMA has to provide each node with a given amount of radio resources. In this paper we present a new distributed algorithm for the allocation of resource blocks in an OFDMA ad hoc network. We are principally interested in allocating resources fairly because the ad hoc networks which we work on are dedicated to be deployed in the areas of natural or man-made disasters and where the guarantee of connectivity is an important issue. Contrary to the commonly applied approach, we consider a resource allocation on the links under a two hop interference distance. The proposed allocation procedure is coupled with our other algorithm which detects and corrects two hop interferences and which has been revised and improved. The performance of our algorithm is evaluated by simulation for different topologies. We observed that simultaneous allocations in large networks allow a constant convergence time to be kept despite the networks size.

I. INTRODUCTION

The radio medium access in ad hoc networks usually relies on a contention mechanism like CSMA/CA (Carrier Sense Multiple Access with Collision Avoidance) used for example in the standards IEEE 802.11 [START_REF]IEEE Standard -Part 11: Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications[END_REF] or IEEE 802.15.4 [2]. With such a protocol, ensuring a minimum data throughput remains a very difficult task because this mecanism does not guarantee a data access rate or any service. On the other hand, the access mechanisms based on reservation like OFDMA (Orthogonal Frequency-Division Multiple Access) [START_REF] Lawrey | Multiuser OFDM[END_REF] ensures to each node a given amount of radio resources. The WiMAX uses this mechanism but only in the infrastructure mode where the allocation is centralized in Access Points. Nevertheless, the IEEE 802.16j Working Group [4] is interested in extending this technology to the ad hoc networks. In OFDMA, the smallest resource allocation unit is defined by a pair (time, frequency) called Resource Block (RB). The allocation of RBs is still an open problem. In this paper, we present a distributed algorithm for allocating RBs in an OFDMA ad hoc network.

There are two ways to allocate RBs in an OFDMA network. The first one, called the node-oriented [START_REF] Bjorklund | Resource optimization of spatial TDMA in ad hoc radio networks: A column generation approach[END_REF] allocation, assigns the RBs to communication equipment. With this solution, RBs are assigned to each node and used to communicate with all its neighbors. The other approach, called the link-oriented [START_REF] Bjorklund | Resource optimization of spatial TDMA in ad hoc radio networks: A column generation approach[END_REF] allocation, assigns resources to communication links. When a resource can be used for both directions of communication, a link is considered as bidirectional. Another approach consists in dedicating resources to a single direction, and links are considered as unidirectional. The choice of the approach depends essentially on the nature of communications. In the case of broadcasts, for example with a handheld transceiver, the node-oriented allocation is more appropriate. However, for point-to-point communications, such as data transfer or video transmission, the link-oriented allocation is preferable. Both approaches also have different properties. The article [START_REF] Gronkvist | Assignment methods for spatial reuse TDMA[END_REF] has shown that allocating resources on the links has a better spatial reuse of resources than the node-oriented allocation and therefore a better bandwidth utilization.

Authors of [START_REF] Rhee | Drand: Distributed randomized tdma scheduling for wireless ad hoc networks[END_REF], [START_REF] Zhu | A Five-Phase Reservation Protocol (FPRP) for Mobile Ad Hoc Networks[END_REF] and [START_REF] Kanzaki | Dynamic TDMA Slot Assignment in Ad Hoc Networks[END_REF] provide a solution for the nodeoriented resources allocation by considering an interference area equal to the communication area. Concerning the link oriented allocation, the authors of [START_REF] Salonidis | Distributed dynamic scheduling for endto-end rate guarantees in wireless ad hoc networks[END_REF] offer a solution for bidirectional links. The unidirectional approach is discussed in papers [START_REF] Djukic | Link scheduling for minimum delay in spatial re-use TDMA[END_REF], [START_REF]Delay aware link scheduling for multi-hop TDMA wireless networks[END_REF] and [START_REF] Gandham | Link scheduling in wireless sensor networks: Distributed edge-coloring revisited[END_REF]. In all these works, an interference area of the same size as the transmission area is considered. In our opinion this assumption is not realistic and leads to numerous interfering allocations as the power of a signal may be too weak to be received, but still strong enough to interfere with another signal. We consider an interference area twice as large as the communication area.

To the best of our knowledge, the RBs allocation on unidirectional links, taking into consideration interference at two hops, has not been addressed beforehand. The contribution of this paper is as follows. It proposes an allocation algorithm assigning RBs to unidirectional links taking into consideration a two hop interference distance for ad hoc networks. An RB cannot be allocated if there is another transmitter for this RB less than two hops away from the receiver.

This assumption makes the allocation algorithm more complex because each node, performing an allocation, impacts its two hop neighborhood whereas it can only communicate with its direct neighbors. In Section 2, we present the model used and assumptions made in this paper. In Section 3, we present our algorithms. In Section 4, we study the performance of our algorithm by simulations on chain, grid and random topologies.

II. MODEL AND ASSUMPTIONS

A. Network modeling

The network is modeled by a directed graph G = (V, E) where V is a set of nodes and E is a set of arcs. This graph is called the connectivity graph in the rest of the paper. Each node is identified by a unique identifier (ID) and models a communication unit. Each arc, (t, r), represents a unidirectional link between a transmitter, t ∈ V and its associated receiver, r ∈ V . In this paper, we consider only connected graphs. We assume that there is a set H of RBs, indexed from 0 to M , assigned to each edge from E.

B. Interference modeling

Assuming that a signal may be too weak to be received, but still strong enough to interfere with another signal, we consider that a transmitter can interfere with all its two-hop receivers. We say that two nodes are two hops away if they are not neighbors but have at least one neighbor in common. To model interference, we define an undirected graph G C = (V C , E C ) which consists of a set of vertices V C and a set of edges E C . Set V C has a one-to-one relation with the edges in set E of graph G, i.e. for each edge e ∈ E, there is a corresponding vertex v c ∈ V C . We call this graph a conflict graph. An edge (e 1 , e 2 ) from E C exists only if the transmitter from an arc associated with one end of this edge and the receiver from the arc of the other end are nodes which are two hops away. Since each resource is supposed to be independent, we defined a set of conflict graphs G Ci = (V Ci , E Ci ), one graph for each channel i of H. The set V Ci represents the subset of elements of V C which uses RB of index i. Set E Ci is obviously the subset of edges from E C between two elements of V Ci . Formally, the solution proposed in the next section aims at realizing ∀e ∈ E, a(e) = 1 with ∀i ∈ H, E Ci = ∅ with a(e) the function representing the number of resources to be allocated to the arc e.

III. RESOURCES ALLOCATION AND CONFLICTS

CORRECTION

In [START_REF] Pomportes | Self-stabilizing algorithm of two-hop conflict resolution[END_REF] we have proposed an algorithm for the correction of two hop allocation conflicts. In this section we describe firstly an algorithm we propose to allocate resources fairly. Next, we explain the algorithm to resolve interference conflicts, which is based on our previous work [START_REF] Pomportes | Self-stabilizing algorithm of two-hop conflict resolution[END_REF]. The modification we introduce here allows the algorithm to make a better decision on which allocation should be broken when conflict occurs.

A. Allocation of Resource Blocks

We introduce a new algorithm to allocate a single RB on each arc of the connectivity graph. This algorithm, tries to maximize the spatial reuse of resources. Under this constraint, our algorithm can operate even with a reduced number of resources. Maximizing the spatial reuse rate means trying to use the same RB as many times as possible. Taking into account our assumption of two hop interferences, the optimal spatial reuse consists in allocating a given RB every three-hop. However, each new allocation has to take into account the already existing ones. This fact makes simultaneous allocations difficult as the optimal spatial reuse may be obtained to the price of a very long convergence time. The solution proposed in this paper offers a tradeoff between the maximization of spatial reuse and the convergence time.

This tradeoff is obtained by controlling the number of simultaneous allocations. Allocations to a link with a high level of interference are made first. We justify this choice by the fact that since a link shares the resources with its interfering links, the difficulty of assigning an available RB for a link increases with the number of interfering links. In our algorithm, this priority is introduced with a weight W n for each node n which is defined by

W n = outDegree(n) + v∈Neighbours(n) (outDegree(v))
, where outDegree(n) is the number of outgoing links of n. taking into account the allocations which have already been made, each node computes a virtual weight, called W ′ n defined as W ′ n = W n -Number of outgoing links of v which has already received an allocation, if all outgoing links of n have not yet been allocated, and 0 otherwise. Each node n of the network which wants to realize an allocation compares its virtual weight with that of all its neighbors. If there is a neighbor with a higher virtual weight, the node does not have priority and it has to wait. In case of virtual weights equality, the differentiation is done on the weight and finally on the ID. When a node n takes the allocation precedence, it allocates a resource on the outgoing link shared with its neighbor, which has a higher virtual weight. To select the resource to be allocated, the algorithm looks among all the free RBs for the one that maximizes spatial reuse. A resource which maximizes spatial reuse is identified by the states used in the algorithm of conflict correction. If such a resource does not exist, the algorithm takes an available free resource. When performing the choice between multiple resources, the algorithm always chooses the one with the smallest identifier. This property forces the algorithm to maximize use of the same resources across the network, thereby increasing the level of spatial reuse.

B. Correction of two hops Allocation Conflicts

The algorithm for the conflict correction is distributed on each node p of the network. This algorithm associates four variables for each RB indexed by i ∈ H. The variable RBState i (p) provides information on the location of the closest transmitters and receivers. This variable can be in one of the eight different states whose meaning is detailed in Table I. The second variable P A i (p), is an allocation pointer. When an RB is used to transmit (to receive, respectively), P A i (p) points to the associated receiver (the transmitter, respectively). The variables P 1,i (p) and P 2,i (p) represent two priorities. They are used to select more efficiently the allocation which should be retained when conflict occurs. The first priority P 1,i (p) is based on weight W p defined in the previous paragraph. The second priority, P 2,i (p), is used as a

; if CandidatesList = ∅ then if ∃ v ∈ N eighbors p tq W ′ v > W ′ p OR (∃ v ∈ N eighbors p tq W ′ v = W ′ p AND W v < W p ) OR (∃ v ∈ N eighbors p tq W ′ v = W ′ p AND W v = W p AND ID v < ID p )
then candidate ← node in the candidatesList with the highest virtual weight; if ∃ a RB r such as r is usable and maximizes the spatial reuse then Allocation of the RB r on the outgoing link from p to candidate; else if ∃ a RB r such as r is usable then Allocation of the RB r on the outgoing link from p to candidate; is unique for each allocation using the same resource. Indeed, two arcs could share the same priority if their transmitters have the same identifier. However each node has a unique identifier. Therefore, it is impossible that two different arcs share the same priority P 2,c for the same resource c. Precisely speaking, when conflict occurs, P 1,i (p) preserves the allocation with the highest priority. When multiple allocations sharing the same value of P 1,i (p) are in conflict, P 2,i (p) is used to remove uncertainty by selecting the allocation which has the higher identifier.

The algorithm is composed of three steps. The first one checks that each RB already assigned is not in conflict with another of higher allocation priority. The second step updates the variables of unused RBs. When a node p detects a transmitter pointing to it, the third step determines whether p has to be associated with that transmitter or whether p has to ignore this request.

IV. SIMULATIONS RESULTS

To validate our algorithm, we developed a simulator in Java. This simulator is available at [15]. We performed numerous simulations for chain, grid and random topologies. To generate random topologies, the simulator places a node inside the a simulation area of size 1000 × 1000, uniformly in both dimensions. When the distance between two nodes is within a threshold called radio range, these nodes may communicate the maximal clique in the conflict graph. The maximal clique corresponds to the maximum number of mutually interfering links in the network. To check whether our solution is scalable, we should increase the size of the network without changing the size of its maximal clique. For chain and grid topologies, the size of the maximal clique reaches a constant value but this is not the case for random topologies. Moreover, it is impossible to extract the random topologies of same maximal clique size because the computation of this size is an NPcomplete problem [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. Nevertheless, the size of the network's maximal cliques partially depends on the nodes' average degree. Therefore, in the simulations, we vary the radio range to maintain a constant nodes degree. Figures 4a, 5a and 6a show the evolution of the convergence time depending on the number of nodes for different topologies. For the chain, we choose to isolate the most unfavorable case where the ID of the nodes are ordered linearly (Fig. 4).

In this case, time increases linearly with the number of nodes. Indeed, we have seen that a node has allocation precedence if it has the highest virtual weight compared with the nodes in its neighborhood. We have also seen that when several candidates share the same virtual weight comparison is made on weight and, finally, on the IDs. In a chain, weight and virtual weight are identical for all nodes except the first two and the last two which have fewer neighbors two hops away than the other nodes. As the final comparison is made according to the ID, only the candidates with the highest IDs can perform it. If the IDs are randomly distributed, simultaneous allocations are possible dividing the convergence time by the number of simultaneous allocations in the network. In Figs. 5a and6a we can observe a tendency of the convergence time depending on network size for the grid and random topologies. In fact, the convergence time increases with the number of links, but it decreases with the number of simultaneous allocations. Thus, 6: Results for random topologies with constant average node degree when the size of a network increases, the increase of the number of links is balanced by the number of simultaneous allocations which also increases.

The simulations show that the resource demand for a random topology is greater than the demand for a grid, which is itself larger than the demand for a chain. The main reason is that the mean degree of nodes follows the same tendency as the RB requirement and indeed, the average degree is two for the chain, four for the clique and five for the random topologies. Concerning the scalable property of our solution, we can observe a convergence of the quantity of RB requirements when the network size increases. This is due to spatial reuse which allows the same resource to be used more than once in the network. As a number of simultaneous allocations is balanced by the increase in the number of links for simulation time, the increase in spatial reuse balances the increase of RB requirements.

V. CONCLUSION

In this paper we presented our work concerning the improvement of the quality of service guaranteed in ad hoc networks. It was motivated by our previous studies of ad hoc networks deployed in disaster areas. We proposed an RB allocation and a conflict correction algorithm which are dedicated to networks with unidirectional links, assuming that the interference range is twice as large as the communication range. We developed a simulator to validate our solution and we evaluated its performance on chain, grid and random topologies. The simulation results show that the convergence time and number of RBs required converge asymptotically when the network size increases. It also shows that the value to which our solution converges depends considerably on the average degree of network nodes. The algorithms proposed here provide a solution to allocate one RB to each arc of the connectivity graph as we were principally interested in the fairness of the proposed reservation. Our future work will focus on creating a new algorithm which allocates groups of RBs to links and which takes into account a traffic matrix.

Fig. 1 :

 1 Fig. 1: Example of chain topology

Fig. 2 :Fig. 3 :Fig. 5 :

 235 Fig. 2: Example of grid topology

  Fig.6: Results for random topologies with constant average node degree

TABLE I :

 I Description of variable RBState

	State	Scope	Meaning
	T	local	slot locally used for transmission
	R	local	slot locally used for reception
	N R1	one hop	at least one neighbor uses the slot for reception
	N T 1	one hop	at least one neighbor uses the slot for transmission
	N T R	one hop	slot used by a transmitter and a receiver in neighborhood
	N T 2	one hop	slot used by a transmitter which is two hops away
	N R2	two hops	slot used by a receiver which is two hops away
	F	-	there is no transmitter and no receiver two hops away

Algorithm 1: Allocation algorithm CandidatesList ←Neighbors p \ {neighbors for which the outgoing link has already received an allocation.}

time is configurable, we decided to represent the convergence time in number of rounds. One round is a time interval required for all nodes to send a message. In further figures, every point is the average of 100 simulations for chain and grid, and 1000 for random topologies. We calculated the 95% confidence interval for these averages.

In our simulations, we are mainly interested in the evolution of the number of RBs required depending on the size of the network. Theoretically, this quantity depends on the size of