Frédéric Aman
email: frederic.aman@imag.fr

Véronique Aubergé
email: veronique.auberge@imag.fr

Michel Vacher
email: michel.vacher@imag.fr

Influence of expressive speech on ASR performances: application to elderly assistance in smart home

Keywords: Expressive speech, distress call, Ambient Assisted Living, Home Automation

Smart homes are discussed as a win-win solution for maintaining the Elderly at home as a better alternative to care homes for dependent elderly people. Such Smart homes are characterized by rich domestic commands devoted to elderly safety and comfort. The vocal command has been identified as an efficient, well accepted, interaction way, it can be directly addressed to the "habitat", or through a robotic interface. In daily use, the challenges of vocal commands recognition are the noisy environment but moreover the reformulation and the expressive change of the strictly authorized commands. This paper focuses (1) to show, on the base of elicited corpus, that expressive speech, in particular distress speech, strongly affects generic state of the art ASR systems (20 to 30%) (2) how interesting improvement thanks to ASR adaptation can regulate (15%) this degradation. We conclude on the necessary adaptation of ASR system to expressive speech when they are designed for person's assistance.

Introduction

Smart homes aim at anticipating and responding to the special needs of elderly people living alone at their own home by assisting them in their daily life. A large variety of sensors are used in this framework [START_REF] Peetoom | Literature review on monitoring technologies and their outcomes in independently living elderly people[END_REF] but one of the best interfaces seems to be the speech interface, that makes possible interaction using natural language and that is well adapted to people with reduced mobility. The purpose of this study is to develop a system able to detect distress calls uttered by elderly people. In a previous study, we reported that elderly voices resulted in degraded performances for generic Automatic Speech Recognition (ASR) systems [START_REF] Aman | Analysing the Performance of Automatic Speech Recognition for Ageing Voice: Does it Correlate with Dependency Level?[END_REF]. The content of test sentences was related to distress but these sentences were read and uttered in a neutral manner without any particular emotion; in this article we are using "neutral" when utterances are not highly perturbated by intentions, attitudes or emotions.

In a real situation and in the framework of home automation, vocal commands can be expressive for 2 reasons:

1. if the user integrates a representation of the home as a communicative entity [START_REF] Clarcke | A Space Odyssey[END_REF], thus prosody and morphosyntax can implement intentions, attitudes and social affects; 2. the command is motivated by an emotional context (panic, happiness to speak with a relative, etc.).

In a distress situation, a cry for help (i.e., "I fell", "Help me") may be charged in emotion, prosodic modifications influencing voice quality [START_REF] Audibert | Prosodie de la parole expressive: dimensionnalité d'énoncés méthodologiquement contrôlés authentiques et actés[END_REF][5] so that emotions could be perceived by the interlocutor. Many studies are related to emotion recognition from voice (distress, fear, stress, etc.) [START_REF] Vidrascu | Analyse et détection des émotions verbales dans les interactions orales[END_REF][7][8] [START_REF] Chastagnol | Reconnaissance automatique des dimensions affectives dans l'interaction orale homme-machine pour des personnes dépendantes[END_REF], in [START_REF] Soury | Stress detection from audio on multiple window analysis size in a public speaking task[END_REF], the authors used energy, pitch, jitter and shimmer. Scherer et al. [START_REF] Scherer | Vocal communication of emotion: A review of research paradigms[END_REF] compare effects of different emotions on prosodic parameters with regard to neutral voice. An immediate and more reachable issue consists in making ASR more robust in presence of emotion instead of recognizing the distress state but few studies were done in this domain excepted the study of Vlasenko et al. [START_REF] Vlasenko | Towards robust spontaneous speech recognition with emotional speech adapted acoustic models[END_REF]. They observed a performance degradation with expressive speech when ASR were trained with neutral voice.

In former studies [START_REF] Aman | How affects can perturbe the automatic speech recognition of domotic interactions[END_REF], the different ASR performances between acted expressive voice and spontaneous expressive voice from a corpus recorded using a Wizard of Oz technique was brought to light; in this study, results were compared to those of neutral voice. The aim of our study is to highlight the lack of robustness of ASR systems with expressive voice. For in-depth study, we recorded the Voix Détresse corpus which is made of distress sentences uttered in expressive and neutral manner in an elicitation protocol. A reduction of performance decrease was observed thanks to the use of an adapted acoustic model. Expressive voice characterisation is presented at the end.

Method

Ideally, this work should have been held on natural dataset. But in the spontaneous corpus of distress, situations (calls to emergency centers, etc.) are: (1) not addressed to a system (human isolation), but to a person devoted to assistance (voice human contact);

(2) after a phone call, that is not in the so serious situation that the elderly is not anymore able to give a call. It is why, in this primary study, we decided to be closer to the natural situations but to control the relevant cues that will be the real situations for which there is yet no available data. Thus we proposed a devoted elicitation protocol. The Voix Détresse French corpus was recorded in our laboratory in the following elicitation protocol:

a list of 20 prototypical utterances was selected from the AD80 corpus [START_REF] Aman | Analysing the Performance of Automatic Speech Recognition for Ageing Voice: Does it Correlate with Dependency Level?[END_REF], that was validated to be representative of distress commands spontaneously produced by subjects in real distress situations at home; to elicitate expressive speech in the specific context of these distress situations, for each utterance is associated a devoted photography showing an elderly in this given distress home situation; the subject had to produce the utterance (as many times as the subject decided) trying to stand in that individual's shoes; before to enter the elicitation procedure, in order to get a control sub-corpus, each speaker had to read the 20 distress sentences in a neutral manner, that is without any elicitation support. The expected emotions for these distress elicitations were mainly negative emotions like fear, anger, sadness. The elderly speakers had some difficulties to enter the elicitation protocol, so that only 5 speakers, all females, had sufficient performances of "quite-natural" distress expressions to be selected for the corpus. It has to be noted that the females are the main users of the Smart homes, because the elderly population over 80 is mainly females, and moreover some French investigations show that the females are more interested by such Smart homes solutions. To complete these low number of elderly subjects, together with collecting a non elderly control sub-set, 20 other subjects, from 23 to 60 years old, were recorded in the same elicitation protocol. Voix Détresse is described in Table 1. The corpus is made of 1481 sentences, its length is 24min 19s.

One half of the corpus was reserved for model training and the remaining part for testing. Since 25 speakers were recorded and regarding each speaker, approximately 10 neutral sentences were reserved for adapting and 10 for testing, and about 20 expressive sentences for training and 20 for testing (see Table 2).

For our study, ASR was operated using Sphinx3 [START_REF] Seymore | The 1997 CMU Sphinx-3 English broadcast news transcription system[END_REF] based on a context dependant triphone acoustic model. Our HMM model was trained with BREF120 corpus [START_REF] Lamel | BREF, a large vocabulary spoken corpus for french[END_REF] which is made of 100 hours of texts read by 120 speakers (65 females, 55 males), this generic model was called BREF120.

For the decoding, a trigram language model (LM) with a 10K lexicon was used. It results from the interpolation of a generic LM (with a 10% weight) and a specialized LM (with a 90% weight). The generic LM was estimated on about 1000M of words from Gigaword, a collection of French newspapers3 , this model is unigram and made of 11018 words. The specialized LM was estimated from the distress and call for help sentences of the AD80 corpus, it is made of 88 unigrams, 193 bigrams et 223 trigrams. This combination has been shown as leading to the best WER for domain specific application [START_REF] Lecouteux | Distant speech recognition in a smart home: Comparison of several multisource ASRs in realistic conditions[END_REF]. The interest of such combination is to bias the recognition towards the domain LM but when the speaker deviates from the domain, the general LM makes it possible to correctly recognize the utterances. Table 3. WER (%) from Sphinx3 decoding (BREF120 acoustic model) for the two groups.

E01F84 E02F85 E03F83 E04F67 E05F73 Y01M31 Y02M30 Y03F60 Y04F26 Y05M26 Y06M32 Y07F23 Y08F25 Y09F29 Y10F24 Y11M25 Y12F53 Y13M41 Y14M59 Y15M28 Y16F28 Y17F25 Y18F30 Y19F24 Y20F29
3 Experiment

Decoding with BREF120 acoustic model

Decoding was operated on the Test part of the Voix Détresse corpus. WER results are presented in Figure 1 for each speaker and as the average for each group in Table 3. A severe performance degradation can be observed: absolute WER difference is 29.95% for younger speakers and 19.6% for elderly speakers. A Welch's t test was operated in order to determine if WER differences between neutral and expressive voices are significant. Test results indicate a significant difference between the corresponding WER (t=-7.2026. df=21.365. p<0.05 (p=3.834e-07) for younger speakers and t=-2.5165. df=7.874. p<0.05 (p=0.03645) for elderly speakers).

As shown in Figure 1, we can observe a WER increase for expressive voice with regard to neutral voice for all speakers, excepted for Y08F25, and these variations are more or less important according to the speaker. One explanation of these differences could be due to the difficulty of the participants, who are not professionnal, to "play the game" of distress during the records.

Overall, persons more comfortable to play distress situations are persons with a higher WER for distress sentences. In real situations, it is feared that very expressive sentences will be badly recognized by ASR, while it is essential to act in this situation.

Decoding with an adapted acoustic model

In order to establish the possibility of improvement, a MLLR adaptation of BREF120 model to each speaker was operated: (1) BREF120_N by using the [START_REF] Soury | Stress detection from audio on multiple window analysis size in a public speaking task[END_REF] sentences uttered by the speaker, (2) BREF120_E by using the 20 expressive sentences, and (3) BREF120_N+E by using neutral and expressive sentences. For each speaker, 3 different models adapted to his voice were obtained. Decoding was operated as in section 3.1 on the 739 neutral and expressive sentences ot the Test part of the corpus using the 3 adapted acoustic models. Results are given in Table 4 with regard to the standard BREF120 model.

ANOVA was operated on the 2 groups neutral and expressive voices (samples are Gaussian and homogeneity of group variances is verified). There is no significant difference inside the samples of the neutral voice group (F(3;96)=1.293; p=0.281) but a significant difference (p-value < 0.05) inside the samples of the expressive voice group (F(3;96)=7.828; p=9.96e-05). Results of the Tukey HSD test are given in Table 5 for neutral voices and in Table 6 for expressive voices. For neutral voices, no significant difference is observed between the WER obtained from the different acoustic models, p>0.05 for each acoustic model pair. By contrast for expressive voices, a significant WER difference exists between the generic BREF120 model and the models adapted to expressive voice BREF120_E, with p<0.05 (p=0.0011), and an absolute difference between average WER equals to 14.78%. Conclusion is the same for BREF120 and BREF120_N+E, with p<0.05 (p=0.0002), and an absolute difference of 16.78%. By contrast, there is no significant difference between the other models.

Therefore, with regard to expressive sentences, speaker adaptation from neutral sentences is not sufficient to improve significantly the WER. Models adapted to expressive voice may be used, but adaptation can be done from mixing neutral and expressive sentences because there is no significant difference between models BREF120_E and BREF120_N+E. In the same way, neutral sentences can be decoded using adapted model BREF120_N+E without important performance decrease.

Expressive voice characterisation

Average values of prosodic parameters, flow, F0, jitter, shimmer et Harmonic to Noise Ratio (HNR), were measured and compared between neutral and elicited expressive voices. Data were the 739 neutral and expressive sentences of the Testing part of Voix Détresse.

We observed for expressive voice, with regard to neutral voice and in average, a decrease of the flow, an increase of F0 and HNR, and an increase of jitter and shimmer. Averages concerning all speakers for each parameter are given in Table 7 with significance tests of the difference between neutral and expressive voices.

In a first step, we considered for our analysis that the distress corpus (i.e., the expressive sentences) of Voix Détresse was homogeneous with regard to expressiveness of acted emotion (eventually with variable expressive intensity), and varied possibly only with the speaker (difference of interpretation or idiosyncratic physiologic variation). However, by more precise hearing of the considered utterances, sentence by sentence, we can suppose that the nature of the distress varies according the situation suggested by the picture submitted for elicitation, and that this nature could be homogeneous by situation (i.e., by sentence) without dominant variability of the subject. For example, the sentence "Je ne me sens pas bien" (I don't fell very well) was associated to a picture showing a very demoralized man taking his arm, supported by an other person, this person was very close and had already started to help the man in distress situation. By contrast, the sentence "A moi" (Help me) is associated to a picture showing an isolated drowning person: it is a question of a request for emergency assistance by an isolated person in a critical condition. The two situations are two opposite situations of our corpus: in the first case a helper took care of the person but in the second case a person in a critical situation was calling but there was no identified helper.

That is what we wanted to verify by evaluating prosodic parameters not by speaker but by context, The 489 expressive sentences were spread into 10 contexts. Jitter, shim- According to the 2D model of Russel [START_REF] Russell | A circumplex model of affect[END_REF] of cognitive psychology, emotions can be represented along 2 axis: valence (positive/negative) and arousal (active/passive). Thus, the context "A moi" could correspond to a very "active" distress, by hearing they are perceived as very stressed. The context "Je ne me sens pas bien" could be more "passive" and is perceived as the less stressed. Measures corroborate these observations: we can observe in Table 8 that the context "A moi" has the lowest jitter and the highest F0 with regard to the 9 other contexts, and the context "Je ne me sens pas bien" has the highest jitter and the lowest F0.

We can see that some analog expressions are characterized by the same parameter values, as "Du secours s'il vous plaît" and "Ne me laissez pas tout seul" with jitter=1.550% and F0=276Hz for these 2 contexts.

Conclusion

After adaptation, WERs of expressive voices remain higher than WERs of neutral voices with the generic model BREF120, the difference is 13.1%. Moreover, the record of spontaneous expressive voices is very difficult and the Voix Détresse corpus is made of acted expressive voice. In all likelihood, performances may be lower in a real case.

Acoustic characteristics of distress are very subtle. It is clear that in the harsh acoustic conditions of a smart home, we have to be extremely vigilant on the physiologic, functional and communicative context of utterances to recognize: a high amount of data must be recorded in order to obtain a robust system. Moreover these data must be finely characterized as regard to their expressiveness. Those conditions are necessary but insufficient ones for the use of speech technologies in real conditions for the assistance of elderly persons.

Fig. 1 .

 1 Fig. 1. WER from Sphinx3 decoding (BREF120 acoustic model) for neutral and expressive voices

Table 1 .

 1 Voix Détresse corpus.

	Composition	Younger group		Elderly group
	All speakers	12F, 8M, 23-60 years old		5F, 67-85 years old
	Neutral voice	400 sentences, 5min 41s (total)	100 sentences, 1min 38s (total)
	Expressive voice	782 sentences, 12min 44s (total)	199 sentences, 4min 15s (total)
	Young voice	Training	Test	Elderly voice	Training	Test
	Neutral	200	200	Neutral		50	50
	Expressive	393	389	Expressive	99	100

Table 2 .

 2 Number of sentences of the Training and Testing partitions of the corpus.

Table 4 .

 4 neutral Training WER (%) for the different voices using generic and adapted models.

	Group	BREF120	BREF120_N	BREF120_E	BREF120_N+E
	Young neutral	9.27	7.80	11.03	7.21
	Young expressive	39.22	28.86	22.45	20.23
	Elderly neutral	18.82	17.06	16.48	15.88
	Elderly expressive	38.42	35.48	31.64	30.50
	Average Neutral	11.18	9.65	12.12	8.94
	Average Expressive	39.06	30.18	24.28	22.28
		BREF120 BREF120_N BREF120_E BREF120_N+E
	BREF120	-	p=0.8298	p=0.9529	p=0.5977
	BREF120_N	-	-	p=0.5174	p=0.9786
	BREF120_E	-	-	-	p=0.2923
	BREF120_N+E -	-	-	-
	WER (%)	11.18	9.65	12.12	8.94

Table 5 .

 5 Neutral voice group: p-value of Tukey HSD test and WER.

Table 6 .

 6 Expressive voice group: p-value of Tukey HSD test and WER.

			BREF120	BREF120_N	BREF120_E	BREF120_N+E
	BREF120		-		p=0.0976	p=0.0011	p=0.0002
	BREF120_N	-		-	p=0.4120	p=0.1683
	BREF120_E	-		-	-	p=0.9526
	BREF120_N+E	-		-	-	-
	WER (%)		39.06		30.18	24.28	22.28
	Parameters	Neutral voice	Expressive voice	Welch's t test (significant difference if p<0.05)
	Flow (Phonemes/s)		13.58		11.71	p<0.05 (p=0.0019)
	F0 (Hz)		162.74		260.43	p<0.05 (p=1.16e-06)
	Jitter (%)		3.07		2.30	p<0.05 (p=0.00015)
	Shimmer (%)		13.63		9.07	p<0.05 (p=1.41e-07)
	HNR (dB)		12.44		14.76	p<0.05 (p=0.00077)

Table 7 .

 7 Average prosodic parameters as a function of kind of voice, and Welch's t test p-values.

Table 8 .

 8 Prosodic parameters for some different French distress sentences.mer, F0 and flow were measured, results are given in Table8. It appears that 2 extreme situations "A moi" (Help me) and "Je ne me sens pas bien" (I don't feel very well) are characterized, throughout all speakers, by extreme values of F0, jitter and shimmer; this fact could confirm our hypothesis on the variability of the nature of the distress.

	Distress sentences	Jitter (%)	Shimmer (%)	F0 (Hz)	Intensity (dB)
	A moi !	1.028	5.962	317	70.27
	Oh la la !	1.367	6.383	224	65.72
	J'ai un malaise !	1.411	6.778	232	63.53
	Aidez-moi !	1.423	6.801	279	66.46
	Du secours s'il vous plaît !	1.550	6.791	276	64.80
	Me laissez pas tout seul !	1.550	7.131	276	64.27
	Qu'est-ce qu'il m'arrive !	1.654	7.470	239	62.21
	Au secours !	1.679	5.414	294	69.79
	Je ne peux plus bouger !	1.708	8.181	282	64.72
	Je ne me sens pas bien !	1.718	8.372	220	58.08

http://catalog.ldc.upenn.edu/LDC2006T17

Acknowledgements This study was supported by the French funding agencies ANR and CNSA through the project CIRDO -Industrial Research (ANR-2010-TECS-012). The authors would like to thank the persons who agreed to participate in recordings.


Evgenios N Kornaropoulos
email: evgenios.kornaropoulos@ecp.fr

Evangelia I Zacharaki
email: evangelia.zacharaki@ecp.fr

Pierre Zerbib
email: pierre.zerbib@aphp.fr

Chieh Lin

Alain Rahmouni
email: alain.rahmouni@aphp.fr

Nikos Paragios
email: nikos.paragios@ecp.fr

Optimal estimation of diffusion in DW-MRI by high-order MRF-based joint deformable registration and diffusion modeling

Over the last years, the apparent diffusion coefficient (ADC), computed from diffusion-weighted magnetic resonance (DW-MR) images, has become an important imaging biomarker for evaluating and managing patients with neoplastic or cerebrovascular disease. Standard methods for the calculation of ADC ignore the presence of noise and motion between successive (in time) DW-MR images acquired by changing the b-value. In order to accurately quantify the diffusion process during image acquisition, we introduce a method based on a high-order Markov Random Field (MRF) formulation that jointly registers the DW-MR images and models the spatiotemporal diffusion. Spatial smoothness on the ADC map, as well as spatiotemporal deformation smoothness, is imposed towards producing anatomically meaningful representations. The high-order dependencies in our MRF model are handled through Dual Decomposition. Performance of registration is compared to a state-of-the art registration approach in terms of obtained fitting error of the diffusion model in the core of the tumor. Preliminary results reveal a marginally better performance of our method when compared against the standard ADC map used in clinical practice, which indicates its potential as a means for extracting imaging biomarkers.

Introduction

Since the proposal of background body signal suppression [START_REF] Takahara | Diffusion weighted whole body imaging with background body signal suppression (dwibs): technical improvement using free breathing, stir and high resolution 3d display[END_REF], DW-MRI has gained clinical impetus for the detection of primary tumors and metastases throughout the whole body [START_REF] Koh | Diffusion-weighted mri in the body: applications and challenges in oncology[END_REF], [START_REF] Padhani | Diffusion-weighted magnetic resonance imaging as a cancer biomarker: consensus and recommendations[END_REF]. Besides the qualitative anatomical information [START_REF] Balbo-Mussetto | Whole-body mri with diffusionweighted imaging: a valuable alternative to contrastenhanced ct for initial staging of aggressive lymphoma[END_REF], [START_REF] Albano | Comparison between whole-body mri with diffusion-weighted imaging and pet/ct in staging newly diag-nosed fdg-avid lymphomas[END_REF], [START_REF] Vilanova | Diffusion-weighted wholebody mr screening[END_REF], [START_REF] Schmidt | Wholebody mri for the staging and follow-up of patients with metastasis[END_REF], quantitative analysis of whole-body DW-MRI is becoming possible by scanning the patient with different b-values (a parameter determining the strength and duration of the diffusion gradients). The apparent diffusion coefficient, which reflects the gradient of water diffusivity in the body, can be subsequently computed at each voxel and displayed as a parametric map (ADC map) [START_REF] Le Bihan | Separation of diffusion and perfusion in intravoxel incoherent motion mr imaging[END_REF]. ADC has been examined not only for the characterization and staging of lesions but also for the response to treatment [START_REF] Lin | Whole-body diffusion magnetic resonance imaging in the assessment of lymphoma[END_REF], [START_REF] Gu | Whole-body diffusion-weighted imaging: the added value to whole-body mri at initial diagnosis of lymphoma[END_REF], [START_REF] Wu | Diffusionweighted mri in early chemotherapy response evaluation of patients with diffuse large b-cell lymphoma-a pi-lot study: comparison with 2-deoxy-2-fluoro-d-glucosepositron emission tomography/computed tomography[END_REF], [START_REF] Lin | Whole-body diffusion-weighted magnetic resonance imaging with apparent diffusion coefficient mapping for staging patients with diffuse large b-cell lymphoma[END_REF]. The calculation of ADC was performed without motion or artifact correction in the previous studies. It is reported though that ADC is a parameter susceptible to artifacts, the most frequent of all being patient's motion and breathing, resulting in misregistration of the images obtained with different b-values [START_REF] Albano | Comparison between whole-body mri with diffusion-weighted imaging and pet/ct in staging newly diag-nosed fdg-avid lymphomas[END_REF], [START_REF] Kwee | Whole-body diffusionweighted imaging for staging malignant lymphoma in children[END_REF]. Therefore, in order to correctly visualize diffusion in biological tissues, image registration should precede the ADC calculation [START_REF] Koh | Diffusion-weighted mri in the body: applications and challenges in oncology[END_REF]. Specifically, due to non-linear motion (cardiac pulses, respiration), non-rigid alignment of the DW-MRI sequence becomes necessary [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF].

In this paper, we propose a joint deformable registration and appearance model computation framework. The deformation model maps each one of the m DW images in a common space (the ADC space), whereas the appearance model describes the signal intensity related to the diffusion of water molecules that takes place during DW-MRI acquisition. The unknown variables of the combined deformation and appearance model correspond to a m + 2-layer graphical model, m layers that represent the m 3D deformation fields, one for each DW image, and two layers for the two parameters of the appearance model, namely the 3D ADC map and the 3D signal intensity without diffusion weighting (explained later with more details). The latent deformation variables are obtained using the discrete formulation introduced in [START_REF] Glocker | Dense image registration through mrfs and efficient linear programming[END_REF]. The m + 2 layers are interconnected in order to achieve consistency between the predicted diffusion (by the appearance model) and the diffusion calculated by the registered DW images. Each deformation field depends thus on the 2-layers appearance model, which amounts to a high-order dependency that requires a high-order inference scheme. Spatial smoothness constraints on the parameters of the appearance model and the deformations are imposed through pairwise interactions. Moreover, the m layers deformation fields have pairwise interconnections (based on proximity of b-values) towards imposing smoothness on the change of deformation through time (the images sequence is acquired with increasing b-values).

Our approach lies between frameworks of concurrent segmentation and registration [START_REF] Wyatt | Map mrf joint segmentation and registration of medical images[END_REF], [START_REF] Yezzi | A variational framework for integrating segmentation and registration through active contours[END_REF], [START_REF] Parisot | Joint tumor segmentation and dense deformable registration of brain mr images[END_REF], [START_REF] Alchatzidis | Local atlas selection for discrete multi-atlas segmentation[END_REF] and image reconstruction and modeling [START_REF] Rajan | Simultaneous estimation of super-resolved scene and depth map from low resolution defocused observations[END_REF], [START_REF] Park | Super-resolution image reconstruction: a technical overview[END_REF]. The main difference of our approach with respect to the concurrent segmentation and registration works is that, our ADC variables are continuous and have unique values on a voxel-wise basis, whereas segmentation labels are discrete and common within regions. In that sense, methods for image reconstruction [START_REF] Rajan | Simultaneous estimation of super-resolved scene and depth map from low resolution defocused observations[END_REF] and denoising [START_REF] Pižurica | Estimating the probability of the presence of a signal of interest in multiresolution singleand multiband image denoising[END_REF] are closer to our approach; those however are defined in a single spatial domain and do not require deformable image registration. On the other hand, methods for recovering appearance change and motion in image sequences target a similar goal, but they account usually for global motion or smooth deformations, as well as systematic changes in appearance exhibiting spatial or temporal structure that can be modeled [START_REF] Black | A framework for modeling appearance change in image sequences[END_REF]. Intensity changes in domain-specific cases are captured by learning models using training examples. Such training images are not available in our case where the tumor appears in different areas of the body.

By elaborating more on concurrent segmentation and registration frameworks, we aim at finding the optimal deformation fields of the DW-MR images together with the optimal ADC values in a one-shot optimization approach. Alchatzidis et al. [START_REF] Alchatzidis | Local atlas selection for discrete multi-atlas segmentation[END_REF] proposed a joint group-wise registration and segmentation method for brain atlas selection, but they optimized the deformation fields iteratively, whereas Sotiras et al. [START_REF] Sotiras | Graphical models and deformable diffeomorphic population registration using global and local metrics[END_REF] didn't include a segmentation layer in their group-wise, one-shot registration method. In our case the one-shot optimization problem imposes additional challenges. To this end, Dual Decomposition (DD) for MRF optimization [START_REF] Komodakis | Mrf energy minimization and beyond via dual decomposition[END_REF] is used to determine the optimal solution to the joint problem, in which linear programming and duality [START_REF] Komodakis | Performance vs computational efficiency for optimizing single and dynamic mrfs: Setting the state of the art with primal-dual strategies[END_REF] are incorporated. It is also worth to mention that our framework is modular with respect to the registration/annotation grid and can be adapted to the application domain. Finally, to the best of our knowledge, we are the first to propose a DW-MRI registration scheme that optimizes explicitly the derived ADC map. Very recently, two other works of DW-MRI registration were published [START_REF] Guyader | Influence of image registration on apparent diffusion coefficient images computed from free-breathing diffusion mr images of the abdomen[END_REF], [START_REF] Veeraraghavan | Simultaneous segmentation and iterative registration method for computing adc with reduced artifacts from dw-mri[END_REF] aiming at computing the ADC without motion artifacts. Guyader at al. [START_REF] Guyader | Influence of image registration on apparent diffusion coefficient images computed from free-breathing diffusion mr images of the abdomen[END_REF] examined the influence of pairwise regis-tration on the ADC calculation, while Veeraraghavan et al. [START_REF] Veeraraghavan | Simultaneous segmentation and iterative registration method for computing adc with reduced artifacts from dw-mri[END_REF] proposed a group-wise DW-MRI registration method that requires structure segmentation, thus is not completely automated. Both methods (i) use mutual information criterion for image similarity which is less appropriate than the fitting error to a diffusion model we are using. (ii) calculate the motion for each image independently without considering temporal consistency (smooth transition in time) and (iii) extract the ADC by voxel-wise curve fitting from the warped images without accounting for spatial constraints.

The remainder of this paper is organized as follows: section 2 describes the combined registration and ADC calculation strategy, while experimental validation is presented in the following section. Discussion and future directions conclude the paper.

Graph-based joint DW image registration and diffusion modelling

Let us consider a sequence of m DW-MR images, each one described by intensity values s t (x), t = 1, ..., m, with x ∈ Ω t , Ω t ∈ R 3 . On top of that, let's consider an extra image z(x) with x ∈ Ω z which represents the ADC template corresponding to the previous image sequence. This ADC template is regarded as the reference frame of an optimal alignment among the DW-MR images. We therefore seek for two results: a set of deformations d = {d : x R = d t (x t ), t = {1, ..., m}} which map mutually corresponding points from the m-image spaces to the same point of the reference frame Ω z and the optimal -in terms of spatial consistency and smoothness -ADC values that can be calculated when all m-images are mutually aligned. This joint mapping and ADC refinement aims at (i) correcting deformations due to (local) organ deformations or motion, (ii) ensure temporal consistency in the diffusion process and (iii) impose spatial consistency on the deformation fields and the derived ADC map. The appearance model used for the calculation of the ADC, as well as the deformation model, are jointly optimized to define the reference pose. In the following section, the problem formulation is presented, first for the case of negligible deformations and then for all possible deformations.

Appearance (diffusion) model and MAP estimation

In DWI, we assume that the intensity values of the images depend mainly on the amount of signal loss occurred during the application of two equally rectangular gradients before and after the 180 • -refocusing pulse [START_REF] Bammer | Basic principles of diffusion-weighted imaging[END_REF]. The relation between the signal attenuation expressed by the signal intensity s and the diffusion of the water molecules in a microscopic level is known to be exponential [START_REF] Le Bihan | Separation of diffusion and perfusion in intravoxel incoherent motion mr imaging[END_REF]. Having acquired the images with high b values (> 50s/mm 2 ), we can assume that the DW images are insensitive to vascular capillary perfusion and diffusion is the only type of motion present, thus the ADC is equal to the water diffusion [START_REF] Koh | Diffusion-weighted mri in the body: applications and challenges in oncology[END_REF]. Therefore the relation between the ADC (referred to with the letter z) and the signal intensity s can be modeled as:

s t = s 0 e -bt•z (1)
where s t and b t the image intensity vector and b-value at instant t respectively, s 0 the signal attenuation for b = 0 and z is the ADC vector (all vectors are linearized 3D matrices).

Denoting with y t the natural logarithm of the image vector s t , the relation between the image vector and the ADC can be expressed as:

y t = -b t z + y 0 + n t (2)
where n t the noise vector. We are assuming noise to be zero mean i.i.d., thus the multivariate pdf of n t is given by:

P (n t ) = 1 (2π) M 2 σ M n exp{- 1 2σ 2 n n T t n t } (3)
where M the size of the linearized noise vector and σ 2 n denotes the variance of the noise process. If all m DW-MR images are perfectly aligned, the standard way of estimating ADC is by calculating the slope of a line fitted on the plot of the logaritm of the relative signal intensities of the tissue against the b-values:

z = m i=1 (b i -b)((y i -ȳ) m i=1 (b i -b) 2 (4)
where b, ȳ the mean b-value and mean of the natural logarithm of the image vector respectively. We approach the calculation of z, as well as y 0 using the Maximum a Posteriori (MAP) technique:

ẑ, ŷ0 = argmax z,y0 P(z, y 0 | y 1 , • • • , y m ) (5)
The random fields z and y 0 are assumed to be statistically independent in this study, as they refer to two independent processes. From Bayes' rule, equation ( 4) can be written as:

ẑ, ŷ0 = argmax z P(y 1 , • • • , y m | z, y 0 )P(z)P(y 0 ) P(y 1 , • • • , y m ) (6)
Since the denominator is not a function of z nor of y 0 , it does not influence their estimation and therefore can be ignored. Taking the log probability of posterior probability:

ẑ, ŷ0 = argmax z,y0 {[log P(y 1 , • • • , y m | z, y 0 ) + log P(z) + log P(y 0 )]} (7)
Next we solve equation ( 7) assuming there is no deformation between the DW images.

From equation [START_REF] Glocker | Dense image registration through mrfs and efficient linear programming[END_REF], since n t 's are assumed to be independent:

ẑ, ŷ0 = argmax z,y0 {[log m t=1 P(y t | z, y 0 ) + log P(z) + log P(y 0 )]} = argmax z,y0 {[ m t=1 log P(y i | z, y 0 ) + log P(z) + log P(y 0 )]} (8)
Since noise is assumed to be i.i.d. Gaussian, from equations ( 2) and ( 3) we have:

P(y 1 , • • • , y m | z, y 0 ) = [ m t=1 log 1 (2πσ 2 n ) M 2 exp{- y t + b t z -y 0 2 2σ 2 n }] = m t=1 -y t + b t z -y 0 2 2σ 2 n - M 3 log(2πσ 2 n ) (9)
Substituting equation ( 9) in equation ( 8) we obtain:

ẑ, ŷ0 = argmax z,y0 {[ m t=1 -y t + b t z -y 0 2 2σ 2 n + log P(z) + log P(y 0 )]} (10)
The derivation of a solution for ADC based on equation [START_REF] Koh | Diffusion-weighted mri in the body: applications and challenges in oncology[END_REF] is the first contribution of our work, in which the ADC can be calculated based on the data (DW-MR images) and some prior knowledge about context dependencies on ADC map and y 0 values, acting as regularization constraints.

MAP estimation for joint appearance and deformation model

Assuming now that some motion has occurred during the acquisition process of the DW-MR images, equation (2) takes the form:

y t • D t = -b t z + y 0 + n t ( 11
)
where D t denotes the deformation field that maps s t to the ADC template space. In this case, we would like to compute the m deformation fields that would bring the m DW-MR images to a common reference space in which the ADC map can be accurately estimated. Let's denote the geometric deformation that maps the DW-MR image s t to the common (ADC) space with D t (x) = x + d t (x), where d t is the displacement field of image s t . The joint posterior probability for the appearance model parameters and the set of

displacement fields [d 1 , • • • , d m ]
is given by:

P(d 1 , • • • , d m , z, y 0 | y 1 , • • • , y m ) = P(y 1 , • • • , y m | d 1 , • • • , d m , z, y 0 )P(d 1 , • • • , d m , z, y 0 ) P(y 1 , • • • , y m ) (12)
The random fields d t ,z and y 0 are once again assumed to be statistically independent, as they refer to three independent processes. Therefore, the joint MAP solution could be derived similarly as in the previous section, resulting in the following optimization problem:

ẑ, ŷ0 , d1 , • • • , dm = argmax z,y0,d1,••• ,dm {[ m t=1 -y t • D t + b t z -y 0 2 2σ 2 n + log P(z) + log P(y 0 ) + log P(d 1 , • • • , d m )]} (13)
The latter equation for computing a joint solution for ADC and deformation fields, using some priors on context dependencies on z, y 0 and deformation fields, is the second contribution of our work.

MRF Formulation of Joint Problem

We formulate our joint DW-MRI registration and ADC modeling problem using Markov Random Fields (MRF). The joint model is parameterized by a set of m + 2 isomorphic grid graphs G = {G 1 , . . . , G m+2 }, each of the first m superimposed onto the corresponding DW-MR image, whereas the last two grid graphs are superimposed onto the ADC map that we want to compute and the y 0 image vector respectively. To this end, we define three sets of labels, corresponding to the three different types of variables in our model. The label set L d is a quantized version of the displacement space. A label assignment l d p ∈ L d (with p ∈ V d ) is equivalent to displacing the control point p by displacement d p . Correspondingly, the labels l z p ∈ L z and l y0 p ∈ L y0 are quantized versions of the ADC and y 0 intensities, respectively.

For every control point in each grid i there is a node p that represents either its displacement if the grid is a deformation grid either the ADC or the y 0 value in case of the last two grids respectively (p ∈ V,with V denoting the set of all nodes on a grid that encode the latent variable). Regarding the deformation, it is evaluated on the deformation grid's control point p as D(x) = x + p∈G η( x -p )d p , where d p is the displacement of the grid's control point p. By moving the grid's control points, we end up deforming the DW-MR images. In the end, the displacement of a voxel x of the image will be determined by the control points displacements and the influence of each control point on x, which is given by the projection function η (a cubic B-splines in our case). Due to the isomorphic nature of the grid, p indexes a common to all grids control point position. Moreover, the nodes in the graph are connected with a set of edges E that encode the interactions between the deformation and the ADC variables. The first type of dependencies between the nodes is a high-order dependency term that relates to the coupling term, which penalises any deviation due to registration of the image intensity value from the attenuation value calculated based on equation [START_REF] Albano | Comparison between whole-body mri with diffusion-weighted imaging and pet/ct in staging newly diag-nosed fdg-avid lymphomas[END_REF]. The second type of dependencies in our energy formulation is related to temporal constistency of the deformation fields, which means that we assume that the deformations occurring on the image query should change smoothly from one image to the other. This term is applicable only on edges connecting nodes on different deformation grids (inter-deformation dependency). The third, fourth and fifth terms relate to the smoothness constraints and are applicable only on edges connecting nodes on the same grid (intra-deformation, intra-ADC and intra-b0 dependency respectively) (see figure 1,left).

Energy terms: The constructed graph is associated with an energy, consisting of four pairwise terms:

E joint (d 1 , • • • , d m , z, y 0 ) = λ c m t=1 p∈Gt p∈Gm+1 p∈Gm+2 V c (l dt p , l z p , l y0 p ) + λ inter m-1 t=1 p∈Gt V inter (l dt p , l dt+1 p ) + λ intra 1 m t=1 p∈Gi q∈(N (p)∩Gi) V intra 1 (l dt p , l dt q ) + λ intra 2 p∈Gm+1 q∈(N (p)∩Gm+1) V intra 2 (l z p , l z q ) + λ intra 3 p∈Gm+2 q∈(N (p)∩Gm+2) V intra 3 (l y0 p , l y0 q ) ( 14
)
where V c is the coupling potential of a third-order dependency, V inter is a pair-wise potential for temporal regularization, Coupling term: It penalizes deformations that lead to an increase in the matching error among z, y 0 and the signal attenuation at a given instant t, based on equation [START_REF] Kwee | Whole-body diffusionweighted imaging for staging malignant lymphoma in children[END_REF]. It is applicable on a third-order clique consisted by the corresponding control points p (same spatial location) on the deformation, the z and the y 0 grid:

V intra n (•, •), n = 1, • • • ,
V c (l dt p , l z p , l y0 p ) = Ωs t η( x t -p ) (y(x t + d l d t p ) + b t z l z p -y l y 0 p 0 ) 2 dx t ( 15
)
This term is our data term as it connects the DW-MR images with all three types of variables. It is defined on the image domain Ω and a function η is used to back project the voxel-wise information on the grids' control points, by determining how much voxel x influences the control point p.

Spatial and temporal deformation smoothness: We assume that the organs' motion evolves smoothly during image acquisition process therefore we expect the deformation on each DW-MR image to not deviate much from the deformation computed on each immediate neighbor, meaning the image preceding it. In the same time, each deformation field should deform the tissues in a way in which the anatomical information is retained. Therefore in order to impose temporal and spatial smoothness on the deformation fields we penalize the magnitude of displacement vector differences. The pairwise potential would then be:

V (l di p , l dj q ) = d l d i p -d l d j q p -q (16)
In case of temporal deformation smoothness (V inter ) p ≡ q and j = i + 1, whereas in case of spatial deformation smoothness (V intra 1 ) i ≡ j. For both cases we are assuming a 8-nodes connectivity. Spatial smoothness on appearance model: Finally the last two terms in equation ( 14) relate to the spatial constraints applied within each of the images of our appearance model (z and y 0 ). The squared error between two different appearance labels was used towards this end:

V intra 2 (l z p , l z q ) = (z l z p -z l z q ) 2 , V intra 3 (l y0 p , l y0 q ) = (y 0 l y 0 p -y 0 l y 0 q ) 2 (17)
Again we are assuming a 8-nodes intra-image connectivity.

Optimization and implementation details

Our discrete high-order MRF-based energy minimization problem was solved using Dual Decomposition (DD), which has been introduced as a framework for MRF optimization, offering global guarantees [START_REF] Komodakis | Mrf energy minimization and beyond via dual decomposition[END_REF]. Dual Decomposition consists in decomposing a global difficult problem into smaller solvable subproblems (referred to as slaves) and then extracting a solution by combining the solutions from these subproblems. In our framework, the graph is decomposed into subgraphs that constitute the set of subproblems, such that they cover (at least once) every node and edge in each graph. One slave problem is defined for each line parallel to a coordinate axis in space and time (corresponding to pairwise interactions), and one slave is defined for the higher order clique (related to the data term). Two different optimization schemes are used for the solution of the subproblems: linear programming and duality [START_REF] Komodakis | Performance vs computational efficiency for optimizing single and dynamic mrfs: Setting the state of the art with primal-dual strategies[END_REF] for pairwise interactions and exhaustive search for the higher order clique. In the end, global optimization is achieved through DD, that ensures agreement among the solutions provided by all slaves.

Our algorithm follows a coarse-to-fine process that first estimates parameters at a coarse spatial image resolution level and then updates the solution to a finer resolution level. We applied a multi-resolution scheme with 3 resolutions levels; two levels in which the joint problem is solved and one at full image resolution in which only the appearance model is optimized, keeping the previously obtained deformation fields fixed. The rationale behind this multiresolution approach is that we want to adopt the grid size to the information content (frequencies) of the modeled variables, i.e. a sparse grid for the deformation model to cope with potential rigid transformations, a mid-level grid for local deformation and global intensity variations and a last dense grid for local intensity changes. Finally, regarding the run-time performance of the method, all the experiments were computed on a 8-cores machine (3.6 GHz processing power) with a memory requirement of 4 GB.

Experimental Validation

In order to validate the proposed registration method, we used a data set provided by the Centre Hospitalier Universitaire Henri-Mondor, France. The data set consisted of 3D images of 17 different patients with lymphoma, scanned with a 3 Tesla MRI scanner using 5 b values: 20,50,150,400 and 800 (s/mm 2 ). Therefore, the sequence of the images to be registered was always of size 5 (m = 5) per patient. Both male and female patients were included. The image size was 130x106x22 for every patient, with a pixel resolution of 2.46x2.46 mm 2 in the axial plane and 5 mm 2 in the z-axis. Each image had been manually annotated by experts to indicate the contour of the tumor (see figure 2).

We compared the proposed joint method against three other approaches of computing an ADC map. The standard way of computing an ADC map is by plotting the logarithm of the relative signal intensity of the tissue on the y-axis against the b values on the x-axis and then fitting a line through the plots with the use of linear regression (see equation ( 4)). Then the ADC is described by the slope of the line. This standard approach used in clinical practice does not take into account any motion correction scheme, therefore we refer to the ADC map derived by this approach as the unregistered case. We also compared the ADC derived by the proposed joint method against an ADC map computed based on deformed images, following the afore-mentioned ADC computation approach as a post-processing to the registration step. An implementation of a state-of-the art group-wise registration approach [START_REF] Sotiras | Graphical models and deformable diffeomorphic population registration using global and local metrics[END_REF] was examined and the ADC computed based on the derived deformed images was calculated. An entropy estimator based on spacings [START_REF] Bouzebda | New entropy estimator with an application to test of normality[END_REF] was used as a global similarity criterion and mutual information (MI) was used in the inter-images comparisons. Finally, we computed an ADC map which values were extracted from the wraped images by our approach.

The ADC map in all cases except the case where it is computed explicitely by the proposed method was computed based on the DW-MR images (deformed or not). Then the accuracy on the ADC maps computed by these four different approached was quantitatively assessed by computing the sum-squared error (SSE) among the predicted image values derived by equation ( 1) and the actual (deformed or not) image values. (see Fig. 1,right). The mean SSE over all voxels belonging to a tumor in the case of the explicitely computed by the proosed method ADC was lower than the corresponding mean SSE computed by all three rest cases.

Finally, in order to qualitatively evaluate the results, the ADC map for the same sagital view of a patient is shown in Fig. 2a, as derived by the four different cases, as well as the s 0 image computed also explicitely by our method. Fig. 2 (first row) illustrates an example of the DW-MR intensity variation for the five different b-valuess. In order to visually assess the results, we plot the ADC map of the same patient as calculated by four different methods. The results are illustrated in Fig. 2, second row, and correspond (from left to right) to the ADC calculated by voxel-wise curve fitting from (i) unregistered images, (ii) warped images using entropy-based similarity (implementation of the method in [START_REF] Sotiras | Graphical models and deformable diffeomorphic population registration using global and local metrics[END_REF]), and (iii) warped images based on diffusion model fit (proposed approach without ADC and y 0 regularization constraints). Finally the AC map an the s 0 image resuted from the propoed joint framework are illustated in the fourh and fith colum respectivly.

Conclusion

In this paper, we propose a novel joint deformable registration based on diffusion modeling that derives an optimal geometrical as well as appearance representation towards more accurately depicting the physiological process taking place during DW image acquisition. At that point we would like to note that convergence to a consenus solution among the slaves in the DD scheme was not achieved in all the experiments. In those cases, the majority voting algorithm was applied to derive the final solution to the energy minimization problem. Extensive study with more data and a parametric analysis of the method are required in order to better explore the predictive power of ADC. Overall, we conclude on the usefulness of DW-MRI registration towards improving the quality of ADC, which is in accordance with the conclusions derived by recent works on ADC refinement through image registration.

3

 3 are the three pairwise potentials imposing spatial regularization within each graph, and N represents the spatial connectivity of the nodes. The constants λ c ,λ inter and λ intra n , n = 1, • • • , 3 are used to weight the different types of potentials.

Figure 1 :

 1 Figure 1: Left figure: The node and the edge system of all the connected graphs. With green, blue and brown color the relationship within the grid nodes belonging to ADC, y 0 and each of the deformation field respectively (intra connections). The black dotted edge represents the third-order dependency among deformation field and appearance model (coupling relationship), whereas pink edges denote the temporal relationship between successive deformatio fields (inter relationship). Right figure: Boxplots of the mean SSE based on ADC maps derived by the four different examined cases.

Figure 2 :

 2 Figure 2: An axial slice of the 3D image of a patient. First row: b values of 20, 50, 150, 400 and 800 s/mm 2 (from left to right). The red contour in the "b 50" image corresponds to the tumor's contour, as it was drawn by the expert. 2nd row: ADC map derived by unregistered images,by deformed entropy-based group-wise registration approach, by curve fitting from the wraped by the proposed framework images and finally the ADC derived by the proposed joint approach (from left to right).The last image on the right is the s 0 derived also explicitly by our framework.

Aknowledgements

The authors of this paper would like to thank Enzo Ferrante and Maxim Berman for their helpful and constructive comments on the development of this work.


Evgenios N Kornaropoulos

Evangelia I Zacharaki

Pierre Zerbib

Chieh Lin

Alain Rahmouni

Nikos Paragios

DEFORMABLE GROUP-WISE REGISTRATION USING A PHYSIOLOGICAL MODEL: APPLICATION TO DIFFUSION-WEIGHTED MRI

Keywords: Group-wise registration, DW-MRI, ADC

Intensity variations can often be described by a physiological or temporal model applied on a voxel-wise basis across a group of images. However the voxel correspondence might be unknown, imposing the need for a group-wise deformable registration coupled with the computation of the model parameters. In this paper we propose a group-wise registration method of medical images that incorporates the temporal dimension (reflecting the change of signal amplitude) of the acquisition process. Consistency on the spatiotemporal physiological model, as well as deformation smoothness, is imposed in order to produce anatomically meaningful representations of the 3D images. The performance of the proposed method is compared to two different group-wise registration approaches; one that penalizes the absolute differences in the intensities and one that penalizes the intensity range among the images on corresponding regions. We chose as an application paradigm the registration of diffusion-weighted magnetic resonance (DW-MR) images for the evaluation of patients with lymphomas. A dataset consisting of 25 patients, each scanned with 3 "b values", was used to evaluate the method's accuracy. The proposed registration method outperfomed the other two registration approaches, making it a very promising method for highlighting the importance of DWI as an imaging biomarker.

INTRODUCTION

Deformable image registration is a very active area of research in medical imaging. It constitutes a very challenging problem, as one has to cope with motion correction due to organs' deformations occuring over time or whole body movement [START_REF] Sotiras | Deformable medical image registration: A survey[END_REF]. Finding anatomical correspondences becomes even more challening when intensity changes due to different acquisition parameters or noise and technical artifacts. In terms of image registration of longitudinal scans, various population registration techniques have been proposed towards finding a homology among a population of more than two im-ages. Based on the defined homology, the images are either registered to a selective or constructed reference frame or in a template-free group-wise manner, the former having the limitation of introducing bias to the process through the selection of the reference, whereas the latter approach is computationally demanding [START_REF] Sotiras | Graphical models and deformable diffeomorphic population registration using global and local metrics[END_REF].

Over the last years, diffusion-weighted magnetic resonance imaging (DWI) has become an important modality for evaluating and managing patients with neoplastic or cerebrovascular disease. Quantitative analysis of DW-MRI is becoming possible by scanning the patient with different values of a parameter that characterizes the applied diffusionsensitizing (bipolar) gradients, known as "b-value" [START_REF] Le Bihan | Separation of diffusion and perfusion in intravoxel incoherent motion mr imaging[END_REF]. The apparent diffusion coefficient (ADC), which depicts the microscopic mobility in time of water molecules in the body, can be then computed for each voxel and displayed as a parametric map (ADC map). Retention of the spatial information derived by the ADC map requires registration of DW-MR images [START_REF] Koh | Diffusion-weighted mri in the body: applications and challenges in oncology[END_REF].

In this paper, we propose a group-wise registration scheme without the need of choosing a reference template while in the same time it is computationally efficient. We aim at finding the optimal deformation fields of the diffusionweighted magnetic resonance (DW-MR) images using a temporal constraint, related to the diffusion process, as well as a smoothness penalty on the deformations. The latent variables of the model are n-deformations (B-spline polynomials) of the images, which are obtained using the discrete formulation introduced in [START_REF] Glocker | Dense image registration through mrfs and efficient linear programming[END_REF]. The deformation variables are connected with the observations towards ensuring meaningful temporal correspondence among the DW images. They are also inter-connected in order to decrease the cost of pairwise comparisons between individual images. Linear programming and duality [START_REF] Komodakis | Performance vs computational efficiency for optimizing single and dynamic mrfs: Setting the state of the art with primal-dual strategies[END_REF] are used to determine the optimal solution of the problem. Our work was inspired by the approach of [START_REF] Sotiras | Graphical models and deformable diffeomorphic population registration using global and local metrics[END_REF]. The main difference between the two approaches is found in the encoding of the global similarity of the population. A statistical compactness criterion has been used in [START_REF] Sotiras | Graphical models and deformable diffeomorphic population registration using global and local metrics[END_REF], whereas we incorporate within the registration process a physiological model representing the temporal intensity change. This inten-sity modeling step removes the ambiguity during the search of anatomic correspondences and thus increases the mapping accuracy.

The application of group-wise registration towards calculating the ADC without motion artifacts has also been examined in [START_REF] Guyader | Influence of image registration on apparent diffusion coefficient images computed from freebreathing diffusion mr images of the abdomen[END_REF] and [START_REF] Veeraraghavan | Simultaneous segmentation and iterative registration method for computing adc with reduced artifacts from dw-mri[END_REF]. However, the registration process in the former is based on pair-wise connections between the images, whereas in the latter one the proposed group-wise registration requires prior structure segmentation, making it a not fullyautomated framework. Finally, both methods use mutual information criterion for image similarity, without taking into account the intrinsic intensity changes in DW-MRI due to the diffusion process.

MATERIALS AND METHODS

Let us consider n images {I 1 ,..,I n }, where each image is described by intensity values I i (x), with x ∈ Ω i , Ω i being the domain of image i (R 3 ). We are looking for a set of transformations T = {T : x R = T i (x i ), i = {1, ..., n}} which map mutually corresponding points from the n-image spaces to the same point of a reference frame Ω R . This mapping aims to correct deformations due to anatomical differences or motion, in order to identify intensity differences in time due to change in the measured quantity. The reference pose would then correspond to the consensus geometrical space, in which the ADC map resides. The models used for the deformation and the diffusion process are described next.

Deformation Model

Similarly to [START_REF] Alchatzidis | Local atlas selection for discrete multi-atlas segmentation[END_REF], we used the popular Free Form Deformations transformation model to parametrize each transformation T (x) by a linear combination of K control points, weighted using cubic B-splines. In the end, we have n uniformly distributed grids of control points over the image domain (one corresponding to each 3D image).

Physiological Model

We assume that the intensity values of the images depend mainly on the amount of signal loss occured during the application of the two equal rectangular gradients before and after the 180 • -refocusing pulse [START_REF] Bammer | Basic principles of diffusion-weighted imaging[END_REF]. The relation between the signal attenuation expressed by the signal intensity Î and the diffusion of the water molecules in a microscopic level is well known to be exponential [START_REF] Le Bihan | Separation of diffusion and perfusion in intravoxel incoherent motion mr imaging[END_REF]. Having acquired the images with high "b-values" (> 50s/mm 2 ), we can assume that our diffusion-weighted data are insensitive to vascular capillary perfusion and diffusion is the only type of motion present, thus the ADC is equal to the water diffusion [START_REF] Koh | Diffusion-weighted mri in the body: applications and challenges in oncology[END_REF]. Therefore the relation between ADC and the the signal intensity at each voxel location can be modeled as:

Îi (x) = e -bi•ADC(x) (1)
where b i is the b value used for the acquisition of image Îi . In case of n "b-values" images, the ADC can be calculated as the slope of a line fitted on the plot of the logarithm of the relative signal intensities of the tissue against the "b-values".

The ADC would then be:

ADC(x) = n i=1 (b i -b)((S i (x) -S(x)) n i=1 (b i -b) 2 (2)
where S i = ln Îi and b, S the mean "b-value" and mean of ln Î respectively.

MRF Formulation

We formulate our deformable registration problem using Markov Random Fields, in which the deformation model is parametrized by a set of n isomorphic grid graphs G = {G 1 , . . . , G n }. For every control point in the deformation grid superimposed onto image i, there is a node p i that represents its displacement (p i ∈ V,with V denoting the set of all nodes that encode the latent variables). Due to the isomorphic nature of the grid, p indexes a common to all grids control point position. Moreover, the nodes in the graph are connected with a set of edges E that encode the interactions between the deformation variables. The first type of dependencies between the nodes relates to the fitting error of the physiological model and is applicable only on edges connecting nodes on different grids (inter-image dependency).

The second type of dependencies relates to the smoothness constraint and is applicable only on edges connecting nodes on the same grid (intra-image dependency). Finally, the solution space around a control point is quantized and indexed by a set of labels L that represents possible control point displacements.

Energy terms: The constructed graph is associated with an energy, consisting of three terms, a global unary term and two pairwise terms:

E M RF (T 1 • G 1 , • • • , T n • G n ) = λ global n i=1 p∈Gi V p (l p ) + λ inter n i=1 p∈Gi q∈(N (p)\Gi) V pq (l p , l q ) + λ intra n i=1 p∈Gi q∈(N (p)∩Gi) V pq (l p , l q ) (3)
where V p (•) are the unary potentials,V pq (•, •) are the pair-wise potentials and N represents the connectivities of the nodes. Three different constants, λ global ,λ inter and λ intra are used to weight the three different types of potentials. The optimal values of the afore-mentioned weights were empirically determined in a brute force approach. The optimization of such an energy is non-convex, however the optimizer we used provides guarantees for an almost optimal solution to such an NP-hard MRF problem [START_REF] Komodakis | Performance vs computational efficiency for optimizing single and dynamic mrfs: Setting the state of the art with primal-dual strategies[END_REF].

Global Unary Term: We approximate a global cost for the deformations by assuming that for each node p of a given deformation field/image i, the rest of the images do not get deformed within the current iteration. This assumption is considered for all nodes and for all deformation fields within a given iteration, an assumption that is common in minimizing graphical models through expansion moves. In this case, the cost of a deformation depends only on the label of this node:

V t pi (l pi ) ≈ • • • Ω1 ••• Ωn ωpi (x i ) g(I 1 (T t-1 1 (x 1 )), • • • , I i (T t i (x i )), • • • , I n (T t-1 n (x n ))) dx 1 • • • dx n (4)
where ω(•) is the support function that determines the contribution of the point x on image i to the unary potential of the control point p of the corresponding grid i. By g(•) we define the metric used to assess a global cost for each possible displacement l of a control point. We propose as an objective criterion to be minimized for the case of DWI-MRI registration the exponential regression fitting error among the images in the sequence, in order to impose consistency on the computation of the physiological model. It should be expected that as the images are jointly aligned, the derived optimal pose would express more accurately the diffusion process, enabling the computation of more accurate ADC maps. More specifically, we propose the sum of squared errors of prediction (SSE) for the global measurement towards global registration in DW-MR images. In that case the metric g would be:

g(I 1 , • • • , I n ) = n i=1 (I i -Îi ) 2 ( 5
)
where Îi is calculated by Eq.1 for the current estimate (fitted values) ADC.

Inter-image Pair-wise Term: We consider a second term related to the obervations, meaning the intensities of the images, that aims to reduce the outlier effect. To this end, local pair-wise comparisons between members of the population are examined too. We assume that two images (for example images i and j) are allowed to be deformed, while the rest n -2 remain static within the current iteration. The SSE criterion is used to calculate the cost of a deformation, as expressed by the following pair-wise term:

V piqj (l pi , l qj ) ≈ • • • Ω1 ••• Ωn ωpi (x i ) [(I 1 (T t-1 1 (x 1 )) -f (x 1 )) 2 + • • • + (I i (T t i (x i )) -f (x i )) 2 + (I j (T t j (x j )) -f (x j )) 2 + • • • + (I i (T t-1 n (x n )) -f (x n )) 2 ] dx 1 • • • dx i dx j • • • dx n (6)
The inter-image pair-wise comparisons are performed only between the immediate neighbors in the temporal domain.

Deformation Smoothness Constraints: In order to impose local smoothness to the deformation fields, we used the same smoothness constraint used in [START_REF] Sotiras | Graphical models and deformable diffeomorphic population registration using global and local metrics[END_REF], which is performed by penalizing the magnitude of displacement vector differences. This constraint was applied to neighboring nodes that belong to the same deformation field assuming a 26-nodes intra-image connectivity.

Implementation details: For all registrations we applied a multi-resolution scheme with 4 resolutions levels for the deformation grid. The spacing was set to 36 mm for the coarse resolution and was reduced in each iteration, ending up to 5.5 mm for the dense resolution. Five iterations were performed at each level. With each iteration the radius of candidate displacements was reduced by a factor of 0.8.

EXPERIMENTAL VALIDATION

In order to validate the proposed registration method, we used a data set provided by the Centre Hospitalier Universitaire Henri-Mondor, France. The data set consisted of 3D images of 25 different patients with lymphoma, scanned with a 3 Tesla MRI scanner using 3 b values: 50,400 and 800 (s/mm 2 ). Therefore, the sequence of the images to be registered was always of size 3 (n = 3) per patient. Both male and female patients were included. The image size was 130x106x22 for every patient, with a pixel resolution of 2.46x2.46 mm 2 in the axial plane and 5 mm 2 in the zaxis. Each image had been manually annotated by experts to indicate the contour of the tumor.

We compared the proposed method against two other modelling criteria: the sum of absolute intensity differences (SAD), which is appropriate when the intensity values of the deforming images come from the same distribution (under the presence of noise), and an entropy-based criterion, applicable in multimodal image registration. These two criteria were implemented in Eq.( 4) and Eq.( 6). Specifically, for the latter criterion, the entropy estimator based on spacings [START_REF] Bouzebda | New entropy estimator with an application to test of normality[END_REF] was used in Eq.( 4), whereas mutual information (MI) was used in the inter-images comparisons (Eq.( 6)). The entropy-based criterion is expected to be more appropriate than the SAD criterion for DW-MRI registration since intensity is not constant over time, i.e. for different b values. Fig. 1 illustrates an example of the DW-MR intensity variation of the three different b value images. The performance of the three registration approaches was quantitatively assessed by computing three different validation criteria: the SSE of the physiological model deriving the ADC, the MI and the normalised cross-correlation (NCC), averaged across pairs of the population (see Fig. 2, first row). The proposed method outperfomed the two other registration implementations in all three validation criteria, irrespectively of the similarity criterion used within the other frameworks. This result can be justified by the fact that the other methods do not use the right data term that models the temporal intensity change; thus a suboptimal solution is reached. Moreover, in order to qualitatively evaluate the results, the ADC map is shown for the same axial view of the patient used in Fig. 1, as derived by the three registration scenarios, as well as without registration (common in clinical practice). It can be seen that the ADC map calculated after registering the DW-MR images with our proposed method reveals anatomical structure and highlights the spatial tissue heterogenity. Part of our current work includes the simulation of diffusion-weighted images, which will allow us to better validate and compare the different approaches.

CONCLUSION

In this paper, we propose a novel deformable registration method that derives an optimal representation for shape and intensity change towards more accurately estimating the underlying pathology. Experimental results on DW-MR images demonstrate that the derived by this approach ADC map preserves structure, as it can be seen inside the big cyst superior to the tumor shown in the figures, and reveals tissues boundaries (tissue heterogenity, anatomical structure). A limitation of the proposed method is that it doesn't examine simultaneously the interaction among all the deformation fields of the population towards a global consistency on the physiological model. Such a global term would require the adaptation of high order inter-image cliques. Therefore, a group-wise registration method that includes high-order dependencies among the images would be a promising research direction for future work.

Fig. 1 :

 1 Fig. 1: An axial slice of the 3D image of a patient, with b values of 50, 400 and 800 s/mm 2 (from left to right). The red contour in the "b 50" image corresponds to the tumor's contour, as it was drawn by the expert.

Fig. 2 :

 2 Fig.2: Results on three different criteria for validating each method, applied on DW-MR images. First row: Boxplots of statistics (mean) of SSE by fitting the physiological model, MI and NCC respectively. Second row: ADC map (same axial slice as in Fig.1) derived without registration and with the examined registration methods.