Tao Li

Alexandre Krupa

Christophe Collewet

A robust parametric active contour based on Fourier descriptors

Keywords: Active contour, Fourier descriptors, contour tracking, ultrasound image

This paper presents an optimal parametric active contour method based on Fourier descriptors to detect and track a contour from a sequence of images in real-time. The advantage of this approach is that not only convex but also concave contours can be considered. We also propose a method to initialize the active contour using only three user's clicks in the first image of the sequence and to automatically re-initialize it when possible topological changes are detected during the tracking. Moreover, the algorithm is implemented on GPU to ensure the real-time tracking performance.

INTRODUCTION

An active contour, usually called snake, is an energy minimizing curve defined in image domain that evolves according to internal forces of the curve itself and external forces generated by the image. The snake was firstly introduced in [START_REF] Kass | Snakes: Active contour models[END_REF] and has been widely applied in computer vision. In general, the snakes can be classified in two main categories: parametric snakes [START_REF] Terzopoulos | Dynamic 3d models with local and global deformations: deformable superquadrics[END_REF][START_REF] Precioso | Robust real-time segmentation of images and videos using a smooth-spline snake-based algorithm[END_REF] and geometric snakes [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Xu | A summary of geometric level-set analogues for a general class of parametric active contour and surface models[END_REF]. In our work we are interested in tracking the contour of an object of interest from a sequence of images in real-time. For that reason, we focus in this paper on parametric snakes instead of geometric ones which are not suitable for real-time application due to the high time-consuming process. In our previous work [START_REF] Collewet | Polar snakes: A fast and robust parametric active contour model[END_REF], such a parametric snake has been proposed. Indeed, by using a polar description we have shown that the snake evolution equations could be simplified to a diagonal system of ordinary differential equations. Unfortunately, in polar frame, only one point can be described by polar description for each angle. Thus the polar-snake can not describe concave contours, as for example the ones depicted in Fig. 1. This can be done by B-spline parametric snakes as proposed in [START_REF] Precioso | Robust real-time segmentation of images and videos using a smooth-spline snake-based algorithm[END_REF] for example, but its evolution equations system is no more a diagonal system. In this paper, we propose to model the active contour by using Fourier descriptors (FD) [START_REF] Impedovo | A fourier descriptor set for recognizing nonstylized numerals[END_REF], which lead to a

The authors acknowledge the support of the ANR project PROSIT of the French National Research Agency and the S2E2 cluster. email: firstname.name@inria.fr diagonal system, and have the advantage to describe concave contours in a cartesian frame.

The main contribution of this paper is the determination of the evolution equations of the FD-based snake in order to detect and track an object of interest from an images sequence. Generally, a parametric active contour has to be initialized as nearest as possible to the real contour of the object of interest in order to converge to the correct position. Therefore, we propose as a second contribution of this paper, an approach to automatically initialize the snake. The method consists in initializing the snake to fit a set of contour points that are coarsely extracted using a fast Freeman chain coding algorithm (FCC) [START_REF] Freeman | A corner-finding algorithm for chain-coded curves[END_REF] applied on a binary image. After this automatic initialization, the FD-snake can finely evolve to reach the exact contour from the grey-level image by minimizing its energy criteria. Moreover, the coarsely extracted contour is also used to detect possible brutal topological changes of the object contour shape that require a FD-snake re-initialization.

To take into account the constraint of real-time contour detection, we also propose in this paper a parallel implementation of the FD-snake on GPU using the CUDA parallel computing architecture [START_REF] Yang | Parallel image processing based on cuda[END_REF] to reduce the processing time.

This paper is structured as follows. First, the description of the parametric snake from FD and its evolution equations are developed in Section 2. Then the method we propose to automatically initialize the snake and to detect the topological changes of the object shape is presented in Section 3. Section 4 provides the principle of the snake implementation on GPU. Finally, experimental results of the contour extraction are presented and discussed in Section 5.

EVOLUTION EQUATIONS OF THE FD-SNAKE

Generic equation of snake

Let us define a generic curve model C q of the parametric snake in the image plane by:

C q (u) : x(u) = x c + n ∑ l=0 q l Φ l (u) (1)
C q is a linear form depending on a n-dimensional parameters vector q and x c is a point inside the area delimited by the contour; u ∈ [0, 2π] and Φ l (l = 0, . . . , n) are 2D vector functions that can be substituted by different descriptor terms, such as polar description, B-spline and Fourier descriptors. As we did in [START_REF] Collewet | Polar snakes: A fast and robust parametric active contour model[END_REF], the evolution of the curve is modeled according to the Lagrangian mechanics formalism by considering q as the generalized coordinates of the system. The Euler-Lagrange equation of each component q i is defined by:

d dt ( ∂T ∂ qi ) - ∂U ∂q i = Q i ( 2
)
where T is the kinetic energy and U is the potential energy of the active contour. Q i describes the generalized forces associated to q i :

Q i = ∫ Cq f ⊤ (u) ∂C q (u) ∂q i du (3)
with f the total of the dissipative forces defined by f = f vis + f im . The forces f vis and f im are respectively due to a viscous friction and the image itself. Here, we consider

f im = -∇E im
where E im is the potential energy term in the image I: E im = -∥∇I∥, with ∥∇I∥ the norm of the image gradient ∇I.

Finally, (2) can be written under a vectorial form (see [START_REF] Collewet | Polar snakes: A fast and robust parametric active contour model[END_REF] for more details):

Mq+C q+Kq = Q im (q) = - ∫ Cq ∇E ⊤ im (u) ∂C q (u) ∂q du (4)
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 (5)
with the mass density µ, the viscosity of the medium γ. The scalar k 1 tends to limit the extension of C q while k 2 tends to limit its curvature. Replacing Φ l (u) of ( 1) by orthogonal functions whose first and second derivatives are also orthogonal can simplify (5) to diagonal matrices. As proposed in [START_REF] Terzopoulos | Dynamic 3d models with local and global deformations: deformable superquadrics[END_REF], we set the mass density µ to zero to cancel the inertial behavior of the active contour in order to improve its dynamic stability. The differential system (4) can then be solved as a discrete optimization problem using explicit Euler's method leading to the following recursive expression of q:

q k+1 = q k + ∆tC -1 (Q im (q k ) -Kq k ) (6)
where q k+1 is the new estimate of the generalized coordinates, q k is the previous estimate and ∆t is the sampling rate.

Parametrization using the Fourier descriptors

Replacing the term Φ l of (1) by Fourier descriptors [START_REF] Impedovo | A fourier descriptor set for recognizing nonstylized numerals[END_REF], the curve model C q becomes:

C q (u) : x(u) = ( a 0 c 0 ) + h ∑ k=1 ( a k b k c k d k ) ( cos ku sin ku ) ( 7
)
where h is the number of harmonics. By developing (7), we have q = (a 0 , a

1 • • • a h , b 1 • • • b h , c 0 , c 1 • • • c h , d 1 • • • d h ) and consequently:                                    0 ≤ l < h : Φ l (u)= ( cos lu 0 ) h+1 ≤ l < 2h : Φ l (u)= ( sin(l -h)u 0 ) 2h+1 ≤ l < 3h+1 : Φ l (u)= ( 0 cos(l -2h -1)u ) 3h+2 ≤ l < 4h+2 : Φ l (u)= ( 0 sin(l -3h -1)u ) (8)
Furthermore, by substituting ( 8) in ( 5), the matrices C and K = K 1 + K 2 involved in (4) become diagonal and constant:

                       C = γ diag(2π, π . . . π, π . . . π, 2π, π . . . π, π . . . π) K 1 = k 1 diag(0, π . . . πh 2 , π . . . πh 2 , 0, π . . . πh 2 , π . . . πh 2 ) K 2 = k 2 diag(0, π . . . πh 4 , π . . . πh 4 , 0, π . . . πh 4 , π . . . πh 4 ) (9)
For that reason, the evolution equations system (4) becomes optimal and the computation time of (6) highly decreases. In addition, to ensure a convergence of the snake toward the true contour, even in the case of an initialization far away from it, we propose to add a surface-based energy term as in [START_REF] Collewet | Polar snakes: A fast and robust parametric active contour model[END_REF]. It can be simply done since the surface S q can be analytically computed. Indeed we have:

S q = π( h ∑ k=1 k(a k d k -b k c k )) (10)

AUTOMATIC INITIALIZATION AND TOPOLOGICAL CHANGE DETECTION

As mentioned in the introduction, the initialization of active contour and possible topological changes of the curve during contour tracking are two important issues. To deal with them, we propose to automatically initialize the snake from a set of contour points that are roughly extracted. This extraction is done by binarizing the current image using an optimal adaptive threshold that we compute thanks to the Otsu's method [START_REF] Otsu | A threshold selection method from gray-level histograms[END_REF]. Only the grey-level pixels included in an image bounding box (Fig. 1(a)) are taken into account for the threshold adaptation.

The coarse contour L p (t) and the surface A(t) of the object of interest in the binary image are then obtained by applying the fast Freeman chain coding algorithm [START_REF] Freeman | A corner-finding algorithm for chain-coded curves[END_REF]. At time t = 0, the snake is then initialized from L p (t = 0) in the first image. Furthermore, the snake will be automatically re-initialized if the object of interest separates in two parts (Fig. 2) or a large difference between A(t) and A(t -∆t) is detected from the following images of the sequence (Fig. 3).

IMPLANTATION ON GPU

In order to perform real time contour tracking from a sequence of ultrasound images, for example in the context of an image-guided robotic medical application [START_REF] Mebarki | 2d ultrasound probe complete guidance by visual servoing using image moments[END_REF], we propose in this section a GPU implementation solution of the snake.

In ( 4), the diagonal matrices C and K are constant. q k is the optimized result obtained from previous iteration. So we are interested to perform parallel computation of the term

Q im (q) = - ∫ Cq ∇E ⊤ im (u) ∂Cq(u)
∂q du. In this paper we consider that ∇E im , = -(∥∇I∥ x , ∥∇I∥ y ) leading to:

Q im (q)= - N ∑ n=1 ∂x(n) ∂q i ∥∇I∥ x (n) - N ∑ n=1 ∂y(n) ∂q i ∥∇I∥ y (n) (11)
with N the points number of the snake. The partial derivation terms ∂x(n) ∂qi and ∂y(n) ∂qi are constant according to FD. To gain time-consuming performance, the following algorithm was implemented on GPU:

Step 1: Copy the constant terms C, K, ∂x(n) ∂qi and ∂y(n) ∂qi from host memory (CPU) to device memory (GPU).

Step 2: If k = 0, copy q k=0 from CPU to GPU.

Step 3: Copy the current image from CPU to GPU.

Step 4: Compute ∥ ∇I ∥ in parallel, the image gradient norm of each pixel is computed by a different GPU thread, therefore all the pixels can be processed on GPU in parallel.

Step 5: Compute the spatial gradient ∇E im in parallel.

Step 6: Compute Q im in parallel.

Step 7: Compute q k+1 in parallel.

Step 8: If |q k+1 -q k | > ϵ, then k = k+1 and go to Step 5.

Step 9: If |q k+1 -q k | ≤ ϵ, copy q k+1 from GPU to CPU.

Step 10: Go to step 3 for the next image of the sequence.

EXPERIMENTAL RESULTS

All the following experiments were performed with a computer equipped with a dual core CPU Intel Xeon 5160, 4Gb memory and a graphic card NVIDIA GeForce GTX 285 providing a 240 cores GPU. Fig. 1(a) shows that only 3 points are required to be defined by the user in the first image of the sequence to initialize the snake. The point P 1 indicates the object of interest, the points P 2 and P 3 define the position of the bounding box (the white rectangle in Fig. 1.(a)) where the adaptive threshold will be computed. The red line around the joystick cover shows the snake result after convergence in the image. Fig. 1(b) and Fig. 1(c) are results of the polar-snake and the FD-snake with an echographic concave object section. In Fig. 1(b) we can see that the polar-snake diverges while the FD-snake perfectly adapts to the true contour of the concave object both for the example of Fig. 1 Fig. 2 shows the results obtained from a sequence of US images where a fission of the object of interest was observed. Two coarsely extracted contours L G1 and L G2 are extracted from two different germs G 1 and G 2 used by the FCC. If the obtained contours L G1 and L G2 are different, we consider that the tracked object separates in two parts (Fig. 2.(c)) and in this case we automatically initialize a second snake (the red one in Fig. 2.(c)).

Fig. 3 shows the results obtained from a sequence of US images where a topological change of the gallbladder was observed. In Fig. 3(a), the gallbladder is seen as an ellipse at the beginning of the US images sequence. During the object tracking, we apply a rotation motion to the US probe in order to see the gallbladder as a tube (Fig. 3(d)). Thanks to the method proposed in Section 3, the FD-snake can track this topological change by automatically re-initialize itself when a large difference of the coarse contour surface A(t) extracted by the FCC is observed between 2 successive images.

Moreover we compare the time-consuming of the CPU and GPU algorithm of the FD-snake during the tracking of the gallbladder (Fig. 3). Table1 reports the processing time of the CPU and GPU algorithm with different harmonics number h and points number N of the active contour. As shown in Table1, the GPU implementation consumes less time than the CPU for all experiments. Moreover, when h and N increase, the consuming time of the GPU implementation slightly increases but still ensures the real-time constraint. Note that this saving time can be used to increase h and / or N .

CONCLUSION

We have presented in this paper an optimal FD-based parametric active contour which has the advantage to track both convex and concave contours. A method is also proposed to semi-automatically initialize the snake in the first image and automatically re-initialize it when topological changes are detected in the current image. Finally, we provided a GPU implementation solution to perform the contour tracking in real-time. In a future work we plan to develop a visual servoing method based on the FD-snake to control a medical robot equipped with an US probe (see [START_REF] Mebarki | 2d ultrasound probe complete guidance by visual servoing using image moments[END_REF]).

Fig. 1 .

 1 Fig. 1. (a): Initialization of the FD-snake from three clicks (with h = 19). (b) Polar-snake divergence on a concave US section (with h = 6). (c) FD-snake convergence on the same concave US section (with h = 6).

Fig. 2 .Fig. 3 .

 23 Fig. 2. Topological change (fission) detection of the object of interest and an automatically initialization of the second snake along a sequence of US images (with h = 6).

Table 1 .

 1 The consuming time of CPU and GPU algorithm in function of harmonic number h and points number N .

			Time
	h N	CPU algorithm GPU algorithm
	3	100 400	75ms 82ms	29ms 35ms
	6	100 400	80ms 95ms	31ms 40ms


Keywords: Tense-Aspect-Mood, V, nonimplosive voiced obstruent, Vb: Verb, VN : Verbal Noun

South-Bauchi West TAM System 1 Bernard CARON LLACAN (UMR 8135 : CNRS, Inalco, Paris 7) caron@vjf.cnrs.fr

The 27 or so lects grouped into the SBW sub-group of Chadic languages [START_REF] Shimizu | The Southern Bauchi Group of Chadic Languages. A survey report[END_REF] are spoken in the South of the Bauchi state (Nigeria) stretching from Bauchi to the south of the state over less than 100 kms. Among them, only Zaar is relatively less unknown [START_REF] Caron | Review of 'Barawa Lexicon: a wordlist of eight South Bauchi (West Chadic) languages: Boghom, Buli, Dott, Geji, Sayanci and Zul[END_REF][START_REF] Schneeberg | Sayanci verb tonology[END_REF][START_REF] Schneeberg | Sayanci Phonology[END_REF] and is often used to characterise the whole SBW group.

These languages or dialects are obviously closely related when their basic vocabulary is compared. However, various elements tend to set Zaar (Saya) and Gùùs (Sigidi) apart :

 From a socio-linguistic and historical point of view : the name ßarawa used by Shimizu to name the group is not acknowledged by the Zaar (Sayawa) and Gùùs (Sigidi) people. If all the ßarawa agree to be closely related and all originate from the same area, viz. the Buli and Zaranda hills near Bauchi, the Zaars say they came from Bornu through Zaria ; they say they settled in the hills South of Tafawa Balewa after some cohabitation with the Angas, Pyem and Jarawa.

 Zaar and Gùùs have three tone levels, the others only two. 2

 Zaar and Gùùs have a very elaborate and profusive TAM system when compared to the other languages.

All these elements single out the SBW group of languages for what Petr Zima has called a laboratory for the study of language evolution. This gives us an opportunity to first describe the extraordinary variation among languages that are so small in terms of number of speakers, and so close in terms of geographical and genetic proximity, and then to try and account for this variation.

In order to do so, we have picked one language from the Geji, Polci and Zeem clusters (resp. Bùù, Nyámzàk and Cààrí) and two from the Dass and Saya clusters : BÈráázÈ (Baraza) and Zò•ì (Dott) for Dass ; Gùùs (Sigidi) and Zaar for Saya.

SBW Classification 3

Language

TAM systems

All SBW's TAM systems conform to the same structure : conjugation is expressed in a pre-verbal morpheme which includes a person mark. The verb itself is not touched by TAM, with two exceptions : floating tones (or tone progagation) and the occasional use of Verbal Nouns for Continuous and Future. When the subject is nominal, the personal mark is dropped and a variant of the TAM marker appears alone.

What can be noted though is the common existence of a 'zero' TAM (H. Jungraithmayr's 'Grundaspekt'), where only the bare Pro is used to conjugate the verb. We have chosen to name this TAM Aorist when it is opposed to the couple Completive + Incompletive, and Completive when it is opposed to an Incompletive only. We have named Continuous the TAM with a locative structure involving Pro's having the form of subjects of non-verbal predicates, and a Verbal Noun. The Bùù TAM system is characterised by the use of the plural marker -ní suffixed to the Verb and Verbal Noun in the 2P. In the subjunctive, -ní can be used for the 1P and 3P as well.

Bùù (Zaranda)

 Aorist. The inflexion follows the pattern « Pro + Ø ».

 Completive. The inflexion follows the pattern « Pro + á ». With a noun subject, the marker is ap. Achievement of the process (perfect) is expressed by the addition of ¿ì at the end of the Verb Phrase. The Negative Completive is formed with the negative marker wèè at the end of the sentence. ¿ì and wèè are in complementary distribution.

 Continuous. The inflexion follows the pattern « Pro + íí », followed by a Verbal Noun. With a noun subject, the marker is Negation is expressed with yáÑ hÚ at the end of the sentence.

 Completive. The inflexion follows the pattern « Pro + Ø ». It is followed by a floating L tone in 1p and 2p. The plural. -nÈ is suffixed to the verb in 2P. (NB : contrary to Bùù, -nÈ appears in this context only.) Achievement of the process is expressed by the addition of wì at the end of the Verb Phrase.

 Incompletive. It is restricted to the expression of habits and future events. The inflexion follows the pattern « Pro + Ø » and bears a H tone. With noun subjects, the form is bÈ/bà. A bÉ variant appears between the subject pronoun and the verb in the 3S.

 Continuous. The inflexion follows the pattern of a locative construction « Pro + dàn + VN » using the pronominal paradigm of non-verbal predicate subjects.

 Future. The inflexion follows the pattern « Pro + rÉkÈ ». With a noun subject, the future marker is bàr. The negation is formed with ràÑ at the end of the sentence.

Nyámzàk (Langas)

 Completive. The inflexion follows the pattern « Pro + Ø ». Its functional values are those of general past, subjunctive, and default TAM for state verbs (i.g. "love"). Perfect is expressed by the addition of wì at the end of the Verb Phrase.

 Incompletive. The inflexion follows the pattern « Pro + aa », except for 1S where the vowel is short. With a noun subject, its form is dÉ.

 Continuous. The inflexion follows the pattern of a locative construction « Pro + pÉ + VN » using the pronominal paradigm of non-verbal predicate subjects.

 Future I is derived from the Incompletive through apotony : H > L. With a noun subject, its form is dÈ. The 1S /n/ in nàà is probably etymological.

 Future II is a more recent derivation from the Incompletive through the insertion of the rÉgÈ morpheme, as can be seen from the 1S à rÉgÈ. With a noun subject, its form is rÈgÈ.

Cààrí (Danche)

The survey of Cààrí was very short and the data collected is not very reliable. However, since Cààrí is the only surviving language of the Zeem subgroup, the other three (Zeem, Tule and Dokshi) being extinct, we have decided to include it here to give some idea of its TAM system.  Recent Past. The inflexion follows the pattern « Pro + am », followed by a floating Hi tone. With a noun subject, the marker is am.

 Remote Past. The inflexion follows the pattern « Pro + sËn ». With a noun subject, the marker is âsÉn.

 Future. The inflexion follows the pattern « Pro + Ø », followed by a floating Hi tone which differenciates it from the Subjunctive. With a noun subject, the marker is vOE.

Mood

Subjunctive

3P

ÿì ÿìs  Subjunctive. The inflexion follows the pattern « Pro + Ø ». With a noun subject, the markers are gu for the singular and gi for the plural.

 Counterfactual. The inflexion follows the pattern « Pro + s ». With a noun subject, the marker is âs. The mark appears both in the conditional and in the main clauses, and it functions like the Hausa dàà…dàà.

Zaar (Saya)

2.8.1.

Verb classes Zaar has 2 verb classes : Middle (M) and Hi (H), which can be observed in the Completive. The M class includes both 1-and 2-syllable verbs. A variant appears with some TAM (sp. the plural Incompletive) and depends on the [+/-V] nature of the initial consonant of the verb : [+V] = nonimplosive voiced obstruent. For short, we call the Completive form 'simple', and the other one 'complex'. This is probably a complex case of tone propagation from the TAM markers to the verb. The verbs appear with a high tone on their first syllable in the 1 st and 2 nd persons of the singular. In the 3 rd person singular and in the plural, they appear in the complex form.

 Completive. The inflexion follows the pattern « Pro + áá». The form of the verb is the simple form.

 Incompletive. The inflexion follows the pattern « Pro + yáá ». The verb appears in the simple form in the singular, and in the complex form in the plural.

 Punctiliar. The meaning is both punctiliar and iterative. The closest equivalent is the periphrastic Hausa 'yi ta'. The inflexion follows the pattern « Pro + yii ». The verb appears in the simple form.

 Concomitant. The inflexion follows the pattern « Pro + yáá + naa ». The verb keeps the simple form.

 Continuous. The inflexion follows the pattern « Pro + yi + ká + VN ».

Tense

Future  Subjunctive. The inflexion follows the pattern « Pro + low tone ». In the plural, a low tone is added to the 1 st syllable of the verb.

 Conditional. The inflexion follows the pattern « pronoun + yáá ». In the 1 st and 2 nd person singular, the verb takes a high tone on the 1 st syllable. In the other persons, the verb appears in the complex form  Counterfactual. The inflexion follows the pattern « pronoun + yí ». The verb appears in the complex form. The Counterfactual is introduced by yáàn, if and is used in contexts that would be translated by a dàà ... dàà Hausa construction.

Composition

Compound aspects

The three basic aspects (Perfect, Imperfect, Aorist) cannot be combined.  Perfect Punctiliar. The inflexion has the structure : « Pro + áá + yi ». All verbs following the yi punctiliar morpheme have a simple form. This is true for all combinations.

 Imperfect Punctiliar. The inflexion has the structure : « Pro + yàá + yi ».

 Imperfect Continuous. The structure of the inflexion is « Pro + yàá + yi + ká ». yi + ká is realised [yi‹á]. The verb appears in the form of a Verbal Noun.

 Punctiliar Continuous. The structure of the inflexion is « Pro + yii + yi + ká ». yi + ká is realised [yi‹á]. The verb appears in the form of a Verbal Noun. This combination is not common and Zaar speakers perfer to use the Imperfect continuous instead.  Remote Past Perfect. The inflexion has the structure : « Pro + áá + tá ». The verb is in the complex form.

 Remote Past Imperfect. The inflexion has the structure : « Pro + ta + yáá ». The verb appears in the simple form in the singular, and in the complex form in the plural. (cf. Imperfect).

 Remote past Punctiliar. The inflexion has the structure : « Pro + ta + yi ». The verb keeps the simple form. This combination is associated to the ßógoró / Ñgòßíyà dialect by the Lusa speakers. They prefer to use the Recent Past Imperfect.

 Remote Past Continuous. The inflexion has the structure : « Pro + ta + yi + ká ». The verb is in the form of a Verbal Noun.

 Recent Past Perfect. The inflexion follows the pattern « Pro +áá + náá ». The verb appears in the complex form.

 Recent Past Imperfect. The inflexion follows the pattern « Pro + náá + yáá ».

The verb appears in the simple form in the singular, and in the complex form in the plural. (cf. Imperfect).

 Recent past Punctiliar. The inflexion follows the pattern « Pro + náá + yi ». The verbs keep the simple form (cf. Imperfect). This combination is associated to the ßógoró / Ñgòßíyà dialect by the Lusa speakers. They prefer to use the Recent Past Imperfect.

 Recent past Continuous. The inflexion follows the pattern « Pro + náá + yi + ká ». The verb is in the form of a verbal noun.

 Immediate Past Perfect. The inflexion follows the pattern « Pro + áá + yi ». The first syllable of the verb receives a low tone.

 Immediate Past Imperfect. The inflexion follows the pattern « Pro + yii + yáá ».

The verb appears in the simple form in the singular, and in the complex form in the plural. (cf. Imperfect).

 Immediate Past Continuous. The inflexion follows the pattern « Pro + yi + yi + gá ». The verb appears in the form of a verbal noun.

Mood and aspect

Comments

Just looking at the number of paradigms, the situation varies from 3 different paradigms (Zo•i) to 13 or more if combinations are included (Zaar). Zaar and Guus, in this respect, side together against the other languages. From now on, all that will be said concerning Zaar holds for Guus as well. Zaar's main innovations lies in the presence of 3 past tenses, and a great number of aspects.

The diversity can be explained by the development from a basic opposition between a Ø TAM, and an Incompletive. The Incompletive is generally formed with -aa, although in Dir, the contrast between Compl and Incompl is tonal. At this stage, the functional value of the opposition is not differenciated from the point of view of tense, aspect and mood. Each morpheme is a mixture of functional values. The Ø TAM combines narrative, past, perfect and subjunctive functions (cf. Nyamzak Completive). The Incompletive is basically non-past (present and future), progressive and habitual. Depending on the language, either can be used as the default TAM to conjugate state verbs.

Then, new TAM's appear to convey a specification that is absent in the original system. The Continuous is formed through a locative construction. Other TAM's are made through the integration of particles that are prefixed (Subj : Bùù, BÈráázÈ) or more generally suffixed (Zaar pasts ; Buu and Dir future). These restrict the functional value of the original TAM's. The Completive for example can be reduced to the expression of the subjunctive (Zaar), or the narrative and conjugation of state verbs such as 'want/like' (Buu). The Incompletive can be reduced to the expression of the habitual, although the negative Incompletive is generally still used to negate the Future, and sometimes alternates with a specific Negative Future TAM.

Conclusion

In conclusion, I would like to go back on Zaar. The morphological and functional status of Zaar TAM's as compared to the others look very odd. The basic process at work in the evolution of the Baraawaa languages has been pushed a step futher in its case : proliferation of TAM's ; possibility of multiple combinations ; development of past tenses. As was mentioned in the introduction, this goes with the appearance of a third tone level (may be an infra-low tone linked to +V consonants). This can be explained by the history of the Zaars. Their population far exceeds that of the other languages : over 125,000 for Zaar ; less than 600 for those that are still well spoken (i.e. B'oto, Zo•i, Zaranda, Dir, Baraza), a handful for those that are dying (Luri, Zeem, Zumbul, Tulai, Lushi). They have mixed with and incorporated alien populations (Angas, Pyem, Jarawa, Boi, etc.), whose existence is still acknowledged in the indentity of clans. But what is remarkable is that, far from becoming a communication language characterized by its morphological pidgin-like simplification, this has lead to a process of phonological, morphological and functional complexification. A case study for R. Nicolaï's « Nouveau programme » 5 ? In any case, Zaar must definitely be set apart from the other SBW languages.

5 Cf. in this volume, R. Nicolaï, "Contact et genèse: ouvertures et perspectives pour un "Nouveau Programme" de recherche sur l'évolution des langues".

•É. The Negative Incompletive is formed with the discontinuous marker •áá … wèè. Ex : kíí •áá VN

  … wèè (2S), etc.

	 Future. The inflexion follows the pattern « Pro + m », followed by a Verbal
	Noun. With a noun subject, the marker is nÉm. The Negative Future has two
	variants : one is expressed through the Negative Incompletive ; the other uses a
	discontinuous marker

kÉm ... wèè : kí kÉm VN … wèè (

	2S), etc.
	 Subjunctive. The form with a noun subject is tÈ. The paradigm ot the Subjunctive
	is very similar to that of the Aorist. The only difference appears in 1S and 1P. This
	can be explained by the incorporation of tÈ. Its position indicates that it may be a
	complementiser. The Negative Subjunctive has two forms. With a noun subject
	they are (1) S

ßáá Vb and (2) S bÈ Vb … wèè. With pronominal subjects, ßàà and bÈ form new paradigms :

	2.2. Dìr
		COMPL	ICOMPL CONT	FUT
	1S	à H	á	ám dàn	âr kÈ
	2S	kÈ H	kÉ	kÉ dàn	kÉrÉ kÈ
	3S	yàà H	yáá bÉ	yàx dàn yààr kÉ yaa bár kÈ
	1P	mù L	mú -R	mì dàn	mûr kÈ
	2P	kÈ L--nÈ kÉ -R	kÈn dàn kÉrÈ kÈ
	3P	wù H	wú	wú dàn	wûr kÈ
		NEG SBJ (1) NEG SBJ (2)
		1S ßâ		bàà … wèè
		2S ßì		bìì … wèè
		3S ßÈ		bÈ … wèè
		1P ßân Vb-ní	bÈn … wèè
		2P ßì Vb-ní	bìì Vb-ní … wèè
		3P ßÈ Vb-ní	bÈs … wèè

níí. 2.7. Gùùs (Sigidi)

  The inflexion follows the pattern « nÉ + Pro ». It takes a Lo tone and is followed by a floating Hi tone. The Subjunctive is substituted for the Completive in [+Foc] contexts.Negation is expressed with •áá … bà, both for the Completive and the Incompletive. The TAM marekers being the same, the difference is in the predicate : a verb for the Negative Completive, and a VN for the Negative Incompletive. The inflexion follows the pattern « Pro + ap ». With a noun subject, the marker is ap. The inflexion follows the pattern « Pro + kap ». With a noun subject, the marker is kap followed by a floating Hi tone. The inflexion follows the pattern « Pro + ka/ak ». With a noun subject, the marker is ak, followed by a floating Hi tone.

						NEG COMPL & ICOMPL
				1S			ám •áá … bà
				2S			kÉ •áá … bà
				3S			yèè •áá … bà
				1P			mÈ 4 •áá … bà
	2.7.2.	AOR Tense 2P 3P IMMEDIATE PST	COMPL kù •áá … bà ICOMPL kì •áá … bà RECENT PST REMOTE PST	FUT FUT
	1S 2S 1S 2.6. Zò•ì (Dott) mÉ ká maa 2S kaa		áá káá mam +H COMPL ICOMPL FUT àmÉ àká mÉsËn kam +H kÉsËn	mîr kîr mOE +H kOE +H
	3S 3S	ÿaa	tÉ	1S a		táá ÿam +H ám	tÉká máá ÿísËn	tîr ÿi +H
	1P 1P	mÈ màà	2S ú	màà màm +H k y áà	mÈká kíí mÈsËn	m y èérì mÈ +H
	2P 2P	kàà	kÈ	3S Ø		kàà kàm +H t y áà	kÈká Ø kÈsËn	k y èérì kÈ +H
	3P 3P	ÿàà	tÈ	1P mà	tàà ÿàm +H màà	tÈká màn ÿìsËn	t y èérì ÿì +H
	2P kò 2.5. BÈráázÈ (Baraza) 3P tì	k w àà tàà	kòn tìn
		COMPL	ICOMPL	HAB	SUBJ
	1S	á			ámáá	ám tÉ	nàà H-
	2S	kÉ		káá		kÉ tÉ	nÉkÈ H-
	3S	yèè		yàá/yàà H-	yèè tÉ	nÉyèè H-
	1P	mù		màá/màà H-	mù tÉ	nÉmÈ H-
	2P	kù		k w àá/kwàà H-	kù tÉ	nÉkù H-
	3P	kì			k y àá/k y àà H-	kì tÉ	nÉkì H-
	2.7.1.		Aspect
		AOR	COMPL ICOMPL	HAB
	1S	ma	map		mOEkap +H mOEka / mak + H
	2S	ka	kap		kOEkap +H kOEka / kak + H
	3S	ÿa	ÿap		ÿikap +H ÿika / ÿak + H
	1P	mà	màp		mÈkàp +H mÈkà / màk + H
	2P	kà	kàp		kÈkàp +H kÈkà / kàk + H
	3P	ÿà	ÿàp		ÿìkàp +H ÿìkà / ÿàk + H
	 Aorist. The inflexion follows the pattern « Pro + a ». With a noun subject, the
	marker is a.

 Completive. The inflexion follows the pattern « Pro + Ø ». The Perfect is expressed by the addition of ßá at the end of the Verb Phrase. It is in complementary distribution with the negation •áá … bà. Iteration is expressed by the addition of the combination róó ßá at the end of the Verb Phrase.

 Incompletive. The inflexion follows the pattern « Pro + aa». In the Continuous, the same inflexion is followed by a Verbal Noun instead : ¨àù, "slaughter" > ¨àù-tÈ, "slaughter-ing".

 Habitual. The inflexion follows the pattern « Pro + tÉ».

 Subjunctive.

 Completive. The inflexion follows the pattern « Pro + Ø ». The Perfect is a combination of the Completive + CaK at the end of the Verb Phrase. The habitual is marked by •a inserted between the Completive subject pronouns and the verb.

 Incompletive. Except for 1S, the inflexion follows the pattern « Pro + àà ». It is followed by a VN. With a noun subject, the Incompletive mark is à.

 Future. The inflexion follows the pattern « Pro + VV » in the singular and « Pro + n » in the plural. With a noun subject, the Future mark is  Completive.

 Incompletive.

 Habitual.

 Immediate Past. The inflexion follows the pattern « Pro + aa ». With a noun subject, the marker is aa.

  The inflexion follows the pattern « Pro + Ø ». It has a narrative functional value.

	Verb class		H	M ; +V	M-2syl ; +V	M ; -V	M-2syl ; -V
	Simple (=lexical)	súú	bwaa	daambár	taar	tuurá
	Complex (=modified) sùú (R) bwàà (L) dààmbár (LH) táàr (F) túùrá (FH)
			love	choose	disturb		clear	push
	2.8.2.	Aspect
	Aorist Compl Icompl 1S mOE máá m y áá 2S kOE káá k y áá	Punctiliar miyii/myii kiyii/kyii	Concomitant Continuous myáánaa mi‹á / m y aa kyáánaa ki‹á / k y aa
	3S á 1P mÉ 2P kÉ	àà màà kàà	ÿáá m y àá k y àá	yii mààyi kààyi		ÿáánaa m y àánaa k y àánaa	ÿi‹á / ÿaa mì‹á / m y àá kì‹á / k y àá
	3P tÉ	tàà	ÿàá	tààyi		ÿàánaa	ÿì‹á / ÿàá
	 Aorist.

  The inflexion follows the pattern « Pro + a ». The verb appears in the simple form in the singular, and the complex form in the plural.

	2.8.4.	Mood
		Subjunctive Conditional Counterfactual 1s mÈ m y áá mí/míyí 2s àà k y áá kí/kíyí
		3s tÈ 1p mÈ +L 2p àà + L	yáá m y aá k y aá		ÿí/yí/ÿíyí/tíyí mì/mìyì kì/kìyì
		3p tÈ + L	yaá/ÿaá	ÿì/ÿìyì/tìyì
		Remote Past Recent Past Immediate past
	1S	ma	mOEtá	mOEnáá	míí
	2S	ka	kOEtá		kOEnáá	kíí
	3S	wò	áta		ánáa	ájí
	1P	má	mÈtà	mÈnàà	mìì
	2P	ká	kÈtà		kÈnàà	kìì
	3P	tá	tÈtà		tÈnàa	ÿìì
	 Future.  Remote Past. The Remote Past refers to events which occurred more than two
	days ago. The inflexion follows the pattern « Pro + ta ». The verb keeps the
	simple form.
	 Recent Past. The Recent Past refers to events which occurred on the previous
	day. The inflexion follows the pattern « Pro + náá ». The verb keeps the simple
	form.

 Immediate past. The Immediate Past refers to events which occurred earlier the same day. The inflexion follows the pattern « Pro + yi ». The first syllable of the verb receives a low tone.

  The inflexion has the structure : « Pro + a + yi ». The verb keeps the simple form. Future Continuous. The inflexion has the structure : « Pro + a + yi + ká ». The verb is in the form of a Verbal Noun.

	Tense and aspect  Future Punctiliar.
	Future	Future	Remote Past	Remote Past	Remote past	Remote Past	Recent Past
	Punctiliar	Continuous	Perfect	Imperfect	Punctiliar		Continuous	Perfect
	1s mayi	mayi‹á	máátá	mOEtáyáá	mOEtáyi		mOEtáyi‹á	máánáá
	2s kayi	kayi‹á	káátá	kOEtáyáá	kOEtáyi		kOEtáyi‹á	káánáá
	3s wòyi	wòyi‹á	ààtá	átayáá	átáyi		átáyi‹á	àànáá
	1p máyi	máyi‹á	mààtá	mÈtàyáá	mÈtàyi		mÈtàyi‹á	màànáá
	2p káyi	káyi‹á	kààtá	kÈtàyáá	kÈtàyi		kÈtàyi‹á	kàànáá
	3p táyi	táyi‹á	tààtá	tÈtàyáá	tÈtàyi		tÈtàyi‹á	tàànáá
	Recent Past	Recent Past	Recent past	Recent past	Immediate Past	Immediate Past
	Perfect	Imperfect	Punctiliar	Continuous	Perfect	Continous
	1s máánáá	mOEnááyáá	mOEnááyi	mOEnááyi‹á	mááyí	mííyì‹á
	2s káánáá	kOEnááyáá	kOEnááyi	kOEnááyi‹á	kááyí	kííyì‹á
	3s àànáá	ánáayáá	ánááyi	ánááyi‹á		ààyí		áyyì‹á
	1p màànáá	mÈnààyáá	mÈnààyi	mÈnààyi‹á	mààyí	mììyì‹á
	2p kàànáá	kÈnààyáá	kÈnààyi	kÈnààyi‹á	kààyí	kììyì‹á
	3p tàànáá	tÈnààyáá	tÈnààyi	tÈnààyi‹á	tààyí	ÿììyì‹á

Shimizu 1978, p. 13

The variation u / OE for the vowel of the 1P inflexion is a phonological phenomenon quite common in SBW languages. i/OE is a weak vowel that replaces the others in conditions that need to be further studied.


Charles Kervrann

Jérôme Boulanger

Pierrick Coupé

Bayesian Non-Local Means Filter, Image Redundancy and Adaptive Dictionaries for Noise Removal

Partial Differential equations (PDE), wavelets-based methods and neighborhood filters were proposed as locally adaptive machines for noise removal. Recently, Buades, Coll and Morel proposed the Non-Local (NL-) means filter for image denoising. This method replaces a noisy pixel by the weighted average of other image pixels with weights reflecting the similarity between local neighborhoods of the pixel being processed and the other pixels. The NL-means filter was proposed as an intuitive neighborhood filter but theoretical connections to diffusion and non-parametric estimation approaches are also given by the authors.

In this paper we propose another bridge, and show that the NL-means filter also emerges from the Bayesian approach with new arguments. Based on this observation, we show how the performance of this filter can be significantly improved by introducing adaptive local dictionaries and a new statistical distance measure to compare patches. The new Bayesian NL-means filter is better parametrized and the amount of smoothing is directly determined by the noise variance (estimated from image data) given the patch size. Experimental results are given for real images with artificial Gaussian noise added, and for images with real imagedependent noise. R 2 G a (t)|z(x + t)-z(y + t)| 2 dt := exp -z(x)-z(y) 2 2,a h 2

(2)

Introduction

Denoising (or restoration) is still a widely studied and an unsolved problem in image processing. Many methods have been suggested in the literature, and a recent outstanding review of them can be found in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. Some of the more advanced methods are based on PDEs [START_REF] Mumford | Optimal approximations by piecewise smooth functions and variational problems[END_REF][START_REF] Perona | Scale space and edge detection using anisotropic diffusion[END_REF][START_REF] Rudin | Nonlinear Total Variation based noise removal algorithms[END_REF][START_REF] Weickert | Anisotropic Diffusion in Image Processing[END_REF] and aim at preserving the image details and local geometries while removing the undesirable noise; in general, an initial image is progressively approximated by filtered versions which are smoother or simpler in some sense. Other methods incorporate a neighborhood of the pixel under consideration and perform some kind of averaging on the gray values. One of the earliest examples for such filters has been presented by Lee [START_REF] Lee | Digital image smoothing and the sigma filter[END_REF] and a recent evolution is the so-called bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] with theoretical connections to local mode filtering [START_REF] Van De Weijer | Local mode filtering[END_REF], non-linear diffusion [START_REF] Barash | A Common framework for nonlinear diffusion, adaptive smoothing, bilateral filtering and mean shift[END_REF][START_REF] Buades | The staircasing effect in neighborhood filters and its solution[END_REF] and nonlocal regularization approaches [START_REF] Mrazek | On robust estimation and smoothing with spatial and tonal kernels[END_REF][START_REF] Elad | On the bilateral filter and ways to improve it[END_REF].

However, natural images often contain many structured patterns which can be misclassified either as details to be preserved or noise, when usual neighborhood filters are applied. Very recently, the so-called NL-means filter has been proposed by Buades et al. [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] that can deal with such a "structured" noise: for a given pixel, the restored gray value is obtained by the weighted average of the gray values of all pixels in the image; each weight is proportional to the similarity between the local neighborhood of the pixel being processed and the neighborhood corresponding to the other image pixels. A similar patch-based regularization approach based on the key idea of iteratively increasing a window at each pixel and adaptively weighting input data has been also proposed in [START_REF] Kervrann | Unsupervised patch-based image regularization and representation[END_REF] with excellent results on a commonly-used image database [START_REF] Portilla | Image denoising using scale mixtures of Gaussians in the wavelet domain[END_REF]. The success of the NL-means filter (see [START_REF] Kervrann | Unsupervised patch-based image regularization and representation[END_REF][START_REF] Kinderman | Deblurring and denoising of images by nonlocal functionals[END_REF][START_REF] Mahmoudi | Fast image and video denoising via nonlocal means of similar neighborhoods[END_REF][START_REF] Lukin | A multiresolution approach for improving quality of image denoising algorithms[END_REF][START_REF] Coupé | Fast non-local means denoising for 3D MR images[END_REF][START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF][START_REF] Azzabou | Random walks, constrained multiple hypothesis and image enhancement[END_REF]), inspired by the Efros and Leung's exemplar-based approach for texture synthesis [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF], is mainly related to image redundancy. A similar idea was early and independently proposed for Gaussian [START_REF] De Bonet | Noise reduction through detection of signal redundancy[END_REF] and impulse noise [START_REF] Wang | Restoration of impulse noise corrupted images using long-range correlation[END_REF][START_REF] Zhang | Image information restoration based on long-range correlation[END_REF] removal in images, and more recently for image inpainting [START_REF] Criminisi | Region filling and object removal by exemplar-based inpainting[END_REF]. Similarities between image patches have been also used in the early 90's for texture segmentation [START_REF] Geman | Boundary detection by constrained optimization[END_REF][START_REF] Kervrann | A Markov random field model-based approach to unsupervised texture segmentation using local and global spatial statistics[END_REF]. More recently, other recent denoising methods demonstrated that representations based on local image patches outperform the best denoising methods of the state-of-the-art [START_REF] Awate | Higher-order image statistics for unsupervised, information-theoretic, adaptive, image filtering[END_REF][START_REF] Kervrann | Unsupervised patch-based image regularization and representation[END_REF][START_REF] Dabov | Image denoising with block-matching and 3D filtering[END_REF][START_REF] Hirakawa | Image denoising using total least squares[END_REF][START_REF] Elad | Image denoising via learned dictionaries and sparse representation[END_REF]; in [START_REF] Roth | Fields of experts: a framework for learning image priors with applications[END_REF], patch-based Markov random field (MRF) models and learning techniques have been introduced to capture non-local pairwise interactions, and were successfully applied in image denoising.

In this paper, we present a new Bayesian motivation for the NL-means filter briefly described in Section 2. In Section 3, we adopt a blockwise (vectorial) representation and introduce spatially adaptive dictionaries in the modeling for better contrast restoration (Section 4). Using the proposed Bayesian framework, we revise the usual Euclidean distance used for patch comparison, yielding to a filter which is better parametrized, and with a higher performance. In Section 4, we also show how smooth parts in the image can be are better recovered if the restored image is "recycled" once. Experimental results on artificial and real images are presented in Section 5, and the performance is very close to the most competitive denoising methods. It is worth noting that the proposed modeling framework is general and could be used to restore images corrupted by non-Gaussian noises in applications such as biomedical imaging (microscopy, ultrasound imagery, ...) or remote sensing.

Image denoising with the NL-means filter

In this section, a brief overview of the NL-means method introduced in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] is presented. Consider a gray-scale image z = (z(x)) x∈Ω defined over a bounded domain Ω ⊂ R 2 , (which is usually a rectangle) and z(x) ∈ R + is the noisy observed intensity at pixel x ∈ Ω. The NL-means filter is defined as

NL z(x) = 1 C(x) y∈Ω w(x, y) z(y) (1)
where NL z(x) at pixel x is the weighted average of all gray values in the image and C(x) is a normalizing factor, i.e. C(x) = y∈Ω w(x, y). The weights w(x, y) defined as

w(x, y) = exp - 1 h 2
express the amount of similarity between the vectorized image patches z(x) and z(y) (or neighborhoods) of each pair of pixels x and y involved in the computation. The decay parameter h ≈ 12σ acts as a filtering parameter. A Gaussian kernel G a (•) of standard deviation a is used to take into account the distance between the central pixel and other pixels in the patch. In (2), the pixel intensities of a √ n × √ n square neighborhood B(x) centered at pixel x, are taken and reordered lexicographically to form a n-dimensional vector z(x) := (z(x k ), x k ∈ B(x)) ∈ R n . In [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], it was shown that 7 × 7 patches are able to take care of the local geometries and textures seen in the image while removing undesirable distortions. The range of the search space in the NL-means algorithm can be as large as the whole image. In practice, it is necessary to reduce the total number of computed weights -|Ω| weights for each pixel -to improve the performance of the algorithm. This can be achieved by selecting patches in a semi-local neighborhood corresponding to a search window ∆(x) of 21×21 pixels. The NL-means filter we will now consider, is then defined as

NL h z(x) = 1 C(x) y∈∆(x) e -z(x)-z(y) 2 /h 2 z(y), C(x) = y∈∆(x) e -z(x)-z(y) 2 /h 2 (3)
where, for the sake of simplicity, • denotes the usual ℓ 2 -norm. In practice, it is just required to set the √ n × √ n patch size, the search space ∆(x) and the filtering parameter h. Buades et al. showed that this filter substantially outperforms the bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] and other iterative approaches [START_REF] Rudin | Nonlinear Total Variation based noise removal algorithms[END_REF].

Since, several accelerated versions of this filter have been proposed [START_REF] Mahmoudi | Fast image and video denoising via nonlocal means of similar neighborhoods[END_REF][START_REF] Coupé | Fast non-local means denoising for 3D MR images[END_REF]. In [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], Buades et al. recommended the vectorial (or block-based) NL-means filter defined as

NL h z(x) = 1 C(x) y∈∆(x) e -z(x)-z(y) 2 /h 2 z(y), C(x) = y∈∆(x) e -z(x)-z(y) 2 /h 2 , (4)
which amounts to simultaneously restore pixels of a whole patch z(x) from nearby patches. The restored value at a pixel x is finally obtained by averaging the different estimators available at that position [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. This filter can been considered as a fast implementation of the NL-means filter, especially if the blocks are picked up on a subsampled grid of pixels. In this paper, the proposed filter is inspired by this intuitive vectorial NL-means filter [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], but also by other recent accelerated versions [START_REF] Mahmoudi | Fast image and video denoising via nonlocal means of similar neighborhoods[END_REF][START_REF] Coupé | Fast non-local means denoising for 3D MR images[END_REF], as explained in the next sections. The related Bayesian framework enables to establish the relationships between these algorithms, to justify some underlying statistical assumptions and to give keys to set the control parameters of the NL-means filter. It is worth noting that this framework could be also used to remove noise in applications for which the noise distribution is assumed to be known and non-Gaussian.

A Bayesian risk formulation

In a probabilistic setting, the image denoising problem is usually solved in a discrete setting. The estimation problem is then to guess a n-dimensional patch u(x) from its noisy version z(x) observed at point x. Typically, the unknown vectorized image patch u(x) is defined as u(x) := (u(x k ), x k ∈ B(x)) ∈ R n where B(x) defines the √ n× √ n square neighborhood of point x and the pixels in u(x) are ordered lexicographically. Let us suppose now that u(x) is unknown but we can observe z(x) = f (u(x), v(x)) where z(x) := (z(x k ), x k ∈ B(x)), v(x) represents noise and f (•) is a linear or a non-linear function related to the image formation process. The noise v(x) is a random vector which components are iid, and u(x) is considered as a stochastic vector with some unknown probability distribution function (pdf).

Conditional mean estimator

To compute the optimal Bayesian estimator for the vector u(x), it is necessary to define an appropriate loss function L(u(x), u(x)) which measures the loss associated with choosing an estimator u(x) when the true vector is u(x). The optimal estimator u opt (x) is found by minimizing the posterior expected loss

E[L(u(x), u(x))] = u(x)∈Λ L(u(x), u(x)) p(u(x)|z(x)),
taken with respect to the posterior distribution p(u(x)|z(x)) and Λ denotes the large space of all configurations of u(x) (e.g |Λ| = 256

n if u(x) ∈ {0, • • • , 255}). The loss function used in most cases is L(u(x), u(x)) = 1 -δ(u(x), u(x)) where the δ function equals 1 if u(x) = u(x) and 0 otherwise. Minimizing E[L(u(x), u(x))] is then equivalent to choose uopt(x) = arg max p(u(x)|z(x))
, with the motivation that this should correspond to the most likely vector given the observation z(x). However, this loss function L(u(x), u(x)) may not be the most appropriate since it assigns 0 cost only to the perfect solution and unit cost to all other estimators. Another thought would be to use a cost function that depends on the number of pixels that are in error such as 2 . Assuming this quadratic loss function, the optimal Bayesian estimator is then

L(u(x), u(x)) = u(x) -u(x)
uopt(x) = arg min b u(x) u(x) u(x) -u(x) 2 p(u(x)|z(x)) = u(x) u(x) p(u(x)|z(x)).
Referred as the conditional mean estimator, uopt(x) can be also written as

uopt(x) = u(x) u(x) p(u(x), z(x)) p(z(x)) = u(x) u(x)p(z(x)|u(x))p(u(x)) u(x) p(z(x)|u(x))p(u(x)) (5)
by using the Bayes' and marginalization rules, and p(z(x)|u(x)) and p(u(x)) respectively denote the distribution of z(x)|u(x) and the prior distribution of u(x).

Bayesian filter and image redundancy Ideally, we would like to know the pdfs p(z(x)|u(x)) and p(u(x)) to compute uopt(x) for each point x in the image from a large number of "repeated" observations (i.e. images). Unfortunately, we have only one image at our disposal, meaning that we have to adopt another way of estimating these pdfs. Due to the fact that the pdfs p(z(x)|u(x)) and p(u(x)) cannot be obtained from a number of observations at the same point x, we choose to use the observations at a number of neighboring points taken in a semi-local neighborhood (or window) ∆(x). The window ∆(x) needs to be not too large since the remote samples are likely less significant and can originate from other spatial contexts. We then assume that this set of nearby samples may be considered as a set of samples from p(u(x)|z(x)).

More formally, we suppose that p(u(x)|z(x)) is unknown, but we have a set {u(x 1 ),

u(x 2 ), • • • , u(x N (x) )} of N (x) posterior samples taken in ∆(x).
In what follows, |∆(x)| is fixed for all the pixels but the size N (x) ≤ |∆(x)| is spatially varying since irrelevant and unlikely samples in ∆(x) are preliminarily discarded. From this set, we start by examining the prior distribution p(u(x)). A first natural ideal would be to introduce MRFs and Gibbs distributions to capture interactions between pixels in the image patch, but the MRF framework involves the computationally intensive estimation of additional hyperparameters which must be likely adapted to each spatial position. Due to the huge domain space Λ, a computational alternative is then to set p(u(x)) to uniform, i.e. p(u(x)) = 1/N (x). This means there is no preference to choose a vector u(x i ) in the set {u(x 1 ), • • • , u(x N (x) )} assumed to be composed of N (x) preliminarily selected "similar" patches. Then, we have the following approximations (see [START_REF] Godtliebsen | A nonlinear Gaussian filter applied to images with discontinuities[END_REF]):

1 N (xi) N(x i ) X j=1 u(xj)p(z(xi)|u(xj)) P → X u(x) u(x)p(z(x)|u(x))p(u(x)), 1 N (xi) N(x i ) X j=1 p(z(xi)|u(xj )) P → X u(x) p(z(x)|u(x))p(u(x)),
and we can propose a reasonable estimator u N (x i ) for uopt(x):

u N (x i ) = 1 N (xi) N (xi) j=1 u(x j )p(z(x i )|u(x j )) 1 N (xi) N (xi) j=1 p(z(x i )|u(x j )) . (6)
Nevertheless, we do not have the set {u(x 1 ),

• • • , u(x N (x) )}, but only a spatially vary- ing dictionary D(x) = {z(x 1 ), • • • , z(x N (x)
)} composed of noisy versions is available.

A way to solve this problem will be then to substitute z(x j ) to u(x j ) in ( 6) as shown in the next section. In a second step, this estimator will be refined by substituting the "aggregated" estimator computed from u N (x j ) (see below) to u(x j ). Indeed, the restored patch at pixel x j is a better approximation of u(x j ) than the noisy patch z(x j ) used as a "pilot" estimator, and the performance will be improved at the second iteration.

Aggregation of estimators

The estimator (6) requires spatially sliding windows over the whole image for image reconstruction. Hence, a set of L (constant for uniform image sub-sampling) concurrent scalar values u N,1

(x i ), u N,2 (x i ), • • • , u N,L (x i
) is calculated for the same pixel x i due to the overlapping between patches. This set of competing estimators must be fused or aggregated into the one final estimator u(x i ) at pixel x i . Actually, when no estimator is a "clear winner", one may prefer to combine the different estimators u N,1

(x i ), u N,2 (x i ), • • • , u N,L (x i
) and a natural approach, well-grounded in statistics [START_REF] Bunea | Sequential procedures for aggregating arbitrary estimators of a conditional mean (under revision)[END_REF], is to use a convex or linear combination of estimators [START_REF] Katkovnik | Adaptive window size image denoising based on intersection of confidence intervals (ICI) rule[END_REF][START_REF] Dabov | Image denoising with block-matching and 3D filtering[END_REF]. Here, our aggregate estimator u N (x i ) is simply defined as the average of competing estimators:

u N (x i ) = 1 L L l=1 u N,l (x i ). (7)
In practice, patches are picked up on a sub-sampled grid (e.g. factor 3) of pixels to speed up the algorithmic procedure (e.g. factor 8), while preserving a good visual quality.

As explained before, to compute u N (x i ), we first substitute z(x) to u(x) in ( 6). This yields the following estimator

1 N (xi) N (xi) j=1 p(z(x i )|z(x j ))z(x j ) 1 N (xi) N (xi) j=1 p(z(x i )|z(x j )) ≈ u N (x i ) (8)
which can be computationally calculated provided the pdfs are known. It is worth noting that p(z(x i )|z(x j )) is not well defined if x i = x j and it could be recommended to set p(z(x i )|u(x i )) ≈ max xj =xi p(z(x i )|z(x j )) in (8) (see [START_REF] Kinderman | Deblurring and denoising of images by nonlocal functionals[END_REF]). The central data point involved in the computation of its own average is then re-weighted to get the higher weight. Actually, the probability to detect an exact copy of z(x i ) corrupted by noise in the neighborhood tends to 0 because the space of √ n × √ n patches is huge and, to be consistent, it is necessary to limit the influence of the central patch.

In the remainder of this section, we shall now consider the usual image model

z(x) = u(x) + v(x) (9)
where v(x) is an additive white Gaussian noise with variance σ 2 . We will further assume that the likelihood can be factorized as p(z(x i )|z(x j )) = n k=1 p(z(x i,k )|z(x j,k )) with x i,k ∈ B(x i ) and x j,k ∈ B(x j ). It follows that z(x i )|z(x j ) follows a multivariate normal distribution, z(x i )|z(x j ) ∼ N (z(x j ), σ 2 I n ) where I n is the n-dimensional identity matrix. From ( 8), the filter adapted to white Gaussian noise is then given by

1 C(x i ) N (xi) j=1 e -z(xi)-z(xj ) 2 /(2σ 2 ) z(x j ) with C(x i ) = N (xi) j=1 e -z(xi)-z(xj) 2 /(2σ 2 ) .( 10
)
If we arbitrarily set N (x i ) = N to a constant value and h 2 = 2σ 2 , this filter is nothing else than the vectorial NL-means filter given in (4). However, it is recommended to set h ≈ 12σ in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] to produce satisfying denoising results. In our experiments, it is also confirmed that h ≈ 5σ is good choice if we use (3) and (4) for denoising. Actually, the filtering parameter h is actually set to a higher value than the expected value √ 2σ in practical imaging. In the next sections, we shall see how this parameter can be better interpreted and theoretically estimated.

Spatially adaptive dictionaries

The filter [START_REF] Dabov | Image denoising with block-matching and 3D filtering[END_REF] can be refined if the adaptive dictionary D(x i ) = {z(x 1 ), • • • , z(x N (xi) )} around x i is reliably obtained using an off-line procedure. Since D(x i ) is assumed to be composed of noisy versions of the more likely samples from the posterior distribution p(u(x i )|z(x i )), the irrelevant image patches in ∆(x i ) must be discarded in a preliminary step. Consequently, the size N (x i ) is adaptive according to local spatial contexts and a simple way to detect these unwanted samples can be based on local statistical features between images patches. In our experiments, we have consider two basic features, that is the mean m(z(x)) = n -1 n k=1 z(x k ) and the variance var(z(

x)) = n -1 n k=1 (z(x k ) -m(z(x))) 2 of a vectorized patch z(x) := (z(x k ), x k ∈ B(x)).
Intuitively, z(x j ) will be discarded from the local dictionary D(x i ) if the mean m(z(x j )) is too 'far" from the mean m(z(x i )) when z(x j ) and z(x i ) are compared. More formally, if |m(z(x j ))m(z(x i ))| > λ α σ/ √ n, where λ α ∈ R + is chosen as a quantile of the standard normal distribution, the hypothesis that the two patches belong to the same "population" is rejected. Hence, setting λ α = 3 given P(|m(z(

x j )) - m(z(x i ))| ≤ λ α σ/ √ n) = 1 -α, yields to α = 2(1 -Φ(λ α / √ 2
)) = 0.034 where Φ means the Laplace distribution.

Similarly, the variance var(z(x j )) is expected to be close to the variance var(z(x i )) for the central patch z(x i ). A F -test4 is used and the ratio F = is compared to a threshold T β,n-1 to determine if the value falls within the zone of acceptance of the hypothesis that the variances are equal. The threshold T β,n-1 is the critical value for the F -distribution with n -1 degrees of freedom for each patch and a significance level β. Typically, when 7 × 7 patches are compared, we have P(F > T 0.05,48 = 1.6) = 0.05. If the ratio F exceeds the value T β,n-1 , the sample z(x j ) is discarded from the local dictionary D(x i ).

This formal description is related to the intuitive approach proposed in [START_REF] Mahmoudi | Fast image and video denoising via nonlocal means of similar neighborhoods[END_REF][START_REF] Coupé | Fast non-local means denoising for 3D MR images[END_REF][START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF] to improve the performance of the NL-means filter.

New statistical distance measure for patch comparison In this section, we propose to revise the distance used for patch comparison, yielding to a NL-means filter which is better parametrized. In ( 3) and ( 4), it is implicitly assumed that z(x i )|z(x j ) ∼ N (z(x j ), 1 2 h 2 I n ). Actually, this hypothesis is valid only for non-overlapping and statistically independent patches, but most of patches overlapped in ∆(x) since ∆(x) is not so large (e.g 21 × 21 pixels). At the opposite, if z(x j ) is horizontally (or vertically) shifted by one pixel from the location of z(x j ), z(x i )|z(x j ) is expected to follow a multivariate Laplacian distribution. However, this statistical hypothesis does not hold true for arbitrary locations of overlapping patches in ∆(x). The adjustment of the decay parameter h ≈ 5σ in (3) to a value higher that the expected value √ 2σ is probably related to the fact that the two compared patches are not independent. Note that some pixel values are in common in the two vectors but at different locations.

Hence, p(z(x i )|z(x j )) must be carefully examined and we propose the following definition for the likelihood: p(z(x i )|z(x j )) ∝ e -φ( z(xi)-z(xj) ) . Typically, we can choose φ(t) = t 2 or φ(t) = |t| (or a scaled version) to compare patches. Here, we examine the distribution of z(x i )z(x j ) from the local dictionary D(x i ) to determine φ. First, it can be observed that E[ z(x i )z(x j ) ] is non-zero in most cases and the probability to find an exact copy of z(x i ) in ∆(x) tends to 0, especially if ∆(x) is large. The maximum of the assumed zero-mean multivariate Gaussian distribution in (3) should be then "shifted" to E[ z(x i )-z(x j ) ]. However, this training step could be hard in practice since it must adapted to each spatial position, and we propose instead to use asymptotic results. Actually, we have already assumed (z(x i,k )-z(x j,k )) ∼ N (0, 2σ 2 ) when two pixels taken in z(x i ) and z(x j ) ∈ D(x i ) are compared. Hence, the normalized distance dist(z(x i ), z(x j )) = z(x i )z(x j ) 2 /(2σ 2 ) follows a chi-square χ 2 n distribution with n degrees of freedom. For n large (n ≥ 25), it can be proved that 2 dist(z(x i ), z(x j )) is approximately normally distributed with mean √ 2n -1 and unit variance:

p 2 dist(z(x i ), z(x j ) ∝ exp - 1 2 2 dist(z(x i ), z(x j )) - √ 2n -1 2 (11) ∝ exp - z(x i ) -z(x j ) 2 2σ 2 - z(x i ) -z(x j ) (σ/ √ 2n -1) + 2n -1 2 .
Accordingly, we define the likelihood as p(z(x i )|z(x j )) ∝ exp -φ( z(x i )z(x j ) ) and choose φ(t) = at 2 -b|t| + c with a = 1/(2σ 2 ), b = √ 2n -1/σ and c = (2n -1)/2 depending only on the patch size n and the noise variance σ 2 . From our experiments, it was confirmed that no additional adjustment parameter is necessary provided the noise variance is robustly estimated and the performance is maximum for the true noise variance as expected. Figure 1 (bottom-right) shows the functions e -t 2 /h 2 and e -(|t|/σ- √ 2n-1) 2 /2 by setting n = 49, σ 2 = 1 and h = 5σ, and then illustrates how data points are currently weighted when the original NL-means filter and the so-called Adaptive NL-means filter (ANL) defined as

ANL σ,n z(x i ) = N (xi) j=1 exp - 1 2 z(x i ) -z(x j ) σ - √ 2n -1 2 z(x j ) N (xi) j=1 exp - 1 2 z(x i ) -z(x j ) σ - √ 2n -1 2 . ( 12
)
where N (x i ) = #{D(x i )}, are applied. Note that the data at point x i should participate significantly to the weighted average. Accordingly, since p(z(x i )|z(x j )) tends to 0 when x i = x j in ( 12), we arbitrarily decide to set p(z

(x i )|z(x i )) ≈ max xj =xi p(z(x i )|z(x j ))
as explained before (see also [START_REF] Kinderman | Deblurring and denoising of images by nonlocal functionals[END_REF]).

Bayesian NL-means filter and plugin estimator In the previous sections, z(x j ) was substituted to u(x j ) in (6) to give [START_REF] Criminisi | Region filling and object removal by exemplar-based inpainting[END_REF] and further [START_REF] Elad | On the bilateral filter and ways to improve it[END_REF]. Now, we are free to substitute the vector u ANL (x j ) of aggregated estimators (computed from the set of restored blocks {ANL σ,n z(x j )} , see [START_REF] Coupé | Fast non-local means denoising for 3D MR images[END_REF]) to u(x j ). This plugin ANL estimator defined as

ANL σ,n u ANL (x i ) = N (xi) j=1 exp - 1 2 2 z(x i ) -u ANL (x j ) σ - √ 2n -1 2 u ANL (x j ) N (xi) j=1 exp - 1 2 2 z(x i ) -u ANL (x j ) σ - √ 2n -1 2 (13)
is expected to improve the restored image since u ANL (x j ) is a better approximation of u(x j ) than z(x j ). In ( 13) the restored image is recycled but the weights is a rescaled function (theoretically by a factor 2) of the distance between the "pilot" estimator u ANL (x j ) given by ( 12) and the input vector z(x i ). The estimators are finally aggregated to produce the final restored image (see [START_REF] Coupé | Fast non-local means denoising for 3D MR images[END_REF]).

In this section, we evaluate the performance of different versions of the Bayesian filter and the original NL-means filter using the peak signal-to-noise ratio (PSNR in db) defined as PSNR = 10 log 10 (255 2 /MSE) with MSE = |Ω| -1 x∈Ω (z 0 (x)u(x)) 2 where z 0 is the noise-free original image. We also use the "method noise" described in [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF] which consists in computing the difference between the input noisy image and its denoised version. The NL-means filter (see ( 3)) was applied with h = 5σ and our experiments have been performed with 7 × 7 patches and 15 × 15 search windows, corresponding to the best visual results and the best PSNR values. For all the presented results, we set T 0.05,n-1 = 1.6 and λ 0.034 = 3 to build spatially adaptive dictionaries.

The potential of the estimation method is mainly illustrated with the 512×512 Lena and Barbara images corrupted by an additive white-Gaussian noise (WGN) (PSNR = 22.13 db, σ = 20). We compared the original NL-means algorithm with the proposed modifications and Fig. 1 shows the denoising results using several filters (n = 7 × 7 and N = 15 × 15): i) the NL-means filter NL h z when the similarity is only measured by the Euclidean distance (see ( 3)); ii) the vectorial NL-means filter NL h z with sliding blocks but no spatial sub-sampling (see ( 4)); iii) our Adaptive NL-means filter ANL σ,n z which includes adaptive dictionaries and a similarity measure based on the new distance (see ( 12)); iv) the plugin Adaptive NL-means filter ANL σ,n u ANL (see [START_REF] Elad | Image denoising via learned dictionaries and sparse representation[END_REF]). In most cases, the PSNR values are slightly affected if a spatial sub-sampling (by a factor 3) is used in the implementation, but the time computing is drastically reduced (speed is 8 times less than before): the implementation of the fast Adaptive NL-means filter took 10 sec on a single-CPU PC 2.0 Ghz running linux, and the full Adaptive NL-means filter (with no spatial sub-sampling) took 75 sec for denoising a 512 × 512 image (see table in Fig. 1). In these practical examples, the use of spatially adaptive dictionaries enables to enhance contrasts. Note that the residual noise in flat zones is more reduced with no additional blur, when ANL σ,n u ANL is applied. In Fig. 2, we modified the estimation of noise variance to assess the sensitivity of this parameter on filtering results. In general, the PSNR value is maximum with the true value (e.g. σ = 20 in Fig. 2) but decreases if this value is under-estimated or over-estimated. In Fig. 2, the estimated noise component is similar to a simulated white Gaussian noise but contains few geometric structures if we under-estimate or over-estimate the noise variance. To evaluate the performance of those filters, we reported the PSNR values for different versions of the NL-means filter. In table 1, the numerical results are improved using our filter, with performance very close to competitive patch-based denoising methods. Note that the best results (PSNR values) were recently obtained by filtering in 3D transform domain and combining sliding-window transform processing with block-matching [START_REF] Dabov | Image denoising with block-matching and 3D filtering[END_REF].

In the second part of experiments, we have applied ANL σ,n u ANL to restore an old photograph (Fig. 3 -left column); in that case, the noise variance is estimated from image data (see [START_REF] Kervrann | Unsupervised patch-based image regularization and representation[END_REF]). Nevertheless, in real digital imaging, images are better described by the following model z(x) = u(x)+u γ (x) ε(x) where the sensor noise u γ (x) ε(x) is defined as a function of u(x) and ε(x) is a zero-mean Gaussian noise of variance σ 2 . Accordingly, the noise in bright regions is higher than in dark regions (see Fig. 2 -second column). From experiments on real digital images [START_REF] Faraji | CCD noise removal in digital images[END_REF], it was confirmed that γ ≈ 0.5 (γ = 0 corresponds to WGN in the previous experiments). Accordingly normalized distance as follows: dist(z(x i ), z(x j )) = z(x i )z(x j ) 2 /(2σ 2 z(x j )) with σ ∈ [1. [START_REF] Buades | The staircasing effect in neighborhood filters and its solution[END_REF][START_REF] Barash | A Common framework for nonlinear diffusion, adaptive smoothing, bilateral filtering and mean shift[END_REF]. Moreover, this model z(x) = u(x) + u(x) ε(x) has already been considered to denoise log-compressed ultrasonic images [START_REF] Loupas | An adaptive weighted median filter for speckle suppression in medical ultrasonic images[END_REF]. Preliminary results of ANL σ,n u ANL with this model is shown in Fig. 3, when applied to two log-compressed ultrasonic images and a cryo-Electronic Microscopy image (cross-section of a microtubule (10-30 nm)) where brights areas are smoother than dark areas.

Conclusion

We have described a Bayesian motivation for the NL-means filter and justify some intuitive adaptations described in previous papers [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF][START_REF] Mahmoudi | Fast image and video denoising via nonlocal means of similar neighborhoods[END_REF][START_REF] Coupé | Fast non-local means denoising for 3D MR images[END_REF]. The proposed framework yields to a filter which is better parametrized: the size of the adaptive dictionary and the noise variance are computed from the input image, and the patch size must be large enough.

Our Bayesian approach has been used to remove image-dependent noise and could be adapted in applications with appropriate noise distributions. A more thoroughly evaluation with other methods [START_REF] Awate | Higher-order image statistics for unsupervised, information-theoretic, adaptive, image filtering[END_REF][START_REF] Azzabou | Random walks, constrained multiple hypothesis and image enhancement[END_REF][START_REF] Tschumperlé | Curvature-preserving regularization of multi-valued images using PDE's[END_REF][START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF], and recent improvements of the NL-means filter described in [START_REF] Buades | The staircasing effect in neighborhood filters and its solution[END_REF], would be also desirable.

  max(var(z(xj )),var(z(xi))) min(var(z(xj )),var(z(xi)))

Fig. 1 .

 1 Fig. 1. Comparisons of different NL-means filters. From top to bottom and from left to right: part of noisy images (Barbara, Lena, σ = 20), NL-means filter, vectorial NL-means filter, Adaptive NL-means filter, plugin Adaptive NL-means filter; numerical results for each filter; comparison of exponential weights for the original NL-means filter (blue dashed line) and for the Adaptive NL-means filter (solid green line) (see text).

Fig. 3 .

 3 Fig. 3. Denoising of real noisy images. From top to bottom: original images, denoised images, "method noise" experiments. The old photograph (left) has been denoised using an usual additive noise and the four other images (right) are denoised using the image-dependent noise model (γ = 0.5, see text).

  , we modify the noisy image (σ = 20) NL h z NL h z ANLσ,nz ANLσ,ne uANL

		Timings	Lena	Barbara
	NL h z	58.1 sec	512 × 512 512 × 512 31.85 db 30.27 db
	NL h z	96.2 sec	32.04 db	30.49 db
	ANLσ,nz	75.2 sec	32.51 db	30.79 db
	ANLσ,n e uANL	173.3 sec	32.63 db	30.88 db
	Fast ANLσ,nz	10.6 sec	32.36 db	30.61 db
	Fast ANLσ,ne uANL 21.2 sec	32.49 db	30.71 db

The F -distribution is used to compare the variance of two independent samples from a normally distributed population.

Acknowledgments. We thank B. Delyon and P. Pérez for fruitful conversations and comments.

 3) and the maximum weighting w(x, y) at x = y is given by w(x, x) = max x =y w(x, y)) ; right: performance of the plugin Adaptive NL-means filter ANLσ,ne uANL for different signal-to-noise ratios (WGN).


José V Manjón

Pierrick Coupe

Antonio Buades

Louis D Collins

Montserrat Robles

Pierrick Coupé

D Louis Collins

MRI denoising based on Sparseness and Self-Similarity

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

In this paper, we present two new approaches for MRI denoising. The first is an extension of the original method proposed by [START_REF] Guleryuz | Weighted averaging for denoising with overcomplete dictionaries[END_REF]. Based on local 3D DCT hard thresholding, our proposed method has been adapted to deal with Rician noise (typical of magnitude MR images) using a pseudo-oracle principle. The second proposed method is a new rotationally invariant 3D version of the Rician-adapted Non Local Means filter [START_REF] Buades | A non-local algorithm for image denoising[END_REF], Coupé 2008a[START_REF] Wiest-Daesslé | Rician noise removal by non-Local Means filtering for low signal-to-noise ratio MRI: applications to DT-MRI[END_REF] that uses a prefiltered image.

Methods

Typically, the observed noisy MR image y is considered to be a linear combination of a noise-free image x and a white noise realization n from the measurement process.

Therefore, the goal of any denoising algorithm is to find a good estimate , given y. This section describes two new methods for image denoising based on two different image properties and compares them to two state-of-the-art denoising methods.

Denoising using sparseness

In [START_REF] Guleryuz | Weighted averaging for denoising with overcomplete dictionaries[END_REF] method, the estimate of the noise-free image is obtained using a hard thresholding method. An overcomplete set of 3D block DCTs (4×4×4 block size) is used (Guleryuz used 8×8 blocks in his 2D method, which yields the same number of coefficients). The local denoised estimate at block j, , is obtained by applying a hard thresholding rule. Finally, all local estimates are combined from all overlapping j blocks at position i. On the other hand, Oracle-based filters assume that, if the null coefficients of the original noise-free image are known, they can be used to improve the denoising. Unfortunately, in practice no noise-free image x is available; hence, this approach is unfeasible. However, if we relax the Oracle condition, we can use a prefiltered image with, for example, the described DCT3D method as an approximation of the noise-free image. With this approach, each block of coefficients of the noisy image y can be better thresholded using the corresponding block coefficients of the prefiltered image. The final image reconstruction is performed using the same approach described above. We will refer to this method as Oracle DCT3D (ODCT3D).

The only parameter of this method is the threshold τ from the prefiltering step. In all of our experiments, we used a τ value of 2.7 σ (where σ is the standard deviation of the noise), which is common in DCT and wavelet thresholding methods [START_REF] Mallat | A wavelet tour of signal processing[END_REF].

Denoising using self-similarity

Originally proposed by [START_REF] Buades | A non-local algorithm for image denoising[END_REF], the NLM filter takes advantage of the high level of pattern redundancy in an image, achieving high-quality image denoising by averaging similar realizations of the noisy signals. Basically, this filter reduces the noise in an image by averaging voxels that originally had the same intensity in the noise-free image. To this end, [START_REF] Buades | A non-local algorithm for image denoising[END_REF] suggested that voxels with similar neighborhoods (small 3D patches in our volumetric case) tend to have similar original values. However, if the images are prefiltered using the described ODCT3D method there is no need to compare 3D patches to measure voxel similarities. In this case, the voxel value and the local mean can be used to efficiently estimate the voxel similarity.

Furthermore, the similarity so defined is rotationally invariant and thus the number similar voxels around each voxel will be increased if compared to the original method.

We will refer to this method as PRI-NLM3D.

Results

The proposed methods were compared with two state of the art MRI denoising methods: WSM from Coupe (2008b) and ONRAD from Krisian (2009).. Example results using simulated data from the well-known Brainweb phantom can be observed in Fig.

1

. An example of the application of the PRI-NLM3D method over real data can be observed in Fig. 2.

Conclusion

To new methods for MRI denoising have been presented. The new methods are able to better reduce the random noise in the images when compared with other state-ofthe-art methods. Both methods run in less than a minute, making them usable in most research and clinical settings.

Fig. 1 Fig. 2 -

 12 Fig.1 Example filtering results for an axial slice of the T1w BrainWeb phantom (Rician noise level of 15%). The third row shows the absolute value of the image residuals for the different methods.


Frédéric Blanc

François Peters

Elisabeth Lemaire

Experimental Signature of the Pair-Trajectories of Rough Spheres in the Shear-Induced Microstructure in Non-Colloidal Suspensions

The shear-induced microstructure in a semi-dilute non-colloidal suspension is studied. A high resolution Pair Distribution Function in the plane of shear is experimentally determined. It is shown to be anisotropic, with a depleted direction close to the velocity axis in the recession quadrant. The influence of roughness on the interaction between particles is quantitatively evidenced. The experimental results compare well with a model from particle pair trajectories.

The influence of the particle interactions on the rheology of non-Brownian suspensions, such as normal stress differences and shear-induced migration, and their link with the shear-induced microstructure, are questions of central interest and remain an experimental challenge. A quantitative measurement of the microstructure is provided by the Pair Distribution Function (PDF). Since the pioneering work of Bossis and Brady [1], in Stokesian Dynamics, the PDF of an unbounded suspension of interacting spheres undergoing a simple shear flow is known to have an angular dependence. In a semi-dilute suspension, experimental [2], numerical [3,4] and analytical [5] studies have shown the existence of a particle depleted area along the velocity direction. In their numerical work, Drazer at al. [4] showed that, besides this anisotropy, a fore-aft asymmetry appears in the presence of strong enough non-hydrodynamic forces acting between the spheres. Repulsion forces between particles eliminate closed orbit trajectories and the particles that were close from each other in the approach quadrant are redistributed in the recession quadrant, leading downstream to a tilt of the depleted area away from the velocity direction. Though only a few experimental measurements of PDF have been performed in non-Brownian suspensions [6,7], the depleted area close to the velocity direction in a simple shear flow has already been accurately studied by Rampall et al. in semi-dilute suspensions (φ=0.10-0.15). They explain their experimental results by a model that computes the PDF from the trajectories of two colliding rough spheres. Besides the depleted area, their theoretical PDF shows fore-aft asymmetry. These predictions are in qualitative agreement with their experimental results. However these early measurements suffer from a quite low resolution that precludes a fully quantitative comparison.

In this letter, we report the measurement of a high resolution PDF in a semi-dilute sheared suspension. The results are in qualitative agreement with the works mentionned above but demonstrate a new property: in the recession quadrant, right above the depleted zone, a taillike high particle concentration zone is apparent, that is not present in the approach quadrant. A trajectorybased model shows that this tail is a direct signature of the pair-trajectories of rough spheres in the shearinduced microstructure. This has been observed with a suspension of PMMA particles (Arkema BS510) dispersed in a Newtonian liquid (η 0 = 0.85 Pa.s) designed by Cargille Laboratories (Immersion liquid Code 11295101160) in order to match both the density (d=1180kg.m -3 ) and the refractive index (n=1.49) of the particles. To improve the index matching, the temperature is controlled (T = 37 • C). The particles are sieved to reduce the initially large size distribution. The final size distribution is determined by direct microscope observation over several hundreds of particles. The particles diameter is 2a=170 ±12µm. The suspension is sheared in a cylindrical Couette cell, R i =19mm in inner radius, R o =24 mm in outer radius and 60 mm in height. The whole apparatus design is precisely discribed in [8,9]. Some fluorescent dye (Nile Blue A) is dissolved in the host liquid and a thin laser sheet (635 nm, 30 µm in thickness), perpendicular to the cell axis, illuminates the suspension. The position of the laser sheet, 11 mm above the bottom of the cell, allows to neglect the bottom end effect on the flow [8]. A CCD camera positioned under the cell and equipped with a long pass filter (650nm) images the interstitial fluid, showing the particles in negative (Fig. 1).

The particle positions and radii are detected using an open source Matlab procedure based on the Circular Hough Transform [10,11]. The position of the particle centers is determined with a subpixel resolution. For the sake of comparison, the mean particle radius a=85µm corresponds to 14.7 pixels on the frame. The inner cylinder angular velocity is controlled by a rheometer (Mars II, Thermofisher) and set to the value Ω =1 rpm. The Reynolds number Re = (a

2 d/η 0 )Ω R i /(R o -R i )≈ 4.10 -6
and the Péclet number P e = (η 0 a 3 /k B T )ΩR i /(R o -R i )≈ 6.10 7 allow to neglect inertia and Brownian motion. The shear rate γ = 2ΩR

2 i R 2 o /((R 2 o -R 2 i )R 2
) depends on the distance to the cylinders center R, and its variation across the gap is approximately 50%. This is not a problem however, since at low concentration, the non-brownian suspension rheology does not depend on the shear rate.

For each frame, a 2D PDF is built in the following way. First, for each particle pair, the relative position vector of the pair in the local flow direction frame of particle 1 is determined (see Fig. 2(a)). The plane of the particles relative position r 12 (ρ,θ) is evenly sampled (∆ρ=1 pixel, ∆θ =2 π/80 rad), and the number of particle pairs (1,2) corresponding to a particular sample N (ρ, θ) is normalized by the total number of origin particles N 1 , by the sample surface, and by the mean surface density of particles n 0S .

g(ρ, θ) = 1 n 0S N (ρ, θ) N 1 ρ∆ρ∆θ (1)
The PDF is then averaged over all images. g(ρ, θ) is thus a surface distribution function, that is related to the real 3D PDF through the integration over an interval [-2∆z max ,2∆z max ] in the z direction, as explained in [2] and shown in Fig. 2(b). To reduce this interval, an histogram of the apparent radii of the detected particles is built, and only particles greater than or equal to 13 pixels (0.88a) are taken into account. From this, we evaluate ∆z max to approximately 0.5a (Fig. 2). To avoid frame boundary effects on the large ρ limit of g(ρ, θ), all origin particles (1) outside a central area are rejected (Fig. 1). The distance of the box edges to the frame edges is larger than 11a. In addition, in order to reduce the hydrodynamic wall effects, no particle (2) closer than an arbitrary distance of 3a to the inner or outer cylinder is considered. Thus, the PDF should not suffer any frame boundary induced spurious decrease for interparticle distances smaller than 8a.

We turn now to the PDF measured for a φ = 5% suspension. 50000 frames are recorded at the rate 0.33 fps during more than 41 hours. We have checked that the PDF computed using the frames 1 to 10000 was not different from the PDF shown here, except for statistical noise, suggesting that any transient state can be neglected. Actually, we evaluate the transient time to the time necessary to break the closed orbits, i.e. the collision time 1/ γ ≈ 3s that is much smaller than the experiment duration. Fig. 3 displays the measured PDF. As already shown by Rampall et al. [2], the PDF is anisotropic, with a depleted area in the approximate di- rection of the flow velocity. To be specific, as shown in Fig. 4, the depleted direction is situated around the angle θ ≈ 6 • in the recession quadrant. This slight tilt angle away from the velocity axis suggests some particle roughness. Away from this depleted direction, at large distance from the origin, the PDF is quite isotropic, with a value close to 1, as expected from the normalization in Eq.1. We note in passing that any strong influence of the walls or of a possible concentration gradient, that would destroy this symmetry, can be ruled out. Another classical property is the large increase of the PDF near the surface ρ=2a. Moreover, thanks to the large number of frames that leads to a high-resolution PDF, we observe for the first time a tail-like high particle concentration zone in the recession quadrant, that is not present in the approach quadrant.

This high concentration tail clearly recalls the relative trajectory of two rough particles in a simple shear flow (Fig. 5). As we will show in the following, this tail allows a good estimation of the particles roughness. Indeed, such roughness-modified trajectories have already been used to compute the shear-induced diffusivity in a dilute suspension [12]. In their paper, a roughness length a is defined. When the distance between the particles equals (2+ )a, the relative radial velocity is cancelled, preventing the particles to get closer from each other. This results in a break of the fore-aft symmetry, the spheres separating on streamlines further apart than on their approach (Fig. 5). Following the idea of Rampall et al. [2], we use the roughness-modified trajectories to compute the PDF. For an introductory purpose, we recall here the classical isotropic expression obtained by Batchelor and Green for the PDF from the motion of two smooth particles in an ambient linear flow [13,14].

g BG (r) = 1 1 -A(r) exp ∞ r 3 s B(s) -A(s) 1 -A(s) ds (2)
where A and B are mobility functions that depend only on the interparticle distance r, and whose expressions can be found in [12]. Eq.( 2) applies everywhere a trajectory comes from infinity. In a simple shear flow however, there exists an infinite region of closed trajectories in which it is impossible to predict the PDF values. In the plane of shear, the particles that undergo such closed orbits approach each other at a distance smaller than approximately 2.00004 a [12,15]. Following Rampall et al. [2], and in agreement with the numerical work of Drazer et al. [3], we suppose that, owing to the roughness that prevents particles to come too close to each other, the bound pairs are eliminated. Thus, only the trajectories from infinity will be taken into account.

We use the expression of the mobility functions A(r) and B(r) given in [12] to compute the relative velocity V of two particles in a simple shear flow V ∞ = γy e x , together with their relative trajectory. A(r) is set to 1 when contact occurs at r=(2+ )a. The integration scheme, namely a variable order Adams-Bashforth-Moulton predictor-corrector solver, was tested as in [12] for =0 and provided the same results. Fig. 5 displays trajectories in the plane of shear for different values of the roughness length. For a given value of , all trajectories that drive the particles in contact concentrate downstream on a single trajectory. Thus, downstream, no particle can locate between the x-axis and this trajectory. We note that the displacement of the trajectories depends strongly on . To determine the PDF, we com- pute the trajectories from evenly spaced origins in the plane (x 0 =-10a), with a sample time ∆t=0.01/ γ. The set of trajectories forms a distribution of particle positions r in the probability space, and each point has a velocity V . Since the number of particles is conserved, this distribution is a discrete form of a function p( r) that obeys Eq.3:

∂p ∂t + ∇. p V = ∇. p V = 0 (3)
We recall [3,14] here that the PDF g( r) is related to the conditional probability P( r| r 0 =0) of finding a sphere with its center at r given that there is a sphere with its center at r 0 =0 and the bulk density of particles n 0 through the equation g( r)=P( r| r 0 =0)/n 0 . Since P( r| r 0 =0) obeys Eq.3 too [14], it is sufficient to impose the limit p→1 as r→∞ to set p( r)=g( r). To this purpose, each trajectory is weighted by the inverse of the density in the plane x=x 0 , namely L y L z V x (x 0 ) ∆t, where L y and L z are the distances between the trajectories origins in that plane along the directions y and z. From this set of points, we compute the relevant 2D PDF that could be compared to our experimental results. The region of interest is defined by

z ∈ [-2∆z max , 2∆z max ], ρ = x 2 + y 2 ∈ [2a, R max ]
where we chose R max =4a. This region is sampled in polar coordinates with the sampling periods ∆ρ=a/15, ∆θ=π/40, and the number of points in each cell divided by its volume 4∆z max ∆(ρ 2 )/2 ∆θ gives the sought PDF.

We restrict the plane (x=x 0 =-10a) of trajectory origins outside the closed orbits region. In order the PDF to be smooth enough, the value of the sampling period in the plane x = x 0 have been chosen: L z /(2∆z max )=1/20, L y /a=1/50. ∆z max =0.5a has been manually to fit the experimental data at best. We note that this value corresponds also to the estimation of ∆z max from the distribution of particle radii (Fig. 2(b)). We have checked that this method yield the isotropic distribution of Eq.2 in the case of a pure strain flow. Concerning the simple shear flow, we obtain the PDF in Fig. 6.a for a roughness a=3.10 -3 a=255nm. A trajectory in the plane of shear for which contact occurs is also shown for the same roughness. It is centered on the tail in the recession quadrant. We have checked that away from the depleted area and the tail, the PDF is identical to the isotropic g BG of Eq.2, in agreement with [3]. Fig. 7.a displays the experimental and computed PDF g(ρ, θ) as a function of θ for two values of the distance ρ. The agreement is very good, especially concerning the position of the tail that appears clearly on the experimental curve. The radial profiles can be compared on Fig. 7.b for three different angles. The model describes accurately the experimental PDF, both in the direction of the velocity gradient (85.5 • < θ <90 • ) and near the depleted direction (18 • < θ <22.5 • ). However, in the depleted direction (0 • < θ <4.5 • ), the model predicts the total lack of particles in contrast to the experimental results. As already mentioned in [2], this is probably the consequence of multibody effects. Finally, we have performed surface roughness measurements on the particles using an Atomic Force Microscope. The main roughness height is small (≈20nm), with isolated bumps, approximately 200 nm in height, as shown on Fig. 6.b. Those bumps are frequent, even if we did not perform any extensive statistical study. These direct measurements are fully consistent with the experimental PDF.

In conclusion, we have measured a high resolution PDF in a sheared dilute suspension. Its fore-aft asymmetry is quantitatively explained by the particle roughness as modeled by a minimum approach distance between particles a [12], which can be determined with reasonable accuracy. Some experiments are currently performed on more concentrated suspensions (φ=0.10-0.55). In the semi-dilute regime, we observe that the tilt angle from the velocity axis hardly depends on the volume fraction. For volume fractions higher than 0.2, the tilt angle increases with the particle concentration, as already reported by Drazer et al. [3]. In addition, the tail turns blurred, partly due to poorer statistics (1000 frames per PDF), but also to increasingly large multibody effects, that modify the pair trajectories.

We are grateful to L. Lobry for fruitful discussions, and A. Zenerino for the AFM images. This work was supported by the ANR (Program No. ANR-08-BLAN-0048-CSD 2).

FIG. 1 :

 1 FIG. 1: Typical image. φ = 5%. Only the origin particles in the white box are considered.

FIG. 2 :FIG. 3 :

 23 FIG. 2: (a) Pair position vector relative to the local azimuthal flow direction (b) aapp is the apparent radius of the particle.

FIG. 4 :

 4 FIG. 4: PDF averaged over the ρ-interval [1.97a 2.04a]. The approach quadrant is θ ∈ [90 • , 180 • ].

3 ! = 3 . 10 - 3 ! = 10 - 2 FIG. 5 :

 331031025 FIG. 5: Trajectory from the point x/a=-10, y/a=0.5, z/a=0 for different values of the roughness a

FIG. 6 :

 6 FIG. 6: (a) PDF from particle trajectories. The roughness is =3.10 -3 . Red line: a trajectory in the plane of shear for which contact occurs. (b) AFM image of the particle surface.

FIG. 7 :

 7 FIG. 7: Experimental (symbols) and computed (lines) PDF. (a) g(ρ,θ) vs θ (•) 1.97 < ρ < 2.04; ( * ) 2.31 < ρ < 2.38; (-) 1.93 < ρ < 2.00; (--) 2.33 < ρ < 2.40; (b) g(ρ,θ) vs ρ/a (•) 0 < θ < 4.5; • ( * ) 18 • < θ < 22.5 • ; ( ) 85.5 • < θ < 90; • ; (-) 0 < θ < 4.5 • ; (--) 18 • < θ < 22.5 • ; (-•) 85.5 • < θ < 90 •