Benoît Guédas
email: benoit.guedas@irccyn.ec-nantes.fr

Philippe Dépincé
email: philippe.depince@irccyn.ec-nantes.fr

A compromise definition in multiobjective multidisciplinary design optimization

Keywords: multidisciplinary, multiobjective, optimization, compromise

The decomposition of a global problem into several sub-problems, that can have antagonist goals, requires to find trade-off solutions. Moreover, the sub-problems are often multi-objective optimization problem : the disciplines have several antagonist objectives to optimize simultaneously. Thus, trade-off solutions have to be found both at the discipline level and at the multidisciplinary level. One way to consider the compromise is to compute all the Pareto efficient solutions of the multi-objective problem involving all the objectives of the disciplines at the same time. The optimal solutions can then be defined as the direct product of each disciplines partial ordered set by Pareto dominance relation. Unfortunately, this definition of the compromise is not satisfying. Indeed, information about efficiency of the solutions inside each discipline is lost during the direct product. In this paper, we propose extensions of the partial ordered set defined by the Pareto dominance relation in each discipline that keeps this information. Another dominance relation over disciplines is also presented. 2.

Introduction

The growing complexity of industrial products have led to specialization and distribution of knowledge, tools and work sites. The designer have now to face two main problems : the optimization of the product quality on a first hand, and the organization of work on a second hand. Multidisciplinary design optimization deals with all the problems that can arise in this context. For example, the optimization of the wing of an aircraft is both an aerodynamic and a structural problem in which each discipline have its own objective and uses specific tools, but a global solution have to be found. The decomposition of a global problem into several sub-problems, that can have antagonist goals, requires to find trade-off solutions. Moreover, the sub-problems are often multi-objective optimization problems: the disciplines have several antagonist objectives to optimize simultaneously. Thus, trade-off solutions have to be found both at the discipline level and at the multidisciplinary level.

In multiobjective multidisciplinary related papers, the compromise is not clearely defined and corresponds in fact to the Pareto optimal solution of the global problem including all the objectives of all the disciplines at the same time. This compromise is relevant if the decomposition of the problem in several disciplines is only due to practical organisational aspects. In the case where the decomposition in disciplines corresponds to a hierarchical organisation of the decision process, the global Pareto front may not represent the expeted compromise. For instance if there exists a subset of solutions that are optimal for all the disciplines at the same time, the global Pareto front will not be this subset of solutions but a larger set including solutions that are not optimal for all the disciplines. We propose a compromise method which in this situation only gives the subset of solutions that are optimal for all the disciplines at the same time. More generally, the proposed compromised method is not a compromise among the objectives but among disciplines.

As we are looking for a compromise method among disciplines which are optimization problems, the fourth section defines what are the solution of an optimization problem -and in particular a multiobjective optimization problem. This section introduces the notations and concepts used along this paper. The solutions are defined using ordered sets theory. The fifth section then defines our compromise method using the orederd set definitions and operations presented in the previous section. Finally, properties and examples of the proposed compromise are given in the sixth section before conclusions.

In some multidisciplinary design optimization problems, some variables are shared bewteen several disciplines while others are local to a given discipline. Moreover, many of these problems are coupled: inputs of one discipline are outputs of another discipline and vice versa. As this paper is a first attempt to define compromise bewteen disciplines, we only consider a simplified problem were all the disciplines share the same set of design variables and have no coupling functions.

Solutions of an optimization problem

Before defining a compromise between optimization problems, we need to define what are the solutions of an optimization problem. The following definitions are common to all optimization problems. According to Matthias Ehrgott [START_REF] Ehrgott | Multicriteria Optimization[END_REF], an optimization can be described by five elements:

• the feasible set X

• the objective space E

• the objective function vector f := (f 1 , . . . , f p ) : X → E

• the ordered set (F, ≤)

• the model map θ : E → F An optimization problem is then caracterized by: (X , f, E)/θ/(F, ≤) Most of the multidisciplinary design optimization methods such as MDF, IDF, AAO, CO, CSSO, BLISS, ATC,. . . only consider disciplines having one objective. In this case, in each discipline E = R, θ = id (the identity function) and (F, ≤) = (R, ≤) the canonical ordered set over the real numbers. Other methods such as E-MMGA, MORDACE and COSMOS are specially designed to handle multiobjective disciplinary problems. Here E = R n with n ≥ 2 and the solutions are compared with the Paretodominance relation . Classical method can also be applied to multiobjective problems if we take a model map θ : R n → R, such as, for instance, the weighted sum : (x 1 , . . . , x n ) → n i=1 ω i x i with the weights ω i ∈ R.

The set of solutions of an optimization problem will depend on the choice of the ordered set in which the solutions are compared.

Partial ordered sets

We give here some definitions about partial ordered sets which will be usefull for the definition of the compromise. More definitions can be found in [START_REF] Schröder | Ordered Sets: An introduction[END_REF]. To define the minimum of a set, we need to compare pairwise all the points of the set in order to find which of them is minimal. This comparison is an order relation and is defined as follows :

Definition An order relation ≤ on a set E is a binary relation verifying the following properties :

i) reflexivity: ∀x ∈ E, x ≤ x ii) antisymetry: ∀x, y ∈ E, x ≤ y ∧ y ≤ x ⇒ x = y ii) transitivity: ∀x, y, z ∈ E, x ≤ y ∧ y ≤ z ⇒ x ≤ z (E, ≤
) is then called an ordered set, partial ordered set or poset.

An ordered set in which the elements are comparable given an order relation. Unfortunately, it does not mean that all the elements of the set are comparable but only a subset of them.

Definition Let e 1 et e 2 be two posets (E, ≤). We say that e 1 and e2 are comparable, and we write e 1 ∼ e 2 et e 2 ∼ e 1 if and only if e 1 ≤ e 2 ou e 2 ≤ e 1 . We write e 1 ≁ e 2 if e 1 and e 2 are not comparable.

A set where all the elements are comparables is said totally ordered. (R, ≤) is a totally ordered set.

Definition Let (E, ≤) be an ordered set. ≤ is said total if ∀(e 1 , e 2 ) ∈ E 2 , e 1 ≤ e 2 ∨ e 2 ≤ e 1 .
To each order relation, we can associate its strict relation by removing the reflexivity property. In such a relation, an element cannot be compared to itself.

Definition A strict order relation < on a set E is a binary relation verifying the following poperties :

i) irreflexivity: ∀x ∈ E, x ≮ x ii) transitivity: ∀x, y, z ∈ E, x < y ∧ y < z ⇒ x < z (E, <) is then called a strict ordered set.
Some operations such as the sum and the product can then be defined on the ordered sets.

Definition Let O 1 = (E 1 , ≤ E1 ) and O n = (E n , ≤ En ) n two ordered sets, et E = n i=1 E i . We define the product O = n i=1 O i such that : e ≤ E e ′ ⇔ ∀i ∈ {1, . . . , n}, e i ≤ Ei e ′ i Definition Let Q the chain C h = y 1 < . . . < y h of h elements. Let h posets P i = (X i , ≤ i ).
We say that P = (X, ≤ P ) is the linear sum P i and we write

P = h i=1 P i with : X = h i=1 X i and a ≤ P b ⇐⇒ a ∈ X i , b ∈ X j and a ≤ i b if i = j a ∈ X i , b ∈ X j and y i < Q y j otherwise
In an optimization problem, we are looking for the minima of a poset.

Definition Let (E, ≤) be a poset. An element e * ∈ E is said to be minimal iff ∀e ∈ E withe * ∼ e, e * ≤ e. The set of all minimal elements of E is then written min E.

The minimal elements of (E, ≤) can be defined with these to equivalent definitions :

min E = {e * ∈ E : ∀e ∈ E s.t. e * ∼ e, e * ≤ e} (1)
min

E = {e * ∈ E : ∀e ∈ E, e ≤ e * ⇒ e * = e} (2)
We can now define the rank of an element in a poset which is an application from E into N that to each element of a partial ordered set assigns an integer.

Definition Let P = (E, ≤) a poset. For x ∈ X, we define the rank r(x) of x such that • r(x) := 0, if x is a minimum.

• r(x) := n, if the elements of rank < n have been assigned and x is a minimum in the poset P \{y ∈ P : r(y) < n}.

The rank can also begin to 1 instead of 0.

Optimal and efficient solutions

Let X be the feasible set of the decision space of an optimization problem, and Y be its image under the objective function mapping f :

X → R n . Y is thus a subset of R n . Y := {f (x) : x ∈ X } (3)
The solutions of a multiobjective optimization problem are defined with the Pareto-dominance relation ≻ which represent a preference. For a minimization problem we have the following definition of the Pareto-dominance relation:

Definition Let a, b ∈ Y ⊆ R n ,
we say that a Pareto-dominates b and we write a ≻ b if a is lower than b on each component and if there is at least one component which is strictly lower :

a ≻ b ⇐⇒ ∀i ∈ {1, . . . , n} a i ≤ b i ∃j ∈ {1, . . . , n} a j < b j
The Pareto-dominance relation allows us to define the non dominated set Y N ⊆ Y which is the set of all the points of the objective space which are not dominated.

Y N := {y * ∈ Y : ∄y ∈ Y s.t. y ≻ y * } (4)
The set of the solutions of the problem is called the efficient set X E ⊆ X and is defined as the set of points of the feasible set for which their image in the objective space is not dominated.

X E := {x ∈ X : f (x) ∈ Y N } (5)
We can notice that the Pareto-dominance relation is similar to the definition of the product of orders but without the reflexivity property (because of the strict order relation). Thus, the Pareto dominance relation is the strict order associated to the direct product. For this reason, we will now us the symbol instead of ≻ for the Pareto-dominance relation.

The solutions of a multiobjective minimization problem can then be interprated as the minima of Y with the order relation of the product (R, ≤) n :

Y N := {y * ∈ Y : ∀y ∈ Y s.t. y * ∼ y, y * ≤ y} (6) Y N := {y * ∈ Y : ∀y ∈ Y, y ≤ y * ⇒ y * = y} (7)
The two above definitions of Y N are equivalent to the previous one. Y N is also called the Pareto front and X E the Pareto set .

Compromise between optimization problems

Naïve definition

Let X 1 E and X 2 E be the respectively the efficients sets of the first and the second discipline, and let Y 2

N and Y 2 N be the corresponding non-dominated sets. Let X E the efficient set of the multidisciplinary problem and X G E the efficient set of the global multiobjective optimization problem including all the objectives. We would like that the efficient soltions of the problem have at least these two properties:

1. a compromise solution belongs to the Pareto front of the global optimization problem, 2. if there are efficient solutions that are common to all disciplines at the same time, then this set of solution is exactly the set of efficient solutions of the multidisciplinary problem.

This can be translated as:

1. X E ⊆ X G E 2. X 1 E ∩ X 2 E = ∅ ⇒ X E = X 1 E ∩ X 2 E

The product for compromise

The compromise solutions between n mono-objetive problems is the Pareto set. In other words, the compromise is defined as the minimal solutions of the product of the natural order of R: (R, ≤) n .

We can define the compromise between p disciplines which have q i , i ∈ {1, . . . , p} objectifs each the same way: by the product of the p orders (R qi , ≤ R q i ) of each disciplines. Unfortunately, if we sum all the objectives of all the disciplines and call n the result: n = p i=1 q i , we have:

(R n , ≤ R n ) = p i=1 (R qi , ≤ R q i ) (8)
This equation means that the result of the product of the p disciplines is the same than the product of the n objectives of the problem.

Extension of the direct product

The previous definition is not satisfying because the information on the efficiency of the solutions are lost during the product. Indeed, we can consider that de non-diminated points of each disciplines are better than the others. But a non-dominated point is not always comparable to a dominated point. Then, during the product, if a non-dominated point and a dominated point are not comparable, they are potentially on the global Pareto front if they are also incoparable in other disciplines. We want the optimal solutions of a discipline be better than all the other solutions. We can make the two following comments:

1. all the optimal solutions of a discipline are equivalent, 2. an optimal solution is prefered to a non-optimal solution.

Let Y D := Y\Y N be the set of dominated points of a discipline. We can define the new order relation ≤ ′ including the two previous comments:

1. ∀(y 1 , y 2 ) ∈ Y, y 1 ≤ y 2 ⇒ y 1 ≤ ′ y 2 2. ∀(y 1 , y 2 ) ∈ Y N × Y D , y 1 ≤ ′ y 2
The second point is useless because if all the non-dominated points are equal, by transitivity, they dominate all the other solutions. If the first point is removed and the second one kept, the global Pareto front after the product would be the union of each Pareto front.

We can also extend this new order relation by applying the same argument to all the ranks of the ordered set. We define the new relation ≤ ′′ as follows:

a ≤ ′′ b ⇐⇒ r(a) ≤ r(b)
Many other definitions of the "rank" have been proposed in the litterature related to multiobjective genetic algorithm. Goldberg [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF] first introduced the rank to bias the selection operator based on the evaluations of the objective function. This idea has then been used by Srinivas and [11] Deb in the NSGA algorithm in order to sort the individuals for a multiobjective problem. Other sorting method have been proposed for multiobjective genetic algorithm. But even if they use the term "rank", they do not correspond to the definition of the rank on a poset. That is why we introduce two new definitions: classification and pseudo-classification to name this mappings.

Definition A classification is a mapping c from the set X into N * such that:

• if X = ∅, ∃x ∈ X s.t. c(x) = 1 • ∀x 1 ∈ X s.t. c(x 1 ) > 1, ∃x 2 ∈ X s.t. c(x 2 ) + 1 = c(x 1 )
Definition A pseudo-classification is a mapping ϕ from the set X into N * such that:

• if X = ∅, ∃x ∈ X t.q. ϕ(x) = 1 • ∀x 1 ∈ X s.t. ϕ(x 1 ) > 1, ∃x 2 ∈ X s.t. ϕ(x 2 ) + 1 ≤ ϕ(x 1 )
The pseudo-classification proposed by Belegund and Salagame [START_REF] Belegund | Optimization of laminated ceramic composites for minimum residual stress and cost[END_REF] assigns the value 0 to each element of the non-dominated set and 1 to each element of the dominated set. The pseudo-classification proposed by Fonseca and Fleming [START_REF] Fonseca | Genetic algorithms for multiobjective optimization: Formulation discussion and generalization[END_REF] assigns to each element of the objective space the number of individuals it dominates plus one. In ordered sets theory, it corresponds to the cardinality of the ideal. Alberto and Ascarate [START_REF] Alberto | Multiobjective evolutionary algorithms. pareto ranking[END_REF] also proposed two pseudo-classifications. The first one is based on the number of times an element dominated others with random weighted sums, and the second one is a mix between the rank and Fonseca and Fleming's pseudo-classification.

All the above pseudo-classifications can be used to define the new order ≤ ′′ . We notice that the first extension ≤ ′ can be interpreted as the second one ≤ ′′ with taking Belegund and Salagame pseudoclassification instead of the rank.

The first compromise bewteen p disciplines is define as follows:

Y N := min p i=1 (Y i N ⊕ Y i D ) (9)
More generally, it can be defined as:

Y N := min p i=1 qi j=1 ψ j (Y i ) (10)
with ψ i a function that gives the part of a set which has all its elements scored i by the peudoclassement ϕ:

ψ i (E) := {e ∈ E : ϕ(e) = i} (11)
Of course, we can choose ϕ = r.

Examples and properties

Let us call c0 the compromise corresponding to the Pareto-optimal solutions of the whole problem, c1 the first compromise proposed and c2 the second one.

Examples

The second and third examples come from [START_REF] Dépincé | Collaborative optimization strategy for multi-objective design[END_REF]. The fourth and fifth example come from [START_REF] Engau | 2D decision-making for multicriteria design optimization[END_REF].

Example #1 Example #2 D 1    min xc,x1 f 1 (x c , x 1 ) = 10(x 1 -0, 1) 2 + (x c -3) 2 min xc,x2 f 2 (x c , x 2 ) = 10(x 2 -0, 3) 2 + (x c -α) 2 D 2    min xc,x3 f 1 (x c , x 3 ) = 10(x 3 -0, 3) 2 + (x c -β) 2 min xc,x4 f 2 (x c , x 4 ) = 10(x 4 -0, 5) 2 + (x c -9) 2 With α, β ∈ R Example #3 Objectives f 1 f 2 a (2,
D 1    min x1,x2 f 1 (x 1 , x 2 ) = (x 1 -2) 2 + (x 2 -1) 2 min x1,x2 f 2 (x 1 , x 2 ) = x 2 1 + (x 2 -3) 2 D 2    min x1,x2 f 1 (x 1 , x 2 ) = (x 1 -1) 2 + (x 2 + 1) 2 min x1,x2 f 1 (x 1 , x 2 ) = (x 1 + 1) 2 + (x 2 -1) 2
Under the following constraints:

     (x 2 1 -x 2 ) ≤ 0 (x 1 + x 2 -2) ≤ 0 (-x 1 ) ≤ 0 Example #4                      V (x) = L(2x 1 + √ 2x 2 + √ 2x 3 + x 4 ) d 1 (x) = F L E ( 2 x 1 + 2 √ 2 x 2 - 2 √ 2 x 3 + 2 x 4 ) d 2 (x) = F L E ( 2 x 1 + 2 √ 2 x 2 + 4 √ 2 x 3 + 6 x 4 ) d 3 (x) = F L E ( 6 √ 2 x 3 + 3 x 4 ) C                      x 1 ≥ F σ x 2 ≥ √ 2 F σ x 3 ≤ 3 F σ x 4 ≤ 3 F σ
with: F = 10kN, L = 200cm, E = 2.10 5 kN/cm 2 and σ = 10kN/cm 2 .The deproblem is decomposed into three subproblems as follows:

f 1 = (d 1 , V ), f 2 = (d 2 , V ), f 3 = (d 3 , V ) under the constraints C.

Inclusion

The order relation ≤ ′′ is an extension of ≤ ′ which is itself an extention of ≤. Thus, the solutions for the compromise c2 are all included in the solutions for the compromise c1 which are included in the compromise c0 (global Pareto front). We have the following inclusions:

min p i=1 qi j=1 ψ j (Y i ) ⊆ min p i=1 (Y i N ⊕ Y i D ) ⊆ min p i=1 Y i N (12)
Exemple Rabeau. N := f 2 (X 1 E ) and Y -2 N := f 1 (X 1 E ) the projection of each Pareto front in the other discipline, we can say that if there is a solution for which its image is above Y -1 N or Y -2 N in respectively the first and second discipline, then this solution is not a compromise according to our definition.

f (1,2) a 1 b 1 c 1 d 1 e 1 f 1 0 10 10 f (2,1) f (2,2) a 2 b 2 c 2 d 2
Proposition Let i and j two disciplines, i, j ∈ {1, 2}, i = j. If a solution x 1 ∈ X is such that ∃x 2 ∈ X j E such that f i (x 2 ) f i (x 1 ), this is not a trade-off solution.

Proof ∃x 2 ∈ X j E such that f i (x 2 ) f i (x 1 ) and x 2 ∈ X j E implies that ∄x 3 ∈ X such that f j (x 2 ) f j (x 3 ). Thus f j (x 2 ) f j (x 1 ) and so there is at least one solution x 2 that dominates x 1 .

But unfortunately, the opposite is not true: a solution for which both of its images are under Y -1 N and Y -2

N is not always a compromise.

Stability

Unfortunately, an order based on a pseudo-classification in a partially ordered set is not stable with the spread of the elements in the set. If some elements are added, the order can change as the Fig. 5 shows.

Figure 1 :

 1 Figure 1: Tous les points de Y N (rouges) sont meilleurs que ceux de Y D (bleus)

Figure 2 :

 2 Figure 2: Our compromise en problem #3

Figure 3 :Figure 4 :

 34 Figure 3: Exemple : tous les points sauf f sont efficaces.

Figure 5 :

 5 Figure 5: The order between e and f has been reversed.

June 1-5, 2009, Lisbon, Portugal

Conclusions


Wafa Rekik
email: wafa.rekik@gmail.fr

Dominique Béréziat
email: dominique.bereziat@upmc.frrapportde

Séverine Dubuisson

Keywords: 3D+t reconstruction, optical flow, spherical parametrisation Reconstruction 3D+t, flot opique, paramétrisation sphérique

published or not. The documents may come

Recovering 3D+t scenes from 3D and 2D+t data by motion compensation and using a prior surface model.

Research report LIP6

Recovering 3D+t scenes from 3D and 2D+t data by motion compensation and using a prior surface model.

Reconstruction de scènes 3D+t à partir de données 2D+t par compensation du mouvement sur un modèle de surface.

Introduction

In some computer vision applications, we deal with volumetric, i.e. 3D, acquisitions framing temporal, i.e. 2D+t, ones (Fig. 1). First type of acquisitions provides a description of the 3D scene geometry, thus, purely spatial information. Second ones observes 2D object motion providing temporal and partial spatial information. A possible approach to mine exhaustively these complementary datasets is to carry out a complete spatio-temporal or a 3D+t scene reconstruction. 3D+t sequences are restored by recovering the 3D original volume from each 2D frame belonging to the 2D+t sequence.

As the matter of fact, 3D reconstruction is an inverse problem. In this particular case, it is also an ill-posed problem since a single 2D frame is hardly sufficient to recover the 3D original volume. Indeed, this frame exhibits a great deal of spatial distortions and in-depth loss of information caused by the projective transformation that reduces the 3D real structure into a 2D image.

To handle this critical lack of information, we introduce a prior data model. It consists in a geometrical constraint defined on a spherical topology. Data observing locally spheric shaped objects can be issued from remote sensing acquisitions and especially in meteorology, from astronomical ones focusing on the Sun and in bio-cellular microscopic imagery. A relevant biologic application is described in [START_REF] Staneva | Phospholipase A2 promotes raft budding and fission from giant liposomes[END_REF] where biologists study dynamic of structures of interest evolving on a spherical cell simulation surface. In this context, for instance, spatial datasets consist in multi-focus acquisitions. Hence, they are not actual 3D images since they are formed by a series of equal sized slices, where signal in each slice is bothered by luminosity coming from adjacent ones. A spatial deconvolution stage is therefore needed when biologists calibrate correctly the video microscopy system to cell wall simulation experiments. Due to this limitation and to respect generality, we use synthetic data to validate our algorithm. Our approach uses this geometrical prior constraint in order to select a relevant reconstruction strategy. In the context of locally spheric shaped data, the information of interest is embedded on the outer curved surface. This allows the 3D-to-2D transformation to be modeled simply with a surface based reconstruction approach. On one hand, inverse 3D restoration is performed gradually by motion compensation involving a spatial regularity criterion that balances the in-depth loss of information. On the other hand, since movement has only two angular degrees of liberty, we opt for a spherical parametrisation to investigate motion in the surface domain instead of building a 3D model. For the latter, computation of spatial derivatives is incoherent and noisy, mainly with reference to the normal to the surface. Besides movement constraints are not relevant. We restore then spherical surfaces, parameterized by angular coordinates, that we display directly by an adapted visualisation tool, called MAPVIS [START_REF] Rekik | MAPVIS: a map-projection based tool for visualizing scalar and vectorial information lying on spheroidal surfaces[END_REF].

State-of-the-art

In the literature, a wide scope of approaches addresses the 3D reconstruction problem. Volume based methods aim at restoring an accurate description of luminosity in each voxel of the 3D reconstructed image, while surface based ones provide a 3D model of the object of interest. This model requires shape and texture recovery of the outer visible side of the observed object. Namely, tomographic reconstruction [START_REF] Natterer | The mathematics of computerized tomography[END_REF] focuses on restoring the original volume from multiple projective views. Stereo-vision [START_REF] Huei-Yung | Computer vision techniques for complete 3D model reconstruction[END_REF] approaches make use in general of the triangulation principle in order to recover the surface shape of selected targets. Moreover, 2D-to-3D registration methods re-construct 3D features from only one projective view. This registration consists in identifying a geometrical transformation that aligns 2D frames with 3D datasets often issued from a different modality of acquisition. This alignment determines their relative position and orientation and then locates them in the coordinate system of the 3D images. It is of growing interest in the medical imaging field mainly to register 3D pre-operative images (like MRI or CT images) into 2D series of intra-operative frames acquired during a surgical intervention. Registration may operate on feature-based methods or intensity based ones. Reader is referred to [START_REF] Gueziec | Anatomy-based registration of CT-scan and intraoperative X-ray images for guiding a surgical robot[END_REF][START_REF] Lemieux | a patient-tocomputedtomography image registration method based on digitally reconstructed radiographs[END_REF]. Tomography and stereovison, namely, require multiple projective views to achieve the reconstruction task. Therefore, approaches aligning 2D projection images to 3D datasets appear to be adapted to the context of our study. However these approaches fit into a registration framework rather than a reconstruction one, unless other projective views are available, which is quite unusual. They can be used alternatively for a preprocessing stage for some aforementioned applications. We then build an adapted reconstruction model involving a prior geometrical constraint.

Reconstruction by motion compensation

Our approach is based on the existence and the knowledge of a functional p that reduces every volumetric structure at each instant into a 2D frame. The central idea is to carry out the 3D+t reconstruction using a couple of 3D data acquired at two different instants as well as the 2D+t sequence acquired meanwhile. Since the intermediary video sequence presents an interesting temporal resolution, we propose to restore the underlying 3D structure of each of its frame. The latter provides only a partial description of the real 3D structure that we aim to estimate at a given instant t. We attempt to compensate this lack of information by the mean of a matching procedure with the 3D volumetric image estimated at t -1. For this purpose, we compute a displacement field matching both 3D structures. We introduce, consequently, a motion hypothesis expressing the assumption that a moving voxel keeps the same gray value over time. It is identical to the well known optical flow constraint [START_REF] Horn | Determining optical flow[END_REF]. It leads to a fundamental equation that we solve gradually in order to rebuild the whole 3D+t sequence, using the couple of 3D data as initial border conditions.

Modeling, assumptions

We carry out the 3D+t reconstruction using a series of 2D+t video sequence acquisitions intersected with 3D volumetric ones. We opt for the following notations for the remainder of the paper:

• X a 3D vector position, i.e. X = (x, y, z);

• x a 2D vector position, i.e. x = (x, y);

• W a 3D vector field, i.e. W = (u, v, w);

• I(X, t) a 3D+t sequence;

• I 2D (x, t) a 2D+t sequence;

• Î(X, t) the estimation of I at (X, t);

Concretely, we have two 3D structures at two different instants: I 1 = I(X, t 1 ) and I 2 = I(X, t 2 ). Between t = t 1 and t = t 2 , we also dispose of a video sequence I 2D (x, t) describing all the 2D reductions of the I(X, t) structures evolving in this duration. Each 2D frame I 2D is a projective sight, at a given instant t, of a 3D original volume. This projection only depends on the data acquisition mechanism. We model it by a linear and stationary transformation, called projection p, given by:

p(I)(x, y) = Ê I(x, y, z)h(z)dz, (1)
where h(z) is a function describing the interaction between observed objects and the input light signal. p integrates contributions of all 3D slices to produce a unique 2D image. It introduces then an important loss of information in the z-direction. In order to reconstruct gradually the 3D+t sequence, we estimate the 3D structure at a given instant t, using the data I 2D (x, t) extracted from the 2D+t sequence. We then minimize the discrepancy between p( Î) and I 2D . A straightforward data model equation is:

p( Î)(x, t) -I 2D (x, t) = 0, ∀t ∈ [t 1 , t 2 ] (2)
Since objects of interest evolve slowly from t to t -1, we can match Î(X, t) with the 3D structure I(X, t -1) known at the previous instant. This matching is quantified in terms of a displacement vector field W between both images. It expresses the assumption that a moving voxel keeps the same gray value over time yielding the following motion constraint, similar to optical flow one:

Î(X, t) = I(X + W, t -1) (3)
We estimate a retrogress vector field defining voxel displacement from t -1 to t for notation conveniences (we represent -W in Section 4). Expanding the left-hand side of Equation (3) in a first order Taylor series leads to:

Î(X, t) = ∇I(X, t -1).W + I(X, t -1), (4)
where ∇ = ( ∂ ∂x , ∂ ∂y , ∂ ∂z ) is the gradient operator. Since p is a linear operator, injection of the Equation (4) in the earlier data model constraint [START_REF] Gueziec | Anatomy-based registration of CT-scan and intraoperative X-ray images for guiding a surgical robot[END_REF] provides the 2D data motion constraint related to W: p(∇I(X, t -1).W) + p(I(X, t -1)) -

I 2D = 0 (5)
We compute the velocity vector field W at each instant t.

Formulation in the spherical context

Let us point out that, in the context of spherical shaped data observations, structures of interest evolve on the surface of a ball of constant ray. Even-though the support is embedded on a 3Dspace, movement has only two degrees of liberty. It is then relevant to estimate the vector field using a spherical surface parametrisation. Let's establish equations formulated in the previous section in the spherical context. For this purpose, we introduce two 2D images I s and I s 2D defined by the following substitutions (see Fig. 2):

I s (θ, ϕ, t) = I(x(θ, ϕ), y(θ, ϕ), z(θ, ϕ), t) (6) I s 2D (θ, ϕ, t) = I 2D (x(θ, ϕ), y(θ, ϕ), t) (7)
with the following parametrisation:

x(θ, ϕ) = R cos(ϕ) sin(θ) y(θ, ϕ) = R sin(ϕ) sin(θ) z(θ, ϕ) = R cos(θ)
The spherical coordinates ϕ and θ respectively stand for the longitude and co-latitude angles. Equation ( 6) expresses merely passage from cartesian coordinates to spherical ones. Equation [START_REF] Natterer | The mathematics of computerized tomography[END_REF] involves the mapping of the 2D frame texture into a hemisphere with same diameter that the 3D known balls. Data model constraint (2) becomes:

p( Îs )(θ, ϕ, t) -I s 2D (θ, ϕ, t) = 0 (8)
Let's now set up p in the sphere reference mark schemed in Fig. 2. Operator p integrates, in a plane including both geographical poles, contributions of frontal and dorsal hemispheres. Assuming that frontal and dorsal interactions with the input light signal are the same, p is defined by:

Y ′ C ϕ X ′ Y ′ M1 M2 π -ϕ P Z ′ Z ′ Y X Z O
p(I s )(θ, ϕ) = I s (θ, ϕ) + I s (θ, π -ϕ) (9)
Matching between Îs (., t) and I s (., t-1) are now expressed as performed in the cartesian parametrisation by Equation (3). According to Fig. 2 notations, the position of a moving point on the spherical surface is determined by the couple of angles (θ, ϕ) and the following equality :

--→ CM = R -→ e r ,
where R is the sphere radius and ( -→ e r , -→ e θ , -→ e ϕ ) the Frenet local coordinate system related to M . Thus:

Îs (R -→ e r , t) = I s (R -→ e r + ∂R -→ e r , t -1), (10)
where -→ e r +∂R -→ e r is the elementary displacement at R -→ e r after the duration ∂t. Since R is constant, ∂ -→ e r is defined by :

∂ -→ e r = ∂θ∂ -→ e θ + sin(θ)∂ϕ∂ -→ e ϕ (11)
Inserting right-hand side of Equation [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF] in the equality (10) leads to, for a unit sphere (R = 1):

Îs ( -→ e r , t) = I s ( -→ e r + ∂θ∂ -→ e θ + sin(θ)∂ϕ∂ -→ e ϕ , t -1) (12)
The right member of Equation ( 12) is linearized using a first order Taylor expansion. It yields the following equation:

Îs (θ, ϕ, t) = I s t-1 (θ, ϕ, t -1) + ∂I s (θ, ϕ, t -1) ∂θ ∂θ ∂t + sin(θ) ∂I s (θ, ϕ, t -1) ∂ϕ ∂ϕ ∂t (13)
Inserting second member of Equation ( 13) in [START_REF] Rekik | MAPVIS: a map-projection based tool for visualizing scalar and vectorial information lying on spheroidal surfaces[END_REF]

Numerical Resolution

We compute an angular velocity field defined by the couple of angular variations ( θ, φ):

θ = ∂θ ∂t φ = ∂ϕ ∂t ( 16
)
For this sake, we have to solve the motion constraint given in (15). Let us introduce the aforementioned p expression in Equation (15):

p ∂I s (., t -1) ∂θ ∂θ ∂t + sin(θ) ∂I s (., t -1) ∂ϕ ∂ϕ ∂t = ∂I s (θ, ϕ, t) ∂θ θ(θ, ϕ) + ∂I s (θ, π -ϕ, t -1) ∂θ θ(θ, π -ϕ)+ sin(θ)( ∂I s (θ, ϕ, t -1) ∂ϕ φ(θ, ϕ)+ ∂I s (θ, π -ϕ, t -1) ∂ϕ φ(θ, π -ϕ)) (17)
where ( θ(θ, ϕ), φ(θ, ϕ)) is the couple of angular velocity fields on the frontal hemisphere and ( θ(θ, πϕ), φ(θ, πϕ)) the couple of fields on the dorsal one. Equation (15) involves four velocity fields. This is due to the transformation p that introduces a positioning ambiguity. Indeed, it is not possible to undoubtedly state if the original position of a given point on the projection plane is in the frontal or the dorsal hemisphere. In order to reduce variable space, we suppose that these velocity fields are related. Let's express relation between velocity fields on the frontal as well as on the dorsal hemisphere. This relation takes into account possible object occlusions. In fact, objects with same latitude and a double-quadrant adjusted longitude are mapped at exactly the same position in the projection plane (see Fig. 3). To circumvent this ambiguity of positioning, we assume that: Under this assumption, number of variables is reduced to two. We estimate then velocities on the upper hemisphere, where θ has a range from 0 to π and ϕ from -π 2 to π 2 . For the remainder of the globe, displacement vectors are issued from the system of equations (18). Data motion constraint is merely respect to v = ( θ, φ):

θ(θ, π -ϕ) = θ(θ, ϕ) φ(θ, π -ϕ) = -φ(θ, ϕ) (18)
Y ′ Z ′ C ϕ X ′ Y ′ Z ′ P M 1 M 2 π -ϕ
Å(v)( θ, φ, t) = p( ∂I s (., t -1) ∂θ ) θ+ sin(θ) ∂I s (., t -1) ∂ϕ (θ, ϕ) - ∂I s (., t -1) ∂ϕ (θ, π -ϕ) φ +I ′ s 2D = 0 (19)
Since this equation involves only two variables, we can build a bi-dimensional model to accurately recover motion. In order to estimate a dense motion field we use a variational formulation. We build a functional whose minimum, with respect to v, corresponds to the solution of Equation ( 19):

E 1 (v) = π 0 π 2 -π 2 Å(v) 2 dϕ dθ
Since Equation ( 19) is under-determined, it is necessary to use an additional constraint to solve it. An obvious way is to constrain the spatial variations of v using a L 2 norm [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF]. This regularization penalizes high spatial deformations of v and implicitly important variations with reference to the data known at a previous time:

E 2 (v) = π 0 π 2 -π 2 ∇ θ 2 + ∇ φ 2 dϕ dθ
Other types of regularization may be used (L 1 norm, image driven norm, ...). Finally, our functional E, is, with respect to v:

E = E 1 + αE 2 ,
where α is a tuning parameter between the two terms that weights the importance of the regularization term E 2 . Differentiation of E, with reference to θ and φ yields a set of Euler-Lagrange equations [START_REF] Glowinski | Numerical Methods for Nonlinear Variational Problems[END_REF]. Discretization of the latter with finite differences leads to a linear system, that we solve using a Gauss-Seidel iterative scheme. Our goal, as schemed in Fig. 1, is to recover a full 3D+t sequence from the 2D+t series of data and the couple of 3D samples framing the temporal interval. The central idea is to estimate gradually velocity field between each 2D frame extracted from the 2D+t at a given instant t and the 3D structure previously estimated at t -1, as detailed in Section 3.1. Under the constant brightness assumption from time t -1 to t, an estimation of Îs (., t) is recovered from the data I s (., t -1) and angular variations ( θ, φ) by the following:

Îs (θ, ϕ) = I s (θ + θ, ϕ + φ, t -1) (20)
However angular coordinates θ + θ and ϕ + φ, computed from discrete 2D couples, do not generally correspond to grid-positions in the output sampled image sequence I s : they yield a set of scattered data. We attempt, therefore, to fit a smooth surface through the non-uniform distribution of these data samples. For this purpose, we represent our surface as a sum of weighted and shifted synthesis functions. B-splines stand apart for a good compromise between quality (namely high fidelity reconstruction and regularity) and computational issues, reader is referred to [START_REF] Thévenaz | Image interpolation and resampling[END_REF]. We then use the method introduced in [START_REF] Lee | Scattered data interpolation with multilevel B-splines[END_REF] and dedicated to recover such samples. The implemented algorithm makes use of a coarse-to-fine hierarchy of control lattices in order to generate a sequence of bi-cubic B-spline functions whose sum approaches the desired interpolation function. In order to restore gradually the complete 3D+t sequence, we adopt a quite simple reconstruction strategy. We use, respectively, a forward progression procedure taking as first 3D available image, I 1 and a backward one using the final 3D available image I 2 . Consequently, we generate two displacement vector fields at the median moment of the temporal interval [t 1 , t 2 ]. The final 3D image is computed merely using the average of both fields. Let us point out that we restore surfaces parametrized by angular coordinates (θ, ϕ). To display results, we use an adapted visualisation tool called, MAPVIS [START_REF] Rekik | MAPVIS: a map-projection based tool for visualizing scalar and vectorial information lying on spheroidal surfaces[END_REF], instead of carrying out inverse transformations in cartesian 3D+t grids. This tool displays complex information (scalar and vectorial) lying on real spheroid surfaces or projected ones. In the 3D case, MAPVIS projects data embedded on the actual globe observed part. Moreover, it improves visualisation of texture in perspective sight views of spheric shaped structures. This tool is based on a suitable planar map projection that unrolls the curved surface around a given origin. Therefore, varying the projection origin around the surface allows to observe different views of the sphere. Since the selected map projection minimizes distortions around the projection origin, the closer this point the more accurate the data recovery. Equatorial aspect of map projection, with reference to an origin lying on the equator displays the frontal hemisphere of the observed globe and polar one, with reference to the north geographical pole displays information lying on the northern hemisphere.

Results

We present, in this Section, a series of experiments in order to assess the performance of the proposed approach. They consist in computing the displacement vector field matching a 2D frame with an anterior 3D structure. As real data are not available, we have designed two sets of simple synthetic data. First one is composed of a couple of spheres I(., t 1 ), I(., t 2 ) holding two squares lying on respectively the frontal and dorsal hemispheres. From t 1 to t 2 , the frontal square center moves towards the right-bottom direction and the dorsal square evolves merely downwards. We simulate the 3D-to-2D transformation of I(., t 2 ) yielding I 2D (., t 2 ) and we compute the spherical representation of I and I 2D using Equations [START_REF] Lemieux | a patient-tocomputedtomography image registration method based on digitally reconstructed radiographs[END_REF][START_REF] Natterer | The mathematics of computerized tomography[END_REF]. We then determine angular variation fields ( θ, φ) matching data I s (., t 1 ) with the projection I s 2D (., t 2 ) in order to estimate Îs (., t 2 ) by motion compensation. Second set of data is quite identical to the first one, except that squares occlude each other in the 2D projection. Results are displayed in cartesian coordinates using MAPVIS. We display in Fig. 4 and 5 motion estimation results for respectively first and second set of synthetic data. Let us point out that motion computation produces a set of correct vector fields (from amplitude and direction point of views), perfectly tangent to the enclosing spherical surface. Moreover, our reconstruction method handles ambiguities, due to the 2D data position, by determining which hemisphere it is evolving in. It is noticeable by visualisation of the difference between the estimation of 3D image and the original one, i.e. I(., t 2 ) -Î(., t 2 ). Performance of our algorithm is not bothered by the more complex case showing an object occlusion (see Fig. 5) which does not ease the aforementioned structure distinction problem. However image difference shows also reconstruction errors, doubtless amplified by the double interpolation procedure involved, first by the computation of a smooth surface from samples recovered by motion compensation, second by the (a) From left to right: visualisation of the frontal then dorsal hemisphere of the set (I(., t 1 ), W). visualization of projective views [START_REF] Rekik | MAPVIS: a map-projection based tool for visualizing scalar and vectorial information lying on spheroidal surfaces[END_REF]. These errors are localized in the border of moving objects. They may be caused by spatial derivative computations. Besides, they are due to the smoothing of motion discontinuities introduced by the spatial regularization of the estimated displacement vector field.

Conclusion

In this paper, we have presented an original approach dedicated to 3D+t scene reconstruction using 3D temporal data and 2D+t temporal sequences. It is based on motion compensation and involves a prior surface model. The latter concerns spherical topology arising in number of computer vision applications. In this context, structures of interest evolve on a spherical support. Although movement is embedded on a 3D domain, it has only two degrees of liberty. Consequently, it is more interesting to investigate motion in the surface domain than to build a 3D model. Consequently, we build a 2D model using a spherical parametrisation in order to compute angular velocity fields. Three dimensional spherical shapes are gradually restored by estimating angular variations. Matching is carried out under the assumption that a moving voxel keeps the same brightness over time. This assumption corresponds to the optical flow constancy constraint. Our approach recovers reliably the 3D structures and handles possible object occlusions. Thanks to the 2D modeling, implemented algorithms are not prohibitive compared to a full 3D approach. Moreover, we use an adapted tool, with relevant visualisation properties, in order to quickly display the restored surfaces parameterized by spherical angular coordinates. Currently, the whole sequence is computed frame by frame, causing errors cumulation along the sequence. This limitation could be circumvented using a spatio-temporal smoothness constraint computed globally on the whole sequence [START_REF] Weickert | Variational optic flow computation with a spatio-temporal smoothness constraint[END_REF]. Discussing other reconstruction issues, it is possible to generalize the prior geometrical constraint to any regular surface. We propose also in another framework to solve the 3D+t scene recovery taken lower assumptions related to the data model. Moreover, we plan to apply our approach based on spherical topology to microscopic cell wall simulation acquisitions in the context described in [START_REF] Staneva | Phospholipase A2 promotes raft budding and fission from giant liposomes[END_REF]. That requires a spatial deconvolution framework to reconstruct each sphere from the original multi-focus image.

Figure 1 :

 1 Figure 1: Recovering a 3D+t sequence from a couple of 3D images framing a series of 2D projections.

Figure 2 :

 2 Figure 2: System of spherical coordinates. The projection plane corresponding to p is P = (C, Y ′ , Z ′ ).

Figure 3 :

 3 Figure 3: Occlusion of objects centered respectively in M 1 and M 2 .

  (b) From left to right: frontal view of projection I 2D (., t 2 ), frontal then dorsal hemisphere view of the estimated image Î(., t 2 ).

  (c) From left to right: frontal then dorsal hemisphere view of the difference between original image and estimated one I(., t 2 ) -Î(., t 2 ).

Figure 4 :

 4 Figure 4: Result on sequence without occluding squares.

  (a) From left to right: visualisation of the frontal then dorsal hemisphere of the set I s (., t 1 ), W.(b) From left to right: frontal view of projection I 2D (., t 2 ), frontal then dorsal hemisphere view of the estimated image Î(., t 2 ).

  (c) From left to right: frontal then dorsal hemisphere of the difference between original image and estimated one I(., t 2 ) -Î(., t 2 ).

Figure 5 :

 5 Figure 5: Result on sequence with occluding squares.


3while the later was showing that NMR spectroscopy was also able to analyze the structure and dynamics of high temperature molten materials 4 . High temperature was one of our main goals and J.P. Coutures had installed in Orleans an aerodynamic levitation device using a CO 2 laser as heating source. This device allowed to reach very high temperature (>2000°C), minimizing heterogeneous nucleation, and was also small enough to fit in a wide bore magnet NMR probe (Fig. 1). With the help and guidance of F.Taulelle in Paris and D.Muller at Bruker in Germany, we were able to show that the concept of laser heating could be adapted to the acquisition of NMR spectra with a setup able to reach temperature up to ~1300°C in boron nitride crucibles and up to more than 2000°C in aerodynamic levitation conditions. The demonstration experiment 5 , carried out in Paris at J.Livage's laboratory, made it possible to gather the necessary financial support and install a 300MHz MSL solid state NMR spectrometer in Orléans allowing our NMR expertise to develop.

From the first results obtained at high temperature it rapidly became clear that the motionally averaged picture of the high temperature molten state had to be complemented by an accurate description of the structure of materials at room temperature from the whole variety of observable nuclei, including quadrupolar nuclei like 27 Al. Aluminum spectra were easily observed at high temperature in the molten state 6,7 but much more complex for solids (static or magic angle spinning) because of the strong second order broadenings (Fig. 2). Based on the work of Samoson 2 I adapted the fitting program that I developed for the analysis of Mossbauer spectra during my PhD to account for the specificity of high resolution MAS NMR : spinning sidebands, second order quadrupolar patterns… This was the beginning of the well known dmfit program 8 that has today a few thousand users worldwide.

Figure 1 : High Temperature laser heating setup allowing heating up to ~2000°C under aerodynamic levitation 5,6,7 .

My first contact with the NMR community took place at the X th ISMAR conference held in Morzine in 1989. At that time we were only acquiring 1D spectra and it was a shock to learn about multiple dimension experiments in liquids and solids involving dipolar recoupling, double orientation rotation (DOR [START_REF] Wu | Double Rotation[END_REF] ) or dynamic angle spinning (DAS 10 )... Sometime later, in 1992, I had the opportunity to visit A.Pines' laboratory in Berkeley and met P.J. Grandinetti who was about to move to Ohio State University. Jointly with A. Llor we made plans for Phil to visit France with his DAS probe and this resulted in a long lasting and fruitful collaboration. Antoine and Phil were combining theoretical and practical skills of Berkeley and Saclay on quadrupolar nuclei and, from them, I understood how phase cycling could be used to select coherence pathways and that the central transition of the complex quadrupolar spin system could be manipulated as a fictitious I=1/2 spin… which still is the source of new developments. Thanks to this introduction we were ready when L.Frydman introduced MQMAS 11,12 in 1995 and among the first group to implement it, combining MQ excitation and central transition selective excitation 13 . I still have a very clear memory of the active discussions at the poster session of the Boston ENC in 1995 where L.Frydman first presented his experiment and the time spent to understand and explain triple quantum filtering. Among the various subject it is worth mentioning : the extension of PASS to the description of QPASS for quadrupolar nuclei 16 , the access to very high field (40T!) with Z. Gan and A. Samoson 17 , the development of J-coupling based approach for crystalline (with R.K. Harris 18 ), or for noncrystalline (with L.Emsley 19 We are now at a time where Solid State NMR is a well established tool for the characterization of materials but still in the situation were new developments and applications rapidly emerge opening new avenues in terms of observation scale, characterization of through bond and through space connectivity 26 , increased sensitivity, measurement on paramagnetic samples... This includes studying order and disorder in complex materials like glasses, energy related materials for battery or hydrogen storage, bioactive or biocompatible materials, zeolites and catalysts with number of pages still to be written and certainly more on subjects that we hardly imagine today.

Biographical Sketch

Figure 2 :

 2 Figure 2 : 2D TOP processing 14 of the 750MHz 1D 27 Al MAS NMR spectrum of YAG (Y 3 Al 5 O 12 ) showing the complexity and the richness of quadrupolar MAS NMR spectroscopy (CT central transition, ST Satellite Transitions). From these premises our activities developed around the description of new methods for application to quadrupolar nuclei and the study of inorganic or hybrid materials with numerous collaborations in France, Europe and worldwide. I would like specifically to acknowledge here the members of the Orleans group : C. Bessada, P. Florian, F. Fayon, B. Alonso, V. Montouillout, N. Pellerin, M. Deschamps, S. Cadars, A. Rakhmatullin, A.L. Rollet, the postdocs and the PhD students that made the group become what it is today with the constant support of J.P.Coutures, director of CRPHT and CRMHT. This was made possible with numerous collaborations involving Taulelle (Versailles France), F. Babonneau, C. Gervais, C. Sanchez and C. Bonhomme (Paris France), B. Bujoli & J.M. Bouler (Nantes France),

  Massiot. b 1957. Ecole Normale Supérieure Paris 1977-84, Thèse de 3 ème Cycle in Geochemistry 1983 (supervisor Pr. Michel Treuil) ; Chargé de Recherche CNRS CRPHT Orléans1984, Directeur de Recherche CRMHT Orléans 1996, Director of CEMHTI UPR3079 CNRS 2008-; Director of TGIR RMN THC FR3050 2007-; Silver Medal CNRS 2003 ; Approx. 215 publications, author of the dmfit program for the analysis of solid state NMR spectra. Current research specialty: multinuclear solid state NMR and application to inorganic and hybrid materials.

  I came to Solid State Nuclear Magnetic Resonance incidentally in the late 80s, attracted by the increasing capacity of this spectroscopy to characterize the local structure in materials lacking from long range order. As a materials science chemist, trained in geology or geochemistry, my main concerns were silicate glasses and melts. This was the time of the mile stone work by G. Engelhardt 1 , E. Lippmaa, A. Samoson 2 , A. Grimmer and coworkers in Eastern Europe or of J.F. Stebbins in Berkeley and Stanford California. The former were able to evidence directly local coordination of silicon or aluminum from 29 Si or 27 Al MAS NMR in silicate and alumino-silicate glasses

	Massiot, Dominique :
	Materials by NMR of
	Quadrupoles from room
	to high temperature
	Dominique Massiot
	CEMHTI-CNRS, Orléans, France

  ) phosphates, the study of alumino-silicate glasses of geological & industrial interest with D.Neuville & L.Cormier 20 , the study of biocompatible materials with application to the treatment of Osteoporosis with B.Bujoli & J.M.Bouler 21 , the study of sol/gel and hybrid materials with F.Babonneau, C.Sanchez, C.Gervais or C.Bonhomme at LCMCP Paris France 22,23,24,25 …

G.Engelhardt,

Engelhardt, Günter: Silicon-29 NMR of Silicates: From Liquids to Solids. An Historical Perspective. Encyclopedia of Magnetic Resonance 2 A.Samoson, 2007. Samoson, Ago: Development of High-Resolution NMR of Quadrupole Nuclei in Solids. Encyclopedia of Magnetic Resonance