Van Minh Nguyen
email: vanminh.nguyen@alcatel-lucent.com

Chung Shue Chen
email: chung-shue.chen@inria.fr

Laurent Thomas
email: laurent.thomas@alcatel-lucent.com

Handover Measurement in Mobile Cellular Networks: Analysis and Applications to LTE

This paper deals with handover measurement in mobile cellular networks. The work is dedicated to network modeling and performance evaluation. The exposition focuses on neighbor cell scanning and addresses its key probabilistic events: (i) suitable handover target found, (ii) service failure, and (iii) scan withdrawal under the interference-limited condition in a multicell system. We derive their expressions and provide a generalized framework for the analysis of handover measurement failure and target cell quality by the best signal quality and minimum duration outage. Results applied to LTE have also shown its effectiveness.

I. INTRODUCTION

In mobile cellular networks, a user may travel across different cells during a service. Handover (HO) which switches the user's connection from one cell to another is therefore an essential function. Technology advancement is expected to minimize service interrupt and provide seamless handover [START_REF] Kim | Advanced handover schemes in IMT-advanced systems [WiMAX/LTE Update[END_REF]. A handover procedure includes two important functions which are handover measurement and handover decision-execution [START_REF] Anas | Performance analysis of handover measurements and layer 3 filtering for UTRAN LTE[END_REF]. The measurement function is responsible for monitoring the service quality from serving cell and finding a suitable neighbor cell for handover. This consists in measuring some parameters such as cell signal quality or strength. Handover decision-execution is made upon the measurement function. It decides whether a handover is needed and in that case coordinates multi-party handshaking among the user and cells to have handover execution fast and transparent.

Mobile-assisted network-controlled handover is highly recommended in emerging cellular network standards (e.g., LTE [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA) and Evolved Universal Terrestrial Radio Access Network (E-UTRAN) -Overall description: Stage 2[END_REF], LTE-Advanced and WirelessMAN-Advanced [START_REF] Kim | Advanced handover schemes in IMT-advanced systems [WiMAX/LTE Update[END_REF]) for its operational scalability and effectiveness in which the mobile measures the signal quality of its neighbor cells and reports the measurement result to the network to make handover decision. To indicate the quality of a neighbor cell, the mobile will synchronize with the cell and measure its signal quality which is best represented by the signal-to-interference-plusnoise ratio (SINR).

In this paper, we focus on handover measurement function. The primary goal of handover measurement is to find a suitable handover target before the signal quality from serving cell becomes unsustainable. A call or service will drop if the received SINR falls below a minimum level over a certain period such that the resulting bit error rate is unacceptable. In digital communications, the threshold corresponds to the point at which reliable error correction is no longer possible although a connection may not be torn down [START_REF] Veeravalli | A locally optimal handoff algorithm for cellular communications[END_REF]. The secondary goal is to at the same time find a high-quality handover target. The neighbor cell with best signal quality among those scanned in each measurement period is in general preferable.

Advanced wireless broadband systems such as 3G and 4G [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA) and Evolved Universal Terrestrial Radio Access Network (E-UTRAN) -Overall description: Stage 2[END_REF] allow adjacent cells operating in a common frequency band. The mobile can simultaneously scan several neighbor cells operating in the same frequency, which are thus called intra-frequency cells. The number of cells that the mobile can scan per measurement period is called mobile's measurement capability, say k. E.g., a LTE-compliant terminal will be able to scan eight intra-frequency cells in each period [5, §8].

It is obvious that higher measurement capability would enhance the performance of the handover measurement and thus that of handover, although this may require more advanced chipset and consume the mobile more resource. Therefore, some basic questions remain to be answered. Most importantly, what is the exact impact of k on the service failure probability and on the quality of target cell? How optimal is the current suggestion (e.g., k = 8 of LTE standard [START_REF]Evolved Universal Terrestrial Radio Access: Requirements for support of Radio Resource Management[END_REF])? Although handover is an important topic and has received a lot of studies, most of them addressed the handover control problem of the handover decision function. The handover measurement function has received less attention and most investigations and analysis are given by simulations. Note that while a handover control problem can be studied conventionally in a simplified model of two cells in which a handover decision is made by assigning the mobile to one of them, the handover measurement problem involves a more complex system in which the signal quality of best cell among a large number of cells needs to be determined. This often incurs modeling and analysis difficulty especially when stochastic parameters are introduced to better describe a wireless network. There lacks a generalized model and analytic study on the performance of the handover measurement which is essential to network design and optimization.

This paper establishes a framework for the system analysis and probabilistic aspect of handover measurement. By observations on stochastic geometry (see, e.g. [START_REF] Baccelli | Stochastic Geometry and Wireless Networks[END_REF], [START_REF] Nguyen | A stochastic geometry model for the best signal quality in a wireless network[END_REF]) and related work on minimum-duration outage [START_REF] Mandayam | Minimum duration outage for CDMA cellular systems: A level crossing analysis[END_REF]- [START_REF] Jiang | Cell-coverage estimation based on duration outage criterion for CDMA cellular systems[END_REF], it brings the following main contributions:

1) A generalized model is developed to characterize the handover measurement to facilitate network analysis.

2) The probability of the following key events are derived in a multicell system: a suitable handover target found, service failure, and scan withdrawal. Analytical results of the above are known fundamental but challenging.

3) By the above, we investigate the performance of LTE intra-frequency handover measurement in both (i) scanning failure probability and (ii) target cell quality.

II. SYSTEM DESCRIPTION

Poisson point process (p.p.) has received increasing interest [START_REF] Baccelli | Stochastic Geometry and Wireless Networks[END_REF] and is especially effective in modeling the spatial location of base stations (BSs) seen by a mobile whose location is often unknown to the network operator [START_REF] Win | A mathematical theory of network interference and its applications[END_REF]. The spatial distribution of BSs is modeled by a Poisson p.p. with intensity λ in R 2 .

In the following, we consider that a cell is covered by a BS with omni-directional antenna. The signal strength of BS i received by a mobile at a distance d i is given by:

P i = P tx d -β i Z i , ( 1
)
where P tx is the signal power of common pilot channel for handover measurement and is assumed constant [START_REF] Anas | Performance analysis of handover measurements and layer 3 filtering for UTRAN LTE[END_REF], β is the path loss exponent, d -β i accounts for the far-field path loss, and Z i = 10 Xi/10 denotes lognormal shadowing, where {X i } are assumed independent and identically distributed according to a Gaussian distribution with zero mean and standard deviation σ X < ∞.

In the time domain, we consider that X i (t) is stationary and has auto-correlation function R X (τ ) with a finite second derivative at the origin (see, e.g., [START_REF] Mandayam | Minimum duration outage for CDMA cellular systems: A level crossing analysis[END_REF]- [START_REF] Jiang | Cell-coverage estimation based on duration outage criterion for CDMA cellular systems[END_REF]). Within this context, fast fading is averaged out as it usually varies much faster than that of handover decision process [START_REF] Anas | Performance analysis of handover measurements and layer 3 filtering for UTRAN LTE[END_REF]. We will then focus on the above shadowing model. However, relaxing the assumption of lognormal Z i is in general possible to account for different fading scenarios. See, e.g., [START_REF] Nguyen | A stochastic geometry model for the best signal quality in a wireless network[END_REF] and [START_REF] Lai | Minimum duration outages in Rayleigh fading channels[END_REF].

The signal quality of a cell i given by SINR, denoted by Q i , is expressible as:

Q i = P i /(N 0 + I i ), (2)
where I i = ∑ j̸ =i P j , and N 0 is the thermal noise average power at the receiver and is assumed constant at all i.

III. HANDOVER MEASUREMENT

A. Overall Description

The details of handover measurement are described below. In principle, the mobile needs a time duration to successfully estimate the signal quality of a neighbor cell. Consequently, the signal quality of neighbor cells are only available at discrete moments. The time spacing between two consecutive moments is called a measurement period, denoted by T meas . E.g., in WCDMA and LTE, the mobile is required to measure eight intra-frequency cells within T meas = 200 ms [5, §8]. By contrast, the signal quality of the serving cell impacts the ongoing service in a time scale as short as one time slot which is usually much smaller than T meas . So, a service failure such as call drop or link outage may occur before the mobile finds a suitable handover target. Therefore, when studying neighbor cell scanning, the signal quality of the serving cell should be During each [m -1, m], if the signal quality of the serving cell becomes good enough, the mobile will withdraw the scanning. In case the SINR of the serving cell is too bad, the scanning will end in failure as a service interrupt occurs. If the scanning identifies a neighbor cell which has a signal quality better than a certain requirement γ req , the cell will be considered as handover target and the cell scanning thus ends in success. Otherwise, the mobile will continue the scanning and keep monitoring the signal quality of the serving BS.

B. Basic Formulation 1) Suitable Handover Target Found:

A neighbor cell is considered as a suitable handover target if its signal quality is better than the required threshold γ req . Note that when one needs to deal with a handover control problem, more criteria may be considered, e.g., whether the signal quality of the target cell is higher than that of the serving cell by a certain margin. Nevertheless, here we deal with the handover measurement function whose role is to find a suitable handover target and prevent service failure. Criteria for handover control problem in handover decision-execution are thus out of the scope here.

As a mobile will scan k cells in each measurement period and the cell with best signal quality is preferable, the event of having a suitable target at moment m can be defined by:

findtarget m (k) {Y k ≥ γ req } , (3)
where

Y k max i=1,••• ,k Q i , (4)
which denotes the best signal quality received from the k cells.

2) Service Failure: In wireless communication, the signal usually undergoes time-varying fading and other impairments such as interference. When the SINR stays below a minimum allowable level, say γ min , for a long time and successive bits are erroneous, redundancy coding does not help much, leading to a service failure. Therefore, it is more appropriate and also general to incorporate a minimum duration τ min when characterizing the event. A service failure during [m -1, m] is thus defined by the following minimum-duration outage:

fail m { Q 0 (t) < γ min , ∀ t ∈ [t 0 , t 0 + τ ] with τ ≥ τ min and t 0 + τ min ∈ [m -1, m] } , (5)
where Q 0 (t) denotes the SINR received from the serving cell at time t. Notice that when τ min = 0, (5) corresponds to an instantaneous SINR outage.

3) Scanning Withdrawal: Handover measurement introduces overheads such as time gaps in data transmission and mobile's resource consumption. It is thus helpful to withdraw scanning when the signal quality of serving cell becomes good enough. However, it is possible that the SINR may rise only for a very short duration and then drops below again. To ensure with high probability that the serving cell's signal quality will be good, the scanning should be canceled only if the serving cell's signal quality is higher than a threshold γ w for a certain period τ w , e.g. time-to-trigger used in LTE [START_REF]Evolved Universal Terrestrial Radio Access Radio Resource Control: Protocol specification[END_REF]. The event of scanning withdrawal during [m -1, m] is thus defined by:

m { Q 0 (t) ≥ γ w , ∀ t ∈ [t 0 , t 0 + τ ] with τ ≥ τ w and t 0 + τ w ∈ [m -1, m] } . ( 6
)

IV. ANALYTICAL RESULTS OF HO MEASUREMENT

A. Probability of Finding a Suitable Handover Target

To determine the probability of findtarget m (k), one needs to define the set of candidate cells from which k cells are taken. By today's standards, there are two different cases: limited candidate set and unlimited candidate set. In the former, a mobile only scans neighbor cells belonging to a pre-defined set which includes a certain number of potential candidates, say N cell cells. E.g., WCDMA and WiMAX networks use a neighbor cell list (NCL) with N cell = 32. In the latter, the mobile is allowed to scan any cell in the network. Note that a network may have a very large number of cells, scanning without restriction could introduce unsupportable overhead. Therefore, new broadband cellular systems use a set of, say N CSID , cell synchronization identities (CSID) which allow a mobile to identify and measure cells more easily. When required to scan k cells, the mobile just picks k out of the total N CSID CSIDs and conduct standardized cell synchronization and measurement. An example of this is LTE that defines 504 physical cell identities (PCI) [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA) and Evolved Universal Terrestrial Radio Access Network (E-UTRAN) -Overall description: Stage 2[END_REF], which serve as CSID. The mobile uses PCIs to scan neighbor cells autonomously [START_REF]Evolved Universal Terrestrial Radio Access Radio Resource Control: Protocol specification[END_REF].

Denote by B the network area of the candidate set. Then B is a bounded area in the case of limited candidate set, while B = R 2 in the case of unlimited candidate set. By [START_REF] Nguyen | A stochastic geometry model for the best signal quality in a wireless network[END_REF], [START_REF] Nguyen | Best signal quality in cellular networks: asymptotic properties and applications to mobility management in small cell networks[END_REF], we determine P{max k i=1 Q i ≥ γ req } for both cases.

1) Case B = R 2 : The set of neighbor cells which are scanned through a selection of k out of the totally N CSID CSIDs can be modeled as a thinning on B = R 2 , denoted by S k , with retention probability ρ k = k/N CSID [START_REF] Baccelli | Stochastic Geometry and Wireless Networks[END_REF]. Hence, (4) is equivalent to:

Y k ≡ max i∈S k Q i . ( 7
)
and

P{findtarget m (k)} = F Y k (γ req ), (8)
where

F Y k (•) is the tail distribution function of Y k .
By the result of best signal quality [7, Corollary 5], we can have:

P{findtarget m (k)} = F Y k (γ req ) = ∫ +∞ γreq ∫ +∞ 0 exp (-C 1 w α -C 2 (w, u)) π [ - 1 + γ req γ req × cos ( C 1 w α tan( πα 2 ) + C 3 (w, u) + C 4 (w, u) ) + cos ( C 1 w α tan( πα 2 ) + C 3 (w, u) -wu )] dwdu, (9)
where

α = 2/β, C 1 = c α (1 -ρ k )Γ(1 -α) cos( πα 2 ), C 2 (w, u) = ρ k c α u -α 1 F 2 (-α 2 ; 1 2 , 1 -α 2 ; -u 2 w 2 4 ), C 3 (w, u) = ρ k c α αw 1 -α 1 F 2 ( 1-α 2 ; 3 2 , 3-α 2 ; -u 2 w 2 4 ) u α-1 , C 4 (w, u) = w(1 -u(1 + γ)/γ), c α = πλ (P tx /N 0 ) α exp ( α 2 σ 2 Z /2 ) with σ Z = σ X log 10
10 , Γ(•) denotes the gamma function, and 1 F 2 denotes the hypergeometric function. Proofs are skipped due to a lack of space. Simulation result in Section V will also show its accuracy.

2) Case B R 2 : We consider two different possibilities: scattered networks like rural macro cellular networks where inter-site distance is large such that λ is small, and dense networks like urban small cell networks where λ is large. Consider B a disk-shaped network area with radius

R B = √ N cell /(πλ). ( 10
)
Under the Poisson p.p. assumption of BS spatial distribution, B has on average N cell BSs. In light of that, we can approximate B, the region of the N cell neighbor cells, by B.

For small λ and N cell is relatively large, we can have R B ≈ ∞, i.e., B can be approximated by R 2 . Similarly, let S k be a thinning on B with retention probability ρ k = k/N cell such that S k has on average k cells. The probability of finding a target cell can be well approximated by [START_REF] Graziosi | Distribution of outage intervals in macrodiversity cellular systems[END_REF].

For large λ (and so the approximation R B ≈ ∞ is not applicable), denoting the distance from the mobile to its closest neighbor BS by R min > 0, then the k cells scanned are located in a ring-shaped network area with inner and outer radii R min and R B . By [14, Theorem 2], we have:

P{findtarget m (k)} = F Y k (γ req ) ≈ ∫ ∞ γreq { f M k (u) ∫ ∞ 0 2 πw e -δw α sin ( w u -γ req 2γ req ) × cos ( wu + w u -γ req 2γ req -δw α tan πα 2 ) dw } du, (11)
where δ = c α Γ(1 -α) cos(πα/2), and f M k is the probability density function (pdf) of random variable

M k = max k i=1 P i : f M k (x) = k • f P (x) • F k-1 P (x), (12)
with F P denoting the cumulative density function (cdf) of the signal strength received from a cell in B is given by

F P (x) = c ( G 1 (x) a α - G 2 (x) b α - e ω G 3 (x) x α + e ω G 4 (x) x α
) ,

where a = (P tx /N 0 )R -β B , b = (P tx /N 0 )R -β min , c = (P tx /N 0 ) α /(R 2 B -R 2 min ), ω = 2σ 2 X /β 2 ,
and G j=1,2,3,4 refers to the cdf of lognormal distribution of parameters (µ j , σ X ), where µ 1 = log a, µ 2 = log b, µ 3 = µ 1 + ασ 2 X , and 11) is a good approximation [START_REF] Nguyen | Best signal quality in cellular networks: asymptotic properties and applications to mobility management in small cell networks[END_REF].

µ 4 = µ 2 + ασ 2 X . f P (x) = dF P (x)/dx. It is observed that (

B. Probability of Service Failure

Recall that X(t) is a Gaussian process for shadowing. By substituting (1) and ( 2), the event in ( 5) can be re-written as:

Q 0 (t) < γ min ⇔ X(t) < γmin (t), (13)
where

γmin (t) = 10 log 10 ( γ min (N 0 + I(t))d(t) β /P tx ) . (14)
By ( 13), the service failure probability P{fail m } is thus given by the minimum-duration outage of X(t) w.r.t. γmin (t) and τ min , considering [(m -1)T meas -τ min , mT meas ]. Note that τ min and T meas are typically about a few hundreds of milliseconds [START_REF]Evolved Universal Terrestrial Radio Access: Requirements for support of Radio Resource Management[END_REF]. The above interval, i.e., T meas + τ min , is very short so that we consider that in which the distance between the mobile and its serving BS is constant, denoted d m .

Following [START_REF] Mandayam | Minimum duration outage for CDMA cellular systems: A level crossing analysis[END_REF]- [START_REF] Jiang | Cell-coverage estimation based on duration outage criterion for CDMA cellular systems[END_REF], given that R X (τ ) has a finite second derivative at the origin, d(t) = d m , and if the interference is constant, say x 0 , we can have the following closed-form expression of the outage probability by the asymptotic properties of the level-crossing of a stationary Gaussian process:

P{fail m | I(t) = x 0 , d(t) = d m } = P{outage w.r.t. γmin } = η min [ τ min exp(- τ 2 min 2σ 2 min ) + √ 2πσ min q( τ min σ min ) ] , (15)
where

q(x) = ∫ ∞ x exp(-t 2 2 ) √ 2π
dt, σ min is the parameter of the Rayleigh distribution describing the fade duration and is given by:

σ min = √ 2/π q(-γ min /σ X )/η min , ( 16
)
and η min is the down-level crossing rate given by:

η min = √ ν/(2πσ X ) exp ( -γ 2 min /(2σ 2 X ) ) , ( 17
)
where ν = -d 2 dτ 2 R X (τ )| τ =0 . Note that in practice, the interference behaves as a random variable. By the stochastic model described in Section II, the interference field I(t) can be modeled as a shot noise on R 2 [START_REF] Baccelli | Stochastic Geometry and Wireless Networks[END_REF] and has the following characteristic function [START_REF] Nguyen | A stochastic geometry model for the best signal quality in a wireless network[END_REF], [START_REF] Win | A mathematical theory of network interference and its applications[END_REF]:

ϕ I (w) = exp ( -δ|w| α [1 -jsign(w) tan(πα/2)] ) , ( 18
)
where δ is defined under [START_REF] Win | A mathematical theory of network interference and its applications[END_REF]. Consider that β > 2, we have 0 < α < 1 such that ϕ I (w) is absolutely integrable. By [START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF]Thm.3,p.509], the pdf of I is thus given by:

f I (x) = 1 π ∫ ∞ 0 e -δw α cos ( δ tan( πα 2 )w α -xw ) dw. (19)
A more general expression of the service failure probability with time-varying random interference I(t) is expressible as:

P{fail m | d = d m } = ∫ ∞ 0 P{fail m |I = x, d = d m }f I (x)dx = ∫ ∞ 0 ∫ ∞ 0 η min π ( τ min e - τ 2 min 2σ 2 min + √ 2πσ min q( τ min σ min ) ) × e -δw α cos ( δ tan( πα 2 )w α -xw ) dwdx. (20)

C. Probability of Scanning Withdrawal

Similarly to ( 13) and ( 14), ( 6) can be re-written as:

Q 0 (t) ≥ γ w ⇔ -X(t) ≤ γw (t), (21)
where -X(t) is still a stationary Gaussian process, and

γw (t) = -10 log 10 ( γ w (N 0 + I(t))d(t) β /P tx ) . ( 22
)
Consequently, the scanning withdrawal event is given by the minimum-duration outage of -X(t) w.r.t. γw (t) and τ w . Since the interval T meas + τ w is also very short, we consider that during [(m -1)T meas -τ w , mT meas ], the distance between the mobile and its serving BS is constant. Provided that R X has a finite second derivative at the origin and

I(t) = x 0 , P{wdraw m |I(t) = x 0 , d(t) = d m } = P{outage w.r.t. γw } = η w [ τ w exp(- τ 2 w 2σ 2 w ) + √ 2πσ w q( τ w σ w ) ] , (23)
where

η w = √ ν 2πσX e -γ 2 w /(2σ 2 X ) and σ w = √ 2 π q(-γw/σX ) ηw
. Similarly to (20), by (19) we obtain a more general expression:

P{wdraw m | d(t) = d m } = ∫ ∞ 0 ∫ ∞ 0 η w π ( τ w e - τ 2 w 2σ 2 w + √ 2πσ w q( τ w σ w ) ) × e -δw α cos ( δ tan( πα 2 )w α -xw ) dwdx. (24)
V. APPLICATION TO LTE INTRA-FREQUENCY HO

In the following, we study the performance of LTE intrafrequency handover measurement by the above results. A mobile in LTE continuously scans intra-frequency cells during its connected mode. This corresponds to the case γ w = +∞ and so P{wdraw m } = 0. The mobile scans neighbor cells autonomously using 504 PCIs without a NCL [START_REF]Evolved Universal Terrestrial Radio Access (E-UTRA) and Evolved Universal Terrestrial Radio Access Network (E-UTRAN) -Overall description: Stage 2[END_REF], [START_REF]Evolved Universal Terrestrial Radio Access: Requirements for support of Radio Resource Management[END_REF]. P{findtarget m (k)} is thus given by ( 9) with retention probability ρ k = k/504. Fig. 2 describes the temporal evolution of the handover measurement. It is composed of Fail, Scan and TargetFound, which correspond to the states of service failure, scanning continuation, and a suitable handover target found, respectively. Note that when being in Fail and TargetFound, the current The transition probabilities from Scan at m -1 to states Fail, Scan, and TargetFound at m are denoted by f m , s m and g m , respectively. It is clear that they are given by:

f m = P{fail m }, ( 25
) g m = (1 -f m ) × P{findtarget m (k)}, (26)
s m = (1 -f m ) × (1 -P{findtarget m (k)}). (27)

A. Performance Metric

We investigate the system performance. As shown in Fig. 2, an occurrence of service failure at m is given by {Scan 1 , . . . , Scan m-1 , Fail m }. Thus,

P{scanning fails at m} = (s 1 × • • • × s m-1 ) × f m . (28)
The system's scanning effectiveness can be indicated by its asymptotic scanning failure probability expressible as:

F(k) = lim M →∞ M ∑ m=1 ( f m m-1 ∏ l=1 s l ) , ( 29
)
where M represents the time duration of the on-going service. Note that F(k) characterizes the possibility that the scanning procedure cannot identify a suitable target cell before a service failure occurs. However, finding target cell with high quality is as well desirable. Therefore, it is necessary to assess the target cell quality when the scanning is in success.

As a suitable target cell is given by the best cell among k cells scanned and provided that its signal quality is better than γ req , the resulting signal quality is thus given by E{Y k |Y k ≥ γ req }. Therefore, the tail distribution F Y k (y) is given by

F Y k |Y k ≥γreq (y) = 1(y < γ req ) + 1(y ≥ γ req ) F Y k (y) F Y k (γ req ) , ( 30
)
where 1(•) is the indicator function. Hence,

E{Y k |Y k ≥ γ req } = γ req + ∫ ∞ γreq F Y k (y)dy F Y k (γ req ) . ( 31
)
The scanning is in success when a suitable target cell is found, i.e., {Scan 1 , . . . , Scan m-1 , TargetFound m }. Thus, It follows that the quality of the scanning is measurable by the asymptotic target cell quality defined below:

P{Target found at m} = (s 1 × • • • × s m-1 ) × g m . (32)
Q(k) = lim M →∞ M ∑ m=1 P{Target found at m}E{Y k |Y k ≥ γ req },
(33) which, by (31) and (32), is expressible as:

Q(k) = ( lim M →∞ M ∑ m=1 g m m-1 ∏ l=1 s l ) ( γ req + ∫ ∞ γreq F Y k (y)dy F Y k (γ req ) ) . (34)

B. Simulation Setting

System parameters are summarized in Table I following 3GPP recommendations [START_REF]Evolved Universal Terrestrial Radio Access: Requirements for support of Radio Resource Management[END_REF], [START_REF]Evolved Universal Terrestrial Radio Access: Radio Frequency System Scenarios[END_REF] for LTE urban macro cellular networks. The network density λ is set corresponding to hexagonal cellular layout such that λ = 2/(3

√ 3R 2 ) BS/m 2 . In lognormal shadowing [8]-[10], R X (τ ) = σ 2 X exp(-(1/2) (vτ /d c ) 2 ), (35)
where v is the user's velocity. So, ν under (17) is given by: ν

= -R ′′ X (τ )| (τ =0) = σ 2 X (v/d c ) 2 . ( 36
)
We study relative and absolute requirement policies; both are considered in LTE [START_REF]Evolved Universal Terrestrial Radio Access Radio Resource Control: Protocol specification[END_REF]. In the former, γ req is required to be higher than γ min by a HO margin ∆ HO such that γ req = γ min + ∆ HO . In the latter, γ req is on a fixed absolute level.

Regarding user mobility, we considered 100 mobile users in the serving cell and each of them moves away from the serving BS at velocity v in a random direction generated according to a uniform distribution on (-π, π]. This scenario has been considered as the most critical circumstance [START_REF]Evolved Universal Terrestrial Radio Access: Radio Frequency System Scenarios[END_REF].

C. Numerical Result

First, Fig. 3 verifies our analytical model against computer simulation. Fig. 3(a) plots P{findtarget k } of ( 9), which corresponds to common LTE model. Its agreement with simulation illustrates the accuracy of modeling the best signal quality Y k defined in (4) by the maximum of SINRs received from the thinning S k proposed in [START_REF] Nguyen | A stochastic geometry model for the best signal quality in a wireless network[END_REF]. Fig. 3(b) checks the derived probability of minimum-duration outage P{fail m }. Results show that the analytical model is close to the simulation in both settings. The lack of space did not allow to show here comparisons of the other probabilities. However, we observed that analytical and simulation results are well matched.

Fig. 4 illustrates that high measurement capability reduces the scanning failure probability F(k). It also indicates F(k)'s dependence on γ min and γ req respectively. Comparing Fig. 4(b) and 4(a), absolute requirement policy shows better robustness than relative requirement policy to different system configurations. Note that the former uses a fixed and generally lower γ min . Fig. 4(b) shows that increasing k to the order of 10 2 improves the performance very marginally. Besides, we can see that the current LTE requirement of k = 8 is insufficient for reliable handover performance (here, no pre-defined NCL).

Fig. 5 shows the resulting target cell quality. In Fig. 5(a), lower γ min (and so lower γ req ) always results in higher target cell quality Q(k) under the relative requirement scheme. On the other hand, Fig. 5(b) shows a crossing point between the curves when k is in-between 10 and 16 under the absolute requirement scheme, in which given large k, Q(k) can be generally enhanced by higher γ req . However, for small k, it is better to maintain low γ req . Fig. 5 allows the network operator to set γ req parameter so as to obtain the best performance.

VI. CONCLUDING REMARKS

We focus on the analytical modeling of handover measurement function in emerging broadband 3G/4G networks with fundamental establishment in the probabilistic aspects of best signal quality and minimum-duration outage. We derive the probability of the key events in a multicell system: (i) suitable handover target found, (ii) service failure, and (iii) scan withdrawal. Particularly, (ii) and (iii) are derived by minimum-duration outage under interference-limited condition which is in fact more general and realistic and can cover the commonly used instantaneous SINR outage as special case. A generalized framework of HO measurement is provided for network design and optimization purposes. Results indicated that the LTE's current requirement for UE intra-frequency measurement capability with k = 8 seems insufficient for reliable HO measurement without using a neighbor cell list.

Fig. 1 .

 1 Fig. 1. Handover measurement and neighbor cell scanning

Fig. 2 .

 2 Fig. 2. Temporal evolution of the handover measurement

Fig. 3 .

 3 Fig. 3. Validation: analytical expression against simulation result

Fig. 4 .Fig. 5 .

 45 Fig. 4. Asymptotic scanning failure probability F (k)

ACKNOWLEDGMENT

This work was done within the framework of the joint laboratory of Alcatel-Lucent Bell Labs and INRIA. Authors would like to thank Franc ¸ois Baccelli for his valuable discussions and continuous supports in this work.


Christophe Collewet

Eric Marchand

Colorimetry-based visual servoing

. Indeed, as we did for the luminance, color attributes are directly used in the control law avoiding therefore any complex images processing as features extraction or matching. We propose in this paper several potential color features and then a way to select a priori the best choice among them with respect to the scene being observed. Experimental results validate as well the interest of using color attributes as visual features as our selection process.

I. INTRODUCTION

Visual servoing is a widely used technique in robot control [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. It is based on visual features extracted from a vision sensor. More precisely, the control law is designed so that the visual features s extracted from the current image at the current pose r, reach a desired value s * acquired at the desired pose r * leading to a correct realization of the task. The control principle relies on the regulation to zero of the following error vector e = ss * . To design this control law, the knowledge of the interaction matrix L s is usually required. For eye-in-hand systems, it links the time variation of s to the camera instantaneous velocity v:

ṡ = L s v (1)
with v = (v, ω) where v is the linear velocity and ω is the angular velocity. In this approach, the choice of the visual features as well as its related interaction matrix play an essential role [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. Most often, geometric visual features as points, straight lines, pose or homography are chosen to control the robot. However, we recently introduced non geometric visual features [START_REF] Collewet | Visual servoing set free from image processing[END_REF], [START_REF] Collewet | Photometry-based visual servoing using light reflexion models[END_REF]. Indeed, this new approach is very interesting since it strongly limits the image processing. More precisely, it avoids the robust extraction, matching (between x(r k ) and x(r * ), where x denotes the visual measurements required to design and compute s) and real-time spatio-temporal tracking (between x(r k-1 ) and x(r k )) which are well known to be difficult tasks. Supplementary advantages are that this approach is not sensitive to partial occlusions and to coarse approximations of the depths required to compute the interaction matrix. In these papers, instead of using geometric visual features we directly used the luminance. In that case, the visual feature vector s is nothing but the image itself while s * is the desired image. The error e is then simply the difference between the current and desired images. The main issue of these works C. Collewet and E. Marchand are with INRIA Rennes -Bretagne Atlantique, IRISA, Lagadic team, Rennes, France. firstname.name@inria.fr.

was to exhibit the analytical form of the interaction matrix related to the luminance. This was therefore a very different approach than previous works based on luminance [START_REF] Deguchi | A direct interpretation of dynamic images with camera and object motions for vision guided robot control[END_REF]- [START_REF] Silveira | Direct visual servoing with respect to rigid objects[END_REF] (see [START_REF] Collewet | Visual servoing set free from image processing[END_REF] for more details).

In this paper we investigate the use of color images instead of using the luminance as we did before. Indeed, using color attributes are more discriminant than using simply the luminance since two image points with a same luminance value can be differentiated from their color attributes. That is why color is widely used in image processing as for example in image segmentation. In contrast, to our knowledge, color has not been used in visual servoing. Of course, it has been used in robot control to segment the object of interest and then to apply the control law (see for example [START_REF] Cubber | Color-based visual servoing under varying illumination conditions[END_REF]). However, let us point out that our approach is radically different. Our goal is to use features based on the R, G and B components of the image directly in the control law. As we did for the luminance, we also here want to avoid any complex image processing as features extraction and matching. Another contribution of this work is to propose a criterion able to select good attributes from a set of interesting potential color attributes. As we shall see, this way to proceed allows to greatly improves the behavior of the control law, especially the 3D trajectory of the camera.

This paper is organized as follows. First, we recall in Section II the way we used the luminance in our previous works [START_REF] Collewet | Visual servoing set free from image processing[END_REF]. Then, we present in Section III potential color attributes and compute their related interaction matrix. In Section IV, a way to select them is discussed. Experimental results on positioning tasks are shown in Section V.

II. LUMINANCE AS A VISUAL FEATURE

To make the paper more readable, before investigated color visual features, we first recall how the luminance has been used in [START_REF] Collewet | Visual servoing set free from image processing[END_REF].

In this work, we have considered the luminance I of each point of the image. Thus, we have

s(r) = I(r) = (I 1• , I 2• , • • • , I N • ) (2)
where I k• is the k-th line of the image. I(r) is then a vector of size N × M where N × M is the size of the image.

The main problem with such visual features is, of course, to derive the analytical form of the interaction matrix related to the luminance. It has been performed in [START_REF] Collewet | Visual servoing set free from image processing[END_REF] for Lambertian scenes while the general case has been performed in [START_REF] Collewet | Photometry-based visual servoing using light reflexion models[END_REF]. However, for the sake of clarity, we only present here the former case.

The basic hypothesis assumes the temporal constancy of the brightness for a physical point between two successive images. Therefore, assuming that the image point x has a displacement dx in the time interval dt, the previous hypothesis leads to

I(x + dx, t + dt) = I(x, t). (3)
Written with a differential form, a first order Taylor series expansion of this equation around x gives

∇I ⊤ ẋ + İ = 0. (4)
with İ = ∂I/∂t. It becomes then straightforward to compute the interaction matrix L I related to I by plugging the interaction matrix L x related to x into (4). We obtain

İ = -∇I ⊤ L x v. (5)
Finally, if we introduce the interaction matrices L x and L y related to the coordinates x and y of x, we obtain

L I = -(∇I x L x + ∇I y L y ) (6)
where ∇I x and ∇I y are the components along x and y of the image gradient ∇I. Let us point out that it is actually the only image processing step necessary to implement this approach.

III. POTENTIAL COLOR VISUAL FEATURES

Instead of considering the luminance, we consider here color attributes. Similarly to the previous section, the simpler visual features could be one of the three R, G, B components of the color vector C = (R, G, B). Note that another choices could be any combinations of 2 of these color planes or even all the three components. In these cases, the interaction matrix is simply the one obtained by stacking all the interaction matrices related to the planes that have been chosen. Note also that the interaction matrix related to one color plane is obtained from [START_REF] Kallem | Kernelbased visual servoing[END_REF] when I has to be substituted by R, G, or B. However, more complex visual features can be used as non-linear functions of C. Interesting non-linear functions are color invariants [START_REF] Healey | Segmenting images using normalized color[END_REF]. Indeed, such functions are invariant to any changes in the scene geometry (that is the relative positions observer / lighting source / surface being observed). They are also invariant to intensity changes of the lighting.

Assuming a white lighting, two classes of color invariants exist: invariants dedicated to Lambertian materials and invariants dedicated to specular materials (see for example [START_REF] Gevers | Object recognition based on photometric color invariants[END_REF]). For Lambertian materials, the following invariants are defined [START_REF] Healey | Segmenting images using normalized color[END_REF]:

L ϕ i = ϕ N i (7)
where ϕ ∈ {R, G, B} and i = 1 or 2. N 1 is defined as

N 1 = R + G + B ( 8
)
while N 2 is defined as

N 2 = R 2 + G 2 + B 2 . ( 9
)
This invariants are known as L 1 or L 2 invariants since they performed a L 1 or L 2 normalization of the (color) image.

For specular materials, the following invariants, known as l 1 l 2 l 3 , are used [START_REF] Gevers | Object recognition based on photometric color invariants[END_REF]:

   l 1 = (R -G) 2 /N 3 l 2 = (R -B) 2 /N 3 l 3 = (G -B) 2 /N 3 (10)
where

N 3 = (R -G) 2 +(R -B) 2 +(G-B) 2 .
Note that this invariant is also normalized. However, this latter invariant may be very noisy or even undefined at any point where the color is not saturated1 That is for the set of points belonging to the achromatic axis (the grey points). In contrast, this problem occurs with the L 1 invariants only for R ≈ G ≈ B ≈ 0. That is for black points only.

In addition, let us point out that the components of all of these invariants [START_REF] Silveira | Direct visual servoing with respect to rigid objects[END_REF] or [START_REF] Gevers | Object recognition based on photometric color invariants[END_REF] are not independent, given two of them the third one can be computed; therefore two components are sufficient to describe completely the color at a point x. That is why we never combine all the components of a normalized invariant.

derive a control law directly based on one of a component of ( 7) or [START_REF] Gevers | Object recognition based on photometric color invariants[END_REF], or on different combinations of two of their components, the related interaction matrix has to be determined. This computation is, however, simple. Considering a component I inv of an invariant given by ( 7) or [START_REF] Gevers | Object recognition based on photometric color invariants[END_REF], the temporal luminance constancy hypothesis (3) holds, leading similarly as [START_REF] Kallem | Kernelbased visual servoing[END_REF] to

L Iinv = -∇I inv ⊤ L x . (11)
Note that the spatial gradient of I inv can be easily computed from [START_REF] Silveira | Direct visual servoing with respect to rigid objects[END_REF] or [START_REF] Gevers | Object recognition based on photometric color invariants[END_REF] from the spatial gradients of R, G and B.

For example, considering [START_REF] Silveira | Direct visual servoing with respect to rigid objects[END_REF], we simply have

L L ϕ i = - 1 N i (∇ϕ -L ϕ i ∇N i ) ⊤ L x . (12)
Recall that to remove the temporal luminance constancy hypothesis the approach described in [START_REF] Collewet | Photometry-based visual servoing using light reflexion models[END_REF] has to be used leading, of course, to a different and more complicated expression than [START_REF] Feddema | Automatic selection of image features for visual servoing of a robot manipulator[END_REF].

A more complicated problem is rather to choose the best color attributes to use from a given color image of the scene being observed. For example, let us consider the simple case where an image only contains the blue and green colors, it is obvious that using the R color plane as the only visual features will lead to a bad behavior of the control law. This problem is the subject of the next section.

IV. SELECTION OF VISUAL FEATURES

The answer to that question must be seen from the visual servoing problem, that is which color attributes will lead to a better behavior of the control law ? that essentially means a convergence from a larger domain and a better 3D camera trajectory. Usually, the condition number (CN) of the interaction matrix is used as a criterion able to compare visual features [START_REF] Feddema | Automatic selection of image features for visual servoing of a robot manipulator[END_REF], [START_REF] Chaumette | Image moments: a general and useful set of features for visual servoing[END_REF]. Using the L 2 norm, it writes as

CN s = max σ i min σ i (13)
with i = 1, • • • , 6 for a 6 ddl robot. The vector σ is obtained from the singular value decomposition of the interaction matrix related to s. A low value of this criterion leads to a better behavior of the closed-loop system with respect to the noise. It means that for small changes in the pose or in the interaction matrix, changes in the visual features remain small too. Let us introduce another criterion. To do that, contrary to the seminal paper [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF] where the control law was designed so that a desired behavior was required for e, we propose to see the design of the control law as an optimization problem as proposed in [START_REF] Collewet | Visual servoing set free from image processing[END_REF], [START_REF] Miura | Robust and uncalibrated visual servoing without Jacobian using a simplex method[END_REF], [START_REF] Malis | Improving vision-based control using efficient secondorder minimization techniques[END_REF]. In this case, the goal is to minimize the following cost function (see Section I).

C(r) = 1 2 e(r) 2 . ( 14
)
Of course, the graal would be to find some visual features so that C(r) would be an hypersphere. Indeed, only a global minimum would exist and a simple steepest descent method would ensure to this minimum. Therefore, we are interested in studying the Hessian of ( 14). It is given by

∇ 2 C(r) = ∂s ∂r ⊤ ∂s ∂r + i=dim s i=1 ∇ 2 s i s i (r) -s i (r * ) . (15
) However, this expression is far too complex to derive some useful results. Thus, we study it around the desired position r * , leading to

∇ 2 C(r * ) = ∂s ∂r ⊤ ∂s ∂r . (16)
Moreover, since we have ṡ = ∂s ∂r ṙ = L s v, we are interested in practice in the following matrix

H * = L s * ⊤ L s * . (17)
This matrix allows us to estimate the cost function around r * . Indeed, a first order Taylor series expansion of the visual features s(r) around r * gives

s(r) = s(r * ) + L s * ∆r (18)
where ∆r denotes the relative pose between r and r * . Therefore, by plugging (18) into ( 14), we obtain an approximation of the cost function in a neighborhood of r *

C(r) = 1 2 ∆r ⊤ H * ∆r. (19)
Consequently, another criterion to select visual features is to choose them so that the eigenvalues of H * are the most similar as possible, therefore we introduce the following criterion

c s = max λ i min λ i ( 20
)
where the λ i 's are the eigenvalues of H * . Of course, the optimal choice for s occurs when c s = 1 since in that case the cost function is an hypersphere. In contrast, if c s is a high value, then the shape of the cost function presents a narrow valley, which is well known to be a complex optimization problem. Therefore, dedicated control laws have to be designed as the one we have proposed in [START_REF] Collewet | Visual servoing set free from image processing[END_REF]. The problem is now to compare the criterion given by ( 13) and the one given by (20). First, let us write the interaction matrix from its singular value decomposition:

L s = UΣV ⊤ (21)
where Σ is a diagonal matrix which diagonal is σ. By multiplying each side of (21) by L s * ⊤ , this relation becomes

L s * ⊤ L s * = VΣ 2 V ⊤ . (22)
However, since L s * ⊤ L s * is a symmetric matrix, it is diagonalizable as

L s * ⊤ L s = WΛW ⊤ . ( 23
)
where W is an orthonormal matrix. Consequently, we have Λ = Σ 2 leading finally from the definitions of CN s and c s to

CN s = √ c s . (24)
This last relation means that the condition number is not only an interesting criterion related to robustness and sensitivity to noise but also related to the shape of the cost function at the desired position.

V. EXPERIMENTAL RESULTS

In all the experiments reported here, the camera is mounted on a 6 degrees of freedom gantry robot. Control law is computed on a Core 2 Duo 3Gz PC running Linux. Image are acquired at 66Hz using an IEEE 1394 camera with a resolution of 320 × 240. Moreover, we used specular and planar objects (some posters) as shown on Fig. 1. The experiments consist in achieving positioning tasks by visual servoing using the following control law

v = -λ (H + µdiag(H)) -1 L s ⊤ (s(r) -s(r * )) (25)
with H = L s ⊤ L s and where λ and µ are positive scalars (see [START_REF] Collewet | Visual servoing set free from image processing[END_REF] where the way to set these parameters is described).

For all the experiments the desired positions was so that objects and CCD planes are parallel. We compared the use of the following different sets of visual features from the values of the criterion (24): set 1: the color planes R, G, B or the luminance; set 2: all the combinations of the color planes R, G or B; set 3: one component of each invariants; set 4: the combinations of two components of a given invariant. For the first experiment, the initial pose was r init = (14 cm, -18 cm, 2.5 cm, -16 • , -10 • , -1.3 • ). The initial and desired images are reported respectively on Fig. 2a and Fig. 1a. The worse results have been obtained when using the l 1 l 2 l 3 invariant as well for the set of visual features 3 as for the set of visual features 4. In those cases, the condition numbers were very bad, between 195 and 274. Therefore, either the control law diverged or converged very slowly. The control law also diverged when using the luminance (the condition number was 148). Better results were obtained using the sets of visual features 1 or 2 (the condition numbers were around 140). However, the best results were obtained for the sets of visual features 3 or 4 when using either the L 1 or the L 2 invariants as can be seen on Fig. 3. More precisely, Fig. 3a depicts the behavior of the different normalized cost functions (a cost function evaluated at time k divided by the cost function evaluated at the first frame) while Fig. 3b depicts the trajectories of the camera in a fixed frame. As can be seen using a L 1 or L 2 invariant leads to a better trajectory than when using visual features from the set 1, it also leads to better decreases of cost functions. Fig. 3 also shows that the differences between the L 1 and L 2 invariants are very low (for both the condition numbers were very similar, around 105). These results validate that a bad condition number leads to a bad behavior of the convergence of the control law (a slow convergence or even a divergence) and to a bad 3D trajectory of the camera. This also validates that the condition number is also related to the form of the cost function.

The second experiment still concerns the same object, only the initial pose has been changed. Thus, the desired image is still the one reported on Fig. 1a. The initial image is represented on Fig. 2b. For this experiment we had r init = (4.3 cm, -18 cm, 0 cm, -16 • , 2.5 • , -12.5 • ). As during the first experiment, using visual features from first and second sets of visual features leads to either a divergence, a slow convergence or a bad camera trajectory. For example, using the R or B color plane in the first experiment yields a very slow convergence (not described in Fig. 3) while in the second experiment their use leads to a divergence. It is also the case for the G color plane where as well during the first as during the second experiment a bad behavior have been obtained. Besides, here again using the luminance failed. Moreover, as during the first experiment, the best results have been obtained using either the invariants L 1 or L 2 (those based on l 1 l 2 l 3 failed); in addition, since the values of their criterion are similar, a similar behavior with respect to the convergence of the control law or with respect to the camera trajectory is observed. Fig. 5 reports the L 1 and L 2 invariants at the desired position for the experiments 1 and 2. As can be seen, the differences between a L ϕ 1 and L ϕ 2 are small. Moreover, except considering ϕ = B, these images are quite similar. Indeed, for ϕ = B, the images are dark which means that these invariants are not discriminant at all. That explains why

. (b) L R 2 . (c) L G 1 . (d) L G 2 . (e) L B 1 . (f) L B 2 .
we obtained among the L ϕ i invariants the worse value for the criterion (120). For ϕ = B, the control law has converged very slowly during the first experiment and has diverged during the second one.

The third experiment concerned a different object. The initial and desired images are reported respectively on Fig. 6a and Fig. 1b. The initial pose was the same as for the first experiment. As for the previous experiments, the l 1 l 2 l 3 invariants have very bad values for their criterion (around 114) and consequently have failed. The sets of visual features 1 and 2, and the luminance, have also some bad values, but slight better (around 109). Thus, as can be seen on Fig. 7, using the R color plane or the luminance have led to a slow convergence and to a bad camera trajectory, but they have converged. As previously, Fig. 8 reports the L 1 and L 2 invariants at the desired position for the experiments 3 and 4. As can be seen, these images are much more different than for the first and second experiments, especially for L R 1 and L R 2 , the L R 2 image is reacher than the L R 1 one. It is also confirmed by the value of the criterion, respectively 105 and 96. Therefore, as expected, using L R 1 leads to a divergence while L R 2 leads to a convergence. On the other hand, as can be observed on Fig. 8c andd, these images are quite homogeneous and thus these invariants are not interesting (recall that the interaction matrix requires non null values for the spatial image gradients, see [START_REF] Feddema | Automatic selection of image features for visual servoing of a robot manipulator[END_REF]). This is validated by bad criterion values (around 107), thus these invariants have led to a divergence. In contrast, using L B i has led to a convergence (the criterion value was 103).

The fourth experiment concerned the same object as for the third experiment but the initial pose was the one of the second experiment. This experiment validates the results of the third experiment, the visual features with a bad criterion value have also led to a divergence or to a bad camera trajectory during this experiment (see Fig. 9).

Remark 1:

The reader may wonder why we do not present any experiment concerning robustness against changes in the scene geometry as mentioned in section III. In fact, the L 1 and L 2 invariants only perform on true Lambertian materials (such as cloths for example) and are not robust at all (concerning changes in the geometry) if the object is not Lambertian as the objects we used in our experiments. In addition, we have tried to use invariants to specular objects (the l 1 l 2 l 3 ones), but these invariants have a very bad criterion value and, consequently, have always led to a bad behavior of the control law. Remark 2: Even if none result are given here, as in [START_REF] Collewet | Visual servoing set free from image processing[END_REF], using color attributes with a low criterion value also leads to a very low positioning error (typically 0.1 mm for the translations and 0.01 • for the rotations). It is because ss * is very sensitive to the pose r. Such visual features also lead to an approach that is not sensitive to partial occlusions (this nice behavior is due to the high redundancy of the visual features we use) and also not sensitive to coarse approximations of the depths required to compute the interaction matrix.

VI. CONCLUSION

It has been shown it this paper that using color attributes greatly improved simply the use of the luminance. As in the luminance case, the complex problems of features extraction and matching are avoided. However, this paper have also shown that a selection of these color attributes is required since their values highly influence the convergence of the control law and the 3D camera trajectory. Especially, the visual features obtained from the L 1 and L 2 invariants have led to the best behaviors.

Fig. 1 .

 1 Fig. 1. Desired image for the objects used in the experiments. (a) First and second experiment. (b) Third and fourth experiment.

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Initial images. (a) For the first experiment. (b) For the second experiment.

Fig. 5 .

 5 Fig. 5. Images of the L ϕ i invariants for the first object. (a) L R 1. (b) L R 2 . (c) L G 1 . (d) L G 2 . (e) L B 1 . (f) L B 2 .

Fig. 6 .Fig. 7 .

 67 Fig. 6. Initial images. (a) For the third experiment. (b) For the fourth experiment.

Fig. 8 .

 8 Fig. 8. Images of the L ϕ i invariants for the second object. (a) L R 1. (b) L R 2 . (c) L G 1 . (d) L G 2 . (e) L B 1 . (f) L B 2 .

Fig. 9 .

 9 Fig. 9. Fourth experiment : comparison of potential color visual features.

The color at a given image point is not saturated if R ≈ G ≈ B.


Igor Yanovsky

Carole Le Guyader

Alex Leow

Arthur Toga

Paul Thompson

Luminita Vese

Unbiased Volumetric Registration via Nonlinear Elastic Regularization

In this paper, we propose a new large-deformation nonlinear image registration model in three dimensions, based on nonlinear elastic regularization and unbiased registration. Both the nonlinear elastic and the unbiased functionals are simplified introducing, in the modeling, a second unknown that mimics the Jacobian matrix of the displacement vector field, reducing the minimization to involve linear differential equations. In contrast to recently proposed unbiased fluid registration method, the new model is written in a unified variational form and is minimized using gradient descent on the corresponding Euler-Lagrange equations. As a result, the new unbiased nonlinear elasticity model is computationally more efficient and easier to implement than the unbiased fluid registration. The model was tested using three-dimensional serial MRI images and shown to have some advantages for computational neuroimaging.

Introduction

Given two images, the source and target, the goal of image registration is to find an optimal diffeomorphic spatial transformation such that the deformed source image is aligned with the target image. In the case of non-parametric registration methods (the class of methods we are interested in), the problem can be phrased as a functional minimization problem whose unknown is the displacement vector field u. Usually, the devised functional consists of a distance measure (intensity-based, correlation-based, mutual-information based [START_REF] Modersitzki | Numerical Methods for Image Registration (Numerical Mathematics and Scientific Computation)[END_REF] or metric-structure-comparison based [START_REF] Lord | Simultaneous registration and parcellation of bilateral hippocampal surface pairs for local asymmetry quantification[END_REF]) and a regularizer that guarantees smoothness of the displacement vector field. Several regularizers have been investigated (see Part II of [START_REF] Modersitzki | Numerical Methods for Image Registration (Numerical Mathematics and Scientific Computation)[END_REF] for a review). Generally, physical arguments motivate the selection of the regularizer. Among those currently used is the linear elasticity smoother first introduced by Broit [3]. The objects to be registered are considered to be observations of the same elastic body at two different times, before and after being subjected to a deformation as mentioned in [START_REF] Modersitzki | Numerical Methods for Image Registration (Numerical Mathematics and Scientific Computation)[END_REF]. The smoother, in this case, is the linearized elastic potential of the displacement vector field. However, this model is unsuitable for problems involving large-magnitude deformations.

In [4], Christensen et al. proposed a viscous fluid model to overcome this issue. The deforming image is considered to be embedded in viscous fluid whose motion is governed by Navier-Stokes equations for conservation of momentum: µ△v(x, t) + (ν + µ)∇(∇ • v(x, t)) = f(x, u(x, t)),

(1) v(x, t) = u t (x, t) + ∇u(x, t) • v(x, t).

(2)

Here, equation [START_REF] Lord | Simultaneous registration and parcellation of bilateral hippocampal surface pairs for local asymmetry quantification[END_REF], defining material derivative of u, nonlinearly relates the velocity and displacement vector fields. One drawback of this method is the computational cost. Numerically, the imagederived force field f(x, u(x, t)) is first computed at time t. Fixing the force field f, linear equation ( 1) is solved for v(x, t) numerically using the successive over-relaxation (SOR) scheme. Then, an explicit Euler scheme is used to advance u in time. Recent works [5-7] applied Riemannian nonlinear elasticity priors to deformation velocity fields. These alternating frameworks, however, are time-consuming, which motivates the search for faster implementations (see for instance [START_REF] Bro-Nielsen | Fast fluid registration of medical images[END_REF] or [START_REF] D'agostino | A viscous fluid model for multimodal non-rigid image registration using mutual information[END_REF] in which the instantaneous velocity v is obtained by convolving f with a Gaussian kernel).

In this paper, which is inspired from related works on segmentation [START_REF] Guyader | A combined segmentation and registration framework with a nonlinear elasticity smoother[END_REF] and on two-dimensional registration [START_REF] Lin | A landmark based nonlinear elasticity model for mouse atlas registration[END_REF], we propose an alternative approach to fluid registration. The new model is derived from a variational problem which is not in the form of a two-step algorithm and which can also produce large-magnitude deformations. For that purpose, a nonlinear elasticity smoother is introduced in three dimensions. As will be seen later, the computation of the Euler-Lagrange equations in this case is cumbersome. We circumvent this issue by introducing a second unknown, a matrix variable V , which approximates the Jacobian matrix of u. The nonlinear elastic regularizer is now applied to V , removing the nonlinearity in the derivatives of the unknown u in the Euler-Lagrange equations. The Euler-Lagrange equations are straightforwardly derived and a gradient descent method is used.

Also, allowing large deformations to occur may yield non-diffeomorphic deformation mappings. In [4], Christensen et al. proposed a regridding technique that resamples the deforming image and re-initializes the process once the value of the deformation Jacobian drops below a certain threshold. In [START_REF] Haber | Numerical methods for volume preserving image registration[END_REF], Haber and Modersitzki introduced an elastic registration model subjected to volume-preserving constraints. To ensure that the transformation g(x) = xu(x) is volume-preserving (that is, for any domain Ω, Ω dx = g(Ω) dx), they proposed the following pointwise constraint: det(I -Du(x)) -1 = 0. Pursuing in the same direction in [START_REF] Haber | Image registration with guaranteed displacement regularity[END_REF], the authors introduced a minimization problem under inequality constraints on the Jacobian.

Here we use an information-theoretic approach previously introduced in [START_REF] Yanovsky | Topology preserving log-unbiased nonlinear image registration: Theory and implementation[END_REF]. In [START_REF] Yanovsky | Topology preserving log-unbiased nonlinear image registration: Theory and implementation[END_REF], the authors considered a smooth deformation g that maps domain Ω bijectively onto itself. Consequently, g and g -1 are bijective and globally volume-preserving. Probability density functions can thus be associated with the deformation g and its inverse g -1 . The authors then proposed to quantify the magnitude of the deformation by means of the symmetric Kullback-Leibler distance between the probability density functions associated with the deformation and the identity mapping. This distance, when rewritten using skew-symmetry properties, is viewed as a cost function and is combined with the viscous fluid model for registration, which leads to an unbiased fluid registration model. Unlike the unbiased fluid registration model, the unbiased nonlinear elasticity method, introduced here, allows the functional to be written "in closed form". The new model also does not require expensive Navier-Stokes solver (or its approximation) at each step as previously mentioned.

Method

Let Ω be an open and bounded domain in R 3 . Without loss of generality, we assume that the volume of Ω is 1, i.e. |Ω| = 1. Let I 1 , I 2 : Ω → R be the two images to be registered. We seek the transformation g : Ω → Ω that maps the source image I 2 into correspondence with the target image I 1 . In this paper, we will restrict this mapping to be differentiable, one-to-one, and onto. We denote the Jacobian matrix of a deformation g to be Dg, with Jacobian denoted by |Dg(x)| = det(Dg(x)) (thus we will use the notation |V | := det(V ) for any 3 × 3 matrix V ). The displacement field u(x) from the position x in the deformed image I 2 • g(x) back to I 2 (x) is defined in terms of the deformation g(x) by the expression g(x) = xu(x) at every point x ∈ Ω. Thus, we consider the problems of finding g and u as equivalent.

In general, nonlinear image registration models may be formulated in a variational framework. The minimization problems often define the energy functional E as a linear combination of image matching term F and the regularizing term R:

inf u {E(u) = F (u) + λ 0 R(u)}.
Here, λ 0 > 0 is a weighting parameter.

Registration metrics

In this paper, the matching functional F takes the form of the L 2 norm (the sum of squared intensity differences), F = F L 2 , and the mutual information, F = F M I . L 2 -norm: The L 2 -norm matching functional is suitable when the images have been acquired through similar sensors (with additive Gaussian noise) and thus are expected to present the same intensity range and distribution. The L 2 distance between the deformed image I 2 • g(x) = I 2 (xu(x)) and target image I 1 (x) is defined as

F L 2 (u) = 1 2 Ω I 2 (x -u(x)) -I 1 (x) 2 dx. (3)
Mutual Information: Mutual information can be used to align images of different modalities, without requiring knowledge of the relationship of the two registered images [START_REF] Collignon | Automated multi-modality image registration based on information theory[END_REF][START_REF] Viola | Alignment by maximization of mutual information[END_REF]. Here, the intensity distributions estimated from I 1 (x) and I 2 (xu(x)) are denoted by p I1 and p I2 u , respectively, and an estimate of their joint intensity distribution by p I1,I2 u . We let

i 1 = I 1 (x), i 2 = I 2 (x -u(x)) denote intensity values at point x ∈ Ω. Given the displacement field u, the mutual information computed from I 1 and I 2 is provided by M I I1,I2 u = R 2 p I1,I2 u (i 1 , i 2 ) log[p I1,I2 u (i 1 , i 2 )/(p I1 (i 1 )p I2 u (i 2 ))] di 1 di 2 .
We seek to maximize the mutual information between I 2 (xu(x)) and I 1 (x), or equivalently, minimize the negative of M I I1,I2 u :

F M I (I 1 , I 2 , u) = -M I I1,I2 u .
(4)

Nonlinear Elastic Regularization

The theory of elasticity is based on the notion of strain. Strain is defined as the amount of deformation an object experiences compared to its original size and shape. In three spatial dimensions, the strain tensor, E = [ε ij ] ∈ R 3×3 , 1 ≤ i, j ≤ 3, is a symmetric tensor used to quantify the strain of an object undergoing a deformation. The nonlinear strain is defined as

ε ij (u) = 1 2 ∂ j u i + ∂ i u j + 3 k=1 ∂ i u k ∂ j u k
, with the nonlinear strain tensor matrix given by

E(u) = 1 2 Du t + Du + Du t Du . (5)
Stored energy (Saint Venant-Kirchhoff material) is defined as

W (E) = ν 2 (trace(E)) 2 + µtrace(E 2 ),
where ν and µ are Lamé elastic material constants. The regularization for nonlinear elasticity becomes

R E (u) = Ω W (E(u))dx.
The regularization term R E (u) can be minimized with respect to u. However, since the regularization term is written in terms of partial derivatives of components of u, the Euler-Lagrange equations become complicated and are computationally expensive to minimize. Instead, following earlier theoretical work [START_REF] Negron-Marrero | A numerical method for detecting singular minimizers of multidimensional problems in nonlinear elasticity[END_REF], we minimize an approximate functional by introducing the matrix variable

V ≈ Du (6)
and thus consider a new form of nonlinear elasticity regularization functional

R E (u, V ) = Ω W ( V ) dx + β 2 Ω ||V -Du|| 2 F dx, (7)
where

V = 1 2 V t + V + V t V , β
is a positive constant, and || • || F denotes the Frobenius norm. For β large enough, R E (u) is well approximated by R E (u, V ). In the limit as β → +∞, we obtain Du ≈ V in the L 2 -topology.

The idea of duplication of variables was also used in previous work [START_REF] Cachier | Iconic feature based nonrigid registration: The PASHA algorithm[END_REF] by Cachier et al. but for a different problem. In their case, the registration energy depends on two variables that are both vector fields. The first vector field C is a set of correspondences between points based on intensity similarity, while the second, denoted by T, is a smooth vector field constrained by the regularization energy and attracted by the set of correspondences C.

Unbiased Registration Constraint

In [START_REF] Yanovsky | Topology preserving log-unbiased nonlinear image registration: Theory and implementation[END_REF], the authors proposed an unbiased fluid image registration approach. Contrary to classical methods for which the term unbiased is used in the sense of symmetric registration, in [START_REF] Yanovsky | Topology preserving log-unbiased nonlinear image registration: Theory and implementation[END_REF], unbiased means that the Jacobian determinants of the deformations recovered between a pair of images follow a log-normal distribution, with zero mean after log-transformation. The authors argued that this distribution is beneficial when recovering change in regions of homogeneous intensity, and in ensuring symmetrical results when the order of two images being registered is switched. As derived in [START_REF] Yanovsky | Topology preserving log-unbiased nonlinear image registration: Theory and implementation[END_REF] using information theory, the unbiased regularization term is given as

R U B (u) = Ω |D(x -u(x))| -1 log |D(x -u(x))|dx. (8)
It is important to note that R U B generates inverse-consistent deformation maps. The inverse-consistent property of the unbiased technique was shown in a validation study of the unbiased fluid registration methods [START_REF] Yanovsky | Validating unbiased registration on longitudinal MRI scans from the Alzheimer's Disease Neuroimaging Initiative (ADNI)[END_REF]. Also, to see why minimizing equation ( 8) leads to unbiased deformation in the logarithmic space, we observe that the integrand is always non-negative, and only evaluates to zero when the deformation g is volumepreserving everywhere (|Dg| = 1 everywhere). Thus, by treating it as a cost, we recover zero-change by minimizing this cost when we compare images differing only in noise.

Given equation ( 6), we have Dg = I -Du ≈ I -V , where I is the 3 × 3 identity matrix. Therefore, as in subsection 2.2, to simplify the discretization, we introduce

R U B (V ) = Ω (|I -V | -1) log |I -V | dx. (9)
Recall that here |I -V | = det(I -V ).

Unbiased Nonlinear Elasticity Registration

The total energy functional employed in this work, is given as a linear combination of the similarity measure F (which is either F L 2 from (3) or F M I from (4)), nonlinear elastic regularization R E in (7), and unbiased regularization R U B in ( 9):

E(u, V ) = F (u) + R E (u, V ) + λR U B (V ). (10)
The explicit weighting parameter is omitted in front of R E (u, V ), since this term is weighted by Lamé constants ν and µ. We solve the Euler-Lagrange equations in u and V using the gradient descent method, parameterizing the descent direction by an artificial time t,

∂u ∂t = -∂E u (u, V ) = -∂ u F (u) -∂ u R E (u, V ), (11)
∂V ∂t = -∂E V (u, V ) = -∂ V R E (u, V ) -λ∂ V R U B (V ), (12)
which gives systems of three and nine equations, respectively. Explicit expressions for the gradients in these equations are given in Section 3.

Remark 1. The regularization on the deformation g proposed in this work can be expressed in a general form R(g) = Ω R 1 (Dg)dx + Ω R 2 (|Dg|)dx, with |Dg| := det(Dg). For the minimization, an auxiliary variable can also be introduced to simplify the numerical calculations, removing the nonlinearity in the derivatives.

3 Implementation

The Energy Gradients

Computing the first variation of functional R E (u, V ), in equation ( 7), with respect to u gives the following components of gradient ∂ u R E (u, V ):

∂ u k R E (u, V ) = β ∂ 1 v k1 + ∂ 2 v k2 + ∂ 3 v k3 -△u k , k = 1, 2, 3. The first variation of R E (u, V ) with respect to V , with V = [v ij ], gives ∂ V R E (u, V ): ∂ v11 R E (u, V ) = β(v 11 -∂ 1 u 1 ) + νc 1 (1 + v 11 ) + µ c 2 (1 + v 11 ) + c 5 v 12 + c 6 v 13 , ∂ v12 R E (u, V ) = β(v 12 -∂ 2 u 1 ) + νc 1 v 12 + µ c 3 v 12 + c 5 (1 + v 11 ) + c 7 v 13 , ∂ v13 R E (u, V ) = β(v 13 -∂ 3 u 1 ) + νc 1 v 13 + µ c 4 v 13 + c 6 (1 + v 11 ) + c 7 v 12 , ∂ v21 R E (u, V ) = β(v 21 -∂ 1 u 2 ) + νc 1 v 21 + µ c 2 v 21 + c 5 (1 + v 22 ) + c 6 v 23 , ∂ v22 R E (u, V ) = β(v 22 -∂ 2 u 2 ) + νc 1 (1 + v 22 ) + µ c 3 (1 + v 22 ) + c 5 v 21 + c 7 v 23 , ∂ v23 R E (u, V ) = β(v 23 -∂ 3 u 2 ) + νc 1 v 23 + µ c 4 v 23 + c 6 v 21 + c 7 (1 + v 22 ) , ∂ v31 R E (u, V ) = β(v 31 -∂ 1 u 3 ) + νc 1 v 31 + µ c 2 v 31 + c 5 v 32 + c 6 (1 + v 33 ) , ∂ v32 R E (u, V ) = β(v 32 -∂ 2 u 3 ) + νc 1 v 32 + µ c 3 v 32 + c 5 v 31 + c 7 (1 + v 33 ) , ∂ v33 R E (u, V ) = β(v 33 -∂ 3 u 3 ) + νc 1 (1 + v 33 ) + µ c 4 (1 + v 33 ) + c 6 v 31 + c 7 v 32 ,
where

c 1 = v 11 + v 22 + v 33 + 1 2 v 2 11 + v 2 21 + v 2 31 + v 2 12 + v 2 22 + v 2 32 + v 2 13 + v 2 23 + v 2 33 , c 2 = 2v 11 + v 2 11 + v 2 21 + v 2 31 , c 5 = v 21 + v 12 + v 11 v 12 + v 21 v 22 + v 31 v 32 , c 3 = 2v 22 + v 2 12 + v 2 22 + v 2 32 , c 6 = v 31 + v 13 + v 11 v 13 + v 21 v 23 + v 31 v 33 , c 4 = 2v 33 + v 2 13 + v 2 23 + v 2 33 , c 7 = v 32 + v 23 + v 12 v 13 + v 22 v 23 + v 32 v 33 .
We can compute the first variation of ( 9), obtaining ∂ V R U B (V ). We first simplify the notation, letting J = |I -V |. Also, denote L(J) = (J -1) log J. Hence, L ′ (J) = dL(J)/dJ = 1 + log J -1/J. Thus,

∂ v11 R U B (V ) = -(1 -v 22 )(1 -v 33 ) -v 32 v 23 L ′ (J), ∂ v12 R U B (V ) = -v 23 v 31 + v 21 (1 -v 33 ) L ′ (J), ∂ v13 R U B (V ) = -v 21 v 32 + (1 -v 22 )v 31 L ′ (J), ∂ v21 R U B (V ) = -v 32 v 13 + v 12 (1 -v 33 ) L ′ (J), ∂ v22 R U B (V ) = -(1 -v 11 )(1 -v 33 ) -v 13 v 31 L ′ (J), ∂ v23 R U B (V ) = -v 12 v 31 + v 32 (1 -v 11 ) L ′ (J), ∂ v31 R U B (V ) = -v 12 v 23 + v 13 (1 -v 22 ) L ′ (J), ∂ v32 R U B (V ) = -v 21 v 13 + v 23 (1 -v 11 ) L ′ (J), ∂ v33 R U B (V ) = -(1 -v 11 )(1 -v 22 ) -v 12 v 21 L ′ (J).

Algorithm

We are now ready to give the algorithm for the unbiased registration via nonlinear elastic regularization.

Algorithm 1 Unbiased Registration via Nonlinear Elastic Regularization 1: Initialize t = 0, u(x, 0) = 0, and V (x, 0) = 0. 2: Calculate V (x, t) using equation [START_REF] Haber | Numerical methods for volume preserving image registration[END_REF].

Steps 3-5 describe the procedure for solving equation [START_REF] Lin | A landmark based nonlinear elasticity model for mouse atlas registration[END_REF]. 3: Calculate the perturbation of the displacement field R(x) = -∂Eu(u, V ). 4: Time step △t is calculated adaptively so that △t • max(||R||2) = δu, where δu is the maximal displacement allowed in one iteration. Results in this work are obtained with δu = 0.1. 5: Advance equation [START_REF] Lin | A landmark based nonlinear elasticity model for mouse atlas registration[END_REF], i.e. ∂u(x, t)/∂t = R(x), in time, with time step from step 4, solving for u(x, t). 6: If the cost functional in [START_REF] Guyader | A combined segmentation and registration framework with a nonlinear elasticity smoother[END_REF] decreases by sufficiently small amount compared to the previous iteration, then stop. 7: Let t := t + △t and go to step 2.

Results and Discussion

We tested the proposed unbiased nonlinear elastic registration model and compared the results to those obtained with the unbiased fluid registration method [START_REF] Yanovsky | Topology preserving log-unbiased nonlinear image registration: Theory and implementation[END_REF], where the unbiased regularization constraint (8) was coupled with the L 2 matching functional (3) and fluid regularization (1), [START_REF] Lord | Simultaneous registration and parcellation of bilateral hippocampal surface pairs for local asymmetry quantification[END_REF]. Here, both methods were coupled with the L 2 and mutual information (MI) based similarity measures. In our experiments, we used a pair of serial MRI images (220 × 220 × 220) from the Alzheimer's Disease Neuroimaging Initiative (ADNI). Since the images were acquired one year apart, from a subject with Alzheimer's disease, real anatomical changes are present, which allows methods to be compared in the presence of true biological changes.

In the tests performed using unbiased nonlinear elasticity coupled with L 2 matching, values of β = 20000 in equation ( 7) and λ = 2000 in equation [START_REF] Guyader | A combined segmentation and registration framework with a nonlinear elasticity smoother[END_REF] were chosen. For MI matching, β = 80 and λ = 8 were used. The values of the Lamé coefficients were chosen to be equal, µ = ν, in all experiments. Bigger values of µ and ν allow for more smoothing. For unbiased fluid registration model, described in [START_REF] Yanovsky | Topology preserving log-unbiased nonlinear image registration: Theory and implementation[END_REF], λ = 500 was chosen for L 2 matching, and λ = 5 for MI matching.

Figure 2 shows the images being registered along with the resulting Jacobian maps. Results generated using the fluid and nonlinear elasticity based unbiased models are similar, both suggesting a mild volume reduction in gray and white matter and ventricular enlargement that is observed in Alzheimer's disease patients. The advantages of the unbiased nonlinear elasticity model is its more locally plausible reproduction of atrophic changes in the brain and its robustness to original misalignment of brain volumes, which is especially noticeable on the brain surface. The unbiased nonlinear elasticity model coupled with L 2 matching generated very similar results to those obtained with the MI similarity measure, partly because difference images typically contain only noise after registration. Unbiased fluid registration method, however, is more effective in modeling the regional neuroanatomical changes, showing more clearly which parts of the volume have undergone largest tissue changes, such as ventricular enlargement as shown in Figure 2.

Figure 3 shows deformed grids generated with unbiased fluid and unbiased nonlinear elastic registration models. Figure 4 shows the energy decrease per iteration for both models.

In Figure 5, we examined the inverse consistency of the mappings [20] generated using unbiased nonlinear elastic registration. Here, the deformation was computed in both directions (time 2 to time 1, and time 1 to time 2) using mutual information matching. The forward and backward Jacobian maps were concatenated (in an ideal situation, this operation should yield the identity), with the products of Jacobians having values close to 1.

The unbiased nonlinear elasticity model does not require expensive Navier-Stokes solver (or its approximation), which is employed at each iteration for fluid flow models. Hence, in our experiments, unbiased nonlinear elasticity iteration (based on explicit scheme) took 15-20% less time than the unbiased fluid step. Convergence was obtained after roughly the same number of iterations for both methods, resulting in better performance for the unbiased nonlinear elasticity model.

To conclude, we have provided an alternative unified minimization approach to the unbiased fluid registration model and have compared both models. The proposed method proves to be easier to implement and is less computationally intensive. Also, a key benefit of the variational framework and of the numerical scheme of the unbiased

Fig. 1 .

 1 Fig. 1. Serial MRI images from the ADNI follow-up dataset (images acquired one year apart) are shown. Volumes I1 (row 1) and I2 (row 2) are depicted as a brain volume (column 1) and from sagittal (column 2), axial (column 3), and coronal (column 4) views. Nonrigid registration aligns volume I2 into correspondence with volume I1.

Unbiased

  Registration via Viscous Fluid Flow coupled with L 2 matching Unbiased Registration via Nonlinear Elasticity coupled with L 2 matching Unbiased Registration via Viscous Fluid Flow coupled with Mutual Information Unbiased Registration via Nonlinear Elasticity coupled with Mutual Information

Fig. 2 .

 2 Fig. 2. Nonrigid registration was performed on the Serial MRI images from the ADNI Followup dataset using unbiased fluid registration and unbiased nonlinear elasticity registration, both coupled with L 2 and MI matching. Jacobian maps are superimposed on the target volume.

Fig. 5 .

 5 Fig. 5. This figure examines the inverse consistency of the unbiased nonlinear elastic registration. Here, the model is coupled with mutual information matching. Jacobian maps of deformations from time 2 to time 1 (column 1) and time 1 to time 2 (column 2) are superimposed on the target volumes. The products of Jacobian maps, shown in column 3, have values close to 1, suggesting inverse consistency.

⋆ This work was funded by the National Institutes of Health through the NIH Roadmap for Medical Research, Grant U54 RR021813 entitled Center for Computational Biology (CCB).

Unbiased models with L 2 matching Unbiased models with MI matching Fluid Nonl.Elasticity Fluid Nonl.Elasticity nonlinear elastic registration model is its robustness to numerical constraints such as CFL conditions. The method allows to remove the nonlinearity in the derivatives of the unknown u in the Euler-Lagrange equations. Future studies will examine the registration accuracy of the different models where ground truth is known, and will compare each model's power for detecting inter-group differences or statistical effects on rates of atrophy.


Christophe Collewet

POLAR SNAKES: A FAST AND ROBUST PARAMETRIC ACTIVE CONTOUR MODEL

Keywords: Active contour model, polar description, segmentation, contour tracking

 to show the benefit of using a polar description to model the contour, in particular to cope with the well-known initialization problem. Indeed, we show that this way to proceed leads to diagonal and constant matrices in the equations of the snake evolution yielding therefore to a faster algorithm. Experimental results on image segmentation and contour tracking validate the efficiency of this new formulation.

INTRODUCTION

Since their first introduction in [START_REF] Kass | Snakes: Active contour models[END_REF], active contours or snakes have proven to be very powerful in various contexts, including image segmentation [START_REF] Tauber | Robust b-spline snakes for ultrasound image segmenta-tion[END_REF], edge detection [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Brigger | B-spline snakes: a flexible tool for parametric contour detection[END_REF] and target tracking [START_REF] Bascle | Tracking complex primitives in an image sequence[END_REF][START_REF] Precioso | Robust real-time segmentation of images and videos using a smoothing-spline snake-based algorithm[END_REF]. An active contour is a deformable curve defined in the image domain that evolves according to internal forces coming from the curve itself and external forces generated from the image (most often gradient-based or regionbased information). Numerous approaches are reported in the literature. Nevertheless, they can be classified into two main categories: geometric snakes [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Malladi | Shape modeling with front propagation: A level set approach[END_REF] where the curves is described from a level set representation, or parametric snakes where the curve is either described as a discrete collection of points [START_REF] Kass | Snakes: Active contour models[END_REF][START_REF] Xu | Snakes, shapes, and gradient vector flow[END_REF] or by an appropriate parametric function such as superquadrics [START_REF] Terzopoulos | Dynamic 3d models with local and global deformations : deformable superquadrics[END_REF] or B-splines [START_REF] Tauber | Robust b-spline snakes for ultrasound image segmenta-tion[END_REF][START_REF] Bascle | Tracking complex primitives in an image sequence[END_REF] for example. In this paper we are interested in the latter case, that is in parametric active contours. Indeed, we will consider here robotics issues such as 3D reconstruction by a moving camera or robot control by visual servoing [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF]. In those cases, a high acquisition rate is required and thus the image processing duration has to be as near as possible to the video rate. Currently, it is well known that geometric approaches are too much time consuming and consequently, they cannot be used in this context. In addition, compared to traditional methods based on points, parametric snakes based on parametric functions are less time consuming. Such works are reported for example in [START_REF] Terzopoulos | Dynamic 3d models with local and global deformations : deformable superquadrics[END_REF][START_REF] Brigger | B-spline snakes: a flexible tool for parametric contour detection[END_REF][START_REF] Bascle | Tracking complex primitives in an image sequence[END_REF][START_REF] Precioso | Robust real-time segmentation of images and videos using a smoothing-spline snake-based algorithm[END_REF].

However, the major drawback of parametric active contours is that they are very sensitive to the initial position of the contour. Indeed, it has to be chosen close to the true contour to ensure its convergence to the correct position. Numerous approaches have been proposed to cope with this problem. A pressure force acts as a force that can inflate or deflate the contour according to the normal at each point belonging to it [START_REF] Cohen | On active contour models and balloons[END_REF]. Another approach is to introduce an external force based on a distance map between a point belonging to the active contour and the closest boundary point [START_REF] Cohen | Finite-element methods for active contour models and balloons for 2-d and 3-d images[END_REF]. Multiresolution methods have also been proposed to cope with this problem [START_REF] Brigger | B-spline snakes: a flexible tool for parametric contour detection[END_REF][START_REF] Leroy | Multi-resolution algorithms for active contour models[END_REF]. A new external force, called GVF, computed as a diffusion of the image gradient map has been presented in [START_REF] Xu | Snakes, shapes, and gradient vector flow[END_REF] to ensure a convergence even from an initial contour far from the true one. However, computing such forces is very time consuming.

The main contribution of this paper is to show the benefit of using a polar description as the parametric curve to describe the contour. In addition, we introduce an external force derived from an energy term based on the area inside the contour. It is very closed in mind to the pressure forces [START_REF] Cohen | On active contour models and balloons[END_REF]. However, since the polar description is based on orthogonal functions, the evolution equations (including the energy term) of this kind of active contour are very simple and, consequently, make the computations faster.

This paper is structured as follow. First, we revisit in Section 2 the evolution equations of a physics-based deformable model as proposed in [START_REF] Terzopoulos | Dynamic 3d models with local and global deformations : deformable superquadrics[END_REF] for a generic family functions. Next, this result is used in Section 3 to show that a polar description highly simplifies these equations. Section 4 points out that it becomes very simple to introduce an area-based energy term. The experimental results are presented in Section 5.

EVOLUTION EQUATIONS OF THE CONTOUR

First, let us consider the generic family curve C q depending on a n-dimensional parameter vector q and defined in the image plane by

C q (u) : x(u) = x c + l=n l=0 q l Φ l (u) (1)
x c being a point inside the contour, u ∈ [0, 2π] and Φ l (l = 0 • • • n) some 2d vector functions. Note that C q writes as a linear form with respect to q.

To derive the evolution equation of C q , we search the evolution of q as proposed in [START_REF] Terzopoulos | Dynamic 3d models with local and global deformations : deformable superquadrics[END_REF]. More precisely, this curve is supposed to be embedded in a viscous medium and supposed to be undergone to internal forces f i and to external forces f e . This problem can be set according to the formalism of Lagrangian mechanics when considering q as the generalized coordinates of the system. The Euler-Lagrange equations give for each component q i of q:

d dt ∂L ∂ qi - ∂L ∂q i = Q i (2)
where Q i describes the generalized force associated to q i :

Q i = Cq f ⊤ (u) ∂C q (u) ∂q i du. (3)
f are the dissipative forces, L the Lagrangian, it is given by L = T -U with T the kinetic energy and U the potential energy. T writes as

T = Cq 1 2 µ C t 2 du (4)
with C t = ∂C q /∂t and µ the mass of each point supposed to constitute C q . Most often, U is described from the Tikhonov regularization terms:

U = Cq 1 2 k 1 C u 2 du + Cq 1 2 k 2 C uu 2 du (5)
where C u = ∂Cq/∂u, C uu = ∂Cq/∂u. The scalar k 1 will tend to limit the extension of C q while k 2 will tend to limit its curvature. According to (1), we have

     C t = l ql Φ l C u = l q l Φ ′ l (with Φ ′ l = ∂Φ l /∂u) C uu = l q l Φ ′′ l (with Φ ′′ l = ∂ 2 Φ l /∂ 2 u). (6)
Therefore, it becomes straightforward to evaluate the contribution of U in (2): we have ∂U/∂ qi = 0 and

∂U ∂q i = Cq l k 1 q l Φ ′ l ⊤ Φ ′ i + l k 2 q l Φ ′′ l ⊤ Φ ′′ i du. (7)
Similarly, we have ∂T /∂q i = 0 and

1 2 ∂ C t 2 ∂ qi = l ql Φ l ⊤ Φ i (8)
leading to

d dt ∂T ∂ qi = Cq µ l ql Φ l ⊤ Φ i du. (9)
Thus, from ( 7) and ( 9), the left part of (2) can be evaluated. Concerning, the right part, we first consider a force due to a viscous friction f vis = -γC t , where γ is the viscosity of the medium, and, second, a force f im due to the image. In addition, we suppose that f im derives from a potential energy term E im : f im = -∇E im where ∇E im denotes the spatial gradient of E im in the image. In this paper, we will use E im = -∇I . Concerning the contribution of vis , we have:

Cq f ⊤ vis ∂C q (u) ∂q i du = - Cq γ l ql Φ l ⊤ Φ i du. (10)
Finally, from ( 7), ( 9) and ( 10), ( 2) can be written under a matrix form:

Mq + C q + Kq = Q(q) = - Cq ∇E im ∂C q (u) ∂q du (11)
where

M = [M mn ], C = [C mn ] and K = [K mn ] with                M mn = µ Cq Φ ⊤ m Φ n du C mn = γ Cq Φ ⊤ m Φ n du K mn = k 1 Cq Φ ′ ⊤ m Φ ′ n du + k 2 Cq Φ ′′ ⊤ m Φ ′′ n du. (12)
To solve the differential system [START_REF] Cohen | On active contour models and balloons[END_REF], it is more convenient to set the mass µ to zero as proposed in [START_REF] Terzopoulos | Dynamic 3d models with local and global deformations : deformable superquadrics[END_REF]. Moreover, we use an explicit Euler's method, leading to

q k+1 = q k + ∆tC -1 Q(q) -Kq ( 13
)
where ∆t is the sampling rate. In practice, its value must be small since E im is only valid in a small domain around q k . Lets us point out that the computation of matrices C and K can be time consuming. Indeed, they are obtained from an integration along C q . Besides, the terms to integrate may be themselves complex to compute. That is why we propose in the next section a way to decrease the computation time of [START_REF] Leroy | Multi-resolution algorithms for active contour models[END_REF]. It is based on a judicious choice for Φ l .

USING A POLAR DESCRIPTION

A way to simplify ( 12) is to use orthogonal functions such that their first and second derivative are also orthogonal in the sense of the inner product defined by the integrals in [START_REF] Cohen | Finite-element methods for active contour models and balloons for 2-d and 3-d images[END_REF]. Note that it cannot be achieved using traditional orthogonal polynomials. Conversely, it can be done using a polar description of the contour C q :

C q (u) : x(u) = x c + ρ q (u) cos u sin u (14)
and by performing a Fourier expansion of the radius ρ q (u):

ρ q (u) = q 0 + k=h k=1 q k cos ku + q k+h sin ku ( 15
)
where h is the number of harmonics. Thus, we have q = (q 0 , • • • , q 2h ) and consequently:

Φ l = cos u cos lu sin u cos lu for 0 ≤ l ≤ h (16)
and

Φ l = cos u sin(l -h)u sin u sin(l -h)u for h < l ≤ 2h. ( 17
)
In that case, it is easy to show that the matrices involved in ( 12) becomes diagonal and constant. Indeed, we simply have

C = γ diag (2π, π, • • • , π) (18) K 1 = k 1 diag (2π, π(1 + 1 2 ), • • • , π(1 + h 2 ), π(1 + 1 2 ), • • • , π(1 + h 2 ) ( 19
)
K 2 = k 2 diag (2π, π(1 + 6×1 2 + 1 4 ), • • • , π(1 + 6h 2 + h 4 ), π(1 + 6×1 2 + 1 4 ), • • • , π(1 + 6h 2 + h 4 ) (20) with K 1 and K 2 defined such that K = K 1 + K 2 .

USING AN AREA-BASED ENERGY TERM

We are interested in this section to add an external force to the active contour model to help the curve to reach the true contour even if its initialization is far from it. However, adding such a force is time consuming since its contribution to (3) requires the computation of an integral of a complex term. It is the case for example with the pressure forces proposed in [START_REF] Cohen | On active contour models and balloons[END_REF] where the normal has to be computed at each point. However, as pointed out in [START_REF] Cohen | Finite-element methods for active contour models and balloons for 2-d and 3-d images[END_REF], this force can also be interpreted as the gradient of an area energy term E area defined by E area = k 3 S q where S q denotes the area of the inside region delimited by the curve and k 3 a scalar. If k 3 is a positive value, the area will decrease until the contour is reached. Thus, an initial estimation outside the contour is required. Conversely, if we want to consider an initialization inside the contour, a negative value has to be chosen.

Thanks to the polar description, this surface is very simple to compute from the well-known Green-Riemann's formula. Indeed, we have

S q = π q 2 0 + 1 2 k q 2 k (21)
that also leads to a very simple way to compute ∂E area /∂q i involved in [START_REF] Kass | Snakes: Active contour models[END_REF]. Indeed, we only have to introduce the following diagonal and constant matrix:

K 3 = k 3 diag (2π, π, • • • , π) (22)
that we have to add to the matrices K 1 and K 2 .

EXPERIMENTAL RESULTS

The first experiment is devoted to show the robustness of the convergence wrt a bad initialization of the active contour. The final goal of this experiment is to allow disable persons to grasp an object with the help of a vision system and a robot arm mounted on their wheel chair [START_REF] Dune | Active rough shape estimation of unknown objects[END_REF]. Therefore, only one click has been used. Fig. 1 shows the behavior of the active contour when objects of the everyday life are considered. The yellow quadrilateral is the initial position of the contour, it has been built from the click (the green cross) and the nearest image border. The curves in blue shows the behavior of the active contour during its convergence; the curve in green depicts the active contour when q no more evolves. As can be seen, using a surface-based energy term is very efficient. For this experiment, the following values have been used for both objects: h = 20, k1 = 0, k2 = 5.10 -4 , k3 = 0.1, γ = 500 (where γ = γ∆t and k i = ki ∆t).

The second experiment (see Fig. 2) is also devoted to the effectiveness of the initialization process. An inside initialization has been considered. Here we have chosen : h = 30, k3 = -0.2 (the other parameters remain unchanged). The left image shows a segmentation of a red muscle of a raw ham. The right one, the segmentation of MRI image of a pork lung. Here again, despite bad initializations, the contour converges to the true one.

The goal of the last experiment is to show the capability of our approach to achieve real-time applications (see Fig. 3). Indeed, we are interested here to control by visual servoing [START_REF] Chaumette | Visual servo control, Part I: Basic approaches[END_REF] the motion of a 2D ultrasound probe held by a robot. The control law is designed so that the current contour (the green one) reaches a learned contour (the red one) leading to a correct position of the probe. More precisely, the control law is based on some 2D images moments (computed from the vector q, see [START_REF] Mebarki | Automatic guidance of an ultrasound probe by visual servoing based on b-mode image moments[END_REF] for more details). Therefore, to achieve such a task the current contour has to be tracked to compute the control law. To do that, q k-1 is used to initialize the contour at time k. Using a polar description allowed to reach the video rate. For this experiment, we have used: h = 6, k1 = 0, k2 = 1.10 -2 , k3 = 0 and γ = 200.

CONCLUSION

We have presented in this paper a physics-based deformable model in the case where the contour is described by a polar description. We have shown that this representation leads to a simple evolution equation of the contour suitable for real-time applications as contour tracking for visual servoing issues. In addition, we have also shown that, with this representation, it is very easy to introduce an additional area-based energy term that efficiently deals with a contour initialization far from the contour to reach. Experiments on various objects have validated our approach as well to perform image segmentation as for contour tracking.

Fig. 1 .

 1 Fig. 1. Robustness in the case of an outside initialization for objects of the everyday life.

Fig. 2 .

 2 Fig. 2. Robustness in the case of an inside initialization for various objects.

Fig. 3 .

 3 Fig. 3. Tracking contour to achieve a visual servoing task.