Jamila Ben Slimane

Ye-Qiong Song

Anis Koubaa

Mounir Frikha

Joint duty cycle Scheduling, resource Allocation and multi-constrained QoS Routing algorithm

Wireless mesh sensor networks (WMSNs) have recently gained a lot of interest due to their communication capability to support various applications with different Quality of Service (QoS) requirements. The most challenging issue is providing a tradeoff between the resource efficiency and the multiconstrained QoS support. For this purpose, we propose a cross-layer algorithm JSAR (Joint duty cycle Scheduling, resource Allocation and multi-constrained QoS Routing algorithm) for WMSNs on based multi-channel multi-time slot MAC. To our best knowledge, JSAR is the first algorithm which simultaneously combines, a duty cycle scheduling scheme for energy saving, a resource allocation scheme for efficient use of frequency channels and time slots, and an heuristic for multi-constrained routing protocol. The performance of JSAR has been evaluated, showing that it is suitable for on-line implementation.

Introduction

Wireless mesh sensor networks (WMSNs) are expected to support various applications with different QoS requirements. According to novel application requirements, QoS constraints become more and more critical in terms of end-to-end delay, data throughput and packet-error-rate. Also, due to energetic constraints at node level, energy saving remains the most challenging issue. Various cross-layer designs are applicable to WM-SNs [START_REF] Luo | Cross-Layer Optimization for Energy-Timeliness Tradeoff in TDMA Based Sensor Networks[END_REF][START_REF] Bouabdallah | Cross-Layer Design for Energy Conservation in Wireless Sensor Networks[END_REF][START_REF] Yuebin | An Energy Optimization Protocol Based on Cross-Layer for Wireless Sensor Networks[END_REF][START_REF] Gang | Minimum latency joint scheduling and routing in wireless sensor networks[END_REF], but the majority of proposed approaches generally focus on WMSNs based single-channel MAC protocol and only [START_REF] Gang | Minimum latency joint scheduling and routing in wireless sensor networks[END_REF] focus on Multi-channel based WMSNs. However dealing with all the WMSN's requirements requires tight collaboration of all network's layers. Although several cross-layer designs have been proposed for WMSNs to ensure energy efficiency or to improve the network's performance or to support QoS guaranties, few approaches have taken into account simultaneously multiple WMSN's requirements. The authors in [START_REF] Luo | Cross-Layer Optimization for Energy-Timeliness Tradeoff in TDMA Based Sensor Networks[END_REF] propose a cross-layer strategy that explores the tradeoff between energy efficiency and packet timeliness in time division multiple access (TDMA) based WSNs, by transmission power allocation and routing path selection schemes. The idea in [START_REF] Bouabdallah | Cross-Layer Design for Energy Conservation in Wireless Sensor Networks[END_REF] is based on a traffic balancing inside the network and a judiciously allocation of the retry limit to each link. [START_REF] Yuebin | An Energy Optimization Protocol Based on Cross-Layer for Wireless Sensor Networks[END_REF] proposes a cross-layer design based on power transmission management, routing, and duty-cycle schedule to optimize the WSNs energy-efficiency. We find that only [START_REF] Gang | Minimum latency joint scheduling and routing in wireless sensor networks[END_REF] proposes a cross-layer design for WM-SNs based on multi-channel access. The goal of the joint wakeup/sleep scheduling and routing algorithm proposed in [START_REF] Gang | Minimum latency joint scheduling and routing in wireless sensor networks[END_REF] is the minimization of the communication latency while providing energy efficiency for nodes in FDMA-based multi-channel WMSNs. The proposed routing algorithm takes into account only delay constraint to minimize communication latency and save energy. We note that the flow differentiation, the link reliability, the residual energy per node and the energy consumption per path have not been taken into account during routing process. None of these approaches combines simultaneously multi-constrained QoS routing, duty cycle scheduling and efficient frequency channels and time slots allocation.

The main problem, we address, is designing an efficient cross-layer algorithm that simultaneously takes into account various WMSN's requirements. For a solution of such NP-hard problem, we propose, in this paper, JSAR, an algorithm that combines simultaneously a network duty cycle scheduling policy for energy saving and so network lifetime maximization, a frequency channels and time slots allocation strategy for an efficient sharing of frequency channels and time slots during routing process, and an heuristic for a QoS support offering guaranteed services in compliance with the application requirements and taking into account the network configuration. In this current work, we only focus on the performance evaluation of JSAR. The network performance evaluation under JSAR is our future work. The rest of the paper is organized as follows: In section 2, we present the system model. Section 3 describes the principle of PMCMTP [START_REF] Ben Slimane | A Prioritized Multi-Channel Multi-Time slot MAC Protocol For Large-Scale Wireless Sensor Networks[END_REF] with the use of JSAR. In section 4, we detail the proposed cross-layer design. In section 5, we evaluate JSAR performance by analyzing and commenting some simulation results.

System model and notations

We consider an WMSN composed of one coordinator located at the center of the network and a set of distributed sensor nodes acting as routers and/or sources organized in fully meshed topology (See second tier of the network architecture proposed in [START_REF] Ben Slimane | A Three-Tiered Architecture for Large-Scale Wireless Hospital Sensor Networks[END_REF]). We can model the network as a weighted directed graph G(V, E, w). We admit that the network supports five data flow priority levels from P 1 to P 5 (P 1 for the highest priority level). The control traffic refers to the network synchronization and resource allocation information. All network nodes are concerned by the beacon frames. So, a simple manner to quickly share such information between all nodes is broadcasting it from the coordinator which is more powerful than others nodes, and has no energetic constraint compared to others nodes. Also, the amount of control traffic is negligible in comparison to the data traffic so there is no need to balance such traffic over the network. A single frequency channel will be enough to ensure the control traffic and the rest of available frequency channels can be used to maximize the number of parallel data communications. For those reasons, we propose one hop control communications using one frequency channel. To balance both load and energy consumption and to minimize power transmission, we propose multi-hop routing with the use of J frequency channels of K time slots per channel for data communication between network's nodes. Each node u ∈ V has a unique identifier (id), and it is characterized by its initial battery capacity (BC u ), residual energy (Re u ) and data queue state (DataQ u ). A node presents three activity states respectively busy (transmitting or receiving data), forced sleep (see section 3.A) and free (in sleep).

Notations:

V: set of vertices (nodes), E: set of edges (link) N: number of nodes, L: number of links, w: link's weight,P i : flow priority levels, J : number of frequency channels, K : number of time slots per channel dc u : refers to a node duty cycle or activity state, dc u = 1 for sleep or free state dc u = 0 for busy and in forced sleep states Each link l uv ∈ E is associated with multiple parameters which can be classified into: additive metric (as delay, link cost, energy consumption, number of hops), multiplicative metric (as packet reception rate P RR) and concave metric (as available bandwidth [START_REF] Dr | Exploring Issues for QoS Based Routing Algorithms[END_REF]). For non-additive metrics, we propose to prune all links that do not satisfy the constraints of concave metrics and for multiplicative metrics we transform it into additive metrics by applying a logarithm operation on it. For the rest of the paper, we consider only additive metrics for path evaluation. Each link l uv ∈ E is associated with a cost parameter cost(u, v) and m additive nonnegative weights w i (u, v) ≥ 0, i = 1..., m which can be denoted by the vector w(u, v) = w 1 (u, v), ..., w m (u, v) and each weight is corresponding to a constraint, and the upper bound of the i th constraint is denoted by C i (i.e vector C = C 1 , ..., C m corresponds to the vector of the m constraints). At a given time t, each link l uv ∈ E is characterized by a set of parameters as follows:

-Link's logarithmic value of packet reception rate prr uv : We define prr uv as the absolute value of the logarithmic value of P RR uv (packet reception rate).

prr uv = |log 10 (P RR uv )| (1)
-Link's delay δ uv : For each data-flow (F ), the wireless interface delay would typically involve the data transfer delay (δ F tr ), the resource allocation and the channel access delay (δ F ch ). δ uv is expressed in terms of time slots as follows:

δ uv = δ F tr + δ F ch (2)
-Link's energy consumption E c uv : We define E c uv as the sum of the energy consumed by node u during radio transmission and data processing and queueing (E c u ) and the energy consumed by node v during radio receiving (E c v ).

E c uv = E c u + E c v (3)
-Link's residual energy Re uv : We define Re uv as the minimum value of residual energy of nodes u and v.

Re uv = min(Re u , Re v ) (4)
-Link's energy availability Ae uv : We define Ae uv as a boolean value reflecting the energy availability for both nodes u and v.

Ae uv = 1 if Re u > E c u and Re v > E c v 0 otherwise -Link's available time slots T S uv :
We assume the use of K equal time slots (τ ) with J frequency channels (i.e. KxJ τ in total). We consider T S j u (resp.T S j v ) to represent the available time slots of the node u (resp.v) of the j th channel and T S j uv to represent the link's available time slots of the j th channel. We model T S j u (resp.T S j v ) as a binary vector of dimension K:

T S j u = tsu j,k 1≤k≤K tsu j,k = 1 if time slot τ j,k is available 0 else (5)
τ j,k represents the k th time slot of the j th channel.

T S j uv = tsu j,k * tsv j,k 1≤k≤K (6)
At a given time t, we characterize a path P s→d , from a source s to a destination d, by a set of parameters as follows:

-Path's packet reception rate prr s→d : We transform this multiplicative metric into a cumulative metric by applying a logarithm operation on it.

P RR s→d = luv∈P s→d P RR uv ⇒ prr s→d = luv∈P s→d prr uv (7)
-Path's delay δ s→d : We define path's delay as the sum of delay expired per each link in the path.

δ s→d = luv∈P s→d δ uv (8)
-Path's energy consumption E c s→d : We define path's energy consumption as the sum of energy consumed per each link in the path.

E c s→d = luv∈P s→d E c uv (9)
-Path's residual energy Re s→d : We define path's residual energy as the minimum of residual energy of the set of all links in the path.

Re s→d = min

luv∈P s→d Re uv (10)
-Path's energy availability Ae s→d : We define path's energy availability as the minimum of energy availability of the set of all links in the path.

Ae s→d = min

luv∈P s→d Ae uv ( 11
)
3 PMCMTP with the use of JSAR

In [START_REF] Ben Slimane | A Prioritized Multi-Channel Multi-Time slot MAC Protocol For Large-Scale Wireless Sensor Networks[END_REF], we have proposed the PMCMTP for frequency channels and time slots allocation inside a personal area network. A key concept in PMCMTP is the elementary active cycle as shown in Figure 1, which is composed of two consecutive superframes, the first for synchronization and collection of resource requests, and the second for the request scheduling algorithm (RSA), the reception of the second beacon and the data communications. In this paper, we propose the extension of the PMCMTP MAC protocol by including the JSAR instead of RSA to ensure the access to the medium, the duty cycle scheduling of network's nodes, the efficient sharing of available frequency channels and time slots and the search of the best routing paths in compliance with multiple application QoS requirements. Let us assume that the coordinator maintains a global vision of the network configuration. Whenever there is a change in the network, the coordinator updates related parameters (matrix and graphs...). The principle of PMCMTP using the JSAR is based on the following four phases: which, the coordinator is listening to the requests of its members. Based on the number of source nodes, the coordinator assigns to each source node a specific time slot according to their id. Just following the reception of the first beacon, each source node waits for its own time slot to send its resource allocation request. The request packet is composed of five fields including the request identifier, the request's priority index, the number of required time slots per hop, the source and destination addresses. -Routing and resource allocation process: According to the JSAR, after reception of all resource allocation requests, the coordinator schedules them according to their priorities. Once the list of requests are scheduled, the coordinator launches the routing process, an heuristic algorithm based on Dijkastra algorithm, to find the best path, then it ensures the effective allocation of time slots and frequency channels. For each request, it tries to find the adequate available time slots per available channels to assign it for the entire selected path from the source to the destination taking into account QoS constraints and the resource availability. After processing all resource allocation requests, the coordinator computes the satisfactory note SN of the current active cycle, then it registers a trace of requests, which were not served to analyze it during the next cycle. Then, it inserts into the next beacon frame the necessary information of the served requests as the index of allocated channel, the index of the first allocated time slot, the number of allocated time slots and source and destination addresses. -Data transmission: After listening to the second beacon, nodes can have a feedback of different requests. Each concerned sensor switches to the suitable channel at the suitable time slot and it begins sending or receiving data frames during the allocated duration.

WSN's cross-layer design

In this section, we propose a new WMSN's cross-layer framework, with joint consideration of the network's duty cycle scheduling policy, the resource allocation management strategy and the network layer QoS routing process to simultaneously ensure energy saving, optimal resource sharing and QoS support. The MAC layer has to support the duty cycle scheduling and the spectral and temporal resource allocation whereas the Network layer has to ensure the multi-constrained QoS routing process.

Duty cycle scheduling

The objective of the network's duty cycle scheduling is to maximize the node's sleep duration in order to prolong the network's lifetime while ensuring QoS requirements in a dynamic manner (i.e by balancing load and enhancing energy saving without affecting to the network performance). The dynamic Sleep/Wake scheduling problem in WMSN can be formulated as follows: how should the activity of nodes be scheduled in a WMSN in order to enhance energy efficiency while ensuring network's QoS requirements?

To determine which subset of nodes to be turned off for a time interval and the manner of sleep interval scheduling are both in joint correlation with routing process and application requirements. In order to minimize and balance energy consumption and to avoid the node's failure due to draining of battery, we suggest forcing certain routers to enter into sleep state. The nodes, whose residual energy is lower than the threshold value T H n and their data queues (DataQ i ) are empty, are turned into forced sleep state during the n th elementary active cycle (See Alg.1). T H n is defined as follows:

T H n = N i=1 Re n i N * SN n-1 , SN n-1 =
N umber of satsif ied routing requests T otal number of routing requests (12)

N and SN n-1 refer respectively to the number of WMSN's nodes and the satisfactory note of the n -1 th active cycle (SN n-1 ∈]0, 1]). The satisfactory note SN n-1 is computed at the end of the n -1 th active cycle.

Let us assume that the duration of the n th data communication phase is equal to K time slots with the use of J frequency channels. We model the duty cycle during data communication phase of the node i by the vector DC n i of dimension K.

DC n i = dc n i,k 1≤k≤K dc n i,k = 0 in f orced sleep or busy 1 if f ree or in sleep (13)
The duty cycle matrix DC n gives an idea about the activity states of the WMSN's nodes. This matrix can be updated during routing process (i.e. by putting the adequate time slots to 0, for each node participating in a path).

DC n =     DC n 1 . . DC n N     N K =     dc n 1,1 .. dc n 1,K . .. . . .. . dc n N,1 .. dc n N,K     N K (14)
Algorithm 1 : Forced Sleep Test, FST

for i = 1 to N { if (Re n i < T H n and DataQi is empty) then DC n i = 0 , update(G(Va, Ea, w) n )}
DC n i : State of the i th node f or the n th elementary active cycle Va,Ea: T he active vertices and edges of the n th active cycle

Resource Allocation

The objective of the resource allocation process is to optimally allocate available resource in terms of time slots per available frequency channels to ensure QoS support and to enhance both energy efficiency and channel utility. As Sleep/Wake scheduling problem, the resource allocation is also in joint correlation with routing process and application requirements. We model the resources of the the n th data communication phase by the binary matrix RS n JK . Initially, RS n JK = I JK (i.e. initially there are in total KxJ available time slots) which is updated during routing process.

RS n JK =     rs n 1,1 .. rs n 1,K . .. . . .. . rs n J,1 .. rs n J,K     JK rs n j,k =
1 if the k th slot of the j th channel is f ree, 0 otherwise.

(15) For a given link l uv ∈ E, we compute the link's available time slots T S uv , for the n th active cycle, as given in Alg.2.

In order to maximize the resource utility and to optimally share such resource between WMSN's nodes, we suggest to implement a centralized resource allocation policy at the network coordinator. The resource allocation decision is taken by the coordinator in response to the set of collected resource allocation requests, according to the network configuration, the duty cycle scheduling and the routing decision. As shown in Fig. 2, first, the coordinator tries to find candidate paths according to the application requirements and network organization. During this process, it ensures the temporary resource allocation (allocation of the most earlier available time slots of available frequency channel per link) to compute path's parameters and cost, then it selects the best path. At this moment, the coordinator makes the effective resource allocation and it updates the resource matrix RS, the duty cycle matrix DC and the energetic graph G e (i.e. by decreasing the residual energy, according to the energy consumption estimation, for each node participating in the path). We propose TRR (Temporary Resource Reservation) algorithm (Alg.3) to ensure the temporary resource reservation of the current path that can be added as extensible or candidate path. This algorithm consists to find the suitable available time slots per available channels to temporary assign it to the new edge l uv . The main objective of this algorithm is to minimize, as possible, the end to end delay.

Multi-constrained QoS routing

The objective of QoS routing is to find a path from a source s to a destination d using the minimum amount of network resources such as energy consumption per route, residual battery power per node, available time slots, while satisfying various QoS constraints, such as delay, reliability, etc. When multiple routing metrics are considered, the problem becomes a multi-constrained path problem, which has been mathematically proven in [START_REF] Wang | QoS routing for supporting resource reservation[END_REF] to be NP-complete. However, various heuristic and approximation algorithms, proposed in the literature, can approximate or solve similar problems in polynomial or pseudo-polynomial time [START_REF] Wendong | An efficient heuristic algorithm for multi-constrained path problems[END_REF][START_REF] Wan | Approach for Multiple Constraints Based Qos Routing Problem of Network[END_REF]. For WMSNs, we propose a centralized multi-constrained QoS reactive routing, which jointly interacts with the proposed duty cycle scheduling and resource allocation methods. The algorithm JSAR, given by Alg.4, describes the proposed routing approach. Let P G and P s→d represent respectively the set of all paths and the set of all paths from a source node s to a destination node d in the network. P G = ∀s,d∈V P s→d .

The basic routing problems can be defined as follows:

Definition 1. The MCP (Multi-Constrained Path) problem is to find paths P s→d from s to d such that:

w(P s→d ) = luv∈P s→d       w 1 (u, v) w 2 (u, v) . . w m (u, v)       ≤       C 1 C 2 . . C m       (16)
The goal is to find the set of feasible paths P f s→d from s to d that satisfy multiple constraints simultaneously.

Definition 2. The MCOP (Multi-Constrained Optimal Path) is one variant of M CP problem which is to find a path P opt s→d from s to d such that:

Cost(P opt s→d ) = min Cost(P s→d )

∀P s→d ∈P f s→d ( 17
)
The goal is to find the optimal feasible path from s to d in such a way that network resources are efficiently utilized. According to Lemma (1), we can use the concept of path domination to reduce the computational complexity by considering only non-dominated feasible paths, given that the set of non-dominated feasible paths from s to d represents a sub set of P f s→d . We define P f nd s→d as the set of all non-dominated feasible paths in P f s→d .

Lemma 2. Given P s→d ∈ P f nd s→d , let P u→v ⊂ P s→d =⇒ P u→v ∈ P f nd u→v (Optimality principle: Every optimal path is formed from optimal sub-paths).

According to Lemma (2), a sub-path of a non-dominated feasible path is even a nondominated feasible path. We consider only non-dominated feasible sub-paths to be extended to reach last destinations. Sub-optimal paths (i.e, paths that are dominated by or are equal to others) are ignored during path's extension process. Let two paths P and P to vertex u and P is dominated by P , if path P can be extended to a path that satisfies Eq. 16, then so also can P . So, there is no need to retain P path extension. The determination of all non-dominated paths (i.e P f nd s→d ) from s to d becomes a very hard task to accomplish when the number of such paths is too high. In this context, we propose to limit the search of non-dominated feasible paths for a set of at most X paths (P Candidate s→d ). The main sub-problems that we must study are as follows:

-How to define the parameter X taking into account its impact on the routing decision and consequently on the network performance, -How to evaluate a path, -How to evaluate and select the set of the X most efficient non-dominated feasible paths if they exist.

Once the set P Candidate s→d is obtained, one still faces the problem of selecting the best final solution. Definition of the parameter X: To improve the efficiency of the path finding algorithm, X should be large, but it is to the detriment of a computation complexity. This parameter must be carefully determined in order to ensure a tradeoff between routing efficiency and both space and computation complexity. We propose two methods to compute X:

-Statically: The coordinator defines X as a fixed number computed as soon as the beginning of the network operation. We can use the definition used in [START_REF] Wendong | An efficient heuristic algorithm for multi-constrained path problems[END_REF]. -Dynamically: For each active cycle, the coordinator can dynamically compute X according to the network configuration, the supported load, and the satisfactory note of the previous active cycle.

Evaluation of a path: First, we propose a link's cost function. For resource allocation requests with the priority's level equals to P 1 or P 2 (i.e. hard real time constraint), we consider the delay as the metric of performance. Then, for resource allocation requests with the priority's level equals to P 3 or P 4 , we propose the energy consumption per link as the metric of performance given that energetic aspect becomes constraining. Finally, for the rest of requests we propose the residual energy per node as the metric of performance to avoid premature death of some nodes.

Cost uv =    δ uv if P = P 1 , P 2 E c uv if P = P 3 , P 4 Re uv if P = P 5 (18)
To ensure jointly QoS support, load balancing and energy saving, we propose the link's cost function (Eq.18) and the path's cost function (Eq.19) based on the resource allocation request's priority. For requests with priority P 1 and P 2 , the path cost represents respectively the end to end delay and the average delay. For requests with priority P 3 Constraints are randomly generated and five different ranges of the constraint vector C (< C 1 , C 2 , C 3 >) are selected for each experiment (Table .2). The result is based on 100 randomly generated resource allocation requests for each range. The first experiment is carried out to test how X value affects the performance of JSAR. We compared JSAR with the exact algorithm where the length of path queue is infinite. According to Table .3, we note that the SN increases as the constraints get looser. Also, the SN increases with the increase of X. When X is equal to 5, JSAR is almost as good as the exact one. According to Fig. 3, we observe that the execution time of JSAR increases with the increase of X, also it increases as the constraints get looser. Globally, the execution time of JSAR is acceptable for typical applications for its on-line implementation. The second experiment is carried out to respectively test how the number of available channels and the number of time slots affect the performance of JSAR. According to Table .4 and Table .5, in the majority of cases, the SN increases with the increase of respectively the number of available channels and the number of time slots. Respectively, the increase of the number of available channels and the number of time slots can significantly increase the performance of the algorithm especially when the constraints get looser. According to the network configuration and the application requirements (e.g. desired SN ), those results can be used in practice to compute the necessary resources.

Conclusion

In this paper, we presented a new cross-layer design for WMSNs based on multichannel MAC protocol. We have decomposed the optimization problem into three subproblems: duty cycle scheduling to save energy and extend the network lifetime, frequency channels and time slots allocation to optimally allocate resource and enhance the network's performance and multi-constrained QoS routing to ensure QoS support and save energy. We solved the three problems jointly and using JSAR that simultaneously satisfies multiple QoS requirements and ensures the optimal resource allocation. The evaluation results demonstrate the performance of JSAR. As future work, we will implement the proposed algorithm in a network simulator to evaluate the performance of the network in terms of end to end delay, throughput, lifetime and resource utility.

Algorithm 2 :Fig. 2 .

 22 Fig. 2. Resource allocation algorithm for WMSN

Elementary active cycle First beacon Resource request Time slot Second beacon + RSA process Data transmission Time slot Resource Request phase Data Transmission phase J data Channels K equal Time slots Using a control Channel First SuperFrame Second SuperFrame

Fig. 1. An elementary active cycle -Forced Sleep test and synchronization: Let us assume that the WMSN's coordinator maintains an updated energetic graph G e of network's nodes. Before the start of each elementary active cycle, the coordinator executes a Force Sleep Test (FTS) to exclude nodes with low levels of energy from participating in the routing process. It must include in the beacon payload the id of nodes that must go to forced sleep state. By listening to the first beacon, WMSN's nodes adjust their wake-up clocks and check the beacon payload to know if they must go to forced sleep state. Only nodes that pass the F ST will be active during the current elementary active cycle and the others must go to forced sleep state to save their batteries. -Requesting resources: This phase represents a set of equal short time slots, during

  Definition 3. Let P s→d and P s→d two different paths from s to d. Path P s→d is dominated by path P s→d if and only if w(P s→d ) ≤ w(P s→d ) (i.e. P s→d ≺ P s→d ) with at least one strict inequality. Lemma 1. if ∃P s→d ∈ P f s→d =⇒ ∃ a non-dominated P s→d ∈ P f s→d .

Table 1 .

 1 Network configuration

	Parameter	Default value
	Topology	ANSNET model
	Number of Vertices -Edges	32-54
	Number of Time Slots per data communication phase	15
	Number of Channels per data communication phase	3
	Length of a path queue or X	5
	Number of routing and resource allocation requests	100

Table 2 .

 2 Range of C1, C2 and C3 for each experimentCaseRange of C =< C1, C2, C3 >)1 C1 ∼ unif orm[1, 2],C2 ∼ unif orm[1, 3],C3 ∼ unif orm[40, 120] 2 C1 ∼ unif orm[2, 3],C2 ∼ unif orm[3, 6],C3 ∼ unif orm[120, 240] 3 C1 ∼ unif orm[3, 4],C2 ∼ unif orm[6, 9],C3 ∼ unif orm[240, 360] 4 C1 ∼ unif orm[4, 5],C2 ∼ unif orm[9, 12],C3 ∼ unif orm[360, 480] 5 C1 ∼ unif orm[START_REF] Gang | Minimum latency joint scheduling and routing in wireless sensor networks[END_REF][START_REF] Dr | Exploring Issues for QoS Based Routing Algorithms[END_REF],C2 ∼ unif orm[12, 15],C3 ∼ unif orm[480, 600]

and P 4 , the path cost represents respectively the average energy consumption per path and the total energy consumption per path. Finally for last requests, in order to avoid network partition, the cost represents the residual energy per path. Evaluation and selection of most efficient non dominated X paths: According to the relax operation (See Alg.4), each new generated path ((current path p) ∪ l uv ) must pass the feasibility test then the non-dominance test. In this case, the new founded path is stored, according to its cost, in the suitable position at the scheduled queues PATH(v) (i.e P f s→v ) and Q (non-dominated paths). The length of queues can not exceed X, otherwise, the last path must be discarded.

N hop : the number of hops per a path P s→d from s to d. Selection of the best path: The best path can be selected among a set of candidate paths. The selection's criterion is the path's cost function. So, the minimum cost path will be selected as the best path. Consequently, the effective resource allocation takes place.

JSAR

As already explained, for each active cycle, the JSAR (See Alg. 4) performs firstly the Forced Sleep Test (line 2 in the main). After the reception of all resource allocation requests, it computes the parameter X (line 3 in the main). For each request it tries to find the set of candidate paths (line 4 in the main) then, from this set it selects the best path and ensures the effective resource allocation (line 5 in the main). Finally it computes the satisfactory note of the current cycle (line 6 in the main).

JSAR evaluation

We have implemented the proposed algorithm in MATLAB environment. The performance of the algorithm is evaluated by experiments on an ANSNET [START_REF] Comer | Internetworking with TCP/IP, 3 rd Edition[END_REF] network (Table.1). We consider a three-constrained path problem. Link weights are as follows:

w 1 (u, v): represents the parameter prr uv randomly selected from uniform [0, 1], -w 2 (u, v): represents the parameter δ uv assumed be equal to 1 Time slot, -w 3 (u, v): represents the parameter E c uv assumed be equal to 40 . --------------X efficient PATHs(X, q, RStemp, DCtemp, Ge temp , Queue)

11 Remove q, RStemp, DCtemp,G etemp from Q 12 if(!=dominated) then 13 q, RStemp, DCtemp, Ge temp = p ∪ luv, RStemp, DCtemp, Ge temp 14 X efficient PATHs(X, q, RStemp, DCtemp, Ge temp , P AT H(v)) 15 X efficient PATHs(X, q, RStemp, DCtemp, Ge temp , Q) --------------- --------------


S F Eskildsen

P Coupé

V Fonov

L R Østergaard

D L Collins

Effect of non-local means denoising on cortical segmentation accuracy with FACE

Introduction

Cortical thickness measurements based on MRI have the potential to accurately detect changes in the cortical gray matter. During the last decade, several methods have been developed to estimate cortical thickness using surface based reconstruction [START_REF] Dale | Cortical surface-based analysis: I. Segmentation and surface reconstruction[END_REF][START_REF] Kim | Automated 3-D extraction and evaluation of the inner and outer cortical surfaces using a Laplacian map and partial volume effect classification[END_REF][START_REF] Eskildsen | Active surface approach for extraction of the human cerebral cortex from MRI[END_REF]). Due to the high sensitivity of reconstruction methods to acquisition artifacts, noise in MRI signal may directly affect the obtained measurements. Usually, this noise is reduced using various smoothing filters. Recently, the non-local means (NLM) filter has demonstrated very high denoising performance compared to these traditional filters (Coupe 2008). In this study, we investigated the sensitivity to noise of cortical surface reconstruction and the ability of the NLM filter to reduce the effect of noise.

Methods

Sixteen repeated T1-weighted scans with isotropic 1mm voxels of the same subject were acquired over a short period of time [START_REF] Holmes | Enhancement of MR images using registration for signal averaging[END_REF]) and were used to create an average image with very high signal-tonoise ratio. The 16 scans were bias field corrected [START_REF] Sled | A nonparametric method for automatic correction of intensity nonuniformity in MRI data[END_REF], intensity normalized, registered to ICBM space and aligned using rigid transformations [START_REF] Collins | Automatic 3D intersubject registration of MR volumetric data in standardized Talairach space[END_REF]) before computing a voxel-by-voxel intensity average dataset with greatly reduced noise (reference image). In order to simulate several levels of noise, the reference image was added Rician noise (1-10%). The optimized blockwise NLM filter (Coupe 2008) was then applied to the 10 simulated noisy images as well as the 16 original scans. FACE (fast accurate cortex extraction, Eskildsen 2006) was applied to the reference image, the simulated noisy scans with and without NLM filtering, and the original 16 images with and without NLM filtering to generate cortical surfaces delineating the WM/GM and the GM/CSF boundary. Cortical surfaces extracted from the reference image were considered as the reference during comparison.

To measure the cortical segmentation error, the root mean square (RMS) of the distances from the vertices of the test surface to the reference surface was calculated for both WM and GM surfaces. Furthermore, the surfaces were voxelated onto a 1mm isotropic voxel grid to enable calculation of the Dice coefficient (DC) [START_REF] Dice | Measures of the amount of ecologic association between species[END_REF]. The WM, GM, and CSF classifications (fuzzy c-means), which are a part of the cortical extraction process, were also evaluated by DC.

Results

Fig. 1 shows input images, tissue classification images, and images with cortical surfaces overlaid for the reference image built on the 16 repeated scans, the reference image with 9% of added Rician noise before and after NLM filtering. Note that the GM surface has difficulties capturing the GM/CSF boundary in the noisy image, while after denoising the delineation is improved. This improvement can be seen on Fig. 2 and3, where the segmentation accuracy on denoised images is better for all the noise levels in terms of DCs and RMS distances. The mean noise level over the 16 repeated scans was estimated at 4.8±0.3% using the noise estimation proposed in (Coupe 2010). This level is indicated on the charts in Fig. 2 and3. Fig. 4 shows the effect of denoising the 16 original scans on the segmentation accuracy compared to the reference image. Denoising improved the accuracy significantly (p<0.001, paired t-test) for both WM and GM surfaces.

Conclusions

The results demonstrate that applying NLM filtering to the MRI scans improves the accuracy of cortical surfaces extracted using FACE. For the experiments with synthetic noise, the RMS distance to the reference was below 0.5mm for noise levels up to 10%, which covers the realistic range of noise present in MRI scans. Segmentation accuracy of the original scans was also improved by denoising, which resulted in RMS distances of 0.44±0.03mm and 0.42±0.03mm for GM surfaces and WM surfaces compared to 0.49±0.03mm and 0.52±0.06mm for the unfiltered scans. The RMS distances obtained are lower than previous experiments using another cortical reconstruction method on the same data [START_REF] Kim | Automated 3-D extraction and evaluation of the inner and outer cortical surfaces using a Laplacian map and partial volume effect classification[END_REF].


Vladimir Fonov

Pierrick Coupé

Simon F Eskildsen

D Louis Collins

Atrophy specific MRI brain template for Alzheimer's disease and Mild Cognitive Impairment

Keywords: anatomical template, atrophy, magnetic resonance imaging, Alzheimer's

Learning objectives:

1. To create a brain template specific to the different levels of atrophy in the brains of patients with Mild Cognitive Impairment and Alzheimer's disease . 2. To study the pattern of brain change correlated

Background

Rapid brain loss is characteristic for the patients with mild cognitive impairment (MCI) and Alzheimer disease (AD) [START_REF] Carlson | Trajectories of brain loss in aging and the development of cognitive impairment[END_REF]. In particular, increase of the lateral ventricular volume is strongly correlated with the progression of the disease. Other parts of the brain (hippocampus, entorhinal cortex, etc.) are also subject to the rapid neuronal loss. Many automated image processing methods rely on usage of an anatomical template for registration into a common stereotaxic space for tissue classification, structure segmentation or voxel-based morphometry (VBM). Unfortunately, most commonly available anatomical templates are based on the scans of healthy individuals: for example the MNI-ICBM152 template used in the publicly available SPM and FSL image analysis tools is based on MRI scans of young healthy adults [START_REF] Fonov | Unbiased average age-appropriate atlases for pediatric studies[END_REF], and thus may be suboptimal for image registration techniques. The use of disease-specific templates have been proposed in the literature [START_REF] Hua | 3D characterization of brain atrophy in Alzheimer's disease and mild cognitive impairment using tensor-based morphometry[END_REF], however high variability in the degree of atrophy for subjects with AD and MCI makes use of a single disease-specific template challenging. Ideally, one would like to use an average MRI template that is matched to the level of atrophy for the subject in question.

In this paper we propose a novel approach to generate a continuous four-dimensional template, where the 4 th dimension is a surrogate measure of overall brain atrophy, thus making it possible to generate the appropriate representative anatomical template from the wide range of possible levels of atrophy such that it is closely matched to the progression of the disease of the subject, while maintaining continuous one-to-one mapping to the common stereotaxic space.

Methods

In this study we used MRI scans obtained from the ADNI database (www.loni.ucla.edu/ADNI). This database contains 1.5T and 3T T1-w MRI scans of different populations such as Cognitively Normal (CN) subjects (n=233), patients with MCI (n=408) and patients with early AD (n=200) acquired at baseline. Automated methods to correct for image intensity non-uniformity [START_REF] Sled | A nonparametric method for automatic correction of intensity nonuniformity in MRI data[END_REF], normalize the image intensity range for 0-100, linearly register all available subjects to the MNI-152 stereotaxic space [START_REF] Collins | Automatic 3D intersubject registration of MR volumetric data in standardized Talairach space[END_REF] and estimate intracranial capacity (ICC) and lateral ventricles volume (LVV) [START_REF] Fonov | Robust automatic segmentation and characterization of lateral ventricle size in the ADNI cohort[END_REF] was applied to all available datasets. A human expert assessed quality of the resulting segmentations, and datasets failing any of the data processing steps were removed from processing. This QC process resulted in 173 NC, 150 AD and 332 MCI datasets that were available to participate in the 4D model creation.

The ratio between LVV and ICC (RLVV) was used as a surrogate measure of overall brain atrophy with mean (standard deviation) value of 2.46 (0.87)%. Subsets from all subjects (CN, MCI and AD) were selected with uniform distribution of RLVV from 1.0 to 6.0%, resulting in a total of 160 subjects.

Our algorithm for creating a Minimum Deformation Template [START_REF] Fonov | Unbiased average age-appropriate atlases for pediatric studies[END_REF] was modified to perform simultaneous 1) creation of the template and 2) linear regression of image intensity and shape versus RLVV. The output of the algorithm is a series of 3D volumes representing both linear model parameters of the intensity change and deformation at each voxel in the brain. These parameters can be used to generate an atlas for any value of RLVV.

Results

The ratio between LVV and ICC yielded values of mean(sd) 2.13(0.72)% for NC, 2.45(0.84)% for MCI and 2.84(0.91)% for AD.

The continuous, four dimensional anatomical template was created. For a given RLVV, an appropriate three dimensional anatomical template may be constructed, reflecting the average shape of the brain and the contrast between different tissue types for the given level of atrophy. Figure 1 shows transverse, sagittal and coronal images through 6 example values of increasing RLVV. As expected, CSF spaces in the lateral ventricles and sulcal spaces increase. Interestingly, the corpus callosum thins, and the superior and lateral surfaces of the brain appear to 'sag' down with increasing atrophy. The cerebellum atrophies as well.

Finally, it is worth noting that explicit, one-to-one mapping is calculated between templates corresponding to different RLVV values, making it possible to co-register MRI scans from subjects with very different levels of atrophy within same stereotaxic space, thus facilitating VBM studies.

Conclusions

The proposed method and resulting template will be useful tools for the development of robust automatic image processing methods targeted to the study of the populations with high degree of variability of atrophy. Furthermore, method presented is not limited to be used only with MCI and AD subjects, but can also be easily adopted for other neurodegenerative studies.


Pierrick Coupé

Vladimir Fonov

Simon F Eskildsen

José V Manjón

Douglas L Arnold

D Louis Collins

Influence of the Training Library Composition on a Patchbased label fusion method: Application to Hippocampus Segmentation on the ADNI dataset

Keywords: hippocampus, biomarkers, early detection, magnetic resonance imaging

Learning objectives:

1. To investigate the accuracy of a patch-based label fusion method for hippocampus segmentation on the ADNI database. 2. To evaluate the influence of the composition of the training library on the segmentation accuracy of a patch-based label fusion method.

Background:

The atrophy of medial temporal lobe structures is may serve as early biomarkers of Alzheimer's disease. The evaluation of hippocampus (HC) atrophy is estimated by volumetric studies requiring a segmentation step that can be very time consuming when done manually. This limitation can be overcome by using automatic segmentation methods. In this study, we propose to validate our nonlocal patch-based method [1] on the Alzheimer's Disease Neuroimaging Initiative (ADNI) database by segmenting the HC of Cognitively Normal (CN) subjects and patients with early Alzheimer's Disease (AD).

Method:

We used images obtained from the ADNI database. For our experiments, we randomly selected 10 MRI of CN and 10 MRI of AD. The first experiment was designed to evaluate the accuracy of our segmentation method on a group constituting of CN and AD (8 CN and 8 AD). The second experiment was designed to investigate the impact of the composition of the training library. Three different training populations were built (i.e., the CN population: 8 CN, the AD population: 8 AD, and the mixed AD / CN population: 4 CN and 4 AD). Through a leave-one-out procedure, our automatic segmentation method was applied on each of the 20 MRI scans. The quality of the obtained automatic segmentations was evaluated by estimating the Dice Kappa similarity index.

Results: For the first experiment, the median Dice Kappa values are presented in Table 1. The segmentation accuracy was significantly better (p-value=0.002) for CN (median k=0.883 for both HC) than for AD (median k=0.838 for both HC). For the second experiment, Figure 1 and Table 2 show the Dice Kappa similarity index for both studied populations according to the training library composition. For both populations, the best median Kappa values were obtained with the mixed CN / AD training library (k=0.875 for CN and k=0.835 AD).

Conclusion:

First, we demonstrated that our patch-based method provides high segmentation accuracy for both CN and AD populations. Second, we showed that the characteristic of the training library has a significant impact on the segmentation accuracy.

1.

Coupe

Figure 1 :

 1 Figure 1: Kappa index distribution according to the composition of the training library for both studied populations.

Table 1 .

 1 , P., et al., Neuroimage, 2011. 54(2): p. 940-54. Results of the first experiment using 16 templates in the training library.

	Median Kappa values	Left HC	Right HC	Both HC
	CN population	0.891	0.866	0.883
	AD population	0.830	0.858	0.838

Table 2 .

 2 Results of the second experiment using different training libraries. The symbol * is used to indicate a significant difference at 95% of confidence between a group and the two others.

	Median Kappa values	8 CN in training	4 AD / 4 CN in	8 AD in training
		library	training library	library
	CN Population	0.874	0.875	0.855*
	AD Population	0.824*	0.835	0.833


Pierrick Coupé

José V Manjón

Montserrat Robles

D Louis Collins

Adaptive Multiresolution Non-Local Means Filter for 3D MR Image Denoising

In this paper, an adaptive multiresolution version of the Blockwise Non-Local (NL-) means filter is presented for 3D Magnetic Resonance (MR) images. Based on an adaptive soft wavelet coefficient mixing, the proposed filter implicitly adapts the amount of denoising according to the spatial and frequency information contained in the image.

Two versions of the filter are described for Gaussian and Rician noise. Quantitative validation was carried out on Brainweb datasets by using several quality metrics. The results show that the proposed multiresolution filter obtained competitive performance compared to recently proposed Rician NL-means filters. Finally, qualitative experiments on anatomical and Diffusion-Weighted MR images show that the proposed filter efficiently removes noise while preserving fine structures in classical and very noisy cases. The impact of the proposed denoising method on fiber tracking is also presented on a HARDI dataset.

I. INTRODUCTION

In magnetic resonance (MR) imaging, denoising is an important preprocessing step used to improve image analysis. Since the noise directly impacts the accuracy of automatic approaches such as segmentation and registration, it is important to improve the Signal-Noise-Ratio (SNR) of images used during quantitative analysis. This is usually done by using a denoising method. However, the challenge is the remove of the noise component while preservation all the image features. In the past, many approaches have been proposed to address the difficult problem of MR image denoising [START_REF] Gerig | Nonlinear anisotropic filtering of MRI data[END_REF]- [START_REF] Awate | Feature-preserving MRI denoising: A nonparametric empirical bayes approach[END_REF].

Recently, in 2D natural image denoising, a new filter has been introduced by Buades et al. [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]: the Non-Local (NL-) means filter. Based on the natural redundancy of the images, this filter uses comparison of image patches in order to achieve the denoising. Due to its simplicity and high denoising performance, this filter is one of the best state-of-the-art methods. Recently, the NL-means or similar patch-based approaches are being adapted to MR denoising [START_REF] Coupé | An Optimized Blockwise NonLocal Means Denoising Filter for 3-D Magnetic Resonance Images[END_REF], [START_REF] Awate | Feature-preserving MRI denoising: A nonparametric empirical bayes approach[END_REF], [START_REF] Manjón | MRI denoising using Non-Local Means[END_REF], [START_REF] He | A Nonlocal Maximum Likelihood Estimation Method for Rician Noise Reduction in MR Images[END_REF]. This paper focuses on a new adaptive multiresolution version of the NL-means filter. As shown in [START_REF] Brox | Efficient nonlocal means for denoising of textural patterns[END_REF], there exists no set of parameters optimal for all the regions of a single image. According to the nature of the underlying structures (edge, texture, constant area...), the filtering parameters should be adapted. The main contribution of this paper is to introduce a new Adaptive Soft wavelet Coefficient Mixing (ASCM) procedure to implicitly optimize the denoising over all the space-frequency resolutions of the image. The main idea is to combine, in the wavelet transform domain, denoised images obtained with different sets of filtering parameters using the noisy input image.

II. METHODS

A. The Non Local Means filter

In the classical formulation of the NL-means filter [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], the restored intensity û(x i ) of the voxel x i , is a weighted average of the voxels intensities u(x j ) in the "search volume" V i of radius M :

û(x i ) = xj∈Vi w(x i , x j )u(x j ) (1)
where w(x i , x j ) is the weight assigned to intensity value u(x j ) to restore the intensity u(x i ) at voxel x i . More precisely, the weight evaluates the similarity between the intensity of the patches N i and N j centered on voxels x i and x j , such that w(x i , x j ) ∈ [0, 1] and xj∈Vi w(x i , x j ) = 1.

For each voxel x j in V i , the computation of the weight is based on the squared Euclidean distance between intensity patches u(N j ) and u(N i ), defined as:

w(x i , x j ) = 1 Z i e -u(N i )-u(N j ) 2 2 h 2 (2)
where Z i is a normalization constant ensuring that j w(x i , x j ) = 1, and h acts as a smoothing parameter. The main limitation of the standard NL-means filter is its computational burden, especially in 3D. Therefore, we used the blockwise approach and the patch preselection strategy proposed in [START_REF] Coupé | An Optimized Blockwise NonLocal Means Denoising Filter for 3-D Magnetic Resonance Images[END_REF] to drastically reduce the computational time. The blockwise approach consists in denoising the entire patch at the same time. By this way, we obtained restored values for all the elements of the patch:

û(N i ) = Nj∈Vi w(x i , x j )u(N j ), (3)
Since the patches overlap, several estimations are obtained for each voxel. The final restored intensity is thus obtained by averaging these different estimations. More details can be found in [START_REF] Coupé | An Optimized Blockwise NonLocal Means Denoising Filter for 3-D Magnetic Resonance Images[END_REF].

B. Adaptation to Rician noise

The raw data of an MR image are complex values corrupted by white additive Gaussian noise with the same variance in real and imaginary parts [START_REF] Gudbjartsson | The Rician distribution of noisy MRI data[END_REF]. In magnitude MR images, the distribution of the noise is transformed into a Rician distribution [START_REF] Gudbjartsson | The Rician distribution of noisy MRI data[END_REF], [START_REF] Sijbers | Estimation of the noise in magnitude MR images[END_REF]:

p(m) = m σ 2 exp(-
m 2 + A 2 2σ 2 )I 0 ( Am σ 2 ). ( 4
)
where σ is the standard deviation of Gaussian noise in the complex domain, A is the amplitude of the signal without noise, m is the value in the magnitude image and I 0 is the zeroth order modified Bessel function. The second-order moment of a Rician distribution can be written as:

E(m 2 ) = A 2 + 2σ 2 (5)
Recently, the nonlocal means filter has been adapted to Rician noise using Conventional Approach (CA) [START_REF] Manjón | MRI denoising using Non-Local Means[END_REF], [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio mri: Applications to DT-MRI[END_REF] and Maximum Likelihood Estimators (MLE) [START_REF] He | A Nonlocal Maximum Likelihood Estimation Method for Rician Noise Reduction in MR Images[END_REF].

In [START_REF] He | A Nonlocal Maximum Likelihood Estimation Method for Rician Noise Reduction in MR Images[END_REF], authors proposed a Non-local Maximum Likelihood (NLML) method considering that the N most similar pixels to the current pixel (in the patch distance sense) comes from the same Rician distribution. Knowing the level of noise, an optimization procedure is used to maximize the Rician-based log likelihood function:

logL = N i=1 log m i σ 2 - N i=1 m 2 i + A 2 2σ 2 + N i=1 log I 0 Am i σ 2 (6)
and the MLE is obtained as:

ÂMLE = arg max A (log L) (7)
The main limitation of the extention to 3D of these methods is the required computational load. Moreover, based on a minimization procedure, these approaches can converge to local minimum and require extra parameters.

In [START_REF] Manjón | MRI denoising using Non-Local Means[END_REF], [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio mri: Applications to DT-MRI[END_REF], authors have proposed Rician-adapted versions of the original NL-means filter using the properties of the second-order moment of Rician distribution called Conventional Approach (CA). The even order moments of a Rician law are usually used due to their simplicity. Indeed, by working on the squared signal, the bias intensity can be removed by subtraction as follows:

ÂCA = max E(m 2 ) -2σ 2 , 0 (8)
In denoising context, E(m 2 ) is usually estimated as the restored value of the squared noisy signal. The CA does not require extra computational time nor minimization procedure. Nevertheless, this method can artificially clamp background values to zero for very high levels of noise. Finally, the CA is more biased than the MLE when a small number of samples is used [START_REF] Sijbers | Estimation of the noise in magnitude MR images[END_REF].

As we will show in our experiments, the NL-means framework is better suited to CA approach since it involves a large number of samples during estimation. Therefore, we chosen to use the CA proposed in [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio mri: Applications to DT-MRI[END_REF] with the blockwise approach defined in [START_REF] Coupé | An Optimized Blockwise NonLocal Means Denoising Filter for 3-D Magnetic Resonance Images[END_REF] as follows:

û(N i ) = max     Nj∈Vi w(x i , x j )u(N j ) 2   -2σ 2 , 0   (9)
C. Multiresolution Framework 1) Hard Wavelet Subbands Mixing: As for all the denoising filters, the choice of the filtering parameters is crucial in NL-means-based restoration. Lukin [START_REF] Lukin | A multiresolution approach for improving quality of image denoising algorithms[END_REF] proposed a multiresolution framework to implicitly adapt the filtering parameters over the different space-frequency resolutions by combining in the transform domain the results obtained with different sets of filtering parameters. We have adapted this framework to NL-means-based restoration of 3D MR images in [START_REF] Coupé | 3D wavelet subbands mixing for image denoising[END_REF]. This Hard Subband Mixing (HSM) can be described as follows:

1) Denoise the original noisy MR image I n using two sets of filtering parameters: S u adapted to the image features preservation (i.e. small patches and/or low h parameter) and S o , adapted to the noise components removal (i.e. large patches and/or high h parameter). This yields two denoised images: Îu with efficient feature preservation (slightly undersmoothed) and Îo with efficient noise removal (slightly oversmoothed).

2) Decompose Îu and Îo into subbands with 3D discrete Wavelet Transform (DWT) at first level decomposition.

By using the usual notation: LLL denotes the low subband whereas LHH, LHL, LLH, HLL, HLH, HHL and HHH denote the high subbands

3) Select the LLL, LLH, LHL and HLL subbands of Îu and the LHH, HHL, HLH and HHH subbands of Îo .

4) Reconstruct the final image Î by an inverse 3D DWT from the combination of the selected subbands.

This procedure can be written as: ĉk = ĉk,u dk = dk,u for LLH, LHL and HLL dk = dk,o for HHL, HLH, LHH and HHH [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] where ĉk , ĉk,u are the scaling coefficients (i.e. coefficients in LLL low band) of Î and Îu decompositions, and dk , dk,u and dk,o are the wavelet coefficients in detailed subbands of Î, Îu and Îo decompositions (more details can be found in [START_REF] Coupé | 3D wavelet subbands mixing for image denoising[END_REF]). This HSM procedure improves the denoising performances of the NL-means filter by selecting the optimal denoised subbands along the image frequencies. Nevertheless, the HSM is not adapted according to noise level, nor to spatial information contained in subbands (i.e. the underlying structure). Moreover, with HSM, the subband selection is empirically optimized instead of being driven by the input signal.

2) Adaptive Soft Wavelet Coefficient Mixing:

In this paper, we propose to automatically adapt the mixing procedure according to the noise level and the spatial information contained in the detailed subbands. To achieve that, we use a similar strategy than in wavelet denoising thresholding method: detect wavelet coefficients corresponding to noise from coefficients corresponding to data by comparison with a threshold. Usually, the coefficients below this threshold are considered as noise and the others correspond to data information [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF]. In our case, the distance between coefficients of noisy input image and the threshold is used to drive the mixing of the coefficients of the denoised images.

First, we decompose I n , Îu and Îo into subbands with a 3D DWT. For each subband of I n , the signed distances between absolute values of noisy wavelet coefficients d k,n obtained from the original noisy image I n and a threshold are computed.

Then, this signed distance is used to perform a soft mixing between denoised wavelet coefficients, dk,u and dk,o , obtained from the two denoised images (see Fig. 1). When signed distance is superior to zero, the corresponding noisy wavelet coefficient is mainly considered as data, thus denoised wavelet coefficient from Îu should be privileged.

Inversely, when signed distance is inferior to zero, the corresponding noisy wavelet coefficient is mainly considered as noise, thus denoised wavelet coefficient from Îo should be privileged. In our approach, the weight assigned to each denoised wavelet coefficient is estimated using a weighting function as described below.

3) Adaptive Threshold: In the wavelet literature several thresholds have been proposed to perform denoising [START_REF] Donoho | Ideal spatial adaptation by wavelet shrinkage[END_REF], [START_REF] Chang | Adaptive wavelet thresholding for image denoising and compression[END_REF]. In our approach the BayesSkrink threshold [START_REF] Chang | Adaptive wavelet thresholding for image denoising and compression[END_REF] has been chosen. This threshold Tb can be computed for a wavelet coefficients subband b as:

Tb (σ) = σ 2 σX ( 11
)
where

σX = max(σ 2 b -σ 2 , 0) (12)
where σ2 b is the variance of the subband b.

4) Soft Coefficient Mixing:

It has been shown that the soft thresholding produces less artifacts and better visual results than hard thresholding [START_REF] Chang | Adaptive wavelet thresholding for image denoising and compression[END_REF]. Thus, we propose to use soft mixing instead of hard mixing such as the final coefficient is a weighted average of the coefficients from the two denoised images (see Fig. 1).

ĉk = ĉk,u dk = φ(d k , Tb (σ)) dk,u + (1 -φ(d k , Tb (σ))) dk,o (13)
where

φ(d k , Tb (σ)) = 1 1 + e (-λ(|dk|-Tb(σ))) (14)

III. MATERIALS

A. Phantom Database

In order to evaluate the performance of the different variants of the NL-means filter on 3D MR images, tests were performed on the BrainWeb database [START_REF] Collins | Design and construction of a realistic digital brain phantom[END_REF]. T1-w phantom MR images with different levels of noise were simulated to validate the performance of the denoising. For each image, ten levels of noise from 3% up to 21%

have been simulated for Gaussian and Rician distributions.

The Rician noise was built from white Gaussian noise in the complex domain:

• I r (x j ) = I(x j ) + η 1 (x j ), η 1 (x j ) ∽ N (0, σ) • I i (x j ) = η 2 (x j ), η 2 (x j ) ∽ N (0, σ)
where I is the noise free "ground truth" from Brainweb, I r is the real component and I i is the imaginary, and σ is the standard deviation of the added white Gaussian noise. Then, the noisy image is computed as:

I n (x j ) = I r (x j ) 2 + I i (x j ) 2 (15)

B. Filters compared

During the experiment, different implementations of the NL-means filter were compared:

• RNLM: voxelwise implementation of the Rician adapted NL-means filter proposed in [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio mri: Applications to DT-MRI[END_REF].

• ONLM/ORNLM: Optimized blockwise implementation of the NL-means filter proposed in [START_REF] Coupé | An Optimized Blockwise NonLocal Means Denoising Filter for 3-D Magnetic Resonance Images[END_REF] 1 without (ONLM) and with (ORNLM) the Rician adaptation described in [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio mri: Applications to DT-MRI[END_REF] .

• ONLM/ORNLM with HSM: ONLM or ORNLM filters with hard subbands mixing as proposed in [START_REF] Coupé | 3D wavelet subbands mixing for image denoising[END_REF].

• ONLM/ORNLM with ASCM: the proposed ONLM or ORNLM filters with adaptive soft coefficients mixing.

• NLML: 3D extention of the NonLocal-Maximum Likelihood approach proposed in [START_REF] He | A Nonlocal Maximum Likelihood Estimation Method for Rician Noise Reduction in MR Images[END_REF].

The selected filtering parameters for the different implementations were:

• For the RNLM filter: h 2 = σ 2 , α = 1, and M = 5 as proposed in [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio mri: Applications to DT-MRI[END_REF].

• For the ONLM/ORNLM filters: h 2 = σ 2 , α = 1, n = 2, M = 5, µ 1 = 0.95 and σ 2 1 = 0.5.

• For the ONLM/ORNLM with HSM/ASCM filters as proposed in [START_REF] Coupé | 3D wavelet subbands mixing for image denoising[END_REF]:

n = 2, S u = (α u , M w , h 2 ) = (1, 3, σ 2 )
for Îu and S o = (α o , M w , h 2 ) = (2, 3, σ 2 ) for Îo . For both we set µ 1 = 0.95 and σ 2 1 = 0.5.

• For the 3D extention of NLML filter: α = 1, N = 25 and M = 4. As in [START_REF] He | A Nonlocal Maximum Likelihood Estimation Method for Rician Noise Reduction in MR Images[END_REF], the number of pixels/voxels N used to maximize the likelihood function was set to 25.

As usually done in nonlocal means litterature, we use a semi-local paradigm instead of fully nonlocal paradigm by restricting the size of search volume. This enables to obtain a practical computational time but also to obtain near-optimal or optimal results excepted for image with repetitive textures [START_REF] Brox | Efficient nonlocal means for denoising of textural patterns[END_REF].

C. Quality metrics

The evaluation of denoising performance depends on the balance between noise removal and structure preservation. Each quality metric tends to privilege the first or second aspect, making it difficult to compare denoising filters. Nevertheless, an optimal denoising should improve the both at the same time. We used six different metrics to evaluate the noise removal and the structure preservation performance of the compared filters.

1) The Signal-to-Noise Ratio (SNR):

SN R = 10log 10 xi∈Ω I(x i ) 2 xi∈Ω (I(x i ) -Î(x i )) 2 ( 16
)
where Ω is the image grid. This metric is based on the noisefree "ground truth" image I and the denoised image Î. Higher SNR values correspond to good denoising results.

2) The MSE on gradient:

MSE on gradient = 1 |Ω| xi∈Ω (∆I(x i ) -∆ Î(x i )) 2 (17)
where the gradient modulus image is computed as: ∆I = ∆ x I 2 + ∆ y I 2 + ∆ z I 2 . This metric is based on the mean L 2 -norm between the first derivate of the images.

3) The Correlation Coefficient [START_REF] Sattar | Image enhancement based on a nonlinear multiscale method[END_REF]:

CoC = xi∈Ω (I(x i ) -Ī)( Î(x i ) -Ī) xi∈Ω (I(x i ) -Ī) 2 xi∈Ω ( Î(x i ) -Ī) 2 (18)
where Ī and Ī are the means of the noisefree and the denoised image.

4) The Edge Preservation Index [START_REF] Sattar | Image enhancement based on a nonlinear multiscale method[END_REF]:

EP I = xi∈Ω (∇I(x i ) -∇I)(∇ Î(x i ) -∇ Î) xi∈Ω (∇I(x i ) -∇I) 2 xi∈Ω ( ∇I(x i ) -∇ Î) 2 (19)
where the Laplacian operator ∇ is obtained with the standard 3 × 3 × 3 voxels approximation. The EPI is computed as the CoC between the second derivatives of the images.

5) The Universal Quality Index [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]:

U QI = σ I Î σ I σ Î . 2I Î (I) 2 + ( Î) 2 . 2σ I σ Î σ 2 I + σ 2 Î ( 20
)
where σ I Î corresponds to the covariance between the noisefree and the denoised image.

6) The Mean Structural SIMilarity (MSSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF]:

SSIM = (2I Î + C 1 )(2σ I Î + C 2 ) ((I) 2 + ( Î) 2 + C 1 )(σ 2 I + σ 2 Î + C 2 ) (21)
The constants C 1 and C 2 are used to stabilize the UQI when the means or the variances of the compared images are close to zero. During our experiment, the SSIM is locally computed within a 5 × 5 × 5 voxels Gaussian kernel. Finally, the MSSIM is estimated as the mean value of the local SSIM values.

IV. RESULTS

A. Comparison of CA and MLE

As in [START_REF] Sijbers | Estimation of the noise in magnitude MR images[END_REF], we propose to compare CA and MLE on a constant signal = 100 corrupted with several Rician noise levels. As in [START_REF] He | A Nonlocal Maximum Likelihood Estimation Method for Rician Noise Reduction in MR Images[END_REF], a MATLAB Nelder-Mead simplex method has been for optimization procedure during MLE estimation. Each experiment was repeated 10,000 times and the mean estimation is used as the result. In [START_REF] Sijbers | Estimation of the noise in magnitude MR images[END_REF] the size of the sample was N = 16 to be close to multiple acquisition number or neighborhood size used in classical denoising filters. Since NL-Means-based approaches use much larger numbers of samples (e.g. 11 3 = 1331 with a search area similar to [START_REF] Coupé | An Optimized Blockwise NonLocal Means Denoising Filter for 3-D Magnetic Resonance Images[END_REF]), we have studied the impact of the size of the sample on CA and MLE with N = 10, N = 100 and N = 1000. The results presented in Fig. 2 show that for small sample size (i.e. N = 10), the MLE obtained slightly better results than CA (around 1% at 100% noise). These results are similar to those provided by Sijbers in [START_REF] Sijbers | Estimation of the noise in magnitude MR images[END_REF]. However, when N increases, the difference between CA and MLE decreases. At N = 1000, CA and MLE provide similar estimation accuracy, while the computational time required by MLE is 1500 times higher than for CA. This computational time explained why the NL-means methods based on MLE use a limited number of samples (in [START_REF] He | A Nonlocal Maximum Likelihood Estimation Method for Rician Noise Reduction in MR Images[END_REF], N = 25 for NLML) and has been proposed only for 2D images. The presented experiment clearly demonstrates that it is better to use a CA with a high number of samples (N > 1000 in case of 3D NLM) than a MLE on a limited number of samples (N = 25 in case of NLML).

B. Impact of the ASCM on Gaussian Noise

We first studied the impact of the ASCM on Gaussian noise on a simple case (i.e. without the Rician bias intensity problem). Fig. 3 shows the results obtained for the different quality metrics. For all the metrics (except EPI) and for all the noise levels, the ASCM-based filter showed the best performances.

C. Comparison of Rician NL-means filters

Figure 4 presents the results obtained for the Rician adapted NL-means filters . For all the metrics and all the noise levels, the proposed ORNLM filter with ASCM obtained the best results. The ORNLM filter and the ORNLM filter with HSM obtained similar results. The empirical subband mixing proposed for Gaussian noise in [START_REF] Coupé | 3D wavelet subbands mixing for image denoising[END_REF] did not improve results on Rician noise, while, without any tuning, the ASCM automatically optimizes the denoising for both Gaussian and Rician noises.

Figure 5 shows the denoising results at 17% Rician noise. The voxelwise implementation of Rician NL-means filter (i.e. RNLM filter) obtained worse results than ORNLM-based filters due to soft structure degradation (see Fig. 5). As shown in [START_REF] Coupé | An Optimized Blockwise NonLocal Means Denoising Filter for 3-D Magnetic Resonance Images[END_REF], the preselection strategy better preserves contrast and edges. The 3D extention of NLML approach proposed in [START_REF] He | A Nonlocal Maximum Likelihood Estimation Method for Rician Noise Reduction in MR Images[END_REF] tended to leave more residual noise (see Fig. 5), thus adversely impacting the quality metrics.

Table I presents the computational time required by the different filters used during our comparison.

V. EXPERIMENTS ON CLINICAL DATASETS

A. T1-weighted anatomical MR image

The proposed method has been applied to a T1-w MR image of 170 × 256 × 256 voxels acquired on a 1.5T

Philips Gyroscan (Best, the Netherlands). Figure 6 shows the denoised image and the absolute residual image.

Visually, the ORNLM with ASCM preserves the anatomical structures in the image while efficiently removing the noise.

B. DW-MRI with high b-value

As assessed by quantitative results on phantom dataset, the proposed multiresolution approach is particularly efficient in the case of high levels of noise. To emphasize the improvement brought by using ASCM, the ORNLM Medical Systems, Erlangen, Germany) using an 8-channel phased-array head coil. The dataset was acquired with 99 diffusion encoding directions and 10 b=0 images, 2mm isotropic voxel size, 63 slices, b=3000 s/mm 2 , TE=121ms and TR=11.1s. As proposed in [START_REF] Aja-Fernandez | Restoration of dwi data using a rician lmmse estimator[END_REF], [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio mri: Applications to DT-MRI[END_REF], each b=0 image and each gradient direction has been denoised separately.

Figure 7 shows the results obtained on one of the gradient directions. The corresponding FA map is presented in Fig. 8. The FA maps computed on the denoised DW-MRI contained less noise and less artifacts than the FA map computed on the original data, especially on the frontal region (see green box in Fig. 8).

Finally, we conclude with an example of fiber tracking to show the qualitative performance of our denoising filter on a DTI application using the same data as the previous example. Streamline tractography was obtained by using the method described in [START_REF] Campbell | Validation and regularization in diffusion MRI tractography[END_REF]. The central slice of the corpus Callosum was used as seed region for the tractography (i.e. all the voxels in this region were used for the initialization).

Figure 9 shows the tractography obtained with the noisy data and the DW-MRI denoised with the proposed filter.

The fibres tracked on denoised image were longer and more dense.

Noisy Undersmoothed

Oversmoothed ORNLM with ASCM two denoised images, where the input noisy image is used to automatically drive the mixing. This new framework demonstrated better denoising capability using 6 quality metrics on Gaussian and Rician noise.

To deal with intensity bias introduced by Rician noise, Conventional Approach (CA) has been used within our NL-means framework as in [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio mri: Applications to DT-MRI[END_REF]. Experiment on bias estimation indicates that using CA for Rician noise removal within a 3D NL-Means framework appears to be the most efficient and practical. The denoising results on the BrainWeb MR phantom shows that the proposed Adaptive Multiresolution NL-means filter outperforms the optimized implementations presented in [START_REF] Coupé | An Optimized Blockwise NonLocal Means Denoising Filter for 3-D Magnetic Resonance Images[END_REF], [START_REF] Coupé | 3D wavelet subbands mixing for image denoising[END_REF] on Gaussian noise. Moreover, the Rician version of the proposed filter obtained competitive results compared to recently proposed Rician NL-means filters [START_REF] He | A Nonlocal Maximum Likelihood Estimation Method for Rician Noise Reduction in MR Images[END_REF], [START_REF] Wiest-Daesslé | Rician noise removal by non-local means filtering for low signal-to-noise ratio mri: Applications to DT-MRI[END_REF] in terms of denoising performance and computational time.

Fig. 1 .

 1 Fig. 1. Adaptive Soft Wavelet Coefficient Mixing Procedure. First, the wavelet transform of the original image is used to compute the threshold Tb for each subband b. Each noisy wavelet coefficient from original image is then compared to Tb in order to drive the soft mixing between the denoised wavelet coefficients. Finally, an invert wavelet transform is used to obtain final image.

Fig. 2 .

 2 Fig. 2. Comparison between the conventional and ML estimators.

Fig. 3 .

 3 Fig. 3. Impact of the ASCM on Gaussian noise. Results obtained with the different implementations of the ONLM filter on T1-w phantom.

Fig. 4 .

 4 Fig. 4. Comparison of the Rician-adapted NL-means filters on Rician noise.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Comparison of the Rician-adapted NL-means filters. Results obtained on T1-w phantom corrupted by 17% Rician noise. The images and the maps have the same scale for all the Rician filters compared.

Fig. 7 .Fig. 8 .

 78 Fig. 7. Gradient direction of the DWI dataset with high b-value . Top. Left: Noisy DW-MRI. Right: Undersmoothed Iu DW-MRI (SN R = 7.12dB estimated with [23]). Bottom. Left: Oversmoothed Io DW-MRI. Right: DW-MRI denoised with ORNLM filter with ASCM.

Fig. 9 .

 9 Fig. 9. Fibres seeded from the central slice of the Corpus Callosum.

This filter can be freely tested at: http://www.irisa.fr/visages/benchmarks

ACKNOWLEDGEMENTS

We would like to thank Ilana Leppert and Prof. Bruce Pike for providing DW-MRI data and help in result analysis.

We also want to thank Nicolas Guizard and Dr Vladimir Fonov for their help during tractography experiment. We thank Dr Luis Concha for his fruitful discussions about DTI processing. This work has been partially supported by the canadian grant Industry Cda (CECR)-Gevas-OE016 and by the Spanish Health Institute Carlos III through the RETICS Combiomed, RD07/0067/2001.